Nuclear quantum effects on an enzyme-catalyzed reaction with reaction path potential: proton transfer in triosephosphate isomerase

Efficient Surface Hopping Approach for Modeling Charge Transport in Organic Semiconductors

The trajectory surface hopping (TSH) method is nowadays widely applied to study the charge/exciton transport process in organic semiconductors (OSCs). In the present study, we systematically examine the performance of two approximations in the fewest switched surface hopping (FSSH) simulations for charge transport (CT) in several representative OSCs. These approximations include (i) the substitution of the nuclear velocity scaling along the nonadiabatic coupling vector (NCV) by rescaling the hopping probability with the Boltzmann factor (Boltzmann correction (BC)) and (ii) a phenomenological approach to treat the quantum feedback from the electronic system to the nuclear system (implicit charge relaxation (IR)) in the OSCs. We find that charge mobilities computed by FSSH-BC-IR are in very good agreement with the mobilities obtained by standard FSSH simulations with explicit charge relaxation (FSSH-ER), however, at reduced computational cost. A key parameter determining the charge carrier mobility is the reorganization energy, which is sensitively dependent on DFT functionals applied. By employing the IR approximation, the FSSH method allows systematic investigation of the effect of the reorganization energies obtained by different DFT functionals like B3LYP or ωB97XD on CT in OSCs. In comparison to the experiments, FSSH-BC-IR using ωB97XD reorganization energy underestimates mobilities in the low-coupling regime, which may indicate the lack of nuclear quantum effects (e.g., zero point energy (ZPE)) in the simulations. The mobilities obtained by FSSH-BC-IR using the B3LYP reorganization energy agree well with experimental values in 3 orders of magnitude. The accidental agreement may be the consequence of the underestimation of the reorganization energy by the B3LYP functional, which compensates for the neglect of nuclear ZPE in the simulations.

Uncovering the Role of Key Active-Site Side Chains in Catalysis: An Extended Brønsted Relationship for Substrate Deprotonation Catalyzed by Wild-Type and Variants of Triosephosphate Isomerase

We report results of detailed empirical valence bond simulations that model the effect of several amino acid substitutions on the thermodynamic (ΔG°) and kinetic activation (ΔG^⧧) barriers to deprotonation of dihydroxyacetone phosphate (DHAP) and d-glyceraldehyde 3-phosphate (GAP) bound to wild-type triosephosphate isomerase (TIM), as well as to the K12G, E97A, E97D, E97Q, K12G/E97A, I170A, L230A, I170A/L230A, and P166A variants of this enzyme. The EVB simulations model the observed effect of the P166A mutation on protein structure. The E97A, E97Q, and E97D mutations of the conserved E97 side chain result in ≤1.0 kcal mol^-1 decreases in the activation barrier for substrate deprotonation. The agreement between experimental and computed activation barriers is within ±1 kcal mol^-1, with a strong linear correlation between ΔG^⧧ and ΔG° for all 11 variants, with slopes β = 0.73 (R² = 0.994) and β = 0.74 (R² = 0.995) for the deprotonation of DHAP and GAP, respectively. These Brønsted-type correlations show that the amino acid side chains examined in this study function to reduce the standard-state Gibbs free energy of reaction for deprotonation of the weak α-carbonyl carbon acid substrate to form the enediolate phosphate reaction intermediate. TIM utilizes the cationic side chain of K12 to provide direct electrostatic stabilization of the enolate oxyanion, and the nonpolar side chains of P166, I170, and L230 are utilized for the construction of an active-site cavity that provides optimal stabilization of the enediolate phosphate intermediate relative to the carbon acid substrate.

Coarse-graining of many-body path integrals: Theory and numerical approximations

Feynman's imaginary time path integral approach to quantum statistical mechanics provides a theoretical formalism for including nuclear quantum effects (NQEs) in simulation of condensed matter systems. Sinitskiy and Voth [J. Chem. Phys. 143, 094104 (2015)] have presented the coarse-grained path integral (CG-PI) theory, which provides a reductionist coarse-grained representation of the imaginary time path integral based on the quantum-classical isomorphism. In this paper, the many-body generalization of the CG-PI theory is presented. It is shown that the N interacting particles obeying quantum Boltzmann statistics can be represented as a system of N pairs of classical-like pseudoparticles coupled to each other analogous to the pseudoparticle pair of the one-body theory. Moreover, we present a numerical CG-PI (n-CG-PI) method applying a simple approximation to the coupling scheme between the pseudoparticles due to numerical challenges of directly implementing the full many-body CG-PI theory. Structural correlations of two liquid systems are investigated to demonstrate the performance of the n-CG-PI method. Both the many-body CG-PI theory and the n-CG-PI method not only present reductionist views of the many-body quantum Boltzmann statistics but also provide theoretical and numerical insight into how to explicitly incorporate NQEs in the representation of condensed matter systems with minimal additional degrees of freedom.

