Phosphate-Dependent Enzyme l-Serine Dehydratase

Sep 24, 2008 - The catalytic mechanism of a pyridoxal 5′-phosphate-dependent enzyme, l-serine dehydratase, has been investigated using ab initio ...
0 downloads 0 Views 1MB Size
Download PDF
J. Phys. Chem. B 2008, 112, 13091–13100

13091

A Quantum Mechanical/Molecular Mechanical Study on the Catalysis of the Pyridoxal 5′-Phosphate-Dependent Enzyme L-Serine Dehydratase Zheng Zhao and Haiyan Liu* School of Life Sciences, and Hefei National Laboratory for Physical Sciences at the Microscale, UniVersity of Science and Technology of China (USTC), Hefei, Anhui, 230027, China ReceiVed: March 14, 2008; ReVised Manuscript ReceiVed: June 19, 2008

The catalytic mechanism of a pyridoxal 5′-phosphate-dependent enzyme, L-serine dehydratase, has been investigated using ab initio quantum mechanical/molecular mechanical (QM/MM) methods. New insights into the chemical steps have been obtained, including the chemical role of the substrate carboxyl group in the Schiff base formation step and a proton-relaying mechanism involving the phosphate of the cofactor in the β-hydroxyl-leaving step. The latter step is of no barrier and follows sequentially after the elimination of the R-proton, leading to a single but sequential R, β-elimination step. The rate-limiting transition state is specifically stabilized by the enzyme environment. At this transition state, charges are localized on the substrate carboxyl group, as well as on the amino group of Lys41. Specific interactions of the enzyme environment with these groups are able to lower the activation barrier significantly. One major difficulty associated with studies of complicated enzymatic reactions using ab initio QM/MM models is the appropriate choices of reaction coordinates. In this study, we have made use of efficient semiempirical models and pathway optimization techniques to overcome this difficulty. Introduction The importance of pyridoxal 5′-phosphate (PLP)-dependent enzymes in living organisms has been long known. Renewed interest in these enzymes has been triggered by recent functional genomics and other studies using improved experimental techniques, and the discovery of their involvements in diseases.1,2 PLP-dependent enzymes can be classified into at least three families, R, β, and γ, based on their amino acid sequences.3 One particularly interesting member of the β-family is the L-serine dehydratase (SDH, EC 4.3.1.17). SDH catalyzes the R, β-elimination of L-serine to yield pyruvate and ammonia (Figure 1).2,4-7 Possible physiological functions of SDH in different organisms have been extensively studied. For example, in microorganisms, two isozymes have been found subjecting to allosteric regulations by AMP as an activator and isoleucine as an inhibitor.8-10 In mammals, SDH has been found with relatively abundant concentrations in rat liver. The enzyme activity can be induced remarkably under various abnormal dietary conditions, including the consumption of high-protein diets or starvation.11,12 Thus it has been suggested that SDH may play important roles in gluconeogenesis.11,12 Recent results that SDH activity changes in response to acidosis suggested that SDH participated in gluconeogenesis also in kidney besides in liver.13,14 In contrast to rat, SDH was found to have low abundance and low activity in human liver, suggesting that this enzyme may have different functions in human.15,16 It was also suggested that SDH was related to tumors in mammals, as its activity has been found to be absent in human colon carcinoma and rat sarcoma.1,17 Besides its physiological functions, the catalytic mechanism of SDH is also an issue of interest. It has been suggested that the R, β-elimination of the hydroxyl of L-serine goes in two steps (Figure 1).2,4,6,7 The first step involves the elimination of * To whom correspondence should be addressed. E-mail: hyliu@ ustc.edu.cn.

Figure 1. The R, β-elimination reaction of L-serine.

the hydroxyl group to produce aminoacrylate, and the aminoacrylate is deaminated in the second step. Only the first step requires enzyme catalysis with the second being a nonenzymatic hydrolysis reaction. The suggested mechanism for the enzymatic step assigns the cofactor PLP an important role. In SDH as well as in other PLPdependent enzymes, PLP is linked to a lysine in the form of a Schiff base.1,2,10 However, it has been suggested that in different enzymes, PLP may be in different protonation states at the N1 atom of the pyridine ring, which may result in two different classes of mechanisms for the R, β-elimination.10 In one class, the pyridine is protonated and can form a quinine-like resonance state after the R hydrogen of the substrate has been extracted. In the other class, the pyridine is not protonated and no quininelike resonance state may occur. Recently, high-resolution crystal structures of SDH from different species have been determined.1,10,16,18 Particularly, the crystal structure of rat liver SDH (rSDH) indicates a hydrogen bond network at the active site. This network is consistent with a model in which the enzyme-bound PLP is not protonated.19-22 Thus the mechanism of rSDH may belong to the second class which lacks the resonance mechanism to stabilize the intermediate. A possible catalytic mechanism for rSDH has been described by Yamada and co-workers (Figure 2) based on crystal structures.10 In this mechanism, the amino group of the substrate serine becomes a strong nucleophile after transferring a proton to the phosphate group of PLP (see step 1 in Figure 2, which is the transfer of proton H1 from N3 of serine to O2 of phosphate). The deprotonated amino group attacks the PLP-Lys41 Schiff

10.1021/jp802262m CCC: $40.75  2008 American Chemical Society Published on Web 09/24/2008

13092 J. Phys. Chem. B, Vol. 112, No. 41, 2008

Figure 2. Schematic representations of possible steps and intermediate species during the catalysis of serine dehydratase. Various atoms have been labeled using numerical subscripts to facilitate discussions. K41 represents the side chain of Lys41 in the active site.

base to form the PLP-Ser Schiff base (see step 2 in Figure 2). Subsequently, the R hydrogen (H3 in Figure 2) and the β-hydroxyl group (O5H in Figure 2) of the substrate serine are eliminated, resulting in a water molecule and a PLP-aminoacrylate intermediate (see step 3 in Figure 2). Finally, a reversed nucleophilic attack by Lys41 on this intermediate takes place, the PLP-Lys41 Schiff base recovered, and the aminoacrylate released (see step 4 in Figure 2). In the above mechanism, the phosphate group participates as a general base. The direct binding environments of the phosphate in PLP-rSDH and in some other PLP binding proteins have been compared using sequence and structure comparisons.10 It has been found that the compositions of charged and hydrogen bond-forming residues in these binding environments differ qualitatively. The unusual lacking of charged residues interacting with the phosphate in PLP-rSDH has been suggested to facilitate the phosphate to play the assumed catalytic role. Besides the role of the phosphate, the proposed mechanism assigns no extra catalytic group for the second proton transfer from the amino group of the substrate serine to the amino group of Lys41, implicating a direct transfer either in concert or in sequential with the nucleophilic attack on the PLP-Lys41 Schiff base. The mechanism also does not specify how the eliminations of the R-hydrogen and the hydroxyl group from the substrate serine take place. Thus the current understanding of the catalytic mechanism is in a sense ambiguous and incomplete. The main issues include whether there are other groups participating in the chemical mechanism, what is the rate limiting step, and how the corresponding transition state is stabilized by the enzyme environment. Given the available experimental data for the rSDH system, combined ab initio QM/MM approach is well suited for addressing these questions.23 Such an approach has produced reliable insights into mechanistic issues for a number of enzyme systems.24-31 The multistep mechanism of rSDH is a typical case of complicated enzymatic reactions involving various bond-forming and bond-breaking steps in unclear orders or even involving uncertain chemical factors. One major technical difficulty in theoretical studies of such processes is the appropriate choice of reaction coordinates (RCs). This is especially true when ab initio QM/MM models are applied and extensive explorations of different RCs in higher-dimensional spaces are too expensive.

