A Comprehensive Understanding of Enzymatic Catalysis by

Research Article pubs.acs.org/acscatalysis

A Comprehensive Understanding of Enzymatic Catalysis by Hydroxynitrile Lyases with S Stereoselectivity from the α/β-Hydrolase Superfamily: Revised Role of the Active-Site Lysine and Kinetic Behavior of Substrate Delivery and Sequential Product Release Yuan Zhao,† Nanhao Chen,‡ Chaojie Wang,*,† and Zexing Cao*,§ †

The Key Laboratory of Natural Medicine and Immuno-Engineering, Henan University, Kaifeng 475004, People’s Republic of China ‡ School of Pharmaceutical Sciences, Sun Yat-sen University, Guangzhou 510006, People’s Republic of China § State Key Laboratory of Physical Chemistry of Solid Surfaces and Fujian Provincial Key Laboratory of Theoretical and Computational Chemistry, College of Chemistry and Chemical Engineering, Xiamen University, Xiamen 360015, People’s Republic of China S Supporting Information *

ABSTRACT: The highly homologous hydroxynitrile lyases from Manihot esculent (MeHNL) and Hevea brasiliensis (HbHNL) both belong to the α/β-hydrolase superfamily, and they convert cyanohydrins into the corresponding ketone (aldehyde) and hydrocyanic acid, which is important for biosynthesis for carbon−carbon formation. On the basis of extensive MM and ab initio QM/MM MD simulations, one-dimensional and twodimensional free energy profiles on the whole enzymatic catalysis by MeHNL have been explored, and the effects of key residues around the channel on the delivery of substrate and product have been discussed. The residue Trp128 plays an important gate-switching role to manipulate the substrate access to the active site and product release. In particular, the release of acetone and HCN has been first detected to follow a stepwise mechanism. The release of HCN is quite facile, while the escape of acetone experiences a barrier of ∼10 kcal/mol. The chemical reaction is an endergonic process with a free energy barrier of ∼17.1 kcal/mol, which dominates the entire enzymatic efficiency. Such energy costs can be compensated by the remarkable energy release during the initial substrate binding. Here the carbon−carbon cleavage is the rate-determining step, which differs from that of HbHNL. The protonation state of Lys237 plays an important role in carbon−carbon bond cleavage by restoring the Ser80Ala mutant system to the wild system, which explains the discrepancy between MeHNL and HbHNL at the molecular or atomic scale. The present results provide a basis for understanding the similarity and difference in the enzymatic catalysis by MeHNL and HbHNL. KEYWORDS: hydroxynitrile lyases, α/β-hydrolase superfamily, QM/MM MD, mechanism, biosynthesis, product release HbHNL of the α/β-hydrolase superfamily has been explored, both experimentally and theoretically,6,7,21,26−33 but the similarity and difference of a variety of HNLs with S stereoselectivity are still less well understood. For MeHNL, various experimental research has focused on crystal structure identification, mechanism prediction, separation methods, and applications to asymmetric cyanohydrination.16,18−25 With regard to the enzymatic mechanism, Wajant et al. identified three crucial residues (Ser80, Asp208, and His236) to be

1. INTRODUCTION Hydroxynitrile lyase (HNL) catalyzes the cleavage of cyanohydrins into the corresponding aldehyde (or ketone) and hydrocyanic acid (HCN), which plays an important role in plant protection systems and asparagine biosynthesis.1−3 This reaction is reversible, and its reverse process is also of great potential significance in the biocatalytic retrosynthesis of organic compounds.4−17 Herein, the HNLs from Manihot esculent (MeHNL) and Hevea brasiliensis (HbHNL) both belong to the α/β-hydrolase superfamily, and they can catalyze the synthesis of chiral molecules with high S stereoselectivity, showing potential values in biosyntheses of pharmaceuticals and agrochemicals.16,18−25 © XXXX American Chemical Society

Received: December 15, 2015 Revised: January 31, 2016

2145

DOI: 10.1021/acscatal.5b02855 ACS Catal. 2016, 6, 2145−2157

Research Article

ACS Catalysis important for the enzyme activity,20 and Lauble et al. obtained the Ser80Ala mutant complex of MeHNL with acetone cyanohydrins (PDB ID: 1E8D), which contains 520 residues in the asymmetric unit composed of chains A and B.18 Previous studies proposed a similar catalytic mechanism, shown in Figure 1.18,20 First, the substrate is deprotonated by Ser80 and

in the biocatalytic retrosynthesis industry, and probably, the delivery of a carbonyl compound in an enzyme is not the same as that of HCN. Furthermore, the suggested distinct catalytic roles of lysine residue experimentally in both HbHNL and MeHNL need to be interpreted. In order to get insight into the whole catalytic process by MeHNL, quite a few crucial issues need to be understood with the aid of QM/MM MD simulations. Examples are as follows. (i) For the proposed different roles of the active-site lysine residues (Lys237 in MeHNL and Lys236 in HbHNL) in the two highly homologous MeHNL and HbHNL, does the crystal structure of the MeHNL-Ser80Ala mutant enzyme provide enough information for the mechanistic speculation? What is the effect of different protonation states of Lys237 on the enzymatic catalysis? (ii) For the dynamic properties of substrate binding and product (acetone and HCN) release, does the release of acetone and HCN follow a sequential or concerted mechanism? What is the role of the key residues around the delivery channel? (iii) What are similarities and differences in the catalytic mechanism, efficiency, and substrate (or product) delivery between MeHNL and HbHNL? Herein extensive combined quantum mechanics and molecular mechanics molecular dynamics (QM/MM MD) simulations and classical molecular dynamics (MD) simulations on the overall catalytic process by MeHNL, including the catalytic mechanism, substrate binding, and acetone and HCN release, have been performed to achieve a more in-depth understanding of the enzymatic catalysis by HNLs of the α/β-hydrolase superfamily.

Figure 1. Experimentally proposed catalytic mechanism based on the crystal structure of the Ser80Ala mutant.18,20

