Origin of the Temperature Dependence of Isotope Effects in Enzymatic

Jun 16, 2007 - The origin of the temperature dependence of kinetic isotope effects (KIEs) in enzyme reactions is a problem of general interest and a m...
1 downloads 8 Views 270KB Size
Download PDF
7852

J. Phys. Chem. B 2007, 111, 7852-7861

Origin of the Temperature Dependence of Isotope Effects in Enzymatic Reactions: The Case of Dihydrofolate Reductase Hanbin Liu and Arieh Warshel* Department of Chemistry, UniVersity of Southern California, 3620 McClintock AVenue, Los Angeles, California 90089-1062 ReceiVed: February 2, 2007; In Final Form: April 30, 2007

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.

I. Introduction In numerous enzymatic reactions, one of the key steps, which is sometimes also the rate-limiting step, is a transfer of proton, hydride, or hydrogen. Since these are comparatively light particles, the transfer is often associated with nuclear quantum mechanical (NQM) effects (e.g., zero-point energy and tunneling corrections). Thus, it is of significant interest to understand the magnitude of NQM effects in enzymatic reactions and to explore their possible contribution to enzyme catalysis.1-8 Furthermore, the temperature dependence of the NQM effects as manifested in the corresponding kinetic isotope effect (KIE) is a topic of significant current interest.1,9-12 A part of this interest is due to the hope that this temperature dependence can provide useful information about possible dynamical contributions to enzymatic reactions.1,10 The development of methods for simulating NQM effects in enzymatic reactions dates back to the early 1990s (see, e.g., refs 6, 8, 13, and 14) and has recently been quite an active field.7,9,15-18 Such calculations can be used in principle in assessing the catalytic contribution of NQM and in analyzing the temperature dependence of KIEs. However, the theoretical challenges are quite difficult in each of these tasks. That is, alhtough theoretical studies of the catalytic contributions of NQM have progressed significantly (e.g., see refs 5, 7-9, and 16-20), there are major open problems and major challenges of reproducing the observed temperature effects of KIE by actual microscopic simulations. Apparently, it is possible to account for the observed temperature dependence by using phenomenological vibronic treatments (e.g., refs 11-13) and adjusting * To whom correspondence should be addressed. E-mail: warshel@ usc.edu.

key parameters (as was done in refs 1 and 5). Unfortunately, the vibronic treatment involves fundamental problems (see refs 17 and 18 and this work), and the parameters obtained by the fitting procedure are unlikely to provide a quantitative description of the reacting system. On the other hand, it is extremely challenging to evaluate the temperature dependence by computer simulations that do not involve any special parametrization. In fact, we are not aware of any successful attempts to reproduce the temperature dependence of KIE of an enzymatic reaction in a quantitative way by a first principle simulation, despite attempts to do so (e.g., refs 9, 19). On the other hand, studies of the temperature dependence of NQM in gas-phase reaction seem to provide more quantitative results.21 Another interesting issue is the relationship between the KIE and the distance between the donor and acceptor in hydride transfer (HT), protein transfer (PT) and hydrogen transfer reaction. It is tempting to assume that enzyme catalyzes reactions by compressing the distance between the donor and acceptor (e.g., ref 22) or the related idea of near attack conformation (NAC) (e.g., ref23). This would mean that the enzyme leads to change in the shape of the barrier (and presumably enhances the contribution from tunneling to catalysis23-29). This idea was also used to rationalize the temperature dependence of observed KIEs (e.g., ref 27). On the other hand, our view has been that tunneling effects do not contribute significantly to catalysis (e.g., ref 9, 17, 30) and that enzymes do not apply strong stress that leads to a significant “compression” relative to the corresponding situation in solution (e.g., ref 18). Furthermore, our studies (e.g., ref 9, 17) and those of others20,31 started to indicate that the KIE and the relative NQM contribution increase, rather than decrease, when the donor-acceptor distance increases. Obvi-

10.1021/jp070938f CCC: $37.00 © 2007 American Chemical Society Published on Web 06/16/2007

Isotope Effects in Enzymatic Reactions

J. Phys. Chem. B, Vol. 111, No. 27, 2007 7853

Figure 1. The reaction system in DHFR. All the atoms shown (except the NHR and R atoms) are included in region I of the EVB treatment.

ously it is important to try to resolve this issue by consistent microscopic simulations. The present work attempts to elucidate the origin of the temperature dependence of the KIE in proteins by simulating this effect in different mutants of dihydrofolate reductase that has been the subject of careful recent studies (e.g., ref 10). The use of the quantum classical path (QCP) and the empirical valance bond (EVB) approaches allows us to reproduce the observed trend while considering explicitly the NQM effects of the protein + substrate + solvent system. II. Methods II.1. The EVB Surface. In order to simulate chemical reactions in proteins, it is essential to have a reasonable potential energy surfaces. In the case of NQM calculations, which require long simulations to obtain convergence of the corresponding free energy, it is also helpful if the relevant potential energy surface is available in an analytical form. In this respect, the EVB approach provides a very powerful tool that will be used in this work. The EVB method has been described extensively elsewhere (e.g., ref 32-34), and we give here only a few key points. The EVB method is a QM/MM method that describes reactions by mixing diabatic states that correspond to classical valence-bond (VB) structures. These states represent the reactant, intermediate (or intermediates), and product states. The potential energies of these diabatic states are represented by classical MM potential functions of the form i Hii ) i ) Rigas + Uiintra(R, Q) + USs (R, Q, r, q) + Uss(r, q) (1)

Here, R and Q represent the atomic coordinates and charges of the diabatic states, and r and q are those of the surrounding protein and solvent. Rigas is the gas-phase energy of the ith diabatic state (where all the fragments are taken to be at infinity), Uintra(R, Q) is the intramolecular potential of the solute system (relative to its minimum), and USs(R, Q, r, q) represents the interaction between the solute (S) atoms and the surrounding (s) solvent and protein atoms. Uss(r, q) represents the potential energy of the protein/solvent system (“ss” designates surrounding-surrounding). The i of eq 1 form the diagonal elements of the EVB Hamiltonian (HEVB). The off-diagonal elements of this Hamiltonian, Hij, are assumed to be constant or to be represented by a simple function, such as an exponential function of the distances between the reacting atoms. These Hij elements are assumed to be the same in the gas phase, in solutions, and in the proteins. The adiabatic ground state energy, Eg, and the corresponding eigenvector, Cg, are obtained by solving the secular equation HEVBCg ) EgCg

