Coupling of hydrogenic tunneling to activesite motion in the hydrogen radical transfer catalyzed by a coenzyme B_{12}dependent mutase
 *Institute of Applied Radiation Chemistry, Department of Chemistry, Technical University of Lodz, Zeromskiego 116, 90924 Lodz, Poland;
 ^{§}Department of Biochemistry, University of Nebraska, Lincoln, NE 685880664; and
 ^{¶}Department of Chemistry and Supercomputing Institute, University of Minnesota, 207 Pleasant Street SE, Minneapolis, MN 554550431
See allHide authors and affiliations

Edited by Richard L. Schowen, University of Kansas, Lawrence, KS, and accepted by the Editorial Board May 14, 2007 (received for review March 9, 2007)
Abstract
Hydrogen transfer reactions catalyzed by coenzyme B_{12}dependent methylmalonylCoA mutase have very large kinetic isotope effects, indicating that they proceed by a highly quantal tunneling mechanism. We explain the kinetic isotope effect by using a combined quantum mechanical/molecular mechanical potential and semiclassical quantum dynamics calculations. Multidimensional tunneling increases the magnitude of the calculated intrinsic hydrogen kinetic isotope effect by a factor of 3.6 from 14 to 51, in excellent agreement with experimental results. These calculations confirm that tunneling contributions can be large enough to explain even a kinetic isotope effect >50, not because the barrier is unusually thin but because cornercutting tunneling decreases the distance over which the system tunnels without a comparable increase in either the effective potential barrier or the effective mass for tunneling.
The determination of the crystal structure of coenzyme B_{12} or 5′deoxyadenosylcobalamin (AdoCbl) by Lenhart and Hodgkin (1) revealed, unexpectedly at the time, a cobaltcarbon bond to C5′ of the deoxyadenosyl ligand. Homolytic rupture of this bond is a key step in the carbonskeletal rearrangements catalyzed by AdoCbldependent enzymes (2, 3), although atomistic details of the mechanism have remained elusive. B_{12}dependent rearrangements play vital roles in amino acid catabolism and fermentation, and methylmalonylCoA mutase (MMCM) is the only B_{12}dependent isomerase that is found in both mammals and bacteria. MMCM catalyzes the rearrangement of methylmalonylCoA to succinylCoA, which then can enter the Krebs cycle. An anomalously large kinetic isotope effect (KIE) has been observed for the hydrogen atom transfer from substrate to the deoxyadenosyl radical (dAdo·) in the step that initiates the isomerization reaction in the enzyme (4). This KIE poses a challenge to theoretical understanding (5), and we use it here as a key to elucidating atomistic details of the reaction catalyzed by MMCM.
Recent advances in computational enzyme kinetics have allowed for a detailed understanding of how many enzymes work, and it now is appreciated that enzymecatalyzed proton and hydride transfer reactions often occur by quantum mechanical tunneling (6, 7). Hydrogen atom transfers are more elusive because they often proceed by protoncoupled electron transfer. However, B_{12}dependent isomerases and reductases are believed to operate by radical translocation (3, 8–11), and various substrate and product radical intermediates have been detected by EPR spectroscopy in these enzymes. The observation of very large hydrogen/deuterium (H/D) KIEs for hydrogen radical transfer to or from dAdo· in several enzymes (5, 12–15) raises the possibility that tunneling may be a common catalytic strategy for such translocation. Confirmation of this interpretation presently is best achieved by atomistic simulation of the KIEs because, despite considerable experimental attention, some key mechanistic details remain unclear (4, 16).
The rate of CoC bond homolysis is accelerated by a factor of ∼10^{11} in the ternary complex of enzyme, AdoCbl, and substrate (S) relative to the uncatalyzed rate in solution (17). We consider the mechanism where reaction 1 is homolysis, reaction −1 is recombination, and reaction 2 is the transfer of H or D from the methyl group of the substrate to the 5′carbon of dAdo·, as illustrated in Scheme 1. Another question regarding the mechanism of B_{12}dependent mutases that brings some confusion is the recently proposed concertedness of the homolysis and the Hatom transfer steps that, until now, usually were considered separate steps. In support of the hypothesis (18) that the productlike radical is stabilized by cob(II)alamin, Kozlowski et al. (19) have found a transition state that comprises homolysis and Hatom transfer. These density functional theory results, in fact, are very similar to those reported earlier on the basis of semiempirical calculations (5). These results, however, are in disagreement with recent calculations that demonstrate the stepwise mechanism with forward and reverse potentialenergy barriers of only ∼10 and ∼8 kcal/mol for reactions 1 and −1, respectively, which makes both reactions much faster than reaction 2, even though step 1 positions the dAdo· radical to facilitate step 2 (20). Therefore, the measured KIE should correspond to step 2. The most probable reason for the differing conclusions seems to be the extreme sensitivity of the system to the size of the model and the details of the dynamical model (5, 16, 21). Our experience indicates that only inclusion of a substantial part of the enzyme environment yields results that are pertinent to the enzymecatalyzed reaction (22).
In the present work, we calculate the quantum tunneling dynamics of reaction 2 and analyze the molecular mechanism of the tunneling process in the enzymecatalyzed reaction, which allows us to understand the tunneling dynamics at a level of detail not previously possible. We previously have shown (23) that Austin Model 1 (AM1) (24) provides accurate (±2 kcal/mol) barrier heights for H radical transfer between carbon centers, and the present calculations employ a combined quantum mechanical and molecular mechanical (QM/MM) (25, 26) potentialenergy surface with unrestricted (27) AM1 for a 45atom QM part and CHARMM22 (28) for the MM part (14,833 atoms).
Rate constants measured under ordinary conditions correspond to an average over fluctuating conformational states of the protein (29). The present calculations include this average by ensembleaveraged variational transitionstate theory (EAVTST). The simulations include 672 aa of Protein Data Bank entry 4REQ (8) plus 299 atoms of the substrate and coenzyme and 1,388 water molecules.
Results
By using methods described previously (30–32), tunneling was included by an ensemble average over eight reaction valleys descending from points selected from the transitionstate ensemble. Including tunneling increases the rate constant for the unsubstituted (perprotio) case by a multiplicative transmission coefficient (33) κ_{H} and for the perdeuterated methyl case by κ_{D}. The KIE increases by κ_{H}/κ_{D}. Tunneling is included by several models corresponding to successively improved descriptions of the multidimensional tunneling process. Zerocurvature tunneling (ZCT) is a model (34) in which tunneling, for each reaction valley, occurs along the minimum energy path (MEP) through that valley, as would be the case if the MEPs were not curved. Each actual MEP curves its way through the 135dimensional primary zone. Tunneling paths on the concave sides of the MEPs increase the tunneling probability because cornercutting shortens the path and decreases the distance over which exponential decay of the nuclearmotion wave function occurs (33, 35–41). Smallcurvature tunneling (SCT) (42) allows cornercutting (illustrated in refs. 33 and 41) but only out to approximately the concaveside turning points of the vibrations transverse to the MEPs. The largecurvature tunneling (LCT) approximation (43) allows tunneling on paths both near and far from the MEP, even extreme cornercutting. Finally, we used microcanonically optimized multidimensional tunneling (OMT) (43), which variationally selects from among both SCT and LCT paths and further increases the calculated transmission coefficients.
The EAVTST calculations, with quantized vibrations but classical motion along the reaction coordinate, give a KIE at 278 K of 14.3, much smaller than the experimental value (4) of 49.9. These values are compared with the results including tunneling in Table 1. Table 1 first shows that ZCT gives only a modest increase in the KIE, to 22.0, but that even modest cornercutting, as in the SCT approximation, approximately doubles the KIE as compared with tunneling along the MEPs. Allowing extreme cornercutting, as in either the LCT or OMT approximations, yields a KIE that agrees with experiment within experimental error.
Discussion
Because the final results of the quantum mechanical atomistic simulation agree with experiment so well, we can analyze them to better understand the nature of the tunneling events. In a fully quantal calculation, the wave packet is delocalized, and interpretation in terms of well defined tunneling paths is not possible. A key issue here is that the tunneling calculations are semiclassical, and they correspond to calculating imaginary action integrals iθ along tunneling paths that can be visualized (35) where the tunneling probability is approximately (1 + e ^{2θ})^{−1}. When the system tunnels at an energy significantly below the effective barrier top, θ is large, and the tunneling probability reduces approximately to e ^{−2θ}. Semiclassically, the magnitude ψ of the wave function decreases as e ^{−θ}, and e ^{−2θ} may be associated with the decay of the probability, ψ^{2}, of finding the tunneling particle on the other side of the barrier. The imaginary action integrals iθ involve three physical quantities: the effective barrier for tunneling, the effective reduced mass for tunneling, and the distance along the tunneling path. A high effective barrier, a high reduced mass, or a long path each decreases the tunneling probability. A key element of the multidimensional tunneling approximations is that they allow cornercutting in the configurationspace path from reactants to products; cornercutting shortens the path and thereby increases the tunneling probability, as illustrated in Fig. 1.
One of the key goals of chemical kinetics is finding the best transition state, which is the hypersurface of maximum free energy of activation (44, 45). The coordinate normal to this hypersurface is the variationally best reaction coordinate. Therefore, we can find the hypersurface of maximum free energy of activation by finding the value of the reaction coordinate that corresponds to the maximum free energy of activation. Thus, instead of taking the reaction coordinate to be the coordinate of the single hydrogen atom that is transferred, the goal is to use this variational principle to identify the participation of individual atomic coordinates in a multidimensional reaction coordinate. When a reaction is dominated by tunneling, this goal becomes to identify the participation of individual atomic coordinates in the tunneling paths. Analysis of the transmission coefficients in Table 1 accomplishes this objective and provides considerable insight into the atomistic details of the tunneling process.
In ZCT, the coupling of the hydrogenic motion to its environment affects the energy available for tunneling but not the path of the tunneling atom. Table 1 shows that the ZCT model qualitatively underestimates the effect of tunneling. Thus, the model of tunneling along the MEP is inadequate, and the participation of the other degrees of freedom (beyond the progress coordinate of the hydrogen atom) in the reaction coordinate for tunneling is essential to explaining the large KIE. We conclude that the large KIE does not result mainly from a particularly thin barrier but rather because cornercutting tunneling decreases the distance over which the system tunnels without a comparable increase in either the effective potential barrier or the effective mass for tunneling.
Why does the multidimensional character of the tunneling raise the KIE? The multidimensional nature of tunneling may be considered to have two components. First is the effect of mass on orthogonal vibrations (i.e., vibrations orthogonal to the tunneling path). Second is the effect of mass on the tunneling path itself.
For a typical reaction, the zero point vibrational energy is smaller at the transition state than for reactants, which lowers the requirement for energy in modes transverse to the reaction path, and thereby it lowers the effective barrier to reaction. Because the lowering tends to be larger where the barrier is higher, this kind of barrierlowering tends to flatten and broaden the effective barrier to tunneling. Furthermore, because zero point vibrational energy tends to be larger for lower masses, the amount of zero point vibrational energylowering, and hence the amount of flattening, tends to be more significant for the lower mass, and thus it tends to lower the amount of tunneling more for the lighter mass. Because the lighter mass appears in the numerator of the KIE, this effect tends to reduce the KIE from a onedimensional picture. Note that even though multidimensional effects lower the KIE in this scenario (so far), the KIE still usually is greater than one because the onedimensional tunneling picture greatly favors H tunneling as the square root of the effective mass for tunneling appears in a negative exponent in the tunneling probability.
Next we consider the fact that we have different paths for H and D transfer, partly because of cornercutting. Although tunneling actually occurs along a distribution of paths, we will discuss the relevant issues in terms of a single representative or dominant tunneling path for each mass case. The tunneling path with the greatest tunneling probability will dominate; hence, the system will cut the corner only to the extent that this cornercutting increases the probability of tunneling. In general, H will take advantage of this more than D, and hence this multidimensional tunneling effect usually will increase the KIE. Why is that?
We can give two reasons. First, as stated above, the H tunneling probability usually is larger, which amplifies any effect that tends to increase tunneling, so the increase is amplified for H. Second, the MEP for H tunneling usually is associated with a smaller skew angle in isoinertial coordinates. (Isoinertial coordinates are any massweighted or massscaled coordinates in which the reduced mass is the same in all directions. The skew angle is the angle between the entrance coordinate and the exit coordinate of the MEP, where the entrance coordinate runs along the entrance valley leading to the saddle region, and the exit coordinate runs along the exit valley leading away from it.) A smaller skew angle means that the shortening of the tunneling distance by a given amount of cornercutting in bond distance space (i.e., by a given amount of potentialenergy barrierraising) is greater for H than for D. Thus, H tends to cut the corner more, even though this leads to a higher barrier (recall that the tunneling path with the greatest tunneling probability will dominate, and hence the system will cut the corner only to the extent that this increases the probability of tunneling; in isoinertial coordinates, this involves the competition between barrierraising and pathshortening, whereas in general coordinates, one would have a third factor, namely, finding a path with a small reduced mass).
Our multidimensional tunneling methods include these competitions as part of the calculation, and we find, in both enzymes and smallmolecule reactions, that cornercutting tends to increase the KIE. A key factor that may appear counterintuitive from some points of view is that although H sees a higher energy along its tunneling path, the tunneling probability is greater because in isoinertial coordinates the path is shorter. If that were not the case, H would not cut the corner more because the tunneling path with the greatest tunneling probability will dominate, and in isoinertial coordinates, this involves the competition between barrierraising and pathshortening. This argument explains why the distribution of H tunneling paths leads to more tunneling than the D distribution does.
It is very difficult to characterize inhomogeneous dynamics in an ensembleaveraged measurement (46), but EAVTST provides information on the distribution of transmission coefficients (33). The standard deviations of the transmission coefficients are given in Table 1 and are small compared with the effect of reactionpath curvature, indicating that the participation of fluctuating coordinates in the tunneling path is dominated by the coordinates in the 45atom primary zone, not the fluctuations of the rest of the substrate, coenzyme, and protein. Because Table 1 shows that the LCT model is in qualitative agreement with the optimized tunneling calculations, we can make further progress in understanding the tunneling paths by examining a representative LCT path for a representative secondary zone of the transitionstate ensemble. By statistical analysis, we found a representative configuration of the secondary zone, and for this configuration, we found the most likely tunneling energy both for H and for D; we then analyzed representative LCT tunneling paths for these energies and secondary zone configurations. The results are presented in Table 2, where R _{aH} denotes the distance from H or D to the acceptor carbon, R _{dH} denotes the distance from H or D to the donor carbon, R _{ad} denotes the acceptor–donor distance, and V − V ^{‡} denotes the potential energy relative to that at the nearest saddle point. The critical configuration along a tunneling path is defined for Table 1 as the point of maximum potential energy along the path. The termini are points on the MEP at the start and end of the tunneling path. Table 2 shows that the tunneling paths start and end at potential energies ≈5–6.5 kcal/mol below the saddle point energy and pass within a few hundredths of an angstrom of the saddle point, at least as far as the distances tabulated are concerned.
Fig. 2 a (for H) and b (for D) show the isotopedependent vibrationally adiabatic groundstate potentials for the representative protein configuration. This is the effective potential for both ZCT and SCT tunneling and for the regions near the path termini for LCT and OMT tunneling. Fig. 2 a and b also show, as horizontal dotted lines, the representative tunneling energies in the OMT approximation. For H, the most typical tunneling energy is 3.4 kcal/mol below the barrier in the effective potential, whereas for D, it is 3.8 kcal/mol below the corresponding barrier top. We emphasize, as a very general point, that in a onedimensional tunneling model the most likely tunneling energy for D would be closer to the top of the barrier than that for H; so this is a clear indication of the inadequacy of onedimensional tunneling approximations.
Figs. 3 and 4 show twodimensional representations of the reaction paths for H and D transfer, respectively. Saddle points are represented by squares, critical configurations are shown by dots, and representative tunneling paths are denoted by dashed lines.
This simulation illustrates the importance of multidimensional quantum dynamics in enzymecatalyzed radical kinetics and confirms that the simulation tools now available for modeling such processes explain the atomistic details. We did not require the explicit inclusion of protein dynamics in individual reaction paths, although coupled protein motion effectively is included in the reaction coordinate by our ensemble averaging procedure, and we did not require abandoning the framework of transitionstate theory, as sometimes has been suggested. Instead, we generalized the idea of a dynamical bottleneck (transition state) for overbarrier reaction also to encompass the critical configurations that dominate the tunneling although a multidimensional barrier, and we allowed the dominant multidimensional reaction paths for tunneling to differ from those that dominate overbarrier kinetics. Representative tunneling paths for both H and D transfer reveal that the tunneling event, which increases the rate by two orders of magnitude, mostly is attributable to cornercutting. These tunneling paths illustrate the nature of tunneling dynamics at an unprecedented level of detail; in particular, we have shown that the large KIE in the hydrogen atom transfer between methylmalonylCoA and dAdo· results from tunneling along paths close to, but not along, the minimumenergy path for reaction. Although the tunneling paths are close to the minimumenergy path, the deviation (along the optimum tunneling paths) of primaryzone atoms from the minimumenergy path greatly enhances the tunneling probabilities, especially for H, because of the exponential dependence of the tunneling probabilities on the magnitude of the imaginary action integral. The KIE calculated from the optimized tunneling calculations, in conjunction with variational transitionstate theory, is in excellent agreement with experiment.
Methods
The simulation system contains 14,878 atoms. This number results from 672 amino acid residues, 1,388 water molecules (in a 30Å equilibrated water sphere), and 299 atoms of the substrate and coenzyme (see Fig. 5). The QM subsystem consisted of 45 atoms, including the entire Ado ligand and part of the substrate (methylmalonylCH_{2}−). The boundary between the QM and MM subsystems was defined by using the generalized hybrid orbital (GHO) approach (25, 26) and placed on the C_{2} atom of the βmercaptoethylamine part of the CoA chain. The total charge of the QM part was −1. The radical species (Ado· and succinyl radicals) were described by using the unrestricted Hartree–Fock (UHF) method (27), which was implemented and used for the first time within the semiempirical module of the CHARMM (47) program. It is known that the corrin ring gets flattened upon binding to MMCM. To model the ring as planar without determining a new set of parameters, the corrin ring was described by using atom types for the heme molecule. All of the nitrogen atoms were assigned as pyrrolic nitrogens, and the carbons in the methylene bridges were assigned as heme alpha and mesocarbons. The cobalt and the carbon atoms from the pyrrolic rings were defined as new atom types. All histidine residues were protonated, with the exception of His610, which is the lower ligand of the cobalt atom. The total charge of the system was −7. Additional CHARMM parameters required for AdoCbl and substrate were based on crystallographic (48) geometries and PM3 (49) partial charges.
Dynamics calculations were carried out with the CHARMM program (47) and the CHARMMRATE (50) module, with quantum effects (vibrational quantization and participation in tunneling) included for a 45atom primary zone by EAVTST (30–32) with OMT (43). The rest of the system was treated by classical mechanics and was called the secondary zone. We used stochastic boundary conditions (51).
For the KIE calculations, we considered a perdeuterated methyl group on the substrate, as in the experimental measurement (4).
Acknowledgments
We thank Antonio Fernandez Ramos, Titus Albu, and Mireia GarciaViloca for helpful assistance. This work was supported in part by State Committee for Scientific Research (KBN), Poland, Grant 3/T09A/70/19 (to P.P.), Fogarty International Research Collaboration Award (FIRCA)–National Institutes of Health Grant TW0121201 (to R.B. and P.P.), National Institutes of Health Grant DK45776 (to R.B.), and National Science Foundation Grant CHE0349122 (to D.G.T.).
Footnotes
 ^{†}To whom correspondence may be addressed. Email: truhlar{at}umn.edu or olaf{at}p.lodz.pl

