Path-Integral Calculations of Nuclear Quantum Effects in Model Systems, Small Molecules, and Enzymes via Gradient-Based Forward Corrector Algorithms

Origin of Catalysis by Nitroalkane Oxidase

The rate of proton abstraction of the carbon acid nitroethane by Asp402 is accelerated by a factor of 10⁸ in the enzyme nitroalkane oxidase (NAO) relative to that by the organic base acetate ion in water. The Cα proton of nitroalkanes is known to exhibit an abnormal correlation between its acidity strength and the rate of deprotonation, with an unusually slow rate of deprotonation in water. This work examines the origin of NAO catalysis, revealing that the rate enhancement by the enzyme is due to transition-state stabilization, restoring the normal behavior of the linear free energy relationship of Bronsted acids. Interestingly, NAO employs the ubiquitous cofactor flavin adenosine diphosphate (FAD) to perform the subsequent oxidation. Does the FAD cofactor also affect the catalytic rate of the initial proton transfer process of the overall nitroalkane oxidation? Classical molecular dynamics and path-integral simulations using a reaction-specific combined quantum mechanics/molecular mechanics (QM/MM) approach were carried out to obtain the free energy reaction profiles, or the potentials of mean force, for the enzymatic reaction and for a model reaction in aqueous solution, as well as for the 2'-deoxy-FAD co-factor-modified NAO. Free energy perturbation calculations suggest that transition-state stabilization of the reactive fragment is the primary cause of the catalytic effect. It is found that the FAD cofactor plays a crucial role in increasing the Cα proton acidity, via specific hydrogen bonding and π-stacking interactions, although these factors have a smaller effect on the enhancement of the rate of deprotonation. Model QM calculations of the π-stacking complexes between the FAD isoalloxazine ring and the neutral and anionic nitroethane, respectively, reveal that the anionic π-stacking complex is more stable than the neutral one by 15.7 kcal/mol, and a net π-stacking energy of 17.3 kcal/mol is obtained. Hence, the isoalloxazine ring, in addition to serving as a very potent oxidizing agent via the formation of covalent intermediate structures, is able to exert a considerable amount of catalytic effect through noncovalent π-stacking interactions.

Generalized fourth-order decompositions of imaginary time path integral: Implications of the harmonic oscillator

The imaginary time path integral formalism offers a powerful numerical tool for simulating thermodynamic properties of realistic systems. We show that, when second-order and fourth-order decompositions are employed, they share a remarkable unified analytic form for the partition function of the harmonic oscillator. We are then able to obtain the expression of the thermodynamic property and the leading error terms as well. In order to obtain reasonably optimal values of the free parameters in the generalized symmetric fourth-order decomposition scheme, we eliminate the leading error terms to achieve the accuracy of desired order for the thermodynamic property of the harmonic system. Such a strategy leads to an efficient fourth-order decomposition that produces third-order accurate thermodynamic properties for general systems.

Temperature-Dependent Kinetic Isotope Effects in R67 Dihydrofolate Reductase from Path-Integral Simulations

Calculation of temperature-dependent kinetic isotope effects (KIE) in enzymes presents a significant theoretical challenge. Additionally, it is not trivial to identify enzymes with available experimental accurate intrinsic KIEs in a range of temperatures. In the current work, we present a theoretical study of KIEs in the primitive R67 dihydrofolate reductase (DHFR) enzyme and compare with experimental work. The advantage of R67 DHFR is its significantly lower kinetic complexity compared to more evolved DHFR isoforms. We employ mass-perturbation-based path-integral simulations in conjunction with umbrella sampling and a hybrid quantum mechanics-molecular mechanics Hamiltonian. We obtain temperature-dependent KIEs in good agreement with experiments and ascribe the temperature-dependent KIEs primarily to zero-point energy effects. The active site in the primitive enzyme is found to be poorly preorganized, which allows excessive water access to the active site and results in loosely bound reacting ligands.

The Impression of a Nonexisting Catalytic Effect: The Role of CotB2 in Guiding the Complex Biosynthesis of Cyclooctat-9-en-7-ol

Terpene synthases generate terpenes employing diversified carbocation chemistry, including highly specific ring formations, proton and hydride transfers, and methyl as well as methylene migrations, followed by reaction quenching. In this enzyme family, the main catalytic challenge is not rate enhancement, but rather structural and reactive control of intrinsically unstable carbocations in order to guide the resulting product distribution. Here we employ multiscale modeling within classical and quantum dynamics frameworks to investigate the reaction mechanism in the diterpene synthase CotB2, commencing with the substrate geranyl geranyl diphosphate and terminating with the carbocation precursor to the final product cyclooctat-9-en-7-ol. The 11-step in-enzyme carbocation cascade is compared with the same reaction in the absence of the enzyme. Remarkably, the free energy profiles in gas phase and in CotB2 are surprisingly similar. This similarity contrasts the multitude of strong π-cation, dipole-cation, and ion-pair interactions between all intermediates in the reaction cascade and the enzyme, suggesting a remarkable balance of interactions in CotB2. We ascribe this balance to the similar magnitude of the interactions between the carbocations along the reaction coordinate and the enzyme environment. The effect of CotB2 mutations is studied using multiscale mechanistic docking, machine learning, and X-ray crystallography, pointing the way for future terpene synthase design.