Carbon Acidity in Enzyme Active Sites

The pK_a values for substrates acting as carbon acids (i.e., C-H deprotonation reactions) in several enzyme active sites are presented. The information needed to calculate them includes the pK_a of the active site acid/base catalyst and the equilibrium constant for the deprotonation step. Carbon acidity is obtained from the relation pK_eq = pKar–pKap = ΔpK_a for a proton transfer reaction. Five enzymatic free energy profiles (FEPs) were calculated to obtain the equilibrium constants for proton transfer from carbon in the active site, and six additional proton transfer equilibrium constants were extracted from data available in the literature, allowing substrate C-H pK_as to be calculated for 11 enzymes. Active site-bound substrate C-H pK_a values range from 5.6 for ketosteroid isomerase to 16 for proline racemase. Compared to values in water, enzymes lower substrate C-H pK_as by up to 23 units, corresponding to 31 kcal/mol of carbanion stabilization energy. Calculation of Marcus intrinsic barriers (ΔG0‡) for pairs of non-enzymatic/enzymatic reactions shows significant reductions in ΔG0‡ for cofactor-independent enzymes, while pyridoxal phosphate dependent enzymes appear to increase ΔG0‡ to a small extent as a consequence of carbanion resonance stabilization. The large increases in carbon acidity found here are central to the large rate enhancements observed in enzymes that catalyze carbon deprotonation.

Neutron scattering in the biological sciences: progress and prospects

The scattering of neutrons can be used to provide information on the structure and dynamics of biological systems on multiple length and time scales. Pursuant to a National Science Foundation-funded workshop in February 2018, recent developments in this field are reviewed here, as well as future prospects that can be expected given recent advances in sources, instrumentation and computational power and methods. Crystallography, solution scattering, dynamics, membranes, labeling and imaging are examined. For the extraction of maximum information, the incorporation of judicious specific deuterium labeling, the integration of several types of experiment, and interpretation using high-performance computer simulation models are often found to be particularly powerful.

The strengths and limitations of effective centroid force models explored by studying isotopic effects in liquid water

The development of effective centroid potentials (ECPs) is explored with both the constrained-centroid and quasi-adiabatic force matching using liquid water as a test system. A trajectory integrated with the ECP is free of statistical noises that would be introduced when the centroid potential is approximated on the fly with a finite number of beads. With the reduced cost of ECP, challenging experimental properties can be studied in the spirit of centroid molecular dynamics. The experimental number density of H2O is 0.38% higher than that of D2O. With the ECP, the H2O number density is predicted to be 0.42% higher, when the dispersion term is not refit. After correction of finite size effects, the diffusion constant of H2O is found to be 21% higher than that of D2O, which is in good agreement with the 29.9% higher diffusivity for H2O observed experimentally. Although the ECP is also able to capture the redshifts of both the OH and OD stretching modes in liquid water, there are a number of properties that a classical simulation with the ECP will not be able to recover. For example, the heat capacities of H2O and D2O are predicted to be almost identical and higher than the experimental values. Such a failure is simply a result of not properly treating quantized vibrational energy levels when the trajectory is propagated with classical mechanics. Several limitations of the ECP based approach without bead population reconstruction are discussed.

Origin of the Non-Arrhenius Behavior of the Rates of Enzymatic Reactions

The origin of the non-Arrhenius behavior of the rate constant for hydride transfer enzymatic reactions has been a puzzling problem since its initial observation. This effect has been used originally to support the idea that enzymes work by dynamical effects and more recently to suggest an entropy funnel model. Our analysis, however, has advanced the idea that the reason for the non-Arrhenius trend reflects the temperature dependence of the rearrangements of the protein polar groups in response to the change in the charge distribution of the reacting system during the transition from the ground state (GS) to the transition state (TS). Here we examine the validity of our early proposal by simulating the catalytic reaction of alcohol dehydrogenase (ADH) and determine the microscopic origin of the entropic and enthalpic contributions to the activation barrier. The corresponding analysis establishes the origin of the non-Arrhenius behaviors and quantifies our original suggestion that the classical effect is due to the entropic contributions of the environment. We also find that the quantum effects reflect in part the temperature dependence of the donor-acceptor distance.

