Tunneling and Classical Paths for Proton Transfer in an Enzyme

3032

J. Phys. Chem. B 2007, 111, 3032-3047

Tunneling and Classical Paths for Proton Transfer in an Enzyme Reaction Dominated by Tunneling: Oxidation of Tryptamine by Aromatic Amine Dehydrogenase Laura Masgrau,†,| Kara E. Ranaghan,‡ Nigel S. Scrutton,§ Adrian J. Mulholland,*,‡ and Michael J. Sutcliffe*,† Manchester Interdisciplinary Biocentre, School of Chemical Engineering and Analytical Science, and Faculty of Life Sciences, UniVersity of Manchester, 131 Princess Street, Manchester M1 7DN, United Kingdom, and Centre for Computational Chemistry, School of Chemistry, UniVersity of Bristol, Cantock’s Close, Bristol BS8 1TS, United Kingdom ReceiVed: NoVember 28, 2006; In Final Form: January 12, 2007

Proton tunneling dominates the oxidative deamination of tryptamine catalyzed by the enzyme aromatic amine dehydrogenase. For reaction with the fast substrate tryptamine, a H/D kinetic isotope effect (KIE) of 55 ( 6 has been reportedsone of the largest observed in an enzyme reaction. We present here a computational analysis of this proton-transfer reaction, applying combined quantum mechanics/molecular mechanics (QM/ MM) methods (PM3-SRP//PM3/CHARMM22). In particular, we extend our previous computational study (Masgrau et al. Science 2006, 312, 237) by using improved energy corrections, high-level QM/MM methods, and an ensemble of paths to estimate the tunneling contributions. We have carried out QM/MM molecular dynamics simulations and variational transition state theory calculations with small-curvature tunneling corrections. The results provide detailed insight into the processes involved in the reaction. Transfer to the O2 oxygen of the catalytic base, Asp128β, is found to be the favored reaction both thermodynamically and kinetically, even though O1 is closer in the reactant complex. Comparison of quantum and classical models of proton transfer allows estimation of the contribution of hydrogen tunneling in lowering the barrier to reaction in the enzyme. A reduction of the activation free energy due to tunneling of 3.1 kcal mol-1 is found, which represents a rate enhancement due to tunneling by 2 orders of magnitude. The calculated KIE of 30 is significantly elevated over the semiclassical limit, in agreement with the experimental observations; a semiclassical value of 6 is obtained when tunneling is omitted. A polarization of the C-H bond to be broken is observed due to the close proximity of the catalytic aspartate and the (formally) positively charged imine nitrogen. A comparison is also made with the related quinoprotein methylamine dehydrogenase (MADH)s the much lower KIE of 11 that we obtain for the MADH/methylamine system is found to arise from a more endothermic potential energy surface for the MADH reaction.

Introduction A lively current debate in enzymology centers on quantum tunneling in enzyme-catalyzed proton, hydride, and hydrogen atom transfers.1-12 Experimental studies of kinetic isotope effects (KIEs) have demonstrated the role of quantum mechanical tunneling in enzymic hydrogen-transfer reactions.2,13-20 Atomic-level elucidation of the factors involved in these reactions is key to understanding a large number of reactions in biology. C-H bond cleavage occurs in many biological reactions, and all of these are likely to involve tunneling to some degree. Atomistic simulations have a vital role to play here. Computational modeling methods are able to analyze fundamental mechanisms of enzymic reactions and can provide * Authors to whom correspondence should be addressed. Phone: +44 117 928 9097 (A.J.M.); +44 161 306 5153 (M.J.S.). Fax: +44 117 925 1295 (A.J.M.); +44 161 306 8918 (M.J.S.). E-mail: adrian.mulholland@ bristol.ac.uk; [email protected]. † Manchester Interdisciplinary Biocentre and School of Chemical Engineering and Analytical Science, University of Manchester. ‡ Centre for Computational Chemistry, School of Chemistry, University of Bristol. § Manchester Interdisciplinary Biocentre and Faculty of Life Sciences, University of Manchester. | Current address: Unite ´ de Bioinformatique Structurale, Institut Pasteur, 25-28 Rue du Dr Roux, Paris F-75724, France.

uniquely detailed insight and analysis.7,21-25 For simulations of H transfer, both quantum mechanics/molecular mechanics (QM/ MM)26,27 and empirical valence bond (EVB)21 potentials have been used successfully within the framework of modern transition state theory (i.e., including quantum nuclear effects).7,8,28,29 The quinoprotein enzymessthe mechanisms of which are generally well-understood2,30,31sare ideally suited to studies of proton transfer during substrate oxidation.3 The enzymes aromatic amine dehydrogenase (AADH) and methylamine dehydrogenase (MADH) are bacterial tryptophan tryptophyl quinone (TTQ)-dependent quinoproteins30,32 that catalyze the oxidative conversion of aromatic and aliphatic primary amines, respectively, to the corresponding aldehyde and ammonia The

RCH2NH2 + H2O f RCHO + NH3 + 2H+ + 2erate-limiting step in the reductive half-reaction14 involves transfer of a proton from a C-H bond of the substrate-derived iminoquinone intermediate to the catalytic aspartate base, making these enzymes excellent model systems for investigating tunneling in proton transfer. Our recent experimental and computational studies of AADH have identified quantum tunneling as the dominant mechanism in the oxidative deamination of tryptamine catalyzed by AADH.2,14 The associated

10.1021/jp067898k CCC: $37.00 © 2007 American Chemical Society Published on Web 03/01/2007

Tunneling and Classical Paths for Proton Transfer

J. Phys. Chem. B, Vol. 111, No. 11, 2007 3033

Figure 1. Schematic overview of the reaction mechanism of the reductive half-reaction of AADH with tryptamine.2 Atoms derived from the substrate are depicted in dark gray, and those derived from water in light gray. All other enzyme-derived atoms are depicted in black; the catalytic base, Asp128β, is labeled in I. The proton transfer (H+ tunneling) step is boxed; III and IV refer to the reactant and product of the proton-transfer step, respectively.

H/D KIE of 55 ( 6 (for substitution from perprotio to dideutero, T ) 5-20 °C),2 one of the largest yet observed in any enzyme reaction, is well in excess of the semiclassical limit of ∼7. In the temperature range studied (5-20 °C), an essentially temperature-independent apparent activation free energy of 12.6 kcal mol-1 was obtained for tryptamine oxidation (which at 300 K would correspond to a rate constant for proton transfer of ∼3500 s-1). KIE values elevated over the semiclassical limit have also been reported for the related TTQ-dependent enzyme methylamine dehydrogenase ( which shares 30% sequence identity with AADH).13,33,34 In our previous work,2 high-resolution crystal structures of substrate-free AADH and of AADH in complex with the substrate tryptamine and with the inhibitor phenylhydrazine were presented. Although the structures of several enzyme-bound intermediates were resolved, neither the reactant nor the product of the proton-transfer step could be isolated (III and IV in Figure 1). Therefore, the tryptamine-derived iminoquinone complex (III), the starting point of all our simulations of the reaction by combined QM/MM methods,27,35 was built, following standard procedures with CHARMM, from the 1.1 Å structure of the Schiff base intermediate V.2 The combination of classical molecular dynamics simulations36 with variational transition state theory calculations including multidimensional tunneling corrections (VTST/MT)37-39 permits separate analysis of the classical and quantum reactions. Thus, the contribution of tunneling in lowering the effective free energy barrier, thereby enhancing the rate of reaction, and the effect on the KIEs are obtained. We extend here our previous computational study2,7 on AADH by using improved energy corrections and high-level QM/MM methods. We analyze in detail both the classical and the tunneling paths to obtain a detailed understanding of the atomic changes accompanying this proton-transfer event. The suitability of the reaction path used in the analysis of the tunneling reaction is investigated by comparison to an ensemble of paths obtained by other approaches. We (and others) have also performed variational transition state theory calculations with small-curvature tunneling corrections (VTST/SCT) calculations on the MADH/methylamine system,33,34,40-42 and we compare here the two related systems AADH and MADH.

Methods The proton-transfer step in the AADH/tryptamine system (boxed step in Figure 1) was modeled using QM/MM methods27,35 available within the CHARMM program (versions c27b2 and c29b2),35 with tunneling modeled using VTST/SCT calculations37,38 available via the CHARMMRATE interface.43 QM/ MM methods allow reactions in enzymes to be modeledsthe reacting groups in the active site of an enzyme are treated quantum mechanically, while including the effects of the rest of the enzyme environment with a simpler force field description. QM/ MM simulations can provide structural and energetic insight into enzyme mechanisms not possible with experiment7,21,25,44,45s as exemplified in this study, structural characterization of the reactant or product structures of this key reaction step for X-ray analysis was not possible. Employing umbrella sampling QM/ MM dynamics and VTST/SCT calculations, the classical and tunneling paths for this proton transfer can be followed. The free energy and tunneling calculations reported here were carried out at a simulation temperature of 300 K. Preparation of the Tryptamine-Derived Iminoquinone Reactant Complex (III). The structure of the tryptaminederived iminoquinone complex (III), the starting point of all our simulations, was modeled from the 1.1 Å resolution X-ray crystal structure of the related Schiff base intermediate V (PDB46 code 2AGY).2 Hydrogen atoms were built into the structure using the H-Build subroutine in CHARMM, and the structure was solvated with a sphere of water molecules. The CHARMM22 all-hydrogen force field47 for the protein residues and TIP3type48 water molecules were used for the molecular mechanics representation. The MM parameters for the TTQ cofactor and the bound substrate were taken from our previous work.2 The system was truncated to a 25 Å sphere centered on the NT atom of the tryptamine-derived iminoquinone in the β subunit (Figure 2) and comprised 8182 atoms, 4090 from the protein and substrate and 1364 water molecules. The QM/MM partition was defined with the QM region consisting of 48 atoms (including 3 HQ-type link atoms;49 Figure 2). HQ link atoms are QM hydrogens without classical van der Waals or bonded force field terms but do interact with MM charges. (Charges of MM host

3034 J. Phys. Chem. B, Vol. 111, No. 11, 2007

Figure 2. Schematic of the active site of the reactant, iminoqinone III, illustrating the QM/MM partition used in the calculations. The QM region is in black, and the surrounding MM region in gray, with link atoms circled; hydrogen bonds between the QM and the MM regions are depicted by dashed gray lines. (The backbones of Asp157β and Asn156β and the carboxylate of Asp84β, which form hydrogen bonds with the tryptamine indole, are not shown for clarity.)

groups must be set to zero.) The overall charge of the QM system was zero, from the -1e of Asp128β and +1e of the bound tryptamine. After initial equilibration of the water molecules by molecular dynamics simulation with fixed protein atoms (T ) 300 K), the full model was minimized by imposing gradually decreasing harmonic restraints on the protein atom positions. For QM/MM simulations, a stochastic boundary approach was used,50 and water was restrained to remain within the simulation system using a spherical deformable potential of radius 25 Å.50 The reaction region was defined as a 21 Å sphere centered on NT: Atoms in this region were not subject to positional restraints. The region between 21 and 25 Å was defined as the buffer region: Protein heavy atoms in this region were restrained using force constants scaled to increase with distance from the center of the system. Atoms further than 25 Å from the center were fixed. QM/MM minimization of III was carried out using the adopted basis Newton-Raphson minimization algorithm to a gradient tolerance of 0.01 kcal mol-1 Å-2. Both AM1 and PM3 semiempirical methods were used to describe the QM region. For comparison with the semiempirical methods, III was also optimized at the B3LYP/ 6-31G(d)/CHARMM22 level of theory using Jaguar51 and Tinker.52 Output and input to these programs are linked by inhouse routines described elsewhere.44,53 The results showed (see the Results section) that the B3LYP/6-31G(d)/CHARMM22 structure agreed best with the PM3/CHARMM22 result, and therefore PM3/CHARMM22 was chosen as the more suitable semiempirical/MM method to study this proton-transfer step. QM/MM Adiabatic Mapping. An adiabatic mapping procedure was followed, starting from the PM3/CHARMM22optimized reactant III; reaction coordinates for the two possible proton transfers were defined as ZO2 ) (d(C1-H1) - d(O2H1)) Å and ZO1 ) (d(C1-H1) - d(O1-H1)) Å. This definition of the reaction coordinate has been used since it has been shown previously to model proton transfers well.33,54-56 The system was optimized at different values of ZO2 or ZO1, passing from

