Computer Simulation of Primary Kinetic Isotope Effects in the Proposed Rate-limiting Step of the Glyoxalase I Catalyzed Reaction*

The proposed rate-limiting step of the glyoxalase I catalyzed reaction is the proton abstraction from the C1 carbon of the substrate by Glu172. Here we examine primary kinetic isotope effects and the influence of quantum dynamics on this process by computer simulations. The calculations utilize the empirical valence bond method in combination with the molecular dynamics free energy perturbation technique and path integral simulations. For the enzyme-catalyzed reaction a H/D kinetic isotope effect of 5.0 ± 1.3 is predicted in reasonable agreement with the experimental result of about 3. Furthermore, the magnitude of quantum mechanical effects is found to be very similar for the enzyme reaction and the corresponding uncatalyzed process in solution, in agreement with other studies. The problems associated with attaining the required accuracy in order for the present approach to be useful as a diagnostic tool for the study of enzyme reactions are also discussed.

The proposed rate-limiting step of the glyoxalase I catalyzed reaction is the proton abstraction from the C1 carbon of the substrate by Glu 172 . Here we examine primary kinetic isotope effects and the influence of quantum dynamics on this process by computer simulations. The calculations utilize the empirical valence bond method in combination with the molecular dynamics free energy perturbation technique and path integral simulations. For the enzyme-catalyzed reaction a H/D kinetic isotope effect of 5.0 ؎ 1.3 is predicted in reasonable agreement with the experimental result of about 3. Furthermore, the magnitude of quantum mechanical effects is found to be very similar for the enzyme reaction and the corresponding uncatalyzed process in solution, in agreement with other studies. The problems associated with attaining the required accuracy in order for the present approach to be useful as a diagnostic tool for the study of enzyme reactions are also discussed.
Glyoxalase I (GlxI) 1 is the first enzyme of the glyoxalase system, which catalyzes the transformation of toxic 2-oxoaldehydes, such as methylglyoxal, into the corresponding nontoxic ␣-hydroxycarboxylic acid. A glutathione (GSH) molecule and a divalent metal ion are required as cofactors in GlxI. Upon binding, methylglyoxal and GSH spontaneously form a hemithioacetal, which is the actual substrate of GlxI. In the reaction, this hemithioacetal is converted to S-D-lactoylglutathione. Subsequently, glyoxalase II hydrolyzes S-D-lactoylglutathione into D-lactic acid and GSH ( Fig. 1) (1).
The reaction mechanism catalyzed by GlxI has been proposed to proceed via an unstable enediol(-ate) intermediate formed through proton abstraction from the C1 carbon of the substrate by an amino acid residue acting as a base (2). This type of mechanism is analogous to that of several other enzymes catalyzing reactions of carbon acids, such as triosephosphate isomerase (3), enolase (4), and mandelate racemase (5). It is thus of interest to elucidate possible common features and principles of enzyme catalysis for this class of reactions. In GlxI a primary deuterium kinetic isotope effect of about 3 with methylglyoxal and phenylglyoxal substrates has been interpreted such that this proton transfer is the rate-limiting step (1). The recently solved crystal structure of GlxI in complex with the transition state analogue S-(N-hydroxy-N-p-iodophenylcarbamoyl)glutathione (6) shows a homodimer with two identical active sites in which a Zn 2ϩ ion is coordinated by two glutamate residues (Glu 172 and Glu 99 ), a histidine (His 26 ), a glutamine (Gln 33 ), and the two carbamoyl oxygens of S-(Nhydroxy-N-p-iodophenylcarbamoyl)glutathione. This structure implies that the enediolate intermediate is in cis configuration with both its oxygens pointing toward the metal. Based on structural analysis and site-directed mutagenesis studies, it seems that Glu 172 is the proton-abstracting base for the reaction with the S-form of the substrate (2, 6).
Recent empirical valence bond (EVB) simulation studies by us on the presumed rate-limiting step of the GlxI reaction with the S enantiomer of the methylglyoxal GSH hemithioacetal substrate and Glu 172 as the proton-abstracting base gave an activation free energy in good agreement with kinetic data (7). The calculations showed, also in agreement with experiments, that the enzyme is not extremely sensitive to the divalent metal ion species bound to the active site. It was found that Mg 2ϩ , Zn 2ϩ , and Ca 2ϩ are able to promote catalysis, although the latter ion yields a ϳ1 kcal/mol higher activation barrier than the two former. Furthermore, these simulations clearly showed that the active site metal ion is the key catalytic factor in GlxI and that it provides the main part of the essential stabilization of the enediolate intermediate (7).
Because H/D/T kinetic isotope effect (KIE) measurements are one of the major experimental probes of reaction mechanisms involving hydrogen transfer (8), it would be very helpful to have reliable theoretical methods for predicting such KIEs in enzyme reactions. Whereas ab initio methods are useful for calculation of KIEs in the gas phase, solution and enzymecatalyzed reactions present major challenges in this respect. However, significant progress in theoretical modeling of enzyme isotope effects in a realistic environment have been made in recent years, in particular using so-called path integral calculations (9, 10) as well as with semiclassical methods combined with variational transition state theory including tunneling corrections (11).
In this work, we employ a version of the path integral simulation technique (10,12,13) together with the EVB method, free energy perturbation, and molecular dynamics (MD) simulations (7,14,15) to investigate the H/D/T KIEs in the first enzymatic step of GlxI. We also examine the influence of quantum mechanical effects (zero-point energies and tunneling) on the rate acceleration relative to the corresponding uncatalyzed process in water.

