Simulation of Tunneling in Enzyme Catalysis by Combining a Biased

Dec 11, 2007 - molecular dynamics simulation of SLO-1. Since the propagations require an extremely long time for treating realistic activation barrier...
0 downloads 0 Views 108KB Size
Download PDF
5950

J. Phys. Chem. B 2008, 112, 5950-5954

Simulation of Tunneling in Enzyme Catalysis by Combining a Biased Propagation Approach and the Quantum Classical Path Method: Application to Lipoxygenase† Janez Mavri,‡,§ Hanbin Liu,§ Mats H. M. Olsson,§ and Arieh Warshel*,§ National Institute of Chemistry, HajdrihoVa 19, 1000 Ljubljana, SloVenia, and Department of Chemistry, UniVersity of Southern California, 3620 South McClintock AVenue, Los Angeles, California 90089-1062 ReceiVed: July 25, 2007; In Final Form: October 8, 2007

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

1. Introduction Elucidating the origin of the large catalytic power of enzymes is a major challenge for computational chemistry. Our previous simulation studies gave strong evidence that the main contribution comes from preorganized electrostatics associated with the protein environment.1 However, the electrostatic preorganization idea has not been uniformly accepted and alternative proposals for enzymatic catalysis are still invoked frequently (for reviews, see, e.g., refs 1-4). One of the prominent proposals involves the idea that the transition state theory (TST) is not valid and that deviations from it contribute to large catalytic effects.5-7 It is also argued frequently that nuclear quantum mechanical effects (NQM) lead to major catalytic contributions.8,9 In our view, the validity of TST in assessing the major contributions to enzyme catalysis has been established repeatedly.1,2,10 Furthermore, many studies that actually compared NQM in enzyme and in solution (and thus examined catalysis) found that NQM do not contribute to catalysis (see review in ref 2). For example, our detailed calculations of the isotope effect of lipoxygenase and related systems found similar NQM in the protein and in solution.1,2,11-13 Nevertheless, it is important to quantify the role of NQM in enzymatic reactions. The quantum classical path (QCP)14-18 has been used in most studies of the catalytic contributions from NQM (the difference between NQM in the enzyme and in water; see ref 2 for a review), but other approaches have been found to provide powerful ways of exploring NQM in enzyme and in solution and of calculating isotope effects in enzymes.19-21 Despite the advances of several effective methods for such calculations, †

Part of the “Attila Szabo Festschrift”. * Corresponding author. E-mail: [email protected]. ‡ National Institute of Chemistry. § University of Southern California.

there were only a few studies that compared the performance of the different methods (e.g., ref 22). In this work, we try to explore the ability of the wave function propagation method to evaluate NQM in lipoxygenase. More specifically, we generate one-dimensional sections of the multidimensional reaction surface in the fluctuating enzyme, propagate the wave functions of H or D by the approach of Kosloff and Kosloff,23 and calculate the rate constant quantum-dynamically from the timedependent wave function. Since the enzyme works on a millisecond time scale, which is far longer than the practically feasible propagation time, we use a biasing potential and evaluate the effect of this potential on the quantum mechanical rate constant. This is done by a new formulation which is based on the QCP approach. Overall, we do not expect the new approach to be particularly useful relative to our standard QCP approach, but it gives us a way to estimate the quantum mechanical transmission factor, which is not provided by the direct probabilistic centroid approach. Our work takes lipoxygenases as a benchmark for the biased propagation approach. This system consists of a structurally related family of non-heme iron containing dioxygenases that catalyze the addition of molecular oxygen to unsaturated fatty acids. The rate determining step of the catalytic reaction of this enzyme involves a C-H bond cleavage, as depicted in Figure 1. Here, we focus on the reaction of soybean lipoxygenase-1 (SLO-1). The nature of the NQM effects of SLO-1 has been explored extensively by recent staudies.5,6,8,11,12,24-26 It was suggested that the enzyme catalyzes its reaction by enhancing tunneling by coupling the enzyme vibrational modes to the CH stretching coordinate and therefore enhancing the reaction rate.5,8 Such proposals have been supported by oversimplified models that have not explored the actual catalytic effect.6 However, our detailed QCP studies11,12 have demonstrated that the tunneling

10.1021/jp0758420 CCC: $40.75 © 2008 American Chemical Society Published on Web 12/11/2007

Simulation of Tunneling in Enzyme Catalysis

J. Phys. Chem. B, Vol. 112, No. 19, 2008 5951

Figure 1. Schematic description of the rate determining step in the reaction of SLO-1.

