# Dynamics and dissipation in enzyme catalysis

See allHide authors and affiliations

Edited by Donald G. Truhlar, University of Minnesota, Minneapolis, MN, and approved August 5, 2011 (received for review April 21, 2011)

## Abstract

We use quantized molecular dynamics simulations to characterize the role of enzyme vibrations in facilitating dihydrofolate reductase hydride transfer. By sampling the full ensemble of reactive trajectories, we are able to quantify and distinguish between statistical and dynamical correlations in the enzyme motion. We demonstrate the existence of nonequilibrium dynamical coupling between protein residues and the hydride tunneling reaction, and we characterize the spatial and temporal extent of these dynamical effects. Unlike statistical correlations, which give rise to nanometer-scale coupling between distal protein residues and the intrinsic reaction, dynamical correlations vanish at distances beyond 4–6 Å from the transferring hydride. This work finds a minimal role for nonlocal vibrational dynamics in enzyme catalysis, and it supports a model in which nanometer-scale protein fluctuations statistically modulate—or gate—the barrier for the intrinsic reaction.

Protein motions are central to enzyme catalysis, with conformational changes on the micro- and millisecond timescale well-established to govern progress along the catalytic cycle (1, 2). Less is known about the role of faster, atomic-scale fluctuations that occur in the protein environment of the active site. The textbook view of enzyme-catalyzed reaction mechanisms neglects the functional role of such fluctuations and describes a static protein environment that both scaffolds the active-site region and reduces the reaction barrier (3). This view has grown controversial amid evidence that active-site chemistry is coupled to motions in the enzyme (4–6), and it has been explicitly challenged by recent proposals that enzyme-catalyzed reactions are driven by vibrational excitations that channel energy into the intrinsic reaction coordinate (7, 8) or promote reactive tunneling (9, 10). In the following, we combine quantized molecular dynamics and rare-event sampling methods to determine the mechanism by which protein motions couple to reactive tunneling in dihydrofolate reductase and to clarify the role of nonequilibrium vibrational dynamics in enzyme catalysis.

Manifestations of enzyme motion include both statistical and dynamical correlations. Statistical correlations are properties of the equilibrium ensemble and describe, for example, the degree to which fluctuations in the spatial position of one atom are influenced by fluctuations in another; these correlations govern the free-energy (FE) landscape and determine the transition state theory kinetics of the system (6). Dynamical correlations are properties of the time-evolution of the system and describe coupling between inertial atomic motions, as in a collective vibrational mode. Compelling evidence for long-ranged (i.e., nanometer-scale) networks of statistical correlations in enzymes emerges from genomic analysis (11), molecular dynamics simulations (11–13), and kinetic studies of double-mutant enzymes (14–16). But the role of dynamical correlations in enzyme catalysis remains unresolved (4, 5, 7, 17, 18), with experimental and theoretical results suggesting that the intrinsic reaction is activated by vibrational modes involving the enzyme active site (9, 19, 20) and more distant protein residues (7, 8, 21). The degree to which enzyme-catalyzed reactions are coupled to the surrounding protein environment, and the lengthscales and timescales over which such couplings persist, are central questions in the understanding, regulation, and de novo design of biological catalysts (22).

*Escherichia coli* dihydrofolate reductase (DHFR) is an extensively studied prototype for protein motions in enzyme catalysis. It catalyzes reduction of the 7,8-dihydrofolate (DHF) substrate via hydride transfer from the nicotinamide adenine dinucleotide phosphate (NADPH) cofactor (Fig. 1*A* and Fig. S1). We investigate this intrinsic reaction using ring polymer molecular dynamics (RPMD) (23, 24), a recently developed path-integral method that enables inclusion of nuclear quantization effects, such as the zero-point energy and tunneling, in the dynamics of the transferring hydride. RPMD simulations with over 14,000 atoms are performed using explicit solvent and using an empirical valence bond (EVB) potential to describe the potential energy surface for the transferring hydride; the EVB potential is obtained from an effective Hamiltonian matrix, with diagonal elements (*V*_{r}(**x**) and *V*_{p}(**x**)) corresponding to the potential energy for the reactant and product bonding connectivities and with the constant off-diagonal matrix element fit to the experimental rate (11, 25). The vector **x** includes the position of the quantized hydride and all classical nuclei in the system. The EVB potential employed here was chosen for consistency with earlier simulation studies of DHFR (11, 26); although the mechanistic issues addressed in this study are not expected to be highly sensitive to the details of the potential energy surface, it should be noted that more accurate electronic structure theory methods for describing the enzyme potential energy surface, including the combined quantum mechanical and classical mechanical (QM/MM) method, are available and widely used in enzyme simulations.

