## Abstract

The idea that enzyme catalysis involves special factors such as coherent fluctuations, quantum mechanical tunnelling and non-equilibrium solvation (NES) effects has gained popularity in recent years. It has also been suggested that transition state theory (TST) cannot be used in studies of enzyme catalysis. The present work uses reliable state of the art simulation approaches to examine the above ideas. We start by demonstrating that we are able to simulate any of the present catalytic proposals using the empirical valence bond (EVB) potential energy surfaces, the dispersed polaron model and the quantized classical path (QCP) approach, as well as the approximate vibronic method. These approaches do not treat the catalytic effects by phenomenological treatments and thus can be considered as first principles approaches (at least their ability to compare enzymatic reaction to the corresponding solution reactions). This work will consider the lipoxygenase reaction, and to lesser extent other enzymes, for specific demonstration. It will be pointed out that our study of the lipoxygenase reaction reproduces the very large observed isotope effect and the observed rate constant while obtaining no catalytic contribution from nuclear quantum mechanical (NQM) effects. Furthermore, it will be clarified that our studies established that the NQM effect decreases rather than increases when the donor–acceptor distance is compressed. The consequences of these findings in terms of the temperature dependence of the kinetic isotope effect and in terms of different catalytic proposals will be discussed.

This paper will also consider briefly the dynamical effects and conclude that such effects do not contribute in a significant way to enzyme catalysis. Furthermore, it will be pointed out that, in contrast to recent suggestions, NES effects are not dynamical effects and should therefore be part of the activation free energy rather than the transmission factor. In view of findings of the present work and our earlier works, it seems that TST provides a quantitative tool for studies of enzyme catalysis and that the key open questions are related to the nature of the factors that lead to transition state stabilization.

## 1. Introduction

Many proposals have been put forward to rationalize the origin of the large catalytic power of enzymes (see Warshel 1991; Villa & Warshel 2001 for a partial list). Unarguably, this is a complex issue where many factors may be assumed to be important. Thus, it is important to find out which proposals account for the major effect in enzyme catalysis. Our previous simulation studies have indicated that the major effect to enzyme catalysis comes from the pre-organization of the protein environment, where the enzyme plays a role of a super solvent with smaller reorganization energy than the corresponding reaction in aqueous solutions (Warshel 1978). However, the pre-organization idea has not yet been uniformly accepted and several alternative views have been suggested. This work will focus on some of the most prominent recent proposals.

Great interest has been placed on ideas that somehow transition state theory (TST) is not valid and that deviations from TST can account for large catalytic effect. These ideas have roughly followed two directions; in one it was implied that tunnelling and other nuclear quantum mechanical (NQM) effects provide a major catalytic factor, and in the other, it has been implied that non-equilibrium solvation (NES) effects and sometimes dynamical effects modify TST and can provide significant catalytic effects.

This work will explore the above proposals and show that they are not supported by careful simulation studies. This will be done by first describing the relevant simulation approaches and second by describing the results from different simulation studies.

## 2. Methods

### (a) Modelling reactions by empirical valence bond and other quantum mechanics/molecular mechanics approaches