(2)

The EVB treatment provides a natural picture of intersecting electronic states, which is useful for exploring environmental effects on chemical reactions in condensed phases.32 The ground-state charge distribution of the reacting species (“solute”) polarizes the surroundings (“solvent”), and the charges of each resonance structure of the solute then interact with the polarized solvent.32 This coupling enables the EVB model to capture the effect of the solvent on the quantum mechanical mixing of the different states of the solute. For example, in cases that ionic and covalent states are describing the solute, the solvent stabilizes the ionic state to a greater extent, and the resulting ground state has more ionic character and more solvation energy. The EVB approach obtains the free energy surface in two steps. In the first step, it accumulates the statistical information about the surface by changing the system adiabatically from one diabatic state to another. In the simple case of two diabatic states, this “mapping” potential, m, can be written as a linear combination of the reactant and product potentials, 1 and 2 m ) (1 - θm) 1 + θm2

(0 e θm e 1) (3)

when θm is changed from 0 to 1 in n + 1 fixed increments (θm ) 0/n, 1/n, 2/n, ..., n/n) potentials. The free energy, ∆Gm, associated with changing θm from 0 to m/n is evaluated by the FEP procedure described elsewhere.33 The second step evaluates the free energy functional that corresponds to the adiabatic ground state surface, Eg, is then obtained by the FEP-umbrella sampling (FEP-US) method,33 which can be written as ∆g(x′) ) ∆Gm - β-1 ln〈δ(x - x′) × exp[-β(Eg(x) - m(x))]〉m (4) where m is the mapping potential that keeps x in the region of x′. If the changes in m are sufficiently gradual, the free energy functional ∆g(x′) obtained with several values of m overlap over a range of x′, and patching together the full set of ∆g(x′) gives the complete free energy curve for the reaction. The simulation system was constructed on the basis of the X-ray crystallographic structure of Escherichia coli DHFR (PDB id: 1RX2). The corresponding reacting system is depicted in Figure 1. The G121V, M42W-G121V mutants were generated by starting from the wild type DHFR structure and running long relaxation runs. The relaxation and the FEP-US simulations were performed by the standard all-atom simulation procedure of the ENZYMIX module of the MOLARIS simulation package.35,36 The simulation system was divided to four regions: region I included the folate and the NADPH (see Figure 1); region II included the protein plus the water molecules in and around the protein up to a radius of 20 Å; and region III included protein atoms and water molecules that were subjected to distance and polarization constraints according to the surface-constrained all-

7854 J. Phys. Chem. B, Vol. 111, No. 27, 2007

Liu and Warshel

atom solvent (SCAAS) boundary condition.37 The rest of the system was represented by a bulk region with a dielectric constant of 80. The long-range electrostatic effects were treated by the local reaction field (LRF) method38 that provides one of the most rigorous ways of treating electrostatic effects in protein and other nonperiodic systems. More details about the classical simulations are given in ref 39. II.2. QCP Simulations. With the analytical EVB surface of the reacting system and its surrounding protein + water system, our task is to obtain the quantum correction to the classical activation free energy. This is done by using the quantum classical path (QCP) centroid path integral approach developed in our previous studies of NQM effects in chemical reactions in solution and proteins.6,8 Note that the QCP has been adopted recently by other research groups (e.g., refs 40-42). In the QCP approach, the nuclear quantum mechanical rate constant is expressed as kqm ) Fqm kBT/h exp(-β∆gqqm)

(5)

where Fqm, kB, T, h, and β are, respectively, the transmission factor, Boltzmann’s constant, the temperature, and Planck’s constant, and β ) 1/kΒT. The quantum mechanical activation barrier, ∆gqqm, includes almost all the nuclear quantum mechanical effects, whereas only small effects come from the preexponential transmission factor in the case of systems with a significant activation barrier.7,43 The quantum mechanical free energy barrier, ∆gqqm, can be evaluated by Feynman’s path integral formulation,44 in which each classical coordinate is replaced by a ring of quasiparticles that are subjected to the effective “quantum mechanical” potential. p

Uqm )

1

∑ 2pMΩ ∆x

2

2

k

k)1

1 + U(xk) p

(6)

Here, ∆xk ) xk+1 - xk (where xp+1 ) x1), Ω ) p/pβ, M is its mass, and U is the actual potential used in the classical simulation. The total quantum mechanical partition function can then be obtained by running classical trajectories of the quasiparticles with the potential Uqm. The probability of being at the transition state is approximated by a probability distribution of the center of mass of the quasiparticles (the centroid) rather than the classical single point. Thus, we follow the strategy of the centroid path integral method.45-47 However, the regular centroid approach requires one to consider the trajectory of all the quasiparticles, and thus, the use of this approach in condensed phase reactions is very challenging and may have major convergence problems. The QCP approach offers an effective and rather simple way for evaluating this probability without changing the simulation program significantly. This is done by propagating classical trajectories on the classical potential energy surface of the reacting system and using the positions of the atom of the system to generate the centroid position for the quantum mechanical partition function. This treatment is based on the finding that the quantum mechanical partition function can be expressed as6,8,14,17 Zq(xj) ) Zcl(xj)〈〈exp{-(β/p)

∑ U(x ) - U(xj)}〉 k

fp〉U

(7)

k

where xj is the centroid position, 〈...〉fp designates an average over the free particle quantum mechanical distribution obtained with the implicit constraint that xj coincides with the current position of the corresponding classical particle, and 〈...〉U

