Scholarly article on topic 'Studying the role of protein dynamics in an SN2 enzyme reaction using free-energy surfaces and solvent coordinates'

Studying the role of protein dynamics in an SN2 enzyme reaction using free-energy surfaces and solvent coordinates Academic research paper on "Chemical sciences"

0
0
Share paper
Academic journal
Nature Chemistry
OECD Field of science
Keywords
{""}

Academic research paper on topic "Studying the role of protein dynamics in an SN2 enzyme reaction using free-energy surfaces and solvent coordinates"

ARTICLES

PUBLISHED ONLINE: 26 MAY 2013 | DOI: 10.1038/NCHEM.1660

nature

chemistry

Studying the role of protein dynamics in an SN2 enzyme reaction using free-energy surfaces and solvent coordinates

Rafael Garcia-Meseguer1, Sergio Marti2, J. Javier Ruiz-Pernia2*, Vicent Moliner2 and Inaki Tunon1*

Conformational changes are known to be able to drive an enzyme through its catalytic cycle, allowing, for example, substrate binding or product release. However, the influence of protein motions on the chemical step is a controversial issue. One proposal is that the simple equilibrium fluctuations incorporated into transition-state theory are insufficient to account for the catalytic effect of enzymes and that protein motions should be treated dynamically. Here, we propose the use of free-energy surfaces, obtained as a function of both a chemical coordinate and an environmental coordinate, as an efficient way to elucidate the role of protein structure and motions during the reaction. We show that the structure of the protein provides an adequate environment for the progress of the reaction, although a certain degree of flexibility is needed to attain the full catalytic effect. However, these motions do not introduce significant dynamical corrections to the rate constant and can be described as equilibrium fluctuations.

One of the most intriguing characteristics of enzymes is their flexibility. It has been stressed that, to function, enzymes must be stable enough to retain their three-dimensional (3D) structure, but flexible enough to permit the evolution of the protein among the different conformations relevant at each step of the full catalytic process1,2. For example, in many cases conformational changes are known to be a requisite for substrate binding and product release. In fact, the multidimensional free-energy surface that corresponds to an enzymatic process is very rugged, and contains multiple minima that appear along the multiple conformational coordinates available for the protein3. Transitions between different conformational substrates of the macromolecule and their population distribution are governed by motions that happen in a broad range of timescales, from milliseconds to femtoseconds2'4.

However, the impact of protein flexibility on the rate constant of the chemical step remains the subject of a long-standing debate in scientific literature5-15. Even if the active site is designed to catalyse the reaction and then to accommodate the charge distribution of the transition state (TS)16, some protein motions are still needed to evolve from the reactant state (RS) to the TS^^17. In other words, a certain degree of protein reorganization is needed to attain the TS15. This reorganization would play a similar role to that of solvent polarization in the Marcus theory for electron transfer in solution18, although much slower conformational components could be involved in enzymatic catalysis12,19. In general, the coordinate that describes the transformation of the system from reactants to products is a collective coordinate that involves the degrees of freedom not only of the solute/substrate, but also of the environment (the solvent and/or the enzyme).

Theoretical descriptions of chemical reactions in enzymatic environments are usually based on the selection of a distinguishing reaction coordinate, typically defined in terms of the substrate or solute coordinates (that is, some valence coordinates related to the bonds to be broken and/or formed). Within this description, some protein motions can be coupled to the progress of the

system along the reaction coordinate. For example, correlated motions within the protein promote tunnelling in hydride-transfer reactions20,21 and compressive local motions within the active site facilitate the approach between the donor and the acceptor atoms in transfer reactions11,22. The participation of protein motions in the progress of the reaction is thus a well-established fact invoked using different terminology: protein reorganization15, coupled motions9 or promoting vibrations8,22. A more controversial question is the way in which these motions must be described and whether the current theoretical frameworks are adequate or not to explain the rate of enzymatic reactions14,22.

As a fundamental approach to describe the reaction rate of chemical reactions, transition state theory (TST) provides the tools for analysing the rate of enzymatic reactions23. The basic assumption in TST is that the selected reaction coordinate (Z) is separable from the rest of the coordinates of the system such that the averaged dynamic behaviour of the trajectories that evolve from the reactants to products can be represented as the equilibrium flux across the dividing surface. In this case the rate constant can be related simply to the free-energy difference between the TS and the reactants (AG*). For a unimolecular process

k-kT k - h

where kB, h and R are the Boltzmann, Planck and gas constants, and T is the absolute temperature. This approach means that the remaining degrees of freedom can be considered at equilibrium at any value of the reaction coordinate, described as a Boltzmann distribution. If the reaction coordinate is defined exclusively in terms of the degrees of freedom of the substrate (or 'solute coordinate'), then the protein degrees of freedom (or 'solvent coordinate') are assumed to be at equilibrium at all the stages of the reaction path. Obviously, protein dynamics spans a hierarchy of timescales2,4,10 and only those motions much faster than the reaction coordinate can be