The rates of chemical reactions in solutions and proteins are determined by the corresponding rate constants (e.g. Warshel 1991).(2.1)where *κ* is the transmission factor, is the average of the absolute value of the velocity along the reaction coordinate at the transition state (TS), and *β*=1/*k*_{B}*T* (where *k*_{B} is the Boltzmann constant and *T* the absolute temperature). The term Δ*G*^{‡} designates the multidimensional activation free energy that reflects the probability for the system to be in the TS region. The free energy is composed of enthalpic and entropic contributions, which also include NES effects (Villa & Warshel 2001) and, as will be shown below, NQM effects. It is also useful to comment here on the common description of the rate constant as(2.2)This Arrhenius expression is very useful for experimental analysis, but it may lead to unnecessary confusion about the factors that determine the rate constant. In other words, as is now recognized by the chemical physics community (e.g. Keck 1966; Bennett 1977; Grimmelmann *et al*. 1981), equation (2.1) provides an accurate description of the rate constant when all the dynamical effects are cast into the transmission factor and all the probabilistic effects are expressed by Δ*G*^{‡}. Indeed, Δ*G*^{‡} includes the activation entropy whereas using equation (2.2) places this important effect in *A* and makes it hard to separate the dynamical and probabilistic factors.

With the above background, we start the discussion by focusing on evaluating Δ*G*^{‡}, which in fact is the most important step. We would also like to emphasize that the ability to calculate Δ*G*^{‡} (and the corresponding free energy profiles) for enzyme reactions and the corresponding reference solution reaction is crucial for any attempt to obtain a quantitative understanding of enzyme reactions.

The common approach to obtain potential energy surfaces for chemical reactions is to use quantum chemical computational methods. Though these methods have become quite effective in treating small molecules in the gas phase (e.g. Pople 1999), we are here interested in chemical reactions in very large systems that at present cannot be explored by *ab initio* methods. Similarly, molecular mechanics (MM) simulations (e.g. Shurki & Warshel 2003), which have been proven to be very effective in exploring protein configurational space, cannot be used to describe bond-breaking and bond-making reactions in proteins or solutions. The generic solution to this problem has been provided by developing hybrid quantum mechanics/molecular mechanics (QM/MM) approaches (Warshel & Levitt 1976). These approaches divide the simulation system (e.g. the enzyme/substrate complex) into two regions. The inner region, region I, contains the reacting fragments, which are represented quantum mechanically, whereas the surrounding protein/solvent region, region II, is represented by a MM force field.

Molecular orbital (MO) QM/MM methods are now widely used in studies of complex systems in general, and enzymatic reactions in particular, and we can only mention several works (e.g. Field *et al*. 1990; Théry *et al*. 1994; Bakowies & Thiel 1996; Gao 1996; Friesner & Beachy 1998; Monard & Merz 1999; Mulholland *et al*. 2000; Zhang *et al*. 2000; Cui *et al*. 2001; Garcia-Viloca *et al*. 2001; Martí *et al*. 2001; Field 2002). Despite these advances, we are not yet at the stage where MO–QM/MM approaches can be used in fully quantitative studies of enzyme catalysis. One of the major problems is associated with the fact that evaluating the potential energy surface quantitatively for the reacting fragment requires *ab initio* electronic structure calculations, and such calculations are too expensive to allow for the configurational averaging needed for proper free energy calculations. Specialized approaches can help one move toward *ab initio* QM/MM free energy calculations (Strajbl *et al*. 2002; Olsson *et al*. 2003), but even these approaches are still in a development stage. Fortunately, one can use approaches that are calibrated on the energetics of the reference solution reaction to obtain reliable results with semi-empirical QM/MM studies, and the most effective and reliable way of doing so is the empirical valence bond (EVB) method described below.

During our search for reliable methods for studies of enzymatic reactions, it becomes apparent that in studies of chemical reactions, it is more physical to calibrate surfaces that reflect bond properties (i.e. valence-bond (VB) based surfaces) than to calibrate surfaces that reflect atomic properties (e.g. MO-based surfaces). Furthermore, it appears to be very advantageous to force the potential surfaces to reproduce the experimental results of the broken fragments at infinite separation in solution, which can effectively be accomplished with the VB picture. Though the resulting EVB method has been discussed extensively elsewhere (Warshel 1991; Aqvist & Warshel 1993), we will outline its main features below.

EVB is a QM/MM method that describes reactions by mixing resonance states (or more precisely diabatic states) that correspond to VB structures, which describe the reactant, intermediate (or intermediates) and product states. The potential energy of these diabatic states is represented by a classical MM force field of the form:(2.3)Here ** R** and

**represent the atomic coordinates and charges of the diabatic states, and**

*Q***and**

*r***are those of the surrounding protein and solvent. is the gas-phase energy of the**

*q**i*th diabatic state (where all the fragments are taken to be at infinite separation),

*U*

_{intra}(

**,**

*R***) is the intra-molecular potential of the solute system (relative to its minimum);**

*Q**U*

_{Ss}(

**,**

*R***,**

*Q***,**

*r***) represents the interaction between the solute (S) atoms and the surrounding (s) solvent and protein atoms.**

*q**U*

_{ss}(

**,**

*r***) represents the potential energy of the protein/solvent system (‘ss’ designates surrounding–surrounding). The**

*q**ε*

_{i}of equation (2.3) forms the diagonal elements of the EVB Hamiltonian (

*H*

_{ii}). The off-diagonal elements of the Hamiltonian,

*H*

_{ij}, are represented typically by simple exponential functions of the distances between the reacting atoms. The

*H*

_{ij}elements are assumed to be the same in the gas phase, in solutions and proteins. The ground state energy,

*E*

_{g}, is obtained by solving(2.4)Here,

*C*_{g}is the ground state eigenvector and

*E*

_{g}provides the EVB potential surface.

The EVB method satisfies some of the main requirements for reliable studies of enzymatic reactions. Among the obvious advantages of the EVB approach is the efficient way to obtain proper configurational sampling and converging free energy calculations. This includes the inherent ability to evaluate NES effects (Villa & Warshel 2001) and the ability to correctly capture the linear relationship between activation free energies and reaction energies observed in many important reactions (e.g. Warshel 1991). Furthermore, the EVB benefits from treating the solute–solvent coupling consistently and conveniently. This feature is essential not only for physically consistent modelling of charge-separation reactions, but also for effective use of experimental information and *ab initio* data in calibrating potential energy surfaces.

### (b) Microscopic modelling of the fluctuations of the environment and nuclear quantum mechanical effects

The EVB method provides an effective way to explore the effect of the fluctuations of the environment on proton transfer (PT) reactions. In other words, 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 emerged from our early studies (Warshel 1982, 1984; Warshel & Chu 1990; Villa & Warshel 2001) and is shown in figure 1 where we show how the fluctuations of the field from the environment change the energy of the ionic state, and thus the potential for a PT at a given fixed environment. The figure also describes the time dependence of the adiabatic ground state and the corresponding barrier (Δ*E*^{‡}) for a PT process. As indicated by the arrow, the actual transfer would occur when the product state is stabilized and Δ*E*^{‡} is reduced. As pointed out in our early studies, the fluctuations of the energy gap between the back and front panels of figure 1 tell us when Δ*E*^{‡} will be reduced (the same point was adopted in Geissler *et al*. 2001 overlooking its origin). Furthermore, these fluctuations can be used to evaluate dynamical effects by considering the autocorrelation of the time-dependent gap between the energies of the reactant and product states (see Discussion in Warshel & Parson 2001).

The fluctuations of the energy gap can also provide an interesting insight about NQM effects. In other words, we can use an approach that is formally similar to our previous treatment of electron transfer (ET) reactions in polar solvents (Warshel 1982) where we considered ET between the solute vibronic channels (e.g. Warshel 1982). Our starting point is the overall rate constant(2.5)where *β*=1/(*k*_{B}*T*) (with *k*_{B} designating the Boltzmann constant) and *E*_{am} is the energy of the *m*th vibronic level of state ‘a’. By equation (2.5), we assume that the reactant well vibrational states are populated according to Boltzmann distribution. The individual vibronic rate constant, , is evaluated by monitoring the energy difference between *E*_{am} and *E*_{bm} as a function of the fluctuations of the rest of the molecule and therefore as a function of time. The time-dependent energy gap can be used to evaluate the probability of surface crossing between the two states by adopting a semi-classical trajectory approach (Tully & Preston 1971) to rate processes in condensed phases (Warshel 1982; Warshel & Hwang 1986). This approach can best be understood and formulated by considering figure 2 and asking what is the probability that a molecule in state ‘a*m*’ will cross to state ‘b*m*′’. Treating the fluctuation of the vibronic state classically, one finds that the time-dependent coefficient for being in state , while starting from state , is given by (Warshel & Chu 1990)(2.6)where the is the Franck–Condon factor for transition from *m* to *m*′ and *H*_{ab} is the off-diagonal electronic matrix element of the EVB Hamiltonian. As shown in figure 2, the energy gap, Δ*ϵ*, is given by(2.7)The corresponding rate constant is given by (Warshel & Hwang 1986)(2.8)where *u* is the electronic energy gap relative to its average values, given by(2.9)and 〈 〉_{a} designates an average obtained over the fluctuations around the minimum of state ‘a’. The above rate constant can be treated using a cumulant expansion (Warshel & Hwang 1986) giving(2.10)In the high-temperature limit one obtains (Warshel & Hwang 1986; Warshel & Chu 1990)(2.11)where *λ* is the ‘solvent reorganization energy’ defined by(2.12)Basically, this expression reflects the probability of vibronic transition from the reactant well to the product well (as determined by the vibrational overlap integrals (the *S*_{mm′}) modulated by the chance that Δ*ϵ* will be zero. This chance is determined by the activation free energy, Δ*g*^{‡}, whose value can be approximated by the activation free energy, Δ*g*^{‡} and can be approximated by(2.13)This relationship is only applicable if the system can be described by the linear response approximation (Warshel & Hwang 1986), but this *does not* require that the system will be harmonic. The above vibronic treatment is similar to the expression developed by Kuznetsov & Ulstrup (1999). However, the treatment that leads to equation (2.13), which was developed by Warshel *et al*. (Warshel 1982; Warshel & Hwang 1986; Warshel & Chu 1990), is based on a more microscopic approach and leads to a more consistent treatment of Δ*g*^{‡} (see also below) where we can use rigorously Δ*G*^{0} rather than Δ*E*. Furthermore, our dispersed polaron (spin boson) treatment (Warshel & Hwang 1986) of equation (2.10) and if needed equation (2.13) gives a clear connection between the spectral distribution of the solvent fluctuations and the low-temperature limit of equation (2.10). It is also useful to note that Borgis & Hynes (1991) and Antoniou & Schwartz (1997) 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 is sufficiently small. Now, in cases of PT processes *H*_{ab} is far too large to allow for a diabatic treatment. However, in some cases the magnitude of is more relevant. Here we note that for 0→0 transitions, *S*_{mn} may be quite small since it is given by where Δ is the dimensionless origin shift for the PT (Δ_{H}≈6.5(Δ*r*/0.6)(*ω*/3000)^{1/2} when Δ*r* is given in Å, *ω* in cm^{−1}, Δ in dimensionless units, and this expression is based on the fact that Δ is 6.5 at Δ*r*=0.6 Å). However, *S*_{mn} approaches unity for hot transitions (vibronic transitions that involve thermally excited vibrational states when *m* is different than zero). At any rate, if the largest contribution to the rate constant comes from *S*_{00}, we may use equation (2.7) as a guide. However, the magnitude of the parameters in equation (2.7) is far from obvious. First, if we consider a collinear PT and treat all the coordinates except the X–H stretching frequency, then we have to deal with large intramolecular reorganization energy. Second, the effective frequency for the X–H stretch can be very different than the typical frequency of about 3000 cm^{−1}, once the X⋯Y distance starts to be shorter than 3.2 Å. In this range H_{12} starts to affect the ground state curvature in a drastic way (see below). Here, one can use the idea introduced by Warshel & Chu (1990) 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 *S*_{00} term it is reasonable to represent this effect by using(2.14)where *α* is approximately 0.5. It is also reasonable to use the empirical *ω*(Δ) of Mikenda (1986). Furthermore, it may also be useful to account for the complications due to the fact that the intra-molecular solute coordinate contributes to the apparent reorganization energy and is also coupled to the PT coordinate. A part of this problem can be reduced by integrating the vibronic rate constant over the ‘soft’ coordinates, and in particular the X⋯Y distance. This can be done by writing(2.15)The effectiveness of the above treatment will be examined in §4, but even if it can guide us with regard to the general trend, it cannot be used in a phenomenological way to estimate actual molecular properties. Perhaps for this reason, microscopic estimates (Olsson *et al*. 2004) of the parameters in equation (2.11) for the reaction catalysed by lipoxygenase were found to be very different than those obtained by fitting equation (2.11) to the observed isotope effects (Knapp *et al*. 2002; Lehnert & Solomon 2003; Smedarchina *et al*. 2003).