designates an average over the classical potential U. Using eq 7, we can obtain the quantum mechanical free energy surface by evaluating the corresponding probability by the same combined FEP-US of eq 4. Now we use the double average of eq 7 rather than an average over a regular classical potential. At any rate, the main point of the QCP is that the quantum mechanical free energy function can be evaluated by a centroid approach that is constrained to move on the classical potential. This provides stable and relatively fast converging results that have been shown to be quite accurate in studies of well-defined test potentials, or which the exact quantum mechanical results are known (e.g., ref 8). In order to obtain the NQM effect, we run classical trajectories of the quasiparticles on the “quantum mechanical” potential given by eq 6. For these trajectories, we have used 15 particles (less than 10 particles seems to give unstable results). In addition, since any reliable protein investigation should average the relevant energies over different protein configurations, the energy of the quasiparticle trajectory has been averaged for every 0.5 fs of a 20 ps protein simulation. In the end, the ∆∆gqclfqm term is averaged over 20 000 protein and quasiparticle configurations. In order to compensate for the comparatively short simulation over the protein configurations (it turns out to be very time-consuming to average over both protein and quasiparticle configurations with a 15 particle system), the system was first equilibrated thoroughly for 500 ps. II.3. A Simplified Vibronic Treatment. The EVB method provides an effective way to explore the effect of the fluctuations of the environment on PT and HT reactions. That is, the electrostatic potential from the fluctuating polar environment interacts with the charge distribution of each resonance structure, and thus, the fluctuation of the environment is directly included in the time dependence of the EVB Hamiltonian. This point has emerged from our early studies13,15,48,49 and discussed in these works. Here, we will only consider its implication with regard to the NQM effects. As discussed elsewhere (e.g., ref 13), we can use the autocorrelation of the fluctuating EVB energy gap and obtain an approximated qusiharmonic rate constant. Briefly, our starting point is the overall rate constant, kab )

∑k

am,bm′

exp{-Eamβ}/

mm′

∑ exp{-E

amβ}

(8)

m

where β ) 1/(kBT) (with kB designating the Boltzmann constant) and Eam is the energy of the mth vibronic level of state a. By eq 8, we assume that vibrational states in the reactant well are populated according to Boltzmann distribution. The individual vibronic rate constant, kam,bm′, is evaluated by monitoring the energy difference between Eam and Ebm as a function of the fluctuations of the rest of the molecule. The use of the corresponding fluctuating energy gap and a cumulant expansion (see ref 50) gives kam,bm′ ) |{Hab}/{p}|2

∑S

2 mm′

∫

∞

-∞

exp[(i/p)〈∆bm′,am〉 + γ(t) dt]

γ(t) ) -(i/p)2

∫ (t - t′)〈∆(0)∆(t′)〉

a

dt′

(9)

where the Sm,m′ is the Franck-Condon factor for transition from m to m′ and Hab is the off-diagonal electronic matrix element of the EVB Hamiltonian. Here, u is the electronic energy gap relative to its average values, given by u ) b - a - 〈∆ba〉a

(10)

Isotope Effects in Enzymatic Reactions

J. Phys. Chem. B, Vol. 111, No. 27, 2007 7855

and 〈 〉a designates an average obtained over the fluctuations around the minimum of state a. In the high-temperature limit, one obtains13,50 kam,bm′ ) |HabSmm′/p|2(πp2/kBTλ)1/2 exp{-∆gqmmβ}

(11)

where λ is the “solvent reorganization energy” defined by λ ) 〈∆ba〉a - ∆G0

∑pω (m′ - m ) + λ] /4λ 2

r

r

(13)

r

(16)

This expression is given for a single ∆r, and a more consistent expression can be obtained by integrating the vibronic rate constant over the “soft” coordinates, and in particular the D‚‚‚A distance. This can be done by writing khab )

(12)

Basically, eq 11 reflects the probability of vibronic transition from the reactant well to the product well (as determined by the vibrational overlap integrals (the Smm′) modulated by the chance that ∆ will be zero. This chance is determined by the activation free energy, ∆gq, whose value can be approximated by q ∆gmm′ ≈ [∆G0 +

kHa0,b0(∆r)/kDa0,b0(∆r) ≈ exp{24∆r2}

∫k

ab(∆r)

exp{-w(∆r)} dr

(17)

where w(∆r) is the potential of mean force (PMF) for the donor-acceptor distance. If the main contribution to kab comes from ka0,b0, then we can approximate the KIE by KIE = khHa0,b0/khDa0,b0 =

∫ F(∆r, T)e

-58∆ 2

r dr /

∫ F(∆r, T)e

-58∆r2x2

dr (18)

where

r

where ω is the vibrational frequency of the quantized mode. This relationship is applicable only if the system can be described by the linear response approximation (see ref 50), but this does not require that the system will be harmonic. The above vibronic treatment is similar to the expression developed by Kuznetsov and Ulstrup.11 However, the treatment that leads to eq 13, which was developed by Warshel and co-workers,13,48,50 is based on a more microscopic approach and leads to much more consistent treatment of ∆gq, where we can use rigorously ∆G0 rather than ∆E. Furthermore, our dispersed polaron (spin boson) treatment50 of eq 9 and, if needed, eq 13, gives a clear connection between the spectral distribution of the solvent fluctuations and the low-temperature limit of eq 9. It is also useful to note that Borgis and Hynes51 and Antoniou and Schwartz52 have used a very similar treatment but considered only the lowest vibrational levels of the proton. Unfortunately, the use of the above vibronic treatment is only valid in the diabatic limit when HabSmn2 is sufficiently small. Now, in cases of PT and HT processes, Hab seems to be far too large to allow for a diabatic treatment. However, it is more significant to ask what is the magnitude of HabSmn2. The FranckCondon factor is determined by the origin shift, ∆r ) R - 2b (where R is the donor-acceptor distance and b is the equilibrium D-H distance, which is typically around 1.05 Å). ∆r is usually converted to a dimensionless origin shift, ∆, with the proper conversion (e.g., see refs 53,54). Here, we note that for 0 f 0 transitions, Smn may be quite small, since it is given by S00 ) exp{-∆2/4}. Now, using the fact that ∆ is 6.5 for ∆r ) 0.6 Å, we obtain for HT, PT, and hydrogen transfer the following: ∆H = 6.5(∆r/0.6)(mHωH/3000)1/2 ∆D = 6.5(∆r/0.6)(mH × 2 × ωH/(2) /3000) 1/2

where ∆r, ω, and M are given in Å, Thus, we have

cm-1,

F(∆r, T) ) [Hab2(∆r)/(λ(∆r))1/2] × exp{-[w(∆r) + ∆gq00 (∆r)]β} (19) If we assume that λ and Hab are independent of ∆r in the region with the largest contribution to k00, we can write KIE =

∫ exp{-w (∆r)β} exp(-58∆r ) dr / ∫ exp{-w (∆r)β} exp{-58∆r x2} dr (20) 2

