9.9 Quantum Mechanical Methods for Enzyme Modeling: Accurate Computation of Kinetic Isotope Effects

9.9 Quantum Mechanical Methods for Enzyme Modeling: Accurate Computation of Kinetic Isotope Effects J Gao, University of Minnesota, Minneapolis, MN, USA r 2012 Elsevier B.V. All rights reserved.

9.9.1 9.9.2 9.9.2.1 9.9.2.2 9.9.2.2.1 9.9.2.2.2 9.9.2.2.3 9.9.2.3 9.9.3 9.9.3.1

Introduction Method Quantum Transition State Theory and Path Integral Simulations Potential Energy Surface Combined QM/MM potentials The X-Pol potential The MOVB method The PI-FEP/UM Method for Accurate Computation of KIEs Illustrative Examples The EH-MOVB Potential Energy Surface for the Hydride Transfer from Trimethylamine to Flavin Cofactor 9.9.3.2 The Decarboxylation of N-Methyl Picolinate in Water 9.9.3.3 Alanine Racemase 9.9.3.4 Nitroalkane Oxidase 9.9.4 Conclusion Acknowledgments References

Abbreviations BLMO CCSDT CHARMM DFT EA-VTST/MT

EH-MOVB FEP GH KIE MC MM

9.9.1

block-localized MO coupled cluster methods chemistry at harvard macromolecular mechanics density functional theory ensemble averaged variational transition state theory with multidimensional tunneling effective Hamiltonian MOVB free energy perturbation Grote-Hynes kinetic isotope effect Michaelis complex molecular mechanics

Introduction

Although the dominant factor in enzyme catalysis is the lowering of the free energy barrier of the enzymatic reaction in comparison with the uncatalyzed process,1–2 quantum mechanical tunneling plays an important role in the calculation of rate constants since proton, hydride, and hydrogen atom transfers are ubiquitous in these processes.3–5 Because of its relatively small mass, zero-point energy and quantum tunneling are significant in determining the free energy reaction barriers for hydrogen transfer reactions.4–8 An intriguing question is whether enzymes have evolved to enhance tunneling to accelerate the reaction rate. However, it is rather difficult to address this question because quantum effects on rate acceleration are much smaller than other

Comprehensive Biophysics, Volume 9

doi:10.1016/B978-0-12-374920-8.00912-7

MO MOVB NAO NDDO NQE PES PI PI-FEP/UM PLP PMF QM QTST VB X-Pol

149 150 150 151 151 152 153 154 155 155 156 156 157 159 159 159

molecular orbital molecular orbital and valence bond nitroalkane oxidase neglect diatomic differential overlap nuclear quantum effect potential energy surface path integral path-integral free energy perturbation and umbrella sampling pyridoxal 50 -phosphate potential of mean force quantum mechanics quantum transition state theory valence bond explicit polarization

factors such as hydrogen bonding and the overall stabilization of the transition state.2,8–9 Nevertheless, a small factor of two in rate enhancement can still have an important physiological impact. The incorporation of nuclear quantum effects (NQEs) in rate calculation is also important for reactions involving heavy atoms since one of the most direct experimental approaches to probe the transition state and the mechanism of an enzymatic reaction is through measurements of kinetic isotope effects (KIE),3,10–11 which is of quantum mechanical origin. Therefore, it is critical to develop computational methods to accurately determine both the primary (a relatively easy problem if hydrogen is involved) and secondary (a more challenging quantity and critical test of computational methods) KIEs for enzymatic reactions.

149

150

Quantum Mechanical Methods for Enzyme Modeling: Accurate Computation of Kinetic Isotope Effects

Each enzyme is unique in that it selectively binds a small number of similar substrates (not the transition state)  sometimes so precise that only a single molecule is recognized by its binding site  along with the remarkable ability to stabilize the transition state, accompanied by protein dynamic conformation change, to accelerate the rate of reaction.12 Electrostatic interactions are undoubtedly the dominant factor in controlling enzyme activity,8,13 but the statement of electrostatic stabilization of transition state provides little insight into the unique catalytic mechanism of a given enzyme. If this were the case, the design of catalytic antibodies would have been much more successful in attaining enzyme activities. Rather, it is the understanding of the specific details that the barrier is lowered and the uniqueness of individual enzymes that are most interesting, and the formulation of general features of catalysis from analyses of individual enzymes that is most useful.8 In the same vein, the notion of contributions from protein dynamics to catalysis cannot be simply characterized by transition-state-theory transmission factor because it is always possible to find a transition state dividing surface, in classical mechanics, such that it is unity.5 Thus, if the question is framed in this fashion, there is no possibility for scientific advance since by definition there is no difference in transmission factor ‘dynamics’ between the catalyzed or uncatalyzed reactions. However, the dynamic fluctuations of different enzymes can be different as demonstrated by numerous experiments,14–19 particularly in the large amplitude motions that are dictated by the specific fold of the enzyme. Thus, the search for a connection between these dynamic motions and the reaction coordinate motion in an enzymatic process is very important.20–21 Quantum mechanics is essential for modeling enzyme mechanism and kinetics, and there are two major challenges in quantum enzymology.5,22 First, electronic structural theory is essential to provide an accurate representation of the potential energy surface (PES) and the description of the transition structure. Depending on the specific biological questions to be addressed, one may use a simple cluster model for the active site, treated by the most sophisticated and accurate quantum chemical models such as density functional theory (DFT) or coupled cluster methods (CCSDT). If dynamics simulations are required, the currently popular approach is to model the active site quantum mechanically (QM) and the rest of the protein and solvent system by molecular mechanics (MM). Such a combined QM/MM approach offers the accuracy of QM that can be systematically improved, and computational efficiency of the force field for large systems. To further increase the capability and reliability of biomolecular computation, it is likely that a full quantum mechanical representation of the entire enzyme system in solution is required. In such a quantal force field, electronic polarization and charge transfer effects are naturally included in the PES.23 The second challenge is to incorporate explicitly NQEs into the time-evolution of macromolecular systems, including zeropoint energy, nonadiabatic coupling and quantum mechanical tunneling. The NQEs are essential in the study of hydrogen transfer processes and the computation of KIEs.5 In this chapter, several novel computational approaches that have been developed in the author’s group for studying enzymatic reactions will be described, including a full

quantum mechanical treatment of solvated proteins and a Feynman path integral simulation technique for accurate computation of both primary and secondary KIEs.

9.9.2 9.9.2.1

Method Quantum Transition State Theory and Path Integral Simulations

The theoretical framework in the present discussion is path integral quantum transition state theory (QTST),24–26 which is derived by writing the rate expression analogous to classical transition state theory. First, in the discrete Feynman path integral (PI) method, each quantized nucleus, atom n, is represented by a ring of P quasiparticles called beads, whose ðnÞ coordinates are denoted as rðnÞ ¼ fri ; i ¼ 1;? ; Pg.27 The ðnÞ ðnÞ discrete paths are circular with rPþ1 ¼ r1 , in which each particle (bead) is connected with its neighbors by a harmonic potential, kðnÞ ¼ 2pP=bl2n , where ln ¼ ð2pb_ 2=Mn Þ1=2 is the de Broglie thermal wavelength of the atom with a mass of Mn. A key concept is the classical correspondence of the atomic coordinates to the centroid variable in path integration,24,26,28–31 defined as the geometrical center of the quasiparticles: rðnÞ ¼

P 1X ðnÞ r P i¼1 i

½1

where the superscript (n) specifies the nth quantized atom. The parameter P is chosen to be sufficiently large such that the numerical result converges to the quantum limit. The QTST rate constant is given by,24–26 ,Z a z 1 a kQTST ¼ /9˙z9Sza ebwðz Þ dz ebwðzÞ ½2 2 N where w(z) is the potential of mean force (PMF) as a function a of the centroid reaction coordinate z[r],  z is the value of z[r] at the maximum of the PMF, and 9˙z9 za ¼ ð2=pbMeff Þ1=2 is a dynamical frequency factor approximated by the velocity for a free particle of effective mass Meff along the reaction coordinate direction. The exact rate constant is obtained by multiplying the QTST rate constant by a correction factor or transmission coefficient gq:25 k ¼ gq  kQTST

½3

Equations [2] and [3] have identical forms to that of the classical rate constant, but unlike transition state theory, there is no variational upper bound in the QTST rate constant because the quantum transmission coefficient gq may be either greater or less than one. There is no practical procedure to compute the quantum transmission coefficient gq. For a model reaction with a parabolic barrier along the reaction coordinate coupled to a bath of harmonic oscillators, the quantum transmission coefficient is the Grote-Hynes (GH) classical transmission coefficient kGH.32 Often, the classical gq is used to approximate the quantum transmission coefficient; however, there is no correspondence between classical and quantum dynamic trajectories and the effects of tunneling may greatly affect reaction dynamics near the barrier top.