Ab initio path-integral calculations of kinetic and equilibrium isotope effects on base-catalyzed RNA transphosphorylation models

Detailed understandings of the reaction mechanisms of RNA catalysis in various environments can have profound importance for many applications, ranging from the design of new biotechnologies to the unraveling of the evolutionary origin of life. An integral step in the nucleolytic RNA catalysis is self-cleavage of RNA strands by 2'-O-transphosphorylation. Key to elucidating a reaction mechanism is determining the molecular structure and bonding characteristics of transition state. A direct and powerful probe of transition state is measuring isotope effects on biochemical reactions, particularly if we can reproduce isotope effect values from quantum calculations. This article significantly extends the scope of our previous joint experimental and theoretical work in examining isotope effects on enzymatic and nonenzymatic 2'-O-transphosphorylation reaction models that mimic reactions catalyzed by RNA enzymes (ribozymes), and protein enzymes such as ribonuclease A (RNase A). Native reactions are studied, as well as reactions with thio substitutions representing chemical modifications often used in experiments to probe mechanism. Here, we report and compare results from eight levels of electronic-structure calculations for constructing the potential energy surfaces in kinetic and equilibrium isotope effects (KIE and EIE) computations, including a "gold-standard" coupled-cluster level of theory [CCSD(T)]. In addition to the widely used Bigeleisen equation for estimating KIE and EIE values, internuclear anharmonicity and quantum tunneling effects were also computed using our recently developed ab initio path-integral method, that is, automated integration-free path-integral method. The results of this work establish an important set of benchmarks that serve to guide calculations of KIE and EIE for RNA catalysis.

Common enzymological experiments allow free energy profile determination

The determination of a complete set of rate constants [free energy profiles (FEPs)] for a complex kinetic mechanism is challenging. Enzymologists have devised a variety of informative steady-state kinetic experiments (e.g., Michaelis-Menten kinetics, viscosity dependence of kinetic parameters, kinetic isotope effects, etc.) that each provide distinct information regarding a particular kinetic system. A simple method for combining steady-state experiments in a single analysis is presented here, which allows microscopic rate constants and intrinsic kinetic isotope effects to be determined. It is first shown that Michaelis-Menten kinetic parameters (kcat and Km values), kinetic isotope efffets, solvent viscosity effects, and intermediate partitioning measurements are sufficient to define the rate constants for a reversible uni-uni mechanism with an intermediate, EZ, between the ES and EP complexes. Global optimization provides the framework for combining the independent experimental measurements, and the search for rate constants is performed using algorithms implemented in the biochemical software COPASI. This method is applied to the determination of FEPs for both alanine racemase and triosephosphate isomerase. The FEPs obtained from global optimization agree with those in the literature, with important exceptions. The method opens the door to routine and large-scale determination of FEPs for enzymes.

Quantum dynamics in continuum for proton transport--generalized correlation

As a key process of many biological reactions such as biological energy transduction or human sensory systems, proton transport has attracted much research attention in biological, biophysical, and mathematical fields. A quantum dynamics in continuum framework has been proposed to study proton permeation through membrane proteins in our earlier work and the present work focuses on the generalized correlation of protons with their environment. Being complementary to electrostatic potentials, generalized correlations consist of proton-proton, proton-ion, proton-protein, and proton-water interactions. In our approach, protons are treated as quantum particles while other components of generalized correlations are described classically and in different levels of approximations upon simulation feasibility and difficulty. Specifically, the membrane protein is modeled as a group of discrete atoms, while ion densities are approximated by Boltzmann distributions, and water molecules are represented as a dielectric continuum. These proton-environment interactions are formulated as convolutions between number densities of species and their corresponding interaction kernels, in which parameters are obtained from experimental data. In the present formulation, generalized correlations are important components in the total Hamiltonian of protons, and thus is seamlessly embedded in the multiscale/multiphysics total variational model of the system. It takes care of non-electrostatic interactions, including the finite size effect, the geometry confinement induced channel barriers, dehydration and hydrogen bond effects, etc. The variational principle or the Euler-Lagrange equation is utilized to minimize the total energy functional, which includes the total Hamiltonian of protons, and obtain a new version of generalized Laplace-Beltrami equation, generalized Poisson-Boltzmann equation and generalized Kohn-Sham equation. A set of numerical algorithms, such as the matched interface and boundary method, the Dirichlet to Neumann mapping, Gummel iteration, and Krylov space techniques, is employed to improve the accuracy, efficiency, and robustness of model simulations. Finally, comparisons between the present model predictions and experimental data of current-voltage curves, as well as current-concentration curves of the Gramicidin A channel, verify our new model.