q

q

2

where wq ) w + ∆gq00. Equation 20 can be used as a qualitative guide, but we must keep in mind that we are dealing with a very problematic approximation. The general approximation that leads to eqs 9 and 11 is the requirement that we are dealing with the diabatic approximation, which means that the Landau Zener approximation is valid for each term in eq 9 (that reflects surface crossing between the am and bm′ vibronic levels). This requires (see ref 50) that mm′ ΓL,Z ) 1 - exp{-2π(HabSmm′)2/pu˘ }

∝ -2π(HabSmm′)2/pu˘ e 1

where u˘ is the change of the vibronic energy gap with time, at the crossing region (namely, ∂∆Ebm′,am/∂t). u˘ can be expressed in terms of the dimensionless solvent coordinate, Q, as (ref 50) u ) pωQ∆, where ω is the effective solvent frequency and ∆ is the solvent origin shift. Using the equal partition relationship, we have Q˙ ) (kBTω/p)1/2. We also have the relationship 2λ ) pω∆2; thus, we can write Γmm' ≈ 2π(HabSmm′)2/(pω(2λkBT/p2)1/2)

1/2

(14)

and amu, respectively.

kHa0,b0 ∝ Hab2 exp(-∆H2/4) × exp(-∆gq00β) = Hab2 exp(-58 × ∆r2) × exp{-∆g00β} kDa0,b0 ∝ Hab2 exp(-∆D2/4) × exp(-∆gq00β) = Hab2 exp(-58 × ∆r2 × x2) × exp{-∆gq00β} (15) where we assumed that ωH is ∼3000 cm-1 (as will be discussed below, this assumption is problematic). With this we can write

(21)

(22)

A similar expression has been estimated by related consideration (e.g., ref 11, 55). For a typical value of Hab ) 20 kcal/ mol, λ ) 20 kcal/mol, ω ) 100 cm-1, and ∆r ) 0.6 Å, we can obtain Γ00 ∼ 0.05. However, for ∆r ) 0.4 Å, we get Γ ∼5 × 103, which is clearly out of the diabatic regions. The estimate of ∆r > 0.5 Å can be an underestimate, since Hab may be much larger than 20 kcal/mol if the diabatic surfaces are fully harmonic (with ωH ) 3000 cm-1), as implied by the vibronic formulation. Thus, for ∆r . 0.6 Å and for zero-zero transition, we satisfy that LZ conditions, but for smaller ∆r or for hot transition (where m q 0), the diabatic approximation is not valid. A similar conclusion about the problem with the diabatic approximation was reached by Kuznetsov, Ulstrup, and co-

7856 J. Phys. Chem. B, Vol. 111, No. 27, 2007

Liu and Warshel

Figure 3. An example of the nuclear quantum contribution to the free energy profile at the reactant state and transition state when the donoracceptor distance is kept at 3.0 Å in native DHFR. Figure 2. A semiclassical vibronic treatment of proton transfer. This model, which is valid only for small H12, treats the donor-hydrogen stretching vibration quantum mechanically and the rest of the system classically. In this way, we monitor the energy gap between the vibronic states 1 + hωH/2(n1 + 1/2) and 2 + hωH/(n2 + 1/2) for trajectories of the system with a fixed X-H bond length. The figure depicts the time dependence of 1, 2, and 2 plus single and double excitations of the X-H bond and also provides the energy levels at two points on the trajectory. A semiclassical surface-hopping treatment of the crossing probability between the vibronic states, due to the fluctuating energy gap, leads to eq 8 (see ref 13).

workers.11,56 Furthermore, as will be discussed below, we face other problems with regard to the parameters in eq 11. First, if we consider a collinear PT or HT and treat all the coordinates except the X-H stretching frequency, then we have to deal with a very large intramolecular reorganization energy (see ref 9). Second, the effective frequency for the X-H stretch can be very different from the typical frequency of ∼3000 cm-1 once the X‚‚‚Y distance starts to be shorter than 3.2 Å. In this range, H12 starts to affect the ground state curvature in a drastic way (see below). Here, one can use the idea introduced by Warshel and Chu13 and modify the diabatic potential to make it close to the adiabatic potential. As long as the main contribution to the rate constant comes from the S00 term, it is reasonable to represent this effect by using ωH(∆) ) ω0H(1 - (4H12/λH)2R)

(23)

where R is ∼0.5. It is also possible to use the empirical ωH(D) of ref 57. Finally, it may also be useful to account for the complications due to the fact that the intramolecular solute coordinate contributes to the apparent reorganization energy and is also coupled to the hydride transfer coordinate (e.g., ref 18). III. Results and Discussion Before studying the NQM effects, it is important to have a reasonable description of the classical activation free energy in the enzyme studied. Here, we started with the results from our recent study39 of DHFR and its mutants, where we reproduced in a semiquantitative way the effect of different mutations on the classical rate constant. Using the same EVB parameters, we performed QCP calculations of the quantized activation free energy for the H and D reactions of the native DHFR58 the G121V mutant59 and the M42W-G121V double mutant.10 The calculations were performed for donor and acceptor distances

Figure 4. The temperature dependence of the KIE in native DHFR for the case when the donor-acceptor distance is kept at 3.0 Å (a) and 3.5 Å (b).

of 3.0 and 3.5 Å and for different temperatures. Typical results for the NQM contribution to the free energy profile at the reactant state equilibrium geometry and at the transition state are given in Figure 3. The corresponding activation barriers were converted to rate constants using eq 5, and the resulting isotope effects were averaged on three independent calculations with different initial conditions. The average results are summarized in Figures 4, 5, and 6. As seen from the figures and as expected from eq 16 and from our previous studies,17,18,60 the KIE increases significantly upon increase of the donor and acceptor distance. Apparently, the change of the KIE with temperature at a given distance is smaller than the changes upon increase in distance. Thus, it is very likely that the observed temperature dependence reflects to a large extent the change in donor-acceptor distance from 3.0 to 3.5 Å. In order to take into account the distance dependence of the KIE, it is useful to follow the spirit of eq 17 and write

∫ k(R) e ) ∫ (k T/h)e

-w(R)β

khqm ≈

qm

dR/

q -(∆gcl(R)β

∫e

-w(R)β

dR

+ w(R)β) -∆∆gqm(R)β q

e

B

dR/