Zhao and Liu To avoid choosing one-dimensional RCs artificially or by trials-and-errors in our ab initio QM/MM modeling, we employed the computationally more efficient semiempirical QM/ MM model32-35 to explore possible reaction paths. This is done in two different ways for different reaction steps. When one or two RCs can be unambiguously identified (for example, steps 1 and 2 in Figure 2), simple restrained optimizations are used to scan the one- or two-dimensional potential energy surfaces (PESs). For the more complicated reaction step (step 3 in Figure 2), we use an nudged elastic band (NEB) method36,37 adapted for enzyme systems38 to optimize the reaction path without any presumed RC. The PESs or MEP energy curves computed using the MNDO-type semiempirical QM/MM model and using single point ab initio QM/MM calculations are of similar shapes, suggesting that appropriate one-dimensional RCs can be derived based on the semiempirical results. Subsequently, these onedimensional RCs are employed to drive restrained optimizations using the ab initio QM/MM model and the iterative optimization approach described by Zhang et al.39 and to obtain the final minimum energy paths (MEP) for the different reaction steps. We note that the computational efficiency of the semiempirical QM/MM model has the potential to allow statistical mechanical sampling of the reaction system. However, for the reaction steps studied here, the MNDO-type semiempirical QM/ MM models, although the produced PESs or energy curves of correct general shapes, give barrier heights differing significantly from those obtained using the ab initio model (see Results and Discussion). Thus the semiempirical QM/MM results have been used only to guide the selections of RCs. In addition, given the focus of the current study, which is on the mechanistic insights generated by the final ab initio QM/MM results and not the relative strengths and weaknesses of different QM/MM models, no extensive testing of other semiempirical models on the particular reactions has been attempted. Despite this, we note that QM/MM approaches based on improved semiempirical models including orthogonalization corrections or d-orbital,40-42 or system-specific parameters43,44 are being developed. For some well-calibrated systems, such improved semiempirical approaches potentially provide useful alternative route, by allowing for extensive statistical mechanics sampling. Materials and Methods 1. The Enzymatic System. The initial structure of rSDH in complex with PLP and L-serine has been constructed from the crystallographic structure of rSDH in complex with PLP-(Omethylserine) aldimine (PDB ID: 1PWH).10 The substrate serine has been modeled by replacing the methyl group of Omethylserine with a hydrogen atom and adjusting the O-H bond length to 1.0 Å. The amino group of the substrate forms a Schiff base with PLP. The resulting model corresponds to IM2 in Figure 2. The phosphate group of PLP is surrounded mainly by neutral amide groups and is in the monoprotonated form. The N1 atom of the pyridine ring of PLP is not protonated, and the O1 is in the oxyanion form to accept two hydrogen bonds from its surroundings as indicated by the crystal structure.10 Except for Lys41, all other titratable side chains of the protein are charged to satisfy salt-bridge and hydrogen-bond requirements. The amino group of Lys41 is uncharged in IM2. To model an explicit solvent environment, the solute complex was solvated in a rectangular box filled with pre-equilibrated waters. The box was centered at the C1 atom (see labels in Figure 2), and had side lengths of 60.8, 74.8, and 83.9 Å, respectively. Water molecules closer than 2.3 Å or further than 8.0 Å from any solute non-hydrogen atom were then removed, resulting in

Pyridoxal 5′-Phosphate-Dependent Enzyme L-Serine Dehydratase a shell of solvent enclosing the entire complex and containing 2725 water molecules. 2. QM/MM Models. The QM part consisted of 51 atoms including PLP, the substrate serine, and the side chain of Lys41. The rest formed the MM part. The QM part has been described either by the semiempirical model AM145 or by ab initio density functional models at various levels of theory, and the MM part by the AMBER force field and the TIP3P model for water.46,47 The QM/MM nonbonded interactions have been treated as described in detail by Field and Karplus et al.,39,48 and the bonded boundaries treated by the pseudobond approach.49-51 We would like to mention that QM/MM calculations based on the PM352 and MNDO53 models have also been performed. These two models performed similarly as AM1 for our particular system, and the corresponding results are not presented. 3. Equilibration and Optimization of the Initial Model. First, the system was optimized until the root-mean-square energy gradient was below 1.0 kcal mol-1 Å-1 with a position restraining force constant of 100.0 kcal mol-1 Å-2 for the solute atoms. Then the solvent shell was equilibrated with 200 ps restrained molecular dynamics simulation at 300 K, during which the positions of all solute heavy atoms presented in the crystal structure had been restrained by harmonic potentials with force constants of 100.0 kcal mol-1 Å-2. The resulting solvent configurations were optimized again until the root-mean-square gradient is below 1.0 kcal mol-1 Å-1. Using the resulting configuration, a set of active atoms were defined which included all solute atoms and solvent molecules within 18 Å from the center of geometry of all QM atoms. In all further calculations, the positions of the shell of atoms outside this active sphere had been fixed. The active sphere was further equilibrated sequentially by four 20 ps molecular dynamics with gradually reducing force constants for the position restraints (60.0, 40.0, 20.0, and 0.0 kcal mol-1 Å-2, respectively). Finally, a 200 ps molecular dynamics simulation was performed without restraining any atoms in the active sphere to sample an equilibrated initial configuration for the mechanistic investigations. This equilibrium simulation has been performed using a pure MM model. The structure of the active site has been relatively well maintained during the 200 ps simulation. In the starting structures, the carboxylate of the substrate serine forms hydrogen bonds with the backbone amide groups of Ala65 and Ala68, with donor-acceptor heavy atom distances of 2.89 and 3.15 Å, respectively. During the simulation, the averaged distances between the O7, O6 atoms on this carboxylate with the respective donor hydrogen atoms on Ala65 and Ala68 are 1.98 and 2.14 Å, respectively. The PLP phosphate is held in position by five hydrogen bonds with the backbone groups of Gly168 to Leu172. These hydrogen bonds have also been relatively well maintained in the simulation with averaged acceptor oxygen-donor hydrogen distances between 1.72 to 2.45 Å. The hydrogen bonds between O1 on PLP and Asn67 has also been kept, with an averaged acceptor oxygen-donor hydrogen distance of 1.94 Å. These interactions between the reaction center and the surrounding groups kept the relative positioning of the reaction groups. The averaged distance between N2 and H3 is 2.94 Å, between C1 and N2 is 2.70 Å, and between O5 and H1 is 2.50 Å. We selected a snapshot from the 200 ps trajectory to start the QM/MM calculations based on three interatomic distances (the ones between C1 of PLP and N2 of Lys41, between C2 of Ser and N2 of Lys41, and between O5 of Ser and O2 of the phosphate group). The root-mean-square deviations of these distances in different recorded conformations from the respective