Masgrau et al. reactant to product in 0.1 Å steps to obtain an indication of the potential energy change associated with the transfer of the proton to O2 or O1 and of the suitability of the reaction coordinate chosen. Simulation of Classical Reaction Path. The PM3/ CHARMM22-optimized reactant III was equilibrated using a 30 ps PM3/CHARMM22 stochastic boundary molecular dynamics simulation. Langevin dynamics were applied to atoms in the buffer region (the list of atoms in the buffer was updated every 50 steps), using friction coefficients of 250 ps-1 for protein heavy atoms and 60 ps-1 for water oxygen atoms, and a 1 fs time step. Nonbonded interactions, calculated within a 13 Å cutoff, were updated every 25 steps. MM/MM interactions were scaled down using an atom-based shifting function, whereas QM/MM electrostatics were scaled down by a group-based switching function. The equilibrated structure was subjected to a further 30 ps of unrestrained dynamics. Molecular dynamics simulation of the AM1/CHARMM22-minimized reactant III illustrated that the differences in geometry observed between the minimized structures produced using the two semiempirical methods persisted. Therefore, subsequent calculations used only the PM3/CHARMM22-modeled system. Umbrella sampling molecular dynamics (PM3/CHARMM22)36,57 was used to simulate the proton transfer from III to IV in 0.1 Å intervals along the reaction coordinate, from ZO2 ) -1.5 to 0.7 Å and ZO1 ) -0.6 to 1.1 Å, using a force constant of 300 kcal mol-1 Å-2. (This force constant was chosen as it resulted in sufficient overlap of the regions sampled between adjacent Z values.) Each subsequent simulation was started from the 1 ps point of the preceding simulation. The value of Z was recorded at each step of the 30 ps dynamics simulation for statistical analysis, and the reaction coordinate statistics were combined using the weighted histogram analysis method (WHAM58) to give the PM3/ CHARMM22 (classical) free energy profile for proton transfer. The resulting classical activation free energy (∆G q,C) was used to estimate the classical rate constant (the rate constant without inclusion of any tunneling or quantized-vibration contribution). Semiempirical methods such as PM3 perform well for many organic species, but their reliability and accuracy can vary significantly for different systems. It is often necessary to correct semiempirical QM/MM energies to obtain satisfactory reaction energetics.59 PM3 has been shown to represent poorly the energy of proton transfer in the similar methylamine MADH/ methylamine system.34,60 It was therefore necessary to correct the energies to provide more realistic values. PM3 can be reparameterized to reproduce experimental data or high-level ab initio data for a particular reaction of interest. Alhambra and co-workers34 developed specific reaction parameters (SRPs)61,62 for the MADH/methylamine system based on higher-level quantum mechanical results for small models. They found that it was only necessary to modify the parameters for the donor carbon atom involved in the proton transfer to reproduce their higher-level results. As the model reactions used in that work also represent the proton transfer in AADH, the SRPs for this carbon atom in MADH were implemented in our versions of CHARMM, and single-point PM3-SRP/CHARMM22 energy calculations were carried out on the structures along the PM3/ CHARMM22 pathways generated by adiabatic mapping (see above). The difference between the PM3/CHARMM22 and the PM3-SRP/CHARMM22 energies (relative to the reactant) was then used to adjust the free energy at each Z value. This approach was tested by performing a PM3-SRP/CHARMM22 adiabatic mapping procedure, optimizing the structures at each pointsthis procedure gave similar results.



Analysis of Classical Reaction Path. Snapshots of the PM3/ CHARMM22 classical reaction path were recorded every 0.5 ps of the molecular dynamics (MD) simulation and analyzed to determine hydrogen-bonding patterns and the contributions of individual residues. A hydrogen bond was defined as a hydrogen to acceptor distance of e2.6 Å (2.8 Å if the donor is sulfur) and a donor to acceptor distance of e3.5 Å (4.0 Å if the donor is sulfur); the angle between atoms forming the hydrogen bond must also be greater than 90°.63 The contributions of different residues (e.g., to QM/MM electrostatic stabilization) were calculated by individually deleting the residues (identified as important in the hydrogen bond analysis) from each structure in the dynamics trajectory at the reactant (III), transition state (TS), and product (IV) complexes and recalculating the energy. The difference between the original QM/MM electrostatic energy and the energy with the residue deleted was then the contribution of that particular residue to the electrostatic interaction. It is important to note that the interaction energies calculated in this manner do not include the effects of dielectric shielding. As a result the energies are not directly comparable to experimental results (e.g., ∆∆Gq for mutated enzymes); however, this type of decomposition analysis is useful for identifying catalytically important interactions.54,64-66 Mulliken charge analysis was also carried out for three snapshots at each Z value along the path, and the average atomic charges were calculated. Charges of the reactant complex were also calculated for the QM region with no MM region present and also with the QM atoms of Asp128β deleted (MM region present) to investigate the effects of the environment on the charges of the QM atoms. This method is limited by the fact that Mu¨lliken charges do not represent the charge distributions of the polarization by the MM region very accurately;67 however, it gives an indication of the general trends. VTST/MT Calculations: Calculation of the KIE. To calculate the KIE for the proton transfer (for doubly deuterated versus doubly protiated tryptamine as the substrate), VTST/ SCT37-39 were carried out with the CHARMMRATE interface.43 For each isotopic reaction, we calculated a representative free energy profile, which included the quantized-vibrational energy for the active zone atoms (atoms that were not restrained during the VTST calculation, see below) and the transmission coefficient, which accounts for quantum nuclear effects (including tunneling) on the reaction coordinate.28,37 The VTST/MT rate ) was obtained using eq 1 constant (kVTST/MT i

kVTST/MT (T) ) κMT i i (T)

kBT exp[-∆G TS,o (T)/RT] i h

(1)

where ∆G TS,o (T) is the standard-state molar free energy of i activation for the isotopic reaction i, with the quantizedvibrational energy included (thus, it is a quasi-classical activation free energy), and κMT is the transmission coefficient, calcui lated in this work within the SCT approximation. T, kB, h, and R are the temperature and Boltzmann’s, Planck’s, and gas constants, respectively. The KIE was then obtained by dividing the calculated perprotio kVTST/SCT rate constant by the one calculated for the doubly deuterated substrate reaction. The starting point for these VTST/SCT calculations was the structure corresponding to the maximum of the PM3/ CHARMM22 adiabatic mapping potential energy profile (see above). The QM/MM system was divided into two zones,28 corresponding to the QM atoms (active zone) and the MM atoms (bath zone). A saddle point, reactant, and product were located on the PM3/CHARMM22 potential energy surface by optimiz-

ing the QM atoms (active zone) embedded in the potential created by the frozen bath zone and characterized by normalmode analysis. From the saddle point, a minimum energy path (MEP) in an isoinertial mass-weighted Cartesian coordinate system was then calculated downhill for each isotopic reaction until the energy of the system approached that of the reactant and the product. The Page-McIver algorithm68 and a step size, δs, of 0.01 bohr amu1/2 were used (where s denotes the massweighted distance along the MEP; it is defined as zero at the saddle point, 0 on the product side). First and second derivatives of the energy were calculated at each step and every 10 steps of the MEP, respectively. The vibrational frequencies for the active zone and the PM3/ CHARMM22 free energy along the MEP, and thereby ∆G TS,o , i were obtained. However, as mentioned above, PM3 does not satisfactorily describe the energetics of this type of proton transfer2,34,60 and tends to underestimate the activation free energy when compared to experiment. Therefore, following the same approach used to correct the classical path, we carried out PM3-SRP/CHARMM22 single-point energy calculations along the PM3/CHARMM22 path, calculating a total of 18 energies along each MEP. The interpolated single-point energy (ISPE) corrections69 algorithm in POLYRATE70 was used to apply the correction to the . With the PM3/CHARMM22 potential energy and to ∆G TS,o i vibrational frequencies (νi(s)) and within the harmonic approximation, the zero-point energy contribution (EZPE(s)) was calculated at each point of the MEP with respect to the value at 3N-7 1 3N-6 1 the reactant (EZPE(s) ) ∑i)1 ( /2)hνi(s) - ∑i)1 ( /2)hνi(sR)) (where sR denotes the reactant of the VTST/SCT calculation). The sum of this zero-point energy contribution and the corrected potential energy gives the effective potential used in the calculation of the tunneling transmission coefficients, κMT i . We will refer to this path as R1. Further analysis of these MEPs was performed to provide insight into the origin of the elevated KIE observed experimentally. Tunneling Contribution to the Activation Free Energy. As explained above, the classical activation free energy for the proton transfer (∆G q,C) was obtained from umbrella sampling simulations, which were sampled over protein configurations. To include the effects of proton tunneling, we used the transmission coefficient obtained from the VTST/SCT calculations based on path R1. The contribution of hydrogen tunneling to the activation free energy (∆G q,tun)37,71 was obtained using eq 2

∆G q,tun(T) ) -RT ln[κMT i (T)]

(2)

This value was combined with the classical activation free energy to give the phenomenological activation free energy

∆G q,act(T) ) ∆G q,C(T) + ∆G q,tun(T)

(3) (∆G q,act)

Note that the phenomenological activation free energy should also include a term for the quantized-vibrational contribution.72 Although we have taken the quantized vibrations into account in the calculation of the KIEs (eq 1), we have not calculated the quantized-vibrational contribution to ∆G q,act. This contribution is usually found to lower the classical activation barrier significantly, and in the present case, it may be expected to be similar to that obtained for the related MADH/methylamine system (-3.2 kcal mol-1).34 The phenomenological activation free energy71,73 can be compared with the experimentally derived activation free energy. It predicts a rate constant for this step of


k(T) )

kB T exp[-∆Gq,act(T)/RT] h

Masgrau et al.

(4)