∫e

-w(R)β

dR (24)

where ∆∆gqqm is the quantum correction and R is the donoracceptor distance. We can also write khqm )

∫ P(R) (k T/h)e q cl

B

-∆∆gqm (R)β q

(25)

Isotope Effects in Enzymatic Reactions

J. Phys. Chem. B, Vol. 111, No. 27, 2007 7857

Figure 5. The temperature dependence of the KIE in the G121V mutant of DHFR for the case when the donor-acceptor distance is kept at 3.0 Å (a) and 3.5 Å (b).

Figure 7. The PMF for the donor-acceptor distance for the systems considered in Figures 3, 4, and 5.

Figure 6. The temperature dependence of the KIE in the M42WG121V mutant of DHFR for the case when the donor-acceptor distance is kept at 3.0 Å (a) and 3.5 Å (b).

where P(R) ) e-(∆

gclq(R)β + w(R)β)

∫e

-w(R)β

dR ) e-∆w (R)β) #

∫e

-w(R)β

dR (26)

is the classical probability of being at the transition state for different values of the donor and acceptor distance and w is the potential of mean force for the distance between the donor and acceptor (Figure 7). Note that this PMF is rather qualitative due to well-known convergence problems61 that increase when one considers large distance changes. The PMF presented in Figure 7 represents more extensive study than the one used in ref 39, but still reflects major uncertainties due to very significant hysteresis and is being used here only as a qualitative guide for our considerations. In order to evaluate eq 24, we calculated the classical free energy profiles for the systems studied while constraining the donor and acceptor distance at 3.0, 3.25, and 3.5 Å. The calculations were performed for a temperature of 300 K, and the profiles were shifted in a way that their minima would follow the PMF of Figure 7. Next, we evaluated the average rate constant of eq 24 while replacing the integral by the average of the contributions from R ) 3.0, 3.25, and 3.5 Å. In doing so, we approximated the ∆∆gqqm from R ) 3.25 Å by the average of the corresponding contributions from R ) 3.0 and 3.5 Å. The calculations were performed for both the H and D reactions at 285, 300, and 315 K. The classical activation barrier was fixed at its value at T ) 300 K, assuming that the temperature dependence of the barrier does not affect significantly the temperature dependence of the KIE (this restricted assumption is removed below). The resulting calculated temperature dependence of the KIEs appears to depend on the change in (∆ gqeff + w) between R)3.0 Å and R)3.5 Å. This key parameter j qeff could not is referred to here as ∆w j qeff. The exact value of ∆w be determined at the present stage due to the following

difficulties. First, the value of this parameter for the solution reaction (which is crucial for the EVB calibration procedure) could not be determined from the ab initio calculations of the potential surface of the solution reaction39 due to the error range of these calculations. The other problem is associated with the above-mentioned uncertainties in the calculated PMF. Thus, we considered this value as a parameter, which was adjusted to get the best fit for the observed temperature dependence of the KIE of the double mutant. However, the other calculated results, including the shape of the free energy profiles and the NQM corrections, were kept unchanged. At any rate, the relationship between the classical surfaces at R ) 3.0, 3.25, and 3.5 Å for the optimal ∆w j qeff are depicted in Figure 8. Finally, the temperature dependence of the KIEs obtained by our analysis is summarized in Figures 9 and 10, together with the corresponding experimental results (taken from refs 10, 58, and 59). Figure 9 gives the calculated results, assuming that the classical activation barrier has the same temperature dependence at any R. As seen from the figure, we obtain a reasonable trend, although for the double mutant, the calculated KIE occurs at a higher average value than the corresponding experimental results. To further explore the origin of the observed effect, we focused on the double mutant, where we have the largest deviation between the calculated and observed results. Here, we explored first the effect of the temperature dependence at a fixed R by using constant KIE for all temperatures. The corresponding result (Figure 10a) indicates that the temperature dependence at a fixed R is likely to play a major role in the case of the double mutant; without this dependence, we obtain an almost constant KIE at all temperatures, despite the effect of the PMF. In other words, it is hard to reproduce the temperature dependence of the KIE by using a reasonable PMF and quantum correction obtained by the vibronic treatment (which does not provide a temperature dependence at a fixed R). Next, we explored the additional possible effect of the change in activation entropy upon change in distance. As seen clear from ref 10, we have a really large change in -T∆S# between the native and double mutant cases. This difference is likely to reflect free energy changes in much a smaller temperature range than 0 to 300 K (which is assumed in deriving ∆S#), and it can

7858 J. Phys. Chem. B, Vol. 111, No. 27, 2007

Figure 8. The classical free energy profiles for the hydride transfer when the donor-acceptor distance is 3.0 Å (in black), 3.25 Å (in red), and 3.5 Å (in green) for the native DHFR (a), the G121V mutant (b) and the M42W-G121V mutant (c). The calculations of the classical surfaces were performed at T ) 300 K. The figure represents the relative behavior of the profiles for the optimal values of ∆w j qeff.

be used to give a rough estimate of the change of T∆S# between R ) 3.0 Å and R ) 3.5 Å (a much more rigorous estimate requires extensive simulation of the type described in ref 62 and is left to subsequent studies). Assuming that T∆S# increases by ∼2 kcal/mol between R ) 3.0 Å and R ) 3.5 Å, we obtained the results depicted in Figure 10b. As seen from this Figure, we obtained a better agreement with the experimental results and also obtain a reasonable temperature dependence when the KIE at R ) 3.5 Å is temperature-independent. IV. Concluding Remarks This work explored the temperature dependence of the KIE in the native DHFR and two mutants by microscopic molecular simulations. The calculations reproduced in a reasonable way the overall trend in the temperature dependence thus provided a reasonable “first principle” calculation of the temperature dependence of KIE in enzymes. The reason for our relative success has been the averaging over different donor-acceptor distances and over different simulations. Before considering in detail the view that emerged from the present study about the temperature dependence of the KIE,