J. Phys. Chem. B, Vol. 112, No. 41, 2008 13093 distances in the crystal structure were computed and compared. The last recorded conformation in the trajectory with a deviation below 0.1 Å (the snapshot at 134.2 ps) was chosen. These distances had been chosen because they are directly involved in the bond-forming or breaking process or in the proposed proton transfer mechanisms. Their respective values in the crystal structure are consistent with the chemical mechanism and provide suiting references for choosing a starting structure in lack of more specific information. The chosen conformation was then optimized using the semiempirical QM/MM model until the root-mean-square energy gradient was below 0.1 kcal mol-1 Å-1. 4. Calculating the MEPs. Starting from the optimized IM2 configuration, we computed the MEPs for the first half of the reaction (from the reactant to IM2) in the backward direction and for the second half (from IM2 to IM3) in the forward direction. Step 4 which is the reversed nucleophilic attack by Lys41 has not been considered here because it can be viewed as a reversing of step 2. To avoid confusion, we will describe all investigated steps in the forward order. Determining the MEP for Step 1. Starting from the reactant, the charged amino group of the substrate serine was activated by transferring a proton to the PLP phosphate, resulting in IM1. For this step (step 1 in Figure 2), we computed a onedimensional PES along the RC ) d[N3-H1]-d[O2-H1], where d[X-Y] stands for the distance between atoms X and Y. We directly used the B3LYP/6-31G*/MM model and the reaction coordinate driving (RCD) and iterative optimization approach as described by Zhang et al.39 Subsequently, single point calculations using diffusion basis functions (B3LYP/6-311+G*/ MM) were performed on structures along the MEP. A Two-Dimensional PES for Step 2. From IM1 to IM2 (step 2 in Figure 2), the amino group of Lys41 was substituted by that of the substrate serine. In addition, a proton is transferred between the respective amino groups. We first considered a mechanism involving no other chemical group to assist this proton transfer. This leads to the definition of two RCs, one for the amino group substitution, RC1 ) d[C1-N2]-d[C1-N3] and the other for the direct proton transfer, RC2 ) d[N3-H2]-d[N2-H2]. As it would be computationally too demanding to determine a fully relaxed two-dimensional PES using the ab initio QM/MM model, we first determined the two-dimensional PES using the AM1/MM model. Then the PES was recomputed applying the B3LYP/6-31G*/MM model to the AM1/MM geometries. Comparing the semiempirical and ab initio PES indicated the qualitative correctness of the AM1/MM surface. The Two-Dimensional PES for an Indirect Proton Transfer Mechanism in Step 2. The previous PES indicated that the above mechanism for step 2 results in a high activation barrier, mainly associated with the direct proton transfer between the two amino groups (see Results and Discussion). The proton transfer takes place between two local minima on the PES, both having RC1 values close to zero (that is, the two amino groups are of approximately the same distances from C1). We denote the minimum with the proton on the serine as IM1′ and the other minimum with the proton on Lys41 as IM2′. Close examination of the structures of IM1′ and IM2′ suggested that the carboxyl group of the substrate serine is at a position to chemically assist the proton transfer. Using IM1′ and IM2′ as ending states, we computed another two-dimensional PES using RC1′ ) d[N2-H2]-d[O6-H2] and RC2′ ) d[N3-H2]-d[O6-H2]. Again complete relaxations of the geometries for each point on the PES were performed at the AM1/MM level.

13094 J. Phys. Chem. B, Vol. 112, No. 41, 2008

Zhao and Liu

Figure 4. The overall potential energy surface along the final minimum energy path. Circles are the B3LYP/6-31G*/MM results, and squares are the B3LYP/6-311 + G*/MM//B3LYP/6-31G*/MM results. The reaction coordinates are d[N3-H1]-d[O2-H1] from R to IM1, d[C1-N2]-d[C1-N3] from IM1 to IM1′, (d[N2-H2]-d[O6-H2]) + (d[N3-H2]-d[O6-H2]) from IM1′ to IM2′, d[C1-N2]-d[C1-N3] from IM2′ to IM2, and d[O2-H1] + d[C3-O5]-d[O5-H1] + d[C2H3]-d[N2-H3] + d[N2-H4]-d[O2-H4] from IM2 to IM3. Here d[X-Y] refers to the distance between two atoms X and Y. Figure 3. Different proton-relaying schemes for the R, β-elimination step (step 3 in Figure 2).

Subsequent B3LYP/6-31G*/MM calculations on the AM1/MM geometries confirmed that the AM1/MM surface is qualitatively correct. Determining the Final MEP for Step 2. The two-dimensional PESs described above suggested that a complete one-dimensional PES for step 2 should consist of three substeps. The first is from IM1 to IM1′, for which the appropriate RC can be d[C1-N2]-d[C1-N3]. The second is from IM1′ to IM2′. For this step, we chose (d[N2-H2]-d[O6-H2]) + (d[N3-H2] d[O6-H2]) as the RC. The last sub step is from IM2′ to IM2, for which d[C1-N2]-d[C1-N3] was again chosen as the RC. Starting from IM2 and going backward until IM1, the RCD and iterative QM/MM optimization approach was employed to obtain a continuous one-dimensional PES and a fully relaxed MEP at the B3LYP/6-31G*/MM level. As for step 1, B3LYP/ 6-311+G*/MM single point calculations were performed on the final MEPs. Exploring Different Proton-Relaying Mechanisms for Step 3 Using NEB. On the basis of the structure of IM2, two possible mechanisms for the abstraction of the R proton and for the protonation of the leaving hydroxyl group in step 3 can be proposed (Figure 3). These mechanisms involve different proton transfer processes that are difficult to explore using PESs of lower dimensions. Thus for initial explorations, we used a NEB method adapted for enzymatic systems to determine the MEPs at the AM1/MM level. For each candidate mechanism shown in Figure 3, a model of the ending IM3 state was obtained using successive optimizations starting from IM2, first with distancerestraining potentials forcing the broken bonds to be longer than 2.0 Å and the newly forming bonds to be within 1.5 Å (or 1.0 Å for bonds involving hydrogen) and then without any restraints. Then an initial path represented by a number of intermediate structures was obtained by linear interpolation in the Cartesian space between the ending IM2 and IM3 structures. This initial path was minimized to obtain the MEP using the adapted NEB approach described by Xie et al.38 The different candidate mechanisms resulted in different MEPs because they defined different correspondences between protons in IM2 and IM3. We then performed single point B3LYP/6-31G*/MM calculations for the intermediate structures on the AM1/MM MEPs.