The thermal reaction rate is calculated from the product of the Boltzmann-weighted activation FE and the reaction transmission coefficient (24), both of which are calculated in terms of the dividing surface *λ*(**x**) = -4.8 kcal/mol where *λ*(**x**) = *V*_{r}(**x**) - *V*_{p}(**x**). The FE surface *F*(*λ*) is obtained using over 120 ns of RPMD sampling (Fig. 1*B*), and the transmission coefficient is obtained from over 5,000 RPMD trajectories that are released from the Boltzmann distribution constrained to the dividing surface (Fig. 1*C*). In contrast to mixed quantum-classical and transition state theory methods, RPMD yields reaction rates and mechanisms that are formally independent of the choice of dividing surface or any other reaction coordinate assumption (24). Furthermore, the RPMD method enables generation of the ensemble of reactive, quantized molecular dynamics trajectories, which is essential for the following analysis of dynamical correlations. Calculation details, including a description of the rare-event sampling methodology used to generate the unbiased ensemble of reactive trajectories (27, 28), are provided in *Materials and Methods* below.

## Results and Discussion

The time-dependence of the transmission coefficient in Fig. 1*C* confirms that reactive trajectories commit to the reactant or product basins within 25 fs. The near-unity value of this transmission coefficient at long times indicates that recrossing of the dividing surface in reactive trajectories is a modest effect, although it is fully accounted for in this study, and it confirms that the collective variable *λ*(**x**) provides a good measure of progress along the intrinsic reaction. We find that quantization of the hydride lowers the FE barrier by approximately 3.5 kcal/mol (Fig. S2), in agreement with earlier work (29, 30).

Statistical correlations among the protein and enzyme active-site coordinates are shown in Fig. 2*A*. The normalized covariance among atom position fluctuations, *c*_{ij} = *C*_{ij}/(*C*_{ii}*C*_{jj})^{1/2} such that [1]is plotted for the Boltzmann distribution in the reactant, dividing surface, and product regions. The figure shows correlations among the protein *α*-carbons and the heavy atoms of the substrate and cofactor; the corresponding all-atom correlation plots are provided in Fig. S3. As has been previously and correctly emphasized (11, 29, 30), structural fluctuations in the active-site and distal protein residues are richly correlated within each region, which contributes to nonadditive effects in the kinetics of DHFR mutants (14, 31). Furthermore, the network of correlations varies among the three ensembles, indicating that fluctuations in distal protein residues respond to the adiabatic progress of the hydride from reactant to product. However, these time-averaged quantities do not address the role of dynamical correlations between the transferring hydride and its environment, which depend on the hierarchy of timescales for motion in the system.

To characterize dynamical correlations in the intrinsic reaction, we introduce a measure of velocity cross-correlations in the reactive trajectories, *d*_{ij}(*t*) = *D*_{ij}(*t*)/(*D*_{ii}(*t*)*D*_{jj}(*t*))^{1/2} such that [2]Here, denotes an average over the nonequilibrium ensemble of phase-space points that lie on reactive trajectories that crossed the dividing surface some time *t* earlier and subsequently terminate in the product basin. This quantity, which vanishes for the equilibrium ensemble, reports on the degree to which atoms move in concert during the intrinsic reaction step. Fig. 2 *B*–*D* show *d*_{ij}(*t*) for several atomic pairs in the active site. Negative dynamical correlations are seen between the donor and acceptor C atoms (Fig. 2*B*), which move in opposite directions (first approaching each other, then moving apart) during the hydride transfer. Similarly, positive correlations are seen between atom pairs on the cofactor (Fig. 2*C*) and on the substrate (Fig. 2*D*) that move in concert as the hydride is transferred. In each case, the primary features of the correlation decay within *τ* = 100 fs.