Quantum Mechanical Methods for Enzyme Modeling: Accurate Computation of Kinetic Isotope Effects

As in classical transition state theory, the PMF, w(z), can be computed from the equilibrium average: 1 hdðz½r  zÞi DFðzÞ ¼ wðzÞ  wðzR Þ ¼  ln b hdðz½r  zR Þi

½4

where DF(z) is the free energy at z relative to that at the reactant state minimum zR, and the ensemble average /? S is obtained by a quantum mechanical effective potential:27 Veff

hn

ðnÞ

ri

N P o i X pP X ðnÞ ðnÞ ;S ¼ ðri  riþ1 Þ2 2 n¼1 bln i

þ

P 1X ð1Þ ðNÞ Uðri ; ?; ri ; SÞ P i

½5

with N being the number of quantized atoms, S representing the coordinates of all classical (solvent) particles, and ð1Þ ðNÞ Uðri ; ?; ri ; SÞ the potential energy of the system. Note that the dynamics generated by the effective potential of eqn [5] has no physical significance; it is used as a procedure to obtain the correct ensemble of configurations.27,33 Studies have shown that the activation free energy in the centroid PI-QTST ‘‘captures most of the tunneling and quantization effects’’; which give rise to deviations from classical TST.34–37 In centroid path integral, the centroid positions, r, are used as the principle variable and the canonical QM partition function of a hybrid quantum and classical system, with one quantized particle for convenience of discussion, can be written as follows:27 qm

QP ¼

!3P=2 Z Z Z 1 P dS dr 2 dRebVeff ðfrg;SÞ O lM

½6

where O is the volume element of classical particles, P is the number of quasiparticles, Veff ðfrg; SÞ is the effective potential (eqn [5]), and, importantly, eqn [6] can be rewritten exactly as a double average,27,38 which is the theoretical basis in the simulation approach of Sprik et al., called the hybrid classical and path integral,39 and of Hwang and Warshel, called QCP.40 Qqm ¼ Qcm //ebDUðr;SÞ SFP;r SU

½7

where the factor Qcm is the classical partition function defined in ref. 27, and the average h?iU is a purely classical ensemble average obtained according the potential U(r,S). In eqn [7], the average differential potential is given by DUðr;SÞ ¼

P 1X fUðri ; SÞ  Uðr;SÞg P i

½8

and the inner average /?SFP;r represents a path-integral freeparticle sampling, carried out without the external potential U(r,S).27,39–40 This double averaging procedure39–41 yields the exact path integral centroid density, which can be used to determine the centroid PMF:42 e

bwðzÞ

¼e

bwcm ðzÞ

/dðz ¼ zÞ/e

bDUðz½r;SÞ

SFP;z SU

151

To enhance convergence of free-particle sampling in centroid path integral simulations, Dan Major developed a bisection sampling technique for a ring of beads, by extending the original approach of Ceperley for free particle sampling in which the initial and final beads are not connected.43–44 Since the free-particle distribution is known exactly at a given temperature, each ring-beads distribution is generated according to this distribution and thus 100% accepted.44 Furthermore, in this construction, each new configuration is created independently, starting from a single initial bead position, allowing the new configuration to move into a completely different region of configurational space. This latter point is especially important in achieving convergence by avoiding being trapped in a local region of the classical potential in the path integral sampling. This method has been thoroughly tested in a number of applications.42,45–48

9.9.2.2

Potential Energy Surface

The potential energy function describes the energetic change as a function of the variations in atomic coordinates, including thermal fluctuations and rearrangements of the chemical bonds. The accuracy of the potential energy function used to carry out molecular dynamics simulations directly affects the reliability of the computed wðzÞ.8,22,49–50 The use of combined QM/MM methods for enzyme reactions will be discussed, and then, the idea is extended to a mixed molecular orbital and valence bond (MOVB) theory to represent the active site.51–53 A quantal explicit polarization (X-Pol) potential in which the entire protein–solvent system is treated by electronic structural methods is also discussed.23,54–57

9.9.2.2.1

Combined QM/MM potentials

The current popular approach to model the PES of enzymatic reactions is combined quantum mechanical and molecular mechanical (QM/MM) methods, in which the system is partitioned into a QM region, typically consisting of the substrate and residues directly participating in the chemical process, and an MM region, including the rest of the system.49,58–59 The QM/MM potential is given by49 o Utot ¼ hCðSÞjHqm ðSÞ þ Hqm=mm ðSÞjCðSÞi þ Umm

½10

o ðSÞ is the Hamiltonian of the QM-subsystem, Umm where Hqm is the classical (MM) potential energy for the MM region, and Hqm/mm(S) is the interaction Hamiltonian between the two regions. The wave function C(S) of the QM-subsystem is optimized to minimize the energy of the electronic Hamiltonian o ðSÞ þ Hqm=mm ðSÞ. Hqm It is most convenient to rewrite eqn [10] as follows:58

Utot ¼ Eoqm ðSÞ þ DEqm=mm ðSÞ þ Umm

½11

where Eoqm ðSÞ is the energy of an isolated QM-subsystem, o Eoqm ðSÞ ¼ hCo ðSÞjHqm ðSÞjCo ðSÞi

½12

½9

where wðzÞ and wcm(z) are the centroid quantum mechanical and the classical mechanical PMF.

and DEqm/mm(S) is the interaction energy corresponding to the energy change of transferring the QM subsystem from the gas phase into the condensed-phase (enzyme active site).

152

Quantum Mechanical Methods for Enzyme Modeling: Accurate Computation of Kinetic Isotope Effects X-Pol wave function:65–66

DEqm/mm(S) is defined by o DEqm=mm ðSÞ ¼ hCðSÞjHqm ðSÞ þ Hqm=mm ðSÞjCðSÞi  Eog ðSÞ ½13

C

XPolX

¼

^ RAx A

M Y

! Fa

½17

a¼1

The QM/MM PES combines the generality of quantum mechanical methods for treating chemical processes with the computational efficiency of molecular mechanics for large molecular systems.60 This is important because the dynamic fluctuations of the enzyme and aqueous system have a major impact on the polarization of the species involved in the chemical reaction, which, in turn, affects the chemical reactivity.58–59,61

9.9.2.2.2

The X-Pol potential

The current force fields for biomolecular systems were established by Lifson in the 1960s,62 which have changed very little in the past half a century. Despite its success in biomacromolecular modeling, there are also many shortcomings, including redundancy of empirical parameters and a lack of unified treatment of electronic polarization. A logical step in the future is to go beyond the molecular mechanical representation of biomolecules. To this end, we have introduced the X-Pol method based on quantum mechanics as a framework for the development of next-generation force fields.23,54–56,63–64 The X-Pol potential is designed with a hierarchy of three approximations.64 The macromolecular system is partitioned into subsystems, which are called fragments or blocks. The first approximation is that the wave function of the full system (CX-Pol) is a Hartree product of the determinant wave functions of individual fragments: CXPol ¼ Rx

M Y

^ a AF

½14

a¼1

where M is the number of blocks or residue fragments, Rx is the normalization constant, Aˆ is an antisymmetrization operator, and Fa is a product of the occupied spin-orbitals fjai g in block a (eqn [14]): Fa ¼ ja1 ja2    jana

½15

The assumption made in eqn [14] neglects the exchange interactions between different fragments. To account for the short-range exchange repulsion as well as the long-range dispersion interactions, the Lennard-Jones potential is used EvdW ab

¼

A X B X a¼1 b¼1

" 4eab

sab Rab

12



sab  Rab

6 # ½16

where A and B are the number of atoms in fragments a and b, and the parameters eab and sab are obtained using standard combining rules such that eab ¼ (eaeb)1/2 and sab ¼ (sa þ sb)/2, in which e and s are atomic empirical parameters. The use of the R12 repulsive terms in eqn [16] to approximate the short-range exchange repulsion significantly reduces the computational cost. Alternatively, the exchange energies can be determined explicitly by antisymmetrizing the

Here, the notation X-Pol-X is used to indicate that the X-Pol potential includes eXchange explicitly. The wave functions of the individual fragments are optimized by the self-consistent field method in the presence of the external electrostatic potential of all other blocks until the energy or electron density of the entire system is converged.54,63 Thus, for fragment a the external potential, Va(r), is ( Z ) B X dr0 rb ðr0 Þ X Zb ðbÞ  Va ðrÞ ¼ þ ½18 9r  r0 9 9r  Rb ðbÞ9 b¼1 ba a where rb(r0 ) is the electron density of molecule b, derived from the wave function, rb(r0 ) ¼ 9Cb(r0 )92. The total potential energy of the system is N   X ^ Etot ¼ F9H9F  Eoa