The treatment of NQM effects can be accomplished on a more quantitative level by including the adiabatic limit and modifying the centroid path integral approach (Gillan 1987; Voth *et al*. 1989; Voth 1996). The centroid path integral represents the unifying approach, which is valid both in the adiabatic and diabatic limits. This is done in a way that allows us to use classical trajectories as a convenient and effective reference for the corresponding centroid calculations. This quantized classical path (QCP) approach (Hwang & Warshel 1993, 1996) will be described briefly below.

In the QCP approach, the NQM rate constant is expressed as(2.16)where *F*_{qm}, *k*_{B}, *T* and *h* are, respectively, the transmission factor, Boltzmann constant, temperature and Planck constant, and *β*=1/*k*_{B}*T*. The quantum mechanical activation barrier, , includes almost all the NQM effects, whereas only small effects come from the pre-exponential transmission factor in the case of systems with a significant activation barrier (Billeter *et al*. 2001; Warshel & Parson 2001).

The quantum mechanical free energy barrier, , can be evaluated by Feynman's path integral formulation (Feynman 1972), where each classical coordinate is replaced by a ring of quasiparticles that are subjected to the effective ‘quantum mechanical’ potential(2.17)Here, (where ), , *M* is the mass of the particle 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 *U*_{qm}. The probability of being at the TS is this way approximated by a probability distribution of the centre of mass of the quasiparticles (the centroid) rather than the classical single point.

Such calculations of centroid probabilities in the condensed-phase reactions are very challenging, since they may involve major convergence problems. The QCP approach offers an effective and rather simple way for evaluating this probability without significantly changing the simulation program. This is done by propagating classical trajectories on the classical potential 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 as (Hwang *et al*. 1991; Hwang & Warshel 1993)(2.18)where is the centroid position, designates an average over the free particle quantum mechanical distribution obtained with the implicit constraint that coincides with the current position of the corresponding classical particle, and designates an average over the classical potential *U*. It is worth stressing that path integral calculations involving computationally expensive quantum chemical evaluation of forces and energies would benefit much from the QCP scheme. In quantum chemical calculations involving quasiparticles, one cannot realize exclusions between quasiparticles and therefore computational effort is proportional to the number of the quasiparticles in the necklace. In the present approach, the quantum connection can be performed *a posteriori*, using the stored trajectory. It would be interesting to apply this approach to a Car-Parrinello path integral scheme and demonstrate almost negligible increases of the CPU time relative to the classical treatment of the nuclei (Tuckerman & Marx 2001).

Using equation (2.17), we can obtain the quantum mechanical free energy surface by evaluating the corresponding probability by the same combined free energy perturbation umbrella sampling (FEP-US) approach that has been repeatedly applied in our classical simulations as well as in our quantum mechanical simulations, but now we use the double average of equation (2.18) rather than an average over a regular classical potential. The actual equations used in our FEP-US calculations are given elsewhere (e.g. Hwang & Warshel 1993, 1996), but the main point of the QCP is that the quantum mechanical free energy function can be evaluated by a centroid approach, which is constrained to move on the classical potential.

## 3. Simulation studies of dynamical effects, tunnelling and related nuclear quantum mechanical effects

Recent studies (see Kohen & Klinman 1999 and references in Villa & Warshel 2001) have suggested that vibrationally enhanced tunnelling (VET) of PT, hydride transfer and hydrogen transfer reactions play a major role in enzyme catalysis. According to this interesting proposal, nature created by evolution protein vibrational modes that are strongly coupled to the hydrogen atom motion. Some workers (e.g. Sutcliffe & Scrutton 2000) assumed that there is an entirely new phenomenon that makes TST inapplicable to enzymatic reactions. However, the VET effect is not new and is common to many chemical reactions in solution (German 1981; Warshel 1982; Borgis & Hynes 1991). Moreover, the VET is strongly related to TST. In other words, when the solvent fluctuates and changes the energy gap (Warshel 1982, 1984) the light atom sees a fluctuating barrier that leads in some cases to a large contribution from nuclear tunnelling. As shown in Warshel (1984), these fluctuations are taken into account in the statistical factor of the classical TST and the same is true when quantum effects are taken into account. Thus, the recent realization that the solvent coordinates should be considered in tunnelling studies is neither new nor means that this effect is important in catalysis.

Warshel & Chu (1990) and Hwang *et al*. (1991) were the first to actually calculate the contribution of tunnelling and other nuclear quantum effects to PT in solution and enzyme catalysis, respectively. Since then, and in particular in the past few years, there has been a significant increase in simulations of NQM effects in enzyme and in solution reactions (Friesner & Beachy 1998). The approaches used range from the QCP (e.g. Hwang & Warshel 1996; Feierberg *et al*. 2000; Villa & Warshel 2001), the centroid path integral approach (Gillan 1987; Voth 1996), and variational TST (Alhambra *et al*. 2000), to the molecular dynamics with quantum transition (MDQT) surface hopping method (Billeter *et al*. 2001) and density-matrix evolution (Berendsen & Mavri 1993; Lensink *et al*. 1999; Mavri & Grdadolnik 2001). Most studies of enzymatic reactions did not examine the reference water reaction, and thus could only evaluate the quantum mechanical contribution to the enzyme rate constant, rather than the corresponding catalytic effect. However, studies that explored the actual catalytic contributions (e.g. Hwang & Warshel 1996; Feierberg *et al*. 2000; Villa & Warshel 2001) concluded that the quantum mechanical contributions are similar for the reaction in the enzyme and in solution, and thus do not contribute to catalysis.

Interestingly, the MDQT approach of Hammes-Schiffer *et al*. (Billeter *et al*. 2001) allowed them to explore the quantum mechanical transmission factor. It was found that even with quantum mechanical considerations, the transmission factor is not so different than unity, and thus, we do not have a large dynamical correction to the TST rate constant.

It is important to clarify here that the description of PT processes by curve crossing formulations is not a new approach nor does it provide new dynamical insight. In other words, the view of PT in solutions and proteins as a curve crossing process has been formulated in early realistic simulation studies (Warshel 1982, 1984; Warshel & Chu 1990) with and without quantum corrections, and the phenomenological formulation of such models has been introduced even earlier by Kuznetsov and others (Kuznetsov & Ulstrup 1999). Furthermore, the fact that the fluctuations of the environment in enzymes and solution modulate the activation barriers of PT reactions has been demonstrated in realistic microscopic simulations of Warshel (1982, 1984). However, as clarified in these works, the time dependence of these fluctuations does not provide a useful way to determine the rate constant. In other words, the electrostatic fluctuations of the environment are determined by the corresponding Boltzmann probability and do not represent a dynamical effect. In other words, the rate constant is determined by the inverse of the time it takes the system to produce a reactive trajectory, multiplied by the time it takes such trajectories to move to the TS. The time needed for generation of a reactive trajectory is determined by the corresponding Boltzmann probability and the actual time it takes the reactive trajectory to reach the TS (on the order of picoseconds), and is more or less constant in different systems.

It is also important to clarify that the solvent reorganization energy, which determines the amplitude of the solvent fluctuations, is not a ‘static dynamical effect’ (as proposed in Knapp & Klinman 2002) but a unique measure of the free energy associated with the reorganization of the solvent from its reactant to its product configuration (see King & Warshel 1990 for a more rigorous definition). In fact, the reorganization energy, *λ*, and the reaction energy, Δ*G*, determine the activation free energy (Δ*g*^{‡}) and the corresponding Boltzmann probability of reaching the TS (see Discussion in Warshel 1984; Villa & Warshel 2001). Now, since the Δ*g*^{‡} is a probability factor, it can be determined by Monte Carlo simulations without any dynamical considerations. Furthermore, since the transmission factor (*κ* in equation (2.1)) is the only dynamical part of the rate constant and since *κ* is close to unity in enzymes and solutions (e.g. Villa & Warshel 2001), the corresponding rate constants do not show significant dynamical effects.

## 4. Lipoxygenase as a test case

At this point, we find it useful to discuss the specific case of lipoxygenase, which has been brought recently (Ball 2004) as an example of the role of NQM effects in catalysis. The catalytic reaction of lipoxygenase involves a very large isotope effect (approx. 80) and thus it is tempting to suggest that the enzyme catalyses this reaction by enhancing the tunnelling effect. This proposal (Knapp *et al*. 2002) has been examined by several theoretical approaches ranging from phenomenological models (Kosloff & Kosloff 1983; Billeter *et al*. 2001; Knapp *et al*. 2002) to continuum representations of the protein (Hatcher *et al*. 2004) and to complete microscopic treatment of the system using the QCP approach (Olsson *et al*. 2004). Here, we will consider both the vibronic and the QCP studies and clarify their limitations. We will also emphasize the conclusions that emerged from these treatments.

### (a) Vibronic treatment

Several research groups (e.g. Knapp *et al*. 2002; Lehnert & Solomon 2003), have used the vibronic treatment of equation (2.11) as a phenomenological tool for fitting the calculated and observed kinetic isotope effect (KIE) of lipoxygenase. The most extensive study has been done by Klinmann *et al*. (Knapp *et al*. 2002) who obtained a reorganization energy of 20 kcal mol^{−1} in their fitting. The problems with the magnitude of this parameter will be discussed below. Hammes-Schiffer *et al*. (Hatcher *et al*. 2004) have progressed beyond the phenomenological treatment using a more realistic microscopic treatment of the solute (on the other hand, the protein effect was modelled macroscopically with an arbitrarily low dielectric constant, making it hard to explore the actual role of the protein). These workers treated all hydrogen motions quantum mechanically (rather than just the stretching mode as was done in Knapp *et al*. 2002). They estimated the intramolecular reorganization energy to be around 20 kcal mol^{−1} from microscopic considerations, but in contrast to the possible implication of Hatcher *et al*. (2004), this reasonable estimate should not be similar to the estimate obtained in Knapp *et al*. (2002) (see below). Hatcher *et al*. (2004) estimated the outer sphere reorganization energy to be around 2 kcal mol^{−1}, representing the protein as a sphere with *ϵ*=4. Now, although this estimate is similar to the *λ*_{out}∼2.5 kcal mol^{−1} obtained by the careful microscopic simulations of Olsson *et al*. (2004), it is somewhat arbitrary since the protein cannot be represented as a continuum and the effective dielectric of such a model is system dependant (Muegge *et al*. 1997).

The inter- and intramolecular reorganization energies of lipoxygenase were examined by dispersed polaron simulation by Olsson *et al*. (Warshel *et al*. 2005). The results of the simulations, which are summarized in figure 3, indicate that we are dealing with a highly coupled system that might require special convolution treatments to extract the effective multidimensional Franck–Condon factors the *S*_{m,m′} of equation (2.6). Although such a treatment is possible, we simply evaluated the inter- and intramolecular reorganization energies using the long-time average of equation (2.12). The corresponding reorganization energies (Olsson *et al*. 2004) were *λ*_{in}∼100 and *λ*_{out}∼2.5 kcal mol^{−1} (please note that the value of *λ*_{in} is completely correlated with H_{12} where larger *λ*_{in} can be accommodated by H_{12}). Before we move further to the actual vibronic treatment, it is useful to consider the reorganization energies obtained by different treatments.

As stated above, Hammes-Schiffer *et al*. (Hatcher *et al*. 2004) obtained reasonable intramolecular reorganization energy (*λ*_{in}∼20 kcal mol^{−1}) by a treatment that described the hydrogen bending motions quantum mechanically and thus did not include their effect in *λ*_{in}. Knapp *et al*. (2002), on the other hand, treated only the hydrogen stretching quantum mechanically, and also obtained *λ*_{in}∼20 kcal mol^{−1} by a phenomenological fitting procedure (which also presumably included *λ*_{out} in the single reorganization parameter). In view of the fundamental difference between both models and the results of our simulation, we believe that the effective *λ* in the model of Knapp *et al*. should be significantly larger than 20 kcal mol^{−1}. We also find it useful to clarify that the *λ* in Knapp *et al*. (2002) should not be confused with the *λ* which has been found by our microscopic treatment to be very small (*λ*_{out}∼2.5 kcal mol^{−1}).

While the above discussion of the reorganization energy can be instructive, we are not interested in a fully vibronic treatment considering the problems with the overall validity of such a treatment. Thus, we limit our analysis to only one quantum mechanical mode (the hydrogen stretch) and treat the rest of the solute and solvent classically. Furthermore, we will focus on the effect of the donor–acceptor distance, which we examine in a parametrical way. The result of our treatment is summarized in figures 4–6. Figure 4 shows different lipoxygenase reaction profiles in water and in the protein for various donor–acceptor distances. As can be seen from figure 4, and which was also briefly discussed in §2*b*, the effective ground state curvature and the corresponding frequency (*ω*) decrease with a shorter donor–acceptor distance. Figures 5 and 6 describe the dependence of *k*_{H} and *k*_{H}/*k*_{D} on the displacement of the two parabolas and on the effective stretching frequency (*ω*). Again, one should keep in mind that the effective *ω* decreases when Δ*R* is decreased due to the effect of H_{12}. As already pointed out in §2*b*, the isotope effect increases upon increasing Δ*R*, reaching a value of 80 at Δ*R*=0.6 Å for *ω*=1500 cm^{−1} (which seems to be a reasonable representation of the lipoxygenase case). Another way to look at this problem is to assume a reasonable potential of mean force (PMF) (*U*(*R*)) in equation (2.15) and *ω*(*R*) according to equation (2.14) or according to the empirical relationship of Mikenda (1986) and then to evaluate an average *k*_{ab}. Such a treatment is summarized in figure 7, where we used the empirical dependence of Mikenda (1986). Indeed, it is rather trivial to reproduce the temperature dependence of the KIE by adjusting *U*(*R*). For example, the observed trend can be reproduced (see figure 8) using *U*(*R*)=0.6((3.4/*R*)^{12}−(3.4/*R*)^{6}). Although the agreement is not perfect, we do not find it useful to improve the agreement by refining *U*(*R*) or the empirical *ω*(*R*). In other words, the vibronic treatment can only be considered as a qualitative guide and figure 8 is used only to demonstrate that the temperature dependence can be reproduced easily by phenomenological parameter fitting approaches.

The take-home message from the above treatment is that the isotope effect and tunnelling correction increase rather than decrease when Δ*R* increases. The drastic dependence of the tunnelling correction on Δ*R* is evidenced from the fact that the NQM correction for the H-transfer reaction reaches a factor of 10^{6}. Indeed, this must be due to tunnelling corrections rather than to zero-point energy effects. The implication from this finding will be discussed below.

### (b) Quantum classical path treatment of the lipoxygenase reaction

Since the vibronic treatment cannot be considered as a quantitative analysis, we also consider here the QCP treatment, which considers the NQM effect and the dynamics of the entire protein in a consistent way. In view of the findings from the above studies, we explored the dependence of the rate constant on Δ*R*. We started by evaluating the EVB surface while fixing the donor–acceptor distance to different Δ*R* values (figure 4) for both the enzyme and solution reaction. Next, we performed QCP simulations for each of these constrained systems. The corresponding changes in activation free energy (ΔΔ*g*^{‡} and its components at the reactant state (RS) and TS regions) and the KIE are shown in figure 9. As seen in figure 9, the KIE increases from a comparatively small value to several hundreds when Δ*R* is increased from 2.5 to 2.9 Å. This primarily reflects an increased curvature and NQM tunnelling correction in the TS region. Qualitatively, the above vibronic treatment shows the same behaviour as this QCP approach, but this approach is a quantitative fully microscopic treatment that takes into account the entire protein with its static and dynamic effects. We are not aware of any related study by other research groups.

### (c) General conclusions about activation barriers and catalysis in lipoxygenase

The above discussion has demonstrated the qualitative insight obtained from the vibronic treatment and the quantitative result obtained by the QCP treatment. With this background, we can address several remaining important issues that are directly related to the overall catalytic effect.

The first point that should be recognized is the finding of the QCP simulations that the NQM effects are very similar in the enzyme and solution reactions. This means that NQM effects do not lead to a significant catalytic effect in lipoxygenase. Here, it is important to point out that the rate constant is determined almost exclusively by the activation barrier of equation (2.16) and that the actual transmission factor is close to unity. A common confusion in this respect emerges from the problematic focus on Δ*E*^{‡} of equation (2.2). Basically the entropic and enthalpic contributions to Δ*g*^{‡} are relatively unstable and can easily change with mutations and temperature. What is stable and thus useful for analysis and predictive purposes is the free energy. In other words, it is far more constructive to exploit the entropy/enthalpy compensation and to focus on the activation free energy rather than to attempt to simulate or analyse its contributions. The discussion and to some extent the confusion about the activation free energy can be traced not only to the use of equation (2.2), but also to the Δ*g*^{‡} of the vibronic treatment (equation (2.11)). In other words, it has been argued that the reason for the temperature dependence of KIE in lipoxygenase is the small activation energy that was deduced from the (Δ*G*+*λ*)^{2}/4*λ* term in equation (2.13) (without the vibrational sum). Unfortunately, this term does not represent the actual Δ*g*^{‡}, but only the contribution to the classical modes to the transition between the vibronic levels. The correct Δ*g*^{‡} includes the barrier associated with the quantum modes (i.e. the hydrogen stretching mode). Treating these modes quantum mechanically does not eliminate the corresponding contribution to the activation barrier. In fact, the activation free energy of the reaction in lipoxygenase is much higher than the corresponding Δ*E*^{‡}, which is about 14 kcal mol^{−1} (table 1).

In concluding this section, it is useful to point out an interesting conclusion that emerged from the use of equations (2.14) and (2.16). Using these equations, we found that the degree of tunnelling (or at least the magnitude of the isotope effect) decreases, rather than increases, when the donor–acceptor distance is reduced. This reflects the reduction in the effective diabatic X–H stretching frequency. This means that the common idea that enzymes can catalyse H-transfer or PT reactions by compressing the donor and acceptor complexes (Ball 2004) is very problematic. Such compression will in fact reduce the tunnelling contributions, since the largest contribution to the KIE comes from , whose value decreases drastically when Δ_{H} and *ω*_{H} decrease. Furthermore, as we have shown repeatedly, enzymes are flexible and unlikely to be able to change drastically the reaction surface (Warshel 1991; Shurki *et al*. 2002). Thus, it is quite likely that the main difference between reactions with large tunnelling and small tunnelling corrections is the intrinsic shape of the potential surface rather than the effect of the environment on this surface. It is thus possible that radical reactions involve steep potential surfaces with relatively small H_{12}, while other reactions involve large H_{12} and shallow surfaces with small effective *ω*_{H}.

### (d) The temperature dependence of the kinetic isotope effect in lipoxygenase and other enzymes

The temperature dependence of the KIE in lipoxygenase has been a subject of significant interest (e.g. Knapp *et al*. 2002), in part because of the assumption that this dependence provides a support to the idea that dynamical effects contribute to catalysis. Our attempt to reproduce the temperature dependence of the rates in lipoxygenase reproduced the observed trend in the H-transfer reaction, but was less successful in the D-transfer reaction (Olsson *et al*. 2004), although more extensive simulations (M. H. M. Olsson & A. Warshel 2005, unpublished work) improved the results for the D-transfer reaction. Unfortunately, the difficulties in simulating the temperature dependence led some workers to assume that there are fundamental problems with our simulations or with our conclusions about the small catalytic contribution of NQM effects. Thus, it is important to clarify several misunderstandings.

First, as shown in §4*a*, since the isotope effect depends very strongly on the donor–acceptor distance, it is trivial to reproduce the temperature dependence by empirically changing the form of *U*(*R*) (figure 8). However, this is very different from actual unbiased simulations that include explicitly the protein in the simulation, where it is quite hard to reproduce the exact shape of the very shallow *U*(*R*). Presently, the only attempt to reproduce the temperature dependence by simulations that include the protein is due to Olsson *et al*. (2004) (note in this respect that the instructive work of Hatcher *et al*. 2004 represented the protein by a phenomenological continuum model and fit rather than obtained the temperature dependence). Thus, the difficulties in reproducing the KIE by microscopic simulations only reflect difficulties in reproducing the exact shallow *U*(*R*) which is not so relevant to catalysis (see below).

Second, our calculations (table 1) showed that the activation free energy is reproduced almost quantitatively even by simulations that do not reproduce the temperature dependence of the KIE and these calculations give a very similar NQM corrections to the reaction in solution and in the enzyme active site. Furthermore, the fact that PMF for the donor–acceptor distance is shallow means that the enzyme does not compress the donor and acceptor to a close distance and does not contribute to catalysis by such a compression effect. In fact, the large tunnelling correction indicates that the enzyme might even pull the donor and acceptor away from each other.

## 5. General comments about dynamical effects

The idea that dynamical effects play a crucial role in enzyme catalysis requires some additional discussion (in addition to the above discussion of quantum mechanical effects) in view of the perception that TST may not be applicable to such studies. Thus, we will consider here several proposals.

### (a) Non-equilibrium solvation is a thermodynamic rather than dynamical effect

Some studies of reactions in solution have suggested that ‘NES’ can affect the transmission factor for a reaction (Bergsma *et al*. 1987; Gertner *et al*. 1989; Truhlar *et al*. 1993; Muller & Warshel 1995; Neria & Karplus 1997; Villa & Warshel 2001). In fact, Garcia-Viloca *et al*. (2004) recently implied that these effects comprise an important correction to TST and provide a key modern consideration for understanding enzymatic reactions. However, NES is not really a dynamical effect or part of the transmission factor, but rather a well-defined free energy factor (a part of Δ*g*^{‡}) that has been recognized for some time (Muller & Warshel 1995; Villa & Warshel 2001; Warshel & Parson 2001) and has been included consistently in almost all EVB studies. When the solute is at its TS, the solvent reorganization energy creates a barrier between the equilibrium configurations of the solvent in the reactant and product states. This barrier constitutes the entire activation free energy in outer sphere ET reactions, and its contribution to Δ*g*^{‡} can be significant in other types of reactions as well (Muller & Warshel 1995).

The contribution of NES to Δ*g*^{‡} is overlooked in PMF treatments of reaction kinetics in which the solvent and solute coordinates are handled separately and the solvent is allowed to relax on its reaction coordinate during each step along the solute coordinate (Warshel & Parson 2001).

Unfortunately, FEP mapping procedures do not reflect the reorganization energy required for changing the solvent coordinate between the mapping steps, and the TS for the solvent coordinate is never actually sampled. The resulting underestimate of Δ*g*^{‡} will be most serious if the coupling of the solute and solvent charges also is neglected.

The contribution of NES to Δ*g*^{‡} can be studied by constraining the EVB solute to its TS geometry (*R*^{‡}) and using an FEP procedure to evaluate the free energy surface for moving along the solvent coordinate from its equilibrium configuration on one side of the TS to its equilibrium configuration on the other side (Warshel & Parson 2001). For example, figure 10 shows calculations of the NES effect for a step in the reaction catalysed by subtilisin, along with similar calculations for the corresponding process in solution (Villa & Warshel 2001). NES contributes about 4 kcal mol^{−1} to the barrier in solution, and about 1 kcal mol^{−1} in the enzyme. The decrease in can be viewed as one feature of the pre-organized environment in the enzyme's active site. Using direct simulations of downhill trajectories gave a transmission factor of approximately 0.6 for both the enzyme and the solution reaction, showing clearly that the NES effect is not a dynamical effect (Warshel & Bentzein 1999).

The EVB/FEP-US method described in §2 incorporates the NES effect automatically by calculating the probability of finding the system at the TS of the solute–solvent system, without having to divide Δ*g*^{‡} explicitly into equilibrium and non-equilibrium components (Hwang *et al*. 1988). Therefore, the EVB Δ*g*^{‡} should differ from by approximately (Bentzien *et al*. 1998; Warshel & Bentzien 1999).

Neria & Karplus (1997) have discussed effects of NES on the PT step in triosephosphate isomerase. They used EVB simulations of a ‘frozen-solvent’ model described by Hynes *et al*. (Gertner *et al*. 1987, 1989) in which fluctuations of the solvent are assumed to be much slower than those of the solute. In this model, the rate constant is written(5.1)where *k*_{PMF} is the TST rate constant obtained by using the PMF as the activation energy of the system and *κ*_{froz} is a transmission factor that represents the NES obtained within the frozen solvent approximation. A careful consistent analysis (Villa & Warshel 2001; Warshel & Parson 2001) has demonstrated that the above treatment is simply unjustified. Because of space limitation, we will not repeat the consideration of Aqvist & Warshel (1993) and Warshel & Parson (2001), but suggest that the reader will examine the relevant discussion.

In conclusion, critical studies of NES effects show that these effects are similar in enzymes and in solution. The main differences appear to be associated with magnitude of the reorganization energy rather than with dynamical effects.

### (b) The nature of reactive vibrational modes is similar in enzymes and solution

The EVB approach allows one to evaluate the projections of protein or solvent motions along the reaction coordinate. In the ‘dispersed-polaron’ or ‘spin-boson’ treatment, this is done by relating the fluctuations of Δ*ϵ*_{12} during an MD trajectory to the fluctuations of an equivalent harmonic system. Using this approach, we were able to show that the distribution of projections of the protein modes on the reaction coordinates of enzymatic reactions are similar to the projections of the solvent modes on the reaction coordinates of the corresponding reference reactions (e.g. Olsson & Warshel 2004). This finding indicates that the nature of the vibrational modes that lead to the TS does not play a significant role in catalysis. This point of view was further strengthened by running downhill dynamics and examining the nature of the motions of the reactive trajectories (e.g. Olsson & Warshel 2004). Such studies have demonstrated that the reactive trajectories are similar in proteins and solution and the main difference is associated with the amplitude of these motions, which are a direct reflection of the corresponding reorganization energy.

The idea that increasing the rate of vibrational equilibration could lead to dynamical effects (Cui & Karplus 2002) merits some discussion. We are not aware of any formal expression that relates the quantum mechanical transmission factor to the classical time of crossing the TS or the rate of energy redistribution. In a classical picture, the speed of a single crossing of the TS is constant; what counts is the time required to dissipate an amount of energy in the order of *k*_{B}*T* and thus to make the barrier crossing effectively irreversible. Nevertheless, we have here a picture of fast solvent relaxation in both the enzyme and solution, which means that changes in the rate of this relaxation are unlikely to contribute significantly to catalysis.

## 6. Concerted motions in an enzyme usually do not make dynamical contributions to catalysis

Nuclear magnetic resonance spectroscopy has been used to study the motions of enzymes over a wide range of time-scales (Palmer *et al*. 2001; Wang *et al*. 2001; Akke 2002; Eisenmesser *et al*. 2002; Palmer 2004). An interesting case is the enzyme cyclophilin A, which catalyses *cis–trans* isomerization of peptidyl proline bonds (Akke 2002; Eisenmesser *et al*. 2002; Falke 2002; Li & Cui 2003). A detailed analysis of the relaxation dynamics indicated that numerous residues of the protein have motions in the millisecond time range, coinciding with the turnover time of the enzyme, and that some of these motions change when substrate is added. In particular, the transverse relaxation of Arg 55, a residue that is hydrogen bonded to the substrate and is essential for catalysis, accelerates in the presence of the substrate and approximately matches the dynamics of forming the TS, suggesting that motions of Arg 55 might play a dynamical role in the catalytic mechanism (Eisenmesser *et al*. 2002). In analysing this proposal, it is important to bear in mind that if Arg 55 or any other residue moves along the reaction coordinate, and if its position changes in the TS, it necessarily moves on the same time-scale as the reaction. For example, the first shell of solvent molecules in a reaction in solution must rearrange during the reaction and so must move at more or less the same rate as the solute atoms. Most of the reorganization of the environment occurs on the same time-scale as the reaction, and the motions of protein residues near the reacting substrate are not fundamentally different in this regard from the motions of the solvent molecules in solution (e.g. Olsson & Warshel 2004). Further, as long as the motions of the protein residues follow Boltzmann's law, they simply reflect probabilistic effects and not *bona fide* dynamical effects. In a QM/MM study of cyclophilin A, Li & Cui (2003) found that Arg 55 underwent only very small displacements between the enzyme–substrate (ES) complex and the TS. They concluded that the Arg residue stabilizes the TS electrostatically, and that its motions do not contribute significantly to the catalysis.

Nunez *et al*. (2004) recently suggested that the catalytic reaction of purine nucleoside phoshorylase involves protein modes that reduce the barrier height by 20% by compressing the reacting fragments. To evaluate the contribution of such modes to catalysis, one must calculate the barrier height from the minimum in the ground state, taking into account the energy associated with the compression. Without such an analysis, it might appear that a sufficiently strong compression would eliminate the barrier completely. In addition, it is important to bear in mind that compression modes similar to those that occur in the protein also occur in the reference solution reaction. In the cases that we have studied, the costs of bringing the reactants to the same distance are similar in solution and in the enzyme, which means that the compression does not contribute significantly to catalysis.

Benkovic *et al*. (Epstein *et al*. 1995; Cameron & Benkovic 1997; Miller & Benkovic 1998; Schnell *et al*. 2004*a*,*b*) have studied the reaction of dihydrofolate reductase by nuclear magnetic resonance (NMR). They found that site-directed mutations of residues in a loop that undergoes relatively large backbone motions had detrimental effects on catalysis, and they suggested that the dynamics of these residues could be important for catalysis. This suggestion was supported by Brooks *et al*. (Radkiewicz & Brooks 2000; Rod *et al*. 2003) who carried out MD simulations of three ternary complexes of the enzyme. Motions of some residues were strongly correlated, and were different in the ES and enzyme–product (EP) complexes. Some of these motions were modified in simulations of mutant enzymes with diminished activity. However, these studies did not examine any of the TSs in the reaction or demonstrate any dynamical effects on the rate constant. The different motions of the ES and EP complexes could just reflect the coupling of ES interactions to interactions of various groups in the protein which are common to all enzymes. In simulations using the MDQT approach, Hammes-Schiffer *et al*. (e.g. Watney *et al*. 2003; Hammes-Schiffer 2004) identified a network of correlated conformational changes with projections on the reaction path, but concluded that these reflect equilibrium structural effects rather than dynamical effects.

It is important to emphasize that identification of correlated motions does not provide a new view of enzyme catalysis, because reorganization of the solvent along the reaction path in solution also involves highly correlated motions (Warshel 1984; Hwang *et al*. 1991). Correlated motions of an enzyme do not necessarily contribute to catalysis, and indeed could be detrimental if they increase the reorganization energy of the reaction. The EVB and dispersed-polaron approaches described above consider the enzyme reorganization explicitly and automatically assess the complete structural changes along the reaction coordinates. A dispersed-polaron analysis of the type presented in figure 3 determines the projection of the protein motion on the reaction coordinate and provides a basis for a quantitative comparison with a reference reaction in solution.

The coupling of protein motions to a reaction in an enzyme involves fluctuating electrostatic interactions of the solute with charged or polar residues and bound water molecules. In solution, it involves reorientation of the dipoles in the solvation shells. Clearly, the reaction coordinate in both cases will involve components along the environment (solvent) coordinate. The real difference is the amplitude of the change in the solvent coordinates during the reaction, which determines the reorganization energy and generally is smaller in the enzyme because of pre-organization of the active site.

## 7. Concluding remarks

Our studies have demonstrated that NQM effects do not contribute significantly since similar effects occur in the reference solution reactions. We have also pointed out that dynamical effects do not play an important role in enzyme catalysis. Nevertheless, one should be able to calculate NQM effects and treat dynamical effects in order to reproduce the actual rate constant in enzymes and to analyse and discriminate between different catalytic proposals. The ability to address such problems has been provided by our approaches.

## Acknowledgements

This work was supported by NIH grants GM-24492 and GM-40283 and NSF grant MCB-0342276. One of us (J.M.) thanks the J. William Fulbright Scholarship Board for the award of a research scholarship. We also gratefully acknowledge the University of Southern California's High Performance Computing and Communications Center for computer.

## Footnotes

One contribution of 16 to a Discussion Meeting Issue ‘Quantum catalysis in enzymes—beyond the transition state theory paradigm’.

- © 2006 The Royal Society