RANK-ORDERING THE BINDING AFFINITY FOR FKBP12 AND H1N1 NEURAMINIDASE INHIBITORS IN THE COMBINATION OF A PROTEIN MODEL WITH DENSITY FUNCTIONAL THEORY

The quantum mechanical interaction energies between FKBP12 as well as H1N1 neuraminidase and their inhibitors were directly calculated with an efficient density functional theory by mimicking the whole protein with a protein model composed of the amino acids surrounding the ligands. It was found that the calculated quantum mechanical interaction energies correlate well with the experimental binding free energies with the correlation coefficients of 0.88, 0.86, and the standard deviation of 0.93 and 1.00 kcal/mol, respectively. To compare with force field approach, the binding free energies with the correlation coefficient R = 0.80 and 0.47 were estimated by AutoDock 4.0 programs. It was indicated that the quantum interaction energy shows a better performance in rank-ordering the binding affinity between FKBP12 and H1N1 neuraminidase inhibitors than those of AutoDock 4.0 program. In combination protein model with density functional theory, the estimated quantum interaction energy could be a good predictor or scoring function in structure-based computer-aided drug design. Finally, five new FKBP12 inhibitors were designed based on calculated quantum mechanical interaction energy. In particular, the theoretical K i value of one compound is as low as 0.05 nM, nearly 8-fold more active than FK506.

SUPPLEMENTING THE PBSA APPROACH WITH QUANTUM MECHANICS TO STUDY THE BINDING BETWEEN CDK2 AND N2-SUBSTITUTED O6-CYCLOHEXYLMETHOXYGUANINE INHIBITORS

Because classical Poisson–Boltzman Surface Area (PBSA) model does not allow re-polarization of charges and does not account for charge transfer when a ligand binds to a protein, we have examined a hybrid approach in which we describe the protein–ligand interface by quantum mechanics and the rest of the system with the classical PBSA model. We found this approach to rank order the binding of five N2 -substituted O6 -cyclohexylmethoxyguanine inhibitors to CDK2 (cyclin-dependent kinase 2) properly. The calculated binding free energy correlated well with experimental Log(IC50) with a correlation coefficient of 0.94. A regression fit between experimental Log(IC50) and calculated binding free energy yielded a root-mean-square error of 0.48 when Log(IC50) spanned a range over three units. In addition, we observed charge transfer between the ligand and the protein at the interface — an effect not accounted for by the classical PBSA model. We also found that the direct interactions between the protein and the ligands provided the dominant factor to distinguish the binding affinity of the five ligands studied here. This hybrid approach can better prioritize derivatives of lead compounds for synthesis and biological evaluation.

Coarse-grained (multiscale) simulations in studies of biophysical and chemical systems

Recent years have witnessed an explosion in computational power, leading to attempts to model ever more complex systems. Nevertheless, there remain cases for which the use of brute-force computer simulations is clearly not the solution. In such cases, great benefit can be obtained from the use of physically sound simplifications. The introduction of such coarse graining can be traced back to the early usage of a simplified model in studies of proteins. Since then, the field has progressed tremendously. In this review, we cover both key developments in the field and potential future directions. Additionally, particular emphasis is given to two general approaches, namely the renormalization and reference potential approaches, which allow one to move back and forth between the coarse-grained (CG) and full models, as these approaches provide the foundation for CG modeling of complex systems.

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.

Multiscale modeling of biological functions

Recent years have witnessed a tremendous explosion in computational power, which in turn has resulted in great progress in the complexity of the biological and chemical problems that can be addressed by means of all-atom simulations. Despite this, however, our computational time is not infinite, and in fact many of the key problems of the field were resolved long before the existence of the current levels of computational power. This review will start by presenting a brief historical overview of the use of multiscale simulations in biology, and then present some key developments in the field, highlighting several cases where the use of a physically sound simplification is clearly superior to a brute-force approach. Finally, some potential future directions will be discussed.