The second shown in Figure 3 was found to have a substantially lower activation barrier than scheme 1. We note that in our NEB calculations, we did not make any assumption about the order of the R-deprotonation and β-hydroxyl leaving processes. For scheme 2, we also performed alternative NEB calculations that treated these two processes as separated steps. The results are the same NEB path and potential energy surface as obtained by not treating them separately. Determining the Final MEP for Step 3. On the basis of the NEB results, a complex one-dimensional RC was chosen based on scheme 2. The RC was defined as RC ) d[O2-H1] + d[C3-O5]-d[O5-H1] + d[C2-H3]-d[N2-H3] + d[N2-H4]-d[O2-H4]. The RCD approach was again employed to obtain the one-dimensional PES along this RC using the B3LYP/631G*/MM model and iterative QM/MM optimizations. Again, B3LYP/6-311+G*/MM single point calculations were performed on the final MEP. 5. Other Computational Details. The semiempirical QM/ MM calculations have been performed using a program54 interfaced with TINKER.55 The ab initio QM/MM calculations have been performed with Gaussian56 and TINKER.39,55 In the molecular dynamics simulations, all bonds have been constrained using SHAKE57 and the integration time step was 2 fs. In semiempirical QM/MM calculations, coupling to external temperature bath has been achieved using the weak coupling approach58 with a relaxation time of 0.1 ps. Unless specified otherwise, convergence criterion for geometry optimizations has been that the root-mean-square energy gradient falls below 0.1 kcal mol-1 Å-1. Results and Discussions We will first describe the chemical mechanisms and the final MEPs for all the three steps. Then, we will present the mechanism-exploring results, including the two-dimensional PESs for step 2 and the NEB results for step 3, to justify the RCs used to obtain the final MEPs. Finally, we will analyze how transition states have been stabilized by the enzyme environments, based on the final MEPs. 1. A Detailed Chemical Mechanism of the Enzymatic Step in rSDH Catalysis. Figure 4 shows the PESs along the final MEPs for all intermediate steps from the reactant to IM3. We

Pyridoxal 5′-Phosphate-Dependent Enzyme L-Serine Dehydratase will focus on the results at the B3LYP/6-311+G*/MM//B3LYP/ 6-31G*/MM level. Step 1 is the proton transfer from the substrate amino group to the phosphate. There is a single transition state (TS1) with a computed barrier height of 7.8 kcal mol-1. The same barrier is 6.4 kcal mol-1 at the B3LYP/6-31G*/MM level. At TS1, d[N3-H1] ) 1.10 Å and d[O2-H1] ) 1.50 Å. The product of step 1, IM1, is computed to be of similar energy as the reactant. Thus once bound to the enzyme-PLP complex, the substrate serine may easily transfer one of its amino protons to the phosphate in a nonrate limiting step. Step 2 (from IM1 to IM2) consists of three sub steps. From IM1, the amino group of the substrate serine attacks the Schiff Base without losing its second proton to produce the germinal diamine intermediate IM1′. The corresponding transition state TS2-1 is associated with a low energy barrier (2.5 kcal mol-1). At this transition state, the C1-N3 bond has not yet formed (d[C1-N3] ) 2.04 Å) and the C1-N2 bond keeps its doublebond character (d[C1-N2] ) 1.35 Å). At IM1′, these distances changed further to produce a tetrahedral configuration at C1 (d[C1-N2] ) 1.45 Å and d[C1-N3] ) 1.53 Å). The second substep from IM1′ to IM2′ involves the proton transfer from the amino group of the substrate serine to the amino group of Lys41, indirectly via the carboxyl group of the substrate. The transition state TS2-2 is associated with a moderate barrier (9.4 kcal mol-1). At this transition state, the C-N bond lengths have not changed much relative to IM1′. The transferred proton is still far from the final acceptor nitrogen (d[N2-H2] ) 2.49 Å). The proton-donor nitrogen bond has almost broken (d[N3-H2] ) 1.37 Å). A new bond between the proton and one of the carboxyl oxygen atoms of the substrate has been formed with a slightly lengthened bond length (d[O6-H2] ) 1.17 Å). After TS2-2 and in IM2′, the proton has been transferred completely to Lys41. IM2′ is of almost the same energy as IM1′. In the last substep from IM2′ to IM2, the protonated amino group of Lys41 leaves the tetrahedral intermediate, resulting in the serine Schiff base. As TS2-1 for the first substep, the corresponding transition state TS2-3 for this substep is also associated with a low activation barrier (3.0 kcal mol-1). At this transition state, the C1-N2 bond has been broken (d[C1-N2] ) 1.94 Å) and the C1-N3 bond length is of double-bond character (d[C1-N3] ) 1.34 Å). Interestingly, starting from TS2-3 until IM2, the energy of the system decreased significantly (by more than 10 kcal mol-1). Step 3 from IM2 to IM3 is a single step combining R-deprotonation and β-hydroxyl-leaving in a sequential order. The transition state takes place after the R-deprotonation and before the β-hydroxyl-leaving, and is associated with the highest activation barrier. This step is determined to be rate-limiting, with a computed barrier of 19.4 kcal mol-1. Enlarging the basis set also has the largest effects on the computed barrier heights for this step: the computed barrier is 31.0 kcal mol-1 at the B3LYP/6-31G*/MM level. The high energy of TS3 relative to IM2 and its strong basis set size-dependence is probably caused by the significant concentration of negative charges on the substrate, which already bears a unit negative charge on the serine carboxyl. Changes in key interatomic distances from IM2 through TS3 to IM3 (see Table 1 and Figure 5) suggest that the R-deprotonation and the β-hydroxyl-leaving take place sequentially. At TS3, the R hydrogen has already been abstracted by Lys41, as indicated by the changes in the C2-H3 and N2-H3 bond lengths from IM2 to TS3 (Table 1). Although the C3-O5 bond has been

J. Phys. Chem. B, Vol. 112, No. 41, 2008 13095

Figure 5. Changes of various interatomic distances along the minimum energy path for step 3. The position of TS3 is indicated by an arrow.