Multiple Paths for the VTST/SCT Calculations. To estimate the robustness of using the R1 path to analyze hydrogen tunneling in the reaction of AADH with tryptamine, we carried out VTST/SCT calculations on an ensemble of MEPs (proton transfer to O2). Following a scheme similar to that of Gao and Truhlar,42 the starting structures for these MEPs were taken from the umbrella sampling molecular dynamics window corresponding to the classical free energy maximum (ZO2 ) 0.0 Å for the present case). As before, the MEPs were calculated with the QM region (active zone) embedded in the potential created by the frozen protein environment; the Page-McIver algorithm68 and a step size, δs, of 0.01 bohr amu1/2 were used. The PM3SRP/CHARMM22 level of theory mentioned above was used here not only for energy calculations but also for geometry optimization and first and second derivatives of the energy. The vibrational frequencies for the active zone and the free energy along the MEP, and thereby ∆G TS,o (eq 1), were obtained for i each MEP. The individual KIEs, tunneling transmission coefficients, and tunneling contributions to the activation free energy (eqs 1 and 2) obtained for the different MEPs were used to calculate the corresponding averaged values. A total of four paths were used, which will be referred to as TS1-TS4. Due to convergence (technical) problems when using the HQtype link atoms in our calculations, a different QM/MM partition from the one used in the previous sections had to be used to calculate these MEPs. Generalized hybrid orbitals (GHOs)47 were used for the boundary atoms separating the QM and MM regions, which were positioned at the CR carbon atoms of Trp160β (second tryptophan residue forming the TTQ cofactor), Trp109β, and Asp128β. A total of 60 atoms formed the new QM region, comprising the full TTQ cofactor, bound tryptamine, and the catalytic aspartate (all of them up to the CR carbons). The following procedure was used to allow the system to accommodate the new QM/MM description. A patch was used to apply the new QM/MM partition to the last structure of the umbrella sampling molecular dynamics window corresponding at ZO2 ) 0.0 Å, followed by 200 of steepest descent and 500 steps of adopted basis Newton-Raphson (ABNR) energy minimization. Stochastic boundary molecular dynamics simulations were carried out as described above for a total of 60 ps. A restraint was applied to keep the ZO2 coordinate value close to that of the classical TS. Four snapshots from the last 40 ps of the simulation were taken as initial structures to localize the corresponding saddle points, reactants, and products and to calculate the four MEPs used for averaging. VTST/MT Calculations on the MADH/Methylamine System. To compare the two TTQ-dependent quinoproteins AADH and MADH, we have also carried out VTST/SCT calculations as described previously for the oxidation of methylamine by MADH.13,33,34,40-42,60 The MEP for abstraction by the O2 oxygen (equivalent to O2 in AADH) will be compared with R1 for the AADH/tryptamine system (see the Discussion section). A QM/MM partition with 31 atoms in the QM region and 3 HQ-type link atoms49 was used. The setup of the system, based on the methylamine-enzyme complex structure prepared in an earlier work,33,74 was done in an analogous way to that presented here for AADH (see ref 74 for a complete description). Single-point PM3-SRP/CHARMM22 energy calculations were carried out at 18 points along the PM3/CHARMM22 MEP to correct the potential energy surface. Rate constants for the

Figure 3. Averaged structure of the active site in (a) the reactant complex (III) (ZO2 ≈ -1.5 Å), (b) the TS (ZO2 ≈ 0.0 Å), and (c) the product (IV) (ZO2 ≈ 0.7 Å).

perprotio and the trideuterated methylamine reactions, and therefore the KIE, were obtained from eq 1. Results Simulation of the Classical Reaction Path. Reactant (III). In the averaged complex III structure resulting from the PM3/ CHARMM22 unrestrained molecular dynamics simulation of this complex (Figure 3a), the catalytic residue Asp128β is



TABLE 1: Average Acceptor-to-Hydrogen Distances (d(A-H)) for Some Key Hydrogen Bonds in the Reactant (III), TS, and Product (IV) Complexes from PM3-CHARMM MD Simulations A-H HE1-O7 HE1-Ala82β O HNT-O7 HNT-Asp84β O HNT-Thr172β OG1 O7-Asp84β HN HNW-Asn156β O O2-Thr172β HG1 O1-Trp160β HN O1-Cys161β HN a

d(A-H)a in reactant (III) ZO2 ≈ -1.5 Å 2.84 2.10 2.55 2.19 3.29 2.45 3.92 1.91 2.04 3.38

(0.13) (0.19) (0.11) (0.24) (0.26) (0.33) (0.34) (0.13) (0.17) (0.29)

d(A-H)a in TS ZO2 ≈ 0.0 Å 2.80 2.13 2.54 2.05 3.56 2.45 3.88 2.00 2.04 2.95

d(A-H)a in product (IV) ZO2 ≈ 0.7 Å

(0.09) (0.18) (0.11) (0.30) (0.32) (0.42) (0.39) (0.19) (0.13) (0.28)

2.80 2.11 2.50 2.10 3.67 2.37 3.82 2.68 2.17 2.90

(0.12) (0.18) (0.12) (0.20) (0.28) (0.34) (0.41) (0.78) (0.23) (0.24)

All distances are in Å. Values in parentheses indicate the standard deviation of the average distance.