contributions are similar in the enzyme and in solution and thus do not contribute to catalysis. The catalytic reaction of SLO-1 was simulated by other workers and, in particular, by HammesSchiffer and co-workers, who provided an interesting insight by using a simplified model that represented the enzyme by a continuums model24 and subsequently included the enzyme explicitly in the simulation.26 An additional study was reported by Garcia-Viloca and co-workers25 who calculated the H/D kinetic isotope effect (KIE) in lipoxygenase by considering the full protein environment, while using a quantum mechanic/ molecular mechanic (QM/MM) approach in conjunction with the variational transition state theory. The classical activation energy obtained by the applied semiempirical MO method was too low leading to underestimated H/D KIEs. This problem was then corrected by introducing a large shift in the original surface in order to reproduce the observed isotope effect, without reproducing, however, the observed rate constant. While this approach is reasonable in some respects, it is important to clarify that, in contrast to the suggestion of ref 25, this approach is quite different than the approach used by Olsson et al.11 That is, the QM/MM empirical valence bond (EVB) surface of ref 11 was calibrated by fitting it to ab initio results and then by introducing a minor refinement in order to reproduce the observed kH, without any adjustment or parametrization that would reproduce the observed KIE. In other words, the systematic EVB approach of ref 11 reproduced quantitatively the observed KIE result without using any prior information about this effect. We believe that obtaining the large observed KIE of SLO-1 without biasing the surface to reproduce this effect is an important step in gaining confidence about the conclusion that emerges from the study of ref 11 (see also concluding remarks). At any rate, the present work is not about the QCP/EVB studies of SLO-1, where we already obtained extremely encouraging results.11 In this paper, we attempt to use a very different approach of direct propagation plus a biasing potential (and a specialized QCP-based evaluation of the effect of this potential) to obtain reasonable KIE results. 2. Calculating the Fluctuating One-Dimensional Hydrogen Potential In order to use a propagation method in studies of hydrogen transfer in a multidimensional system, it is essential to capture the time dependent fluctuations of the hydrogen transfer potential because of its coupling to the rest of the system. Thus, we used EVB model to simulate the time dependence of the entire enzyme substrate complex and used the simulated results to evaluate the fluctuations of the one-dimensional hydrogen transfer potential. This was done by the following strategy: (i) the reactive system was modeled by the EVB QM/MM method, while including in the quantum subsystem the methylene group, the OH group, and the iron ion by a two state model described in ref 11. The rest of the classical environment was incorporated into the quantum region according to the standard EVB procedure.10,29 (ii) EVB MD runs were propagated on the reactant state of the reacting system for 20 ps restraining the CHO angle to 180° with a force constant of 10 kcal/(mol rad2).

Figure 2. Schematic representation of the evaluation of the 1-D hydrogen potentials. The calculations involve running trajectories on the reactant state of the EVB surface and “cutting” the surface along the hydrogen transfer vector in increments of equal time step.

Every 20 fs, we “cut” a one-dimensional potential for a transfer of the hydrogen from the reactant to product state along the line connecting the carbon and oxygen atoms. These potentials are described schematically in Figure 2 (see also Figure 2 of ref 32). In the present case, we accumulated 1000 potentials. The initial CH distance for the calculation of the hydrogen potential was 0.6 Å, and its final value was set to 2.40 Å. For each CH distance, all other degrees of freedom were kept fixed, implying that H or D asymmetric stretching dynamics is fast relative to the other degrees of freedom. This corresponds to Marcus’ picture of the electron-transfer reactions where the electrons feel frozen nuclear degrees of freedom. Each hydrogen potential was evaluated at 64 points along the C‚‚‚O vector for 64 CH distance values. The hydrogen potentials between the points were calculated using a spline interpolation procedure. The series of hydrogen potentials evaluated by the above approach were subjected to the quantum calculations described below, and the corresponding results were averages to give the final rate constants. The simulation system was described by the surface constrained all atoms solvent (SCAAS)30 model, and the local reaction field (LRF)31 model was used to handle the long-range effect. The molecular dynamics (MD) simulations were done with a 0.5 fs time step at 300 K. More details about the specific EVB treatment used can be found in ref 11. All atom simulations of the essential reaction step of lipoxygenase were performed with the MOLARIS suite of programs.27,28 3. Propagating the H Wave Function The grid propagator method of Kosloff and Kosloff23 was used for solving the time dependent Schro¨dinger equation for the hydrogen potentials described above. This method is based on using a Fourier transform of the wave function in the calculation of the kinetic energy term. The corresponding formulation is robust and can be coded conveniently. The hydrogen potentials were interpolated as a function of time by using cubical splines. The initial wave function was Gaussian shaped, centered in the reactant well. The initial width of the wave function did not have a significant effect on the calculated rate constant. The dividing surface between the reactants and the products was set at 1.36 Å from the carbon atom. One hundred random starting points were used in each calculation,