Fig. 2*E* summarizes the extent of dynamical correlations throughout the enzyme system in terms of . Only atoms in the substrate and cofactor regions (Fig. 2*E*, lower triangle) and a small number of protein atoms in the active-site region exhibit appreciable signal. The same conclusions are reached upon integrating the absolute value of the *d*_{ij}(*t*) (Fig. S4), emphasizing that this lack of signal in the protein residues is not simply due to the time integration. Instead, Fig. 2*E* reveals that the dynamical correlations between distal protein residues and the intrinsic reaction do not exist on any timescale. We also provide measures for dynamical correlations that are nonlocal in time (Fig. S5) and for dynamical correlations among perpendicular motions (Fig. S6), but the following conclusion is unchanged. The extensive network of statistical correlations (Fig. 2*A*) is neither indicative of, nor accompanied by, an extensive network of dynamical correlations during the intrinsic reaction (Fig. 2*E*).

A combined measure of the dynamical correlation between a given atom and the intrinsic reaction event can be obtained from the nonequilibrium ensemble average of velocities in the reactive trajectories. Specifically, we consider , where *ξ*∈{x,y,z} indicates the component of the velocity, the filter selects configurations in the region of the dividing surface, and is the average magnitude of *λ*(**x**) in the reactant and product regions. Each component of **f**_{i}(*t*) vanishes trivially at equilibrium. Fig. 3 *A*–*C* presents the measure for various atoms in the active-site region. The donor and acceptor C atoms (Fig. 3 *A* and *B*) are both strongly correlated with the dynamics of the intrinsic reaction, whereas the O atom in the Y100 residue of the active site (Fig. 3*C*) reveals smaller, but nonzero, signatures of dynamical correlation. Fig. 3*D* presents for each atom, summarizing the degree to which all atoms in the active site exhibit dynamical correlations, and Fig. 3*E* compares the correlation lengthscales in the enzyme. The main panel in Fig. 3*E* presents *f*_{i} as a function of the distance of heavy atoms from the midpoint of the hydride donor and acceptor, and the inset similarly presents the distance dependence of the statistical correlation measure , where *c*_{ij} is defined previously and where indices *μ* and *ν* label the donor and acceptor carbon atoms, respectively. Whereas the statistical correlations reach the nanometer lengthscale and involve the protein environment, dynamical correlations are extremely local in nature and primarily confined to the enzyme substrate and cofactor.

Fig. 4 illustrates that dynamical correlations in the intrinsic reaction are limited by disparities in the relative timescales for enzyme motion. The figure presents two-dimensional projections of the FE surface, *F*(*λ*,Θ_{α}), where *α*∈{1,2}, Θ_{1}(**x**) is the distance between hydride donor and acceptor atoms, and Θ_{2}(**x**) is the separation between active-site protein atoms I14 *C*_{δ} and Y100 O (side chain). Overlaid on the surfaces are the minimum FE pathway between the reactant and product basins, *s*, and the time-parameterized pathway followed by the ensemble of reactive trajectories, . Nonzero slope in *s* indicates statistical correlation of Θ_{α} with *λ*, whereas the same feature in *σ* indicates that the dynamics of Θ_{α} and *λ* are dynamically correlated. Fig. 4*A* confirms that the donor-acceptor distance is both statistically and dynamically correlated with the intrinsic reaction. In contrast, Fig. 4*B* reveals significant statistical correlation between Θ_{2} and the intrinsic reaction, but the reactive trajectories traverse the dividing surface region on a timescale that is too fast to dynamically couple to the protein coordinate.