His236; meanwhile, the carbon−carbon bond is broken. Afterward, the His236 provides a proton for the hydrocyanic acid formation. It should be noted that Lauble et al.18 ruled out the effect of active-site lysine residue (Lys237) on the basis of three-dimensional X-ray structural data of the Ser80Ala mutant complex. Three years later, Gruber et al.21 identified the crystal complex of HNL from Hevea brasiliensis (HbHNL) with the biological substrate acetone cyanohydrin by a freeze−quench method. They found that HbHNL and MeHNL are highly homologous (77% sequence identity), which may lead to the similar binding style in the active site and molecular mechanism. However, they proposed that the active-site residue (Lys236) plays an important role in HbHNL but not in MeHNL. With regard to the substrate binding and product release, the X-ray crystal structure of wild-type MeHNL reveals a narrow channel through the surface of the protein to the active site.19 Bühler et al. has identified that the Trp128Ala mutant system can enlarge the access channel to the active site of MeHNL and facilitate the entrance of sterically hindered substrates to the active site.23,25 For HbHNL, similar consequences are proposed;26,27 moreover, Wagner et al.13 suggested that the ketone and HCN would enter (or leave) in a sequential fashion for HbHNL on the basis of the crystal structure. However, for these HNLs, detailed mechanistic aspects for substrate binding and product (acetone and HCN) release, as well as structural and kinetic effects of related residues, are still open. Theoretical knowledge about the catalytic mechanism and substrate (or product) delivery for MeHNL has been lacking up to now. Very recently, we performed extensive QM/MM and MM MD simulations on the catalytic process by HbHNL. Our calculated results show good agreement with the experimental data,31 and the importance of protonated Lys236 in the catalytic chemical step and the nontrivial role of Trp128 in the substrate delivery have been established, on the basis of calculations. Such results of HbHNL provide a basis for comparison with MeHNL, whereas further investigations are highly required for a comprehensive understanding of the whole enzymatic catalysis. For example, the release of product acetone was not considered in previous calculations, which is definitely a matter

2. COMPUTATIONAL DETAILS 2.1. Preparation of the Enzyme−Substrate Complex Model. The initial structure of the enzyme−substrate complex was built on the basis of the crystal structure of the Ser80Ala mutated MeHNL complex with acetone cyanohydrin (PDB ID: 1E8D)18 by restoring alanine to serine. The protonation states of ionizable residues were determined via the H++ program34−37 and the neighboring hydrogen bond networks. The acetone cyanohydrin and the protein were described by the AMBER GAFF force field (GAFF)38 and the AMBER99SB force field,39,40 respectively. The partial atomic charges of ligand were calculated by the restrained electrostatic potential (RESP) charge41 at the HF/6-31G(d) level with the GAUSSIAN 09 program.42 The whole system was solvated into a 94 × 93 × 105 Å rectangular box of TIP3P water43 with a 15 Å buffer distance between the box edge and the nearest solute atoms. The system was neutralized by adding Na+ ions. The protons were added automatically, and topology parameters as well as initial coordinates were generated by the tleap Amber tool. After the energy minimization, the system was heated gradually from 0 to 300 K for 50 ps, and another 50 ps MD simulation was carried out to relax the system density to about 1.0 g cm−3. Finally, the standard 10 ns MD simulation under NVT ensemble was performed with an integration time step of 1 fs via use of the periodic boundary condition. The cutoff value was set to 12 Å for van der Waals and electrostatic interaction calculations. Long-range electrostatic interactions were dealt with by the particle mesh Ewald (PME) method.44,45 The Langevin method was used to maintain the temperature at 300 K. All the bonds involving hydrogen atoms were constrained by using the SHAKE scheme.46 The root-mean-square deviation (RMSD) was used to estimate the stability of the backbone for the enzyme, and the last 2 ns trajectories were utilized for the ligand−residue interaction decomposition of the 2146


Research Article

ACS Catalysis two subunits. Herein, the MMPBSA.py module47 has been employed, which can calculate streamlining end-state free energy by using ensembles of representative structures, derived from molecular dynamics (MD) or Monte Carlo (MC) simulations. Several implicit solvation models are available. Since Hou et al.48−50 found that the MM-GBSA method using the generalized born model is suitable for both binding-pose predictions and binding free energy estimations, in this work the method has been used to analyze the binding free energy for the protein−ligand interactions. The per-residue decomposition was employed to estimate the contribution of each residue to the total free energy, which can provide reliable results in comparison with experimental alanine scanning data.51 The whole MD simulations are accomplished by applying the AMBER 12 software.52 2.2. QM/MM MD Simulations. In order to explore the distinct effect of the active-site lysine residue in MeHNL and HbHNL, the performances of different protonation states of Lys237 for MeHNL have been discussed. Two QM/MM models with different protonation states of the active-site lysine residue are constructed by deleting solvent molecules beyond 24 Å of the C1 atom of substrate on the basis of the snapshot obtained from the classical MD trajectory (models A and B refer to the protonated and unprotonated Lys237, respectively), as shown in Figure 2. The QM subsystem includes the side

Figure 2. QM/MM models for the MeHNL enzyme with different protonation states of Lys237.

Figure 3. (a) Binding mode of the active site for chain A. (b) Free energy decomposition based on the per-residue type for chain A. (c) Energy decomposition of the key residues for chain A.

chains of Thr11, Ser80, Asp208, His236, Lys237, and acetone cyanohydrin, which are described by using the B3LYP functional and the 6-31G(d) basis set. This approach has been successfully applied to other enzymatic reactions,53−60 including the HbHNL system.31 The MM part was treated with the AMBER99SB force field,39,40 which is the same as that used in the above MD simulations. The TIP3P model43 was used for water molecules. The QM/MM boundary was handled by the pseudobond approach61−63 with improved parameters.64 The spherical boundary condition was applied, in which all atoms within a radius of 20 Å from the spherical center atom were allowed to move freely. The cutoffs of 12 and 18 Å were used for van der Waals and electrostatic interactions among MM atoms, respectively. Following the interactive minimization procedure, the minimum reaction energy path was mapped out

by using the reaction coordinate driving method (RCD).65 Then a 500 ps MD simulation was performed for each structure identified along the path with the QM subsystem frozen. The obtained snapshot was used for the subsequent QM/MM MD simulations with umbrella sampling.66 For each window, 20 ps QM/MM MD simulation was carried out with a time step of 1 fs by Beeman algorithm.67 The Berendsen thermostat68 was used to control the temperature at around 300 K. The last 5 ps was used to obtain the potential of mean force (PMF) by the weighted histogram analysis method (WHAM).69,70 Herein the appropriate step size and the force of the biasing harmonic potential were chosen to guarantee the sufficient sampling. All QM/MM calculations were carried out by using modified Q-CHEM 4.071 and TINKER 4.2 programs.72 2147


Research Article

ACS Catalysis

Figure 4. Structures of the reactant, intermediate, and transition states for model A. The average distances for selected bonds are given in Å.