'Departament de Química Física, Universität de Valencia, 46100 Burjassot, Spain, 2Departament de Química Física i Analítica; Universität Jaume I, 12071 Castellón, Spain. *e-mail: ignacio.tunon@uv.es; pernia@qfa.uji.es

considered to be at equilibrium. Thus, this approach breaks down for motions coupled to the solute coordinate and that take place in similar or slower timescales. In such a case an explicit treatment of vibrational dynamics in the enzyme may be required. This is suggested for some cases covering both fast vibrational motions22 and slow conformational changes19. However, other analyses stress that the role of protein motions can be incorporated satisfactorily in the description of the chemical step as equilibrium fluctuations and thus explicit dynamic treatments of protein motions are not really necessary for modelling enzymatic catalysis5,9,13,14,24.

TST offers a convenient framework to incorporate non-statistical effects of protein motions. A transmission coefficient (k) can be calculated from time-dependent simulations and incorporated into the TST expression of the rate constant25:

kBT _DG*(0

As the transmission coefficient takes values lower than unity the equilibrium approach (equation (1)) provides an upper limit to the rate constant. There is no unambiguous way to separate the effects of the enzyme on the activation free energy (or on the equilibrium) from those on the transmission coefficient (non-equilibrium effects) because both depend on the choice of the reaction coordinate23. If this is defined exclusively on the basis of the coordinates of the substrate, deviations from the equilibrium distribution of protein motions should be reflected in the transmission coefficient. Assuming that these effects are important in the vicinity of the dynamic bottleneck, the transmission coefficient can be evaluated from the frictional force exerted on the reaction coordinate at the TS26 or by means of rare-event simulations that count the number of recrossings across the dividing surface27. To date, simulations performed on enzymatic reactions show that the transmission coefficient deviates only modestly from unity for reasonable choices of the reaction coordinate. Typical values usually range between 0.5 and 0.9 (refs 23,24,28), which reflects a modest participation of protein motions in the passage of the system over the barrier top.

A convenient way to analyse quantitatively the role of protein motions during the chemical transformation is to project the multidimensional free-energy surface (FES) of the enzymatic reaction in a 2D model obtained as a function of a solute coordinate and a solvent coordinate29,30. Such a FES allows one to estimate, at a quantitative level, the timing and coupling between the solute and solvent motions along the reaction path. Furthermore, the transmission coefficient can be derived from the differences between the dividing surface defined on the FES obtained under the equilibrium solvation approximation and that obtained on the 2D FES (which reflects the coupling between the solute and solvent coordinates).

A prototypical example in which solvent effects can play an important role is the SN2 reaction. Studies in aqueous solution have already established that the reaction proceeds with a large rearrangement of water molecules around the nucleophile and the leaving group30,31. Thus, the change in the electrostatic potential created by the environment on the nucleophile and the leaving group is a good coordinate on which to follow the evolution of the environment along the reaction process32, and the distances associated with the bonds to be broken and formed provide an adequate solute coordinate. In this work, we studied the SN2 nucleophilic reaction between dichloroethane and Asp124 in DhlA, a haloalkane dehalogenase from Xanthobacter autotrophicus (Fig. 1)33,34, an enzymatic reaction analysed theoretically by several groups35-38. We also modelled the counterpart SN2 reaction in aqueous solution, in which a molecule of acetate is employed as the nucleophile. The comparison between the enzymatic and in-solution process offers an excellent opportunity to analyse the

Trp175

GC'......h-nA

Trp125

H Asp124

Figure 1 | Schematic representation of the SN2 reaction catalysed by DhlA.

The QM subsystem is depicted by the atoms labelled in blue.

role of structure, flexibility and dynamics of the environment on a fundamental class of chemical reactions.

Results

The PM3/MM FESs that correspond to nucleophilic attack in the aqueous solution and in the active site of DhlA traced as a function of the solute (j) and solvent (s) coordinates are presented in Fig. 2a,b, respectively. The solute coordinate (j — d(CCl) — d(CO)) evolves from negative to positive values as the reaction proceeds. The solvent coordinate is obtained from the electrostatic potential created by the environment on the leaving group and the nucleo-phile (s — VCl — VO; see Methods) and evolves from negative values at the RS (where the nucleophile atom bearing the negative charge is stabilized by electrostatic interactions with the environment) to positive values at the product state (where the electrostatic potential takes larger positive values on the leaving group that now is negatively charged). This figure also displays, as continuous lines on the FESs, the minimum free-energy paths (MFEPs) obtained from the free-energy gradient. The free-energy paths obtained using the solute coordinate and assuming that the solvent coordinate equilibrates at each value of the former (equilibrium free-energy paths (EFEPs)) are also shown in Fig. 2 (dashed lines).

Analysis of FESs. The free-energy differences between the saddle points and the reactant minima located on the FESs are 27.4 and 37.1 kcal mol 1 for the catalysed and uncatalysed processes, respectively. The activation free energy deduced from the experimental rate constant of the enzymatic reaction at 298 K is 15.3 kcal mol 1 (ref. 34) and the barrier estimated for the process in aqueous solution is 26 kcal mol 1 (ref. 36) (a value of 29.9 kcal mol 1 was reported for the reaction in solution at 373K)35,39. The overestimation observed in our theoretical values results from the use of the PM3 Hamiltonian35,37. The M06-2X-corrected free-energy barriers (see Methods) for the enzymatic and in-solution processes are 16.5 and 27.4 kcal mol 1, respectively, in better quantitative agreement with the experimental values. In any case, the PM3/MM calculations provide a correct estimation of the catalytic effect, defined as the difference between the in-solution and the enzymatic free-energy barriers. The PM3/MM difference is 9.7 kcal mol—1, in good agreement with the difference derived from the M06-2X values, 10.9 kcal mol 1, and from the experimental values, 10.7 kcal mol 1.

The FESs show noticeable differences between the reaction in the solution and in the enzyme. For the RS the protein structure provides a much more adequate environment for the progress of the reaction than does the solution. In solution, the RS occurs at a value for the solvent coordinate of about — 100 kcal mol 1 |e| 1, and the enzymatic RS occurs at about — 25 kcal mol 1 |e| 1. This latter value of the solvent coordinate is much closer to the value needed to reach the TS

ARTICLES

NATURE CHEMISTRY doi: 10.1038/NCHEM.1660

-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0

a (A) a (A)

Figure 2 | FESs corresponding to the SN2 reaction between an acetate anion and dichloroethane. a,b, FESs for the reaction in aqueous solution (a) and in

the active site of DhlA (b). For a and b the solute distinguishing coordinate is the antisymmetric combination of the bond-breaking and -forming distances (j — d(ClC) — d(OC)); the solvent coordinate is the antisymmetric combination of the electrostatic potential created by the environment on the leaving group and the nucleophilic oxygen (s — Va — VO). The isoenergetic free-energy lines represent steps of 3 kcal mol—1. The continuous red and blue lines represent the minimum free-energy path on the FESs and the dashed ones correspond to the path obtained assuming equilibrium solvation along the solute coordinate. The white arrow represents the reaction path from the enzymatic Michaelis complex in a completely rigid protein environment.

(s ~ —10 kcal mol 1 |e| 1 in both environments). For the Michaelis complex the protein is already organized, from the electrostatic point of view, to favour the reaction, but in aqueous solution the environment largely needs to be reorganized to facilitate the reaction30'31'38.

Nevertheless, the protein structure does not behave as a rigid scaffold in which the reaction takes place. The reaction would be significantly more difficult in a frozen-protein environment in which the solvent coordinate remains unchanged from RS to TS. A straight line from the enzymatic RS at a constant value of the solvent coordinate represents the rigid-environment path. This path (white arrow in Fig. 2b) shows a PM3/MM free-energy barrier of 31.5 kcal mol 1, ~4 kcal mol 1 higher than the value observed for the flexible protein. A rate constant that corresponds to a hypothetical rigid protein would be about 103 times smaller, at room temperature, than that in the real enzyme. This means that protein flexibility plays a role in catalysis and that the environment needs to be rearranged when going from the Michaelis complex to the TS to reach a maximum reduction in the activation free energy.

Flexibility can be quantified by means of the force constant associated with the solvent coordinate. The force constants obtained from a parabolic fit of the free-energy change along the solvent coordinate in solution and in DhlA are given in Table 1. The force constant obtained in the enzyme is about 4.2 times larger than that obtained in the solution. The protein structure is stiffer than the structure of water, which is related to the existence of a network of covalent bonds in the former. However, as stated above, the change needed in the solvent coordinate to reach the TS from the RS is much smaller in the enzyme than in the solution. The final result is that the work to be done on the solvent coordinate to reach the TS is significantly smaller in the enzyme than in the solution. According to the MFEP traced on the FESs, the free-energy difference between the TS and the RS can be written as approximately

AG} « AGs,(sfS ^ s};jRS) + DGj ^ jnz;sf)