½19

a¼1

where Hˆ is the Hamiltonian of the system defined below, Eoa is the energy of an isolated fragment a whose wave function is Coa . The individual wave functions of the subsystems can be obtained at any level of theory – ab initio Hartree-Fock, semiempirical molecular orbital theory, correlated wave function theory, or Kohn-Sham DFT.67 Without further approximation, it is necessary to compute the two-electron integrals arising from different molecules, which would be too expensive for a force field designed for condensed phase simulations. Consequently, we introduce the second approximation in the X-Pol method by approximating the ‘‘quantum mechanical’’ electrostatic potential of eqn [18] by a classical multipole representation. We have used the monopole term, that is, partial atomic charges, in early applications, in which the formally two-electron integrals are reduced to one-electron integrals, which are computationally efficient.54–55,64 The third approximation is the specific quantum mechanical model to be used for a given problem and a specific purpose. To this end, we adopt semiempirical Hamiltonians based on the neglect diatomic differential overlap (NDDO) approximation.68–69 The X-Pol potential has been tested and applied to the simulation of liquid water,55 and has it been extended to molecular dynamics simulations of polypeptide in solution.23 Analytic gradient technique has been developed and implemented into the program CHARMM (Chemistry at Harvard Macromolecular Mechanics),70 and the X-Pol method has been illustrated to be feasible for extensive molecular dynamics simulations of fully solvated proteins under periodic boundary conditions.23 In that study, 3,200 (3.2 ps) full energy and gradient evaluations can be performed in molecular dynamics simulations of a BPTI protein in water, consisting of nearly 15 000 atoms and 30 000 basis functions, on a single 2.6 GHz processor in 1 day.

Quantum Mechanical Methods for Enzyme Modeling: Accurate Computation of Kinetic Isotope Effects 9.9.2.2.3

The MOVB method

Often, it is desirable to specifically define the electronic state corresponding to the reactant and product configurations, in context of valence bond (VB) theory, to model the PES for chemical and enzymatic reactions. In principle, this provides a procedure to estimate the solvent and protein reorganization energies along the diabatic state coordinates for the reaction and product configurations. This can be conveniently done using empirical potential functions.71 However, it is difficult to use delocalized molecular orbital and density functional theory to model charge localized diabatic states. To this end, we have developed a mixed molecular orbital and valence bond (MOVB) theory, in which the effective diabatic states are constructed by the X-Pol-X method (also known as block-localized wave function).51–53 Here, molecular orbitals (MOs) are strictly block-localized within individual fragments of a molecular system based on the Lewis structures of the reactant or product states. The potential energy surface of the adiabatic ground state is obtained by avoided crossing that couples the two or more diabatic states. Key features of the MOVB model include (1) that the MOs within each fragment are orthogonal, which makes computation efficient; and (2) that the block-localized MOs (BLMOs) between different fragments are nonorthogonal, which retains important characteristics of VB theory. In the following we use the hydride transfer reaction between trimethylamine, (CH3)3N (TMA-H), and a flavin cofactor (Nf þ ) model (hereafter simply called flavin) to illustrate the MOVB method (Scheme 1). First, we define the reactant diabatic state (Scheme 1), Cr(R), as follows: ^ TMA-H wNf g Cr ¼ Afw r r þ

þ -Nf denote the products of occupied and wH where wTMA p p BLMOs expanded over basis orbitals in fragments TMA þ and H-Nf, respectively. The MOVB wave function for the reactive system is written as a linear combination of the diabatic states.

Fg ðRÞ ¼ ar Cr ðRÞ þ ap Cp ðRÞ

þ

^ TMA wH-Nf g Cp ðRÞ ¼ Afw p p þ

½21

CH3

H

N

N N

N O

½22

where ar and ap are the configurational coefficients for the reactant and product diabatic state, respectively.53,66,72 The potential energy of the adiabatic ground state, Vg(R), can be obtained either by consistently optimizing both the configuration coefficient and the orbital coefficients in the diabatic configurations (called CDC),53 or by using the lower energy root of the equation:    Hrr ðRÞ  VðRÞ Hrp ðRÞ  Srp ðRÞVðRÞ   ½23 ¼0    Hpr ðRÞ  Spr ðRÞVðRÞ Hpp ðRÞ  VðRÞ where Hrr(R) and Hpp(R) are the variationally optimized Hamiltonian matrix elements for the reactant and product diabatic states (called VDC),53 respectively, Hrp(R) ¼ Hpr(R) is the exchange integral (off-diagonal matrix element), and Srp(R) ¼ Spr(R) is the overlap integral between the two diabatic states. The Hamiltonian matrix elements in eqn [23] are given as follows:52,67  h i i 1 h Hab ¼ Sab Tr ðDab ÞT h þ Tr ðDab ÞT JDab 2 i 1 h  Tr ðDab ÞT KDab þ Enuc 4

½20

where wrTMA-H and wNf specify the products of occupied r BLMOs in fragments TMA-H and Nf þ , respectively (Scheme 1). Similarly, the wave function of the product state (Scheme 1), Cp(R), is expressed as

153

 ½24

where the subscripts a and b specify either the reactant (r) or the product (p) state or both, Enuc is the nuclear Coulomb energy, Sab and Dab are the overlap integral and density matrix over nonorthogonal determinant wave functions, and h, J, and K are the standard one-electron, Coulomb and exchange matrices. Equation [24] is a general formula, valid both for molecular orbital theory and for Kohn-Sham DFT.67 In the

CH3

H

O

N

N

H

N

O N

H

Nf+

O

H

H-Nf

H CH3

N

CH2

CH3

(CH3)2N

CH2 TMA+

TMA-H

Scheme 1 Schematic illustration of the definition of the reactant (left) and product (right) diabatic states in the hydride transfer reaction between trimethylamine (TMA-H) and a model for the flavin cofactor (Nf þ ). Cembran, A.; Payak, A.; Lin, Y. L.; Xie, W.; Song, L.; Mo, Y.; Gao, J. A non-orthogonal block-localized effective hamiltonian approach for chemical and enzymatic reactions. J. Chem. Theory Comput. 2010, 6, 2242. Copyright by American Chemical Society.

154

Quantum Mechanical Methods for Enzyme Modeling: Accurate Computation of Kinetic Isotope Effects

latter case, the exchange integral K is replaced by the exchangecorrelation potential.67 An attractive MOVB approach is to employ the semiempirical models based on the NDDO approximation,68–69 and this approach, called effective Hamiltonian MOVB (EHMOVB), has been implemented in the Chemistry at Harvard Macromolecular Mechanics (CHARMM) package.73 Typically, semiempirical models, such as AM1,74 PM375 or the selfconsistent charge tight-bonding density functional algorithm,76 can adequately describe the qualitative features of the potential surface. However, the quantitative results are often not satisfactory. If we treat the MOVB matrix elements as effective Hamiltonian,72 then, the quantitative errors in semiempirical methods can be eliminated in a way similar to that in empirical71,77 or semiempirical VB models.78–80 Specifically, we introduce a parameter in the off-diagonal Hamiltonian matrix element Hrp, aimed to reproduce the barrier height for a given chemical reaction: EH ¼ Hrp þ grp Hrp

½25

where grp is a parameter that dominantly affects the computed EH is the total effective Hamiltonian barrier height, and Hrp resonance integral. Another formalism is to scale the off-diagonal matrix element as follows:72 EH Hrp

¼ zrp Hrp

½27

where the parameter De ¼ DEexpt  Hpp(Rp) þ Hrr(Rr). The procedure outlined above is identical to that used in the parameter ‘calibration’ of empirical valence bond models, such as that in refs77 and71, or of semiempirical valence bond,78–79 which allows the barrier height and reaction energy to be readily fitted to their target values exactly. In general, however, it is more challenging to ‘calibrate’ the change of molecular structure along the entire reaction path, especially the precise geometry of the transition state.81 Since the Hamiltonian matrix elements are computed explicitly using electronic structural methods and are dependent on all degrees of freedom of the system, the PES can be adequately modeled using the MOVB method.53,82

9.9.2.3

" #" # R QLqm ðzLa Þ QH kL qm ðzH Þ bfFR ðzRL ÞFR ðzRH Þg L H KIE ¼ H ¼ e a k QH QLqm ðzRL Þ qm ðzH Þ

½28

½26

Both options can be useful, depending on the performance of the semiempirical model and the specific reaction considered. For the hydride transfer reaction between trimethylamine and flavin, eqn [25] was employed. The second parameter of the EH-MOVB model is the adjustment of the relative energy between the reactant and product diabatic states by shifting the product state energy, Hpp, by an amount of De to yield the desired reaction energy: EH ¼ Hpp þ De Hpp