Enzymatic control of product distribution in terpene synthases: insights from multiscale simulations

In this opinion, we review some recent work on terpene biosynthesis using multiscale simulation approaches, with special focus on contributions from our group. Terpene synthases generate terpenes employing rich carbocation chemistry, including highly specific ring formations, proton, hydride, methyl, and methylene migrations, followed by reaction quenching. In these enzymes, the main catalytic challenge is not rate enhancement, but rather control of intrinsically reactive carbocations and the resulting product distribution. Herein, we review multiscale simulations of selected mono-, sesqui-, and diterpene synthases. We point to the many tools adopted by terpene synthases to achieve correct substrate fold, carbocation formation, carbocation reaction environment, and reaction quenching. A better understanding of the toolbox employed by terpene synthases is expected to aid in the search for new enzymatic and biomimetic synthetic routes to natural and unnatural terpenes.

Role of Microsolvation and Quantum Effects in the Accurate Prediction of Kinetic Isotope Effects: The Case of Hydrogen Atom Abstraction in Ethanol by Atomic Hydrogen in Aqueous Solution

Hydrogen abstraction from ethanol by atomic hydrogen in aqueous solution is studied using two theoretical approaches: the multipath variational transition state theory (MP-VTST) and a path-integral formalism in combination with free-energy perturbation and umbrella sampling (PI-FEP/UM). The performance of the models is compared to experimental values of H kinetic isotope effects (KIE). Solvation models used in this study ranged from purely implicit, via mixed–microsolvation treated quantum mechanically via the density functional theory (DFT) to fully explicit representation of the solvent, which was incorporated using a combined quantum mechanical-molecular mechanical (QM/MM) potential. The effects of the transition state conformation and the position of microsolvating water molecules interacting with the solute on the KIE are discussed. The KIEs are in good agreement with experiment when MP-VTST is used together with a model that includes microsolvation of the polar part of ethanol by five or six water molecules, emphasizing the importance of explicit solvation in KIE calculations. Both, MP-VTST and PI-FEP/UM enable detailed characterization of nuclear quantum effects accompanying the hydrogen atom transfer reaction in aqueous solution.

Accelerating equilibrium isotope effect calculations. II. Stochastic implementation of direct estimators

Path integral calculations of equilibrium isotope effects and isotopic fractionation are expensive due to the presence of path integral discretization errors, statistical errors, and thermodynamic integration errors. Whereas the discretization errors can be reduced by high-order factorization of the path integral and statistical errors by using centroid virial estimators, two recent papers proposed alternative ways to completely remove the thermodynamic integration errors: Cheng and Ceriotti [J. Chem. Phys. 141, 244112 (2015)] employed a variant of free-energy perturbation called "direct estimators," while Karandashev and Vaníček [J. Chem. Phys. 143, 194104 (2017)] combined the thermodynamic integration with a stochastic change of mass and piecewise-linear umbrella biasing potential. Here, we combine the former approach with the stochastic change in mass in order to decrease its statistical errors when applied to larger isotope effects and perform a thorough comparison of different methods by computing isotope effects first on a harmonic model and then on methane and methanium, where we evaluate all isotope effects of the form CH_4-xD_x/CH₄ and CH_5-xD_x ⁺/CH₅ ⁺, respectively. We discuss the reasons for a surprising behavior of the original method of direct estimators, which performed well for a much larger range of isotope effects than what had been expected previously, as well as some implications of our work for the more general problem of free energy difference calculations.

Direct Computation of the Quantum Partition Function by Path-Integral Nested Sampling

In the present work we introduce a computational approach to the absolute rovibrational quantum partition function using the path-integral formalism of quantum mechanics in combination with the nested sampling technique. The numerical applicability of path-integral nested sampling is demonstrated for small molecules of spectroscopic interest. The computational cost of the method is determined by the evaluation time of a point on the potential energy surface (PES). For efficient PES implementations, the path-integral nested sampling method can be a viable alternative to the direct-Boltzmann-summation technique of variationally computed rovibrational energies, especially for medium-sized molecules and at elevated temperatures.

Effect of Asp122 Mutation on the Hydride Transfer in E. coli DHFR Demonstrates the Goldilocks of Enzyme Flexibility