Proton-transport mechanisms in cytochrome c oxidase revealed by studies of kinetic isotope effects

Cytochrome c oxidase (CytcO) is a membrane-bound enzyme, which catalyzes the reduction of di-oxygen to water and uses a major part of the free energy released in this reaction to pump protons across the membrane. In the Rhodobacter sphaeroides aa₃ CytcO all protons that are pumped across the membrane, as well as one half of the protons that are used for O₂ reduction, are transferred through one specific intraprotein proton pathway, which holds a highly conserved Glu286 residue. Key questions that need to be addressed in order to understand the function of CytcO at a molecular level are related to the timing of proton transfers from Glu286 to a "pump site" and the catalytic site, respectively. Here, we have investigated the temperature dependencies of the H/D kinetic-isotope effects of intramolecular proton-transfer reactions in the wild-type CytcO as well as in two structural CytcO variants, one in which proton uptake from solution is delayed and one in which proton pumping is uncoupled from O₂ reduction. These processes were studied for two specific reaction steps linked to transmembrane proton pumping, one that involves only proton transfer (peroxy-ferryl, P→F, transition) and one in which the same sequence of proton transfers is also linked to electron transfer to the catalytic site (ferryl-oxidized, F→O, transition). An analysis of these reactions in the framework of theory indicates that that the simpler, P→F reaction is rate-limited by proton transfer from Glu286 to the catalytic site. When the same proton-transfer events are also linked to electron transfer to the catalytic site (F→O), the proton-transfer reactions might well be gated by a protein structural change, which presumably ensures that the proton-pumping stoichiometry is maintained also in the presence of a transmembrane electrochemical gradient. Furthermore, the present study indicates that a careful analysis of the temperature dependence of the isotope effect should help us in gaining mechanistic insights about CytcO.

Approximate inclusion of quantum effects in transition path sampling

We propose a method for incorporating nuclear quantum effects in transition path sampling studies of systems that consist of a few degrees of freedom that must be treated quantum mechanically, while the rest are classical-like. We used the normal mode centroid method to describe the quantum subsystem, which is a method that is not CPU intensive but still reasonably accurate. We applied this mixed centroid/classical transition path sampling method to a model system that has nontrivial quantum behavior, and showed that it can capture the correct quantum dynamical features.

Theoretical analysis of kinetic isotope effects on proton transfer reactions between substituted alpha-methoxystyrenes and substituted acetic acids

Primary kinetic isotope effects (KIEs) on a series of carboxylic acid-catalyzed protonation reactions of aryl-substituted alpha-methoxystyrenes (X-1) to form oxocarbenium ions have been computed using the second-order Kleinert variational perturbation theory (KP2) in the framework of Feynman path integrals (PI) along with the potential energy surface obtained at the B3LYP/6-31+G(d,p) level. Good agreement with the experimental data was obtained, demonstrating that this novel computational approach for computing KIEs of organic reactions is a viable alternative to the traditional method employing Bigeleisen equation and harmonic vibrational frequencies. Although tunneling makes relatively small contributions to the lowering of the free energy barriers for the carboxylic acid catalyzed protonation reaction, it is necessary to include tunneling contributions to obtain quantitative estimates of the KIEs. Consideration of anharmonicity can further improve the calculated KIEs for the protonation of substituted alpha-methoxystyrenes by chloroacetic acid, but for the reactions of the parent and 4-NO(2) substituted alpha-methoxystyrene with substituted carboxylic acids, the correction of anharmonicity overestimates the computed KIEs for strong acid catalysts. In agreement with experimental findings, the largest KIEs are found in nearly ergoneutral reactions, DeltaG(o) approximately 0, where the transition structures are nearly symmetric and the reaction barriers are relatively low. Furthermore, the optimized transition structures are strongly dependent on the free energy for the formation of the carbocation intermediate, that is, the driving force DeltaG(o), along with a good correlation of Hammond shift in the transition state structure.

Calculating solution redox free energies with ab initio quantum mechanical/molecular mechanical minimum free energy path method