positioned to abstract the proton from the bound tryptaminederived iminoquinone. The hydrogen bond between the proton to be abstracted (H1) and the O2 of the catalytic base was present throughout the simulation. The average donor-protonacceptor distances were d(C1-H1) ) 1.17 ( 0.04 Å, d(C1O1) ) 2.82 ( 0.07 Å, and d(O1-H1) ) 1.74 ( 0.07 Å. The other oxygen atom of the Asp128β carboxylate, O2, was also in close proximity to the proton to be abstracted throughout the simulation (d(C1-O2) ) 3.52 ( 0.20 Å and d(O2-H1) ) 2.68 ( 0.22 Å), and a hydrogen bond between H1 and O2 was observed 50% of the time. As mentioned in the Methods section, the AM1/CHARMM22 level of theory was also tested for modeling the reactant structure. In comparison to PM3, AM1 predicts a slightly different orientation of Asp128β with H1 located more centrally between the two carboxylate oxygens (d(O2-H1) ) 2.56 ( 0.33 Å and d(O1-H1) ) 2.07 ( 0.14 Å). According to the AM1/CHARMM22 simulation, the interaction of H1 with O1 would also be considered a hydrogen bond 100% of the time, but the interaction with O2 would be a hydrogen bond for ∼80% of the simulation. The difference in these distances predicted by the two semiempirical methods was also observed in their corresponding minimized structures, which showed distance values very close to the MD-averaged distances (d(O2-H1) ) 2.70 and 2.59 Å and d(O1-H1) ) 1.73 and 2.06 Å for PM3 and AM1, respectively). At the B3LYP/6-31G(d)/CHARMM22 level of theory, the geometry of the minimized reactant III (note that molecular dynamics of the modeled system at this level of theory were not computationally affordable) shows O1 to be significantly closer to H1 (d(O2-H1) ) 2.65 Å and d(O1H1) ) 1.87 Å). As the PM3/CHARMM22 geometry agrees better with the B3LYP/6-31G(d)/CHARMM22 one, the PM3/ CHARMM22 level was chosen to model the proton transfer. These differences obtained in the O1/2-H1 distances in the modeled complex III are most likely due to differences between the AM1 and the PM3 semiempirical methods, as AM1 favors bifurcated hydrogen bonds.75,76 A similar observationsthat d(O1-H1) is shorter than d(O2-H1)shas been made for the MADH/methylamine system.33,34,40-42,74,77 Several additional hydrogen bonds were observed between the substrate and the surrounding residues in complex III (Table 1). The backbone of Asp84β can interact with the substrate in two different waysseither donating a hydrogen bond (HNO7) or accepting a hydrogen bond (O-HNT). The hydrogen bond between HNT and Asp84β O has an average acceptorto-hydrogen (A-H) distance of 2.19 ( 0.24 Å, and the hydrogen bond between O7 and Asp84β HN has an average A-H distance of 2.45 ( 0.33 Å. A hydrogen bond is present between HE1 of TTQ and Ala82β O with an average A-H distance of 2.10 ( 0.19 Å. The carboxylate oxygens of Asp128β (O2 and O1) form

TABLE 2: Kinetic and Thermodynamica Parameters Calculated for the Protonb Transfer to O2 and O1 at the PM3-SRP//PM3/CHARMM22 Level O2 O1

∆G q,C

∆rG

12.5 16.4

-8.1 184 4.0 52

κH

∆G q,tun ∆G q,act H H -3.1 -2.4

9.4 14.0

kH

κD KIE

8.93 × 105 36 399 21

30 12

a ∆Gq,C, classical free energy barrier; ∆ G, classical free energy of r reaction; κH/D, tunneling transmission coefficient; ∆G q,tun H , tunneling contribution to the free energy barrier; ∆G q,act H , phenomenological free energy barrier; kH, rate constant (s-1); KIE, calculated kinetic isotope effect including the effect of proton tunneling. All energies are in kcal mol-1. b The subscripts H and D refer to the perprotiated and the perdeuterated reactions, respectively.

hydrogen bonds with Thr172β and Trp160β that are present in all structures of complex III with average A-H distances of d(O2-Thr172β HG1) ) 1.91 ( 0.13 Å and d(O1-Trp160β HN) ) 2.04 ( 0.17 Å. Analysis of the modeled reactant complex III shows that, although the abstraction by O1 of Asp128β would seem to be more favorable based on the proton-acceptor distance, ignoring the possible role of O2 as a proton acceptor could be misleading. The proximity between the donor carbon (C1) and both oxygen atoms of Asp128β in the modeled reactant complex (d(C1O2) ) 3.52 ( 0.20 Å and d(C1-O1) ) 2.82 ( 0.07 Å) is in agreement with the crystal structure of AADH with the inhibitor phenylhydrazine2 in which C1 is substituted by a nitrogen atom; thus some differences are expected (d(N-O2) ) 3.2 Å and d(N-O1) ) 2.8 Å). Moreover, while the positions of the hydrogen atoms can be deduced from the crystal structures, there is some ambiguity as to where they are positioned. Therefore, we investigated both possibilities, i.e., transfer to either oxygen atom, and the results are reported below. Classical Free Energy Barrier and Reaction Path. Umbrella sampling molecular dynamics simulations were used to calculate the classical free energy barrier and study the thermally averaged structural changes in and around the active site associated with proton transfer. Although initial adiabatic mapping profiles suggested transfer to O1 to be lower in energy, the PM3-SRP//PM3/CHARMM22 activation free energies (Table 2 and Figure 4) show that transfer to O2 is in fact energetically more favorable. The classical free energy barrier at 300 K for transfer to O2 is 12.5 kcal mol-1, which corresponds to a rate constant of ∼5000 s-1. For transfer to O1, a higher barrier of 16.4 kcal mol-1 (corresponding to a lower rate constant of 7 s-1) was obtained. The lower free energy barrier and more exothermic reaction (-8.1 kcal mol-1 for O2 vs 4.0 kcal mol-1 for O1) are observed for abstraction by O2, the oxygen more distant from the proton in the reactant complex (III). This preference for transfer to O2 did not change with the inclusion


Figure 4. Free energy profiles (relative to the reactant) for abstraction of H1 by either O2 or O1 of Asp128β at the PM3-SRP//PM3/ CHARMM22 level of theory. (Data from simulations without any restraints imposed on the reaction coordinate are also included in the WHAM analysis.)

of hydrogen tunneling in the activation free energy (see below). This shows that approaches based on substrate complexes alone (e.g., near attack conformations (NACs)) would be misleading in this case.78-80 While O1 is closer to the proton to be abstracted, the reorientation of Asp128β to shorten the O2H1 distance occurs at a relatively low energetic cost (Figure 4). This is possibly due to our observation that the interaction with O1 can be maintained while the interaction with O2 is established. Additionally, it appears that protonation of O2s but importantly not O1sis required for tryptamine oxidation to proceed to completion.2 Therefore, given that the energetics of proton abstraction by O2 are significantly more favorable than abstraction by O1, the results for abstraction by O2 are discussed below. Analysis of the simulations shows no major structural changes associated with the proton transfer to O2 (Figure 3). Figure 3b shows the averaged structure of the complex at ZO2 ≈ 0.0 Å, the location of the free energy maximum (TS). The average donor-proton-acceptor distances for this classical TS are d(C1-H1) ) 1.42 ( 0.11 Å, d(C1-O2) ) 2.77 ( 0.10 Å, and d(O2-H1) ) 1.40 ( 0.08 Å. Note that the reaction coordinate value of the simulation does not guarantee that each snapshot has this structure, as the force constant for the restraint is chosen so that there will be some overlap between adjacent simulations. The (classical) TS is analyzed not because it necessarily exists (tunneling will bypass this classical free energy maximum) but rather to give an idea of the structure and properties of representative reactive configurations. The classical TS also gives information about the potential energy surface, which determines the characteristics of the barrier through which the proton tunnels and thereby the tunneling probabilities. Moreover, computational studies of hydrogen tunneling have shown that tunneling is likely to occur once the donor-acceptor motions necessary for the (classical) transfer have taken place,33,34,37 i.e., when donor and acceptor have approached close enough so that the reaction coordinate becomes mainly hydrogen motion. This donor-acceptor distance is usually very close to the donoracceptor distance at the classical TS, even for cases where tunneling occurs at energies well below the classical TS.2 Figure 3c shows the averaged structure corresponding to the free energy minimum window at the product side of the proton transfer, which has average donor-proton-acceptor distances of d(C1-H1) ) 1.75 ( 0.04 Å, d(C1-O2) ) 2.73 ( 0.05 Å,

Masgrau et al. and d(O2-H1) ) 1.01 ( 0.03 Å and a reaction coordinate value of ZO2 ≈ 0.7 Å. According to these distances, the three atoms (C1, H1, and O2) remain at a separation within the limits of a weak C-H‚‚‚O hydrogen bond. VTST/SCT Calculations: Tunneling and Kinetic Isotope Effects. Rate Constant and KIE Calculations. Tunneling increases the rate of proton abstraction by O2 by a factor of 184 to 8.93 × 105 s-1 (Table 2 at the PM3-SRP//PM3/ CHARMM22 level and 300 K). In terms of activation free energy, this is equivalent to a reduction of 3.1 kcal mol-1 in the barrier height, giving a phenomenological activation free energy (∆G q,act in eq 3) of 9.4 kcal mol-1. For the perdeuterH ated reaction the rate enhancement due to tunneling is much smaller (36), which represents a lowering of the free energy barrier by 2.1 kcal mol-1. Note that this lowering of the free energy barrier provides a rough estimate of the maximum contribution of tunnelingsan equivalent reaction in solution might be expected to show some tunneling, thereby lowering (and possibly even eliminating) the contribution of tunneling to catalysis.81 The importance of tunneling in this reaction is also reflected in the KIE calculated using VTST/SCT; a value of 30 is obtained for proton transfer to O2 (27 for monodeuterated tryptamine). This value is 5 times higher than the calculated semiclassical (i.e., no tunneling included) KIE of 6 and thus clearly elevated significantly beyond the semiclassical limit. Other approaches have been used in the literature to calculate transmission factors.37 Application of the zero-curvature tunneling (ZCT) approximation to the AADH/tryptamine system gives a KIE of 26, whereas the Wigner correction82 gives a value of only 10. The ZCT approximation differs from the SCT approximation in that it does not consider the coupling to the reaction coordinate of vibrational normal modes transverse to it in the determination of the tunneling path, although these multidimensional vibrations are included in the effective potential for motion along the path. This coupling effect on the path itself has been found to be important in most hydrogen transfers,56,83 for which the ZCT treatment will underestimate the tunneling contribution. The Wigner correction omits not only the multidimensional effects (both kinds of coupling of vibrational modes) but also the actual shape of the potential energy surface used to calculate the tunneling probabilities, which is instead modeled by an inverse parabola. Furthermore, the Wigner correction approximates the tunneling through this model parabola by only the leading term in an expansion in powers of p.37 The adiabatic energy profiles, i.e., potential energy including zero-point energy corrections, used in the calculation of these tunneling probabilities and KIEs are shown in Figure 5a together with the representative tunneling energy (RTE; the energy of the system dominating the tunneling transfer).37 The RTE values for the protiated and deuterated reaction are 9.1 and 9.4 kcal mol-1, respectively, above the corresponding adiabatic energy of the reactant. The intersections of the RTE and the adiabatic energy barrier (classical turning points, Figure 5a) are indicative of the structures where the system changes from a classical to a tunneling regime. At the reactant side, this occurs at s ≈ -0.4 bohr amu1/2 (where d(H1-O2) ≈ 1.51 Å) for proton transfer and at s ≈ -0.5 bohr amu1/2 (where d(H1-O2) ≈ 1.52 Å) for deuteron transfer. The proton tunnels a distance of ∼0.46 Å (compared to ∼0.42 Å for the deuteron) between the two turning points. For proton transfer to the other carboxylate oxygen of Asp128β, O1, the calculated phenomenological free energy of



Figure 5. (a) Effective potentials (left-hand axis, solid lines with squares) along the reaction coordinate, s, used in the calculation of the tunneling transmission coefficients at the PM3-SRP//PM3/CHARMM22 level and curvature along the reaction path (right-hand axis, in bohr-1 amu-1/2, solid lines with crosses). s denotes the mass-weighted distance along the minimum energy path (defined as zero at the saddle point, 0 on the product side). Black and gray lines are used for the perprotio and the perdeutero reactions, respectively. The corresponding dominant tunneling energies at 300 K (representative tunneling energy (RTE), left-hand axis) are also shown as horizontal lines. (b) Zero-point energy contribution change with respect to the reactant zero-point energy contribution (right-hand axis, line with crosses) and variation of the C-O (triangles), C-H (circles), and H-O (squares) distances (in Å, left-hand axis) along the reaction path for protium transfer to O2. -1 activation (∆G q,act H , Table 2) is 14.0 kcal mol . The corre-1 sponding rate constant (400 s ) is therefore much lower than that for transfer to O2 (8.93 × 105 s-1). Thus, the inclusion of -1 tunneling (∆G q,tun H (O1) ) -2.4 kcal mol ) does not change the preference for transfer to O2 seen in the classical treatment of the reaction. The KIE is also smaller than that obtained for O2 (KIE ) 5 and 12 without and with tunneling, respectively). Therefore, our results show proton transfer to O2 to be kinetically and thermodynamically favored and the calculated KIE to be in closer agreement with experiment than transfer to O1. Analysis of the KIEs. Insight into the origin of the large calculated KIEs can be obtained from analysis of the MEP used to calculate the KIEs (Figure 5). The variationally optimized transition state structure for transfer to O2 occurs at s ) 0.06 bohr amu1/2, very close to the potential energy maximum (s ) 0.04 bohr amu1/2). The corresponding donor-proton-acceptor distances are d(C1-H1) ) 1.45 Å, d(C1-O2) ) 2.73 Å, and d(O2-H1) ) 1.29 Å, values in reasonable agreement with the range of distances obtained in the sampling of transition state configurations reported above. From the normal-mode analysis, imaginary frequencies of 2057i and 1514i cm-1 are obtained for the protiated and doubly deuterated tryptamine transition

Figure 6. Snapshots of the QM region along the calculated PM3/ CHARMM22 minimum energy path showing the change in the dominant movement associated with the reaction coordinate as the reaction proceeds: (a) s ) -2.0 bohr amu1/2, (b) s ) -0.7 bohr amu1/2, and (c) transition state.

states, respectively. The associated transition vector (Figure 6) illustrates the atomic motions that define the reaction coordinate at this structure; it is dominated by the motion of the proton being transferred (Figure 6c), but contributions from the atoms involved in the redistribution of the π electronic system that accompanies the reaction are also observed (i.e., C1, NT, CH2, and C7). Also depicted in Figure 6 are the eigenvectors associated with the reaction coordinate for two other points of the MEP. These enable visualization of the essential movement of the QM atoms in going from the reactant structure to the transition state. The first point was chosen to be an early structure in the reaction path (s ) -2.0 bohr amu1/2, Figure 6a); the second was chosen to correspond to the maxima of the reaction path curvature (see below, s ) -0.7 bohr amu1/2, Figure 6b), as this corresponds to a region of important changes in the motions that describe the reaction coordinate. Consistent with the changes of the

3040 J. Phys. Chem. B, Vol. 111, No. 11, 2007 donor-proton-acceptor distances (Figure 5) and the changes observed along the umbrella sampling reaction path, the MEP begins (from the reactant complex) with the approach of the two heavy atoms and the orientation of H1 toward O2. Rearrangement of the tryptamine indole, the NT atom (nitrogen of the iminoquinone) adjacent to the donor carbon, and the Asp128β carboxylate are observed in this initial part of the reaction path. Once the C1-O2 separation has reached a value of ∼2.75 Å (at s ≈ -0.6 bohr amu1/2), it remains nearly constant, while the proton would be (classically) transferred. The separation of C1 and H1 increases from 1.16 Å at s ≈ -0.6 bohr amu1/2 to 1.71 Å at s ≈ +0.6 bohr amu1/2, concomitant with d(O2-H1) decreasing from 1.61 to 1.01 Å. The observed change in the dominant motion from heavy-atom motion to predominantly proton motion is associated with a pronounced decrease of the vibrational frequency of mode 21 (from 2587 cm-1 at s ) -0.6 bohr amu1/2 to 1922 cm-1 at s ) -0.4 bohr amu1/2) This effects a rapid decrease in the zeropoint energy contribution (EZPE) for the protium reaction, with a change of -1.8 kcal mol-1 from s ) -0.6 to -0.4 bohr amu1/2 (Figure 5), with an overall difference in the zero-point energy contribution between the transition state and the reactant of -3.4 kcal mol-1. For the deuterated reaction the change in this frequency is only from 1970 to 1920 cm-1 (s ) -0.9 to -0.6 bohr amu1/2), and the associated zero-point energy contribution change from the VTST/SCT reactant to transition state is -2.3 kcal mol-1. This difference between the two isotopes, reflected in the adiabatic energy barriers depicted in Figure 5a, is the main factor contributing to the calculated semiclassical KIE of 6. Associated with this sudden change of the dominant motion, maxima in the reaction path curvature37 are also observed at s ) -0.7 and -1.0 bohr amu1/2 for protium and deuterium, respectively (Figure 5a). These maxima indicate the coupling of the reaction coordinate to vibrational modes transverse to it, a common feature in hydrogen-transfer reactions that has been shown to be very important for the correct treatment of tunneling.56,83,84 This multidimensional effect represents a shortening of the tunneling path and thus enhances proton tunneling. The effect here is larger for the proton transfer than for the deuterium transfer, which results in a further elevation of the KIE. An analysis of the vibrational modes at the proton curvature maxima shows that the reaction coordinate at this point is still dominated by donor-acceptor (heavy-atom) motion and that the major coupling is with the C1-H1 stretching normal mode (ν ) 2667 cm-1 at s ) -0.7 bohr amu1/2). Soon after this point and as the system approaches the semiclassical transition state structure, the C1-O2 stretching couples to the reaction coordinate, which is now dominated by proton motion. AADH/Tryptamine Results from Multiple Paths. The tunneling transmission coefficients (κH), tunneling contributions to the activation free energy (∆G q,tun H ), and KIEs calculated with the TS1-TS4 reaction paths are given in Table 3. Values within the ranges from 478 to 184 are obtained for κH, from -3.7 to -3.1 kcal mol-1 for ∆G q,tun H , and from 52 to 28 for the KIE. The average values are 295, -3.4 kcal mol-1, and 37, respectively, in good agreement with the values obtained with the R1 reaction path (184, -3.1 kcal mol-1, and 30, respectively). The adiabatic energy profiles for these MEPs are depicted in Figure 7. Paths TS1-TS4 are quite similar except for the energies of reaction (ranging from -6.4 to -1.9 kcal mol-1, exothermic in all cases) and the barrier heights (10.3 to 12.3 kcal mol-1). They all attain their corresponding reactant