MATERIALS AND METHODS
The EVB method (14,15) was used in earlier work to describe the potential energy surface for proton transfer from the C1 carbon atom of the hemithioacetal substrate to the catalytic base Glu 172 (7). The reaction in this case is represented by the two valence bond resonance structures (⌽ 1 and ⌽ 2 ) shown in Fig. 2, whose energies are given by the corresponding diagonal elements of the EVB Hamiltonian.
Here the first three terms denote intramolecular potential energies (Morse bond potentials, angles, dihedral, and improper torsion potentials), and the fourth term describes the nonbonded electrostatic and van der Waals interactions between the reacting fragments. V rs (i) represents the interactions between the reacting groups and their surroundings (enzyme and water), whereas V ss is the potential energy of the surrounding protein-solvent system. The last term is the intrinsic gas phase energy of valence bond state i with the fragments at infinite separation. The gas phase free energy difference ⌬␣ ϭ ␣ (2) Ϫ ␣ (1) as well as the off-diagonal Hamiltonian elements H 12 describing the mixing of the pure valence bond states were calibrated as described in Ref. 7 to reproduce the experimental reaction free energy and activation barrier of the corresponding uncatalyzed reaction in solution. The reaction free energy was calculated from the difference in pK a between Glu (4.1) and the substrate (13.5) (7), and the activation free energy was obtained from the linear free energy relationship in Fig. 3 of Ref. 16. This yielded reaction and activation free energies of ⌬G wat 0 ϭ 12.8 kcal/mol and ⌬G ‡ wat ϭ 22.1 kcal/mol, respectively. The effects of quantum mechanical nuclear motion were examined by treating the substrate C1 atom, the transferring hydrogen, and the acceptor oxygen as quantized particles employing the centroid path integral technique (10,12). In this formulation, each quantized particle is represented by a closed ring (or necklace) of quasiparticles, which are sequentially connected by harmonic springs and each experiencing a fraction of the external potential acting on the real particle. The effective quantum mechanical potential of such a quantized particle is given by, where P is the number of quasiparticles (with coordinates x j ) in the necklace, m is the mass of the real atom, and V cl is the classical potential on the atom. For the interaction between the quantized atoms, quasiparticle j in one necklace interacts with the corresponding j:th bead in the other. As one can see from Equation 2, the dispersion of the quantized atom depends on its mass and the temperature such that it approaches the classical limit as these quantities increase.
Although path integral simulations can be carried out by direct MD on all (quantum and classical) particles, such calculations will suffer from convergence problems, at least as far as small energy differences between different isotopes are concerned, because of the regular noise and divergence of trajectories. However, the fact that the motion of the center of mass (centroid, denoted by x ) of the quantized particle can be rigorously separated from the fluctuations of the quasiparticles around the centroid allows much more efficient computations (10). The two types of trajectories can thus be generated independently. Here we use a standard Monte Carlo (MC) procedure to generate configurations of the different (free particle) quantized necklaces, which are then used together with MD simulations of the full classical system including the corresponding centroids (which are then treated as regular classical particles). By this computational scheme, the usual perturbation formula for calculating reaction free energy profiles (7,14,15) now takes the form, where X is the reaction coordinate, ␤ ϭ 1/kT, V m is the free energy perturbation mapping potential determined by the mapping parameter , and E g is the EVB ground-state energy (see Ref. 15 for details), and the sum runs over all quasiparticles in each necklace. The two types of averages involved are thus taken over the free particle necklace distributions ( fp ) and over the classical MD trajectories on the mapping potential ( Vm ). All simulations of the enzyme-catalyzed reaction used the crystal structure of GlxI in complex with S-(N-hydroxy-N-p-iodophenylcarbamoyl)glutathione (6) as the starting point, where the ligand coordinates were used to dock the S-hemithioacetal of glutathione and methylglyoxal. The active site Zn 2ϩ ion was replaced by Mg 2ϩ , because the two ions give essentially the same catalytic effect (7) and the latter is simpler to treat computationally. For the uncatalyzed reaction, the starting coordinates of the reacting groups were directly taken from the above crystal structure. The reaction center was surrounded by an 18 Å sphere of SPC water (17)  same size, containing both protein and water, in the enzyme simulations. The protein atoms outside this sphere were restrained to their x-ray coordinates and interacted only via bonds across the boundary. A nonbonded cutoff of 10 Å was used together with multipole expansion for long range electrostatics (18). The water surface was subjected to radial and polarization restraints as described elsewhere (19,20). Net charges were assigned to Arg 37 , Glu 99 , Lys 150 , and Glu 172 because these amino acid residues are close to the reaction center, whereas distant ionizable groups were treated as neutral. The simulations were initiated from structures which had been equilibrated for 50 ps.
All perturbations were carried out using 51 discrete free energy perturbation steps (s). An MD timestep of 1 fs and a constant temperature of 300 K was used in all calculations. Equilibrium bond lengths and dissociation energies for the reacting groups were taken from standard Morse potentials and AM1 calculations (7). The GROMOS87 force field was otherwise used in the calculations (17) together with the ion parameters in Refs. 21 and 22. The EVB Hamiltonian (⌬␣ and H ij ) was parameterized against available data on the activation barrier and reaction free energy for the uncatalyzed reference reaction in solution as described in Ref. 7. In the path integral simulations each quantized particle was represented by 20 quasiparticles and 6 ϫ 10 6 MC configurations (after equilibration) were generated for the quantized 1 H, 2 H, 3 H, 16 O, and 12 C necklaces at 300 K. All MD simulations and free energy calculations were carried out using the program Q (20).
Three sets of simulations with different ratios between classical and quantized particle sampling were performed for the reaction with each hydrogen isotope in the enzyme and in water. The first set of simulations had a classical trajectory length (excluding equilibration) of 255 ps with 5 MC configurations/classical step. The second also used 255 ps MD trajectories, but with only 1 MC configuration/step. The third set used 127.5 ps of classical MD with 10 MC configurations/step. Energy data were collected every 5 fs in all cases. Because of the separate calculations of classical trajectories and the necklace configurations, each such trajectory could be combined with different isotope quasiparticle distributions so that isotope effects were evaluated on the same classical trajectory. Likewise, the same MD trajectory could be used with and without considering quantized atoms to examine the magnitude of the quantum effects. As an additional check the enzyme calculations were repeated with an extended set of six quantized atoms (forming angles involving the transferring hydrogen), which yielded virtually identical results.