significantly lengthened (from 1.47 to 1.57 Å), the leaving hydroxyl group is not yet protonated. After TS3 the transfer of H1 and H4 and the leaving of the hydroxyl group take place in a concerted, barrier-less way along the MEP (Figure 5). We note that the sequential order of R-deprotonation and β-hydroxyl-leaving is consistent with experimental studies on the base-catalyzed β-elimination reactions in aqueous solution.59 On the MEP computed here, the carbanion species has been identified as a transition state instead of a meta-stable intermediate. Although uncertainties in our theoretical approach do not allow us to completely exclude the existence of any intermediate in step 3, our results suggest that if such an intermediate exists, it should be only marginally stable and be very close to the true TS3 in both structure and energy. 2. The Three Substeps and an Indirect Proton Transfer Mechanism in Step 2. In the above chemical mechanism, we have suggested for step 2 an indirect proton transfer path in which the substrate carboxyl group chemically assisted the proton transfer. This has been based on the two-dimensional PESs constructed as described in the method section. To construct the first PES, the two RCs are (RC1, RC2) ) (d[C1-N3]-d[C1-N2], d[N2-H2]-d[N3-H2]). The results are shown as Figure 6. Figure 6a shows the AM1/MM results, and 6b shows the B3LYP/6-31G*/MM//AM1/MM results. The positions of IM1 and IM2 are (1.3, 1.8) (in Å, the same below) and (-1.6, -1.8), respectively. Along the MEP in Figure 6a, there are two intermediates, IM1′ at (0.2, 1.5) and IM2′ at (-0.2, -1.5). There are also three transition states, TS2-1 at (0.8, 1.8) separating IM1 and IM1′, TS2-3 at (-0.8, -1.5) separating IM2′ and IM2, and TS2-4 at (0.0, 0.0) separating IM1′ and IM2′. The notation TS2-4 has been used to differentiate this transition state from the transition state TS2-3 which also separates IM1′ and IM2′ but along the indirect proton transfer pathway (see below). Obviously, step 2 can be represented as a sequential process during which N3 approaches C1 first (nucleophilic attack), then the proton H2 is transferred, and finally N2 leaves C1. The B3LYP/6-31G*/MM//AM1/MM PES shown as Figure 6b are qualitatively the same as the AM1/ MM one, suggesting that we can safely base our definitions of one-dimensional RCs for the final ab initio QM/MM optimizations on these approximate PESs. Figure 6,a,b indicates that direct proton transfer from N3 to N2 is energetically unfeasible. The activation barrier associated with TS2-4 is 39.0 kcal mol-1 at the AM1/MM level and 34.0 kcal mol-1 on the single point ab initio QM/MM PES. The MEP in Figure 6a that has been obtained using d[N2-H2]-d[N3-H2] as one of the RCs may not be the one associated with the lowest barrier. Closer examinations of the structure of IM1′ and IM2′

13096 J. Phys. Chem. B, Vol. 112, No. 41, 2008

Zhao and Liu

Figure 6. Two dimensional potential energy surfaces for step 2. RC1 ) d[C1-N2]-d[C1-N3] and RC2 ) d[N3-H2]-d[N2-H2]. (a) The AM1/ MM results. (b) The B3LYP/6-31G*/MM//AM1/MM results. Positions of various intermediates and transition states are indicated.

Figure 7. Two dimensional potential energy surfaces for the indirect proton transfer sub step in step 2. (RC1′, RC2′) ) (d[N2-H2]-d[O6-H2], d[O6-H2]-d[N3-H2]). (a) The AM1/MM results. (b) The B3LYP/6-31G*/MM//AM1/MM results. Positions of various intermediates and transition states are indicated.

TABLE 1: Comparisons of Various Interatomic Distances (In Angstrom) at IM2, TS3, and IM3 IM2 TS3 IM3

O2-H1

O5-H1

C3-O5

C2-H3

N2-H3

N2-H4

O2-H4

1.00 1.03 1.67

1.83 1.59 0.98

1.47 1.57 3.01

1.10 1.98 2.57

2.14 1.05 1.02

1.02 1.05 1.67

2.15 1.64 1.02

suggested that the carboxyl group of the substrate is at an ideal position to assist the proton transfer. As has been described in the method section, we constructed another two-dimensional PES linking IM1′ and IM2′ to find an alternative RC for the proton transfer sub step. The PESs using (RC1′, RC2′) ) (d[N2-H2]-d[O6-H2], d[O6-H2]-d[N3-H2]) are shown as Figure 7. Figure 7a shows the AM1/MM results, and Figure 7b shows the B3LYP/6-31G*// AM1/MM results. Both PESs suggested that the participation of the substrate carboxyl group may significantly lower the activation barrier relative to the direct proton transfer path. On the AM1/MM surface (Figure 7a), there is a minimum at (1.4, -1.3), corresponding to completed proton transfer to the

carboxyl group. The highest barrier exists between IM1′ and this local minimum with a height of 22.9 kcal mol-1, which is much lower than the 39.0 kcal mol-1 value computed for the direct proton transfer mechanism. The corresponding transition state TS2-2 is located at (1.2, -0.1), corresponding to a transition structure for the proton transfer between the amino and carboxyl groups of the substrate. The single point ab initio QM/MM PES (Figure 7b) is in qualitative agreement with the AM1/MM one, except that on the ab initio QM/MM PES, the carboxyl-protonated local minimum and TS2-2 merge into a single TS2-2. In Figure 7b, this transition state is located at (1.3, -0.3) and is associated with an activation barrier of 14.4 kcal mol-1, which is also much lower than the 34.0 kcal mol-1

Pyridoxal 5′-Phosphate-Dependent Enzyme L-Serine Dehydratase

J. Phys. Chem. B, Vol. 112, No. 41, 2008 13097

Figure 9. AM1/MM potential energy curves obtained using different snapshots as starting geometries for scheme 2 of step 3. The reaction coordinates are d[O2-H1] + d[C3-O5]-d[O5-H1] + d[C2-H3]-d[N2-H3] + d[N2-H4]-d[O2-H4] from IM2 to IM3.

Figure 8. Potential energy surfaces along the NEB paths for the two proton-relaying schemes in step 3, shown in Figure 3. Plots (a) and (b) are results for schemes 1 and 2, respectively. Circles are the AM1/MM results and squares are the B3LYP/6-31G*/MM// AM1/MM results.

barrier for the direct proton transfer mechanism computed at the same level of theory. In summary, applying the computationally efficient semiempirical QM/MM model has allowed for the construction of multiple two-dimensional PESs, which suggested a mechanism consisted of three sub steps for step 2 and a chemical role for the substrate carboxyl group in the second sub step. These results in turn led to the definition of a single RC for each sub step, which was then used to determine the final ab initio QM/MM MEPs described above. 3. The Proton-Relaying Mechanism in the Rate-Limiting Step 3. As has been described in the Materials and Methods, we have used the adapted NEB approach to explore two possible proton transfer schemes (Figure 3) for step 3. The resulting PES along the NEB paths is presented in Figure 8a for scheme 1 and 8b for scheme 2. At the AM1/MM level, these schemes show comparable activation barriers (54.4 kcal mol-1 for scheme 1 and 59.7 kcal mol-1 for scheme 2). However, at the B3LYP/ 6-31G*/MM//AM1/MM level, scheme 2 has a substantially lower barrier (33.4 kcal mol-1) as compared with scheme 1 (47.1 kcal mol-1). For scheme 1, the semiempirical and abinitio PESs are qualitatively similar. For scheme 2, the AM1/MM model produced a significantly higher barrier at a later stage, relative to the single point abinitio QM/MM results. This may be caused by an underestimation of the stability of the leaving hydroxyl group by the semiempirical AM1 model. For the abinitio QM/ MM, the barrier takes place before the leaving of the hydroxyl group. As the single point abinitio calculations have produced a much lower barrier for scheme 2 relative to scheme 1, we have defined the RC for step 3 based on scheme 2 to determine the final MEP. It turns out that the RCD results computed using the abinitio QM/MM model (Figure 4) are in good agreement with the single point ab initio QM/MM calculations on the AM1/ MM geometries (Figure 8b), in support of our approach defining RCs based on the NEB results. We have started from a single snapshot at 134.2 ps from the 200 ps MM simulations and performed the MEP optimizations.