Liu and Warshel we would like to point out the fundamental difference between a parameter fitting to a phenomenological formula (ref 1) and actual microscopic simulations used here and in our previous studies. Although it is almost trivial to reproduce the observed effect by arbitrary fitting parameters, it is extremely hard to do so by microscopic simulations. A part of the difficulty is associated with the difficulty of obtaining the exact PMF, and a part is associated with the general difficulties of capturing temperature dependence effects (e.g., entropic effects) by current simulation approaches (see the discussion in ref 61). Here, however, one can use the comparison of the calculated and observed temperature dependences to refine the calculated PMF on the basis of experimental measurements. Doing this while using eq 11 or eq 20 may not be so useful, since these equations do not reflect a quantitative physical picture. The use of a full microscopic treatment in modeling the KIE and its temperature dependence has reduced significantly the uncertainties associated with the use of arbitrary parameters and, more importantly, the problems associated with the use of the vibronic treatment. With this in mind, we explored several general limits that were not considered in previous studies. First, it was found that the temperature dependence of the KIE of the double mutant cannot be reproduced with temperature independence KIE at each distance and with a temperature independent of the w# defined in eq 26. Here, we note that that the vibronic treatment that is the basis of many of the studies in the field (e.g., see refs 1, 20, 63) assumes temperature independence at each R. Of course, it is still possible to reproduce the observed KIE with these conditions using unrealistic fast dependence of the KIE on the distance (i.e., eq 20) that would be produced by the vibronic treatment and an arbitrary w#(r). However, the current treatment is aimed at reducing the arbitrariness of the analysis. Thus, it is likely that some of the temperature dependence of the KIE reflects the dependence of ∆∆gqqm on the temperature at fixed R. We also found that if the activation entropy depends on the distance, we can reproduce the temperature dependence of the KIE even with a fixed KIE at each R. The most important finding is probably the establishment of a clear correlation between the increase in the average donoracceptor distance and the isotope effect. This finding, which is consistent with the qualitative prediction of the vibronic treatment, is apparently opposite to the traditional expectation from NQM corrections. That is, the traditional attribution for NQM contribution to catalysis23-29,64 implied that the protein compresses the reaction fragments and leads to a narrow potential and, thus, to larger tunneling. (Figure 11a). In this case, one would expect a relationship of the form (KIE)tun ∝ exp{-KxV0L}

(27)

On the other hand, it has been found here that the isotope effect increases due to the sharp distance dependent on the zerozero vibrational overlap (Figure 11b). Where the tunneling contribution follows the relationship (see eqs 14-15), (KIE)tun ∝ exp{K′L2(x2 - 1)}

(28)

Thus, in the more recent view, the NQM effects decrease rather than increase upon compression. This has, of course, a major implication with regard to the idea that NQM effects help in enzyme catalysis. In fact, the effects that lead to an increase in the NQM contributions appear to be anticatalytic (the rate constant is smaller for larger donor-acceptor distance).

Isotope Effects in Enzymatic Reactions

J. Phys. Chem. B, Vol. 111, No. 27, 2007 7859

Figure 9. The calculated (in black) and observed (in red) temperature dependence of the KIE for the native DHFR (a), the G121V mutant (b), and the M42W-G121V mutant (c).

Figure 10. The calculated temperature dependence of the KIE of the M42W-G121V mutant for the following limits: (a) the case when ∆gqcl (R) is assumed to be temperature-independent: while the KIE at each R is temperature-independent (9) and while the KIE at each R is temperature-dependent (b). The experimental results are designated by 2. (b) The same as in (a), but with distant dependent -T∆S# (see text).

It should be clarified here that the prediction of the vibronic formula about the relationship between the KIE and the donoracceptor distance has been found before our recent studies (e.g., ref 20, 65). However, with the exception of the work of Hammes-Schiffer,20 we could not find an early realization that the increase in distance of the KIE leads to a picture that is different from the traditional picture. Furthermore, we are not aware of any study (except our own; e.g., ref 18) that actually established the distance dependence of the KIE by proper microscopic simulations without the use of the problematic diabatic approximation and the vibronic formula. Moreover, we could not find a realization that this finding makes it hard to attribute catalysis to the increase in tunneling (see below). Although the qualitative prediction of the vibronic treatment is confirmed by the quantitative path integral simulations, it is important to emphasize that the vibronic treatment may be a very poor approximation in the range of regular donor-acceptor distance, which should be described by the adiabatic approximations (see Section II.3). In fact, the predicted distance depen-

dence of the isotope effect on R is much smaller in the actual path integral simulation than in eq 16. Nevertheless, at long distances when the electronic coupling is small, the vibronic formula provides a reasonable approximation. Thus, the overall trend at a decrease in KIE upon a decrease in distance follows the prediction of the vibronic formula. As stated in Section III, our treatment involves the use of ∆w j qeff as an adjustable parameter because of uncertainties in the calculated ab initio potential surface for the solution reaction and the uncertainties in the PMF. In this respect, one may wonder why such treatment was not needed in the recent treatment of the temperature dependence of the KIE of lypoxygenase.63 This success might reflect the fact that the authors of ref 63 used the phenomenological vibronic expression (although this was done with microscopically derived parameters) and did not explore the performance of such a model on different mutants. Furthermore, they used an average donor-acceptor distance (the equivalent of the averaging procedure of eq 20) and determined this distance by simulations. Although this procedure is very reasonable, it is unlikely that the average distance can be determined with certainty by simulations in cases of large changes in the PMFs (i.e., the double mutant case), since it will be subjected to the same major convergence problems encountered in our PMF calculations. It is sometimes stated that the HT reaction in DHFR and other enzymes is dominated by nuclear tunneling (e.g., ref 19). This view is based on the fact that the NQM contribution to the rate constant is more than 50%. Although this view may be technically reasonable, it can lead to a major confusion. That is, the actual NQM correction of the rate constant is only 1 kcal/mol, as compared to the 15 kcal/mol total barrier. This means that transition state theory provides a very good description of the rate constant and an excellent description as much as catalysis concerned (the NQM corrections are similar in the enzyme and in solution). When the donor and acceptor distance increases, the relative NQM effects increase, but the overall rate constant becomes very small and does not contribute too much to the average rate constant of eq 24. Moreover, the idea that the reaction is dominated by tunneling may give the impression that we have two independent channels, one of which is thought to be the barrier and the other is over the barrier, and that the detailed understanding of both channels is crucial for understanding catalysis. In fact, even if the main contribution to the rate constant were due to contributions from large Rs (in which the diabatic limit is valid), this contribution would involve a relatively large classical barrier (the ∆gq00 of eq 13 includes about 5 kcal/mol contribution from the solvent reorganization and additional contribution from the solute classical modes). However, the fact that the tunneling contribution is relatively small in DHFR and other findings of the present study indicate that the main contribution to the rate constant come from