A quantum mechanical/molecular mechanical minimum free energy path (QM/MM-MFEP) method was developed to calculate the redox free energies of large systems in solution with greatly enhanced efficiency for conformation sampling. The QM/MM-MFEP method describes the thermodynamics of a system on the potential of mean force surface of the solute degrees of freedom. The molecular dynamics (MD) sampling is only carried out with the QM subsystem fixed. It thus avoids "on-the-fly" QM calculations and thus overcomes the high computational cost in the direct QM/MM MD sampling. In the applications to two metal complexes in aqueous solution, the new QM/MM-MFEP method yielded redox free energies in good agreement with those calculated from the direct QM/MM MD method. Two larger biologically important redox molecules, lumichrome and riboflavin, were further investigated to demonstrate the efficiency of the method. The enhanced efficiency and uncompromised accuracy are especially significant for biochemical systems. The QM/MM-MFEP method thus provides an efficient approach to free energy simulation of complex electron transfer reactions.

A Systematic Approach for Computing Zero-Point Energy, Quantum Partition Function, and Tunneling Effect Based on Kleinert's Variational Perturbation Theory

In this paper, we describe an automated integration-free path-integral (AIF-PI) method, based on Kleinert's variational perturbation (KP) theory, to treat internuclear quantum-statistical effects in molecular systems. We have developed an analytical method to obtain the centroid potential as a function of the variational parameter in the KP theory, which avoids numerical difficulties in path-integral Monte Carlo or molecular dynamics simulations, especially at the limit of zero-temperature. Consequently, the variational calculations using the KP theory can be efficiently carried out beyond the first order, i.e., the Giachetti-Tognetti-Feynman-Kleinert variational approach, for realistic chemical applications. By making use of the approximation of independent instantaneous normal modes (INM), the AIF-PI method can readily be applied to many-body systems. Previously, we have shown that in the INM approximation, the AIF-PI method is accurate for computing the quantum partition function of a water molecule (3 degrees of freedom) and the quantum correction factor for the collinear H(3) reaction rate (2 degrees of freedom). In this work, the accuracy and properties of the KP theory are further investigated by using the first three order perturbations on an asymmetric double-well potential, the bond vibrations of H(2), HF, and HCl represented by the Morse potential, and a proton-transfer barrier modeled by the Eckart potential. The zero-point energy, quantum partition function, and tunneling factor for these systems have been determined and are found to be in excellent agreement with the exact quantum results. Using our new analytical results at the zero-temperature limit, we show that the minimum value of the computed centroid potential in the KP theory is in excellent agreement with the ground state energy (zero-point energy) and the position of the centroid potential minimum is the expectation value of particle position in wave mechanics. The fast convergent property of the KP theory is further examined in comparison with results from the traditional Rayleigh-Ritz variational approach and Rayleigh-Schrödinger perturbation theory in wave mechanics. The present method can be used for thermodynamic and quantum dynamic calculations, including to systematically determine the exact value of zero-point energy and to study kinetic isotope effects for chemical reactions in solution and in enzymes.

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.

Free energies of chemical reactions in solution and in enzymes with ab initio quantum mechanics/molecular mechanics methods

Combined quantum mechanics/molecular mechanics (QM/MM) methods provide an accurate and efficient energetic description of complex chemical and biological systems, leading to significant advances in the understanding of chemical reactions in solution and in enzymes. Here we review progress in QM/MM methodology and applications, focusing on ab initio QM-based approaches. Ab initio QM/MM methods capitalize on the accuracy and reliability of the associated quantum-mechanical approaches, however, at a much higher computational cost compared with semiempirical quantum-mechanical approaches. Thus reaction-path and activation free-energy calculations based on ab initio QM/MM methods encounter unique challenges in simulation timescales and phase-space sampling. This review features recent developments overcoming these challenges and enabling accurate free-energy determination for reaction processes in solution and in enzymes, along with applications.

Origin of the Temperature Dependence of Isotope Effects in Enzymatic Reactions: The Case of Dihydrofolate Reductase

The origin of the temperature dependence of kinetic isotope effects (KIEs) in enzyme reactions is a problem of general interest and a major challenge for computational chemistry. The present work simulates the nuclear quantum mechanical (NQM) effects and the corresponding KIE in dihydrofolate reductase (DHFR) and two of its mutants by using the empirical valence bond (EVB) and the quantum classical path (QCP) centroid path integral approach. Our simulations reproduce the overall observed trend while using a fully microscopic rather than a phenomenological picture and provide an interesting insight. It appears that the KIE increases when the distance between the donor and acceptor increases, in a somewhat counter intuitive way. The temperature dependence of the KIE appears to reflect mainly the temperature dependence of the distance between the donor and acceptor. This trend is also obtained from a simplified vibronic treatment, but as demonstrated here, the vibronic treatment is not valid at short and medium distances, where it is essential to use the path integral or other approaches capable of moving seamlessly from the adiabatic to the diabatic limits. It is pointed out that although the NQM effects do not contribute to catalysis in DHFR, the observed temperature dependence can be used to refine the potential of mean force for the donor and acceptor distance and its change due to distanced mutations.