We considered the second proton transfer scheme for step 3 as an example to investigate how sensitive the MEP derived using the semiempirical QM/MM would be to different starting structures. Two other snapshots approximately 70 ps before and after the 134.2 ps one, one at 60 ps and the other at 200 ps, were used separately to recompute the MEP paths using the NEB approach. The potential energy curves obtained using the three different starting structures are shown in Figure 9. The curves are qualitatively the same. The barrier heights are between 59.7 to 62.5 kcal mol-1. The energy differences between the product and the reactant states are between 7.1 to 8.7 kcal mol-1. We note that a number of previous QM/MM studies on other enzymatic systems have also compared MEPs computed using different starting structures sampled by simulations.26,60-62 Qualitatively similar MEPs were also reported, and the reported standard variations in energy barriers were between 1 to 3 kcal mol-1. We would also like to acknowledge that influences of environmental fluctuations on the MEPs have not been considered in this work. Our ab initio QM/MM results produced reasonable barrier heights, suggesting that the MEP calculations have captured the major effects of environmental reorganizations associated with the reaction process. In several earlier studies, the free energy perturbation QM/MM approach has been employed to estimate the contributions of MM environmental fluctuations.24-26,39,63,64 In the majority of cases, MEP provided reasonable first-order potential energy surfaces, with environmental fluctuations contributed less than a few kcal mol-1. These earlier examples and the results in Figure 9 suggest that the major mechanical insights obtained here would not be influenced by the approximation of neglecting environmental fluctuations. 4. Strategy for Transition State Stabilization Utilized by the Enzyme. Figure 10 shows various energy components of the final PES. The intrinsic QM energies (EQM in Figure 10) have been computed as the single point energies of the isolated QM subsystem at the B3LYP/6-31G* level. EMM are the total energies of the MM subsystem. The QM/MM interaction energies include both the van der Waals and the electrostatic part; the electrostatic part has been approximately determined using electrostatic potential-derived point charges on QM atoms and the MM charge. This approximation has been adopted only in analyzing the results. In QM/MM minimizations and computing the total QM/MM energies, the QM/MM electrostatic interactions have been considered as part of the electronic Hamiltonian. The energies are plotted as values relative to those of the reactant.

13098 J. Phys. Chem. B, Vol. 112, No. 41, 2008

Figure 10. Various energy components along the final minimum energy paths. Squares are the results of intrinsic QM subsystem, EQM; triangles are the MM subsystem, EMM; and circles are EQM/MM. The EQM/MM includes the vdw term and the electrostatic term between the QM subsystem and MM subsystem.

Figure 11. Decompositions of the QM/MM interaction energies at different transition states into contributions of individual residues. Values are relative to those of the reactant state. (a) For TS1; (b) for TS2-2; (c) for TS3. Labels beside spikes are residue numbers of major contributing residues.

Figure 10 suggests that the enzyme environment does not contribute much to stabilize the transition states involved in step 1 and step 2. However, the enzyme significantly stabilizes the transition state TS3 for step 3, by as large as 23.2 kcal mol-1. As stated before, step 3 is the rate-limiting step associated with the highest activation barrier. The observation that the computed TS3 is most strongly stabilized as compared with other structures along the path supports that TS3 closely approximate the actual rate-limiting TS, as this is in consistent with the theory that enzymes achieve efficient catalysis through stabilizing the ratelimiting TS. As first order approximations, the transition state stabilization effects of the enzyme environment can be decomposed into contributions of individual residues.65 In Figure , we plot the contributions of individual residues to stabilize TS1, TS2-2, and TS3, respectively. The sparse spikes in these plots indicated

Zhao and Liu that only a small number of residues contribute to transition state stabilization (negative spikes in Figure 11). They mostly locate within a few sequence fragments surrounding the active site. A general picture of how the enzyme environment achieved maximum stabilization of TS3 may emerge from detailed examinations of the interactions of these residues with the reacting QM part. The major contributing residues indicated in Figure 11 fall into different groups that locate within different sequence fragments, and most of them interact specifically with different parts of the QM subsystem (Figure 12). The first group of residues, residues 64-68, interacts specifically with the carboxyl group of the serine. In addition, Asn67 also interacts with O1 of PLP. We note that this pentapeptide sequence fragment is well-conserved in rSDH and in the SDHs of other species.1,18,66-71 The residues have stabilizing effects on all three major transition states (around -1.0 kcal mol-1 for TS1, -2.5∼-3.5 kcal mol-1 for TS2-2 and -5∼-8 kcal mol-1 for TS3), except that Ser64 and Asn67 have small (1.4∼1.5 kcal mol-1) destabilization effects on TS2-2. The moderate transition state-stabilizing effects on TS1 can be caused by the small increases in the negative charges on the carboxyl group, resulted from proton transfer from serine to the phosphate. The stabilizing effects on TS2-2 can be mostly attributed to positional rearrangements of the substrate serine accompanying its nucleophilic attack on the PLP-Lys41 Schiff base. This causes the serine to enter deeper into the active site, leading to shorter interaction distances between the substrate carboxyl and the backbone amides of Ala65, Gly66, and Ala68. The small destabilization effects of Ser64 and Asn67 on TS2-2 are consistent with the charge redistribution on the carboxyl group and the PLP ring at this state. At TS3, the negative charges on O6, O7, and O1 increase significantly, due to the elimination of the R proton of serine. Thus the residues interacting specifically with these groups contribute significantly to stabilize TS3. The second group of residues interacts specifically with the phosphate group. These include residues 168-172, which form the binding pocket for the phosphate group, providing five hydrogen bond donors. These interactions disfavor the protonation of the phosphate, thus destabilize the transition states relative to the reactant. However, these groups may increase the acidity of the protonated-phosphate group, making it a better balanced general acid-general base catalyst in step 3. Another important residue interacting with the phosphate is His145. Its positively charged side chain also disfavors protonation of the phosphate, thus destabilizing TS1 relative to the reactant. The large stabilizing effects on TS2-2 and TS3 may again be attributed to shorter interaction distances between His145 and the phosphate at these states. His145 also facilitates the phosphate to serve as a general acid to protonate the leaving hydroxyl group in step 3. The third group of residues interacts with the amino group of Lys41. These mainly include residues 135-138. At TS2-2 and TS3, the backbone carboxyl of Pro135 forms a hydrogen bond with the amino group of Lys41. In addition, one of its side chain CH2 is in close contact with the leaving hydroxyl group of the substrate serine. The distances from the negative carboxyl groups of Asp137 and Asp138 to the amino group of Lys41 shortened significantly in TS2-2 and TS3 relative to the reactant. The changes are from 9.7 to 8.4 Å for Asp137, 11.9 to 10.5 Å for Asp138 in TS2-2, 9.7 to 6.9 Å for Asp137, and 11.9 to 9.1 Å for Asp138 in TS3. They have large stabilization effects on the TS2-2 and TS3.