Dihydrofolate reductase (DHFR) catalyzes the reduction of dihydrofolate (DHF) to tetrahydrofolate (THF) in the presence of NADPH. The key hydride transfer step in the reaction is facilitated by a combination of enzyme active site preorganization and correlated protein motions in the Michaelis-Menten (E:NADPH:DHF) complex. The present theoretical study employs mutagenesis to examine the relation between structural and functional properties of the enzyme. We mutate Asp122 in Escherichia coli DHFR, which is a conserved amino acid in the DHFR family. The consequent effect of the mutation on enzyme catalysis is examined from an energetic, structural and short-time dynamic perspective. Our investigations suggest that the structural and short-time dynamic perturbations caused by Asp122X mutations (X = Asn, Ser, Ala) are along the reaction coordinate and lower the rate of hydride transfer. Importantly, analysis of the correlated and principle component motions in the enzyme suggest that the mutation alters the coupled motions that are present in the wild-type enzyme. In the case of D122N and D122S, the mutations inhibit coupled motion, whereas in the case of D122A, the mutation enhances coupled motion, although all mutations result in similar rate reduction. These results emphasize a Goldilocks principle of enzyme flexibility, that is, enzymes should neither be too rigid nor too flexible.

Path integrals with higher order actions: Application to realistic chemical systems

Quantum thermodynamic parameters can be determined using path integral Monte Carlo (PIMC) simulations. These simulations, however, become computationally demanding as the quantum nature of the system increases, although their efficiency can be improved by using higher order approximations to the thermal density matrix, specifically the action. Here we compare the standard, primitive approximation to the action (PA) and three higher order approximations, the Takahashi-Imada action (TIA), the Suzuki-Chin action (SCA) and the Chin action (CA). The resulting PIMC methods are applied to two realistic potential energy surfaces, for H₂O and HCN-HNC, both of which are spectroscopically accurate and contain three-body interactions. We further numerically optimise, for each potential, the SCA parameter and the two free parameters in the CA, obtaining more significant improvements in efficiency than seen previously in the literature. For both H₂O and HCN-HNC, accounting for all required potential and force evaluations, the optimised CA formalism is approximately twice as efficient as the TIA formalism and approximately an order of magnitude more efficient than the PA. The optimised SCA formalism shows similar efficiency gains to the CA for HCN-HNC but has similar efficiency to the TIA for H₂O at low temperature. In H₂O and HCN-HNC systems, the optimal value of the a₁ CA parameter is approximately 13, corresponding to an equal weighting of all force terms in the thermal density matrix, and similar to previous studies, the optimal α parameter in the SCA was ∼0.31. Importantly, poor choice of parameter significantly degrades the performance of the SCA and CA methods. In particular, for the CA, setting a₁ = 0 is not efficient: the reduction in convergence efficiency is not offset by the lower number of force evaluations. We also find that the harmonic approximation to the CA parameters, whilst providing a fourth order approximation to the action, is not optimal for these realistic potentials: numerical optimisation leads to better approximate cancellation of the fifth order terms, with deviation between the harmonic and numerically optimised parameters more marked in the more quantum H₂O system. This suggests that numerically optimising the CA or SCA parameters, which can be done at high temperature, will be important in fully realising the efficiency gains of these formalisms for realistic potentials.

Quantum and classical simulations of orotidine monophosphate decarboxylase: support for a direct decarboxylation mechanism

Orotidine 5'-monophosphate (OMP) decarboxylase (ODCase) catalyzes the decarboxylation of OMP to uridine 5'-monophosphate (UMP). Numerous studies of this reaction have suggested a plethora of mechanisms including covalent addition, ylide or carbene formation, and concerted or stepwise protonation. Recent experiments and simulations present strong evidence for a direct decarboxylation mechanism, although direct comparison between experiment and theory is still lacking. In the current work we present hybrid quantum mechanics-molecular mechanics simulations that address the detailed decarboxylation mechanisms for OMP and 5-fluoro-OMP by ODCase. Multidimensional potentials of mean force are computed as functions of structural progress coordinates for the Methanobacterium thermoautotrophicum ODCase reaction: the decarboxylation reaction coordinate, an orbital rehybridization coordinate, and the proton transfer coordinate between Lys72 and the substrate. The computed free energy profiles are in accord with the available experimental data. To facilitate further direct comparison with experiment, we compute the kinetic isotope effects (KIEs) for the enzyme-catalyzed reactions using a mass-perturbation-based path-integral method. The computed KIE provide further support for a direct decarboxylation mechanism. In agreement with experiment, the data suggest a role for Lys72 in stabilizing the transition state in the catalysis of OMP and, to a somewhat lesser extent, in 5-fluoro-OMP.

Electrostatically guided dynamics--the root of fidelity in a promiscuous terpene synthase?

Terpene cyclases are responsible for the initial cyclization cascade in the multistep synthesis of more than 60,000 known natural products. This abundance of compounds is generated using a very limited pool of substrates based on linear isoprenoids. The astounding chemodiversity obtained by terpene cyclases suggests a tremendous catalytic challenge to these often promiscuous enzymes. In the current study we present a detailed mechanistic view of the biosynthesis of the monoterpene bornyl diphosphate (BPP) from geranyl diphosphate by BPP synthase using state of the art simulation methods. We identify the bornyl cation as an enzyme-induced bifurcation point on the multidimensional free energy surface, connecting between the product BPP and the side product camphene. Chemical dynamics simulations suggest that the active site diphosphate moiety steers reaction trajectories toward product formation. Nonetheless, chemical dynamics is not precise enough for exclusive product formation, providing a rationale for the lack of fidelity in this promiscuous terpene cyclase.