Masgrau et al. TABLE 3: Tunneling Transmission Coefficients (KH), Tunneling Contributions to the Activation Free Energy -1 (∆G q,tun H , in kcal mol ) and KIEs Obtained from the TS1-TS4 MEPs (PM3-SRP/CHARMM22 Level of Theory) for the AADH/Tryptamine Systema TS1 TS2 TS3 TS4 average R1

κH

∆G q,tun H

KIE

478 317 202 184 295 184

-3.7 -3.4 -3.2 -3.1 -3.4 -3.1

52 38 31 28 37 30

a The respective averaged values and the corresponding values for the R1 MEP are also given.

Figure 7. Ensemble of minimum energy paths calculated for transfer to O2 in the AADH/tryptamine system (see the Methods section for details).

Figure 8. Comparison of the AADH/tryptmaine (black lines) and MADH/methylamine (gray lines) abstraction reactions by O2. Effective potentials (left-hand axis, solid lines with squares) along the reaction coordinate, s, used in the calculation of the tunneling transmission coefficients at the PM3-SRP//PM3/CHARMM22 level and curvature along the reaction path (right-hand axis, in bohr-1 amu-1/2, solid lines with crosses) are shown. The corresponding dominant tunneling energies at 300 K (representative tunneling energy (RTE), left-hand axis) are also depicted.

structures at s ≈ -3 bohr amu1/2, where d(O2-H1) ≈ 1.76 Å and d(O1-H1) ≈ 2.40 Å. For comparison, the R1 and the R1t paths are also shown in Figure 7; the R1t path was obtained by subtracting 4 kcal mol-1 from R1 at each point along the MEP and is only depicted to facilitate comparison with the TS1TS4 paths. KIE and Potential Energy Profile for MADH/Methylamine. The PM3-SRP//PM3/CHARMM22 potential energy profile (including zero-point energy contributions) for protium transfer to O2 in MADH (equivalent to O2 in AADH) is depicted in Figure 8 together with that for the reaction in AADH. The adiabatic energy maximum occurs at s ) +0.1 bohr amu1/2 and



TABLE 4: Average Energiesa and Residue Contributions to the QM/MM Electrostatic Interaction from Analysis of Structures from PM3-CHARMM MD Simulations reactant (III) (ZO2 ≈ -1.5 Å)

TS (ZO2 ≈ 0.0 Å)

product (IV) (ZO2 ≈ 0.7 Å)

QM QM/MMelec

39.1 -58.9

(4.6) (5.9)

0.0 0.0

34.6 -46.7

(6.8) (4.9)

-4.5 12.2

5.6 -31.5

(5.9) (3.7)

-33.5 27.4

Ala82β Asp84β HNb Asp84β CdOb Asp84β CO2- b Asn156β Asp157β Trp160β Cys161β Phe169β Thr172β

-4.1 -0.1 -5.5 -10.5 -1.1 -9.9 -13.1 -1.0 -2.5 -13.9

(0.9) (0.8) (1.3) (2.4) (0.9) (1.1) (2.1) (0.9) (0.5) (1.7)

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

-3.3 -0.6 -5.4 -9.6 -1.4 -7.7 -10.1 -1.8 -1.4 -8.8

(0.9) (0.9) (1.0) (1.6) (1.0) (1.1) (2.0) (0.7) (0.4) (1.6)

0.8 -0.4 0.2 0.9 -0.3 2.2 3.0 -0.8 1.0 5.1

-2.0 -1.7 -5.2 -6.7 -3.1 -5.3 -5.9 -1.6 -0.3 -3.4

(0.8) (1.2) (1.0) (1.4) (0.7) (1.0) (1.8) (0.7) (0.3) (1.6)

2.1 -1.5 0.3 3.9 -2.0 4.5 7.2 -0.6 2.1 10.5

a The QM energy is the energy of the QM region only at the PM3 level, QM/MMelec is the electrostatic interaction between the QM and MM regions, and the energies associated with a particular residue are the contributions of that residue to the QM/MMelec energy. Negative energies indicate a stabilizing contribution, and positive energies a destabilizing contribution. All energies are in kcal mol-1. Numbers in parentheses indicate the standard deviations, and numbers in bold show the energies relative to that in the reactant complex. b Residue Asp84β interacts with both the tryptamine substrate and the TTQ cofactor (see text). To provide a more detailed description of the interactions of this residue along the reaction coordinate, this reside is divided into three separate electrostatic groups: the carboxylate side chain (CO2-), the NH group of the backbone (NH), and the carbonyl group of the backbone (CdO).