2.3. RAMD-MD and MM MD Simulations. In order to get the possible delivery channels for the reactant and product, the combined random acceleration molecular dynamics and MD simulations (RAMD-MD)73,74 were performed by using NAMD 2.9 software.75 The protein and the ligand were treated with the AMBER99SB39,40 and GAFF38 force fields, respectively. The initial structure of the enzyme−substrate (acetone cyanohydrin) system was constructed by using the snapshot from MM MD simulations. The enzyme−product (acetone and HCN) structure was obtained from the standard MD simulations (the same approach as used in section 2.1), on the basis of the sphere product state from QM/MM MD simulations (see Figure S1 in the Supporting Information). In RAMD-MD simulations, a random force was added to the center of mass. After a certain period of time, when the ligand did not reach the threshold parameter, a new random direction was selected; otherwise, the force direction was maintained. When the ligand fled away from the initial position, the classic MD simulations were performed to equilibrate the sampling, which avoided the error arising from the higher random force. Herein, random accelerations of 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, and 0.20 kcal Å−1 g−1 and a threshold distance value of 0.005 Å were added on the C1 atom of acetone cyanohydrins, the C1 atom of acetone, and the C atom of HCN, and 28 RAMD MD trajectories for each model were identified. On the basis of the results of RAMD MD simulations, the most probable pathway for each model (the enzyme−substrate and the enzyme− product) could be analyzed. By using the umbrella sampling technique,66 one-dimensional and the two-dimensional free energy profiles could be identified for substrate access and product release, respectively. The initial models were built on the basis of the snapshot from conventional MD simulations. Afterward, at least 8 ns MD simulations with the appropriate biasing harmonic potential were carried out for each window along the path of the substrate movement. By using the WHAM technology69,70 with the probability distribution, the free energy profiles were subsequently produced.

3. RESULTS AND DISCUSSION 3.1. Catalytic Mechanism. 3.1.1. Substrate−Enzyme Complex. In order to simulate the real protein environment, the mutated MeHNL−substrate system has been restored to the complex of the wild-type enzyme with substrate, and longtime MM MD simulations have been carried out. On the basis of the root-mean-square deviation (RMSD) results of the backbone atoms, the system approaches equilibrium after ∼6 ns (see Figure S2 in the Supporting Information). The last 2 ns sampling snapshots have been used for the following analysis. Since the two chains in MeHNL are essentially symmetrical, here only one chain will be used to evaluate the binding free energy and the ligand−residue interaction decomposition. As we can see in Figure 3, Thr11, Ile12, Ser80, Leu158, His236, and Lys237 play important roles in stabilizing the substrate in the active site. For the polar amino acid residues, such as Thr11, Ser80, His236, and Lys237, the electrostatic interactions play a predominant role in substrate binding, while for the nonpolar amino acid residues, such as Ile12 and Leu158, the van der Waals interactions are more important. In comparison with the wild system, the binding free energies of Ser80Ala and Lys237Ala mutant systems decrease by 7.1 and 19.1 kcal/mol, respectively, which show that the two residues have crucial effects on the substrate binding. It is worth noting that when the Ser80Ala mutant complex was repaired to the wide-type system, the distance between −NH3+ of Lys237 and −CN of the substrate was shortened notably from 3.5 to 3.1 Å after minimization and approximately remained unchanged during the long equilibrium simulations (see Figure S3 in the Supporting Information), suggesting that there are standing electrostatic interactions between them accordingly. This implies that Lys237 is probably involved in the catalytic process, and thus the proposed catalytic mechanism based on the X-ray structure of the mutant system cannot be universalized. 3.1.2. Comparison of the Catalytic Mechanisms for MeHNL and HbHNL. Experimentally, the active-site lysine in HbHNL is proved to be protonated by site-directed mutagenesis and difference Poisson−Boltzmann calculations.21 2148


Research Article

ACS Catalysis

Figure 6. Statistics of O1−H1, N1−H2, and C1−C2 bond lengths from snapshots along RC1.

used for discussion of the corresponding catalytic mechanism. By use of QM/MM MD simulations and umbrella sampling, the free energy profiles and corresponding configurations of reactant, transition state, intermediate, and product involved in the enzymatic catalysis by MeHNL have been determined. As shown in Figures 4 and 5, the catalytic reaction can be divided into three steps: (i) substrate deprotonation by Ser80 and His236, (ii) carbon−carbon bond cleavage, and (iii) HCN formation with the aid of Ser80. Different reaction coordinates have been tested (see Figure S5 in the Supporting Information). Considering the computational cost and accuracy, the corresponding reaction coordinates are finally defined as −(dO1−H1 + dN1−H2) (RC1), dC1−C2 (RC2), and dO1−H1 − dC2−H1 (RC3), respectively. During the first step, the substrate deprotonation occurs synchronously with the proton transfer from Ser80 to His236. The C1−O2 bond decreases from 1.42 ± 0.03 to 1.37 ± 0.03 Å, and O2δ−−C1δ+ is detected at the IM1 state. Asp208 is close to His236 apparently, which strengthens the electron-withdrawing ability of His236 and Ser80. The strong hydrogen bond interaction and electrostatic attraction can facilitate the proton transfer. The distance between His236 and Asp208 is smaller than that in HbHNL, and both His+···Asp− and His···Asp interactions exist here, which may speed up the rate of catalysis and stabilize the intermediate. As Figure 5a shows, the free energy barrier is ∼4.7 kcal/mol, which is lower than that in HbHNL by ∼2.0 kcal/mol, as expected above. For both HbHNL and MeHNL, from the reactant state B to IM1, multiple proton transfers are involved, while the C1−C2 bond is slightly changed. To clarify the sequence of proton transfer and carbon−carbon bond cleavage in catalysis by MeHNL, the entire snapshots along RC1 have been collected. As Figure 6 shows, the C1−C2 bond distance remains at around 1.5 Å along RC1, while the O1−H1 and N1−H2 bond distances shorten noticeably. It is clear that the substrate deprotonation and carbon−carbon bond cleavage follow a stepwise mechanism. However, there is a quasi-concerted mechanism for this step in HbHNL, as shown in our previous study, since such statistics analysis has not been performed. Herein, we consider that both HNLs should obey similar mechanisms for the proton transfer and carbon−carbon cleavage stage, while their main differences are that there are a more stable IM1 and a lower free energy barrier for the MeHNL.

Figure 5. (a) Free energy profile of the proton transfer for model A. (b) Free energy profile of the carbon−carbon cleavage for model A. (c) Free energy profile of the hydrocyanic acid formation for model A.