RESULTS AND DISCUSSION
One of the main advantages with the EVB method is the possibility to calibrate the reaction free energy surface against experimental data for relevant reference reactions in solution. Because observed rate and equilibrium constants for a given solution reaction include the effects of quantum mechanical zero-point energy and tunneling as well as barrier recrossings, these effects are implicitly taken into account by the standard EVB approach. When calculating enzyme reaction free energy profiles on the basis of such aqueous solution calibrations, it is assumed that the magnitude of these quantum and dynamical effects is similar in the catalyzed and the uncatalyzed reactions. That is, if these "corrections" to the classical transition state theory are similar in the two environments then it is reasonable to expect that classical MD calculations for the enzyme and solution reactions will capture most of the catalytic effect.
That the above assumption is, in fact, a very good approximation is immediately evident from Fig. 3, which shows the results from classical and path integral simulations of the GlxI reaction in the enzyme and in solution. The free energy profiles denoted by diamonds in Fig. 3 were obtained from classical trajectories where the solution profile has been calculated as described in Ref. 7 against experimental data. If instead the quantized path integral trajectories are used, again with parameterization of the water profile against experiment, we obtain the free energy curves denoted by circles in Fig. 3. The fact that the two different enzyme profiles, obtained with and without explicit simulation of quantum mechanical nuclear motion, essentially coincide demonstrates that the quantized motions contribute very little to the catalytic effect. The actual difference in activation free energy between the classical and quantum treatments for this reaction step is only 0.2 Ϯ 0.3 kcal/mol. The free energy profiles denoted by triangles in Fig. 3 are evaluated from the classical trajectories but using the calibration from the quantized path integral simulations of the solution reaction. The deviations of these classical profiles from the quantum ones (illustrated by vertical arrows in Fig. 3) therefore directly depict the magnitudes of quantum effects in the two microenvironments. Hence, it is clear that zero-point energy and tunneling effects contribute significantly to the rate of proton transfer, (⌬⌬G ‡ qu ϳ2.5 kcal/mol), but because the contributions are almost identical in the enzyme and water reactions the overall contribution to the catalytic effect is negligible. The same conclusion was reached in a recent study of proton transfer of carbonic anhydrase, although an error range of 2 kcal/mol was reported there (10). Our present error bar for the difference in ⌬⌬G ‡ qu between enzyme and water is only 0.2 kcal/mol, so we consider this conclusion strongly supported by the present data. Fig. 4 shows a snapshot from the simulations