where i stands for the environment (enzyme or the water solution). The first term of the right-hand side of equation (3) represents the work to be done on the solvent coordinate and the second one the work to be done on the solute coordinate. The values that correspond to the solvent coordinate are about 11 and 3 kcal mol 1 in aqueous solution and in DhlA, respectively. Thus, the free-energy cost associated with the change along the solvent coordinate in the enzyme is substantially smaller than that in the solution, and this difference

represents 80% of the catalytic effect. It is interesting that, in the enzyme, the probability of sampling configurations of the environment conductive to the reaction (s*) is much larger than the probability of sampling adequate values of the solute coordinate (j*). The work associated with the reorganization of the environment in the enzymatic process represents only 11% of the total free-energy barrier.

Another important aspect to be analysed is the timing between the solute and solvent motions. According to the MFEPs obtained in the two environments, the change in the solvent coordinate precedes the change in the solute coordinate. From a dynamic point of view this means that the solvent coordinate is slower than the solute coordinate. A significant difference in the timescales of these two motions could result in important dynamical effects because of the delay between them. To characterize the time evolution of the environment we computed the characteristic frequencies associated with the motion along the solvent coordinate in the two environments using force constants and effective masses deduced from the equipartition principle. These frequencies are provided in Table 1. As observed, both the force constant and the effective mass associated with the solvent coordinate are larger in the enzyme than in the solution. Both effects cancel out and, as a result, the frequency associated with the motion along the solvent coordinate in the enzyme (410 cm 1) is very similar to the value obtained in the solution (480 cm 1). These values essentially correspond to the reorientation of hydrogen-bond donors around the nucleophile and the leaving group, motions that occur in picoseconds or faster38. So, from the dynamical point of view there are no significant differences in the participation of the environment during the reaction progress in aqueous solution or in DhlA. According to the MFEPs, the motions along the solvent coordinate, breaking of hydrogen bonds to the nucleophile and forming of new ones to the leaving group occur before or after the motion along the solute coordinate and the timescales associated