However, for MeHNL, the role of lysine was ruled out on the basis of the Ser80Ala mutant crystal, and the different protonated state of lysine received less attention.18,19 By employing the H++ program,34−37 we obtained the titration curve of lysine for both HNLs (see Figure S4 in the Supporting Information) and found that the results were very similar, which implies that the lysine is protonated for both HNLs. Thus, model A can be 2149


Research Article

ACS Catalysis In the second step, the carbon−carbon bond is completely broken to generate acetone and cyanide anion. The hydrogen bond between the acetone moiety and Thr11 survives, and meanwhile, the newly formed cyanide anion is stabilized by Lys237 and Ser80 via hydrogen bond interactions. It is worth mentioning that the distance between −NH3+ of Lys237 and −CN of the substrate slightly fluctuates around 2.50 Å from B to IM1, and its remarkable change from 2.46 ± 0.45 to 1.74 ± 0.11 Å appears at the IM1 → IM2 step, which is consistent with the evolution of the carbon−carbon bond length. Hence, the protonation of Lys237 plays an extremely crucial role in the carbon−carbon bond cleavage, which supplies the driving force by strong hydrogen bond and electrostatic polarization interactions. What is noteworthy is that the (CN)−···(NH3)+− Lys236 zwitterionic intermediate has been formed at the IM2 state, as found in HbHNL. However, there are different orientations of CN− in MeHNL and HbHNL. For MeHNL, the C2 atom points toward the hydrogen atom of Ser80, while for HbHNL, it points toward the hydrogen atom of His235, which may modify the following catalytic mechanisms to some extent. The free energy barrier for the carbon−carbon cleavage is ∼10.5 kcal/mol (see Figure 5b), which is higher than that of HbHNL by ∼3.9 kcal/mol. Such a barrier difference may be caused by the electrostatic interactions between −NH3+ of Lys237 and −CN of the substrate. Their distance is 2.46 Å in IM1 for MeHNL, which is longer than that in HbHNL by 0.25 Å. The shorter distance in HbHNL suggests a relatively stronger polarization interaction, in comparison to MeHNL, which may activate the C1−CN bond to facilitate the bond cleavage. In a previous study, we found that the longer the distance between −NH3+ of Lys236 and −CN of the substrate, the higher the barrier for the mutant systems of HbHNL.31 In the third stage, the CN− is close to Lys237 and Ser80, and thus it may attract protons from both residues to generate HNC and HCN, respectively. For HNC formation, we defined dLysN‑LysH−dN3‑LysH as the reaction coordinate. The QM/MM scan calculations have been carried out, and no stable intermediate appears with a continuing rise in energy, which means that the reaction to yield HNC hardly occurs (see Figure S6 in the Supporting Information). With regard to the HCN formation, the relative energy profiles exhibit one intermediate and a low barrier (see Figure S5 in the Supporting Information), and thus the QM/MM MD simulations combined with umbrella sampling has been used to discuss the detailed mechanism and corresponding free energy profile. As we can see in Figure 4, the CN− abstracts the proton of Ser80 to generate HCN, and His236 delivers its own proton-N1 to Ser80, together with the proton abstraction of Asp208. To illustrate the sequence of the proton transfers, the entire snapshots along RC3 have been collected to make an analysis. As shown in Figure 7, all of the proton transfers occur simultaneously. The His+···Asp− and His···Asp interactions become His···Asp−; moreover, the hydrogen bond distance between the cyanide anion and Lys237 changes from 1.74 ± 0.11 to 2.82 ± 0.40 Å. We should note that the free energy barrier for this step is ∼8.7 kcal/mol (see Figure 5c), which is lower than that of HbHNL by ∼4.4 kcal/mol. Both of the HNLs should overcome the strong ion-pair bonding interaction of (CN)−···(NH3)+−Lys236 during the formation of HCN, and their main difference is the proton source. For HbHNL, the cyanide anion abstracts the proton of histidine, along with the loss of the hydrogen bond between His235 and Ser80, while for MeHNL, the cyanide

Figure 7. (a) Statistics of C2−H1 and O1−H1 distances from snapshots along RC3. (b) Statistics of N1−H2 and O1−H2 distances from snapshots along RC3. (c) Statistics of O3−H3 and N2−H3 distances from snapshots along RC3.

anion will extract the proton of Ser80 without any hydrogenbond breaking. Therefore, the barrier for this step in MeHNL is lower than that in HbHNL. In the catalytic reaction by MeHNL, the rate-limiting step is carbon−carbon cleavage, differing from HbHNL with the rate-determining step of HCN formation. Here all three steps are endothermic, and these energy costs may be compensated by the energy release during the substrate binding. What is more, the product release can 2150


Research Article

ACS Catalysis

Figure 8. Structures of the reactant, intermediate, and transition states for model B. The average distances for selected bonds are given in Å.

not generated, which can be ascribed to the absence of the ion pair (CN)−···(NH3)+−Lys236. The His+···Asp− interaction is converted into His···Asp− along with proton transfer via His236 → Ser80 → − CN, and HCN is formed. The predicted free energy barrier is ∼22.8 kcal/mol (see Figure 9), much higher than that in model A by ∼8.4 kcal/mol. Clearly, the different protonation states of Lys237 may bring about different mechanisms and influence the catalytic efficiency to some extent, showing that Lys237 has catalytic roles similar to those of Lys236 in HbHNL. However, in comparison with HbHNL, the transition-state platform in the second step is not found for model B of MeHNL, which may arise from the different proton sources for HCN formation. For HbHNL, the transition-state configurations vary greatly during the sampling, and the proton may come from either His236 or Ser80, whereas for MeHNL, the transition-state configurations are relatively stable and the proton comes from Ser80. 3.2. Delivery Pathways of Substrate and Products. 3.2.1. Substrate Binding. In many enzymatic catalysis processes, the rate-limiting step is related to the substrate delivery to the active site or product release.59,76 For MeHNL, the transportation of aldehyde (ketone), HCN, and the corresponding cyanohydrins are very important in the biosynthesis of pharmaceuticals and agrochemicals. Herein, possible transportation pathways and their thermodynamic and dynamic properties are discussed. For the substrate binding, three possible channels (defined as paths A−C) are detected by means of random acceleration molecular dynamics simulations, which are displayed in Figure 10, as well as the residues involved. The active site is deeply buried inside the enzyme, and the distances for paths A−C are about 12, 12, and 15 Å, respectively. Moreover, the pocket width of path A is larger than those of paths B and C; thus, path A may be the most favorable channel. The probabilities for the delivery of substrate and product through each channel are summarized in Table 1. As Table 1 shows, path A, located between residues