Author contributions: A.D.D., P.P., R.B., and D.G.T. designed research; A.D.D., P.P., R.B., and D.G.T. performed research; A.D.D., P.P., R.B., and D.G.T. analyzed data; and A.D.D., P.P., R.B., and D.G.T. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission. R.L.S. is a guest editor invited by the Editorial Board.
 Abbreviations:
 AdoCbl,
 5′deoxyadenosylcobalamin;
 MMCM,
 methylmalonylCoA mutase;
 KIE,
 kinetic isotope effect;
 dAdo·,
 deoxyadenosyl radical;
 QM,
 quantum mechanical;
 MM,
 molecular mechanical;
 EAVTST,
 ensembleaveraged variational transitionstate theory;
 MEP,
 minimum energy path;
 SCT,
 smallcurvature tunneling;
 LCT,
 largecurvature tunneling;
 ZCT,
 zerocurvature tunneling;
 OMT,
 optimized multidimensional tunneling.
 © 2007 by The National Academy of Sciences of the USA
References
 ↵

↵
 Corey EJ ,
 Cooper NJ ,
 Green MLH
 ↵
 ↵
 ↵

↵
 GarciaViloca M ,
 Gao J ,
 Karplus M ,
 Truhlar DG

↵
 Kohen A
 Kohen A ,
 Limbach HH
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
 Weisblat DA ,
 Babior BM
 ↵

↵
 Banerjee R ,
 Truhlar DG ,
 DybalaDefratyka A ,
 Paneth P
 Hynes JT ,
 Klinman JP ,
 Limbach HH ,
 Schowen RL
 ↵
 ↵
 ↵
 ↵

↵
 DybalaDefratyka A ,
 Kwiecien RA ,
 Sicinska D ,
 Paneth P
 Brady J ,
 Gao J ,
 Naidoo KJ ,
 Hann M ,
 Field M

↵
 Banerjee R ,
 DybalaDefratyka A ,
 Paneth P
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
 Alhambra C ,
 Corchado J ,
 Sánchez ML ,
 GarciaViloca M ,
 Gao J ,
 Truhlar DG
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
 Kuppermann A ,
 Adams JT ,
 Truhlar DG
 Cobiand BC ,
 Kurepa MV
 ↵
 ↵

↵
 Truhlar DG ,
 Gordon MS
 ↵
 ↵

↵
 Evans MG
 ↵
 ↵
 ↵

↵
 Evans PR ,
 Mancia F
 ↵

↵
 Alhambra C ,
 Corchado JC ,
 Sanchez ML ,
 Villa J ,
 Gao J ,
 Truhlar DG
 ↵