Table 1 | Force constants, effective masses and characteristic frequencies associated with the solvent coordinate (s) for the Sn2 reaction in solution and in DhlA.

Aqueous solution DhlA

Ks (kcal 2 1 mol |e|2) 3.3 x 10 23 ms (kcal 21 mol |e|2 s2) 0.40 x 10230 ns (cm21) 480 1.4x 1022 2.30 x 10230 410

Figure 3 | The TS ensembles obtained from equilibrium and non-equilibrium pictures differ slightly, as reflected in the respective dividing surfaces.

a, Dividing surfaces for the reaction in the solution. The arrow represents a hypothetical trajectory recrossing the equilibrium dividing surface. b, Dividing surfaces for the reaction in the active site of DhlA. The angle between the equilibrium and non-equilibrium plane is smaller than that in the solution. c, Schematic representation of the TS wells projected on the equilibrium and non-equilibrium dividing surfaces with their characteristic frequencies (vj and vs j). For a, b and c the continuous lines correspond to the dividing surface defined according to the FES traced along the solute and solvent coordinates. The dashed lines correspond to the dividing surface obtained when the solvent coordinate is assumed to be at equilibrium (the equilibrium plane is then defined just as j — j).

with these movements are very similar in the two media. Interestingly, we did not find any evidence that slow conformational motions of the protein affected the chemical step. Obviously, these motions can exist, but either they do not have consequences on the electrostatic coordinate (and thus on the energetics of the reaction) or they do not happen in the close neighbourhood of the Michaelis complex in the free-energy landscape.

Once analysed, we can compare the results obtained using both the solute and the solvent coordinates with those obtained assuming equilibrium solvation at any value of the solute coordinate. This is the usual approximation employed to analyse chemical reactions in condensed environments. The EFEPs, also presented in Fig. 2, go from the reactant to the product minima through the saddle point. Then, because the free energy is a state function, 1D profiles traced along the solute coordinate provide almost the same activation free energies as those of 2D surfaces obtained as a function of the solute and solvent coordinates (the origin of the small differences is discussed below). However, although the free-energy differences are correct, the equilibrium solvation approach is unable to describe properly the timing between the solute and solvent coordinates. Effectively, as observed in Fig. 2, in the equilibrium treatment one pulls along the solute coordinate and the solvent coordinate abruptly changes in the vicinity of the TS, but in the MFEPs, both in the solution and in the enzyme, solvent motions precede the changes along the solute coordinate. In any case, this limitation does not affect the estimation of the

reaction-rate constant, because this is determined mostly by the free-energy difference between the TS and the RS.