also facilitate the reaction process. We will discuss both of the processes in the following section. 3.1.3. Role of the Active-Site Lysine Residue. In previous experimental studies, the active-site lysine residue (Lys236 in HbHNL and Lys237 in MeHNL) has been proposed to have different roles in catalysis by HbHNL and MeHNL. The activity of HbHNL depends on the protonation state of Lys236, while for MeHNL, Lys237 was assumed to be less important for catalysis on the basis of the three-dimensional crystal structure of the Ser80Ala mutated enzyme−substrate complex. The catalytic role of Lys236 in HbHNL has been confirmed in our previous calculations. Herein, we restore the Ser80Ala mutated MeHNL into the wild enzyme at first and support that the protonation state of Lys237 plays a curial role in the enzymatic catalysis by QM/MM MD simulations. To further understand the performance of the active-site lysine Lys237, its deprotonation state (model B) has been also investigated. The reaction coordinates of dO1−H1 + dN1−H2 (RC4) and dC1−C2 (RC5) are used to depict the catalytic process. The structures along with corresponding static bond lengths of key states and the free energy profiles are given in Figures 8 and 9. The catalytic reaction can be divided into two steps: (i) substrate deprotonation by Ser80 and His236 and (ii) carbon−carbon cleavage and HCN formation. In the first step, the proton transfers take place through substrate → Ser80 → Asp208, and the C1−O2 bond is slightly shortened. The charged configuration of the His235···Asp208 dyad is changed from His···Asp− to His+···Asp−. The distance between His235 and Asp208 is longer than that in Model A, which suggests a weaker driving force for the proton transfer. Moreover, there is no decent hydrogen bond between −NH2 of Lys237 and −CN of the substrate. Therefore, the formation of the intermediate and the proton transfer are more difficult than those in model A with the protonated Lys237. In the second step, the carbon−carbon bond is entirely broken, and acetone is formed. The cyanide anion intermediate as found in model A is 2151


Research Article

ACS Catalysis

is closed through the gate movement of Trp128 and Met147, which blocks access of the substrate. In the second stage (11.0 Å < RC ≤ 13.5 Å), Trp128 is flipping gradually, and the channel opens to let the substrate approach the active site. Herein, on the basis of the distance and angle of the −CN group and ring of Trp128, the lone pair···π and π···π interactions can be found between the conjugate ring of Trp128 and the −CN group of the substrate. The electrostatic interaction can be also found between Met147 and the substrate. There is a barrier at this stage, which may come from the flip of Trp128 and the breaking of the hydrogen bond network. In the third stage (8.5 Å < RC ≤ 11.0 Å), the substrate will pass through Trp128 and Met147. −CH3···π interactions are found between the −CH3 group of the substrate and the ring of Trp128, which drives the Trp128 to recover gradually. A hydrogen bond network is formed between the −OH group of Thr11 and the −CN group of the substrate, which helps the substrate access the active site. In the final stage (6.5 Å ≤ RC ≤ 8.5 Å), the substrate is binding at the original active site by Thr11, Ser80, and Lys237 via hydrogen bond interactions. The side-chain residues Ile12, Leu149, Leu158, and Ile210 constitute a hydrophobic pocket to accommodate the two −CH3 groups of the substrate, and the channel is blocked by the residues Trp128 and Met147. For HbHNL, a similar process and the gating role of Trp128 have been found in previous studies.31 3.2.2. Sequential Release of Products. For the acetone release, three possible channels have been discovered, which are similar to those of the substrate entrance. For hydrocyanic acid, the same channels are detected. The trajectory numbers and corresponding possibilities are given in Table 1. Obviously, for the release of both products, path A, located between residues 180−190 and Trp128, is the predominant pathway, and it will be discussed in detail. What is more, in order to get insight into their release sequence, MM MD simulations in combination with umbrella sampling technology for the two-dimensional (2D) free energy surface have been performed, and the free energy profile is shown in Figure 13. For the acetone release, the reaction coordinate is defined as the distance between the CD atom of Lys237 and the carbonyl carbon of acetone (RC7), and for the hydrocyanic acid release, the reaction coordinate is chosen as the distance between the CD atom of Lys237 and the carbon atom of hydrocyanic acid (RC8) (see Figure S8 in the Supporting Information). As shown in Figure 13, three mechanisms will be discussed. (I) HCN releases first, and then acetone moves away from the active site. We note that when acetone stays at the active site, the HCN undergoes a long minimum platform of ∼3 kcal/mol and then easily escapes from the protein by overcoming a barrier of ∼2 kcal/mol. After HCN departs from the protein, the release of acetone experiences a free energy barrier of about 10 kcal/mol, and no transition-state platform appears. Overall, this mechanism is feasible. (II) The acetone releases before HCN. In this process of acetone escape, the free energy barrier is higher and higher, and there is no occurrence of any local minimum point on the free energy surface, showing that the mechanism is not reasonable. (III) The acetone and HCN release simultaneously. We note that that a long transition platform shows up in the 2D free energy profiles, which is around 10 kcal/mol. The mechanism may be probable, although it is less favorable than the first mechanism. The stepwise mechanism I for the product release will be thus discussed in detail, and the corresponding key states are shown in Figure 14.

Figure 9. (a) Free energy profile of the proton transfer for model B. (b) Free energy profile of carbon−carbon cleavage and hydrocyanic acid formation for model B.

180−190 and Trp128, has a dominant share of trajectories (∼82.1%). Only four and one trajectories arise from paths B (∼14.3%) and C (∼3.6%), respectively. The predominant channel, path A, will be discussed in detail in the following part. According to the escaping direction of acetone cyanohydrins, the distance between the C1 atom of the reactant and the CE atom of Lys237 is selected as the reaction coordinate (RC6) (see Figure S7 in the Supporting Information). Classical MD simulations with the umbrella sampling technology have been employed to determine the free energy. As shown in Figure 11, the free energy profiles for different sampling time durations are very similar, showing reliable convergence of MD simulations for the free energy profiles. The reactant binding is a notably exothermic and low-barrier process, which drives the subsequent chemical reaction. On the basis of the conformational changes, the process can be approximately divided into four stages, as shown in Figure 12. In the first stage (RC > 13.5 Å), the substrate approaches the rim of the protein channel. With the moving of the substrate, a hydrogen bond network may be formed between the −OH group of substrate and the −NH2 group of Gln118; moreover, a noncovalent interaction may be detected between the ring of Trp128 and the −CN group of the substrate. Both of these interactions will promote the substrate entrance into the protein channel. At this moment, the channel 2152