5952 J. Phys. Chem. B, Vol. 112, No. 19, 2008

Mavri et al.

TABLE 1: H/D Kinetic Isotope Effect Evaluated with Different Biasing Potentials and Corrected by QCP Calculationsa Al (kcal/mol)

KIEpropl

NQM correction of bias potential

70 65 60 55

7(1 15 ( 2 22 ( 10 42 ( 26

1.59

a

KIEg

66.8 l

The propagation approach was used for all values of A , and the QCP correction term was only evaluated for the smallest Al.

and the time-dependent Schro¨dinger equation was integrated in the forward direction with a time step of 0.001 fs for 5 000 000 steps giving rise to 5 ps of real time with a grid of 64 points (vide infra). The time-dependent population of the reactant well was averaged over 100 independent runs. The rate constant was calculated from the initial decay of the population of the reactant well. For the fast Fourier transform (FFT) scheme, we examined several numbers of grid points of 512, 256, 128, and 64. The latter number of grid points was sufficient to achieve the convergence in terms of the rate constant while simultaneously keeping the low spectral range of the propagator, allowing for a relatively large time step. Our preliminary calculations indicated that it is very hard to obtain meaningful quantum results for the unbiased EVB surface due to the need for an extremely long integration time. Thus, we tried to explore the effect of a biased sampling by elevating the reactant well with a Gaussian wave function u ) A exp(-2.77(1.13 - x)2), where x is the C-H distance, 1.13 Å represents the center of the reactant well, and A is the preexponential factor. Different values of pre-exponential factor A were applied. A higher value of A means a lower effective activation energy and higher rate constant. Four different values of A were used, giving rise to different values of the classical activation energy. The results are summarized in Table 1. The H/D KIE for the given biasing potential was calculated as a ratio of the rate constant for reaction with H and D. Each rate constant was calculated by linear regression of the reactant

population. The uncertainty was estimated as the absolute value of the difference between the slope from linear regression and the slope estimated by visual examination. A typical time averaged course of the population of the reactant well is shown in Figure 3. Decreasing the value of A resulted in the increase of the effective barrier and a larger KIE. However, this led to a decrease of the rate constants for both H and D. However, the uncertainty of the calculated KIE also decreased, due mainly to the low value of the rat constant for D. 4. Sampling the Propagation Results Using the QCP Approach Considering the difficulties in using direct simulations, we looked for a way to use the direct simulations with the biasing potential in evaluation of the actual rate constant on the unbiased reaction surface. In searching for an optimal approach, we felt that the application of advanced propagators33-35 will not solve the problem of biased sampling in the case of large barriers. We note in this respect that biased sampling was applied to simulate much smaller H/D KIE in the case of carboxylic ester hydrolysis,36 but this approach failed to consider the crucial quantum contribution to the biasing free energy (see below). Thus, we looked for a different strategy. The starting point of our derivation is the observation that the centroid path integral approach gives a quantized rate constant in the form37,38

kg ) κqg(kBT/h) exp{-∆g*(Ug)qβ}

(1)

where q designates quantum mechanics, κqg, kB, T, h, and β are respectively the transmission factor on the ground state surface, Boltzmann’s constant, the temperature, Planck’s constant, and β ) 1/kBT. The quantum mechanical activation free energy, ∆gqq, includes almost all of the nuclear quantum mechanical effects, whereas only small effects come from the preexponential transmission factor in the case of systems with a

Figure 3. Time-averaged population of the reactant well for H and D, for calculations that involve the biasing potential ulb ) 55 exp(-2.77(1.1 - x)2).

Simulation of Tunneling in Enzyme Catalysis

J. Phys. Chem. B, Vol. 112, No. 19, 2008 5953 The term exp{-β(g(Ua)x0,q,H - g(Ua)x0,q,D)} can be expressed as

exp{-β(g((Ua)x0,q,H - g(Ua)x0,cl) - (g(Ua)x0,q,D g(Ua)x0,cl))} where cl designates the classical potential surface. Now, we can write

