A full picture of enzymatic catalysis by hydroxynitrile lyases from Hevea brasiliensis: protonation dependent reaction steps and residue-gated movement of the substrate and the product

Abstract

Hydroxynitrile lyases (HNLs) defend plants from herbivores and microbial attack by releasing cyanide from hydroxynitriles. The reverse process has been productively applied to bioorganic syntheses of pharmaceuticals and agrochemicals. To improve our understanding of the catalytic mechanism of HNLs, extensive ab initio QM/MM and classical MM molecular dynamics simulations have been performed to explore the catalytic conversion of cyanohydrins into aldehyde (or ketone) and HCN by hydroxynitrile lyases from Hevea brasiliensis (HbHNLs). It was found that the catalytic reaction approximately follows a two-stage mechanism. The first stage involves two fast processes including the proton abstraction of the substrate through a double-proton transfer and the C–CN bond cleavage, while the second stage concerns HCN formation and is rate-determining. The complete free energy profile exhibits a peak of ∼18 kcal mol⁻¹. Interestingly, the protonation state of Lys236 influences the efficiency of the enzyme only to some extent, but it changes the entire catalytic mechanism. The dynamical behaviors of substrate delivery and HCN release are basically modulated by the gate movement of Trp128. The remarkable exothermicity of substrate binding and the facile release of HCN may drive the enzyme-catalyzed reaction to proceed along the substrate decomposition efficiently. Computational mutagenesis reveals the key residues which play an important role in the catalytic process.

---

1 Introduction

Hydroxynitrile lyases (HNLs) are important members of the α/β-hydrolase superfamily, and they catalyze the cleavage of cyanohydrins into aldehyde (or ketone) and hydrocyanic acid (HCN). The release of HCN not only protects plant systems from herbivores and microbial attack^1–3 but also provides a nitrogen source for the biosynthesis of asparagine.^4,5 The reverse in vivo reactions may also occur under different conditions, with the efficient enantioselectivity for the synthesis of chiral compounds.^6–12 In recent years, HNLs have been successfully utilized in the production of pharmaceuticals and agrochemicals owing to their importance in biocatalytic retrosynthesis.^13–20 Certainly, an improved understanding of the enzymatic catalysis mechanism of HNLs can potentially further optimize the production processes in industry and help the rational design of biocatalysts.

Herein, the hydroxynitrile lyase from Hevea brasiliensis (HbHNL) has been considered. This enzyme can catalyze the formation of carbon–carbon bonds with high (S)-stereoselectivity during the chiral molecule synthesis.^19,22 Up to now, a number of HbHNL crystal structures, including the apo state and the complex state, have been determined (see Table S1 in the ESI†).^21,23–28 These static structures provide an opportunity for subsequent computational simulations which can reveal the dynamical conformational evolution along the enzyme-catalyzed reactions.

Based on the crystal structures of the wild-type enzyme and its complexes with the inhibitor or the substrate in combination with the activity analyses of site-directed mutants, a general acid/base catalysis has been experimentally proposed as the most probable mechanism (see Fig. 1),²¹ where Ser80, His235 and Asp207 serve as the catalytic triad to initiate the reaction. It should be noted that all members of the α/β hydrolases contain a conserved catalytic triad (nucleophile–histidine–aspartate),²⁹ and in HNLs the nucleophile is serine. In consideration of the enzymatic environment, it was hypothesized that the deprotonation and the C–C bond cleavage of the substrate occur via a concerted mechanism, followed by HCN formation and release. We note that previous density functional calculations with an active-site model proposed a three-step mechanism (where the deprotonation and the C–C bond cleavage occur separately; see Fig. 1) and the predicted relative energies of the reaction severely depend on the choice of cluster models.³⁰

Fig. 1 The proposed catalytic mechanisms by experiment²¹ and by the QM calculation.³⁰

The kinetic characterization of the enzyme mutants indicates that the protonation state of Lys236 plays a critical role in substrate binding, together with Thr11 and Ser80.²¹ Moreover, the positively charged –NH₃⁺ group of the protonated Lys236 can stabilize the nascent cyanide ion which acts as a better leaving group. The fact that the cyanide ion is a leaving group is very important in the recycling of the enzyme back to its initial state. Based on the determined X-ray crystal structures of HbHNL complexes with non-natural chiral substrates, Gartler et al.²⁵ proposed that Lys236 may be involved in determining the enantioselectivity of the enzyme. They also claimed that the binding modes of the chiral substrates were identical with acetone cyanohydrin. This means that chiral and achiral substrates may involve the very same mechanism. However, the detailed elucidation of the effects of key residues on both substrate binding and the catalytic process is still elusive.

It has been well known that the accessibility of the active site for substrate binding and product release makes important contributions to the overall enzymatic efficiency. The crystal structure of HbHNL indicates that the active site is deeply buried inside the enzyme and there is only one narrow channel for substrate delivery.²⁴ In the channel for the ligand transportation, the flexibility of the side-chain residue Trp128 can be expected to play a nontrivial role.²⁵ However, plausible mechanisms and dynamical properties of acetone cyanohydrin delivery and hydrocyanic acid release also remain unknown.

It is thus expected that computational studies can provide a microscopic description of the catalytic process and shed light on the exact roles of individual residues and the details of substrate delivery and product release. Our computational investigations will basically focus on the key issues related to the whole enzymatic catalysis accordingly. The questions to be addressed include: (i) what is the most likely catalytic mechanism for the enzymatic reaction in the protein environment? Are the proposed reaction steps of double-proton transfer, C–C bond cleavage, and HCN formation stepwise or concerted? And which is the rate-determining step? (ii) What is the role of the protonation state of Lys236 in these catalytic processes and how does this protonation state influence the overall catalytic reaction? (iii) What are the precise roles of the key residues in the whole catalytic process? (iv) What are the transportation mechanisms for the substrate acetone cyanohydrin and the product hydrocyanic acid? In an attempt to clarify these issues, we have conducted both classical molecular dynamics (MD) and combined quantum mechanics and molecular mechanics (QM/MM) MD simulations, which are expected to generate a full picture of the overall enzymatic catalysis and provide detailed information on enzyme engineering for biosynthesis of organic molecules.

2 Computational methods