Evaluation of the transmission coefficient. Although the MFEPs and EFEPs coincide at the RS and TS, there is a small difference in the activation free energies estimated from the 2D (or non-equilibrium) and the 1D (or equilibrium) treatments. The origin of this difference is in the definition of the TS ensemble obtained in each treatment. As shown in Fig. 3, the dividing surface in the 2D description contains the saddle point and goes through the ridges that separate the reactant and product valleys. In the 1D treatment, only the solute coordinate is employed and thus the dividing surface is defined simply as j — j* (represented as dashed lines in Fig. 3). The difference corresponds to a rotation of the dividing surface in variational TST40. The TS further along the non-equilibrium dividing surface is narrower than the TS further along the equilibrium dividing surface (see Fig. 3c) and so the frequencies associated with the motion of the TS along the former dividing surface (vs j) are larger than those for the latter (vj). This can be translated into an entropic difference in the TS ensembles and consequently in the activation free energies29:

AAGJ = DGt(s,j)-DGt(j) = RT ln ^ (4)

From our FESs we estimated that the ratios between frequencies are about 1.3 and 1.1; using these values we estimated that AAG*

ARTICLES

NATURE CHEMISTRY doi: io.io38/nchem.i66o

at 298 K is about 0.2 kcal mol2 and 0.1 kcal mol2 in the solution and in the enzyme, respectively. These free-energy differences are quite small, below the statistical uncertainty of typical free-energy simulations, and thus non-equilibrium effects make a very small contribution to the activation free energies.

Obviously, a smaller activation free energy is translated into a larger rate constant. This difference should then be compensated by the consideration of a transmission coefficient smaller than unity (see equation (2)). The origin of this value can be understood in terms of the dividing surfaces obtained in the non-equilibrium and the equilibrium descriptions. The dividing surface obtained from the non-equilibrium treatment is defined such that any trajectory arriving at that surface from the reactant side will continue to the product region because the free energy continuously decreases in that direction, and so the transmission coefficient is equal to unity for that surface. However, using the equilibrium dividing surface, some trajectories that go from reactants to products find a free-energy barrier after crossing this surface and they could return to the reactant side (see Fig. 3a). This means that, in this description, the transmission coefficient would be lower than unity. Obviously, using equation (2) the final estimation of the rate constant obtained from the non-equilibrium and the equilibrium approaches would be the same if the transmission coefficient is obtained as29

AAG} k = e- RT

We believe that our analysis of this reaction can be extended easily to other enzymatic systems and so provide an adequate framework for a quantitative discussion on the role of protein motions in catalysis.

Methods

We used a quantum mechanics/molecular mechanics (QM/MM) computational scheme in which dichloroethane and the side chain of residue Asp 124 are described using the PM3 semi-empirical Hamiltonian41. The rest of the enzyme and/or water molecules form the MM subsystem described by means of the all-atoms optimized potential for liquid simulation (OPLS) for the enzyme42 and a flexible TIP3P potential for water molecules43. The Lennard-Jones parameters for the QM/MM interactions are also taken from the OPLS potential, except those for the chlorine atoms, which are taken from Gao and Xia44.

For the system in aqueous solution we placed dichloroethane and acetate in a pre-equilibrated box of water molecules of side 55.8 A, deleting all those water molecules with oxygen atoms found at <2.8 A from any non-hydrogen atom of the solute fragments. For the enzyme-substrate system the X-ray crystal structure coordinates were taken from the Protein Data Bank (code 2DHC)45. The protonation state of titratable residues was determined using the PropKa program46. The whole system was placed in a pre-equilibrated cubic box of water molecules of side 79.5 A. To neutralize the charge of the protein, 16 sodium ions were added so that both in the solution and in the enzymatic system the total charge was —1. The initial coordinates for the TSs in both environments are taken from our previous work37,38. The systems were equilibrated further by means of 200 ps of molecular dynamics simulation in the NVT ensemble at the reference temperature of 298 K using the Langevin integrator with a time step of 1 fs and periodic boundary conditions. A cutoff radii switched between 12.5 and 15 A was applied for all types of interaction.

In this work we obtained FESs using two different coordinates, a solute coordinate (j) and a solvent coordinate (s). The FES can be expressed as

Using the free-energy differences given above, the transmission coefficients for the reaction in the solution and in the enzyme obtained using the solute coordinate as the distinguishing reaction coordinate are 0.8 and 0.9, respectively. Obviously, this procedure leads to a very crude estimate of the transmission coefficient because of the statistical errors associated with the free energies and of the non-explicit treatment of all the degrees of freedom. A more accurate estimation shows that the transmission coefficients for this reaction are about 0.6 and 0.8 in aqueous solution and in the enzyme, respectively38. The transmission coefficients, evaluated at the TS, depend on the participation of the solvent coordinate in the barrier-crossing event, which is quite small for this reaction. However, it must be stressed again that this does not mean that the protein remains unchanged during the chemical step. Protein motions occur first and facilitate the motion along the solute coordinate, but at the TS they can be considered to be in equilibrium. In other words, although the reactions involve protein motions coupled to the chemical transformation, the probability that these motions take the system to the TS is determined mainly by the activation free energy1'5'15.

Discussion