exp{-β[g(Ua)x0,q,R - g(Ua)x0,cl]} ) Figure 4. Illustrating the effect of the biasing potential, ulb, on the reaction surface.

significant activation barrier,20,39 when the NQM probabilistic effects are included in ∆gqq. Since κqg depends mainly on the behavior at the TS region, we will assume that κqg ≈ κql in cases when the potentials Ul and Ug (see Figure 4) are similar at the TS region. Considering the addition of the biasing potential, uib, to the reaction surface (Figure 4), we can write

∆g* (Ug)q ) g(Ug)x*,q - g(Ug)x0,q

(2)

= (g(Ul)x*,q - g(Ul)x0,q) - (g(Ul)x*,q - g(Ug)x*,q) + (g(Ul)x0,q - g(Ug)x0,q) ) ∆g* (Ul)q - ∆g(Ug f Ul)x*,q + ∆g(Ug f Ul)x0,q (3) where Ul ) Ug + ulb, ulb ) Al exp{-2.77(1.13 - Rdonor-H)2}, and where we assume that x0g ≈ x0l . Here, ∆g(Ua f Ub)xi designates the change of the given ∆g at xi when Ua changes to Ub. Before evaluating the above ∆g’s, we note that eq 3 can be used to write

kg = κqg(kBT/h) exp{- ∆g* (Ul)qβ} × exp{- β(∆g(Ug f Ul)x0,q - ∆g(Ug f Ul)x*,q)} (4) Using this approximated relationship and assuming that κg ≈ κl, we can write

kg = k(Ul)q exp{-β(∆g(Ug f Ul)x0,q - ∆g(Ug f Ul)x*,q)} (5) Now if ulb(x* ) f 0, then the correction contribution from x# becomes small and the main contribution comes from x0. This situation is satisfied when Al becomes small, and in this case, we can write

kg = kprop(Ul) exp{-∆g(Ug f Ul)x0,q β}

(6)

Here, kprop(Ul) is the rate constant obtained from the direct propagation with the given Ul. The term exp{-∆g(Ug f Ul)x0,qβ} can be evaluated reliably and effectively by the QCP approach14,15 and is likely to account for most of the quantum effects associated with the change of Ul to Ug and thus for the corresponding correction of the isotope effect. Now our task is to determine the KIE and the ∆g’s correction terms.

KIEg ≈ KIEprop exp{-β(∆g(Ug f Ul)x0,q,H ∆g(Ug f Ul)x0,q,D)} ≈ KIEprop exp{-β(g(Ul)x0,q,H g(Ug)x0,q,H - g(Ul)x0,q,D + g(Ug)x0,q,D)} (7)

a - Ua(xc)])〉fp〉U 〈〈δ(x - x′) exp(-β[Uq,R

cl

a

(8)

where R is H or D, Uq,R is the nuclear quantum mechanic effective potential (see ref 14), and a is g or l. The expression 〈‚‚‚〉fp is the free particle quantum average, 〈‚‚‚〉Ucla is an average over the classical potential, and xc is the coordinate of the classical trajectory. See ref 14 for details. Thus, we calculated our corrections using eqs 7 and 8. The above approach has been used to convert the results of the propagation with a biased potential to the results that would be obtained without a bias. The calculations were done for the case with the smallest Al where the approximation of eq 10 is expected to be valid. As seen from the results reported in Table 1, the calculated KIE increases from the value obtained with Al ) 55 (KIEpropl ) 42) to about 67 with the QCP correction. This value agrees reasonably well with the observed value (KIEobs ≈ 80). 5. Concluding Remarks The study of NQM effects in complex systems presents a significant challenge for computational chemistry. Although there are several approaches that provide a reasonable way of obtaining such results, it is important to explore alternative models. The present work examines the ability of propagation methods to reproduce NQM in enzymes. We start our study by generating one-dimensional hydrogen potentials from EVB simulations of lipoxygenase and propagating the wave function of the hydrogen or the deuterium. It is found that in the case of high barriers (which is the case in lipoxygenase) it is important to apply a biasing potential. In view of this problem, it is essential to use a combination of direct propagation and a biasing potential and then evaluate of the effect of changing the potential from the actual potential to the biasing potential by the QCP approach. The present approach reproduced in a reasonable way the large NQM effect in SLO-1, thus adding further support to the validity of our EVB potential for this system. However, it is important to clarify that our most systematic studies of SLO-1 where done in refs 11 and 12 with the EVB/QCP approach (rather than with the propagation approach) and that the present work is not our deepest study of the relationship between tunneling and catalysis. In fact, since this issue continues to involve significant misunderstandings, we will re-emphasis here a few key issues. The first point is that our studies2 repeatedly show that, despite the exciting finding of large NQM in enzymes, these effects do not contribute to enzyme catalysis. That is, as demonstrated in ref 2, any attempt to actually simulate the NQM both in an enzyme and in the relevant reference solution reaction gave very similar effects. This indicated that nuclear tunneling effects do not contribute to catalysis. Unfortunately most studies that promoted the importance of NQM in catalysis have not examined the reference solution reaction and thus could not validate our findings. Second, in contrast to a long held belief that enzymes can “compress” the reacting