Since the achiral and chiral substrates share the same binding mode and follow similar catalytic mechanisms, the initial computational model was built based on the X-ray crystal structures of HbHNLs complexed with the achiral acetone cyanohydrin (PDB: 1SC9).²¹ The protonation states of ionizable residues were determined at pH 7.5 via the program PROPKA 3.1^31–34 as well as the previous Poisson–Boltzmann calculations.²¹ The Amber99SB force field^35–37 and the TIP3P model³⁸ were employed for the protein and water molecules, respectively. As to the ligand, the force field parameters were obtained by using the Amber GAFF force field (GAFF),³⁹ and the charge parameters were determined by the restrained electrostatic potential (RESP) method⁴⁰ at the HF/6-31G(d) level using Gaussian 03 package.⁴¹ The whole system was solvated in a ∼82 × 83 × 75 Å cubic water box with 15 Å buffer solvents on each side, and neutralized by sodium ions. Then, multistep optimizations were performed to remove poor interatomic contacts. After that, the system was heated up gradually from 0 to 300 K over a period of 50 ps, and another 50 ps MD simulation was carried out for further relaxation of the system. Finally, 10 ns MD simulation for equilibration was performed. All the MD simulations were employed under the NVT ensemble by using the periodic boundary condition with 12 Å cutoff distances for van der Waals and electrostatic interaction calculations. The Langevin method was utilized to control the temperature at 300 K. The bonds involving hydrogen atoms were constrained by using the SHAKE scheme.^42,43 The last 3 ns trajectories were used for ligand–residue interaction decomposition and binding free energy analysis by employing the MM-GBSA method.^44,45 All the molecular dynamics simulations were performed using AMBER12 software.⁴⁶

Based on the MD simulations, two QM/MM models (A and B) were prepared by deleting the water molecules beyond a 24 Å radius from the sulfur atom of Cys81. With respect to Model A, the Lys236 residue took a protonated state, while for Model B, the Lys236 residue was neutral (unprotonated). As shown in Fig. 2, the QM subsystem includes the substrate and the side chains of Ser80, His235, Thr11, Asp207 and Lys236. The QM region was treated by the DFT method with the B3LYP functional^47–50 and the 6-31G(d) basis set, which has been successfully used in many enzymatic catalytic systems.^51–60 As for the MM subsystem, the Amber99SB force field^35–37 was employed as in the above classical MD simulations. The QM/MM boundary was handled by the improved pseudo-bond approach.^61–63 The spherical boundary condition was applied, and the atoms beyond 20 Å from the spherical center were fixed. The cutoff values of 18 and 12 Å were utilized for electrostatic and van der Waals interactions among MM atoms, respectively.

Fig. 2 QM/MM models for the HbHNL enzyme with a protonated Lys236 (Model A) and an unprotonated Lys236 (Model B).

The QM/MM optimization was carried out to derive a minimum energy path using the reaction coordinate driving (RCD) method,⁶⁴ and for the electronic structures along the reaction path, the MM region was equilibrated for 500 ps by MM MD simulations. The snapshots abstracted from these MM MD simulations served as the initial structures for subsequent QM(DFT)/MM MD simulations with the umbrella sampling. Each window was calculated for 20 ps with the time step of 1 fs by using the Beeman algorithm.⁶⁵ The Berendsen thermostat⁶⁶ was employed to maintain the system temperature at 300 K. Afterwards, free energy profiles were determined by the weighted histogram analysis method (WHAM)^67–70 with the probability distribution for each window of the last 15 ps umbrella sampling. All ab initio QM/MM calculations were performed using the modified Q-Chem 4.0⁷¹ and Tinker 4.2⁷² programs.

To figure out the plausible entry and exit channels for acetone cyanohydrin and hydrocyanic acid transportation, the combined random acceleration molecular dynamics (RAMD) and MD simulations (RAMD-MD)^73,74 have been carried out using the NAMD 2.9 software.⁷⁵ The Amber99SB force field^35–37 and GAFF have been utilized for the protein and ligand, respectively. For the complex of the enzyme with achiral acetone cyanohydrin, the initial structure was obtained from the MM MD simulations. For the complex of the enzyme with hydrocyanic acid, the initial structure was set up by the sphere model of the HbHNL-product state (containing the acetone molecule and hydrocyanic acid) from QM/MM MD simulations by removing the acetone molecule.

In RAMD simulations, an additional force with random orientation is added to the center of mass (or other pre-defined points) of the ligand to identify the possible pathways in the binding pocket for the ligand fleeing away from the protein in a computationally feasible time. Within a certain period of time, when the ligand moves away beyond the threshold distance, the direction is maintained. Otherwise, a new random direction will be chosen. However, once the ligand escapes from the initial position, it may move towards the wrong channels because the random force is generally higher than the resistance of the protein. This issue can be avoided by using the combined RAMD-MD approach. Upon the escape of the ligand from the initial position, the conventional MD simulations would be switched on and the equilibration sampling is recovered. Herein, the random accelerations of 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, and 0.15 kcal Å⁻¹ g⁻¹ are applied to the C1 atom of achiral acetone cyanohydrin and the C atom of HCN. 24 RAMD MD trajectories for each model have been explored. In what follows, one or two most favorable channels will be discussed in detail by mapping out their free energy profiles along dynamic pathways using the umbrella sampling technique. A series of biasing harmonic potentials along the defined reaction coordinate have been tested. Based on the most appropriate biasing harmonic potential, at least 8 ns MD simulations are performed for each window. The free energy profiles are generated by using WHAM^67–70 with the probability distribution for all windows.

3 Results and discussion

3.1 Catalytic mechanism of HbHNL

The catalytic mechanism would be explored based on Model A at first, in which Lys236 takes a protonated state as suggested from experiments for the real enzymatic environment. Herein, d_C1–C2 + d_O1–H1 (RC1) and d_H2–Nδ − d_C2–H2 (RC2) are chosen as reaction coordinates. Whereas RC2 obviously refers to the proton transfer from His235 to cyanide ion, the choice of RC1 is intended to examine whether the C–C breaking within acetone cyanohydrin and the proton transfer from acetone cyanohydrin to Ser80 are stepwise or concerted. As we can see from the relative free energy profiles (see Fig. 3a) and the structures of the corresponding key states (see Fig. 4), the whole reaction can be envisioned as a quasi-two-step reaction: the first is the substrate deprotonation mediated by Ser80 and His235, which occurs approximately synchronously with the subsequent C1–C2 bond cleavage (R → TS1 → IM1 → TS2 → IM2) owing to the presence of a metastable proton-transfer intermediate (IM1), and the second is the protonation of the cyanide anion by His235 to yield HCN (IM2 → TS3 → IM3).