functions for two different isotopes through free energy perturbation (FEP) theory. In essence, only one simulation of a given isotopic reaction, for example, the light isotope, is performed, while the ratio of the partition function, that is, the KIE, to a different isotopic reaction, is obtained through FEP by perturbing the mass from the light isotope to the heavy isotope.42 This is in contrast to other approaches that have been reported in the literature employing centroid path integral simulations,40,83–84 in which two separate simulations are performed. The use of ‘mass’ perturbation in free-particle bisection sampling schemes results in a major improvement in computation accuracy for KIE calculations such that secondary kinetic isotope effects and heavy atom isotope effects can be reliably obtained.42,85 To our knowledge, the coupled pathintegral free energy perturbation and umbrella sampling (PI-FEP/UM) method is the only practical approach to yield computed secondary KIEs sufficiently accurate to be compared with experiments.86 Specifically, the light-to-heavy kinetic isotope effect is rewritten in terms of the ratio of the partial partition functions at the centroid reactant and transition state and is given by:42

The PI-FEP/UM Method for Accurate Computation of KIEs

The centroid path integral method described in section 9.9.2.1 enables us to conveniently determine kinetic isotope effects by directly computing the ratio of the quantum partition

The ratio of the partition function in eqn [28] can be written as follows: QH qm ðzÞ QLqm ðzÞ

b

¼

/dðz  zÞ/eP

P i

DUiL-H bDU L

e

SFP;L SL

o

/dðz  zÞeb½FL ðz;SÞFFP  SL

½29

where the subscript L specifies that the ensemble averages are o is done using the light isotope, DU L is defined by eqn [8], FFP the free energy of the free particle reference state for the quantized particles,27 and DUiL-H ¼ Uðri;H Þ  Uðri;L Þ represents the difference in ‘‘classical’’ potential energy at the heavy and light bead positions ri,H and ri,L. In the bisection sampling scheme,45 the perturbed heavy isotope positions are related to the lighter ones by ri;L lM hi ¼ L ¼ ri;H lMH hi

rffiffiffiffiffiffiffiffi MH ; ML

i ¼ 1; 2; ?; P

½30

where ri,L and ri,H are the coordinates for bead i of the corresponding light and heavy isotopes, lML and lMH are the de Broglie wave lengths for the light and heavy nuclei, and hi is the position vector in the bisection sampling scheme which depends on the previous sequence of directions and has been fully described in reference45. In eqn [29], we obtain the free energy (inner average) difference between the heavy and light isotopes by carrying out the bisection path integral sampling with the light atom and then perturbing the heavy isotope positions according to eqn [30]. Then, the free energy difference between the light and heavy isotope ensembles is weighted by a Boltzmann factor for each quantized configuration.42

Quantum Mechanical Methods for Enzyme Modeling: Accurate Computation of Kinetic Isotope Effects

9.9.3

Illustrative Examples

9.9.3.1

The EH-MOVB Potential Energy Surface for the Hydride Transfer from Trimethylamine to Flavin Cofactor

The adiabatic ground state potential energy surfaces determined using DFT with the B3LYP and M06-2X functionals are compared with the standard semiempirical AM1 model and the EHMOVB method in Figure 1 as a function of the reaction coordinate Rc, defined as the difference between the bond distances for the breaking (C–H) and forming (H–N) bonds, for the hydride transfer reaction between trimethylamine and a flavin cofactor (see Scheme 1), which is used as a model reaction for the histone lysine specific demethylase one (LSD1) enzymatic process. The 6-31 þ G(d) basis set was used in DFT calculations. The M06-2X density functional method is currently one of the most accurate DFT methods, which yields an estimated barrier height of 17.4 kcal mol1 and a relative energy of –6.1 kcal mol1 between the product and reactant state. The popular hybrid B3LYP method slightly underestimates the hydride transfer barrier at 15.1 kcal mol1. The semiempirical AM1 energy profile is qualitatively correct, but it contains two main problems; the computed energy of activation is 22.6 kcal mol1, about 5 kcal mol1 too high, and the predicted energy of reaction is too endothermic by 6.4 kcal mol1. Using EH-MOVB,82 we obtained an activation energy of 18.1 kcal mol1 for the hydride transfer between trimethylamine and flavin and an energy of reaction of –4.4 kcal mol1. The minimum energy path for the hydride transfer from trimethylamine to flavin has been optimized as a function of the

155

reaction coordinate Rc defined above. In the present study, we have constrained the hydride migration to be collinear with the donor (C) and acceptor (N) atoms, whereas all other degrees of freedom are fully minimized using the ABNR algorithm in CHARMM.73 The potential energy curves for the reactant and product diabatic states, as well as that of the adiabatic ground state, are shown in Figure 2 as a function of the diabatic energy difference, or the energy-gap reaction coordinate: DE ¼ Hrr  Hpp

½31

Here, the energy gap is projected from the energies of the diabatic states along the minimum energy path. The main objective of Figure 2 is to illustrate that the minimum energy potential surface for the adiabatic ground state and those for the diabatic states can be fully represented with the use of either a geometrical or an energy-gap reaction coordinate. The reactant state potential shows a steady increase as the reaction coordinate changes from the reactant to the product side, which is mirrored by the product state energy in the reverse direction. The trend of the two diabatic potential energy curves is consistent with heterolytic bond cleavages of the reactant (C–H) and the product (N–H) species. At the diabatic state crossing point, which corresponds roughly to the location of the transition state of the hydride transfer reaction, the diabatic state is circa 40 kcal mol1 in energy above the adiabatic ground state, suggesting that there is significant electronic coupling between the reactant and the product states. The coupling energy is similar to values determined for proton transfer and nucleophilic substitution reactions using ab initio WFT and DFT.51–53,67,72,87

25 EH-MOVB (AM1)

175

AM1

Vg

B3LYP

Hrr

M06-2x

Energy (kcal mol−1)

ΔE (kcal mol−1)

15

5

−5

125 Hpp

75

25

−15 −2.0

−1.0

0.0

1.0

2.0

Rc (Å) Figure 1 Computed potential energy profile along the minimum energy path (Rc ¼ R[C  H]  R[H  N]) for the hydride transfer reaction between trimethylamine and flavin model using EHMOVB(AM1) in red, AM1 in light blue, B3LYP/6-31G(d) in navy blue, and M06-2X/6-31G(d) in green. Cembran, A.; Payak, A.; Lin, Y. L.; Xie, W.; Song, L.; Mo, Y.; Gao, J. A non-orthogonal block-localized effective hamiltonian approach for chemical and enzymatic reactions. J. Chem. Theory Comput. 2010, 6, 2242. Copyright by American Chemical Society.

−25 −160

−80

0 ΔE (kcal

80

160

mol−1)

Figure 2 Potential energy profiles for the reactant state (blue), the product state (green), and the adiabatic ground state (red) as a function of the energy difference between the reactant and product diabatic states (i.e., the energy-gap reaction coordinate). Cembran, A.; Payak, A.; Lin, Y. L.; Xie, W.; Song, L.; Mo, Y.; Gao, J. A nonorthogonal block-localized effective hamiltonian approach for chemical and enzymatic reactions. J. Chem. Theory Comput. 2010, 6, 2242. Copyright by American Chemical Society.

156

Quantum Mechanical Methods for Enzyme Modeling: Accurate Computation of Kinetic Isotope Effects

C2

O

C2

N

C

N

CH3

O

CH3

+ CO2

Scheme 2 Decarboxylation of N-methyl Picolinate.

The Decarboxylation of N-Methyl Picolinate in Water

The decarboxylation of N-methyl picolinate (Scheme 2) was used as a model to probe the mechanism of the uncatalyzed reaction in water for orotidine 50 -monophosphate decarboxylase. The primary and secondary heavy atom kinetic isotope effects were determined by Rishavy and Cleland, extrapolated to 25 1C.88 Classical molecular dynamics simulations and umbrella sampling were first carried out for a system containing the N-methyl picolinate, treated by the AM1 Hamiltonian, in a cubic box (ca. 30  30  30 A˚3) of 888 water molecules, described by the TIP3P potential.89 In all calculations, long-range electrostatic interactions were treated by the particle mesh-Ewald method for QM/MM potentials.90 The classical mechanical PMF as a function of the C2-CO2 distance was obtained at 25 1C and 1 atm through a total of 2 ns of dynamics simulations. Then, the coordinates from the classical trajectory were used in PI-FEP/UM simulations, employing a total of 97 800 classical configurations for each isotope (12C, and 14N) for the decarboxylation reaction. For each classical configuration, 10 steps of path-integral freeparticle sampling were performed. Each quantized atom was represented by 32 quasiparticles. Solvent effects are significant, increasing the free energy barrier by 15.2 kcal mol1 to a value of 26.8 kcal mol1, which is accompanied by a net free energy of reaction of 24.7 kcal mol1. The large solvent effect is due to the presence of a positive charge on the pyridine nitrogen, which is annihilated in the decarboxylation reaction. Both the 12C/13C primary KIE and the 14N/15N secondary KIE have been determined, with the immediate adjacent atoms about the isotopic substitution site also quantized. Figure 3 depicts the ratio of the 12C/13C partition functions along the C2–CO2 distance to illustrate the computational sensitivity. First, we note that the NQEs are non-negligible even for bond cleavage involving two carbon atoms, which reduce the free energy barrier by 0.45 kcal mol1. The computed intrinsic 13C primary KIE is 1.0318 7 0.0028 at 25 1C for the decarboxylation of N-methyl picolinate in water. For comparison, the experimental value is 1.0281 7 0.0003 at 25 1C. For the secondary 15 N KIE, the PI-FEP/UM simulation yields an average value of 1.008370.0016, which may be compared with experiment (1.007070.0003).88 The agreement between theory and experiment is excellent, which provides support for a unimolecular decarboxylation mechanism in this model reaction.