5954 J. Phys. Chem. B, Vol. 112, No. 19, 2008 system while making the potential surface narrower and increasing the NQM effects, we found that compressions actually reduce the KIE rather than increase it. Thus, the observation of large KIE reflects the increase (rather than decrease) in the distance between the donor and the acceptor. This seemingly counter-intuitive point is clarified and demonstrated in ref 40, and we will only consider a recent misunderstanding of this issue here. More specifically, a recent interesting study of SLO-125 explored the effect of changing the width of the gas-phase potential on the calculated KIE in a model of SLO-1. It was concluded that “the decompression of the active site will result in a lower contribution to tunneling ... thus indicating the importance of the proposed gating promoting modes”. However, this idea is very problematic since, as we found recently, decompression (or increasing the donor acceptor distance) increases, rather than decreases, the tunneling contribution.2 Unfortunately, the interesting analysis of ref 25 involved changes of the width of the intrinsic potential surface (the one of the isolated substrate without the enzyme) instead of using a single potential surface for the substrate and changing the donor-acceptor distance in the presence of the enzyme (as was done in ref 40 and 12) or actually reproducing the observed effect of mutations that change the donor-acceptor distance (as was done in ref 40). Furthermore, an arbitrary change of the width of the intrinsic potential does not address or represent the proper physics of hydrogen transfer processes since it does not converge to the proper diabatic limit at large distances. More importantly, it does not explore the effect of the protein since the active site strain is not likely to cause a significant narrowing of the potential surface of the reacting atoms (proteins are far too flexible to exert such an effect). We bring this discussion here not in order to criticize the work of ref 25, but rather to point out the risk of following widely held proposals by theoretical studies that model the experimentalist view of such proposals, instead of modeling the given experimental results (e.g., the observed effect of mutations). Although the present approach may provide an interesting way of simulating NQM effects, and despite the availability of other powerful options, we believe that the QCP method (and not the direct propagation approach) is probably the method of choice when it comes to simulations of rate constants of enzymatic reactions with large nuclear quantum effects.11 This is true, in particular, when one is interested in consistent microscopic simulations of the contributions of enzyme fluctuations to NQM effects. The QCP13,41 offers an effective and simple way to evaluate the probability of being at the transition state and at the reactant well and thus provides the quantized activation free energy. In contrast to the difficulties of using propagators on a grid, it is trivial to quantize a large number of degrees of freedom but one. For example, our QCP study of LSO-111,12 quantized the motions of five atoms (the methylene group and OH), which is equivalent to 15 degrees of freedom (larger number of degrees of freedom can be easily accommodated). This would be practically impossible to do with the use of numerical propagators on a grid. Despite the apparent advantage of the QCP over the propagation method, we find the current method to be potentially useful in cases when one is interested in quantum dynamical effects and in the transmission factor. We also note that the current QCP-based biasing approach may be useful in combination with other time-dependent approaches.

