1758
J. Phys. Chem. B 2002, 106, 1758-1767
Enzymatic Reactions of Triosephosphate Isomerase: A Theoretical Calibration Study C. Lennartz,† A. Scha1 fer,† F. Terstegen,‡ and W. Thiel*,‡ Department of Scientific Computing/Computational Chemistry, ZDF/C, BASF Aktiengesellschaft, D-67056 Ludwigshafen, Germany, and Department of Theory, Max-Planck-Institut fu¨ r Kohlenforschung, Kaiser-Wilhelm-Platz 1, D-45470 Mu¨ lheim an der Ruhr, Germany ReceiVed: July 11, 2001; In Final Form: NoVember 6, 2001
Combined quantum mechanical (QM) and molecular mechanical (MM) calculations are reported for the triosephosphate isomerase-catalyzed conversion of dihydroxyacetone phosphate into glyceraldehyde 3-phosphate. The minima and transition states for the relevant proton-transfer reactions have been located on QM/ MM potential surfaces. The primary objective of this work is to study the sensitivity of optimized structures and relative energies toward variations in the QM/MM model, including the choice of the QM method, the size of the QM region, the size of the optimized MM region, and the treatment of the QM/MM boundary. The QM methods that have been applied in combination with the CHARMm force field range from semiempirical (AM1) to density functional (BP86, B3LYP) and ab initio (MP2) methods, the most extensive QM calculations involving 275 atoms and 2162 basis functions at the density functional level. Implications of the different choices of QM/MM options on the energy profile are discussed. From a mechanistic point of view, the present QM/MM results support a four-step proton-transfer pathway via an enediol, with involvement of neutral His95 acting as a proton donor, since the alternative direct intramolecular proton transfer in the enediolate is disfavored by the protein environment.
1. Introduction A detailed understanding of the mechanisms of enzyme catalysis as well as enzyme inhibition is crucial in modern life science. Computational approaches can contribute to such an understanding by providing detailed information on structures, energetics, and dynamics at the molecular level. Quantum mechanical (QM) methods are capable of describing all chemical processes, in principle, but normally they cannot be applied to whole enzymes due to the prohibitively large computational effort. On the other hand, standard molecular mechanical (MM) methods can be used for the simulation of large biochemical systems, but they cannot treat the breaking and making of bonds in enzyme reactions. Therefore, hybrid QM/MM methods offer a natural approach to model enzyme catalysis, since the active site and its protein environment can be described at appropriate QM and MM levels, respectively. The application of QM/MM methods requires a number of choices (e.g., regarding the actual QM and MM methods used, the QM/MM division, the representation of the QM/MM interaction terms, the treatment of the QM/MM boundary, and the approach to geometry optimization or conformational sampling). The numerical results will clearly depend on these choices. While studies on model systems have provided some insight and experience, it would still seem desirable to investigate the sensitivity of the results toward these factors in a real enzymatic reaction. We have chosen the triosephosphate isomerase- (TIM-) catalyzed conversion of dihydroxyacetone phosphate (DHAP) into glyceraldehyde 3-phosphate (GAP) for this purpose. The corresponding reactions have been the subject of several reviews1-3 and numerous detailed experimental4-13 and theoretical14-22 studies. Being a key step in the glycolytic pathway, † ‡
BASF Aktiengesellschaft. Max-Planck-Institut fu¨r Kohlenforschung.
Figure 1. Reaction mechanism for the catalytic interconversion of dihydroxyactone phosphate (DHAP, I) to D-glyceraldehyde 3-phosphate (GAP, V).
the conversion proceeds about 109 times faster in TIM than in aqueous solution under acetate ion catalysis.23,24 In fact, the reaction rate in TIM is limited by diffusion rather than by a chemical transformation,12 and TIM has therefore been characterized as a “perfect enzyme”.1 Isotope labeling experiments5 have established a free energy profile and an overall barrier of approximately 13 kcal/mol for the chemical transformations. As discussed in a recent theoretical paper,15 the precise chemical mechanism of the TIM-catalyzed reactions is not yet fully elucidated. The most widely accepted mechanism14,25 is depicted in Figure 1. The first step involves the abstraction of the pro-R proton of DHAP by the Glu165 side-chain carboxylic group to form an enediolate intermediate (II). A subsequent proton transfer from the pyrrolic nitrogen of the imidazole ring in the side chain of the His95 residue then yields an enediol
10.1021/jp012658k CCC: $22.00 © 2002 American Chemical Society Published on Web 01/11/2002
Enzymatic Reactions of Triosephosphate Isomerase intermediate (III). The imidazolate can be reprotonated by the terminal hydroxyl group of DHAP, thus forming another enediolate intermediate (IV) with the oxy and hydroxy positions exchanged in comparison with II. In the final fourth step, the enediolate (IV) is reprotonated by the glutamic acid side chain to form the product GAP (V). The first and last step of this mechanism are generally accepted, but there is some debate on the role of the enediol intermediate III, which has not yet been observed experimentally. As an alternative to the above two-step transformation II f III f IV (labeled as path A),15 a direct intramolecular proton transfer II f IV between the two enediolate intermediates has been proposed17-19 (path B).15 Another mechanistic variant (path C)15 considers a different two-step transformation, II f III f IV, where the glutamic acid from Glu165 replaces the imidazole from His95 in the role of proton donor and acceptor.10,11 In the present paper, we do not attempt to resolve these mechanistic issues. We shall focus on a QM/MM calibration study for the reactions of path A (see Figure 1) and present some results for path B, but we shall not consider path C. The reader is referred to other recent theoretical papers for further mechanistic discussions.15,16 In the comparisons with previous theoretical work, we shall concentrate mostly on related QM/ MM studies.14,15 This paper is structured as follows: After defining the chosen computational procedures (section 2) we discuss first the semiempirical QM/MM results (section 3) and thereafter density functional and ab initio QM/MM results (section 4). In both cases, small and large QM regions have been employed. The last part of the paper summarizes our conclusions (section 5). 2. Methods The X-ray structure of a TIM inhibitor (phosphoglycolohydroxamate) complex7 was taken from the PDB database (7TIM) and edited by QUANTA9826 to assign hydrogens and MSI-CHARMm27,28 25.2 force-field parameters and to convert the inhibitor into DHAP. All ionizable amino acid residues were assumed to be in their charged states and were not neutralized by counterions in the simulations. The prepared structure contains 8331 atoms, including 247 crystal water molecules, and consists of two TIM monomer subunits. The description as a dimer is important since the interconnecting loops extend to the vicinity of the active site, and one residue (Thr75) from the other subunit has indeed been found to influence the barriers through hydrogen-bonding interactions with His95.15 Two different QM regions were employed. The smaller one consists of 37 atoms including the substrate, the His95 side chain, the Glu165 side chain, and two hydrogen link atoms to saturate the dangling bond at each side chain (for I, DHAP, 3-methylimidazole, and propionate). The larger one contains 275 atoms and incorporates all residues and crystal water molecules within a sphere of about 10 Å centered on the substrate (see Supporting Information for further information and an explicit residue list) with 14 hydrogen link atoms to saturate the broken covalent bonds. The total charge of the QM region is -3 in both cases. Several QM approaches were used to describe the QM region, in particular the semiempirical AM1 method29,30 with standard parameters,30,31 the BP86 density functional,32,33 the B3LYP hybrid functional,33-35 and the correlated ab initio MP2 method.36 Resolution of identity (RI) techniques37 were applied for BP86 and MP238,39 to speed up the calculations. QM calculations for the small QM region were carried out with each
J. Phys. Chem. B, Vol. 106, No. 7, 2002 1759 of the above methods, while those for the large QM region were restricted to AM1 and BP86. The BP86 and B3LYP computations of the small QM region employed a split-valence basis set with polarization functions on all atoms (SVP),40 augmented with diffuse s and p functions for carbon, nitrogen, and oxygen. The exponents of these diffuse functions were derived from the s and p exponents of the SVP basis in an even-tempered manner. This basis will be denoted as SVP+ and comprises 458 basis functions. The RI-MP2 calculations of the small QM region used a valence triple-ζ basis set, polarized by two p and one d function on hydrogen and two d and one f function on all other atoms (TZVPP)41 and augmented by diffuse functions as described above. The resulting basis will be labeled as TZVPP+ and involves 961 basis functions. The BP86 treatment of the large QM region employed the SVP+ basis for the inner 37 atoms of the small QM region (see above) and an SV(P) basis set for the remaining QM atoms,40 yielding a total of 2162 basis functions. In the QM/MM calculations, the above QM approaches were combined with a CHARMm force field description27,28 of the MM region (parameters from MSI-CHARMm 25.2). A standard electronic embedding scheme was chosen (model B):42 the fixed MM atomic charges were included in the one-electron part of the QM Hamiltonian, thus leading to polarization of the QM region, and the QM/MM electrostatic interactions were evaluated from the QM electrostatic potential and the MM atomic charges. For the AM1 electrostatic potential, we adopted a special parametrization43 designed to mimic ab initio (RHF/6-31G*) electrostatic potentials,44 and used the corresponding oneparametric function with the published parameters (Table 2 of ref 43), whereas the electrostatic potentials for the other QM methods were computed analytically. No cutoffs were applied for the nonbonded interactions in the MM region. To estimate the errors introduced by truncating long-range nonbonded QM/MM interactions, we performed some calculations without any QM/MM cutoff and compared the results to calculations with a neutral group based cutoff of 15 Å. At a given optimized geometry, single-point energies with and without such a QM/MM cutoff may differ by up to 3 kcal/ mol. Structures are less sensitive: when starting from a structure optimized with the above cutoff, only a few cycles are needed for reoptimization without cutoffs, and the associated energy changes are generally below 0.35 kcal/mol. Therefore, all energies reported in this paper are calculated without QM/MM cutoffs, while a neutral group based cutoff of 15 Å was employed during geometry optimization. Choosing from the large number of methods available for treating the QM/MM boundary,42,45-53 we adopted hydrogen link atoms in combination with the charge shift model52 as a standard that can be applied for any QM approach. In this model, the link atom interacts with the MM atoms, but spurious polarization effects are minimized as follows: the charge of the MM atom being replaced by the link atom is shifted to the MM atoms directly bonded to this atom, and dipoles are added to these atoms to account for the polarity of the original bond. There are alternative link atom schemes where the link atom does not interact with any MM atom46 (denoted as L150 or QQ53) or where the link atom does not interact with the MM atoms of the adjacent neutral charge group (denoted as L250 or HQ53). At the semiempirical level, link atoms may be avoided by using adjusted connection atoms (ACA) that have been parametrized to mimic the structural and electronic properties of a methyl group.51 To assess the relative merits of these treatments, the
1760 J. Phys. Chem. B, Vol. 106, No. 7, 2002
Lennartz et al.
four approaches outlined above were compared in semiempirical QM/MM calculations with a small QM region. Due to the high conformational flexibility of an enzyme, computational studies of enzymatic reactions should include a sufficient sampling of configurational space, e.g., through molecular dynamics or Monte Carlo calculations. Such studies have been done at the semiempirical QM/MM level (for a recent example see ref 54), and special methods have been devised for QM/MM methods with an efficient sampling of the MM region only.22 In practice, however, the computationally less demanding approach of defining a reaction path and optimizing the relevant structures along this path is more common. For the present calibration study, this alternative approach was adopted. We selected residue-based spheres of different diameters centered on the substrate to be fully optimized, while the remaining atoms were fixed to the X-ray positions (see above). Our standard optimizations involved 576 atoms (14 Å sphere), and more extensive optimizations ranged from 971 atoms (18 Å sphere) to 2748 atoms (32 Å sphere). Explicit lists of the optimized residues are given in the Supporting Information. All QM/MM calculations were performed with the ChemShell software package55 using MNDO9756 for the semiempirical calculations and TURBOMOLE57 for the density functional and ab initio computations, respectively. The CHARMm force field was run through the parallel DL_POLY program,58 which is part of ChemShell. The HDLC optimizer59 included in ChemShell was used to carry out all geometry optimizations. Minima were found through the L-BFGS algorithm.60 Transition states (TS) were located by a microiterative procedure in which the reaction core is optimized by the P-RFO algorithm61 using an explicit numerical Hessian for eigenmode following, while the rest of the system is treated by the L-BFGS scheme. The default convergence criteria were applied,59 and it was verified in the case of transition states that the available Hessian (see above) had exactly one negative eigenvalue and an appropriate transition vector. Further details on the optimizer and the chosen optimization techniques can be found in the literature.59 3. Semiempirical Calculations In this and the following section, energy profiles will be presented for the TIM catalyzed reactions shown in Figure 1. Relative energies are given for the five minima I-V and the four intervening transition states, the zero value being assigned arbitrarily to the initial TIM-DHAP complex (I). These relative energies are derived from QM/MM total energies without zeropoint vibrational or other corrections. Details of the optimized geometries will only be discussed for the active site, which is depicted schematically for I-V in Figure 2. Deviations between optimized geometries are given as root-mean-square (rms) values, either for all optimized atoms or for certain subsets (e.g., crystal water, protein, or QM atoms), with the obvious convention that the fixed outer atoms (see above) are superimposed in such comparisons. For the qualitative characterization of structures through the number of hydrogen bonds, we employ the geometric criterion that a hydrogen bond X-H‚‚‚Y exists between two electronegative atoms (X, Y ) N or O) if the distance H‚‚‚Y is less than 2.5 Å and the angle X-H‚‚‚Y is greater than 90°. Finally, in this and the following section, the type of calculation will normally be indicated by the notation: QM-method(number of QM atoms)/CHARMm(diameter of optimized sphere). A. Link Atom Schemes. The results from AM1(37)/ CHARMm(14 Å) calculations with the charge shift option for the link atoms serve as reference. The corresponding energy
Figure 2. Important internal coordinates of the active site. IIIa and IIIb correspond to AM1 and DFT structures, respectively (see text). Dotted lines indicate the most important nonbonded interactions. Values of geometrical parameters are given in Tables 4 and 5.
profile was generated by following the individual proton-transfer reactions (see Figure 1) sequentially to make sure that the geometries change smoothly along the reaction path and exhibit a consistent hydrogen-bonding pattern. The optimized structures of the minima and transition states were then fully reoptimized at the same level of theory but with the L1, L2, and ACA schemes instead of the charge shift option. The resulting energies are given in Table 1 and shown in Figure 3. The structures obtained with these different procedures turn out to be very similar, with rms deviations of less than 0.1 Å for the QM region compared to the reference structure. The average absolute deviations from the reference energies are 0.1, 1.0, and 1.2 kcal/mol for L1, L2, and ACA, respectively, with corresponding maximum deviations of 0.2, 2.1, and 2.6 kcal/ mol. When assessing these results, one should keep in mind that the charge groups generated by QUANTA/MSI-CHARMm extend over whole residues (i.e., His95 and Glu165 are each treated as a single charge group). According to the definition of the L2 and ACA schemes, all charges in the groups next to
Enzymatic Reactions of Triosephosphate Isomerase
J. Phys. Chem. B, Vol. 106, No. 7, 2002 1761
TABLE 1: Relative Energiesa from AM1(37)/CHARMm (14 Å) for Different Treatments of the QM/MM Boundary structures treatment
I
TS
II
TS
III
charge shift L1 L1b L2 ACA
0.0 0.0 0.0 0.0 0.0
11.6 11.6 9.0 10.7 10.1
6.0 6.1 2.0 5.0 3.7
19.8 19.7 17.5 20.7 19.6
13.8 13.7 13.0 15.9 14.9
TS
IV
TS
V
20.4 4.4 14.2 6.2 20.3 4.5 14.2 6.1 18.0 -1.0 16.5 11.5 21.2 3.1 13.4 5.5 20.1 1.8 12.8 5.1
a Relative energies are given in kilocalories per mole. b Reference 14; AM1/CHARMM energies estimated from Figure 1 of that study.
Figure 4. Comparison of present AM1(37)/CHARMm and B3LYP(37)/CHARMm energy profiles with previous AM1/CHARMM14 and B3LYP/CHARMM15 results (see text).
Figure 3. AM1(37)/CHARMm(14 Å) energy profiles for different treatments of the QM/MM boundary (see text).
the link atoms (i.e., His95 and Glu165) are excluded from the electrostatic interactions with the QM and link atoms, while in the L1 case these charges are present for the QM atoms but invisible for the link atoms. In the charge shift model, the charges close to the link atoms are not deleted (as in L2 or ACA) but shifted to neighboring atoms where they interact with the QM and link atoms. It is obvious from this description that the QM/MM electrostatic interactions involving His95 and Glu165 are covered better by L1 and the charge shift option, and therefore it is not surprising that the corresponding results are almost identical and differ somewhat from the L2 and ACA results. Since previous studies have shown that the L1 scheme is less adequate for density functional and ab initio calculations,50 we shall adopt the charge shift model as a standard throughout this paper. It should be stressed, however, that the charge groups from QUANTA/MSI-CHARMm are untypically large and not representative for those from other protein force fields. They represent a worst-case scenario for the application of the L2 and ACA schemes, which will perform better for smaller charge groups where they neglect much fewer interactions. Given the fact that all four procedures applied yield very similar structures and that the deviations between the relative energies remain small on an absolute scale even for this unfavorable charge group definition, our main conclusion from these comparisons is that the present results are not overly sensitive to the exact treatment of the QM/MM boundary. B. Small QM Region. In an early QM/MM study14 the energy profile of the TIM-catalyzed reactions has been calculated using the same QM and MM methods (AM1/CHARMM), the same QM/MM division (small QM region), the same initial geometry (X-ray), an analogous electronic embedding (model B), and an analogous link atom scheme (L1) as in our present work. Minor technical differences concern the size of the system studied (1585 vs 8331 atoms), the diameter of the sphere being optimized (16 vs 14 Å), the type of optimization performed
(stochastic boundary conditions vs fixed outer atoms), the optimization techniques employed, the version of the force field used (academic CHARMM19 vs MSI-CHARMm 25.2), and the exact definition of the AM1 electrostatic potential.14,46 Furthermore, the energy used previously is given as the sum of the expectation value of the polarized QM Hamiltonian and the van der Waals QM/MM interactions,14 while we use the total QM/MM energy.42 The results from the previous AM1/CHARMM study14 are included in Table 1 and Figure 4. The corresponding energy profile is qualitatively similar to the present ones, with typical deviations of the order of 3 kcal/mol and maximum deviations around 5 kcal/mol for IV and V. Relative to the initial TIMDHAP complex (I), the current energies tend to be slightly higher than those obtained previously, except for the final transition state and the product (V), where they are somewhat lower (and thus more consistent with the experimental finding that the overall reaction is reversible and almost thermoneutral when small entropy differences are assumed). We believe that the observed deviations between these two independent AM1/ CHARMM studies are due to the technical differences outlined above, which prevent complete agreement. While such deviations of the order of 3 kcal/mol may be disconcerting, it is reassuring that both studies yield qualitatively the same energy profile, indicating that the AM1/CHARMM approach is fairly robust. To provide a more detailed understanding, we have analyzed the present energy profile in terms of individual energy contributions.42 The total QM/MM energy EQM/MM can be decomposed into the expectation value EQM of the polarized QM Hamiltonian (including the electrostatic QM/MM interactions) and the remainder EQM/MM - EQM (collecting the pure MM terms and the bonded and van der Waals QM/MM interactions). In Figure 5 the energy terms EQM and EQM/MM - EQM are plotted along the reaction coordinate. Obviously the overall shape of the profile is dominated by EQM, but there are also substantial contributions from EQM/MM - EQM (up to 15 kcal/mol). Further decomposition of EQM/MM - EQM (see above) shows that the electrostatic interactions within the MM region provide most of the corresponding variations (which are mainly due to subtle reorganizations of hydrogen bonds in the MM region in response to reorganizations within the active site). We have also identified a transition state on the AM1/ CHARMm potential surface for an intramolecular proton transfer from intermediate II directly to IV without an active role of His95 (path B). It contains a hydrogen atom that is
1762 J. Phys. Chem. B, Vol. 106, No. 7, 2002
Lennartz et al. TABLE 2: Relative Energiesa from QM(number of QM atoms)/CHARMm(14 Å) structures method AM1(37)b AM1(275)c AM1(275)b BP86/SVP+(37)d MP2/TZVPP+(37)e B3LYP/SVP+(37)e B3LYP/SVP+(37)b B3LYP/6-31+G(d,p)g BP86/SVP+(275)e BP86/SVP+(275)b
Figure 5. Decomposition of the AM1(37)/CHARMm(14 Å) energy profiles (see text).
Figure 6. AM1(number of QM atoms)/CHARMm(14 Å) energy profiles; // denotes single-point calculations (see text).
located almost symmetrically between the two enediolate oxygen atoms. The transition state lies 42 kcal/mol above the intermediate (II) and is thus highly unfavorable. It is well-known, however, that AM1 strongly overestimates the barrier to such proton transfers between oxygen atoms, giving a value of 34 kcal/mol for the simplest gas-phase enediolate model system (HO-CHdCH-O-). It is therefore more meaningful to point out that the protein environment increases the barrier for the direct proton transfer by 8 kcal/mol in this case, which is close to the shift of 5.5 kcal/mol reported at the B3LYP/6-31+G(d,p) level.15 Hence, path B is disfavored due to the presence of His95, according to both AM1/CHARMm and B3LYP/CHARMM.15 C. Large QM Region. As discussed in section 2, the large QM region comprises 275 atoms and includes all residues within a sphere of about 10 Å around the substrate. Application of the energy decomposition scheme described above shows that the AM1(275)/CHARMm energy profile for the TIM-catalyzed reactions is now almost entirely determined by EQM since the differential contributions from EQM/MM - EQM remain below 3 kcal/mol. When such a large QM region is used, the barriers should thus depend directly only on the chosen QM method, although the MM method can exert an indirect influence through the geometries. Single-point AM1(275)/CHARMm energies were calculated at the optimized AM1(37)/CHARMm reference geometries. In addition, starting from these reference geometries, full reoptimizations were carried out at the AM1(275)/CHARMm level. The corresponding relative energies are shown in Table 2 and Figure 6, while selected geometrical data are listed in Table 4 with the atom numbering from Figure 2.
I 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
TS
II
TS
III
TS
IV
11.6 6.0 19.8 13.8 20.4 4.4 11.2 9.1 22.7 19.9 28.1 10.3 20.1 12.4 21.9 14.6 21.9 11.9 6.8 3.4 16.2 12.3 16.3 4.7 13.9 5.8 16.2 9.3 17.0 7.5 12.3 6.8 21.0 16.0 21.2 8.8 9.2 6.4 19.1 15.8 (17)f 8.5 11.3 6.2 14.6 4.5 13.5 9.7 20.9 17.1 21.7 11.5 9.2 15.4 8.7
TS
V
14.2 6.2 18.7 9.1 14.2 2.7 9.0 1.9 16.2 5.0 15.3 4.8 10.7 3.9 7.4 2.3 18.1 10.7 2.5
a Relative energies are given in kilocalories per mole. b Optimized geometries. c Single-point energies at AM1(37)/CHARMm optimized geometries. d Optimized intermediates; transition states located only approximately from reaction path calculations (see text). e Single-point energies at BP86(37)/CHARMm optimized geometries. f Approximate value (see text). g Reference 15; B3LYP/CHARMM energies taken from Figure 1 of that study.
The energy profile from the single-point AM1(275)/CHARMm calculations is qualitatively similar to the AM1(37)/CHARMm reference profile, but the relative energies are higher by ca. 3 kcal/mol around II and V and even by ca. 6 kcal/mol between III and IV. Reoptimization at the AM1(275)/CHARMm level preserves all the stationary points but shifts their relative energies in a nonuniform manner (see, e.g., the increase for II and the transition state leading to II and the decrease for V). These discrepancies indicate that CHARMm and AM1 differ in the treatment of hydrogen bonds or/and the protein backbone, since the transition from the small to the large QM region effectively changes the description of a number of water molecules and residues from CHARMm to AM1. It should be noted in this connection that AM1 shortcomings in these areas have been documented.62 For a more detailed comparison of the AM1(37)/CHARMm and AM1(275)/CHARMm structures, Table 4 reports selected geometrical parameters in the active site (comprising DHAP, His95, and Glu165; see Figure 2). At both levels, the same type of hydrogen-bonding network is found for each intermediate. However, there are some notable differences in the optimized bond lengths, which exceed 0.2 Å in the case of III [see R(O4H3), R(C1-H1), and R(N-H2)] and are on the order of 0.1 Å or less otherwise. Similar trends are found for the bond angles, which again show the largest deviations for III. Nevertheless, in an overall assessment, the active-site geometries are rather similar at both levels. Table 4 also lists the rms deviations in the positions of the protein residues and the water molecules within the optimized 14 Å sphere. The former range from 0.30 Å (I) to 0.38 Å (V) and the latter from 0.54 Å (I) to 0.74 Å (III). These deviations are mostly due to different hydrogen-bonding patterns in the region around the active site. Most of the water molecules solvate the highly polar phosphate group of the substrate. Counting the number of hydrogen bonds formed between the water molecules and the substrate only, we find 14/15 hydrogen bonds for I/III in the AM1(37)/CHARMm structures and 21/24 hydrogen bonds in the AM1(275)/CHARMm structures, respectively. Thus, the shell of solvating water molecules around the phosphate group is more compact in the case of the large QM system. On the other hand, counting the hydrogen bonds in the protein subsystem of I/III yields 36/37 H-bonds for AM1(37)/CHARMm and 27/29 hydrogen bonds for AM1(275)/ CHARMm. These different distributions of hydrogen bonds for
Enzymatic Reactions of Triosephosphate Isomerase
J. Phys. Chem. B, Vol. 106, No. 7, 2002 1763
TABLE 3: Relative Energiesa from AM1(37)/CHARMm (size of optimized MM region) MM sphere
I
TS
II
TS
14 Å 18 Åb 18 Åc 32 Å
0.0 0.0 0.0 0.0
11.6 7.3 13.0
6.0 3.3 8.3 13.2
19.8 14.4 20.1
structures III TS 13.8 7.7 12.0 27.9
20.4 11.2 16.5
IV
TS
V
4.4 3.0 8.2 6.5
14.2 6.5 17.8
6.2 0.2 5.9 15.6
a Relative energies are given in kilocalories per mole. b Geometries (I) obtained by direct optimization (see text). c Geometries (II) obtained after revision of the hydrogen-bond network (see text).
easy to arrive at local minima that differ in the hydrogenbonding network or with regard to other conformational changes. This can then cause significant energy variations that are unrelated to the issues under investigation and therefore completely misleading. It is thus essential to compare only optimized structures along a consistent reaction path, i.e., in our case AM1(37)/CHARMm(18 Å)-II. It may be advisable to restrict the optimization to the smallest sphere that is still acceptable (e.g., 14 Å) to facilitate the control over such factors. Alternatively, molecular dynamics simulations of sufficient length could be used to sample conformational space, but this would be more costly. The optimizations within a 32 Å sphere represent another computational experiment. The resulting relative energies of the intermediates are mostly unrealistic (see Table 3). This is as expected, since almost a complete TIM monomer is optimized when a 32 Å sphere is used, so that part of the enzyme is exposed to vacuum boundary conditions. This leads to large geometry changes for the surface exposed to the vacuum and to concomitant energy changes. Therefore, these results are only included for documentation and will not be discussed further. 4. Density Functional and Ab Initio Calculations
Figure 7. AM1(37)/CHARMm(size of relaxed MM region) energy profiles (see text).
different QM/MM divisions influence the active site generally only in an indirect manner, but there is also a more direct influence by a thiol group of the Cys126 residue: while this group does not exhibit any hydrogen bonding to the active site in the case of the small QM region, it causes changes in the positions in the Glu165 residues of I and II by hydrogen bonding with the carboxylate group in the case of the large QM region. D. Size of the Optimized MM Region. As described in section 2, our standard geometry optimizations involve 576 atoms within a sphere of about 14 Å around the substrate. At the AM1(37)/CHARMm level, we have investigated the effects of increasing the diameters of the relaxed spheres to 18 and 32 Å, respectively. These optimizations were started from the optimized AM1(37)/CHARMm(14 Å) reference geometries. The results are shown in Table 3 and Figure 7 under the label AM1(37)/CHARMm(18 Å)-I. The corresponding relative energies differ appreciably from the 14 Å reference energies: they are generally lower, by up to 9 kcal/mol for the transition state III f IV. These changes are much larger than expected for a mere extension of the relaxed part of the MM region. Closer inspection of the results reveals that the hydrogenbonding patterns between the QM region and the adjacent enzyme residues are not always the same in the 14 and 18 Å structures: there is an additional hydrogen bond between the phosphate group of the substrate and the water residue HOH788 that is present in only some of the 18 Å structures. To correct for this inconsistency, HOH788 was manually reoriented in the respective 18 Å structures, which were then reoptimized. This generated a more consistent reaction profile labeled as AM1(37)/ CHARMm(18 Å)-II, which deviates much less from the reference profile, with maximum deviations of up to 4 kcal/mol in relative energies and up to 2 kcal/mol in individual barriers. The preceding discussion highlights a typical difficulty encountered in automatic geometry optimizations of large systems. Due to the huge number of degrees of freedom, it is
A. Small QM Region. Density functional theory (DFT) and correlated ab initio methods have been applied to check how the QM/CHARMm results change when AM1 is replaced by higher-level QM methods. The chosen computational procedures have been described in section 2, and geometry optimizations have been restricted to a 14 Å sphere. We first discuss the results obtained with a BP86(37)/SVP+ description of the QM region. Starting from the optimized AM1(37)/CHARMm(14 Å) reference geometries, the reaction intermediates I-V were reoptimized successfully. The resulting structures were qualitatively similar to the AM1/CHARMm reference structures, except for the enediol (III), which transformed to the enediolate (IV) during the optimization. After testing several conformations of the two hydroxyl groups in III, we located a chelate minimum structure which is sketched in Figure 2 as IIIb. Hence, two hydrogen bonds are needed to stabilize the negatively charged deprotonated nitrogen of the imidazolate in III with BP86(37)/ SVP+ (in contrast to AM1, where only one such bond is present). Turning to the transition states, the automatic search by eigenmode following did not succeed in any case when starting from the AM1(37)/CHARMm reference structures, despite the correct shape of the starting BP86/CHARMm Hessian of the reaction core (one negative eigenvalue, correct mode). The transition modes were found to be extremely soft with very small eigenvalues, which changed into modes with positive eigenvalues after a few optimization steps, and the search then finally converged to an adjacent minimum (intermediate). To get approximate barrier heights with BP86(37)/SVP+, we performed reaction path calculations. Keeping the distance between the corresponding donor and acceptor atoms fixed at the reactant value, the donor-hydrogen distance was increased in steps of 0.1 Å until the barrier was crossed. The maximum energy values were taken as approximate transition-state energies that are given in Table 2 and in Figure 8 together with the relative energies of the intermediates. Using the corresponding geometries as starting points for the eigenmode following search, the transition states leading from I to II and from IV to V could be located, but the problems outlined above persisted for the other two transition states. In the two successful refinements, the barriers were reduced somewhat (by 1.2 and 3.5 kcal/mol,
1764 J. Phys. Chem. B, Vol. 106, No. 7, 2002
Lennartz et al.
TABLE 4: Comparison of Optimized Structures from AM1(37)/CHARMm and AM1(275)/CHARMma structure parameterb
Ic
IIc
IIIc
IVc
Vc
R(O1-H2) R(O1-H3) R(O2-H3) R(O3-H1) R(O4-H2) R(O4-H3) R(C1-H1) R(C2-H1) R(N-H2) R(N-H3) θ(H2-O1-C1-C2) θ(H1-C1-C2-O2) θ(H3-O2-C2-C1)
0.981, 0.977 2.205, 2.283 1.990, 2.061 2.143, 2.050 2.076, 1.974
0.971, 0.971 2.261, 2.233 1.921, 1.936 0.988, 0.989 2.426, 2.332
0.977, 0.989 2.663, 2.826 0.973, 0.973 0.977, 0.976
2.061, 2.012 3.174, 3.219 0.970, 0.969 0.984, 0.985
2.129, 2.123 3.414, 3.345 0.975, 0.975 2.034, 1.926 2.077, 2.146
2.099, 2.159
2.357, 2.185 2.605, 2.339
2.137, 2.146
1.127, 1.133
2.473, 2.361 1.008, 1.004
1.143, 1.151 0.997, 0.993
94.8, 96.1 134.9, 122.2
1.007, 1.004 82.1, 83.4 90.2, 95.8
2.121, 1.883 2.401, 2.544 12.4, 2.1 76.8, 97.8 45, 62.8
17.7, 4.8 83.8, 93.3 95.9, 98.1
5.7, 21.2
rms(H2O) rms(protein)
0.54 0.30
0.66 0.33
0.74 0.34
0.65 0.33
0.66 0.38
134.2, 125.2
a Important internal coordinates of the active site and rms deviations of the water and protein environment are given. Atom numbering is from Figure 2. b Distances and rms deviations are given in angstroms; angles are given in degrees. c AM1(275) is the first number listed; AM1(37) is the second number listed.
TABLE 5: Comparison of Optimized Structures from BP86(37)/CHARMm and BP86(275)/CHARMma structure parameterb
Ic
IIc
IIIc
IVc
Vc
R(O1-H2) R(O1-H3) R(O2-H3) R(O3-H1) R(O4-H2) R(O4-H3) R(C1-H1) R(C2-H1) R(N-H2) R(N-H3) θ(H2-O1-C1-C2) θ(H1-C1-C2-O2) θ(H3-O2-C2-C1)
1.022, 1.052 2.421, 2.312 1.814, 1.916 2.438, 2.277 1.690, 1.537
0.998, 0.991 2.550, 2.587 1.630, 1.645 1.067, 1.055 2.085, 2.332
1.032, 1.022 2.483, 2.466 1.005, 1.023 1.001, 0.996
1.764, 1.723 3.189, 3.221 0.994, 0.998 1.049, 1.051
2.520, 2.221 3.431, 3.352 1.036, 1.061 1.881, 1.718 1.540, 1.496
1.721, 1.802
2.968, 2.971 2.255, 2.345
1.783, 1.748
1.112, 1.116
1.841, 1.867 1.057, 1.063
1.145, 1.166 1.050, 1.040
68.9, 85.8 149.6, 136.5
1.074, 1.069 82.5, 90.6 106.8, 100.5
1.624, 1.706 1.788, 1.707 0.1, 4.3 103.2, 99.5 12.2, 11.6
8.1, 3.7 107.7, 94.6 102.2, 104.6
3.5, 0.4
rms(H2O) rms(protein)
0.73 0.26
0.74 0.26
0.75 0.26
0.74 0.26
0.72 0.29
127.1, 119.1
a Important internal coordinates of the active site and rms deviations of the water and protein environment are given. Atom numbering is from Figure 2. b Distances and rms deviations are given in angstroms; angles are given in degrees. c BP86(275)is the first number listed; BP86(37) is the second number listed.
Figure 8. BP86,B3LYP,MP2/CHARMm(14 Å) energy profiles; // denotes single-point calculations (see text).
respectively). For the sake of consistency, we shall only discuss the reaction path values in the following. Compared with the AM1(37)/CHARMm results, the BP86(37)/ CHARMm relative energies are generally slightly lower (see Table 2), except for the enediolate IV, which is destabilized by 0.5 kcal/mol, whereas the other intermediates are stabilized by
1-4 kcal/mol relative to I. The barriers for the four proton transfers are also consistently lower with BP86, typically by about 4 kcal/mol. The rate-determining steps are the transformations II f III f IV with overall barriers of ca. 16 kcal/mol (BP86) and 20 kcal/mol (AM1). When the optimized AM1(37)/CHARMm and BP86(37)/ CHARMm structures are compared, the rms deviations of the optimized protein region are of the order of 0.1 Å and below. This similarity is not surprising since this region is described essentially at the same level (CHARMm) in both cases. Small deviations arise from differences in the QM electrostatic potentials used to evaluate the QM/MM Coulomb interactions. These deviations are largest for the most polar interactions between QM and MM regions, i.e., the interactions between water molecules and the charged phosphate group of the substrate where the corresponding rms values range from 0.16 up to 0.30 Å. In the various optimized intermediates, there are typically 15-17 hydrogen bonds crossing the boundary between QM and MM region for BP86/CHARMm, compared with 1718 for AM1/CHARMm. In most cases the atoms involved in these hydrogen bonds are identical. However, AM1/CHARMm predicts interactions between the positively charged Lys12 side chain and the phosphate ester oxygen atom in some structures,
Enzymatic Reactions of Triosephosphate Isomerase which are not present for BP86/CHARMm. Concerning the internal structure of the active site, it has already been pointed out that the optimized AM1/CHARMm and BP86/CHARMm geometries are qualitatively similar in the sense that they show the same hydrogen bonded network (except for the enediol III, see above). Despite this qualitative agreement, substantial quantitative discrepancies can occur for the nonbonded distances (up to 0.5 Å, see Tables 4 and 5). For further energetic comparisons, single-point QM/CHARMm calculations were performed at the BP86/CHARMm optimized geometries describing the QM region by the B3LYP/SVP+ and MP2/TZVPP+ levels of theory (see Table 2). The resulting profiles are compared with the BP86/SVP+ results in Figure 8. It is obvious that the B3LYP-based relative energies are consistently higher than the BP86-based values, by 3-6 kcal/ mol, whereas the energy profile for MP2 resembles the B3LYP profile for the transformations I f TS f II and IV f TS f V and the BP86 profile for the remaining transformations. It is also remarkable that the relative energies from AM1/CHARMm are actually quite close to those from the B3LYP/CHARMm and MP2/CHARMm single-point calculations, with mean absolute deviations of 1.5 and 2.5 kcal/mol, respectively (see Table 2). The potential well for the enediol intermediate III appears to be somewhat deeper with B3LYP and MP2 than with BP86 (see Table 2). Since the two transition states enclosing III could not be located precisely at the BP86/CHARMm level, we reoptimized the available BP86/CHARMm structures at the B3LYP/CHARMm level to see whether these missing transition states could be optimized properly. This was indeed possible for the reaction II f III but not for III f IV, where the search produced low-gradient structures about 1-2 kcal/mol above III with one negative Hessian eigenvalue, which, however, failed to satisfy the standard convergence criteria. For the sake of completeness, we also reoptimized all other BP86/CHARMm structures at the B3LYP/CHARMm level to generate the corresponding energy profile (see Table 2 and Figure 8). Compared with the single-point B3LYP/CHARMm results, the reoptimization lowers all energies relative to I slightly (by 0.23.6 kcal/mol, on average less than 2 kcal/mol). Our best estimated barriers for the proton transfers from the enediol III to II and IV are very small (3.3 and 1-2 kcal/mol, respectively). Given the limited accuracy of our current computational approach, we can therefore not be certain whether the enediol III actually exists as a shallow minimum on the potential surface. In fact, another recent B3LYP/CHARMM study,15 which used a different basis set [6-31+G(d,p)], a smaller QM region during optimization (substrate and His95), other optimization techniques, and a different link atom scheme, did not locate such a minimum for III (see Figure 4 for a visual comparison of the B3LYP-based profiles). Our best B3LYP/CHARMm estimate for the transition structure III f IV resembles the AM1/CHARMm geometry of the enediol with regard to the hydrogen-bonding pattern (see Figure 2, IIIa). According to B3LYP/CHARMm, the proton transfer IIIb f IV involves the breaking of the hydrogen bond between H3 and O4 and the formation of a new hydrogen bond between H3 and N, whereas according to AM1/CHARMm this rearrangement of hydrogen bonds occurs already in the preceding proton transfer II f IIIa (see Figure 2). More generally speaking, the transformation II f IV with involvement of His95 requires the reorganization of at least three hydrogen bonds, so that we may expect considerable energetic compensation from the concomitant breaking and making of these bonds. It may
J. Phys. Chem. B, Vol. 106, No. 7, 2002 1765
Figure 9. BP86/SVP+(number of QM atoms)/CHARMm(14 Å) energy profiles (see text).
therefore be most appropriate to stress that the potential energy surface is rather flat around the enediol (III) and that the proton transfers involving His95 thus proceed almost on an isoenergetic plateau. As in the AM1/CHARMm case, we have also identified a transition state on the BP86/CHARMm potential surface for an intramolecular proton transfer from intermediate II directly to IV without an active role of His95 (path B). Its geometry is qualitatively similar to the one obtained from AM1/CHARMm concerning the location of the migrating proton. The transition state lies 18.5 kcal/mol above the initial TIM-DHAP complex (I), and the barrier for the direct intramolecular proton transfer from II to IV is 15.1 kcal/mol, which is close to the recently reported B3LYP/6-31+G(d,p)/CHARMM value of 14.0 kcal/ mol.15 A comparison with the simplest gas-phase model system (HO-CHdCH-O-, BP86/SVP+ barrier of 5.0 kcal/mol) shows again that the protein environment disfavors the direct proton transfer considerably (here by about 10 kcal/mol). Compared with the two-step process involving His95 (path A), the BP86/CHARMm barrier for the one-step process (path B) is only about 2 kcal/mol higher (see Table 2). The preference for path A is thus less pronounced than at the B3LYP/631+G(d,p)/CHARMM level,15 where the relevant barriers differ by 5.6 kcal/mol. B. Large QM Region. To check the internal consistency of the QM/MM coupling model also at the density functional level, calculations with a large QM region were carried out in analogy to the semiempirical QM/MM computations described in section 3C. The BP86/SVP+ approach was chosen for efficiency reasons (see section 2). Single-point BP86(275)/CHARMm energies were calculated at the available BP86(37)/CHARMm geometries. Furthermore, starting from these geometries, full reoptimizations were performed for all intermediates at the BP86(275)/CHARMm level. The corresponding relative energies are given in Table 2 and Figure 9, while selected structural data are listed in Table 5 with the atom numbering from Figure 2. The relative energies from the single-point BP86(275)/ CHARMm calculations are consistently higher than the BP86(37)/ CHARMm values, by 5-9 kcal/mol. This implies that the largest change occurs for I, whose relative energy has arbitrarily been assigned to zero: shifting the energy scale such that II becomes the origin, the variation in the relative energies due to the extension of the QM region remains below 3 kcal/mol for all species except I, indicating that most of these changes are rather uniform. Compared with the single-point values, reoptimization of the intermediates lowers the relative energies
1766 J. Phys. Chem. B, Vol. 106, No. 7, 2002 slightly for II-IV (0.5-2.8 kcal/mol) and appreciably for V (8.2 kcal/mol). The optimized BP86(37)/CHARMm and BP86(275)/CHARMm geometries of the active site are qualitatively similar in the sense that they show the same hydrogen-bonding network, although substantial deviations in the nonbonded distances may occur [up to 0.3 Å for R(O1-H2) in V]. The rms deviations in the positions of the protein residues and the water molecules within the optimized 14 Å sphere are rather uniform. The former range from 0.26 to 0.29 Å and are thus somewhat smaller than in the AM1/CHARMm case, while the latter lie between 0.72 and 0.76 Å and thus tend to be larger than for AM1/CHARMm (see Table 4). As discussed in section 3C, the major differences between the AM1/CHARMm results for small and large QM regions come from a more compact solvation shell of water molecules around the phosphate group. This is not the case for BP86/ CHARMm since the number of the corresponding hydrogen bonds does not change much: for I/III, there are 14/15 and 14/14 hydrogen bonds between water molecules and substrate with a small and large BP86 region, respectively. However, the positions of the water molecules do change upon reoptimization, both around the phosphate group and elsewhere. In addition, the number of hydrogen bonds in the protein subsystem for I/III decreases from 39/39 to 33/34 when going from BP86(37) to BP86(275). In an ideal QM/MM combination, the potential surfaces generated by the chosen QM and MM methods would match perfectly, and increasing the size of the QM region at the expense of the MM region would not affect the results. Both for the AM1/CHARMm and BP86/CHARMm approaches, this is clearly not the case in our present application. The energy profiles for the TIM-catalyzed reactions remain qualitatively similar, but there are nonnegligible changes when the QM region is extended from 37 to 275 atoms, both in single-point energies and after reoptimization. Likewise, active-site geometries remain qualitatively similar, but there are structural changes in the surroundings, particularly in the positions of water molecules. A more quantitative assessment shows that typical and maximum energy changes as well as structural rms deviations are of comparable magnitudes for AM1/CHARMm and BP86/ CHARMm. However, in an overall assessment, the changes upon extension of the QM region tend to be more uniform in the case of BP86/CHARMm. 5. Discussion and Conclusions The present QM/MM validation study on TIM-catalyzed proton transfers has focused on the question of how much the optimized structures and relative energies are influenced by variations in the QM/MM model, in particular with regard to the choice of the QM method, the size of the QM region, the size of the optimized MM region, and the treatment of the QM/ MM boundary. In a qualitative sense, the computed energy profiles for the four-step conversion of DHAP into GAP (path A) and the active-site geometries of the corresponding intermediates and transition states are not too sensitive to such variations since important qualitative features remain unchanged. The quantitative results are, of course, affected by changes in the QM/MM options. To be more specific, the choice for the description of the QM/MM boundary in semiempirical AM1/CHARMm calculations has not been critical in the present case since different link atom and connection atom schemes lead to similar results (almost identical geometries and typical variations in relative energies of 1 kcal/mol or less). Likewise, an extension of the
Lennartz et al. size of the optimized MM region in semiempirical AM1/ CHARMm computations causes only minor changes in relative energies and barriers (typically around 2 kcal/mol), provided that the different optimized structures refer to a consistent reaction path with analogous conformations and analogous hydrogen-bonding networks. Much care is required to satisfy this condition in large-scale optimizations. Concerning the QM/MM division, both a small QM region (37 atoms) and a large QM region (275 atoms) were used with semiempirical and density functional QM methods. The decomposition of the QM/MM energy expression for the AM1(37)/ CHARMm calculations shows that the shape of the reaction profile mostly reflects the QM energy (which includes the electrostatic QM/MM interactions in the chosen electronic embedding); however, the other MM-type terms (especially the nonbonded Coulomb terms) provide substantial differential contributions (up to 15 kcal/mol in our case). Their influence diminishes in the AM1(275)/CHARMm calculations with the large QM region (differential effects up to 3 kcal/mol only). QM/MM studies with such a large QM region and QMdominated energetics are desirable but normally impractical. When a smaller QM region is used, the consistency between the QM and MM potentials becomes important, especially for the immediate environment of the active site. Comparisons of the present QM/MM results for small and large QM regions indicate that neither AM1 nor BP86 matches CHARMm perfectly in the description of the protein and water environment around the active site. Even though the basic features of the energy profile are preserved, the transition to a larger QM region causes shifts in the relative energies that are commonly around 3 kcal/mol but may also be higher (up to 9 kcal/mol). Moreover, there are large rms deviations for the water coordinates that have been traced to a more compact hydration of the phosphate group in the case of AM1(275)/CHARMm, while the rms deviations for the protein coordinates are smaller, especially for BP86(275)/CHARMm. The results with small and large QM regions generally tend to be more uniform for BP86 than for AM1. In our opinion, the AM1/CHARMm combination seems to be better balanced with a small QM region, and we see no clear evidence that the BP86/CHARMm results are improved significantly upon extension of the QM region. Concerning the choice of the QM method, the expected accuracy and the computational effort increase in the sequence AM1 < BP86 < B3LYP < MP2. It is reassuring that all QM/ CHARMm approaches applied presently lead to qualitatively similar reaction profiles. In all cases the conversion between the enediolates II and IV contains the rate-limiting step of path A. For all types of QM methods, the highest QM/CHARMm barrier on this path is computed to lie between 15 and 22 kcal/ mol (with one exception, see Table 2), which may be compared with the experimental free energy barrier of 13 kcal/mol.5 As pointed out before, the active-site geometries are also qualitatively similar at different QM/CHARMm levels, with the exception of the enediol III, where a chelate structure with two hydrogen bonds at the deprotonated nitrogen atom of the imidazolate ring is found by BP86/CHARMm and B3LYP/ CHARMm, whereas AM1/CHARMm predicts a minimum with a single hydrogen bond (see Figure 2). From a mechanistic point of view, it is debatable whether the enediol III exists at all on the potential surface.15 All QM/ MM methods applied presently (see above) have located a proper enediol minimum that satisfies our standard convergence criteria for energy minimization.59 On the other hand, this minimum is only very shallow, and it has been difficult to find
Enzymatic Reactions of Triosephosphate Isomerase the transition states that separate III from II and IV. Our best estimates (B3LYP/CHARMm) are around 3 and 1-2 kcal/mol, respectively, so that we prefer to stress that the surface around III is quite flat and therefore hard to characterize unambiguously. Another mechanistic issue is the feasibility of a direct intramolecular proton transfer from II to IV without any involvement of His95 (path B).17-19 Here, our current AM1/ CHARMm and BP86/CHARMm results support the notion15 that this path requires more activation since it is disfavored by the protein environment (by about 8-10 kcal/mol). According to our calculations, the alternative two-step conversion with involvement of neutral His95 acting as a proton donor (path A) is the more likely mechanism. Acknowledgment. This work was supported by the European commission (ESPRIT Project 25047, QUASI). We thank Professor M. Karplus for helpful discussions. Supporting Information Available: Lists of residues included in the large QM region and in the different optimized MM regions. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Knowles, J. R.; Albery, W. J. Acc. Chem. Res. 1977, 10, 105. (2) Knowles, J. R. Nature 1991, 350, 121. (3) Gerlt, J. A. Bioorganic Chemistry: Peptides and Proteins; Hecht, S. M., Ed.; Oxford University Press: New York, 1998; pp 279-311 and references therein. (4) Herlihy, J. M.; Maister, S. G.; Albery, W. J.; Knowles, J. R. Biochemistry 1976, 15, 5607. (5) Albery, W. J.; Knowles, J. R. Biochemistry 1976, 15, 5627. (6) Albery, W. J.; Knowles, J. R. Biochemistry 1976, 15, 5588. (7) Davenport, R. C.; Bash, P. A.; Seaton, B. A.; Karplus, M.; Petsko, G. A.; Ringe, D. Biochemistry 1991, 30, 5821. (8) Lolis, E.; Alber, T.; Davenport, R. C.; Rose, D.; Hartman, F. C.; Petsko, G. A. Biochemistry 1990, 29, 6609. (9) Lodi, P. J.; Knowles, J. R. Biochemistry 1991, 30, 6948. (10) Harris, T. K.; Abeygunawardana, C.; Mildvan, A. S. Biochemistry 1997, 36, 14661. (11) Harris, T. K.; Cole, R. N.; Comer, F. I.; Mildvan, A. S. Biochemistry 1998, 37, 16828. (12) Nickbarg, E. B.; Knowles, J. R. Biochemistry 1988, 27, 5939. (13) Nickbarg, E. B.; Davenport, R. C.; Petsko, G. A.; Knowles, J. R. Biochemistry 1988, 27, 5948. (14) Bash, P. A.; Field, M. J.; Davenport, R. C.; Petsko, G. A.; Ringe, D.; Karplus, M. Biochemistry 1991, 30, 5826. (15) Cui, Q.; Karplus, M. J. Am. Chem. Soc. 2001, 123, 2284. (16) Aquist, J.; Fothergill, M. J. Biol. Chem. 1996, 271, 10010. (17) Alagona, G.; Desmeules, P.; Ghio, C.; Kollman, P. A. J. Am. Chem. Soc. 1984, 106, 3623. (18) Alagona, G.; Ghio, C.; Kollman, P. A. J. Mol. Biol. 1986, 191, 23. (19) Alagona, G.; Ghio, C.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 9855. (20) Pera¨kyla¨, M.; Pakkanen, T. A. Proteins: Struct., Funct., Genet. 1996, 25, 225.
J. Phys. Chem. B, Vol. 106, No. 7, 2002 1767 (21) Pera¨kyla¨, M. J. Chem. Soc., Perkin Trans. 1997, 2185. (22) Zhang, Y.; Liu, H.; Yang, W. J. Chem. Phys. 2000, 112, 3483. (23) Hall, A.; Knowles, J. R. Biochemistry 1975, 14, 4348. (24) Richard, J. P. J. Am. Chem. Soc. 1984, 106, 4926. (25) Stryer, L. Biochemistry, 4th ed.; Freeman: New York, 1995. (26) Program QUANTA98, Molecular Simulations Inc., 9685 Scranton Rd., San Diego, CA 92121-3752. (27) CHARMm Program: CHARMm 25.2, Molecular Simulations Inc., 9685 Scranton Rd., San Diego, CA 92121-3752. (28) CHARMM: Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (29) Dewar, M. J. S.; Thiel, W. J. Am. Chem. Soc. 1977, 99, 4899. (30) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. A.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902. (31) Dewar, M. J. S.; Jie, C. J. Mol. Struct. 1989, 187, 1. (32) Perdew, J. P. Phys. ReV. B 1986, 33, 8822. (33) Becke, A. D. Phys. ReV. A 1988, 36, 3098. (34) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (35) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (36) Moller, C.; Plesset, M. S. Phys. ReV. 1934, 46, 618. (37) Eichkorn, K.; Treutler, O.; Oehm, H.; Ha¨ser, M.; Ahlrichs, R. Chem. Phys. Lett. 1995, 240, 283. (38) Weigend, F.; Ha¨ser, M. Theor. Chem. Acc. 1997, 97, 331. (39) Weigend, F.; Ha¨ser, M.; Patzelt, H.; Ahlrichs, R. Chem. Phys. Lett. 1998, 294, 143. (40) Scha¨fer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 2571. (41) Scha¨fer, A.; Huber, C.; Ahlrichs, R. J. Chem. Phys. 1994, 100, 5829. (42) Bakowies, D.; Thiel, W. J. Phys. Chem. 1996, 100, 10580. (43) Bakowies, D.; Thiel, W. J. Comput. Chem. 1996, 17, 87. (44) Ford, G. P.; Wang, B. J. J. Comput. Chem. 1993, 14, 1101. (45) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1986, 7, 718. (46) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11, 700. (47) Thery, V.; Rinaldi, D.; Rivail, J.-L.; Maigret, B.; Ferenczy, G. J. Comput. Chem. 1994, 15, 269. (48) Gao, J.; Amara, P.; Alhambra, C.; Field, M. J. J. Phys. Chem. A 1998, 102, 4714. (49) Zhang, Y.; Lee, T.-S.; Yang, W. J. Chem. Phys. 1999, 110, 46. (50) Antes, I.; Thiel, W. Hybrid Quantum Mechanical and Molecular Mechanical Methods; Gao, J., Ed.; ACS Symposium Series 712; American Chemical Society: Washington, DC, 1998; pp 50-65. (51) Antes, I.; Thiel, W. J. Phys. Chem. 1999, 103, 9290. (52) Vries, A. H.; Sherwood, P.; Collins, S. J.; Rigby, A. M.; Rigutto, M.; Kramer, G. J. J. Phys. Chem. 1999, 103, 6133. (53) Reuter, N.; Dejaegere, A.; Maigret, B.; Karplus, M. J. Phys. Chem. B 2000, 104, 1720. (54) Billeter, S. R.; Hanser, C. F. W.; Mordasini, T. Z.; Scholten, M.; Thiel, W.; van Gunsteren, W. F. Phys. Chem. Chem. Phys. 2001, 3, 688. (55) ChemShell is a modular QM/MM program developed in the European QUASI project under the coordination of P. Sherwood; see http:// www.cse.clrc.ac.uk/Activity/ChemShell. (56) Thiel, W. Program MNDO97, version 5.0, University of Zu¨rich, 1998. (57) Ahlrichs, R.; Ba¨r, M.; Ha¨ser, M.; Horn, H.; Ko¨lmel, C. Chem. Phys. Lett. 1989, 162, 165. (58) Smith, W.; Forester, T. J. Mol. Graph. 1996, 14, 136. (59) Billeter, S. R.; Turner, A. J.; Thiel, W. Phys. Chem. Chem. Phys. 2000, 2, 2177. (60) Liu, D. C.; Nocedal, J. Math. Prog. 1989, 45, 503. (61) Banerjee, A.; Adams, N.; Simons, J.; Shepard, R. J. Phys. Chem. 1985, 89, 52. (62) Mo¨hle, K.; Hofmann, H.-J.; Thiel, W. J. Comput. Chem. 2001, 22, 509.