The use of an explicit solvent coordinate can be a useful strategy to dissect the role of structure, flexibility and dynamics in catalysis. Comparison of the FESs obtained in the solution and in the enzyme for a prototypic SN2 reaction shows that the origin of the catalytic efficiency of the enzyme is due to a protein structure that provides an adequate environment for the progress of the reaction. However, the protein is not completely rigid and some motions contribute to reduce the free-energy barrier. These motions, which occur on a picosecond timescale or faster, dynamically precede the change in the solute coordinate. We also show that these protein motions can be treated as equilibrium fluctuations with the rate constant determined mainly by the free-energy difference between the TS and the RS. The activation free energies estimated under the assumption that the solvent coordinate is in equilibrium with the solute coordinate are only slightly underestimated with respect to those obtained in the non-equilibrium description. These differences can be accounted for by the inclusion of a transmission coefficient in the TST rate constant.

W(j, s) = C - kTln J p(xN)S(j(xN - j -

- s0)dxN

where p(xN) is the probability density of finding the system at the configuration xN. In this case, the solute coordinate (j) is the antisymmetric combination of the distances of the outgoing chloride and the incoming oxygen to the carbon atom (j — d(ClC) — d(OC)). The solvent coordinate selected was the antisymmetric combination of the electrostatic potential created by the environment on the outgoing chlorine atom and the incoming oxygen atom, the electrostatic potential created by the MM environment on the chlorine atom:

s = Vci(xN) - Vo(xN) = £ T

j=1 Ix;

;=i - xo|

where the sums run on the M sites of the environment with point charges This is a collective coordinate that involves all the MM atoms with an electrostatic influence on the donor or acceptor atoms.

The FESs that correspond to the reactions in aqueous solution and the enzyme were obtained using umbrella sampling47, applying parabolic constraints to the solute and the solvent coordinates:

V = 1 K(j-

1 , Vs = 2 *s(S

Molecular dynamics preferentially explores the most probable configurations of the system around the reference values j0 and s0. The joint probability distributions of the two coordinates were obtained by means of the weighted histogram analysis method (WHAM)48. To save computational cost, simulations were performed with any atom beyond 25 A of dichloroethane kept frozen. A total of 5,454 simulation windows consisting of 5 ps of equilibration and 50 ps of production were employed to trace the FESs in the solution, and for the enzyme 3,100 windows were needed. The force constants used to keep the system at the reference values of the solute and solvent coordinates were 2,500 kJ mol—1 A—2 and 0.01 kJ —1 mol |e|2, which provided a good control of the coordinates32.

Importantly, the PM3/MM Hamiltonian results in systematically overestimated energy barriers, but the geometries obtained for the RS and TS are good enough for reasonable estimations of kinetic isotope effects35. We then corrected the systematic error in the activation free energies by means of single-point calculations at higher theoretical levels. With this purpose we optimized ten TS structures starting from different configurations selected from the corresponding simulation. After intrinsic reaction coordinate calculation, the energy barrier was obtained at the PM3/MM level and by means of single-point calculations at the M06-2X/6-311 + G(2df,2p)/MM level49. The correction energy term was evaluated as the averaged difference between the semi-empirical and M06-2X energy barriers.

Received 15 January 2013; accepted 18 April 2013; published online 26 May 2013

r)ô(s(xN )

References

1. Hammes, G. G., Benkovic, S. J. & Hammes-Schiffer, S. Flexibility, diversity, and cooperativity: pillars of enzyme catalysis. Biochemistry 50, 10422-10430 (2011).

2. Henzler-Wildman, K. A. et al. A hierarchy of timescales in protein dynamics is linked to enzyme catalysis. Nature 450, 913-916 (2007).

3. Benkovic, S. J., Hammes, G. G. & Hammes-Schiffer, S. Free-energy landscape of enzyme catalysis. Biochemistry 47, 3317-3321 (2008).

4. Henzler-Wildman, K. & Kern, D. Dynamic personalities of proteins. Nature 450, 964-972 (2007).

5. Garcia-Viloca, M., Gao, J., Karplus, M. & Truhlar, D. G. How enzymes work: analysis by modern rate theory and computer simulations. Science 303, 186-195 (2004).

6. Gao, J. et al. Mechanisms and free energies of enzymatic reactions. Chem. Rev. 106, 3188-3209 (2006).

7. Olsson, M. H. M., Parson, W. W. & Warshel, A. Dynamical contributions to enzyme catalysis: critical tests of a popular hypothesis. Chem. Rev. 106, 1737-1756 (2006).

8. Antoniou, D., Basner, J., Nunez, S. & Schwartz, S. D. Computational and theoretical methods to explore the relation between enzyme dynamics and catalysis. Chem. Rev. 106, 3170-3187 (2006).

9. Nashine, V. C., Hammes-Schiffer, S. & Benkovic, S. J. Coupled motions in enzyme catalysis. Curr. Opin. Chem. Biol. 14, 644-651 (2010).

10. Ramanathan, A. & Agarwal, P. K. Evolutionarily conserved linkage between enzyme fold, flexibility, and catalysis. PLoS Biol 9, e1001193 (2011).

11. Zhang, J. & Klinman, J. P. Enzymatic methyl transfer: role of an active site residue in generating active site compaction that correlates with catalytic efficiency. J. Am. Chem. Soc. 133, 17134-17137 (2011).

12. Bhabha, G. et al. A dynamic knockout reveals that conformational fluctuations influence the chemical step of enzyme catalysis. Science 332, 234-238 (2011).

13. Adamczyk, A. J., Cao, J., Kamerlin, S. C. L. & Warshel, A. Catalysis by dihydrofolate reductase and other enzymes arises from electrostatic preorganization, not conformational motions. Proc. Natl Acad. Sci. USA 108, 14115-14120 (2011).

14. Glowacki, D. R., Harvey, J. N. & Mulholland, A. J. Taking Ockham's razor to enzyme dynamics and catalysis. Nature Chem. 4, 169-176 (2012).

15. Kamerlin, S. C. L. & Warshel, A. At the dawn of the 21st century: is dynamics the missing link for understanding enzyme catalysis? Proteins 78,1339-1375 (2010).

16. Warshel, A. et al. Electrostatic basis for enzyme catalysis. Chem. Rev. 106, 3210-3235 (2006).

17. Kurplus, M & McCammon, J A. Dynamics of proteins: elements and function. Annu. Rev. Biochem. 52, 263-300 (1983).

18. Marcus, R. A. Chemical and electrochemical electron-transfer theory. Annu. Rev. Phys. Chem. 15, 155-196 (1964).

19. Kosugi, T. & Hayashi, S. Crucial role of protein flexibility in formation of a stable reaction transition state in an a-amylase catalysis. J. Am. Chem. Soc. 134, 7045-7055 (2012).

20. Pang, J., Pu, J., Gao, J., Truhlar, D. G. & Allemann, R. K. Hydride Transfer reaction catalyzed by hyperthermophilic dihydrofolate reductase is dominated by quantum mechanical tunneling and is promoted by both inter- and intramonomeric correlated motions. J. Am. Chem. Soc. 128, 8015-8023 (2006).

21. Kanaan, N. et al. Temperature dependence of the kinetic isotope effects in thymidylate synthase. A theoretical study. J. Am. Chem. Soc. 133, 6692-6702 (2011).

22. Hay, S. & Scrutton, N. S. Good vibrations in enzyme-catalysed reactions. Nature Chem. 4, 161-168 (2012).

23. Pu, J., Gao, J. & Truhlar, D. G. Multidimensional tunneling, recrossing, and the transmission coefficient for enzymatic reactions. Chem. Rev. 106, 3140-3169 (2006).

24. Boekelheide, N., Salomon-Ferrer, R. & Miller, T. F. Dynamics and dissipation in enzyme catalysis. Proc. Natl Acad. Sci. USA 108, 16159-16163 (2011).

25. Garrett, B. C. & Truhlar, D. G. in Theory and Applications of Computational Chemistry (eds Dykstra, C. E., Frenking, G., Kim, K. S. & Scuseria, G. E.) 67-87 (Elsevier, 2005).

26. Grote, R. F. & Hynes, J. T. The stable states picture of chemical reactions. II. Rate constants for condensed and gas phase reaction models. J. Chem. Phys. 73, 2715-2732 (1980).

27. Bergsma, J. P., Gertner, B. J., Wilson, K. R. & Hynes, J. T. Molecular dynamics of a model SN2 reaction in water. J. Chem. Phys. 86, 1356-1376 (1987).

28. Ruiz-Pernia, J. J., Tunon, I., Moliner, V., Hynes, J. T. & Roca, M. Dynamic effects on reaction rates in a Michael addition catalyzed by chalcone isomerase. Beyond the frozen environment approach. J. Am. Chem. Soc. 130, 7477-7488 (2008).

29. Gertner, B. J., Bergsma, J. P., Wilson, K. R., Lee, S. & Hynes, J. T. Nonadiabatic solvation model for SN2 reactions in polar solvents. J. Chem. Phys. 86, 1377-1386 (1987). N

30. Hwang, J. K., King, G., Creighton, S. & Warshel, A. Simulation of free energy relationships and dynamics of Sn2 reactions in aqueous solution. J. Am. Chem. Soc. 110, 5297-5311 (1988).

31. Gertner, B. J., Wilson, K. R. & Hynes, J. T. Nonequilibrium solvation effects on reaction rates for model Sn2 reactions in water. J. Chem. Phys. 90, 3537-3558 (1989).

32. Ruiz-Pernia, J. J., Martí, S., Moliner, V. & Tuñón, I. A novel strategy to study electrostatic effects in chemical reactions: differences between the role of solvent and the active site of chalcone isomerase in a Michael addition. J. Chem. Theory Comput. 8, 1532-1535 (2012).

33. Janssen, D. B., Scheper, A., Dijkhuizen, L. & Witholt, B. Degradation of halogenated aliphatic-compounds by Xanthobacter autotrophicus GJ10. Appl. Environ. Microbiol. 49, 673-677 (1985).

34. Schanstra, J. P., Kingma, J. & Janssen, D. B. Specificity and kinetics of haloalkane dehalogenase. J. Biol. Chem. 271, 14747-14753 (1996).

35. Devi-Kesavan, L. S. & Gao, J. Combined QM/MM study of the mechanism and kinetic isotope effect of the nucleophilic substitution reaction in haloalkane dehalogenase. J. Am. Chem. Soc. 125, 1532-1540 (2003).

36. Shurki, A., Strajbl, M., Villa, J. & Warshel, A. How much do enzymes really gain by restraining their reacting fragments? J. Am. Chem. Soc. 124, 4097-4107 (2002).

37. Soriano, A. et al. Electrostatic effects in enzyme catalysis: a quantum mechanics/molecular mechanics study of the nucleophilic substitution reaction in haloalkane dehalogenase. Theor. Chem. Accounts 112, 327-334 (2004).

38. Soriano, A., Silla, E., Tunon, I. & Ruiz-Lopez, M. F. Dynamic and electrostatic effects in enzymatic processes. An analysis of the nucleophilic substitution reaction in haloalkane dehalogenase. J. Am. Chem. Soc. 127, 1946-1957 (2005).

39. Okamoto, K., Kita, T., Araki, K. & Shingu, H. Kinetic studies of bimolecular nucleophilic substitution. 4. Rates of Sn2 and E2 reactions of beta-substituted ethyl chlorides with sodium acetate in aqueous solutions. B. Chem. Soc. Jpn 40, 1913 (1967).

40. Truhlar, D. G. & Garrett, B. C. Variational transition-state theory. Acc. Chem. Res. 13, 440-448 (1980).

41. Stewart, J. J. P. Optimization of parameters for semiempirical methods I. Method. J. Comput. Chem. 10, 209-220 (1989).

42. Jorgensen, W. L., Maxwell, D. S. & Tirado-Rives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 118, 11225-11236 (1996).

43. Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W. & Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926-935 (1983).

44. Gao, J. & Xia, X. A two-dimensional energy surface for a type II Sn2 reaction in aqueous solution. J. Am. Chem. Soc. 115, 9667-9675 (1993).

45. Verschueren, K. H., Seljée, F., Rozeboom, H. J., Kalk, K. H. & Dijkstra, B. W. Crystallographic analysis of the catalytic mechanism of haloalkane dehalogenase. Nature 363, 693-698 (1993).

46. Olsson, M. H. M., Sondergaard, C. R., Rostkowski, M. & Jensen, J. H. PROPKA3: consistent treatment of internal and surface residues in empirical pKa predictions. J. Chem. Theory Comput. 7, 525-537 (2011).

47. Torrie, G. M. & Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: umbrella sampling. J. Comput. Phys. 23, 187-199 (1977).

48. Kumar, S., Rosenberg, J. M., Bouzida, D., Swendsen, R. H. & Kollman, P. A. The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method. J. Comput. Chem. 13, 1011-1021 (1992).

49. Zhao, Y. & Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Accounts 120, 215-241 (2008).

Acknowledgements

The authors acknowledge financial support from the Ministerio de Economía y Competitividad (MEC) through project CTQ2012-36253-C03. J.J.R-P. thanks a Juan de la Cierva contract and R.G-M. a FPU fellowship of the Ministerio de Economía y Competitividad. I.T. acknowledges helpful discussions held with D. Laage and J. T. Hynes during his sabbatical stay at the Eé cole Normale Supéerieure, France. The authors acknowledge computational facilities of the Servei d'Informatica de la Universitat de Valencia on the 'Tirant' supercomputer.

Author contributions

I.T., V.M. and J.J.R-P. designed the computational experiments. S.M. wrote the code and R.G-M. performed the calculations. I.T., V.M. and J.J.R-P. co-wrote the first version ofthe paper. All the authors commented and discussed the results and the final version ofthe manuscript.

Additional information

Reprints and permissions information is available online at www.nature.com/reprints. Correspondence and requests for materials should be addressed to J.J.R-P. and I.T.

Competing financial interests

The authors declare no competing financial interests.