9.9.3.3

Alanine Racemase

Alanine racemase catalyzes the interconversion of L- and D-alanine, the latter of which is an essential component in the peptidoglycan layer of the bacterial cell wall. The chemical

0.85

Q qm(12C)/Q qm(13C)

9.9.3.2

0.80

0.75 1

2

3

z(C2−CO2) (Å) Figure 3 Computed ratio of the quantum mechanical partition functions between 12C and 13C isotopic substitutions at the carboxyl carbon position for the decarboxylation reaction of N-methyl picolinate in aqueous solution from the PI-FEP/UM simulations.

transformation is illustrated in Scheme 3, which has been modeled by a combined QM/MM potential in molecular dynamics simulations. Employing the PI-FEP/UM method, Dan Major computed the primary KIEs both for the forward and reverse processes.48,91 In these calculations, the semiempirical AM1 formalism was used to describe the active site, which includes the pyridoxal 50 -phosphate (PLP) cofactor bound substrate and the acid and base residues Lys39 and Tyr2650 , where the prime indicates a residue from the second subunit of the dimeric enzyme. In this study, the AM1 model was reparametrized to model the AlaR-catalyzed racemization, yielding a highly accurate Hamiltonian that is comparable to mPW1PW91/ 6-311 þ þ G(3df, 2p) results. Stochastic boundary molecular dynamics simulations were carried out for a system of 30 A˚ sphere about the center of the active site, and a series of umbrella sampling simulations were performed, cumulating a total of 24 ns statistical sampling (with 1 fs integration step) to yield the classical potential of mean force (the outer average in eqn [7]). Then, the transferring proton, the donor and acceptor heavy atoms for each process are quantized by a centroid path integral with 32 beads for each particle. About 15 000 configurations saved from these trajectories in regions corresponding to the Michaelis complex reactant state, transition state, and product state were extracted, each of which was subjected to 10 bisection free-particle sampling to yield the centroid PI quantum corrections to the classical free energy profile. The quantum mechanical PMF are displayed in Figure 4, which incorporate the quantum mechanical corrections to the classical PMF from molecular dynamics simulations. For the

Quantum Mechanical Methods for Enzyme Modeling: Accurate Computation of Kinetic Isotope Effects

Tyr265′ O

H3N Lys39

O2C L-Ala + E-PLP

H O H3C

H

H3N Lys39

Tyr265′ OH

CH3

O2C k1

N OPO32−

H O

k−1

N

O2C k2

N

H3C

Arg219

H2N Lys39

Tyr265′ OH

CH3

OPO32− N

H O

k−2

157

H

CH3

N OPO32−

H3C

D-Ala + E-PLP

N Arg219

Arg219

Scheme 3 Forward and reverse proton transfer reactions in the active site of alanine racemase. Major, D. T.; Nam, K.; Gao, J. Transition state stabilization and a-amino carbon acidity in alanine racemase. J. Am. Chem. Soc. 2006, 128, 8114–8115. Copyright by American Chemical Society.

going from the reactant state to the transition state, and that hydrogen tunneling is negligible with an average transmission factor of 1.14 and 1.31 for the hydrogen and deuterium. Furthermore, the net quantum effects, from EA-VTST/MT calculations, lower the classical barrier by 2.7 and 1.9 kcal mol1 for the hydrogen and deuterium transfer reactions, in accord with the PI-FEP/UM calculations. For the D-L conversion of alanine by AlaR, the actual proton transfer step is not rate limiting. If we neglect the complexity of reaction steps involving internal and external aldimine exchange, substrate binding and product release, that is, we only consider the proton abstraction and reprotonation steps, the overall rate constant for the ‘chemical step’ can be expressed as follows:48

20 18.7

cl

Free energy (kcal mol−1)

17.0 15

15.4

qm (H)

16.1

qm (D) 13.2

10

12.1 6.6

5

0

−2.1 −2.8

−2.7 −5 −4

−3

−2

−1

0

1

2

3

4

keff ¼

Reaction coordinate (Å)

Figure 4 Classical and quantum potentials of mean force for the proton and deuteron abstraction of Ala-PLP by Tyr2650 , and the reprotonation of Ala-PLP carbanion intermediate by Lys39 in the active site of alanine racemase.

L-D alanine racemization in AlaR, the first proton abstraction step by Tyr2650 is rate-limiting, and thus the observed rate constant is directly related to this reaction step, kobs ¼ k1, and the KIE is computed using this rate constant. Figure 4 shows that inclusion of NQEs to the computed classical PMF, lowers the free energy barrier by 2.6 and 1.7 kcal/mol for proton and deuteron transfer to Tyr2650 phenolate ion in AlaR. This leads to a computed intrinsic KIE of 4.2 for the a-proton abstraction in the L-D alanine conversion. However, the computational result is much greater than the experimental value of 1.9.92 The origin of this discrepancy is still unclear, which may be a result of one or a combination of factors, including the complexity involved in the analysis of experimental kinetic data, the possibility that the experimental KIE was not exactly the intrinsic value, and computational uncertainty. For comparison, the primary intrinsic KIE was also computed using the ensemble averaged variational transition state theory with multidimensional tunneling, or EA-VTST/MT, method,5,93 and a value of 4.0 for the Tyr2650 proton abstraction was obtained, in close agreement with the PI-FEP/ UM results. Interestingly, the EA-VTST/MT method allows the separation of the total NQEs into zero-point energies and tunneling. We found that the dominant quantum mechanical contribution is due to the change in zero-point energy in

k2 k1 k1 þ k2

½32

Since our simulations (Figure 4) show that k2ck1, we obtain the following rate expression: keff ¼

k2 k1 k1 ¼ k2 K2

½33

where K2 is the equilibrium constant for the proton transfer reaction from L-Ala-PLP to a neutral Lys39. Based on the relative free energies in Figure 4, we obtained an estimated primary KIE of 3.1 from path integral simulations, and 2.6 from EA-VTST/MT. Again, these values are greater than that from a multicomponent analysis of experimental kinetic data (1.4370.20).92 The difference between computational results and values from analysis of experimental data suggest that further studies are needed.

9.9.3.4

Nitroalkane Oxidase

The flavoenzyme nitroalkane oxidase (NAO) catalyzes the conversion of nitroalkanes to nitrite and aldehydes or ketones.94 The a-proton abstraction of the small substrate nitroethane by Asp402 is rate-limiting. The enzymatic rate enhancement is 109 over the uncatalyzed reaction in water.95 Of particular interest is the finding that the deuterium KIEs are noticeably greater for the enzymatic reaction (9.2) than that in water (7.8),95–96 suggesting that tunneling contribution may play a greater role in the enzyme than that in water. The proton transfer reaction between nitroethane and Asp402 catalyzed by nitroalkane oxidase provides an excellent

158

Quantum Mechanical Methods for Enzyme Modeling: Accurate Computation of Kinetic Isotope Effects

O CH3

C O-

H CH3

C H

O

O N+

CH3

O−

CH3

C

C OH

H

O− N+ O−

Scheme 4 Proton transfer from nitroethane to acetate (Asp402) ion. Major, D. T.; Garcia-Viloca, M.; Gao, J. Path integral simulations of proton transfer reactions in aqueous solution using combined QM/MM potentials. J. Chem. Theory Comput. 2006, 2, 236–245.

30

Free energy (kcal mol−1)

25

Table 1 Computed primary (11) and secondary (21) kinetic isotope effects, and computed and experimental total deuterium isotope effects for the proton transfer from nitroethane to Asp402 in NAO and to acetate ion in water. Subscripts and superscripts are used to specify the rate constant for isotope substitutions at the primary and secondary position, respectively