Pyridoxal 5′-Phosphate-Dependent Enzyme L-Serine Dehydratase

J. Phys. Chem. B, Vol. 112, No. 41, 2008 13099

Figure 12. (a) The active site structure of TS3. The QM part is represented as sticks. The substrate serine is behind the PLP moiety. Some of the major contributing residues to the QM/MM interactions are represented as lines and labeled by residue numbers. (b) Specific interactions between the QM part and the enzyme environment. Residues contribute significantly to stabilize TS3 are represented as thicker boxes.

Another outstanding residue implicated by Figure 11 is Ala222. This residue interacts specifically with the amino group of the substrate serine and has varying effects on the three major transition states. Besides residues interacting specifically with the QM subsystem, several residues at a distance from the QM system were also found to contribute non-negligibly to transition state stabilization, including Arg98, Asp114, and Glu194 (Figure 11). Their effects are mainly on TS2-2 and TS3 and through favorite electrostatic interactions.

mechanisms of enzyme catalysis may eventually advance the design of specific inhibitors and the reengineering of molecular catalysts. Acknowledgment. We thank Dr. Yingkai Zhang for help with the QM/MM pseudobond model. Financial support from the Chinese Natural Science Foundation (Grants 90403120 and 30670485) and from the Chinese Ministry of Science and Technology (Grants 2006AA02Z303) are acknowledged. References and Notes

Conclusions We have provided a detailed, energetically feasible chemical mechanism for the serine dehydration reaction catalyzed by rSDH. The following new insights into the chemical mechanism have been obtained: the carboxyl group chemically participates the Schiff base substitution step (step 2) to lower the proton transfer barrier; the elimination of the R-proton and the leaving of the hydroxyl group take place sequentially in a single step (step 3). In this step, the rate-limiting barrier arises because of the elimination of the R-proton. The subsequent leaving of the β-hydroxyl group is energetically a downhill process, assisted by a proton-relay mechanism via the PLP phosphate. Besides chemical catalysis, rSDH provides a well-organized environment to specifically stabilize the rate-limiting transition state TS3. At this state, the extracting of the R proton leads to localization of positive charges on the amino group of Lys41 and negative charges on the carboxyl group of the substrate serine. The transition state stabilization strategy of the enzyme seems to have taken advantage of such distinctive charge redistributions. Most of the enzymatic groups contributing nonchemically to catalysis provide specific interactions with these groups. A few distant residues also contribute to stabilize the redistributed charges. In summary, we have employed a strategy to study enzyme reactions involving complex reaction coordinates combining the physical accuracy of ab initio QM/MM models, the computational efficiency of semiempirical QM/MM models, and pathway optimization techniques. New insights into the chemical mechanism and transition state stabilization strategies of a PLPdependent enzyme have been obtained at the ab initio QM/MM level. Although there have been relatively few reported quantitative kinetic studies on this particular enzyme which would allow us to directly compare the computed barriers with experimental rate constants, our predictions about the roles of various residues should be testable by experiments such as site-directed mutagenesis. The ability to determine and test detailed catalytic

(1) Sun, L.; Bartlam, M.; Liu, Y. W.; Pang, H.; Rao, Z. H. Protein Sci. 2005, 14, 791. (2) Mozzarelli, A.; Bettati, S. Chem. ReV. 2006, 6, 275. (3) Alexander, F. W.; Sandmeier, E.; Mehta, P. K.; Christen, P. Eur. J. Biochem. 1994, 219, 953. (4) Holzer, H.; Cennamo, C.; Boll, M. Biochem. Biophys. Res. Commun. 1964, 14, 478. (5) Nakagawa, H.; Kimura, H. J. Biochem. 1969, 66, 669. (6) Reiber, H. Z. Physiol.Chem. 1974, 335, 1240. (7) Mehta, P. K.; Christen, P. AdV. Enzymol. Relat. Areas Mol. Biol. 2000, 74, 129. (8) Umbarger, H. E. Protein Sci. 1992, 1, 1392. (9) Changeux, J. P. BioEssays 1993, 15, 625. (10) Yamada, T.; Komoto, J.; Takata, Y.; Ogawa, H.; Pitot, H. C.; Takusagawa, F. Biochemistry 2003, 42, 12854. (11) Snell, K. AdV. Enzyme Regul. 1984, 22, 325. (12) Ebara, S.; Toyoshima, S.; Matsumura, T.; Adachi, S.; Takenaka, S.; Yamaji, R.; Watanabe, F.; Miyatake, K.; Inui, H.; Nakano, Y. Biochim. Biophys. Acta 2001, 1568, 111. (13) Imai, S.; Kanamoto, R.; Yagi, I.; Kotaru, M.; Saeki, T.; Iwami, K. Biosci. Biotechnol. Biochem. 2003, 67, 383. (14) Lopez-Flores, I.; Peragon, J.; Valderrama, R.; Esteban, F. J.; Luque, F.; Peinado, M. A.; Aranda, F.; Lupianez, J. A.; Barroso, J. B. Am. J. Physiol.: Regul. Integr. Comp. Physiol 2006, 291, R1295. (15) Yoshida, T.; Kikuchi, G. Biochem. Biophys. Res. Commun. 1969, 35, 577. (16) Kashii, T.; Gomi, T.; Oya, T.; Ishii, Y.; Oda, H.; Maruyama, M.; Kobayashi, M.; Masuda, T.; Yamazaki, M.; Nagata, T.; Tsukada, K.; Nakajima, A.; Tatsu, K.; Mori, H.; Takusagawa, F.; Ogawa, H.; Pitot, H. C. Int. J. Biochem. Cell Biol. 2005, 37, 574. (17) Snell, K.; Natsumeda, Y.; Eble, J. N.; Glover, J. L.; Weber, G. Br. J. Cancer 1988, 57, 87. (18) Ogawa, H.; Gomi, T.; Nishizawa, M.; Hayakawa, Y.; Endo, S.; Hayashi, K.; Ochiai, H.; Takusagawa, F.; Pitot, H. C.; Mori, H.; Sakurai, H.; Koizumi, K.; Saiki, I.; Oda, H.; Fujishita, T.; Miwa, T.; Maruyama, M.; Kobayashi, M. Biochim. Biophys. Acta 2006, 1764, 961. (19) Hyde, C. C.; Ahmed, S. A.; Padlan, E. A.; Miles, E. W.; Davies, D. R. J. Biol. Chem. 1988, 263, 17857. (20) Gallagher, D. T.; Gilliland, G. L.; Xiao, G.; Zondlo, J.; Fisher, K. E.; Chinchilla, D.; Eisenstein, E. Structure 1998, 15, 465. (21) Burkhard, P.; Rao, G. S.; Hohenester, E.; Schnackerz, K. D.; Cook, P. F.; Jansonius, J. N. J. Mol. Biol. 1998, 283, 121. (22) Meier, M.; Janosik, M.; Kery, V.; Kraus, J. P.; Burkhard, P. EMBO J. 2001, 20, 3910. (23) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227. (24) Liu, H. Y.; Zhang, Y. K.; Yang, W. T. J. Am. Chem. Soc. 2000, 122, 6560.