Mavri et al. Acknowledgment. This work was supported by NIH Grant GM24492. One of us (J.M.) thanks the J. William Fulbright Scholarship Board for the award of a research scholarship. Financial support from Slovenian Ministry of Higher Education, Science and Technology (Grant P1-0012) is gratefully acknowledged. We thank the University of Southern California High Performance Computing and Communication Center (HPCC) for computer time. References and Notes (1) Warshel, A.; Sharma, P. K.; Kato, M.; Xiang, Y.; Liu, H.; Olsson, M. H. M. Chem. ReV. 2006, 106, 3210. (2) Olsson, M. H. M.; Parson, W. W.; Warshel, A. Chem. ReV. 2006, 106, 1737. (3) Field, M. J. J. Comput. Chem. 2002, 23, 48. (4) Mulholland, A. J.; Grant, G. H.; Richards, W. G. Protein Eng. 1993, 6, 133. (5) Kohen, A.; Klinman, J. P. J. Am. Chem. Soc. 2000, 122, 10738. (6) Nunez, S.; Antoniou, D.; Schramm, V. L.; Schwartz, S. D. J. Am. Chem. Soc. 2004, 126, 15720. (7) Scrutton, N. S.; Basran, J.; Sutcliffe, M. J. Eur. J. Biochem. 1999, 264, 666. (8) Knapp, M. J.; Rickert, K.; Klinman, J. P. J. Am. Chem. Soc. 2002, 124, 3865. (9) Ball, P. Nature 2004, 431, 396. (10) Warshel, A. Computer Modeling of Chemical Reactions in Enzymes and Solutions, 1997 ed.; John Wiley & Sons: New York, 1991. (11) Olsson, M. H. M.; Siegbahn, P. E. M.; Warshel, A. J. Am. Chem. Soc. 2004, 126, 2820. (12) Olsson, M. H. M.; Mavri, J.; Warshel, A. Philos. Trans. R. Soc. London, Ser. B 2006, 361, 1417. (13) Hwang, J.-K.; Warshel, A. J. Am. Chem. Soc. 1996, 118, 11745. (14) Hwang, J.-K.; Warshel, A. J. Phys. Chem. 1993, 97, 10053. (15) Hwang, J.-K.; Chu, Z. T.; Yadav, A.; Warshel, A. J. Phys. Chem. 1991, 95, 8445. (16) Wang, M. L.; Lu, Z. Y.; Yang, W. T. J. Chem. Phys. 2006, 124, 124516. (17) Wang, Q.; Hammes-Schiffer, S. J. Chem. Phys. 2006, 125, 184102. (18) Major, D. T.; Gao, J. L. J. Mol. Graphics 2005, 24, 121. (19) Gao, J. L.; Truhlar, D. G. Annu. ReV. Phys. Chem. 2002, 53, 467. (20) Billeter, S. R.; Webb, S. P.; Agarwal, P. K.; Iordanov, T.; HammesSchiffer, S. J. Am. Chem. Soc. 2001, 123, 11262. (21) Kuznetsov, A. M.; Ulstrup, J. Can. Chem. J. 1999, 77, 1085. (22) Yamamoto, T.; Miller, W. H. J. Chem. Phys. 2005, 122. (23) Kosloff, D.; Kosloff, R. J. Comput. Phys. 1983, 52, 35. (24) Hatcher, E.; Soudackov, A. V.; Hammes-Schiffer, S. J. Am. Chem. Soc. 2004, 126, 5763. (25) Tejero, I.; Garcia-Viloca, M.; Gonzalez-Lafont, A.; Lluch, J. M.; York, D. M. J. Phys. Chem. B 2006, 110, 24708. (26) Hatcher, E.; Soudackov, A. V.; Hammes-Schiffer, S. J. Am. Chem. Soc. 2007, 129, 187. (27) Lee, F. S.; Chu, Z. T.; Warshel, A. J. Comput. Chem. 1993, 14, 161. (28) Chu, Z. T.; Villa, J.; Strajbl, M.; Schutz, C.; Shurki, A.; Warshel, A. Molaris, 9.05, b; Los Angeles, 2002. (29) Aqvist, J.; Warshel, A. Chem. ReV. 1993, 93, 2523. (30) King, G.; Warshel, A. J. Chem. Phys. 1989, 91, 3647. (31) Lee, F. S.; Warshel, A. J. Chem. Phys. 1992, 97, 3100. (32) Villa, J.; Warshel, A. J. Phys. Chem. B 2001, 105, 7887. (33) Berendsen, H. J. C.; Mavri, J. J. Phys. Chem. 1993, 97, 13464. (34) Billeter, S. R.; van Gunsteren, W. F. Mol. Simul. 1995, 15, 301. (35) Iyengar, S. S.; Jakowski, J. J. Chem. Phys. 2005, 122. (36) Lensink, M. F.; Mavri, J.; Berendsen, H. J. C. J. Comput. Chem. 1999, 20, 886. (37) Gillan, M. J. J. Phys. C. Solid State Phys. 1987, 20, 3621. (38) Voth, G. A.; Chandler, D.; Miller, W. H. J. Chem. Phys. 1989, 91, 7749. (39) Warshel, A.; Parson, W. W. Q. ReV. Biophys. 2001, 34, 563. (40) Liu, H. B.; Warshel, A. J. Phys. Chem. B 2007, 111, 7852. (41) Feierberg, I.; Luzhkov, V.; Aqvist, J. J. Biol. Chem. 2000, 275, 22657.