Research Article

ACS Catalysis

Figure 10. Three plausible channels for the substrate delivery predicted by RAMD MD simulations.

between the closed and open states. The HCN molecule approaches Trp128, and when RC8 = 12 Å, their distance is the closest, which is around 4 Å (see in Figure S9a in the Supporting Information). This means that a noncovalent interaction may exist. At this moment, the indole ring of Trp128 evolves to an open state. The noncovalent interactions between the aromatic π regions with HCN have been explored by previous electrostatic potential calculations,77 which are responsible for formation of the counterintuitive complexes. Here our QM calculations also show the dependence of the HCN orientation on such noncovalent interactions (see Figure S9b), in which the most favorable conformation between HCN and the ring of Trp128 is at 34° and the distance distribution between related atoms from MD simulations (refer to Figure S9a) reveals that the strongest attraction arises from π···π and C−H···π interactions. We note that a lone pair···π interaction at 79° is negligible. In the second stage (6 Å ≤ RC7 < 8 Å, 14 Å ≤ RC8 < 16 Å), HCN passes through the “door” with a barrier of ∼2 kcal/mol, and the conformation of Trp128 changes from the open state to the closed state gradually. In the third stage (6 Å ≤ RC7 < 8 Å, RC8 ≥ 16 Å), the hydrocyanic acid is released from the protein, and the system arrives at a local minimum point energetically. Herein, Trp128 is closed entirely, and the acetone still stays at the active site. In the fourth stage (8 Å ≤ RC7 < 13 Å, RC8 ≥ 16 Å), the acetone is struggling to escape from the active site. Through its move, −CH3···π, lone pair···π, and π···π interactions occur between acetone and the indole ring of Trp128, which may attract the acetone move and induce the conformational change of Trp128 from the closed state to an open state. In the fifth stage (13 Å ≤ RC7 < 15 Å, RC8 ≥ 16 Å), the acetone is still moving toward the outside of the protein, and Trp128 gradually shuts down the channel with the aid of the −CH3···π interaction between acetone and the indole ring. In the final stage (RC7 ≥ 15 Å, RC8 ≥ 16 Å), the acetone is almost exposed to the water environment, and

Table 1. Number of Trajectories and Corresponding Possibility of Paths A−C for the Delivery of Substrate and Product on the Basis of RAMD MD Simulations reactant binding path A path B path C

acetone release

HCN release

no.

possibility, %

no.

possibility, %

no.

possibility, %

23 4 1

82.1 14.3 3.6

25 1 2

89.3 3.6 7.1

24 2 2

85.7 7.1 7.1

Figure 11. Free energy profile for the substrate transportation along the RC6 (the distance between CE of Lys237 and C1 of the reactant).

The product escape from the protein can be divided into six stages on the basis of the free energy profile. In the first stage (6 Å ≤ RC7 < 8 Å, 5 Å ≤ RC8 < 14 Å), the acetone binds at the active site, while HCN leaves from the binding site and moves toward to the “door” of the channel. The Trp128 sways 2153


Research Article

ACS Catalysis

Figure 12. Key residues and their conformational changes at the representative states in the main channel.

which provides the driving force for the following chemical reaction. The hydrophobic pocket accommodates the achiral acetone cyanohydrin, which is also responsible for the chiral substrate binding. The second stage of the chemical step dominates the efficiency of the entire enzymatic catalysis. This rate-determining process is composed of three steps: (i) deprotonation of the substrate by the catalytic triad of Ser80His236-Asp208, (ii) cleavage of the carbon−carbon bond and formation of the cyanide anion (the rate-limiting step in the enzymatic catalysis by MeHNL), and (iii) formation of the hydrocyanic acid arising from proton transfer via His236 → Ser80 → CN−. In the third stage, the hydrocyanic acid and acetone are released in proper sequence, on the basis of the two-dimensional free energy profile. The release of HCN proceeds first, followed by the release of acetone, and their corresponding free energy barriers are ∼2 and ∼10 kcal/mol, respectively. Accordingly, the release of HCN is quite facile, while the release of acetone may influence the catalytic efficiency. Two key residues of Trp128 and Lys237 are found to be important for the enzymatic catalysis by both HNLs. The residue Trp128 serves as the gating switch to regulate the substrate binding and product release, in which the π effects between substrate (or product) and Trp128 can manipulate the substrate movement. The key residue Lys237 plays a crucial role in the C−C bond breaking and the cyanide anion stabilization, whereas it was assumed to be trivial for MeHNL experimentally, on the basis of the X-ray structure of the Ser80Ala mutated enzyme. Moreover, our QM/MM MD simulations show that the protonation state of Lys237 will influence both the reaction mechanism and the catalytic efficiency. The HNL enzyme with the protonated Lys237 is more favorable than that with the unprotonated species, both thermodynamically and dynamically, and this is probably a general character for HNLs of the α/β-hydrolase superfamily with S stereoselectivity. MeHNL has been predicted to follow a mechanism similar to that of HbHNL for the chemical step, whereas the rate-limiting steps are quite different. For HbHNL, the HCN formation

Figure 13. Two-dimensional free energy profiles for the product release along RC7 (the distance between the CD of Lys237 and C1 of the acetone) and RC8 (the distance between the CD of Lys237 and C of the HCN).

nothing can block it; meanwhile the door of the channel is closed (see Figure S10 in the Supporting Information). We assume that the products have released completely. The calculated results agree well with the experimental speculation13 about the release order of products. To directly acquire the energy information on the whole enzymatic process, Figure 15 summarizes the relative free energy profiles for the entire process, including the chemical step and substrate and product transportation, which is useful for a better comprehension of the enzymatic catalysis of α/β-hydrolase superfamily with S stereoselectivity for HNLs.

4. CONCLUSIONS Extensive QM(B3LYP/6-31G*)/MM MD and MM MD simulations have been used to explore the full picture of the enzymatic process by MeHNL. Here the whole catalysis is composed of three major stages, including substrate binding, chemical reaction, and product release. The initial substrate binding is remarkably exothermic with a low free energy barrier, 2154


Research Article

ACS Catalysis