Fig. 3 Predicted free energy profiles of the whole catalytic reaction by ab initio QM/MM MD simulations for Model A (a) and Model B (b).

Fig. 4 The structures of the reactant, transition states, and intermediates involved in the catalytic reaction for Model A (the average distances for selected bonds in Å) from the QM/MM-MD sampling.

During the process of transition from the reactant state R to IM1, the H1 atom of the substrate is abstracted by the O2 atom of Ser80, which is concomitant with the proton transfer from Ser80 to the Nδ-position of His235. In other words, the two proton transfers are concerted. At the same time, the C1–C2 bond is stretched slightly, and the single C1–O1 bond has a notable tendency to form a polarized C1^δ+–O1^δ− bond (see the population analyses in Table S2, ESI†). Thus, IM1 mostly corresponds to the structure with proton transfers via substrate → Ser80 → His235, but the C1–C2 bond retained. In the following reaction step from IM1 to IM2, the double protons are completely transferred together with the C1–C2 bond cleavage. Simultaneously, the single C1–O1 bond is shortened from 1.36 ± 0.02 Å to 1.23 ± 0.01 Å and thus it is converted to a typical carbonyl group. Moreover, the Mulliken charge distributions compiled in Table S2 (ESI†) show that the negative charges on the C2 atom remarkably increase from 0.21 ± 0.09 to −0.53 ± 0.13 a.u., suggesting that the cyanide anion intermediate is formed.

At this first stage, the His235–Asp207 dyad accommodates one proton and its charged configuration is changed from His–Asp⁻ to His⁺–Asp⁻. Clearly, the hydrogen bond between them becomes stronger due to the favorable electrostatic attraction, implying that the two residues play a crucial role in stabilizing the nascent intermediate in catalysis. Here they serve as a general base to accept the proton from Ser80 to facilitate the deprotonation of the substrate by Ser80. Most noteworthy is the distance between the –NH₃⁺ group of Lys236 and the N1 atom, which is significantly shortened from 3.04 ± 0.33 to 1.75 ± 0.13 Å. Such strong hydrogen bond interactions can stabilize the labile cyanide anion, showing that the protonated Lys236 can assist the cleavage of the C–C bond through strong electrostatic (field) interactions as the driving force. The energy barriers for double-proton transfer and the C–C bond cleavage are 6.7 and 6.6 kcal mol⁻¹, respectively, suggesting that both processes are extremely fast and experimentally difficult to differentiate. This may be the reason why experiments support the hypothesis that the deprotonation and the C–C bond cleavage of the substrate are concerted. Overall, for the first stage from the reactant state to the stable IM2, the free energy change ΔG is 4.8 kcal mol⁻¹.

At the second stage, the cyanide anion gradually approaches His235 to form cyanic acid. We note that the charge on the C2 atom remarkably changes from −0.53 ± 0.13 to 0.35 ± 0.09 a.u., while the charge on the Nδ atom goes towards the opposite direction, lending support to the proton transfer from His235 to the cyanide anion. Meanwhile, the distance between the –NH₃⁺ group of Lys236 and CN⁻ is stretched from 1.75 ± 0.13 to 3.05 ± 0.59 Å, and the hydrogen bond between Ser80 and His235 dissolves. The hydrogen bond distance between His235 and Asp207 slightly increases from 1.62 ± 0.12 to 1.83 ± 0.12 Å. All these structural changes point to the formation of HCN at this stage.

What is noteworthy is the relatively high free energy barrier (13.1 kcal mol⁻¹) for HCN formation, indicating that this step in the protein environment is rate-determining and crucial for the overall enzymatic catalysis. This finding is quite different from the previous study without considering the protein environment and its dynamics,³⁰ which claimed the cleavage of the C–C bond as the rate-limiting step. The comparison of the equilibrium configurations of IM2 and IM3 suggests that the loss of strong ion-pair bonding interaction of (CN)⁻⋯(NH₃)⁺–Lys236 and the N–H bond cleavage in the process of HCN formation cannot be sufficiently compensated by the H–C bond formation, resulting in a relatively high energy barrier and a metastable intermediate configuration IM3. Furthermore, the remarkable configuration evolution and the environmental effect of proteins during the proton abstraction of His235 by CN⁻ are also responsible for these thermodynamic and dynamical properties to a certain extent.

Knowing that the rate-determining step has a barrier of 13.1 kcal mol⁻¹, and the two earlier steps have barriers of 6.7 and 6.6 kcal mol⁻¹, we anticipate that the complex reaction barrier would be very close to the value of 16.1 kcal mol⁻¹ estimated by the transition state theory from the experimental data (k_cat = 10 s⁻¹ at 303 K. Though the kinetic parameters were measured with a partially purified enzyme, they can still provide useful information about the kinetic properties of the enzyme).⁷⁶ It should be pointed out that the IM3 intermediate is very unstable compared to IM2 and thus the forward catalytic reaction (rate-determining step) seems to be quite unfavorable. How does then the reaction take place? Our further simulations (see the following) demonstrate that the facile release of the newly formed HCN plays a key role in facilitating this forward reaction, and the remarkable energy release in the substrate binding may further compensate this energy requirement (vide infra).

3.2 Effect of the protonation state of Lys236 on catalysis

To evaluate the role of Lys236 with different protonation states in catalysis, we performed computational studies similar to the above to probe the possible catalytic mechanism based on Model B where Lys236 is unprotonated and neutral in HbHNL. In contrast to Model A, the energy scanning calculations along RC1 (d_C1–C2 + d_O1–H1) were unsuccessful with very high energy barriers. Alternatively, new reaction coordinates of d_O1–H1 − d_O2–H1 (RC3) and d_H2–N − d_O2–H2 (RC4) were defined to follow the proton transfer. Fig. 3b shows the energy profile with Model B, where the largest energy barrier is 15.1 kcal mol⁻¹, higher than the rate-limiting step with Model A (13.1 kcal mol⁻¹). The present results suggest that, while a neutral Lys236 does slow down the enzyme-catalyzed decomposition of cyanohydrins, the process is still efficient enough as the overall complex barrier is less than 22.4 kcal mol⁻¹. What we see from the comparison between Fig. 3a and b is the dramatic difference in the reaction mechanisms.