corresponds to a value of 15.3 kcal mol-1, very close to that obtained for the AADH path (15.5 kcal mol-1). Despite the similarities in the barrier height, the curvature of the reaction path, and the energy profile of the two systems (Figure 8), the calculated KIE for MADH is much lower than that of AADH (11 vs 30). The semiclassical KIE is 6. This difference arises from the fact that the MADH MEP is more endothermic, and as a result, the proton cannot tunnel until higher energies are attained as shown by the representative tunneling energies (Figure 8). The reactant-side turning point is found at s ≈ -0.4 bohr amu1/2. The transmission coefficients for MADH are 35 and 17 for protium and deuterium, respectively. This corresponds to tunneling contributions to the activation free energies (eq 2) of -2.1 and -1.7 kcal mol-1, respectively. Analysis of the Energetics and Interactions along the Reaction Path. Hydrogen Bond Analysis. Analysis of the hydrogen bonding along the ZO2 ) [d(C1-H1)-d(O2-H1)] abstraction path, i.e., for abstraction by O2 of Asp128β, suggests that no major changes in conformation are associated with proton transfer (Table 1). The hydrogen bond between HE1 and Ala82β O is present in all structures at each reaction coordinate value (d(HE1-Ala82β O) ) 2.10 ( 0.19 Å at ZO2 ≈ -1.5 Å and d(HE1-Ala82β O) ) 2.11 ( 0.18 Å at ZO2 ≈ 0.7 Å). Both hydrogen bonds formed between Asp84β and atoms of the QM region are present throughout the path and are on average slightly shorter in the product (IV) than in the reactant (III) (d(HNT-Asp84β O) ) 2.19 ( 0.24 Å and d(O7-Asp84β HN) ) 2.45 ( 0.33 Å at ZO2 ≈ -1.5 Å and d(HNT-Asp84β O) ) 2.10 ( 0.20 Å and d(O7-Asp84β HN) ) 2.37 ( 0.34 Å at ZO2 ≈ 0.7 Å). One might expect the hydrogen bonds involving O2 and O1 of Asp128β to be significantly affected by the protonation of O2 in this pathway. However, the hydrogen bond between O1 and the Trp160β HN backbone is present in all structures along the reaction coordinate. The acceptor-to-hydrogen distance lengthens slightly upon proton transfer, changing from d(O1Trp160β HN) ) 2.04 ( 0.17 Å in the reactant complex (III) to d(O1-Trp160β HN) ) 2.17 ( 0.23 Å in the product (IV). The hydrogen bond between O2 and Thr172β HG1 observed in the reactant complex is affected by the protonation of O2, but the hydrogen bond is still present in the majority of the structures of the product. d(O2-Thr172β HG1) increases from

an average of 1.91 ( 0.13 Å in the reactant (III) to 2.68 ( 0.78 Å in the product (IV). The larger standard deviation of this distance in the product (compared to that in the reactant complex) shows that Thr172β can move out of hydrogenbonding range with O2. This is due to rotation of the Thr172β hydroxyl toward the TTQ ring and O7 in some structures of the product. Effect of the Enzyme on the Energetics of the Reaction. An energy decomposition analysis was carried out for the structures from the simulations of reactant (III) without any restraints imposed on the reaction coordinate (ZO2 ≈ -1.5 Å), the free energy maximum (ZO2 ≈ 0.0 Å), and the product (ZO2 ≈ 0.7 Å) simulations (Table 4). On the basis of the QM component of the total energy, proton transfer to O2 in the active site of AADH is an essentially barrierless and exothermic process (∆QMq ) -4.5 ( 6.8 kcal mol-1 and ∆rQM ) -33.5 ( 5.9 kcal mol-1 (relative to the reactant complex)). Note that these energies use the standard PM3/CHARMM22 method without the correction being applied. We would expect the correction to increase the QM barrier by ∼4 kcal mol-1. A significant energy barrier for the proton transfer only arises when the QM/ MM electrostatic interactions are considered (∆QM/MMqelec ) +12.2 ( 4.9 kcal mol-1 and ∆rQM/MMelec ) +27.4 ( 3.7 kcal mol-1 (relative to the reactant complex)). This reflects the fact that the maximum electrostatic stabilization of the QM atoms by the surrounding protein occurs at the reactant structure. As mentioned earlier, since we have not also studied the reaction in aqueous solution we cannot assess the contribution to catalysis of these interactions. The main electrostatic interactions in the reactant complex (III) are hydrogen bonds to the Asp128β carboxylate from the side chain of Thr172β and the backbone of Trp160β, with average (absolute) QM/MMelec interaction energies of -13.9 ( 1.7 and -13.1 ( 2.1 kcal mol-1, respectively. These interactions will be affected the most by the reaction since the electron density of the Asp128β carboxylate decreases as the proton is transferred. This stabilization of the carboxylate through hydrogen bonds likely prevents the (already elevated) pKa (7.1 in the free enzyme and 6.0 in the enzyme-substrate complex85) of Asp128β from being elevated even further, thus ensuring the presence of the charged (deprotonated) form of Asp128β.

3042 J. Phys. Chem. B, Vol. 111, No. 11, 2007 There are two other aspartate residues in the active site, Asp84β and Asp157β, which also have a significant electrostatic interaction with the QM region in the reactant. Asp157β forms a hydrogen bond with HNW of the tryptamine adduct in some structures along the reaction coordinate. The absolute QM/ MMelec energy of Asp157β in the reactant (III) is -9.9 ( 1.1 kcal mol-1. The interaction of Asp84β is split into separate groups as this residue forms two hydrogen bonds with the QM atoms in the reactant complex: HNT-Asp84β O and O7Asp84β HN. The first part contains the electrostatic group of the backbone nitrogen (Asp84β HN), the second part contains the electrostatic group of the backbone carbonyl (Asp84β CdO), and the third group contains the carboxylate side chain (Asp84β CO2-). The carboxylate side chain of Asp84β has the largest absolute QM/MMelec energy of the different parts of Asp84β (absolute QM/MMelec ) -10.5 ( 2.4 kcal mol-1), despite the fact that it is not involved in any hydrogen bonds with atoms of the QM region. The nitrogen and carbonyl backbones have absolute electrostatic interaction energies of -0.1 ( 0.8 and -5.5 ( 1.3 kcal mol-1 in the reactant (III), respectively. In general, the interactions that stabilize the reactant (III) diminish upon proton transfer. The most significant reduction in electrostatic interaction is observed for residues Thr172β and Trp160β that are hydrogen-bonded to O2 and O1 of Asp128β, respectively. O2 is protonated during the reaction, and consequently, the contribution of Thr172β shows the most significant change of any residue upon proton transfer with a relative QM/ MMelec energy of +10.5 ( 1.6 kcal mol-1 in the product (IV). The decrease in interaction energy is consistent with the significant (0.77 Å) increase in acceptor-to-hydrogen (A-H) distance upon proton transfer observed for this hydrogen bond. The relative QM/MMelec energy of Trp160β is +7.2 ( 1.8 kcal mol-1 in the product (IV). The hydrogen bond between O1 and Trp160β HN is also longer in the product, but the change in A-H distance (0.13 Å) for this hydrogen bond is not as large as that observed in the hydrogen bond between O2 and Thr172β. For the TS and the product, only Cys161β, Asn156β, and HN-Asp84β provide more QM/MM stabilization than for the reactant (III) (on average a relative stabilization of 99% of the proton transfer goes via a combination of classical activation and proton tunneling instead of via a purely classical route. The resulting phenomenological activation free energy for transfer to O2, ∆GHq,act, is 9.4 kcal mol-1. The contribution of tunneling to catalysis is most usefully addressed by comparison to an equivalent uncatalyzed reaction, preferably in solution.22 An experimental comparison of this sort would prove challenging for the proton transfer studied here. Computational analysis is therefore particularly useful. In principle, the energetic differences found between the two environments could lead to a different role for tunneling. As discussed in ref 7, tunneling can be more important in the enzyme than in solution because the slow step in the uncatalyzed reaction in solution is often very endothermic. Alternatively, in the cases where the reaction energy is not affected and the barrier height is lower in the enzyme, the tunneling contribution could even be smaller in the enzyme than that in solution.7 So far, for the systems where such a comparison has been made, both environments led to similar tunneling contributions, suggesting no catalytic role for hydrogen tunneling in those cases.81 While we have not examined an equivalent solution reaction for the present system (note that others have unsuccessfully tried for the MADH system),41 our calculations show that the activation barrier is reduced by 3.1 kcal mol-1 relative to a purely classical proton transfer in the enzyme. This value provides a rough estimate of the maximum catalytic contribution of tunneling (i.e., as discussed above, an equivalent reaction in solution might be expected to also show tunneling effects, thereby lowering (and possibly even eliminating) the contribution to catalysis of tunneling81). If the reaction in solution showed no tunneling, then the maximum contribution to catalysis of tunneling would be significant (a maximum increase in rate of 2 orders of magnitude) but comparable and indeed less than other catalytic effects7,92 found in other types of enzyme reactions (e.g., transition state stabilization by electrostatic interactions in chorismate mutase).44 It is clear, however, that the well-organized active site of AADH, set up for efficient proton transfer in this reactive complex, also has a structure capable of promoting tunneling.2 The analysis of the representative MEP used to calculate the tunneling contributions shows that the proton tunnels ∼0.46 Å through the barrier from a H1-O2 distance of 1.51 Å. The significant reaction path curvature at this region of the MEP reveals the coupling of the C1-O2 stretching mode to the reaction coordinate as the latter becomes dominated by proton motion. This coupling is translated into an effective reduction of the length of the tunneling path, which enhances the tunneling contribution to the rate constant. (Although it is generally associated with the promotion of tunneling,23 it should be noted that two recent studies with lipoxygenase suggest that such a shortening of the proton-acceptor distance can decrease tunneling93 and that there is an intermediate tunneling path length that corresponds to the maximum for the KIE value.94) If the coupling is omitted, as in the ZCT approximation, then the calculated transmission coefficient is reduced from 184 (with SCT) to 100 (with ZCT). The importance of accounting for this multidimensional effect, termed corner cutting, in the calculation of rate constants and KIEs is well-known for gas-phase

Masgrau et al. reactions,83,95 and it has now been shown to also be important for enzyme-catalyzed proton transfer.33,39,56 The significance of proton tunneling in this reaction is also apparent from the calculated KIE, which, with a value of 30 with the SCT approach, is 5 times higher than the semiclassical calculated value of 6. The effect of the coupling between the reaction coordinate and vibrational modes transverse to it is also observed in the value of the calculated KIEs: If the coupling is omitted (ZCT approach), then the calculated KIE is 13% smaller than the SCT value. It is also instructive to note the (expected) much poorer performance of the Wigner tunneling correction, which gives a KIE 67% smaller than the SCT approach. Therefore, although the Wigner correction is readily calculated, its use is not recommended, especially when tunneling is significant and occurs from energies considerably below the top of the energy barrier, as is the present case. VTST/SCT Results from Different MEPs. The averaged KIEs and tunneling parameters obtained from the TS1-TS4 paths (κH ) 295, ∆G q,tun ) -3.4 kcal mol-1, and KIE ) 37) are in H good agreement with those obtained with the R1 MEP (184, -3.1 kcal mol-1, and 30, respectively). This suggests that R1 is a representative path and that the analysis described above can provide insight into tunneling in this system. The structures for the four paths TS1-TS4 were selected at times of 20, 30, 37, and 58 ps, respectively, of the molecular dynamics simulation. Although four snapshots are not enough to draw a definitive conclusion, the results (Table 3) could indicate a lack of equilibration of the system at the earlier times and a convergence of the KIEs to ∼30. It is important to note that for the R1 MEP a protein (MM region) environment frozen at a configuration close to that at the reactant complex was used, while the TS1-TS4 MEPs had the MM region frozen at a configuration corresponding to the transition state ensemble. The use of protein configurations taken from different parts of the reaction coordinate has been shown to significantly affect the calculated KIEs.40 Although this is not the case here, some differences exist. The use of a more reactant-like configuration in R1 is most probably responsible for an overestimation of the endothermicity, whereas the use of a TS configuration in TS1-TS4 leads to more exothermic paths than R1. Another interesting effect appears as a result of the different environmentspath R1 portrays the reaction with H1 initially closer to O1 than to O2, proceeding from the modeled reactant complex (III) to the product (IV). However, when the protein is frozen at a TS configuration, the MEP leads to a reactant structure with H1 closer to O2 than to O1 (with d(O2-H1) ≈ 1.76 Å and d(O1-H1) ≈ 2.40 Å). Due to the geometric constraints imposed by the frozen environment, this structure is not able to relax into the modeled reactant complex (III) structure. This leads to the apparent significant difference between the R1 and the TS1-TS4 MEPs seen in Figure 7, especially for the barrier heights and the reactant side of the energy profiles. Overlaying R1 (i.e., R1t) and TS1-TS4 illustrates that the regions from s ≈ -0.5 to 0.5 bohr amu1/2 are quite similar. This region comprises the zone where tunneling is possible and most likely to occur, which explains why no large differences are obtained between the R1 and the TS1-TS4 results, with the effects of differences in barrier height, width and shape, and exo/endothermicity essentially cancelling out. Energy Corrections. The standard PM3 method does not model the energetics of this proton transfer very accurately. The activation barrier is too low, when compared with experiment,

Tunneling and Classical Paths for Proton Transfer and the product is very exothermic in terms of QM energy. As a result it is necessary to correct the energies in some way to account for this. Our approach here was to use SRPs developed by Alhambra and co-workers for the analogous MADH/ methylamine system.34 The energy difference between PM3/ CHARMM22 and PM3-SRP/CHARMM22 values for singlepoint energy calculations was used to adjust the free energy at each point along the ZO2 or ZO1 reaction coordinate. The ∆G q,C of 12.5 kcal mol-1 for transfer to O2 represents an increase of ∼+4 kcal mol-1 compared to PM3/CHARMM22, and the ∆rG of -8.1 kcal mol-1 represents a reduction in the exothermicity predicted by PM3/CHARMM22 of ∼6 kcal mol-1. The application of these corrections does not alter the preference of the two oxygenssfor O1 this increase is in fact slightly larger, at ∼+5 and ∼+6.5 kcal mol-1 for the barrier height and reaction energy, respectively. The use of the SRPs makes transfer to O1 endothermic rather than approximately thermoneutral. Similar corrections were obtained for the tunneling calculations when single-point energy calculations were carried out on the MEP structures; in particular, for transfer to O2, corrections were obtained of +4.4 kcal mol-1 (+5.0 kcal mol-1 for O1) for the barrier height and +5.7 kcal mol-1 (+7.0 kcal mol-1 for O1) for the reaction energy. For the preferred transfer to O2, the contribution of tunneling to the activation free energy (eqs 2 and 3) is increased by 19% when this correction is used, with the standard PM3 method giving a ∆G q,tun value of -2.6 H kcal mol-1. The effect on the KIE is even larger, with an increase of 25% over the already high standard PM3 value of 24. However, the improvement in the calculated KIE is still not enough to reproduce the extremely high experimental value of 55 ( 6. The phenomenological activation free energy (a ∆G q,act value of 9.4 kcal mol-1) obtained also underestiH mates the experimental value of 12.6 kcal mol-1. Taken together, these results seem to indicate that the correction used for the barrier height may be too small. Another less rigorous method of correction applied in previous work2 to give an indication of the activation barrier was to use corrections obtained from the energy differences between PM3 and PM3-SRP calculations reported for MADH/ methylamine (without carrying out any further calculations). This gave an increase in barrier height of 8.7 kcal mol-1 and an increase in reaction energy of 6.2 kcal mol-1, leading to a ∆Gq,C value of 17.2 kcal mol-1 for O2 and 20.0 kcal mol-1 for O1. When these values were corrected to include tunneling value of 13.2 kcal mol-1 was obtained for effects, a ∆G q,act H O2, in reasonable agreement with the experimental value. It should be noted that the effect of quantizing the vibrational energy,72 which would probably lower the activation free energy further, is not explicitly included in the value that we report for ∆G q,act H . The use of a much higher correction for the energy barrier than (and a similar correction for the reaction energy to) that used in the present study led to larger tunneling transmission coefficients and a higher KIE.2 Our previous estimate of the KIE for the proton transfer to O2 was 93, thus overestimating the experimental value. The use of different corrections has allowed us to test the sensitivity of the results on the potential energy surface. The different models used seem reasonable and lead to similar conclusions, in particular the importance of tunneling in this enzymic proton-transfer reaction. Large KIEs have been obtained in all of the calculations, elevated well above the semiclassical limit, in agreement with the experimental observations. Comparison with the MADH/Methylamine System. It is informative to compare the results for AADH/tryptamine with

J. Phys. Chem. B, Vol. 111, No. 11, 2007 3045 those for the related MADH/methylamine system. The earliest studies of MADH,33,60 based on the crystal structure of the free enzyme, showed one oxygen atom of the catalytic Asp428β of MADH (equivalent to Asp128β in AADH) to be hydrogenbonded to Tr2 of TTQ (equivalent to Trp160β in AADH) and a water molecule, while the other carboxyl oxygen was not hydrogen-bonded. This led to consideration of only the unbonded carboxyl oxygen (equivalent to O2 in AADH) as the proton acceptor in methylamine and ethanolamine oxidation. In a subsequent study, Tresadern et al.40 found that the water molecule moved away during a molecular dynamics simulation based on the modeled reactant complex, with a subsequent repositioning of Asp428β to form a hydrogen bond with Thr474β, equivalent to the hydrogen bond formed between Asp128β and Thr172β in the crystal structure of AADH/ tryptamine.2 In AADH, a water molecule is present in the same position in the crystal structure of the Michaelis complex2 (I in Figure 1) but has moved away before the proton transfer (Figure 1).2 This highlights the care that must be taken when constructing a model for simulations, especially when basing the model on coordinates from a ligand-free enzyme. The semiclassical KIE obtained in all of the MADH/ methylamine studies (∼5-6),33,34,41,60 including our own recent study,74 is very similar to that reported here for the AADH/ tryptamine reaction. Thus, the much more highly inflated KIE for the AADH/tryptamine system appears to arise mainly from the tunneling contribution and not from differences in the zeropoint energy contributions. When tunneling is included in the reaction, the general picture emerges of heavy atoms approaching close enough for the proton to tunnel while they remain nearly fixed applies in all cases. We have obtained similar values of the reaction coordinate, s, corresponding to the turning points for MADH and AADH; these are also similar to those reported by Tresadern et al.60 However, the adiabatic energy barrier beyond (i.e., to the product side of) these points seems to be sharper for the AADH/tryptamine system, which could explain the larger tunneling contribution obtained for AADH/tryptamine. Comparison of the MEPs for AADH and MADH calculated in our study shows an almost identical energy profile from s ≈ -0.4 to +0.5 bohr amu1/2. The difference between the two systems is the higher endothermicity for MADH, which, in turn, reduces the tunneling region available for this enzyme, resulting in tunneling through a barrier of similar width but lower height. The effect of the reaction energy on the tunneling probabilities is well-known and has been highlighted recently41 in the context of the MADH/methylamine system. Note that the differences in the reactant side of the reaction in Figure 8 are below the energy of the product, which is the energy at which tunneling can commence in an endothermic reaction, and do not impact greatly on the tunneling region. In a recent study of MADH and in agreement with our potential energy surfaces for AADH, Nun˜ez et al.41 obtained lower potential energy barriers for transfer to the closer oxygen atom (equivalent to O1 in AADH).41 Also in accord with AADH,2 the authors find that the tunneling contribution and the KIE are higher and in better agreement with experiment for transfer to the oxygen equivalent to O2 in AADH. Thus, Nun˜ez et al. conclude that while significant proton transfer to the oxygen equivalent to O2 occurs by tunneling in the MADH/ methylamine system it is likely that transfer via a classical overthe-barrier mechanism would also occur to the other oxygen atom (equivalent to O1). The study of MADH by Alhambra et al.34 appears to consider the oxygen atom equivalent to O2 in AADH, although it is not specified explicitly in the text and no

3046 J. Phys. Chem. B, Vol. 111, No. 11, 2007 hydrogen-bonded water molecule is mentioned. Another recent study77sbased on the analysis of a molecular dynamics simulation of the MADH/methylamine reactant complexsconcludes that the oxygen atom of the catalytic aspartate equivalent to O1 in AADH is the only proton acceptor. It should be noted that for AADH the umbrella sampling molecular dynamics calculations predict a lower free energy barrier for transfer to O2 and thus show O2 to be the favored acceptor even before accounting for hydrogen tunneling. Umbrella sampling molecular dynamics results obtained by us74 for the MADH/methylamine system at the same PM3-SRP//PM3/CHARMM22 level give similar classical free energy barriers for the two possible proton transfers. Thus, the question of whether the equivalent of O2 and/or O1 is the proton acceptor in the MADH/ methylamine system remains to be answered; studies to address this, which are currently underway in our groups, will include both dynamics of the enzyme system and the reaction itself, incorporating tunneling effects. In most computationally studied proton transfers in which a carboxylate group acts as the catalytic base only one oxygen is considered; this assumption is made based on geometrical grounds (shortest distance) underpinned by the assumption that both oxygens are otherwise equivalent. We have shown here that this is not necessarily always the case. A different substrate or a slightly different active site could alter this preference. For the AADH/tryptamine system it seems that abstraction by O2 would be required sterically for the reductive half-reaction to proceed to completion.2 This restriction is not necessarily present in MADH for the much smaller methylamine substrate. We could then speculate that the different interplay between the two oxygens abstracting the proton, which present significantly different KIEs in both the MADH and the AADH systems, could be in part responsible for the differences observed in the KIEs. However, and as suggested for MADH,41 the possible competition between the two oxygens and how this changes with temperature would add kinetic complexity to a system with an already high number of degrees of freedom. Conclusions We have used QM/MM molecular dynamics simulations and VTST/SCT calculations to study in detail both the classical and the tunneling paths for the proton transfer catalyzed by AADH with tryptamine as a substrate, for which an experimental KIE of 55 ( 6 (T ) 5-20 °C) has been reported. After the recent publication of the enzyme crystal structure, we have extended our previous computational study on AADH2 by using improved energy corrections, high-level QM/MM methods, and VTST/ SCT results averaged over several MEPs. The modeled AADH/tryptamine reactant complex (iminoquinone adduct) III is configured with the proton much closer to the O1 oxygen of Asp128β than to O2. The ab initio QM/ MM characterization of this reactant complex agrees with this result. Despite this, and contrary to what would be initially expected, transfer to the O2 oxygen of Asp128β is found to be the favored reaction both thermodynamically and kinetically. Therefore, we have focused our analysis on this transfer. Analysis of the Mulliken charges and the enzyme-substrate interactions has revealed that the C-H bond to be broken is polarized as a result of the close presence of the catalytic base (Asp128β) and the (formally) positively charged imine nitrogen. The interactions of the QM region (Asp128β-tryptamine-TTQ) with the active site residues tend to favor the reactant more than the product, mainly as a result of the two strong hydrogen bonds made by the carboxylate of Asp128β. An important driving force

Masgrau et al. for the reaction comes from the energy change in the QM region itself as the proton is transferred and the electron delocalized into the TTQ rings. A calculated KIE of 30 has been obtained (KIE ) 6 when tunneling is omitted), elevated well above the semiclassical limit, in agreement with the experimental observations. Comparison of quantum and classical models of proton transfer allows estimation of the contribution of hydrogen tunneling in lowering the barrier to reaction in the enzyme. A contribution of -3.1 kcal mol-1 has been obtained, which represents a rate enhancement due to tunneling by 2 orders of magnitude. This value provides a rough estimate of the maximum catalytic contribution of tunnelingsan equivalent reaction in solution could invoke tunneling, thereby lowering (and possibly even eliminating) the contribution of tunneling to catalysis. Therefore, although the maximum contribution to catalysis of tunneling calculated in this work is quite significant, it is comparable tosand indeed somewhat less thanssome other catalytic effects found in other types of enzyme reactions. For comparison, VTST/SCT calculations have also been performed on the proton tunneling reaction in the related MADH/methylamine system, for which an experimental KIE of ∼17 has been reported. The much lower calculated KIE that we obtain for the MADH/methylamine system (KIE ) 11) is attributed to a more endothermic potential energy surface for this reaction, although further studies of the MADH/methylamine system are required to confirm which oxygen(s) of the catalytic aspartate act as acceptor(s) for the tunneling proton. Acknowledgment. L.M. and K.E.R. contributed equally to this work. The authors thank the Engineering and Physical Sciences Research Council (EPSRC) for support (Grant Nos. GR/R86102/01 and GR/R86119/01). A.J.M. also thanks the IBM High Performance Computing Life Sciences Outreach Programme and the Biotechnology and Biological Sciences Research Council (BBSRC). M.J.S. also thanks the BBSRC Research Equipment Initiative for funding a high-performance computer cluster, and M.J.S. and L.M. gratefully acknowledge use of the University of Leicester Mathematical Modelling Centre’s supercomputer, which was purchased through the EPSRC Strategic Equipment Initiative. N.S.S. is a BBSRC Professorial Research Fellow. We are grateful to Linus Johannissen for useful discussions. References and Notes (1) Ball, P. Nature 2004, 431, 396. (2) Masgrau, L.; Roujeinikova, A.; Johannissen, L. O.; Hothi, P.; Basran, J.; Ranaghan, K. E.; Mulholland, A. J.; Sutcliffe, M. J.; Scrutton, N. S.; Leys, D. Science 2006, 312, 237. (3) Masgrau, L.; Basran, J.; Hothi, P.; Sutcliffe, M. J.; Scrutton, N. S. Arch. Biochem. Biophys. 2004, 428, 41. (4) Liang, Z. X.; Lee, T.; Resing, K. A.; Ahn, N. G.; Klinman, J. P. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 9556. (5) Liang, Z. X.; Klinman, J. P. Curr. Opin. Struct. Biol. 2004, 14, 648. (6) Kohen, A.; Klinman, J. P. Chem. Biol. 1999, 6, R191. (7) Garcia-Viloca, M.; Gao, J.; Karplus, M.; Truhlar, D. G. Science 2004, 303, 186. (8) Olsson, M. H.; Siegbahn, P. E.; Warshel, A. J. Am. Chem. Soc. 2004, 126, 2820. (9) Hammes-Schiffer, S. Biochemistry 2002, 41, 13335. (10) Antoniou, D.; Caratzoulas, S.; Kalyanaraman, C.; Mincer, J. S.; Schwartz, S. D. Eur. J. Biochem. 2002, 269, 3103. (11) Sutcliffe, M. J.; Scrutton, N. S. Eur. J. Biochem. 2002, 269, 3096. (12) Hwang, J. K.; Warshel, A. J. Am. Chem. Soc. 1996, 118, 11745. (13) Basran, J.; Sutcliffe, M. J.; Scrutton, N. S. Biochemistry 1999, 38, 3218. (14) Basran, J.; Patel, S.; Sutcliffe, M. J.; Scrutton, N. S. J. Biol. Chem. 2001, 276, 6234. (15) Harris, R. J.; Meskys, R.; Sutcliffe, M. J.; Scrutton, N. S. Biochemistry 2000, 39, 1189.

Tunneling and Classical Paths for Proton Transfer (16) Basran, J.; Harris, R. J.; Sutcliffe, M. J.; Scrutton, N. S. J. Biol. Chem. 2003, 278, 43973. (17) Kohen, A.; Cannio, R.; Bartolucci, S.; Klinman, J. P. Nature 1999, 399, 496. (18) Knapp, M. J.; Rickert, K.; Klinman, J. P. J. Am. Chem. Soc. 2002, 124, 3865. (19) Francisco, W. A.; Knapp, M. J.; Blackburn, N. J.; Klinman, J. P. J. Am. Chem. Soc. 2002, 124, 8194. (20) Agrawal, N.; Hong, B.; Mihai, C.; Kohen, A. Biochemistry 2004, 43, 1998. (21) Warshel, A. Computer Modeling of Chemical Reactions in Enzymes and Solutions; John Wiley & Sons: New York, 1991. (22) Villa, J.; Warshel, A. J. Phys. Chem. B 2001, 105, 7887. (23) Benkovic, S. J.; Hammes-Schiffer, S. Science 2003, 301, 1196. (24) Mulholland, A. J. Drug DiscoVery Today 2005, 10, 1393. (25) Claeyssens, F.; Harvey, J. N.; Manby, F. R.; Mata, R. A.; Mulholland, A. J.; Ranaghan, K. E.; Schu¨tz, M.; Thiel, S.; Thiel, W.; Werner, H.-J. Angew. Chem., Int. Ed. 2006, 45, 6856. (26) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 227. (27) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11, 700. (28) Truhlar, D. G.; Gao, J.; Alhambra, C.; Garcia-Viloca, M.; Corchado, J.; Sa´nchez, M. L.; Villa`, J. Acc. Chem. Res. 2001, 35, 341. (29) Hammes-Schiffer, S. Curr. Opin. Struct. Biol. 2004, 14, 192. (30) Mure, M.; Mills, S. A.; Klinman, J. P. Biochemistry 2002, 41, 9269. (31) Hyun, Y. L.; Davidson, V. L. Biochemistry 1995, 34, 816. (32) Govindaraj, S.; Eisenstein, E.; Jones, L. H.; Sanders-Loehr, J.; Chistoserdov, A. Y.; Davidson, V. L.; Edwards, S. L. J. Bacteriol. 1994, 176, 2922. (33) Faulder, P. F.; Tresadern, G.; Chohan, K. K.; Scrutton, N. S.; Sutcliffe, M. J.; Hillier, I. H.; Burton, N. A. J. Am. Chem. Soc. 2001, 123, 8604. (34) (a) Alhambra, C.; Sanchez, M. L.; Corchado, J.; Gao, J.; Truhlar, D. G. Chem. Phys. Lett. 2001, 347, 512. (b) Alhambra, C.; Sanchez, M. L.; Corchado, J.; Gao, J.; Truhlar, D. G. Chem. Phys. Lett. 2002 355 388. (35) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (36) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187. (37) Truhlar, D. G.; Isaacson, A. D.; Garrett, B. C. Generalized Transition State Theory. In Theory of Chemical Reaction Dynamics; Baer, M., Ed.; CRC Press: Boca Raton, FL, 1985; Vol. IV, p 65. (38) Gao, J.; Truhlar, D. G. Annu. ReV. Phys. Chem. 2002, 53, 467. (39) Truhlar, D. G. Variational Transition State Theory and Multidimensional Tunnelling for Simple and Complex Reactions in the Gas Phase, Solids, Liquids and Enzymes. In Isotope Effects in Chemistry and Biology; Limbach, H.-H., Kohen, A., Eds.; Taylor & Francis: Boca Raton, FL, 2006; p 579. (40) Tresadern, G.; Nunez, S.; Faulder, P. F.; Wang, H.; Hillier, I. H.; Burton, N. A. Faraday Discuss. 2003, 122, 223. (41) Nunez, S.; Tresadern, G.; Hillier, I. H.; Burton, N. A. Philos. Trans. R. Soc. B 2006, 361, 1387. (42) Truhlar, D. G.; Gao, J. L.; Garcia-Viloca, M.; Alhambra, C.; Corchado, J.; Sanchez, M. L.; Poulsen, T. D. Int. J. Quantum Chem. 2004, 100, 1136. (43) Garcia-Viloca, M.; Alhambra, C.; Corchado, J. C.; Sanchez, M. L.; Villa, J.; Gao, J.; Truhlar, D. G. CHARMMRATE, version 2.0; University of Minnesota: Minneapolis, MN, 2002. (44) Claeyssens, F.; Ranaghan, K. E.; Manby, F. R.; Harvey, J. N.; Mulholland, A. J. Chem. Commun. 2005, 5068. (45) Warshel, A.; Sharma, P. K.; Kato, M.; Xiang, Y.; Liu, H. B.; Olsson, M. H. M. Chem. ReV. 2006, 106, 3210. (46) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. Nucleic Acids Res. 2000, 28, 235. (47) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; JosephMcCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586. (48) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (49) Reuter, N.; Dejaegere, A.; Maigret, B.; Karplus, M. J. Phys. Chem. A 2000, 104, 1720. (50) Brooks, C. L., III.; Karplus, M. J. Chem. Phys. 1983, 79, 6312. (51) Jaguar, version 4.1; Schrodinger: Portland, OR, 2000.

J. Phys. Chem. B, Vol. 111, No. 11, 2007 3047 (52) Ponder, J. W. TINKER: Software Tools for Molecular Design, version 4.0; Saint Louis, MO, 2003. (53) Harvey, J. N. Faraday Discuss. 2004, 127, 165. (54) Mulholland, A. J.; Lyne, P. D.; Karplus, M. J. Am. Chem. Soc. 2000, 122, 534. (55) Alhambra, C.; Sanchez, M. L.; Corchado, J.; Gao, J. L.; Truhlar, D. G. Chem. Phys. Lett. 2001, 347, 512. (56) Alhambra, C.; Corchado, J.; Sanchez, M.; Gao, J.; Truhlar, D. J. Am. Chem. Soc. 2000, 122, 8197. (57) Proust-De Martin, F.; Dumas, R.; Field, M. J. J. Am. Chem. Soc. 2000, 122, 7688. (58) Bartels, C.; Karplus, M. J. Comput. Chem. 1997, 18, 1450. (59) Hermann, J. C.; Hensen, C.; Ridder, L.; Mulholland, A. J.; Holtje, H. D. J. Am. Chem. Soc. 2005, 127, 4454. (60) Tresadern, G.; Wang, H.; Faulder, P. F.; Burton, N. A.; Hillier, I. H. Mol. Phys. 2003, 101, 2775-2784. (61) Stewart, J. J. P. J. Comput. Phys. 1989, 10, 209. (62) Gonza´lez-Lafont, A.; Truong, T. N.; Truhlar, D. G. J. Phys. Chem. 1991, 95, 4618. (63) Hadfield, A. T.; Mulholland, A. J. Int. J. Quantum Chem. 1999, 73, 137. (64) Mulholland, A. J.; Richards, W. G. Proteins 1997, 27, 9. (65) Ridder, L.; Mulholland, A. J.; Rietjens, I.; Vervoort, J. J. Am. Chem. Soc. 2000, 122, 8728. (66) Ridder, L.; Mulholland, A. J.; Rietjens, I. M. C. M.; Vervoort, J. J. Mol. Graphics Modell. 1999, 17, 163. (67) Cramer, C. J. Essentials of Computational Chemistry: Theories and Models; John Wiley & Sons: Chichester, U. K., 2002. (68) Page, M.; McIver, J. W. J. Chem. Phys. 1988, 88, 922. (69) Chuang, Y. Y.; Corchado, J. C.; Truhlar, D. G. J. Phys. Chem. A 1999, 103, 1140. (70) Corchado, J. C.; Chuang, Y. Y.; Fast, P. L.; Villa, J.; Hu, W. P.; Liu, Y. P.; Lynch, G. C.; Nguyen, K. A.; Jackels, C. F.; Melissas, V. S.; Lynch, I.; Rossi, I.; Coitino, E. L.; Fernandez-Ramos, A.; Pu, J. Z.; Albu, T. V.; Steckler, R.; Garrett, B. C.; Isaacson, A. D.; Truhlar, D. G. POLYRATE, version 9.0; University of Minnesota: Minneapolis, MN, 2002. (71) Kreevoy, M. M.; Truhlar, D. G. Transition State Theory. In InVestigation of Rates and Mechanisms of Reactions; Bernasconi, D., Ed.; John Wiley & Sons: New York, 1986; Vol. 1; p 13. (72) Garcia-Viloca, M.; Alhambra, C.; Truhlar, D. G.; Gao, J. J. Chem. Phys. 2001, 9953. (73) Masgrau, L.; Gonzalez-Lafont, A.; Lluch, J. M. Chem. Phys. Lett. 2002, 353, 154. (74) Ranaghan, K. E.; Masgrau, L.; Scrutton, N. S.; Sutcliffe, M. J.; Mulholland, A. J. ChemPhysChem, submitted. (75) Zheng, Y. J.; Merz, K. M. J. Comput. Chem. 1992, 13, 1151. (76) Jurema, M. W.; Shields, G. C. J. Comput. Chem. 1993, 14, 89. (77) Pierdominici-Sottile, G.; Echave, J.; Palma, J. J. Phys. Chem. B 2006, 110, 11592. (78) Hur, S.; Bruice, T. C. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 1176. (79) Hur, S.; Bruice, T. C. J. Am. Chem. Soc. 2003, 125, 1472. (80) Hur, S.; Bruice, T. C. J. Am. Chem. Soc. 2003, 125, 5964. (81) Olsson, M. H.; Parson, W. W.; Warshel, A. Chem. ReV. 2006, 106, 1737. (82) Wigner, E. P. Z. Phys. Chem., Abt. B 1932, 19, 203. (83) Kreevoy, M. M.; Ostovic, D.; Truhlar, D. G.; Garrett, B. C. J. Phys. Chem. 1986, 90, 3766. (84) Pu, J. Z.; Ma, S. H.; Gao, J. L.; Truhlar, D. G. J. Phys. Chem. B 2005, 109, 8551. (85) Hothi, P.; Khadra, K. A.; Combe, J. P.; Leys, D.; Scrutton, N. S. FEBS J. 2005, 272, 5894. (86) Greatbanks, S. P.; Gready, J. E.; Limaye, A. C.; Rendell, A. P. J. Comput. Chem. 2000, 21, 788. (87) Garcia-Viloca, M.; Truhlar, D. G.; Gao, J. L. J. Mol. Biol. 2003, 327, 549. (88) Siebrand, W.; Smedarchina, Z. J. Phys. Chem. B 2004, 108, 4185. (89) Strajbl, M.; Shurki, A.; Kato, M.; Warshel, A. J. Am. Chem. Soc. 2003, 125, 10228. (90) Shurki, A.; Sˇ trajbl, M.; Villa´, J.; Warshel, A. J. Am. Chem. Soc. 2002, 124, 4097. (91) Ranaghan, K. E.; Mulholland, A. J. Chem. Commun. 2004, 1238. (92) Warshel, A.; Parson, W. W. Q. ReV. Biophys. 2001, 34, 563. (93) Olsson, M. H.; Mavri, J.; Warshel, A. Philos. Trans. R. Soc. B. 2006, 361, 1417. (94) Tejero, I.; Garcia-Viloca, M.; Gonzalez-Lafont, A.; Lluch, J. M.; York, D. M. J. Phys. Chem. B 2006, 110, 24708. (95) Masgrau, L.; Gonzalez-Lafont, A.; Lluch, J. M. J. Phys. Chem. A 1999, 103, 1044.

Tunneling and Classical Paths for Proton Transfer in an Enzyme

Recommend Documents