7860 J. Phys. Chem. B, Vol. 111, No. 27, 2007

Liu and Warshel to the NAC configuration (e.g., ref 23)) reduce the tunneling rather than increase it. Similarly, the concept of environmentally coupled quantum mechanical tunneling, which is presumably linked to catalysis,66,67 is also problematic in its implications. That is, if the issue is the role of the modes orthogonal to the donor and acceptor distance (or the protein reorganization energy), then we are not dealing with anything that is unique to enzymes, and in fact, the coupling to the environment is significantly smaller in the enzyme than in solution (e.g., ref 17). If the issue is the donor-acceptor distance, then we have a similar distance in the native enzyme and the reference reaction in water and similar motional effects. The current use of the term “Marcus-like” model (e.g., refs 10, 65) for the description of the vibronic formula may be somewhat confusing as much as enzyme catalysis is concerned. The most serious problem is the possible implication that the reorganization energy term is closely related to tunneling effects. In fact, the reduction of the reorganization energy by the preorganized enzyme active site is probably the most important factor in enzyme catalysis,68,69 and it has enormous effects in completely adiabatic reactions without any significant tunneling effect. This work and our previous studies (e.g., ref 17) argued that the temperature dependent on the KIE does not provide a support to the proposal that dynamical effects contribute to catalysis. However, this does not mean that studies of the temperature dependent on KIE are not very valuable. In particular, by “inverting” the observed isotope effect obtained from the rate constants of eq 24 and by combining theory with experiment, we may use temperature-dependent experiments of the type reported in ref 10 to construct experimentally determined PMF for the native and mutant enzymes. Such an approach can clearly help in refining simulation approaches and in elucidating the effects of different mutations on the donor-acceptor distance.

Figure 11. The catalytic (anticatalytic effect) of the donor-acceptor distance dependence of the NQM in hydride transfer reactions. The traditional idea that the NQM increase upon compression of the donoracceptor distance is described in Figure 11a. This view assumes that enzyme catalysis involves an increase in tunneling due to the compression effects. On the other hand, as shown here, the NQM effects increase when the donor-acceptor distance decreases (Figure 11b), and thus, it is hard to rationalize a relationship between NQM and catalysis.