H2O, classical Quantum (D) Quantum (H)

20 15

NAO, classical Quantum (D) Quantum (H)

10 5 0 −5 −4

−2

0 2 Reaction coordinate (Å)

4

Figure 5 Classical (blue) and quantum mechanical potential of mean force for the proton (red) and deuteron (green) transfer from nitroethane to acetate ion in water (dashed curves), and to Asp402 (solid curves) in the active site of nitroalkane oxidase. Major, D. T.; Heroux, A.; Orville, A. M.; Valley, M. P.; Fitzpatrick, P. F.; Gao, J. Proc. Natl. Acad. Sci. 2009, 106, 20734. Copyright by PNAS.

opportunity for a comparative study because the enzymatic process can be directly modeled by the reaction of nitroethane and acetate ion in water. The PI-FEP/UM simulation technique was used to model the NQEs.42 Moreover, the PES was described by a combined QM/MM method.47,58 Thus, both the electronic structure of the reacting system and the nuclear dynamics are treated quantum mechanically. The crystal structure 2C12 with the inhibitor spermine was used to construct the Michaelis complex (MC) model through a combined minimization and molecular dynamics simulation technique.97 As it turns out, the computational MC structure was in excellent accord with a parallel study of the D402NNitroethane X-ray structure.86 The centroid PMFs from PI-FEP/UM simulations for the nitroethane deprotonation reaction in NAO and in water (Scheme 4) are shown in Figure 5. The reaction coordinate is defined as the difference between the breaking (donor–proton) and forming (acceptor–proton) bond distances (Scheme 4). The computed free energies of activation are 15.9 and 24.4 kcal mol1 for the enzymatic and the uncatalyzed proton transfer reaction in water, respectively, in accord with experimental results (14.0 and 24.8 kcal mol1).95 The computed primary and secondary KIEs for the nitroethane deprotonation reaction in water and in the enzyme are listed in Table 1 along with the total KIEs that have been determined for the perdeuterated substrate nitroethane at the Ca position.95 In PI-FEP/UM simulations, the computed

11

kHH =kDH

kHH =kTH

kDD =kTD

H2O NAO 21 H2O NAO kHH =kDD H2O NAO

6.6370.31 8.3670.58 kHH =kHD 1.34070.132 1.21370.150 Calc. 8.371.1 10.171.4

13.071.0 18.172.4 kHH =kHT 1.37570.183 1.22970.209 Expt.95 7.870.1 9.270.4

2.1770.04 2.3870.05 kDD =kDT 1.09670.048 1.05070.025

Major, D. T.; Garcia-Viloca, M.; Gao, J. Path integral simulations of proton transfer reactions in aqueous solution using combined QM/MM potentials. J. Chem. Theory Comput. 2006, 2, 236–245.

nuclear quantum effects include both zero-point energies and nuclear tunneling and their contributions to the computed KIEs are not separable.42 Thus, it is useful to use semiclassical methods98 to estimate the tunneling transmission coefficient to gain an understanding of the origin of NQEs in NAO catalysis. Consequently, we have utilized the EA-VTST/MT approach5,22,93 to separate tunneling from the overall quantum effects (Table 1). The results in Figure 5 and Table 1 present overwhelming support that nuclear quantum effects make greater contributions to the enzymatic reaction in nitroalkane oxidase than the uncatalyze process in water. First, there are greater NQEs to lower the classical barrier of the enzymatic process than that of the uncatalyzed model reaction. Second, the difference in quantum effects between proton and deuteron transfer is more significant for the enzymatic reaction than the uncatalyzed one. Third, both experiments and computations show enhanced primary KIEs in the enzyme over that in water. Fourth, the mixed Swain-Schaad exponent for the enzymatic reaction is greater than the semiclassical limit without tunneling and it is increased for the enzyme process. Finally, the tunneling transmission coefficient was found to be about three times greater for the enzyme reaction than the model reaction in water from EA-VTST/MT analysis, an entirely different approach from the PI-FEP/UM simulation method. Analyses of the tunneling paths in EA-VTST reveal that the origin of the difference may be attributed to a narrowing effect in the effective potentials for tunneling in the enzyme than that in aqueous solution.86 These studies demonstrate that differential quantum tunneling contributions are utilized in certain enzymatic catalysis as illustrated by nitroalkane oxidase.

Quantum Mechanical Methods for Enzyme Modeling: Accurate Computation of Kinetic Isotope Effects

9.9.4

Conclusion

In this chapter we presented a number of quantum mechanical methods developed for enzyme kinetics modeling; these include the treatment of the potential energy surface for reactive system and the incorporation of nuclear quantum effects in dynamics simulations. Two aspects are emphasized: (1) the potential energy surface is represented by combined QM/MM method and by a fully quantal X-Pol potential, in which electronic structure theory is used to describe bond forming and breaking processes; and (2) a free energy perturbation approach to transform a light (heavy) isotope into an heavy (light) counterpart in Feynman centroid path integral simulations have been developed to incorporate nuclear quantum effects along the centroid potential of mean force obtained from umbrella sampling simulations. The latter method is called the PI-FEP/UM approach, which coupled with a novel application of the bisection sampling technique allows primary and secondary kinetic isotope effects to be precisely computed for enzymatic reactions. For constructing reactive potential energy surfaces, a novel mixed MOVB theory has been developed, which can be used in context of combined QM/MM potential as the QM representation, and fully coupled with the quantal X-Pol force field as an integral part of the block-localization that includes nonadiabatic coupling. These methods have been illustrated by examining the potential surface for a hydride transfer in the gas phase, the decarboxylation reaction of N-methyl picolinate in water, the proton abstraction and reprotonation process catalyzed by alanine racemase, and the enhanced nuclear quantum effects in nitroalkane oxidase catalysis. These examples show that the incorporation of quantum mechanical effects is essential for enzyme kinetics simulations and offers a great opportunity to more accurately and reliably model the mechanism and free energies of enzymatic reactions.

Acknowledgments I wish to thank many coworkers and collaborators for their invaluable contributions to the project. Their names are listed in the references cited in this article. This work has been supported in part by the National Institutes of Health (GM46736) and by the National Science Foundation (CHE09-57162).

References [1] Wolfenden, R. Degrees of difficulty of water-consuming reactions in the absence of enzymes. Chem. Rev. 2006, 106, 3379–3396. [2] Garcia-Viloca, M.; Gao, J.; Karplus, M.; Truhlar, D. G. How enzymes work: analysis by modern rate theory and computer simulations. Science 2004, 303, 186–195. [3] Kohen, A.; Limbach, H. H., Eds. Isotope Effects in Chemistry and Biology; Taylor & Francis Group, CRC Press: New York, 2005. [4] Nagel, Z. D.; Klinman, J. P. Tunneling and dynamics in enzymatic hydride transfe. Chem. Rev. 2006, 106, 3095–3118. [5] Pu, J.; Gao, J.; Truhlar, D. G. Multidimensional tunneling, recrossing, and the transmission coefficient for enzymatic reactions. Chem. Rev. 2006, 106, 3140–3169.

159