clarified in detail (see Figure 12), which provides a basis for understanding the functional roles of these residues in the biological synthesis of chiral molecules through the inverse reaction. Interestingly, the release of acetone and hydrocyanic acid in MeHNL follows a stepwise mechanism, which may be appropriate for HbHNL and other members of the α/β-hydrolase superfamily with S stereoselectivity. These new findings give guidance regarding site-directed mutagenesis for further experimental studies and the biosynthesis of new higher efficiency and better binding kinetic behavior pharmaceuticals.

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acscatal.5b02855. Building process of initial structure for RAMD MD simulations, root-mean-square deviation (RMSD) results of the backbone, changes in distance between −NH3+ of Lys237 and −CN of the substrate, titration curve of active-site lysine for HbHNL and MeHNL by the H++ program, relative energy profile of the proton transfer and hydrocyanic acid formation for model A by using a QM/MM scan, relative energy profile of the HNC formation that involves proton transfer between Lys237 and CN− for model A by using aQM/MM scan, defined reaction coordinate for delivery of acetone cyanohydrins, defined reaction coordinate for release of acetone and HCN, distance between atoms of HCN molecules and the aromatic ring of Trp128 during the last 1 ns in MD simulations for RC7 = 8 Å and RC8 = 12 Å, QM optimized structures based on different conformations between HCN and the aromatic ring of Trp128 which are extracted from snapshots and corresponding interaction energy, side view of the third and fifth stages of product release (PDF)

■

Figure 14. Representative states of the product release in the main channel from classical MD simulations with umbrella sampling.

AUTHOR INFORMATION

Corresponding Authors

*E-mail for C.W.: [email protected]. *E-mail for Z.C.: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work has been supported by the Ministry of Science and Technology (2011CB808504 and 2012CB214900) and the National Science Foundation of China (21133007, 21373164, and 21172053). We also thank Prof. Ruibo Wu at Sun Yat-sen University, Prof. Yirong Mo at Western Michigan University, and Dr. Shenglong Wang at New York University for their assistance in computation. We thank the National Supercomputing Centers in Hunan and Shenzhen for providing the computational resources.

Figure 15. Relative free energy profile for the entire enzymatic catalytic process.

dominates the whole reaction rate, while for MeHNL the carbon−carbon bond cleavage is the rate-determining step. Such a discrepancy may arise from the different strengths of the hydrogen bond network around acetone cyanohydrin and different proton sources for HCN formation. The free energy barrier for the chemical step of MeHNL is lower than that of HbHNL, which means that MeHNL probably has higher catalytic efficiency. Herein the effect of key residues around the channel on the delivery of substrate and product has been

■

REFERENCES

(1) Conn, E. E. The biochemistry of plants. A comprehensive treatise. Secondary plant products; Academic Press: New York, 1981; Vol. 7; pp 479−500. (2) Hickel, A.; Hasslacher, M.; Griengl, H. Physiol. Plant. 1996, 98, 891−898. (3) Wajant, H.; Effenberger, F. Biol. Chem. 1996, 377, 611−617. 2155


Research Article

ACS Catalysis (4) Sharma, M.; Sharma, N. N.; Bhalla, T. C. Enzyme Microb. Technol. 2005, 37, 279−294. (5) Purkarthofer, T.; Skranc, W.; Schuster, C.; Griengl, H. Appl. Microbiol. Biotechnol. 2007, 76, 309−320. (6) Andexer, J.; Langermann, J.; Kragl, U.; Pohl, M. Trends Biotechnol. 2009, 27, 599−607. (7) Dadashipour, M.; Asano, Y. ACS Catal. 2011, 1, 1121−1149. (8) Paravidino, M.; Sorgedrager, M. J.; Orru, R. V.; Hanefeld, U. Chem. - Eur. J. 2010, 16, 7596−7604. (9) Oroz-Guinea, I.; García-Junceda, E. Curr. Opin. Chem. Biol. 2013, 17, 236−249. (10) Müller, M. ChemBioEng Rev. 2014, 1, 14−26. (11) van Rantwijk, F.; Stolz, A. J. Mol. Catal. B: Enzym. 2015, 114, 25−30. (12) Kawahara, N.; Asano, Y. ChemBioChem 2015, 16, 1891−1895. (13) Wagner, U.; Hasslacher, M.; Griengl, H.; Schwab, H.; Kratky, C. Structure 1996, 4, 811−822. (14) Muschiol, J.; Peters, C.; Oberleitner, N.; Mihovilovic, M. D.; Bornscheuer, U. T.; Rudroff, F. Chem. Commun. 2015, 51, 5798−5811. (15) Fesko, K.; Gruber-Khadjawi, M. ChemCatChem 2013, 5, 1248− 1272. (16) Diebler, J.; von Langermann, J.; Mell, A.; Hein, M.; Langer, P.; Kragl, U. ChemCatChem 2014, 6, 987−991. (17) Kassim, M. A.; Rumbold, K. Biotechnol. Lett. 2014, 36, 223−228. (18) Lauble, H.; Miehlich, B.; Förster, S.; Wajant, H.; Effenberger, F. Protein Sci. 2001, 10, 1015−1022. (19) Lauble, H.; Förster, S.; Miehlich, B.; Wajant, H.; Effenberger, F. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2001, 57, 194−200. (20) Wajant, H.; Pfizenmaier, K. J. Biol. Chem. 1996, 271, 25830− 25834. (21) Gruber, K.; Gartler, G.; Krammer, B.; Schwab, H.; Kratky, C. J. Biol. Chem. 2004, 279, 20501−20510. (22) Chmura, A.; van der Kraan, G. M.; Kielar, F.; van Langen, L. M.; van Rantwijk, F.; Sheldon, R. A. Adv. Synth. Catal. 2006, 348, 1655− 1661. (23) Bühler, H.; Effenberger, F.; Förster, S.; Roos, J.; Wajant, H. ChemBioChem 2003, 4, 211−216. (24) Lauble, H.; Miehlich, B.; Förster, S.; Kobler, C.; Wajant, H.; Effenberger, F. Protein Sci. 2002, 11, 65−71. (25) Bühler, H.; Miehlich, B.; Effenberger, F. ChemBioChem 2005, 6, 711−717. (26) Gruber, K.; Gugganig, M.; Wagner, U. G.; Kratky, C. Biol. Chem. 1999, 380, 993−1000. (27) Gartler, G.; Kratky, C.; Gruber, K. J. Biotechnol. 2007, 129, 87− 97. (28) Mueller-Dieckmann, C.; Panjikar, S.; Schmidt, A.; Mueller, S.; Kuper, J.; Geerlof, A.; Wilmanns, M.; Singh, R. K.; Tucker, P. A.; Weiss, M. S. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2007, 63, 366− 380. (29) Zuegg, J.; Gruber, K.; Gugganig, M.; Wagner, U.; Kratky, C. Protein Sci. 1999, 8, 1990−2000. (30) Schmidt, A.; Gruber, K.; Kratky, C.; Lamzin, V. S. J. Biol. Chem. 2008, 283, 21827−21836. (31) Zhao, Y.; Chen, N.; Mo, Y.; Cao, Z. Phys. Chem. Chem. Phys. 2014, 16, 26864−26875. (32) Cui, F.-C.; Pan, X.-L.; Liu, J.-Y. J. Phys. Chem. B 2010, 114, 9622−9628. (33) Dreveny, I.; Andryushkova, A. S.; Glieder, A.; Gruber, K.; Kratky, C. Biochemistry 2009, 48, 3370−3377. (34) Anandakrishnan, R.; Aguilar, B.; Onufriev, A. V. Nucleic Acids Res. 2012, 40, W537−W541. (35) Myers, J.; Grothaus, G.; Narayanan, S.; Onufriev, A. Proteins: Struct., Funct., Genet. 2006, 63, 928−938. (36) Gordon, J. C.; Myers, J. B.; Folta, T.; Shoja, V.; Heath, L. S.; Onufriev, A. Nucleic Acids Res. 2005, 33, W368−W371. (37) Bashford, D.; Karplus, M. Biochemistry 1990, 29, 10219−10225. (38) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157−1174.