intermediate distances where the vibronic treatment is not valid. In this case, one has to look at the NQM effect as a correction to the activation barrier. Thus, it is not so useful to argue that if the NQM correction reduces a barrier of, say, 15 kcal/mol by 1 kcal mol (and reduces the 1011 effect of the activation barrier by a trivial factor of 10) that NQM effects dominate the reaction. With this view, one would expect that the effect of mutations that change the rate constant by, say, 3 orders of magnitude will be due to tunneling (that presumably dominates the reaction). However, the change in rate constant is entirely due to the change in the classical activation barrier. It is also important to consider here the idea that DHFR increases the reaction rate “through optimization of the coupling of the environmental dynamics and quantum tunneling (see ref 66 or related ideas27) as well as the implication that quantum mechanical tunneling is promoted by correlated motions (e.g., ref 19) and is optimized by the wild type enzymes.10 In fact, the tunneling is smaller in the native enzyme (where the catalysis is maximal) than in the mutant enzyme, and this meant that the motions that move the reacting system in the mutant enzyme to the native configuration (which can be described as moving

Acknowledgment. This work was supported by NIH grant GM24492 and by computer time from the University of Southern High Performance Computing and Communication Center (HPCC). H. Liu thanks Dr. Mats Olsson for helpful discussions. References and Notes (1) Knapp, M. J.; Rickert, K.; Klinman, J. P. J. Am. Chem. Soc. 2002, 124, 3865. (2) Kemsley, J. Chem. Eng. News 2003, 81, 29. (3) Doll, K. M.; Bender, B. R.; Finke, R. G. J. Am. Chem. Soc. 2003, 125, 10877. (4) Wang, S. X.; Mure, M.; Medzihradszky, K. F.; Burlingame, A. L.; Brown, D. E.; Dooley, D. M.; Smith, A. J.; Kagan, H. M.; Klinman, J. P. Science 1996, 273, 1078. (5) Lehnert, N.; Solomon, E. I. J. Biol. Inorg. Chem. 2003, 8, 294. (6) Hwang, J.-K.; Warshel, A. J. Am. Chem. Soc. 1996, 118, 11745. (7) Billeter, S. R.; Webb, S. P.; Agarwal, P. K.; Iordanov, T.; HammesSchiffer, S. J. Am. Chem. Soc. 2001, 123, 11262. (8) Hwang, J.-K.; Warshel, A. J. Phys. Chem. 1993, 97, 10053. (9) Olsson, M. H. M.; Siegbahn, P. E. M.; Warshel, A. J. Am. Chem. Soc. 2004, 126, 2820. (10) Wang, L.; Goodey, N. M.; Benkovic, S. J.; Kohen, A. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 15753. (11) Kuznetsov, A. M.; Ulstrup, J. Can. J. Chem. 1999, 77, 1085. (12) German, E. D.; Kuznetsov, A. M. J. Chem. Soc. Faraday Trans. 1 1981, 77, 397. (13) Warshel, A.; Chu, Z. T. J. Chem. Phys. 1990, 93, 4003. (14) Hwang, J.-K.; Chu, Z. T.; Yadav, A.; Warshel, A. J. Phys. Chem. 1991, 95, 8445. (15) Villa, J.; Warshel, A. J. Phys. Chem. B 2001, 105, 7887. (16) Gao, J. L.; Truhlar, D. G. Annu. ReV. Phys. Chem. 2002, 53, 467. (17) Olsson, M. H. M.; Parson, W. W.; Warshel, A. Chem. ReV. 2006, 106, 1737. (18) Olsson, M. H. M.; Mavri, J.; Warshel, A. Philos. Trans. R. Soc. London, Ser. B 2006, 361, 1417.

Isotope Effects in Enzymatic Reactions (19) Pang, J. Y.; Pu, J. Z.; Gao, J. L.; Truhlar, D. G.; Allemann, R. K. J. Am. Chem. Soc. 2006, 128, 8015. (20) Hatcher, E.; Soudackov, A. V.; Hammes-Schiffer, S. J. Am. Chem. Soc. 2004, 126, 5763. (21) Villa, J.; Gonzalez-Lafont, A.; Lluch, J. M.; Truhlar, D. G. J. Am. Chem. Soc. 1998, 120, 5559. (22) Rodgers, J.; Femec, D. A.; Schowen, R. L. J. Am. Chem. Soc. 1982, 104, 3263. (23) Luo, J.; Kahn, K.; Bruice, T. C. Bioorg. Chem. 1999, 27, 289. (24) Ball, P. Nature 2004, 431, 396. (25) Bahnson, B. J.; Colby, T. D.; Chin, J. K.; Goldstein, B. M.; Klinman, J. P. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 12797. (26) Klinman, J. P. Trends Biochem. Sci. 1989, 14, 368. (27) Mincer, J. S.; Schwartz, S. D. J. Phys. Chem. B 2003, 107, 366. (28) Sutcliffe, M. J.; Scrutton, M. S. Philos. Trans. R. Soc. London, Ser. A 2000, 358, 367. (29) Schwartz, S. D. Vibrationally Enhanced Tunneling and Kinetic Isotope Effects in Enzymatic Reactions. In Isotope Effects in Chemistry and Biology; Kohen, A., Limbach, H. H., Eds.; CRC: Boca Raton, FL 2004; p 475. (30) Warshel, A. J. Biol. Chem. 1998, 273, 27035. (31) Ohta, Y.; Soudackov, A. V.; Hammes-Schiffer, S. J. Chem. Phys. 2006, 125, 144522. (32) Hwang, J.-K.; King, G.; Creighton, S.; Warshel, A. J. Am. Chem. Soc. 1988, 110, 5297. (33) Warshel, A. Computer Modeling of Chemical Reactions in Enzymes and Solutions; John Wiley & Sons: New York, 1991. (34) Aqvist, J.; Warshel, A. Chem. ReV. 1993, 93, 2523. (35) Chu, Z. T.; Villa, J.; Strajbl, M.; Schutz, C.; Shurki, A.; Warshel, A. Molaris; version 9.05, b.; Los Angeles 2002. (36) Lee, F. S.; Chu, Z. T.; Warshel, A. J. Comput. Chem. 1993, 14, 161. (37) King, G.; Warshel, A. J. Chem. Phys. 1989, 91, 3647. (38) Lee, F. S.; Warshel, A. J. Chem. Phys. 1992, 97, 3100. (39) Liu, H.; Warshel, A. Biochemistry 2007, 46, 6011. (40) Wang, M. L.; Lu, Z. Y.; Yang, W. T. J. Chem. Phys. 2006, 124, 124516. (41) Wang, Q.; Hammes-Schiffer, S. J. Chem. Phys. 2006, 125, 184102. (42) Major, D. T.; Gao, J. L. J. Am. Chem. Soc. 2006, 128, 16345. (43) Warshel, A.; Parson, W. W. Q. ReV. Biophys. 2001, 34, 563. (44) Feynman, R. P. Statistical Mechanics; Benjamin: New York, 1972.

J. Phys. Chem. B, Vol. 111, No. 27, 2007 7861 (45) Gillan, M. J. J. Phys. C: Solid-state Phys. 1987, 20, 3621. (46) Voth, G. A. AdV. Chem. Phys. 1996, 93, 135. (47) Voth, G. A.; Chandler, D.; Miller, W. H. J. Chem. Phys. 1989, 91, 7749. (48) Warshel, A. J. Phys. Chem. 1982, 86, 2218. (49) Warshel, A. Proc. Natl. Acad. Sci. U.S.A. 1984, 81, 444. (50) Warshel, A.; Hwang, J.-K. J. Chem. Phys. 1986, 84, 4938. (51) Borgis, D.; Hynes, J. T. J. Chem. Phys. 1991, 94, 3619. (52) Antoniou, D.; Schwartz, S. D. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 12360. (53) Manneback, C. Physica 1951, XVII, 1001. (54) Warshel, A. Annu. ReV. Biophys. Bioeng. 1977, 6, 273. (55) Suhnel, J.; Gustav, K. Chem. Phys. 1984, 87, 179. (56) German, E. D.; Kuznetsov, A. M. J. Chem. Soc. Faraday Trans. 2 1981, 77, 2203. (57) Mikenda, W. J. Mol. Struct. 1986, 147, 1. (58) Sikorski, R. S.; Wang, L.; Markham, K. A.; Rajagopalan, P. T. R.; Benkovic, S. J.; Kohen, A. J. Am. Chem. Soc. 2004, 126, 4778. (59) Wang, L.; Tharp, S.; Selzer, T.; Benkovic, S. J.; Kohen, A. Biochemistry 2006, 45, 1383. (60) Warshel, A.; Olsson, M. H. M.; Villa-Freixa, J. Computer Simulations of Isotope Effect in Enzyme Catalysis. In Isotope Effects in Chemistry and Biology; Kohen, A., Limbach, H. H., Eds.; CRC: Boca Raton, FL 2006. (61) Kato, M.; Warshel, A. J. Phys. Chem. B 2005, 109, 19516. (62) Villa, J.; Strajbl, M.; Glennon, T. M.; Sham, Y. Y.; Chu, Z. T.; Warshel, A. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 11899. (63) Hatcher, E.; Soudackov, A. V.; Hammes-Schiffer, S. J. Am. Chem. Soc. 2007, 129, 187. (64) Tejero, I.; Garcia-Viloca, M.; Gonzalez-Lafont, A.; Lluch, J. M.; York, D. M. J. Phys. Chem. B 2006, 110, 24708. (65) Knapp, M. J.; Klinman, J. P. Eur. J. Biochem. 2002, 269, 3113. (66) Swanwick, R. S.; Maglia, G.; Tey, L.; Allemann, R. K. Biochem. J. 2006, 394, 259. (67) Maglia, G.; Allemann, R. K. J. Am. Chem. Soc. 2003, 125, 13372. (68) Warshel, A. Proc. Natl. Acad. Sci. U.S.A. 1978, 75, 2558. (69) Warshel, A.; Sharma, P. K.; Kato, M.; Xiang, Y.; Liu, H.; Olsson, M. H. M. Chem ReV 2006, 106, 3210.