The results presented here complement previous theoretical efforts to illuminate the role of protein motions in enzyme catalysis. For example, Neria and Karplus (32) used transmission coefficient calculations and constrained molecular dynamics (MD) trajectories to determine that the protein environment in triosephosphate isomerase (TIM) is essentially rigid (i.e., dynamically unresponsive) on the timescale of the intrinsic reaction dynamics; this finding is consistent with the lack of long-lengthscale dynamical correlations reported in the current study. Furthermore, Truhlar and coworkers (33, 34) and Karplus and Cui (35) both demonstrated that quasi-classical tunneling coefficients for hydrogen transfer evaluated at instantaneous enzyme configurations in the transition state region fluctuate significantly with donor-acceptor motions and other local active-site vibrations, which is likely consistent with the direct observation of short-lengthscale dynamical correlations reported here. However, by using quantized molecular dynamics to sample the ensemble of reactive trajectories in DHFR catalysis and to perform nonequilibrium ensemble averages that directly probe dynamical correlation, we provide a framework for strengthening and generalizing these earlier analyses. In particular, the current approach avoids transition state theory approximations by providing a rigorous statistical mechanical treatment of the ensemble of reactive trajectories, it allows for the natural characterization of lengthscales and timescales over which dynamical correlations persist, and it seamlessly incorporates dynamical effects due to both nuclear quantization and trajectory recrossing. We expect this approach to prove useful in future studies of dynamics in other enzymes, which will be necessary to confirm the generality of the conclusions drawn here.

## Concluding Remarks

The physical picture that emerges from this analysis is one in which the intrinsic reaction involves a small, localized group of atoms that are dynamically uncoupled from motions in the surrounding protein environment. As in the canonical theories for electron and proton transfer in the condensed phase (36, 37), the strongly dissipative protein environment merely gates the fast dynamics in the active site. This view reconciles the observed kinetic effects of distal enzyme mutations (14, 15) with evidence for short-ranged dynamical correlations in the active site (9, 19, 20), and it supports a unifying theoretical perspective in which slow, thermal fluctuations in the protein modulate the instantaneous rate for the intrinsic reaction (4, 7, 17, 29, 30, 32, 33, 34, 35, 38, 39, 40, 41). Furthermore, the work presented here provides clear evidence against the proposed role of nonlocal, rate-promoting vibrational dynamics in the enzyme (8), and it reveals the strikingly short lengthscales over which nonequilibrium protein dynamics couples to the intrinsic reaction. The combination of quantized molecular dynamics methods (24) with trajectory sampling methods (28) provides a useful approach for characterizing mechanistic features of reactions involving significant quantum tunneling effects. These findings for the case of the DHFR enzyme, a good candidate for dynamical correlations because of its small size and strong network of statistical correlations, suggest that nonlocal dynamical correlations are neither a critical feature of enzyme catalysis nor an essential consideration in de novo enzyme design.

## Materials and Methods

### Calculation Details.

All simulations are performed using a modified version of the Gromacs-4.0.7 molecular dynamics package (42). Further calculation details regarding the potential energy surface, the system initialization and equilibration protocol, free-energy sampling, and dividing surface sampling are provided in *SI Text*, Figs. S7–S9, and Table S1.

### Ring Polymer Molecular Dynamics.