As Fig. 3b and Fig. S1 in ESI† show, the reaction mechanism with Model B essentially comprises only two steps: one is the deprotonation of the substrate through the double-proton transfer mediated by Ser80 and His235 (R → TS1 → IM1), and the other is the concurrent carbon–carbon bond cleavage and the HCN formation (IM1 → TS platform → IM2). This mechanism is totally different from Model A. Here no C–C bond cleavage intermediate with the CN⁻ moiety (IM2 in Fig. 3a) was identified and a quite flat plateau of transition state appears in the relative free energy profile. Therefore, the protonation state of Lys236 controls the catalytic mechanism.

In Model A, the –NH₃⁺ group of the protonated Lys236 significantly stabilizes the nascent CN⁻ through strong ion-pair bonding interactions and thereby can drive the C–C bond cleavage, while there are no such strong electrostatic interactions in Model B with the unprotonated Lys236. Without the electrostatic field from the protonated Lys236, the intermediate configuration with a CN⁻ species is unlikely formed in Model B. The structures of possible reaction configurations are presented in Fig. S1 and S2 in ESI.†

We note that here the C–C bond cleavage and the HCN formation require a long journey in the active domain, which may account for the occurrence of the flat transition-state plateau in Model B. As a consequence, in Model B transition-state configurations may be very different during our sampling simulations, though they have comparable energies. Here two kinds of reaction tendencies (see Fig. S2 in ESI†) have been discussed. One may proceed towards the reactant-analogue state while the other shows a tendency towards the final HCN formation, where the hydrogen atom in the newly formed HCN may come from either of the two transferred protons.

3.3 The role of key residues in the catalytic process

Experimentally,²¹ it has been assumed that the residues Ser80–His235–Asp207 form a catalytic triad as a general base to take part in the whole reaction. Among them, Ser80 and His235 mediate the double-proton transfer directly, and Asp207 retains the hydrogen bonding with the Nδ atom of His235 to keep the latter at a proper position towards the hydrogen abstraction and strengthen the basicity of Ser80–His235. Apart from the catalytic triad, other important residues include Lys236 and Thr11. As mentioned above, the protonation state of Lys236 plays a role in the catalytic process. The residue Thr11, however, is involved in the hydrogen bonding network around the substrate.

To further elucidate the roles of these key residues, QM(B3LYP)/MM calculations have been performed on the Asp207Ala, Lys236Ala, and Thr11Ala mutant systems. Based on the QM/MM optimized structures, QM/MM energy scans along the reaction coordinate of RC1 were performed. As we can see from Fig. 5, the barrier for the carbon–carbon bond cleavage in the Asp207Ala mutant system is higher than that of the wild system by 4.4 kcal mol⁻¹. Quite interestingly, the T11A mutant system exhibits a much lower barrier of ∼10 kcal mol⁻¹, compared to the wild HbHNL enzyme. However, the T11A mutation may destroy the hydrogen-bonding network around the substrate and make the substrate less likely for the initial proton transfer, resulting in a loss of enzymatic activity as observed experimentally.²¹ As for the Lys236Ala system, the C–C bond breaking cannot be achieved as it needs to couple other bond cleavage processes and thus experiences much high barriers which makes the reaction unrealistic under ambient conditions. In conclusion, mutations of the key residues of Asp207, Lys236, and Thr11 may have a remarkable effect on the C1–C2 bond cleavage.

Fig. 5 The predicted relative energies for the carbon–carbon cleavage in HbHNL, ASP207ALA, and THR11ALA systems by the QM/MM scan.

The QM/MM-optimized configurations indicate that the hydrogen bond distances between N1 of the substrate and –NH₃⁺ of Lys236 are 2.49, 2.74, and 1.99 Å for wild-type HbHNL, mutant systems Asp207Ala and Thr11Ala, respectively. These results show that the lack of hydrogen bonds between the substrate and Thr11 may remarkably enhance the hydrogen bonds between the substrate and (NH₃⁺)–Lys236, which will facilitate the C–CN bond cleavage through strong electrostatic interactions, although this near-attack conformer (NAC) is unlikely to be accessible in such a situation.

3.4 The role of key residues in substrate binding

In previous studies, Gruber et al.²¹ proposed that the residues Lys236, Ser80, and Thr11 play an important role in substrate binding. To get a deeper insight into the individual contributions of these three residues to substrate binding, the binding free energies were calculated for the Lys236Ala, Ser80Ala, Thr11Ala, and original HbHNL systems. The major results are summarized in Fig. 6a and Table S3 (ESI†). We note that the van der Waals, electrostatic, and non-polar solvation interactions are the major driving forces for substrate binding, while the polar solvation energy and the entropy effect are unfavorable (i.e., positive) terms. Although strong electrostatic interactions approximately counteract the polar solvation energy here, they should be important for the C–CN bond activation and breaking. In comparison with the wild HbHNL system, the binding free energies of the Lys236Ala, Ser80Ala, and Thr11Ala mutant systems decrease by 2.96, 4.74, and 1.23 kcal mol⁻¹, respectively, which mainly arise from the electrostatic interactions. On the basis of per-residue type, the free energy decomposition shows that the residues Ile12, Cys81, Leu157, and His235 in the active site also have notable contributions to substrate binding, as shown in Fig. 6b, in which the van der Waals and non-polar solvation energies play a primary role in binding the substrate.

Fig. 6 (a) The binding free energy components for HbHNL, Lys236Ala, Ser80Ala, and Thr11Ala systems. (b) The free energy decomposition on the basis of per-residue type along with the binding mode. (c) The energy decomposition for the key residues.

3.5 The acetone cyanohydrin access to the active site

Most of the studies on enzymatic catalysis focus on chemical reaction processes. However, it is equally important to investigate the delivery processes of the substrate and the product particularly when the active site of an enzyme is buried inside and both the substrate and the product need to go through a channel. To have an insight into the thermodynamic and dynamical properties of HbHNL to accommodate the substrate, we explored the possible pathways for substrate delivery to the active domain. However, in view of the computational efficiency and convenience, its reverse process, i.e., the release of acetone cyanohydrin, in the complexation of the enzyme with acetone cyanohydrin, has actually been used to characterize the substrate access to the active site of HbHNL.