[6] Klinman, J. P. Linking protein structure and dynamics to catalysis: the role of hydrogen tunnelling. Phil. Trans. Roy. Soc., B 2006, 361, 1323–1331. [7] Garcia-Viloca, M.; Alhambra, C.; Truhlar, D. G.; Gao, J. Inclusion of quantummechanical vibrational energy in reactive potentials of mean force. J. Chem. Phys. 2001, 114, 9953–9958. [8] Gao, J.; Ma, S.; Major, D. T.; Nam, K.; Pu, J.; Truhlar, D. G. Mechanisms and free energies of enzymatic reactions. Chem. Rev. 2006, 106, 3188–3209. [9] Schowen, R. L. In Transition States of Biochemical Processes; Gandour, R. D.; Schowen, R. L., Eds.; Plenum Press: New York, 1978; pp 77–114. [10] Devi-Kesavan, L. S.; Gao, J. Combined QM/MM study of the mechanism and kinetic isotope effect of the nucleophilic substitution reaction in haloalkane dehalogenase. J. Am. Chem. Soc. 2003, 125, 1532–1540. [11] Schramm, V. L. Enzymatic transition state poise and transition state analogues. Acc. Chem. Res. 2003, 36, 588–596. [12] Wolfenden, R.; Snider, M. J. The depth of chemical time and the power of enzymes as catalysts. Acc. Chem. Res. 2001, 34, 938–945. [13] Warshel, A.; Sharma, P. K.; Kato, M.; Xiang, Y.; Liu, H.; Olsson, M. H. M. Electrostatic basis for enzyme catalysis. Chem. Rev. 2006, 106, 3210–3235. [14] Eisenmesser, E. Z.; Millet, O.; Labeikovsky, W.; Korzhnev, D. M.; Wolf-Watz, M.; Bosco, D. A.; Skalicky, J. J.; Kay, L. E.; Kern, D. Intrinsic dynamics of an enzyme underlies catalysis. Nature 2005, 438, 117–121. [15] Kern, D.; Eisenmesser, E. Z.; Wolf-Watz, M. Enzyme dynamics during catalysis measured by NMR spectroscopy. Methods in Enzymology 2005, 394, 507–524. [16] Labeikovsky, W.; Eisenmesser, E. Z.; Bosco, D. A.; Kern, D. Structure and dynamics of Pin1 during catalysis by NMR. J. Mol. Bio. 2007, 367, 1370–1381. [17] Wolf-Watz, M.; Thai, V.; Henzler-Wildman, K.; Hadjipavlou, G.; Eisenmesser, E. Z.; Kern, D. Linkage between dynamics and catalysis in a thermophilicmesophilic enzyme pair. Nature Struct. Mol. Bio. 2004, 11, 945–949. [18] Boehr, D. D.; Dyson, H. J.; Wright, P. E. An NMR perspective on enzyme dynamics. Chem. Rev. 2006, 106, 3055–3079. [19] Boehr, D. D.; McElheny, D.; Dyson, H. J.; Wright, P. E. The dynamic energy landscape of dihydrofolate reductase catalysis. Science 2006, 313, 1638–1642. [20] Karplus, M. Aspects of protein reaction dynamics: deviations from simple behavior. J. Phys. Chem. B 2000, 104, 11–27. [21] Antoniou, D.; Basner, J.; Nunez, S.; Schwartz, S. D. Computational and theoretical methods to explore the relation between enzyme dynamics and catalysis. Chem. Rev. 2006, 106, 3170–3187. [22] Gao, J.; Truhlar, D. G. Quantum mechanical methods for enzyme kinetics. Ann. Rev. Phys. Chem. 2002, 53, 467–505. [23] Xie, W.; Orozco, M.; Truhlar, D. G.; Gao, J. X-Pol potential: an electronic structure-based force field for molecular dynamics simulation of a solvated protein in wate. J. Chem. Theory Comput. 2009, 5, 459–467. [24] Voth, G. A.; Chandler, D.; Miller, W. H. Rigorous formulation of quantum transition state theory and its dynamical corrections. J. Chem. Phys. 1989, 91, 7749–7760. [25] Voth, G. A. Feynman path integral formulation of quantum mechanical transition-state theory. J. Phys. Chem. 1993, 97, 8365–8377. [26] Messina, M.; Schenter, G. K.; Garrett, B. C. Centroid-density, quantum rate theory: variational optimization of the dividing surface. J. Chem. Phys. 1993, 98, 8525–8536. [27] Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965. [28] Gillan, M. J. Quantum simulation of hydrogen in metals. J. Phys. Rev. Lett. 1987, 58, 563–566. [29] Voth, G. A. Path-integral centroid methods in quantum statistical mechanics and dynamics. Adv. Chem. Phys. 1996, 93, 135–218. [30] Cao, J.; Voth, G. A. The formulation of quantum statistical mechanics based on the Feynman path centroid density. V. Quantum instantaneous normal mode theory of liquids. J. Chem. Phys. 1994, 101, 6184–6192. [31] Gehlen, J. N.; Chandler, D.; Kim, H. J.; Hynes, J. T. Free energies of electron transfer. J. Phys. Chem. 1992, 96, 1748–1753. [32] Grote, R. F.; Hynes, J. T. The stable states picture of chemical reactions. II. Rate constants for condensed and gas phase reaction models. J. Chem. Phys. 1980, 73, 2715–2732. [33] Chandler, D.; Wolynes, P. G. Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids. J. Chem. Phys. 1981, 74, 4078–4095.

160

Quantum Mechanical Methods for Enzyme Modeling: Accurate Computation of Kinetic Isotope Effects

[34] Makarov, D. E.; Topaler, M. Quantum transition-state theory below the crossover temperature. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1995, 52, 178–188. [35] Messina, M.; Schenter, G. K.; Garrett, B. C. A variational centroid density procedure for the calculation of transmission coefficients for asymmetric barriers at low temperature. J. Chem. Phys. 1995, 103, 3430–3435. [36] Mills, G.; Schenter, G. K.; Makarov, D. E.; Jonsson, H. Generalized path integral based quantum transition state theory. Chem. Phys. Lett. 1997, 278, 91–96. [37] Jang, S.; Voth, G. A. A relationship between centroid dynamics and path integral quantum transition state theory. J. Chem. Phys. 2000, 112, 8747–8757, Erratum: 2001. 8114: 1944 [38] Feynman, R. P.; Kleinert, H. Effective classical partition functions. Phys. Rev. A 1986, 34, 5080. [39] Sprik, M.; Klein, M. L.; Chandler, D. Staging: a sampling technique for the Monte Carlo evaluation of path integrals. Phys. Rev. B 1985, 31, 4234–4244. [40] Hwang, J. K.; Warshel, A. A quantized classical path approach for calculations of quantum mechanical rate constants. J. Phys. Chem. 1993, 97, 10053–10058. [41] Villa, J.; Warshel, A. Energetics and dynamics of enzymatic reactions. J. Phys. Chem. B 2001, 105, 7887–7907. [42] Major, D. T.; Gao, J. An integrated path integral and free-energy perturbationumbrella sampling method for computing kinetic isotope effects of chemical reactions in solution and in enzymes. J. Chem. Theory Comput. 2007, 3, 949–960. [43] Pollock, E. L.; Ceperley, D. M. Simulation of quantum many-body systems by path-integral methods. Phys. Rev. B 1984, 30, 2555–2568. [44] Ceperley, D. M. Path integrals in the theory of condensed helium. Rev. Mod. Phys. 1995, 67, 279–355. [45] Major, D. T.; Gao, J. Implementation of the bisection sampling method in path integral simulations. J. Mol. Graph. Model. 2005, 24, 121–127. [46] Major, D. T.; Garcia-Viloca, M.; Gao, J. Path integral simulations of proton transfer reactions in aqueous solution using combined QM/MM potentials. J. Chem. Theory Comput. 2006, 2, 236–245. [47] Major, D. T.; York, D. M.; Gao, J. Solvent polarization and kinetic isotope effects in nitroethane deprotonation and implications to the nitroalkane oxidase reaction. J. Am. Chem. Soc. 2005, 127, 16374–16375. [48] Major, D. T.; Gao, J. A combined quantum mechanical and molecular mechanical study of the reaction mechanism and a-amino acidity in alanine racemase. J. Am. Chem. Soc. 2006, 128, 16345–16357. [49] Gao, J. Methods and applications of combined quantum mechanical and molecular mechanical methods. In Reviews in Computers & Chemistry; Lipkowitz, K. B.; Boyd, D. B., Eds.; VCH: New York, 1995; Vol. 7, pp 119–185. [50] Truhlar, D. G.; Gao, J.; Alhambra, C.; Garcia-Viloca, M.; Corchado, J.; Sanchez, M. L.; Villa, J. The incorporation of quantum effects in enzyme kinetics modeling. Acc. Chem. Res. 2002, 35, 341–349. [51] Mo, Y.; Gao, J. Ab initio QM/MM simulations with a molecular orbital-valence bond (MOVB) method: application to an SN2 reaction in water. J. Comput. Chem. 2000, 21, 1458–1469. [52] Mo, Y.; Gao, J. An ab initio molecular orbital-valence bond (MOVB) method for simulating chemical reactions in solution. J. Phys. Chem. A 2000, 104, 3012–3020. [53] Song, L.; Gao, J. On the construction of diabatic and adiabatic potential energy surfaces based on ab initio valence bond theory. J. Phys. Chem. A 2008, 112, 12925–12935. [54] Gao, J. Toward a molecular orbital derived empirical potential for liquid simulations. J. Phys. Chem. B 1997, 101, 657–663. [55] Gao, J. A molecular-orbital derived polarization potential for liquid water. J. Chem. Phys. 1998, 109, 2346–2354. [56] Xie, W.; Gao, J. Design of a next generation force field: The X-POL potential. J. Chem. Theory Comput. 2007, 3, 1890–1900. [57] Xie, W.; Song, L.; Truhlar, D. G.; Gao, J. Incorporation of QM/MM buffer zone in the variational double self-consistent field method. J. Phys. Chem. B 2008, 112, 14124–14131. [58] Gao, J.; Xia, X. A prior evaluation of aqueous polarization effects through Monte Carlo QM-MM simulations. Science 1992, 258, 631–635. [59] Senn, H. M.; Thiel, W. QM/MM methods for biomolecular systems. Angew. Chem., Int. Ed. 2009, 48, 1198–1229. [60] Gao, J.; Furlani, T. R. Simulating solvent effects in organic chemistry: combining quantum and molecular mechanics. IEEE Comput. Sci. Eng. 1995, Fall Issue, 24–33.