Ab initio calculation of proton-coupled electron transfer rates using the external-potential representation: A ubiquinol complex in solution

In quantum-mechanical/molecular-mechanical (QM/MM) treatment of chemical reactions in condensed phases, one solves the electronic Schrodinger equation for the solute (or an active site) under the electrostatic field from the environment. This Schrodinger equation depends parametrically on the solute nuclear coordinates R and the external electrostatic potential V. This fact suggests that one may use R and V as natural collective coordinates for describing the entire system, where V plays the role of collective solvent variables. In this paper such an (R,V) representation of the QM/MM canonical ensemble is described, with particular focus on how to treat charge transfer processes in this representation. As an example, the above method is applied to the proton-coupled electron transfer of a ubiquinol analog with phenoxyl radical in acetonitrile solvent. Ab initio free-energy surfaces are calculated as functions of R and V using the reference interaction site model self-consistent field method, the equilibrium points and the minimum free-energy crossing point are located in the (R,V) space, and then the kinetic isotope effects (KIEs) are evaluated approximately. The results suggest that a stiffer proton potential at the transition state may be responsible for unusual KIEs observed experimentally for related systems.

QM/MM Minimum Free Energy Path: Methodology and Application to Triosephosphate Isomerase

Structural and energetic changes are two important characteristic properties of a chemical reaction process. In the condensed phase, studying these two properties is very challenging because of the great computational cost associated with the quantum mechanical calculations and phase space sampling. Although the combined quantum mechanics/molecular mechanics (QM/MM) approach significantly reduces the amount of the quantum mechanical calculations and facilitates the simulation of solution phase and enzyme catalyzed reactions, the required quantum mechanical calculations remain quite expensive and extensive sampling can be achieved routinely only with semiempirical quantum mechanical methods. QM/MM simulations with ab initio QM methods, therefore, are often restricted to narrow regions of the potential energy surface such as the reactant, product and transition state, or the minimum energy path. Such ab initio QM/MM calculations have previously been performed with the QM/MM-Free Energy (QM/MM-FE) method of Zhang et al.1 to generate the free energy profile along the reaction coordinate using free energy perturbation calculations at fixed structures of the QM subsystems. Results obtained with the QM/MM-FE method depend on the determination of the minimum energy reaction path, which is based on local conformations of the protein/solvent environment and can be difficult to obtain in practice. To overcome the difficulties associated with the QM/MM-FE method and to further enhance the sampling of the MM environment conformations, we develop here a new method to determine the QM/MM minimum free energy path (QM/MM-MFEP) for chemical reaction processes in solution and in enzymes. Within the QM/MM framework, we express the free energy of the system as a function of the QM conformation, thus leading to a simplified potential of mean force (PMF) description for the thermodynamics of the system. The free energy difference between two QM conformations is evaluated by the QM/MM free energy perturbation method. The free energy gradients with respect to the QM degrees of freedom are calculated from molecular dynamics simulations at given QM conformations. With the free energy and free energy gradients in hand, we further implement chain-of-conformation optimization algorithms in the search for the reaction path on the free energy surface without specifying a reaction coordinate. This method thus efficiently provides a unique minimum free energy path for solution and enzyme reactions, with structural and energetic properties being determined simultaneously. To further incorporate the dynamic contributions of the QM subsystem into the simulations, we develop the reaction path potential of Lu, et al.2 for the minimum free energy path. The combination of the methods developed here presents a comprehensive and accurate treatment for the simulation of reaction processes in solution and in enzymes with ab initio QM/MM methods. The method has been demonstrated on the first step of the reaction of the enzyme triosephosphate isomerase with good agreement with previous studies.

Fingerprints of Delocalized Transition States in Quantum Dynamics

Reactions with delocalized transition states (plateau reactions) can be characterized statically by their energy profile along the reaction path, where they exhibit a broad, flat region instead of one or several well-defined saddle points on the potential energy surface. Employing our new, highly flexible quantum dynamics code to perform two-dimensional and effective four-dimensional quantum wave packet propagations on ab initio based model potentials, we show that plateau reactions can also be discerned from the other standard reaction types by their dynamics.