Fig. 7 shows four possible pathways (Pa-1, Pa-2, Pb, and Pc) for acetone cyanohydrin release identified by RAMD-MD simulations, which are defined as Pa-1 (between helix D and Trp128), Pa-2 (between Trp128 and helix D1′), Pb (between helix D1′ and helix E), and Pc (between helix D2′ and D3′). 24 RAMD-MD trajectories are summarized in Table 1. We note that 17 trajectories among them follow the Pa-1 channel and Pa-1 is thus the predominant pathway for substrate release (or substrate docking). The possibilities of the other three channels are almost the same. For Pa-1, according to the escaping direction, the distance between the Cδ atom of Lys236 and the C1 atom of acetone cyanohydrin is chosen as the release coordinate (named RC5, see Fig. S3a, ESI†). MD simulations along this coordinate (RC5) from 6.5 to 16.0 Å with a biasing harmonic potential of 30 kcal mol⁻¹ were performed. The potential of the mean force (PMF) profile is depicted in Fig. 8a. Clearly, the relative free energy profiles for different sampling time durations are very similar, showing reliable convergence of MD simulations for the PMF profiles. The predicted relative free energies in Fig. 8a reveal that the access of the substrate to the active site is quite favorable thermodynamically with a negligible barrier (∼1 kcal mol⁻¹), and the binding free energies ΔG are about −15 kcal mol⁻¹. Such remarkable exothermicity may drive the subsequent catalytic process.

Fig. 7 The possible channels for acetone cyanohydrin delivery predicted by RAMD-MD simulations. The acetone cyanohydrin in the active site is shown as a sphere, and the red part of the surface in the upper right figure refers to Trp128.

Table 1 Statistics of 24 trajectories for substrate delivery and HCN release^a

Substrate	Share	Possibility (%)	HCN	Share	Possibility (%)
a Random accelerations of 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, and 0.15 kcal Å^−1 g^−1 are applied to the C1 atom of achiral acetone cyanohydrin and the C atom of hydrocyanic acid, respectively.
Pa-1	17	70.8	Pa	18	75.0
Pa-2	3	12.5	Pb	2	8.0
Pb	2	8.3	Pc	4	17.0
Pc	2	8.3

---

Fig. 8 (a) Free energy profiles for the release of acetone cyanohydrin along RC5 (the distance between the Cδ atom of Lys236 and the C1 atom of acetone cyanohydrin). (b) Free energy profiles for the release of HCN along RC7 (the distance between the CZ atom of Tyr158 and the C atom of HCN).

The PMF profiles for the Pa-2 channel have also been characterized by RAMD-MD simulations (see Fig. S4, ESI†). The distance between the CB atom of Thr11 and the C1 atom of the substrate has been taken as the release coordinate (named RC6, see Fig. S3a, ESI†). MD simulations suggest that the free energies of release ΔG are about 20 kcal mol⁻¹, and thus it should be facile for the transport of the substrate to the active site through the Pa-2 channel energetically. For the reverse process, i.e., the substrate release, such relatively high energy requirement may result in the occurrence of only a few trajectories for the Pa-2 channel during MD simulations. Presumably, the steric effects along the channels are also important for substrate delivery. In what follows, we will discuss the substrate transportation along the most probable channel Pa-1 in detail, in order to study the dynamical features and key residue effects on substrate movement.

We approximately classify the process of acetone cyanohydrin release into four stages based on the conformational evolutions along the reaction coordinate RC5 (Fig. 8a). At the first stage (7.0 Å ≤ RC < 9.0 Å), the channel is totally closed, and Trp128 blocks the release of acetone cyanohydrin. Thr11 and Ser80 may stabilize acetone cyanohydrin at the original active site via the hydrogen bonding interactions. At the second stage (9.0 Å ≤ RC ≤ 11.0 Å), the residue Trp128 turns on one side gradually as a switch, and the “door” (indole ring) swings between an opening state and a closing state. Acetone cyanohydrin tends to move away from the original active site, subsequently inducing a conformational change of helix D from an α-helix to a loop structure. At the third stage (11.0 Å < RC ≤ 14.0 Å), the system overcomes a barrier of ∼16.0 kcal mol⁻¹. Trp128 no longer blocks the ligand from moving out of the channel, as the pocket almost remains open, and helix D recovers its original α-helix structure gradually. Moreover, residues Tyr133, Trp128, and Gln180 provide the driving force for the release of cyanohydrin through the direct and stabilized hydrogen bonding, the σ–π interaction, and the hydrogen-bond network with adjacent water molecules, respectively. At the last stage (RC > 14.0 Å), the acetone cyanohydrin leaves the protein completely, and the Trp128 returns to the original “closed” state once again. The surfaces of the binding pocket from different windows are shown in Fig. 9a, and the secondary structures and key residues related to the release process are shown in Fig. 9b. The flexibility of Trp128 agrees well with the previous experimental study,²⁵ and similar dynamical behaviours have been found in other biosystems.⁷⁷

Fig. 9 (a) Surfaces of the binding pocket from different windows along RC5 (the distance between the Cδ atom of Lys236 and the C1 atom of the ligand) for the release of acetone cyanohydrin. (b) The key residues in the active site of HbHNL from different windows along RC5 (the distance between the Cδ atom of Lys236 and the C1 atom of the ligand) for the release of acetone cyanohydrin. The purplish red part represents helix D, and obviously its conformer changes along with the flip of Trp128. Trp128 is shown in stick representation and the reaction path corresponds to Pa-1 in Fig. 7.

It is worth mentioning that the HNL from Manihot esculenta (MeHNL) has a similar structure to HbHNL. In MeHNL, Trp128 was mutated to Ala, resulting in an improved conversion for larger aromatic cyanohydrins.⁷⁸ In addition, the residue also influences the stereoselectivity.⁷⁹ Despite these testified roles of Trp128 in the catalytic process in MeHNL, they are still required to be elucidated for HbHNL. Therefore, we mutated Trp128 to alanine in silico to quantify its function in the substrate release (or binding) process. We obtained the free energy profile for the Trp128Ala mutant system by the umbrella sampling approach (see Fig. S5, ESI†). The energy required for the substrate release in the Trp128Ala mutant is about 8 kcal mol⁻¹, only half of the value (∼16 kcal mol⁻¹) in the wild system. Clearly, the steric effect of the large indole ring and the conformational dynamics of Trp128 may dominate the substrate delivery, though other hydrogen-bond networks and hydrophobic residues are also involved in this process to some extent.