13100 J. Phys. Chem. B, Vol. 112, No. 41, 2008 (25) Cisneros, G. A.; Liu, H. Y.; Zhang, Y. K.; Yang, W. T. J. Am. Chem. Soc. 2003, 125, 10384. (26) Hu, P.; Zhang, Y. K. J. Am. Chem. Soc. 2006, 128, 1272. (27) Cheng, Y.; Cheng, X.; Radiac, Z.; McCammon, J. A. J. Am. Chem. Soc. 2007, 129, 6562. (28) Altun, A.; Shaik, S.; Thiel, W. J. Am. Chem. Soc. 2007, 129, 8978. (29) Allemann, R. K.; Young, N. J.; Ma, S. H.; Truhlar, D. G.; Gao, J. L. J. Am. Chem. Soc. 2007, 129, 13008. (30) Ma, S.; Devi-Kesavan, L. S.; Gao, J. J. Am. Chem. Soc. 2007, 129, 13633. (31) Sharma, P. K.; Chu, Z. T.; Olsson, M. H. M.; Warshel, A. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 9661. (32) Warshel, A. Computer Modeling of Chemical Reactions in Enzymes and Solutions; Wiley: New York, 1991. (33) Bash, P.; Ho, L.; MacKerell, A.; Levine, D.; Hallstrom, P. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 3698. (34) Thiel, W. J. Mol. Struct.:THEOCHEM 1997, 398, 1. (35) Hermann, J. C.; Hensen, C.; Ridder, L.; Mulholland, A. J.; Holtje, H. D. J. Am. Chem. Soc. 2005, 127, 4454. (36) Mills, G.; Jo’nsson, H. Phys. ReV. Lett. 1994, 72, 1124. (37) Jo’nsson, H.; Mills, G.; Jacobsen, K. W. Nudged elastic band method for finding minimum energy paths of transitions. In Classical and Quantum Dynamics in Condensed Phase Simulations; Berne, B. J., Ciccotti, G., Coker, D. F., Eds.; World Scientific: Singapore, 1998; pp 385. (38) Xie, L.; Liu, H.; Yang, W. J. Chem. Phys. 2004, 120, 8039. (39) Zhang, Y.; Liu, H.; Yang, W. J. Chem. Phys. 2000, 112, 3483. (40) Winget, P.; Selcuki, C.; Horn, A. H. C.; Martin, B.; Clark, T. Theo. Chem. Acc. 2003, 110, 254. (41) Bredow, T.; Jug, K. Theo. Chem. Acc. 2005, 113, 1. (42) Gregersen, B.; Lopez, X.; York, D. J. Am. Chem. Soc. 2004, 126, 7504. (43) Rossi, I.; Truhlar, D. G. Chem. Phys. Lett. 1995, 233, 231. (44) Nam, K.; Cui, Q.; Gao, J.; York, D. M. J. Chem. Theo. Compu. 2007, 3, 486. (45) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902. (46) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (47) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1996, 118, 2309. (48) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11, 700. (49) Zhang, Y. K.; Lee, T. S.; Yang, W. T. J. Chem. Phys. 1999, 110, 46. (50) Antes, I.; Thiel, W. J. Phys. Chem. A 1999, 103, 9290. (51) Zhang, Y. K. J. Chem. Phys. 2005, 122, 024114. (52) Stewart, J. J. P. J. Comput. Chem. 1989, 10, 221. (53) Dewar, M. J. S.; Thiel, W. J. Am. Chem. Soc. 1977, 99, 4899.

Zhao and Liu (54) Liu, H.; Elstner, M.; Kaxiras, E.; Franenheim, T.; Hermans, J.; Yang, W. Proteins: Struct., Funct., Genet. 2001, 44, 484. (55) Ponder, J. W. TINKER, a Software Tools for Molecular Design; Version 3.6 ed.; (the updated version for the TINKER program can be obtained from the webpage at http://dasher.wustl.edu/tinker, maintained by J. W. Ponder.), accessed February, 1998. (56) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; AlLaham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03; ReVision C.02 ed.; Gaussian, Inc.: Wallingford, CT, 2004. (57) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (58) Berendsen, H. J. C.; Postma, J. P. M.; Gunstereno, W. F. V.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (59) Cavestri, R. C.; Fedor, L. R. J. Am. Chem. Soc. 1970, 92, 4610. (60) Zhang, Y.; Kua, J.; McCammon, J. A. J. Phys. Chem. B 2003, 107, 4459. (61) Claeyssens, F.; Ranaghan, K. E.; Manby, F. R.; Harvey, J. N.; Mulholland, A. J. Chem. Commum. 2005, 5068. (62) Klahn, M.; Braun-Sand, S.; Rosta, E.; Warshel, A. J. Phys. Chem. B 2005, 109, 15645. (63) Jorgensen, W. L. Acc. Chem. Res. 1989, 22, 184. (64) Kuhn, B.; Kollman, P. A. J. Am. Chem. Soc. 2000, 122, 2586. (65) Boresch, S.; Archontis, G.; Karplus, M. Proteins: Struct., Funct., Genet. 1994, 20, 25. (66) Datta, P.; Goss, T. J.; Omnaas, J. R.; Patil, R. V. Proc. Natl. Acad. Sci. U.S.A. 1987, 84, 393. (67) Ogawa, H.; Konishi, K.; Fujioka, M. Biochim. Biophys. Acta 1989, 996, 139. (68) Samach, A.; Hareven, D.; Gutfinger, T.; Ken-Dror, S.; Lifschitz, E. Proc. Natl. Acad. Sci. U.S.A. 1991, 88, 2678. (69) Bornaes, C.; Petersen, J. G. L.; Holmberg, S. Genetics 1992, 131, 531. (70) Eisenstein, E.; Yu, H. D.; Fisher, K. E.; Iacuzio, D. A.; Ducote, K. R.; Schwarz, F. P. Biochemistry 1995, 34, 9403. (71) Ogawa, H. Trends Comp. Biochem. Physiol. 2000, 6, 1.

JP802262M