The RPMD equations of motion (23) used to simulate the dynamics of DHFR are [3][4]where *U*(**q**,**Q**_{1},…,**Q**_{N}) is the potential energy function for the system, *n* = 32 is the number of ring polymer beads used to quantize the hydride, **q**_{j} and *m*_{n} are the position and mass of the *j*th ring polymer bead, and **q**_{0} = **q**_{n}. Similarly, *N* is the number of classical nuclei in the system, and **Q**_{k} and *M*_{k} are the position and mass of the *k*th classical atom, respectively. The interbead force constant is *k*_{n} = *m*_{H}*n*^{2}/(*βℏ*)^{2}, where *m*_{H} = 1.008 amu is the mass of the hydride and *β* = (*k*_{B}*T*)^{-1} is the reciprocal temperature; a temperature of *T* = 300 K is used throughout the study. For dynamical trajectories, RPMD prescribes that *m*_{n} = *m*_{H}/*n*.

### Calculating the Statistical Correlation Functions, *c*_{ij}.

In Fig. 3*A*, equilibrium ensemble averages are presented for the system in the reactant region, the dividing surface region, and the product region. These ensemble averages are strictly defined using [5]where *P*(*λ*) = exp(-*βF*(*λ*))/∫*dλ*^{′} exp(-*βF*(*λ*^{′})), and *F*(*λ*) is calculated using umbrella sampling, as described in *SI Text*. For the ensembles in the reactant, dividing surface, and product regions, we employ *λ*^{∗} = -181 kcal/mol, -7 kcal/mol, and 169 kcal/mol, respectively, and *δλ* = 2.5 kcal/mol.

### The Transition Path Ensemble.

Reactive trajectories are generated through forward- and backward-integration of initial configurations drawn from the dividing surface ensemble with initial velocities drawn from the Maxwell–Boltzmann distribution. Reactive trajectories correspond to those for which forward- and backward-integrated half-trajectories terminate on opposite sides of the dividing surface. From the 10,500 half-trajectories that are initialized on the dividing surface (i.e., 5,250 possible reactive trajectories), over 3,000 reactive RPMD trajectories are obtained. For analysis purposes, the integration of these reactive trajectories was continued for a total length of 1 ps in both the forward and backward trajectories.

The reactive trajectories that are initialized from the equilibrium Boltzmann distribution on the dividing surface must be reweighted to obtain the unbiased ensemble of reactive trajectories (i.e., the transition path ensemble) (27, 28, 43). A weighting term is applied to each trajectory α, correctly accounting for the recrossing and for the fact that the trajectories are performed in the microcanonical ensemble (27), [6]where the sum includes all instances in which trajectory α crosses the dividing surface, and is the velocity in the collective variable at crossing event *i*. We find that the relative statistical weight of all reactive trajectories that recross the dividing surface is 1.6%, emphasizing that recrossing does not play a large role in the current study.

## Acknowledgments

This work was supported by the National Science Foundation (NSF) CAREER Award (CHE-1057112) and computing resources at the National Energy Research Scientific Computing Center. Additionally, N.B. acknowledges an NSF graduate research fellowship, and T.F.M. acknowledges an Alfred P. Sloan Foundation fellowship.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: tfm{at}caltech.edu.

Author contributions: N.B., R.S.-F., and T.F.M. designed research, performed research, contributed new reagents/analytic tools, analyzed data, and wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1106397108/-/DCSupplemental.

Freely available online through the PNAS open access option.

## References

- ↵
- ↵
- Boehr DD,
- McElheny D,
- Dyson HJ,
- Wright PE

- ↵
- Cook PF,
- Cleland WW

- ↵
- ↵
- ↵
- Garcia-Viloca M,
- Gao J,
- Karplus M,
- Truhlar DG

- ↵
- ↵
- ↵
- Masgrau L,
- et al.

- ↵
- ↵
- Agarwal PK,
- Billeter SR,
- Rajagopalan PTR,
- Benkovic SJ,
- Hammes-Schiffer S

- ↵
- Pang JY,
- Pu JZ,
- Gao JL,
- Truhlar DG,
- Allemann RK

- ↵
- Rod TH,
- Radkiewicz JL,
- Brooks CL

- ↵
- Wong KF,
- Selzer T,
- Benkovic SJ,
- Hammes-Schiffer S

- ↵
- ↵
- Gira B,
- et al.

- ↵
- Benkovic SJ,
- Hammes-Schiffer S

- ↵
- ↵
- ↵
- ↵
- Saen-Oon S,
- Quaytman-Machleder S,
- Schramm VL,
- Schwartz SD

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Marcus RA,
- Sutin N

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Chemistry

- Biological Sciences
- Biophysics and Computational Biology