3.6 Release of HCN

Needless to say, the release of HCN is significant indicating the whole purpose of HNLs for the plant defense system and biosynthesis of asparagine. From the perspective of kinetics, the removal of HCN from the active site may also drive the catalytic reaction to move along the forward reaction. Thus, by employing the RAMD-MD method, we similarly explored its possible release channels. The distance between the CZ atom of Tyr128 and the C atom of HCN was chosen as the reaction coordinate (named RC7, see Fig. S3b, ESI†), and 23 windows were adopted from 4.0 Å to 15.0 Å. Our simulations revealed that HCN may be released through three pathways, including Pa (between helix D and Trp128), Pb (between helix D₁′ and helix E), and Pc (between helix D₂′ and D₃′), with possibilities of 75%, 8%, and 17%, respectively (see Table 1). Based on the umbrella sampling method, the PMF profiles for Pa have been obtained and shown in Fig. 8b. The predicted free energy barrier is only ∼4.0 kcal mol⁻¹, and thus HCN can easily escape from the protein. This is understandable due to its much small molecular size. The facile release of HCN can drive the catalytic reaction to proceed efficiently, though the product state is higher in energy than the reactant state in the lytic reaction of acetone cyanohydrin. As mentioned above, the remarkable energy release in the substrate binding may compensate for the energy requirement in the catalytic step.

The key residues and secondary structures from different windows are shown in Fig. S6 (ESI†). In the beginning (4.0 Å ≤ RC ≤ 8.0 Å), HCN moves dynamically among the residues His14, Ile12, Lys236, Thr11, Ser80, and His235, with a negligible barrier of about 0.5 kcal mol⁻¹. Afterwards (8.0 Å < RC ≤ 12.0 Å), the energy increases gradually, probably due to the flip of Trp128. As we can see from Fig. S6 (ESI†), HCN sways on both sides of the residue, which may cause the barrier of release to increase. At last (RC > 12.0 Å), HCN moves away from the enzyme completely. At this point, Trp128 returns to the closed state as well, and its gating mechanism plays a role in the HCN release, but overall its steric effect is limited due to the small size of HCN.

4 Conclusions

On the basis of extensive QM/MM MD and RAMD MD simulations, the catalytic mechanism of HbHNL and its substrate delivery and product (HCN) release have been fully explored. Approximately the chemical step in the low-energy enzymatic process involves two stages. The first stage includes the substrate deprotonation by Ser80 and His235 and the C–CN bond cleavage. These two processes are fast and often considered together. The second stage is the cyanic acid formation which is the rate-determining step with a barrier of 13.1 kcal mol⁻¹. Combining these three chemical steps together, we predict that the complex reaction barrier would be pretty much close to the experimental estimate of 16.1 kcal mol⁻¹. Calculations show that the catalytic mechanism strongly depends on the protonation state of Lys236, as the protonated Lys236 can facilitate the C–CN bond breaking and stabilize the nascent CN⁻, resulting in a relatively favorable mechanism dynamically. With a protonated Lys236, the rate-determining step appears to be the HCN formation, owing to the loss of strong ion-pair bonding interactions of (CN)⁻⋯(NH₃)⁺–Lys236, the N–H bond cleavage, and a remarkable conformational change. In contrast, with a neutral Lys236, the rate-limiting step is the concurrent C–CN bond cleavage and the HCN formation. Energetically, however, hydroxynitrile lyases with either protonated or neutral Lys236 can efficiently catalyze the decomposition of cyanohydrins.

Based on RAMD MD simulations, plausible channels for substrate delivery and HCN release have been explored. In the most probable channels, the key residue Trp128 shows a dynamical gating mechanism, which may dominate the transportation of the substrate and HCN. The binding of acetone cyanohydrin is predicted to be remarkably exothermic, while the release of HCN requires an energy of only 4.0 kcal mol⁻¹. The roles of residues in catalysis and delivery of the substrate and HCN are also analyzed based on the MD simulations and free energy decomposition. The mutations of key residues can alter the enzymatic activity remarkably. In comparison with the achiral acetone cyanohydrin, a chiral substrate such as mandelonitrile most probably behaves differently. The aromatic residue may be bound to the hydrophobic pocket and thus has a preferred orientation. However, as was deduced from the X-ray structural analysis, the HbHNL complex with symmetric substrates and asymmetric substrates probably follows similar catalytic mechanisms, binding modes, and release channels.²⁵ Thus we expect that a thorough understanding of the whole enzymatic catalysis as presented in this work can provide meaningful clues for the biocatalytic retrosynthesis of chiral compounds such as mandelonitrile and so on.

Acknowledgements

This work was supported by the Ministry of Science and Technology (2011CB808504 and 2012CB214900), the National Science Foundation of China (21133007 and 21373164), and the US National Science Foundation under the grant CHE-1055310 and Western Michigan University (FRACAA). We thank the National Supercomputing Center in Shenzhen and Guangzhou for providing the computational resource. We also thank Prof. Ruibo Wu at Sun Yat-sen University, Dr Shenglong Wang at New York University and Dr Vlad Cojocaru at Max Planck Institute for Molecular Biomedicine for their kind help.

Notes and references