[61] Garcia-Viloca, M.; Truhlar, D. G.; Gao, J. Importance of substrate and cofactor polarization in the active site of dihydrofolate reductase. J. Mol. Biol. 2003, 327, 549–560. [62] Bixon, M.; Lifson, S. Potential functions and conformations in cycloalkanes. Tetrahedron 1967, 23, 769–784. [63] Xie, W.; Song, L.; Truhlar, D. G.; Gao, J. The variational explicit polarization potential and analytical first derivative of energy: towards a next generation force field. J. Chem. Phys. 2008, 128, 234108/234101–234108/234109. [64] Song, L.; Han, J.; Lin, Y. L.; Xie, W.; Gao, J. Explicit polarization (X-Pol) potential using ab initio molecular orbital theory and density functional theory. J. Phys. Chem. A 2009, 113, 11656–11664, doi:10.1021/jp902710a. [65] Mo, Y.; Peyerimhoff, S. D. Theoretical analysis of electronic delocalization. J. Chem. Phys. 1998, 109, 1687–1697. [66] Mo, Y.; Gao, J.; Peyerimhoff, S. D. Energy decomposition analysis of intermolecular interactions using a block-localized wave function approach. J. Chem. Phys. 2000, 112, 5530–5538. [67] Cembran, A.; Song, L.; Mo, Y.; Gao, J. Block-localized density functional theory (BLDFT), diabatic coupling, and its use in valence bond theory for representing reactive potential energy surfaces. J. Chem. Theory Comput. 2009, 5, 2702–2716. [68] Pople, J. A.; Santry, D. P.; Segal, G. A. Approximate self-consistent molecular orbital theory. I. Invariant procedures. J. Chem. Phys. 1965, 43, S129–S135. [69] Pople, J. A.; Segal, G. A. Approximate self-consistent molecular orbital theory. II. Calculations with complete neglect of differential overlaop. J. Chem. Phys. 1965, 43, S136. [70] Brooks, B. R.; Brooks, C. L.; Mackerell, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, D. M.; Karplus, M. J. CHARMM: the biomolecular simulation program. J. Comput. Chem. 2009, 30, 1545–1614, doi:10.1002/ Jcc.21287. [71] Aqvist, J.; Warshel, A. Simulation of enzyme reactions using valence bond force fields and other hybrid quantum/classical approaches. Chem. Rev. 1993, 93, 2523–2544. [72] Song, L.; Mo, Y.; Gao, J. An effective hamiltonian molecular orbital-valence bond (MOVB) approach for chemical reactions as applied to the nucleophilic substitution reaction of hydrosulfide ion and chloromethane. J. Chem. Theory Comput. 2009, 5, 174–185. [73] Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: a program for macromolecular energy, minimization and dynamics calculations. J. Comput. Chem. 1983, 4, 187. [74] Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. Development and use of quantum mechanical molecular models. 76. AM1: a new general purpose quantum mechanical molecular model. J. Am. Chem. Soc. 1985, 107, 3902–3909. [75] Stewart, J. J. P. Optimization of parameters for semiempirical methods I. Method. J. Comp. Chem. 1989, 10, 209–220. [76] Elstner, M.; Porezag, D.; Juugnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Sukai, S.; Seifect, G. Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties. Phys. Rev. B 1998, 58, 7260–7268. [77] Warshel, A.; Weiss, R. M. An empirical valence bond approach for comparing reactions in solutions and in enzymes. J. Am. Chem. Soc. 1980, 102, 6218–6226. [78] Sato, S. On a new method of drawing the potential energy surface. J. Chem. Phys. 1955, 23, 592. [79] Raff, L. M. Theoretical investigations of the reaction dynamics of polyatomic systems: chemistry of the hot atom (T* + CH4) and (T + CD4) systems. J. Chem. Phys. 1974, 60, 2220. [80] Silver, D. M.; Brown, N. J. Valence bond model potential energy surface for H4. J. Chem. Phys. 1980, 72, 3859. [81] Schlegel, H. B.; Sonnenberg, J. L. Empirical valence-bond models for reactive potential energy surfaces using distributed Gaussians. J. Chem. Theory Comput. 2006, 2, 905–911. [82] Cembran, A.; Payak, A.; Lin, Y. L.; Xie, W.; Song, L.; Mo, Y.; Gao, J. A nonorthogonal block-localized effective Hamiltonian approach for chemical and enzymatic reactions. J. Chem. Theory Comput. 2010, 6, 2242. [83] Hwang, J.-K.; Warshel, A. How important are quantum mechanical nuclear motions in enzyme catalysis? J. Am. Chem. Soc. 1996, 118, 11745–11751.

Quantum Mechanical Methods for Enzyme Modeling: Accurate Computation of Kinetic Isotope Effects

[84] Hinsen, K.; Roux, B. Potential of mean force and reaction rates for proton transfer in acetylacetone. J. Chem. Phys. 1997, 106, 3567–3577. [85] Gao, J.; Wong, K.-Y.; Major, D. T. Combined QM/MM and path integral simulations of kinetic isotope effects in the proton transfer reaction between nitroethane and acetate ion in water. J. Comput. Chem. 2008, 29, 514–522. [86] Major, D. T.; Heroux, A.; Orville, A. M.; Valley, M. P.; Fitzpatrick, P. F.; Gao, J. Differential quantum tunneling contributions in nitroalkane oxidase catalyzed and the uncatalyzed proton transfer reaction. Proc. Natl. Acad. Sci. 2009, 106, 20734. [87] Gao, J.; Garcia-Viloca, M.; Poulsen, T. D.; Mo, Y. Solvent effects, reaction coordinates, and reorganization energies on nucleophilic substitution reactions in aqueous solution. Adv. Phys. Org. Chem. 2003, 38, 161–181. [88] Rishavy, M. A.; Cleland, W. W. Determination of the mechanism of orotidine 50 -monophosphate decarboxylase by isotope effects. Biochemistry 2000, 39, 4569–4574. [89] Major, D. T.; Gao, J. An integrated path integral and free-energy perturbation– umbrella sampling method for computing kinetic isotope effects of chemical reactions in solution and in enzymes. J. Chem. Theory Comput. 2007, 3, 949. [90] Nam, K.; Gao, J.; York, D. M. An efficient linear-scaling ewald method for long-range electrostatic interactions in combined QM/MM calculations. J. Chem. Theory Comput. 2005, 1, 2–13. [91] Major, D. T.; Nam, K.; Gao, J. Transition state stabilization and a-amino carbon acidity in alanine racemase. J. Am. Chem. Soc. 2006, 128, 8114–8115.

161

[92] Spies, M. A.; Woodward, J. J.; Watnik, M. R.; Toney, M. D. Alanine racemase free energy profiles from global analyses of progress curves. J. Am. Chem. Soc. 2004, 126, 7464–7475. [93] Alhambra, C.; Corchado, J.; Sanchez, M. L.; Garcia-Viloca, M.; Gao, J.; Truhlar, D. G. Canonical variational theory for enzyme kinetics with the protein mean force and multidimensional quantum mechanical tunneling dynamics. Theory and application to liver alcohol dehydrogenase. J. Phys. Chem. B 2001, 105, 11326–11340. [94] Fitzpatrick, P. F.; Orville, A. M.; Nagpal, A.; Valley, M. P. Nitroalkane oxidase, a carbanion-forming flavoprotein homologous to acyl-CoA dehydrogenase. Arch. Biochem. Biophys. 2005, 433, 157–165. [95] Valley, M. P.; Fitzpatrick, P. F. Comparison of enzymatic and non-enzymatic nitroethane anion formation: thermodynamics and contribution of tunneling. J. Am. Chem. Soc. 2004, 126, 6244–6245. [96] Valley, M. P.; Tichy, S. E.; Fitzpatrick, P. F. Establishing the kinetic competency of the cationic imine intermediate in nitroalkane oxidase. J. Am. Chem. Soc. 2005, 127, 2062–2066. [97] Nagpal, A.; Valley, M. P.; Fitzpatrick, P. P.; Orville, A. M. Crystal structures of nitroalkane oxidase: Insights into the reaction mechanism from a covalent complex of the flavoenzyme trapped during turnover. Biochemistry 2006, 45, 1138–1150. [98] Fernandez-Ramos, A.; Miller, J. A.; Klippenstein, S. J.; Truhlar, D. G. Modeling the kinetics of bimolecular reactions. Chem. Rev. 2006, 106, 4518–4584.