FIG. 3. Free energy profiles for the proton abstraction reaction (A) in GlxI and (B) in water.
The reaction coordinate ⌬⑀ 12 is the energy gap between the valence bond states ⌽ 1 and ⌽ 2 (15). The curves plotted with diamonds are classical simulations where the profile in water was calibrated as in Ref. 7. The curves denoted by circles are quantized path integral simulations, also with the water profile calibrated against experimental data. The curves plotted with triangles are free energy profiles from the classical simulations obtained with the calibration from the quantized simulations in water (see "Results"). with 1 H in the transition state region, where the arrangement of the catalytically important residues can be seen. In addition, the essence of the path integral treatment can be appreciated from the spatial dispersion of the quantized particles.
The fact that the present EVB model, both without explicit treatment of quantum effects and combined with path integral simulations, yields an activation barrier for the first reaction step in GlxI that agrees well with kinetic data (23), makes it particularly interesting to examine the performance of the latter approach for predicting KIEs. Fig. 5 shows the results of path integral simulations of H/D/T abstraction from the hemithioacetal substrate by Glu 172 in GlxI and by glutamate in aqueous solution. The free energy profiles in Fig. 5 are given relative to the same free energy perturbation mapping (or reference) potential so that the combined quantum mechanical effects because of zero-point energy and tunneling along the reaction path can be appreciated. The calculated average KIEs in the enzyme and in water from three different simulations are given in Table I. We obtain here a H/D KIE of 5.0 Ϯ 1.3 for the enzyme reaction which is in reasonable agreement with the experimental value of about 3, considering the small free energy differences involved (see below). The corresponding calculated H/D KIE for the uncatalyzed reaction in water is 3.6 Ϯ 0.7. The calculated H/T and D/T KIEs are also given in Table I but no experimental values for these have been reported so far. Of course, because we are dealing with very small activation free energy differences between the isotopes, i.e. a KIE of 3 corresponds to a ⌬⌬G ‡ of only 0.6 kcal/mol at room temperature, it is essential to gauge the accuracy of these results. Thus, we carried out several independent simulations with different classical centroid trajectories and different MC samples of the necklace quasiparticle distributions. Analysis of these MD/MC runs shows that the present convergence error range for the KIE is only about 0.2 kcal/mol. To reach this level of accuracy, however, it seems necessary to use fairly long MD trajectories, which together with the MC treatment of the quasiparticles make the calculations rather time consuming but not intractable. Even with relatively long MD trajectories on the order of several hundred ps it is well known that the absolute activation free energies are associated with larger errors, typically ϳ1 kcal/mol. Here, it is the separability of the classical (centroid) and quantized (necklace) motions (10) that allow a significantly higher accuracy to be attained for the relative activation free energies between different isotopes. That is, each KIE calcula-tion can be carried out using the same centroid trajectory for H, D, and T, and it appears that the KIEs do not depend strongly on the classical trajectories. On the other hand, as noted above, different classical trajectories yield slightly different absolute free energy profiles with convergence errors in ⌬G ‡ of about 1 kcal/mol.
For comparison with the above calculations of KIEs in GlxI and in water, it is of interest to examine the intrinsic KIEs originating from zero-point energies of the corresponding gas phase reaction. We thus calculated the ground and transition state structures for the elementary proton transfer step of a model system where Glu 172 was replaced by acetate and the substrate truncated so that the -S-CH 2 -linker was replaced by S-CH 3 (cf. Fig. 2) using AM1 and ab initio 6 -31ϩG* calcula-  tions (24). With a standard scaling of Hartree-Fock frequencies by 0.89 (to compensate for the neglect of electron correlation) and only considering the zero-point energy contributions we obtained KIEs of 4.9, 9.5, and 2.0 for H/D, H/T, and D/T, respectively (the calculated bare energy barrier is 27.8 kcal/ mol). As a comparison, the unscaled AM1 frequencies give KIEs of 5.5, 11.6, and 2.1 ( Table I). The calculated values of the isotope effects are thus rather similar in the three environments (gas phase, water, and enzyme), although there is some indication of larger H/D and H/T effects in the enzyme and gas phase compared with solution. The geometrical arrangement of the donor hydrogen acceptor atoms of the EVB and quantum mechanical gas phase transition states are also quite similar. The rate of hydrogen tunneling in enzyme reactions has attracted considerable attention in recent years, largely because of the experimental efforts of Klinman and co-workers (25,26). In the semiclassical picture, where KIEs only originate from the differences in zero-point energy, the ratio between the H/D and D/T isotope effects are simply predicted from the masses of the isotopes and one obtains ln(H/T)/ln(D/T) ϭ 3.3 (27). Measurements of the ratios of these isotope effects in enzymes have therefore been used as a probe of hydrogen tunneling (25,26). In particular, if the ratio ln(H/T)/ln(D/T) is significantly larger than 3.3 this would indicate that transfer of the lighter proton has a substantial contribution from tunneling. Another useful probe of tunneling effects is the temperature dependence of KIEs (25, 26,28). As far as the ratio ln(H/ T)/ln(D/T) is concerned, we obtain here values of 4.8 Ϯ 1.0, 3.2 Ϯ 0.3, and 3.4 for the enzyme, solution, and ab initio gas phase reactions, respectively. Whereas the calculated KIE ratio in the enzyme appears to differ somewhat from the solution and gas phase reactions, where the ratio is close to that predicted semiclassically, the present accuracy is unfortunately not enough to judge the significance of this difference. On the other hand, our comparison above (Fig. 3) of the reaction free energy profiles with and without simulation of quantum effects clearly shows that the overall magnitude of these effects is very similar in the enzyme and solution reactions. A better test case for examining the issue of nuclear tunneling contributions would therefore probably be the alcohol dehydrogenases studied by Klinman and co-workers (25, 26), where tunneling is said to be significant, and the temperature dependence is unexpected.
In summary, we have used a combination of EVB and path integral simulations in this work to examine kinetic H/D/T isotope effects on the initial proton abstraction step in GlxI. This reaction step is believed to be rate-limiting in the enzyme and is associated with an observed H/D KIE of about 3. Our calculated reaction free energy profiles both from classical and path integral simulations agree well with the experimental activation barrier for the proton abstraction step. Furthermore, we find that the contribution of quantum mechanical effects to catalysis are very small, i.e. the magnitude of these effects are essentially the same in the enzyme and in the uncatalyzed solution reaction.
Nevertheless, it is important to treat quantum and dynamical effects in a consistent way in computational modeling studies of enzyme reactions. In the standard EVB approach (14,15) this problem is solved by calibration of the potential surface against experimental data, which implicitly include such effects. Most other so-called quantum mechanics/molecular mechanics methods are not parametrized against rate constants so that zero-point energy, nuclear tunneling, and barrier recrossing corrections, in principle, need to be considered explicitly. However, it is important to realize that zero-point energy corrections are dominating and that most semiempiri-cal methods, which are often used in quantum mechanics/ molecular mechanics calculations, are parametrized against experimental heats of formation and also implicitly include the zero-point energies. To get direct information regarding the magnitude of the above mentioned corrections to classical transition state theory, it appears that both path integral simulations of the type described herein and the approach of Alhambra et al. (11) can be very useful. Regarding dynamical effects in enzyme reactions, which can be expressed in terms of a transmission factor in the transition state theory, reactive flux simulations of the type recently reported by Neria and Karplus (29) can provide useful information. However, so far there is no convincing evidence that dynamical effects are important in enzyme catalysis (14,29,31). This should not be interpreted in such a way that the enzyme flexibility and conformational changes are unimportant for their function. But there is no evidence that the actual dynamics, or time dependence of motions, contribute significantly to the catalytic rate enhancement compared with uncatalyzed solution reactions. In other words, it is the regular Boltzmann probability of being at the transition state (that could be calculated by time-independent Monte Carlo simulation as well) that is by far the most important determinant of the rate acceleration and not the detailed nature of the trajectories leading to the transition state.
Our calculated value of the H/D KIE for proton abstraction from the hemithioacetal substrate in GlxI is about 5.0 Ϯ 1.3. Within the accuracy of the simulations this value is compatible with the experimentally measured effect of around 3 (23). The corresponding KIE calculated for the uncatalyzed solution reaction is 3.6 Ϯ 0.7, whereas the semiclassical KIE (excluding tunneling) for the gas phase reaction is found to be 4.9 from ab initio calculations. The fact that the KIEs obtained for the different environments and with different methods are rather similar would suggest that this range of H/D KIE is intrinsic to the reaction and not largely affected by the surroundings. In the enzyme reaction it is expected that once the H/D/T has been initially abstracted by the general base it will rapidly be exchanged with solution (30). Therefore, it indeed seems reasonable to interpret the observed KIE in terms of the initial proton transfer being rate-limiting. It can be noted that several of the other enzymes mentioned above that catalyze reactions of carbon acids also show H/D KIEs of around 3 (e.g. triosephosphate isomerase (3), yeast enolase (4), and mandelate racemase (5)). This may indicate that, despite the differing substrates and details of the catalytic machinery, the transition states for enzyme-catalyzed proton abstraction from carbon acids are rather similar.