1. E. E. Conn, in The Biochemistry of Plants: A Comprehensive Treatise, ed. P. K. Stumpf and E. E. Conn, Academic Press, New York, 1981, vol. 7, pp. 479–500.
2. A. Hickel, M. Hasslacher and H. Griengl, Physiol. Plant., 1996, 98, 891–898.
3. H. Wajant and F. Effenberger, Biol. Chem., 1996, 377, 611–617.
4. R. Lieberei, D. Selmar and B. Biehl, Plant Syst. Evol., 1985, 150, 49–63.
5. D. Selmar, ACS Symp. Ser., 1993, 13, 191–204.
6. M. Gruber-Khadjawi, T. Purkarthofer, W. Skranc and H. Griengl, Adv. Synth. Catal., 2007, 349, 1445–1450.
7. J. von Langermann, J.-K. Guterl, M. Pohl, H. Wajant and U. Kragl, Bioprocess Biosyst. Eng., 2008, 31, 155–161.
8. M. Paravidino, M. J. Sorgedrager, V. A. Romano and U. Hanefeld, Chem. – Eur. J., 2010, 16, 7596–7604.
9. S. E. Milner, T. S. Moody and A. R. Maguire, Eur. J. Org. Chem., 2012, 3059–3067.
10. I. Oroz-Guinea and E. García-Junceda, Curr. Opin. Chem. Biol., 2013, 17, 236–249.
11. M. Müller, ChemBioEng Reviews, 2014, 1, 14–26.
12. J. Diebler, J. von Langermann, A. Mell, M. Hein, P. Langer and U. Kragl, ChemCatChem, 2014, 6, 987–991.
13. M. Fechter and H. Griengl, in Enzyme catalysis in organic synthesis, ed. K. Drauz and H. Waldmann, Wiley-VCH Verlag GmbH, Weinheim, 2nd edn, 2002, vol. II, pp. 974–989.
14. T. Purkarthofer, W. Skranc, C. Schuster and H. Griengl, Appl. Microbiol. Biotechnol., 2007, 76, 309–320.
15. J. Holt and U. Hanefeld, Curr. Org. Synth., 2009, 6, 15–37.
16. N. J. Turner and E. O'Reilly, Nat. Chem. Biol., 2013, 9, 285–288.
17. U. Hanefeld, Chem. Soc. Rev., 2013, 42, 6308–6321.
18. M. Sharma, N. N. Sharma and T. C. Bhalla, Enzyme Microb. Technol., 2005, 37, 279–294.
19. M. Dadashipour and Y. Asano, ACS Catal., 2011, 1, 1121–1149.
20. I. Dreveny, A. S. Andryushkova, A. Glieder, K. Gruber and C. Kratky, Biochemistry, 2009, 48, 3370–3377.
21. K. Gruber, G. Gartler, B. Krammer, H. Schwab and C. Kratky, J. Biol. Chem., 2004, 279, 20501–20510.
22. N. Andexer, N. Staunig, T. Eggert, C. Kratky, M. Pohl and K. Gruber, ChemBioChem, 2012, 13, 1932–1939.
23. K. Gruber, M. Gugganig, U. G. Wagner and C. Kratky, Biol. Chem., 1999, 380, 993–1000.
24. U. G. Wagner, M. Hasslacher, H. Griengl, H. Schwab and C. Kratky, Structure, 1996, 4, 811–822.
25. G. Gartler, C. Kratky and K. Gruber, J. Biotechnol., 2007, 129, 87–97.
26. C. Mueller-Dieckmann, S. Panjikar, A. Schmidt, S. Mueller, J. Kuper, A. Geerlof, M. Wilmanns, R. K. Singh, P. A. Tucker and M. S. Weiss, Acta Crystallogr., Sect. D: Biol. Crystallogr., 2007, 63, 366–380.
27. J. Zuegg, K. Gruber, M. Gugganig, U. G. Wagner and C. Kratky, Protein Sci., 1999, 8, 1990–2000.
28. A. Schmidt, K. Gruber, C. Kratky and V. S. Lamzin, J. Biol. Chem., 2008, 283, 21827–21836.
29. D. L. Ollis, E. Cheah, M. Cygler, B. Dijkstra, F. Frolow, S. M. Franken, M. Harel, S. J. Remington, I. Silman, J. Schrag, J. L. Sussman, K. H. G. Verschueren and A. Goldman, Protein Eng., 1992, 5, 197–211.
30. F. Cui, X. Pan and J. Liu, J. Phys. Chem. B, 2010, 114, 9622–9628.
31. H. Li, A. D. Robertson and J. H. Jensen, Proteins, 2005, 61, 704–721.
32. D. C. Bas, D. M. Rogers and J. H. Jensen, Proteins, 2008, 73, 765–783.
33. H. M. Olsson, C. R. Søndergard, M. Rostkowski and J. H. Jensen, J. Chem. Theory Comput., 2011, 7, 525–537.
34. C. R. Søndergaard, M. H. M. Olsson, M. Rostkowski and J. H. Jensen, J. Chem. Theory Comput., 2011, 7, 2284–2295.
35. W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell and P. A. Kollman, J. Am. Chem. Soc., 1995, 117, 5179–5197.
36. V. Hornak, R. Abel, A. Okur, B. Strockbine, A. Roitberg and C. Simmerling, Proteins, 2006, 65, 712–725.
37. J. Wang, P. Cieplak and P. A. Kollman, J. Comput. Chem., 2000, 21, 1049–1074.
38. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
39. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
40. D. A. Pearlman, D. A. Case, J. W. Caldwell, W. S. Ross, T. E. Cheatham, S. Debolt, D. Ferguson, G. Seibel and P. Kollman, Comput. Phys. Commun., 1995, 91, 1–41.
41. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09, Revision A.1, Gaussian Inc., Wallingford, CT, 2009.
42. J.-P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, J. Comput. Phys., 1977, 23, 327–341.
43. S. Miyamoto and P. A. Kollman, J. Comput. Phys., 1992, 13, 952–962.
44. A. Kollman, I. Massova, C. Reyes, B. Kuhn, S. Huo, L. Chong, M. Lee, T. Lee, Y. Duan, W. Wang, O. Donini, P. Cieplak, J. Srinivasan, D. A. Case and T. E. Cheatham, Acc. Chem. Res., 2000, 33, 889–897.
45. D. A. Case, T. A. Darden, T. E. Cheatham III, C. L. Simmerling, J. Wang, R. E. Duke, R. Luo, R. C. Walker, W. Zhang, K. M. Merz, B. Roberts, S. Hayik, A. Roitberg, G. Seabra, J. Swails, A. W. Goetz, I. Kolossváry, K. F. Wong, F. Paesani, J. Vanicek, R. M. Wolf, J. Liu, X. Wu, S. R. Brozell, T. Steinbrecher, H. Gohlke, Q. Cai, X. Ye, J. Wang, M. J. Hsieh, G. Cui, D. R. Roe, D. H. Mathews, M. G. Seetin, R. Salomon-Ferrer, C. Sagui, V. Babin, T. Luchko, S. Gusarov, A. Kovalenko and P. A. Kollman, AMBER 12, University of California, San Francisco, CA, 2012.
46. T. Hou, W. Zhang, D. A. Case and W. Wang, J. Mol. Biol., 2008, 376, 1201–1214.
47. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
48. C. T. Lee, W. T. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
49. S. H. Vosko, L. Wilk and M. Nusair, Can. J. Phys., 1980, 58, 1200–1211.
50. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627.
51. D. W. Rooklin, M. Lu and Y. Zhang, J. Am. Chem. Soc., 2012, 134, 15595–15603.
52. Z. Ke, G. K. Smith, Y. Zhang and H. Guo, J. Am. Chem. Soc., 2011, 133, 11103–11105.
53. R. Wu, Z. Lu, Z. Cao and Y. Zhang, J. Am. Chem. Soc., 2011, 133, 6110–6113.
54. R. Wu, S. Wang, N. Zhou, Z. Cao and Y. Zhang, J. Am. Chem. Soc., 2010, 132, 9471–9479.
55. R. Wu, H. Xie, Z. Cao and Y. Mo, J. Am. Chem. Soc., 2008, 130, 7022–7031.
56. Y. Zhou and Y. Zhang, Chem. Commun., 2011, 47, 1577–1579.
57. Y. Shi, Y. Zhou, S. Wang and Y. Zhang, J. Phys. Chem. Lett., 2013, 4, 491–495.
58. R. Wu, W. Gong, T. Liu, Y. Zhang and Z. Cao, J. Phys. Chem. B, 2012, 116, 1984–1991.
59. Z. Ke, H. Guo, D. Xie, S. Wang and Y. Zhang, J. Phys. Chem. B, 2011, 115, 3725–3733.
60. N. Chen, J. Zhou, J. Li, J. Xu and R. Wu, J. Chem. Theory Comput., 2014, 10, 1109–1120.
61. Y. Zhang, J. Chem. Phys., 2005, 122, 024114.
62. Y. Zhang, Theor. Chem. Acc., 2006, 116, 43–50.
63. Y. Zhang, T. S. Lee and W. Yang, J. Chem. Phys., 1999, 110, 46–54.
64. Y. Zhang, H. Liu and W. Yang, J. Chem. Phys., 2000, 112, 3483–3492.
65. D. Beeman, J. Comput. Phys., 1976, 20, 131–139.
66. H. J. C. Berendsen, J. P. M. Postma, W. F. Van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
67. S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman and J. M. Rosenberg, J. Comput. Chem., 1992, 13, 1011–1021.
68. M. Souaille and B. Roux, Comput. Phys. Commun., 2001, 135, 40–57.
69. A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett., 1988, 61, 2635–2638.
70. A. Grossfield, WHAM: the Weighted Histogram Analysis Method, Version 2.0.6, http://membrane.urmc.rochester.edu/content/wham.
71. Y. Shao, L. F. Molnar, Y. Jung, J. Kussmann, C. Ochsenfeld, S. T. Brown, A. T. Gilbert, L. V. Slipchenko, S. V. Levchenko, D. P. O'Neill, R. A. DiStasio Jr., R. C. Lochan, T. Wang, G. J. Beran, N. A. Besley, J. M. Herbert, C. Y. Lin, T. Van Voorhis, S. H. Chien, A. Sodt, R. P. Steele, V. A. Rassolov, P. E. Maslen, P. P. Korambath, R. D. Adamson, B. Austin, J. Baker, E. F. Byrd, H. Dachsel, R. J. Doerksen, A. Dreuw, B. D. Dunietz, A. D. Dutoi, T. R. Furlani, S. R. Gwaltney, A. Heyden, S. Hirata, C. P. Hsu, G. Kedziora, R. Z. Khalliulin, P. Klunzinger, A. M. Lee, M. S. Lee, W. Liang, I. Lotan, N. Nair, B. Peters, E. I. Proynov, P. A. Pieniazek, Y. M. Rhee, J. Ritchie, E. Rosta, C. D. Sherrill, A. C. Simmonett, J. E. Subotnik, H. L. Woodcock III, W. Zhang, A. T. Bell, A. K. Chakraborty, D. M. Chipman, F. J. Keil, A. Warshel, W. J. Hehre, H. F. Schaefer III, J. Kong, A. I. Krylov, P. M. Gill and M. Head-Gordon, Phys. Chem. Chem. Phys., 2006, 8, 3172–3191.
72. J. W. Ponder, TINKER, Software Tools for Molecular Design, Version 4.2, Washington University, Saint Louis, MO, 2004.
73. S. K. Lüdemann, V. Lounnas and R. C. Wade, J. Mol. Biol., 2000, 303, 797–811.
74. H. Vashisth and C. F. Abrams, Biophys. J., 2008, 95, 4193–4204.
75. J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kale and K. Schulten, J. Comput. Chem., 2005, 26, 1781–1802.
76. R. Selmar, R. Lieberei, B. Biehl and E. E. Conn, Physiol. Plant., 1989, 75, 97–101.
77. Y. Lin, Z. Cao and Y. Mo, J. Am. Chem. Soc., 2006, 128, 10876–10884.
78. H. Bühler, F. Effenberger, S. Förster, J. Roos and H. Wajant, ChemBioChem, 2003, 4, 211–216.
79. H. Bühler, B. Miehlich and F. Effenberger, ChemBioChem, 2005, 6, 711–717.

---

Footnote

† Electronic supplementary information (ESI) available: The crystal structures of HbHNLs, Mulliken charge populations of selected atoms that are relevant to the reaction, the binding free energy and the difference in the binding free energy of HbHNL and mutant systems, the structures of the reactant, the transition state, the intermediate for model B along the reaction path, wherein HCN formed by the deprotonation of His235, the different structures of the transition state during the “platform” for model B, the defined reaction coordinate for the release of acetone cyanohydrin and HCN, free energy profiles along the RC6 coordinate for the release of acetone cyanohydrin and RC5 for the release of acetone cyanohydrin in the Trp128Ala system, the key residues and the secondary structures from different windows along the RC7 coordinate for the release of HCN. See DOI: 10.1039/c4cp04032e

---