(39) 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. 1995, 117, 5179−5197. (40) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Proteins: Struct., Funct., Genet. 2006, 65, 712−725. (41) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269−10280. (42) Frisch, M.; Trucks, G.; Schlegel, H.; Scuseria, G.; Robb, M.; Cheeseman, J.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. Gaussian 09; Gaussian Inc., Wallingford, CT, 2010. (43) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926−935. (44) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089−10092. (45) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577−8593. (46) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. J. Comput. Phys. 1977, 23, 327−341. (47) Miller, B. R., III; McGee, T. D., Jr; Swails, J. M.; Homeyer, N.; Gohlke, H.; Roitberg, A. E. J. Chem. Theory Comput. 2012, 8, 3314− 3321. (48) Hou, T.; Wang, J.; Li, Y.; Wang, W. J. Comput. Chem. 2011, 32, 866−877. (49) Hou, T.; Wang, J.; Li, Y.; Wang, W. J. Chem. Inf. Model. 2011, 51, 69−82. (50) Wang, J.; Hou, T.; Xu, X. Curr. Comput.-Aided Drug Des. 2006, 2, 287−306. (51) Zoete, V.; Meuwly, M.; Karplus, M. Proteins: Struct., Funct., Genet. 2005, 61, 79−93. (52) Case, D.; Darden, T.; Cheatham III, T. E.; Simmerling, C.; Wang, J.; Duke, R.; Luo, R.; Walker, R.; Zhang, W.; Merz, K. AMBER 12; University of California, San Francisco, CA, 2012. (53) Wu, R.; Lu, Z.; Cao, Z.; Zhang, Y. J. Am. Chem. Soc. 2011, 133, 6110−6113. (54) Wu, R.; Wang, S.; Zhou, N.; Cao, Z.; Zhang, Y. J. Am. Chem. Soc. 2010, 132, 9471−9479. (55) Wu, R.; Xie, H.; Cao, Z.; Mo, Y. J. Am. Chem. Soc. 2008, 130, 7022−7031. (56) Shi, Y.; Zhou, Y.; Wang, S.; Zhang, Y. J. Phys. Chem. Lett. 2013, 4, 491−495. (57) Ke, Z.; Smith, G. K.; Zhang, Y.; Guo, H. J. Am. Chem. Soc. 2011, 133, 11103−11105. (58) Rooklin, D. W.; Lu, M.; Zhang, Y. J. Am. Chem. Soc. 2012, 134, 15595−15603. (59) Chen, N.; Zhao, Y.; Lu, J.; Wu, R.; Cao, Z. J. Chem. Theory Comput. 2015, 11, 3180−3188. (60) Zhao, Y.; Chen, N.; Wu, R.; Cao, Z. Phys. Chem. Chem. Phys. 2014, 16, 18406−18417. (61) Zhang, Y.; Lee, T.-S.; Yang, W. J. Chem. Phys. 1999, 110, 46−54. (62) Chen, X.; Zhang, Y.; Zhang, J. Z. J. Chem. Phys. 2005, 122, 184105. (63) Zhang, Y. Theor. Chem. Acc. 2006, 116, 43−50. (64) Zhang, Y. J. Chem. Phys. 2005, 122, 024114. (65) Zhang, Y.; Liu, H.; Yang, W. J. Chem. Phys. 2000, 112, 3483− 3492. (66) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187−199. (67) Beeman, D. J. Comput. Phys. 1976, 20, 130−139. (68) Berendsen, H. J.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. J. Chem. Phys. 1984, 81, 3684−3690. (69) Souaille, M.; Roux, B. t. Comput. Phys. Commun. 2001, 135, 40− 57. (70) Kumar, S.; Rosenberg, J. M. J. Comput. Chem. 1992, 13, 1011− 1021. (71) Shao, Y.; Molnar, L. F.; Jung, Y.; Kussmann, J.; Ochsenfeld, C.; Brown, S. T.; Gilbert, A. T.; Slipchenko, L. V.; Levchenko, S. V.; O’ Neill, D. P. Phys. Chem. Chem. Phys. 2006, 8, 3172−3191. (72) Ponder, J. W. TINKER, Software Tools for Molecular Design, Version 4.2; Washington University School of Medicine, Saint Louis, MO, 2004. 2156


Research Article

ACS Catalysis (73) Lüdemann, S. K.; Lounnas, V.; Wade, R. C. J. Mol. Biol. 2000, 303, 797−811. (74) Vashisth, H.; Abrams, C. F. Biophys. J. 2008, 95, 4193−4204. (75) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. J. Comput. Chem. 2005, 26, 1781−1802. (76) Yao, L.; Li, Y.; Wu, Y.; Liu, A.; Yan, H. Biochemistry 2005, 44, 5940−5947. (77) Murray, J. S.; Shields, Z. P.-I.; Seybold, P. G.; Politzer, P. J. Comput. Sci. 2015, 10, 209−216.

2157


A Comprehensive Understanding of Enzymatic Catalysis by

Recommend Documents