Theoretical study of the hydrolysis mechanism of dihydrocoumarin catalyzed by serum paraoxonase 1 (PON1): different roles of Glu53 and His115 for catalysis

Abstract

Serum paraoxonase 1 (PON1) is a calcium-dependent enzyme that can catalyze the hydrolysis of multiple substrates, including lactones, thiolactones, carbonates, esters and phosphotriesters, as well as the formation of a variety of lactones. To better understand the lactonase mechanism of PON1, the hydrolysis of dihydrocoumarin, which is considered as a native substrate of PON1, has been investigated by using a combined quantum mechanics and molecular mechanics (QM/MM) approach. Two possible reaction pathways with either Glu53 or His115 acts as the general base have been considered. On the basis of our calculations, these two pathways correspond to the overall energy barriers of 12.5 and 9.0 kcal mol⁻¹, respectively. During the catalytic reaction, if one of the two residues (Glu53 and His115) acts as the catalytic base, the other one forms strong hydrogen bonding interaction with the attacking hydroxide to facilitate the hydrolysis. However, mutation studies reveal that Glu53 is necessary for hydrolysis, whereas His115 is not essential but can promote the activity of PON1. Natural population analysis indicates that the catalytic Ca²⁺ does not act as Lewis acid but plays structural role in fixing the orientations of the substrate and related residues. In addition, Asp269 is found to coordinate with Ca²⁺ cation and facilitate the protonation of the alkoxide leaving group by forming hydrogen bond with lactone. These results can explain the fact that mutation of Glu53 results in the loss of activity of PON1, and the hydrolysis of dihydrocoumarin is unaffected by mutation of H115.

---

1. Introduction

Serum paraoxonase 1 (PON1) is a calcium-dependent enzyme that catalyzes the hydrolysis of multiple substrates, such as lactones, thiolactones, carbonates, esters and phosphotriesters including the nerve agents sarin and soman, and insecticides.^1,2 Nerve agents have been used as military weapons in historical warfare and insecticides have been widely used to protect agricultural crops.^3,4 However, they are harmful to both health and environment and thereby many efforts have been made to minimize this effect.⁵ Besides, PON1 also plays important roles in the alleviation of atherosclerosis and has been shown to be involved in drug metabolism.⁶ Recently, it has been shown that PON1 can catalyze the formation of lactones. Therefore, it is interesting and meaningful to explore the catalytic mechanism of PON1.

PON1 belongs to the mammalian PONs which include three subfamilies, while PON1 is the most investigated member owing to its paraoxonase activity.^7,8 Since multiple substrates can be catalyzed by PON1, it has been used to explore the catalytic promiscuity in both experimental and theoretical studies to demonstrate the mechanisms of protein functional evolution.^9–11 In addition, the detoxication of organophosphates by PON1,¹² the nerve agent VX hydrolysis by PON1 and the influence of temperature on the kinetics of PON1 have been studied in the past decades.^13,14 Besides, various mechanisms of PON1 for catalyzing the hydrolysis of esters have been proposed on the basis of kinetic studies and crystal structure.^10,14–16 For example, a proton shuttle mechanism for ester substrates has been suggested,^15,16 in which a histidine-dyad first deprotonates a hydrolytic water molecule to generate a hydroxide anion, which acts as a nucleophile to attack the ester carbonyl to form an oxyanionic tetrahedral intermediate; then the intermediate collapses to generate the acetate ion and phenol or alcohol. In 2014, Debord et al. proposed a similar mechanism for arylesterase and paraoxonase with the substrate phenyl acetate,¹⁴ and this proton shuttle mechanism was verified by Le et al. via QM/MM calculation.¹⁷ It should be pointed out that the His115, which was ever regarded as an important residue for the bound of the substrate, does not directly participate in catalysis.¹⁸ Further experiment by Khersonsky et al. reveals that the hydrolysis of phosphotriesters and dihydrocoumarin was unaffected by His115, but appears to participate in lactone and aryl ester hydrolysis.¹⁶ Besides, Ben-David et al. proposed a modified mechanism for the hydrolysis of lactones by rePON1. In this mechanism, the hydroxide nucleophile is generated via general base catalysis by both His115 and Glu53, in which the residue Glu53 was emphasized compared with the proton shuttle mechanism.¹⁰ Nevertheless, whether His115 or Glu53 is necessity for the hydrolysis reaction is still in dispute. In addition, the role of Asp269 has not been verified, which has been proposed by Ben-David.¹⁰ In this work, we adopt molecular dynamics (MD) simulations and combined quantum mechanical/molecular mechanical (QM/MM) method to explore these issues using dihydrocoumarin as the substrate. This methodology has been successful applied for exploring many enzymatic reaction mechanisms.^19–23

The three-dimensional structures of PON1 with and without substrate-analogue inhibitors have been obtained. The earliest crystal structure (PDB ID: 1VO4) was identified at pH 4.5, which contains a asymmetric unit, at 2.2 Å resolution.¹⁵ Later, the crystal structures of rePON1 (recombinant PON1 variant, PDB ID: 3SRE) and 2HQ/rePON1 (the complex with 2-hydroxyquinoline which is a well known competitive inhibitor of PON1, PDB ID: 3SRG) were determined at pH 6.5 and diffracted to 2.0 Å resolution.¹⁰ As native substrates, the lactones have relatively higher rates (k_cat/k_M ≈ 10⁶ M⁻¹ S⁻¹) than other substrates.²⁴ In particular, dihydrocoumarin has the closest structure compared to the lactam inhibitor 2HQ.¹⁰ Although dihydrocoumarin is not a typical lactone substrate of PON1,^25–27 its behaviors on hydrolysis do reveal a significant difference with lactone substrates but similar to paraoxon.^10,16 Thus, we chose the substrate dihydrocoumarin to study the hydrolysis mechanism, as shown in the Scheme 1.

Scheme 1 Catalytic reaction of dihydrocoumarin by PON1.

In this study, we attempted two possible reaction paths (path_a and path_b) for the hydrolysis of dihydrocoumarin catalyzed by PON1. In path_a, Glu53 acts as the catalytic base, whereas in path_b His115 is the catalytic base (Scheme 2). However, it is worth noting that both Glu53 and His115 form strong hydrogen bonds with the hydrolytic water molecule, and therefore both of the residues participate in the reaction in the two pathways. Besides, in order to explore the different catalytic roles of Glu53 and His115, we also performed two mutations in which Glu53 was mutated to Ala (E53A) or H115 to Gly (H115G). From these studies, we can identify the most likely reaction pathway and the functions of key residues.

Scheme 2 The proposed hydrolysis mechanism of dihydrocoumarin catalyzed by PON1.

2. Computational details

2.1 System pretreatment

The initial enzyme–substrate model was constructed on the basis of the X-ray crystal structure of mammalian serum PON1 (PDB ID: 3SRG).¹⁰ In this crystal structure, the competitive inhibitor 2-hydroxyquinoline (2HQ) is bound to the catalytic calcium ion, as shown in Fig. 1b. And it has been found that the docking mode of dihydrocoumarin with both open and closed structures is overall similar to that of 2HQ. Thus, to construct the enzyme–substrate model, the NH of 2HQ was firstly replaced by an oxygen atom, and the obtained dihydrocoumarin (signed as OCH) was preliminarily optimized through Gaussian 03 program package,²⁸ and then was docked into the active site. As a result, the catalytic calcium ion completed seven-coordination and this structure was used as the starting point for molecular dynamics simulations.

Fig. 1 (a) Crystal structure of PON1 (PDB ID: 3SRG); (b) active site derived from the crystal structure; (c) active site structure of the PON1-substrate complex derived from a series of minimizations and MD simulation, in which the inhibitor 2HQ has been replaced by the substrate OCH. Calcium coordination bonds were shown in dash lines. All the distances are shown in angstroms.

2.2 Molecular dynamics simulations

Firstly, the protonation states of all residues were determined at pH = 8.0 using the PDB2PQR suite of programs^29–31 in combination with actual reaction conditions and checked by the VMD program.³² According to experiment conditions and the pKa values, His115 and Glu53 were set to deprotonated and Asp269 to the protonated. It should be pointed that the aspartic acid residue could act as a proton donor in some calcium-dependent hydrolases (nucleoside hydrolase and α-mannosidase) although the carboxyl group of aspartic acid is usually deprotonated.^33,34 Thus, Asp269 has been suggested to be protonated by comparing with the structurally and functionally related enzyme Staphylococcus aureus.^10,35 The missing hydrogen atoms were added by the HBUILD program of the CHARMM package.³⁶ Subsequently, the computational model was solvated into a pre-equilibrated sphere which contains 4499 TIP3P waters, and 13 Na⁺ were added to neutralize the system (Fig. S1†). Afterwards, 20 ns molecular dynamics simulation was carried out with the CHARMM22/CMAP all-atom force field^36–38 for the protein after a series of minimizations. The system was basically equilibrated after 6 ns with the root-mean-squared deviation (RMSD) value of the backbone atoms of 2.0 Å, as shown in the ESI (Fig. S2†).

2.3 QM/MM calculations

The last snapshot derived from the 20 ns molecular dynamics simulation was taken as the reactant complex for the subsequent QM/MM calculations. The whole system was divided into two regions, the QM and MM regions. The QM region was treated by DFT with B3LYP functional^39–41 in Turbomole⁴² and the MM region was characterized by the CHARMM22 force field³⁸ in DL-POLY.⁴³ In our model, the QM region contains all the calcium-coordinated residues (Glu53, His115, Asp269, Asn168, Asn224, Asn270 and W1), W2 (the proposed hydrolytic water), the substrate OCH and the Ca²⁺ cation, as shown in Fig. 2a. Hydrogen link atoms were used to saturate the valency of atoms across the QM/MM boundary.⁴⁴ All protein residues and water molecules within 15 Å from the carbonyl carbon of substrate OCH were allowed to move, whereas the rest of the system was kept frozen. During the QM/MM calculation, geometry optimizations were implemented by the hybrid delocalized internal coordinates (HDLC) optimizer.⁴⁵ The quasi-Newton limited memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) method^46,47 was used for locating minima and the partitioned rational function optimization (P-RFO) algorithm⁴⁸ was used for transition state searches. To recognize the transition states, we firstly scanned the potential energy profile along the reaction coordinate (for example Fig. S3†), then the highest point on the energy curve was chosen to optimize the transition state. The last structure from the PES scan was optimized to find the intermediate, which was employed to start another scan to find the consecutive transition state and intermediate on the reaction path. The geometry optimizations was calculated at B3LYP/6-31G(d,p) level. Then the single-point energy calculations were further performed at B3LYP/6-311++G(2d,2p) level which includes diffuse functions and double polarization functions on each atom. Besides, the natural population atomic charges of QM region atoms were calculated at 6-311++G(2d,2p) level.^49–51 All the QM/MM calculations were performed with ChemShell program,^52,53 which combines Turbomole module⁴² and DL-POLY program⁴³ for QM and MM regions, respectively.

Fig. 2 (a) The selected QM region in our QM/MM calculations. Possible hydrogen bonds are shown in red dash lines and calcium coordination bonds in black dash lines; (b) QM/MM optimized active site structure (R_ab) of the enzyme–substrate complex. All the distances are shown in angstroms.

3. Results and discussion

3.1 Structure of enzyme–substrate complex

As shown in Fig. 1a, PON1 is a six-bladed β-propeller with a single active site with two Ca²⁺ ions in its central tunnel. According to previous studies,^10,15 the crystal structure of rePON1 (PDB ID: 3SRG) is ∼90% identical in sequence to both rabbit and human PON1s and shows the same enzymatic specificity. Therefore, the enzyme–substrate complex was obtained by changing the competitive inhibitor (2HQ) in the active-site of 2HQ/rePON1 to the substrate dihydrocoumarin (OCH), as shown in Fig. 1b and c. After MD and QM002FMM optimization, there are some slight deviations on the positions of the residues, which are supposed to be reasonable considering the structural difference between the inhibitor and the substrate. In Fig. 2b, the nucleophile W2 forms hydrogen bonds to Glu53 and His115 with lengths of 1.83 Å and 1.99 Å, respectively. Compared to the QM/MM calculation results on lactone hydrolysis by Le et al.,¹⁷ the active site structure of the enzyme–substrate complex (Fig. 2b) has some slight different from their structure: a hydrogen bond exists between Asp269 and OCH with length of 1.88 Å, indicating its possible role in facilitating the protonation of the alkoxide leaving group. Since path_a and path_b start from the same reactant, the finally optimized enzyme–substrate complex is named as R_ab (Fig. 2b).

3.2 Reaction pathways

Path_a. The optimized structures of transition states, intermediates and product are shown in Fig. 3, and the energy profile is displayed in Fig. 4. To clearly show the structural changes in Fig. 3, we use blue dashed lines to show coordination bonds and blue font to distinguish the three residues (Asn168, Asn224 and Asn269) which coordinate with the calcium ion but are not directly involved in the hydrolytic reaction. The whole catalytic reaction contains two elementary steps: the nucleophilic attack of hydrolytic water molecule to the carbonyl carbon and the cleavage of ester linkage. In the first step, the deprotonated Glu53 acts as the catalytic base to abstract a proton from water molecule, and simultaneously the generated hydroxyl group attacks on the carbonyl carbon atom of the substrate, generating the oxyanionic intermediate (IM_a). We also note that, during the reaction, His115 always forms strong hydrogen bond to the hydrolytic water, which facilitates the nucleophilic attack of hydrolytic water molecule to the carbonyl carbon of the substrate. The structure of TS1_a shows that the proton transfers from the hydrolytic water to Glu53 is much later than the nucleophilic attack. Along with the protonation of Glu53, the coordination between Ca²⁺ and Glu53 gradually weakens with the Ca²⁺–O distance increase from 2.42 Å in R_ab to 2.62 Å in IM_a. The first step is calculated to be quite easy, corresponding to an energy barrier of 5.6 kcal mol⁻¹. In the second step, the ester bond is cleaved, generating the final product (P_a). During this step, Asp269 donates a proton to the substrate, forming an alcoholic hydroxyl. This step corresponds to an energy barrier of 10.1 kcal mol⁻¹. However, on the basis of energy profile shown in Fig. 4, the overall energy barrier of path_a is suggested to be 12.5 kcal mol⁻¹. The free energy barrier estimated from the experimental kinetic constant (24) is 14.4 kcal mol⁻¹. Therefore, the calculated barrier is basically consistent with the experiments.

Fig. 3 Optimized geometries for various species in path_a. For clarity, the distances of calcium coordination bonds are shown in blue and the other key bond distances in black. All the distances are in angstroms.

Fig. 4 Energy profiles of dihydrocoumarin hydrolysis reaction catalyzed by PON1 at the B3LYP/6-31++G(2d,2p) level. All the energies of the four reactant models are set to zero.

Path_b. In path_b, His115 acts as the catalytic base which abstracts a proton from the hydrolytic water. During the catalytic reaction, similar cases as in path_a were found. For instance, Glu53 always forms a strong hydrogen bond to the water molecule, and Asp269 donates a proton to the substrate, etc., as shown in Fig. 5. However, some differences can also be found in these two pathways. One is energy barriers of the two elementary steps. They are 4.2 and 9.0 kcal mol⁻¹ in path_b for the first and second steps, respectively, indicating His115 to be the most likely catalytic base.

Fig. 5 Optimized geometries for various species in path_b. For clarity, the distances of Ca²⁺ coordination bonds are shown in blue and the other key distances in black. All the distances are in angstroms.

Over all, the calculation results on path_a and path_b show that either Glu53 or His115 can act as the catalytic base. During the catalytic reaction, if one of the two residues functions as the catalytic base, another residue facilitates the reaction by forming hydrogen bond with the hydrolytic water molecule. Furthermore, these two pathways may compete with each other.

It should be not that, similar to the proposed mechanisms in the references,^10,15, both in path_a and path_b, the protonation states of the two catalytic bases (H115 and E53) and Asp269 are not recovered after the hydrolysis of the substrate. Thus, to restart another catalytic cycle, a series of proton shuttle processes are needed, and the calculated energy barriers in path_a and path_b may do not present the real cases.

3.3 Mutations of His115 and Glu53

To explore the catalytic role of Glu53 and His115, we further constructed two computational models based on the crystal structure of mammalian serum PON1 (PDB ID: 3SRG).¹⁰ Firstly, the residue Glu53 or His115 was manually mutated to Ala (E53A) or Gly (H115G) by VMD program,³² then several nanoseconds of MD simulations were performed to equilibrate the systems. Finally, a series of QM/MM calculations were carried out on the two mutated enzymes. The corresponding reaction pathways are denoted as path_c and path_d, and the optimized structures involved in these two pathways are shown in Fig. 6 and 7, respectively. After mutation of His115 by Gly (path_c), Glu53 still acts as the catalytic base, but the calculated energy barrier for the first step increases to 12.2 kcal mol⁻¹ from 5.6 kcal mol⁻¹ in path_a, which suggests the importance of His115 for the first step. However, the energy barrier of the second step of path_c decreases to 4.2 kcal mol⁻¹ from 10.1 kcal mol⁻¹ in path_a, which means His115 is not conducive to the cleavage of the ester bond. By checking the structures in Fig. 6, one can see that the mutation of His115 by Gly decreases the steric exclusion in that region, thereby facilitating the cleavage of ester bond. In general, after mutation of His115 to Gly, the first step is calculated to be rate limiting, but the overall energy barrier is still comparable to the experimentally free energy barrier.

Fig. 6 The selected QM region and optimized geometries for various species in path_c (His115G). For clarity, the distances of calcium coordination bonds are shown in blue and the key bond distances in black. All the distances are in angstroms.

Fig. 7 The selected QM region and optimized geometries for various species in path_d (Glu53A). For clarity, the distances of Ca²⁺ coordination bond are shown in blue and the key bond distances in black. All the distances are in angstroms.

However, the mutation of Glu53 to Ala (path_d) increases the energy barrier of the first step, as shown in Fig. 4. After this mutation, in the first step, the catalytic base is still His115, which is similar to that in path_b, but the energy barrier increases to 16.7 kcal mol⁻¹, suggesting that the activation of hydrolytic water becomes more difficult without the assistance of Glu53. For the second step, the energy barrier (9.9 kcal mol⁻¹) is almost the same as in path_b (9.0 kcal mol⁻¹). It is reasonable, because Ala53 is far from the ester bond of the substrate and should have no obvious influence on its cleavage. Since the overall energy barrier of path_d is 22.8 kcal mol⁻¹, the mutation of Glu53 to Ala greatly decreases the activity of PON1, which is in consistent with the observation that mutation of Glu53 by all other amino acids results in the loss of activity.¹⁰

During the hydrolysis of dihydrocoumarin, the hydrolytic water forms strong hydrogen bonds with both His115 and Glu53. Our calculations reveal that both of the residues can function as the catalytic base. But owing to the different positions in the active site, and Glu53 also coordinates with the Ca²⁺, His115 and Glu53 have different influence on the catalytic efficiency.

The energy barrier in path_c shows that the hydrolysis can still occur without the assistance of His115, indicating His115 is not essential for catalysis, but it is important for promoting the activity of PON1 by forming hydrogen bond with hydrolytic water molecule.

As for Glu53, the calculation results of path_d reveal that, without the assistance of Glu53, the activation of hydrolytic water by His115 alone becomes more difficult. Besides, we also found that, when Glu53 acts as the catalytic base, the electrostatic interaction between the Ca²⁺ cation and the carboxyl oxygen of Glu53 gradually weakens along with the catalytic reaction. It is understandable, because a proton has been transferred from the hydrolytic water to the carboxyl group of Glu53. The distances between the Ca²⁺ cation and the carboxyl oxygen of Glu53 are displayed in Fig. S4,† which show a gradual increase of this distance in path_a and path_c. Based on the above results, we come to a conclusion that Glu53 not only interacts with the attacking water molecule but also coordinates with Ca²⁺ during the whole reaction, which agrees with the result obtained by Le et al.¹⁷ In addition, Asp269 is found to coordinate with Ca²⁺ cation and forms hydrogen bonding interaction with the lactone to facilitate the protonation of the alkoxide leaving group.

To exam the role of Ca²⁺ in the catalytic reaction, we further calculated the Natural Population Atomic (NPA) charges of Ca²⁺ and some key atoms, which are listed in the ESI (Table S1 to S4†). As shown in these tables, the charge of Ca²⁺ almost remains unchanged during the catalytic reaction. Thus, the Ca²⁺ cation does not act as Lewis acid but coordinates with the key residues and substrate, which mainly plays role in fixing the positions of the substrate and catalytic base (Glu53).

4. Conclusion

In this paper, the detailed mechanism for the hydrolysis of dihydrocoumarin catalyzed by PON1 has been studied by using QM/MM approach. Two possible reaction pathways with either Glu53 or His115 acts as the general base have been considered. Our calculation results reveal that these two pathways are coextensive, corresponding to overall energy barriers of 12.5 and 9.0 kcal mol⁻¹, respectively. During the catalytic reaction, if one of the two residues acts as the catalytic base, the other one forms strong hydrogen bonding interaction with the attacking hydroxide to facilitate the reaction. Calculations on mutant of H115G reveal that, without the assistance of His115, Glu53 alone can still act as the general base to activate the hydrolytic water molecule, and hydrolysis of dihydrocoumarin corresponds to an overall energy barrier of 12.2 kcal mol⁻¹. However, calculations on mutant of E53A indicate that, without the assistance of Glu53, the activation of hydrolytic water by His115 alone becomes more difficult, and the overall energy barrier increases to 22.8 kcal mol⁻¹. Thus, we come to the conclusion that Glu53 is necessary for hydrolysis, whereas His115 is not essential but can promote the activity of PON1. These results can well explain the experimental results that mutation of Glu53 results in the loss of activity of PON1, and the hydrolysis of dihydrocoumarin is nearly unaffected by mutation of H115. Natural population analysis indicates that the catalytic Ca²⁺ does not act as Lewis acid but mainly plays role in fixing the positions of the substrate and catalytic base (Glu53). Our results may provide useful information for further understanding the catalytic mechanism of PON1 and the roles of pocket residues. It also provides new perspective for the regulation of enzyme activity.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (21373125, 21573127, 21403201).

References

1. D. I. Draganov and B. N. La Du, Naunyn-Schmiedeberg's Arch. Pharmacol., 2004, 369, 78–88.
2. B. N. La Du, M. Aviram, S. Billecke, M. Navab, S. Primo-Parmo, R. C. Sorenson and T. J. Standiford, Chem.-Biol. Interact., 1999, 119, 379–388.
3. N. B. Munro, K. R. Ambrose and A. P. Watson, Environ. Health Perspect., 1994, 102, 18–38.
4. F. M. Raushel, Curr. Opin. Microbiol., 2002, 5, 288–295.
5. A. N. Bigley and F. M. Raushel, Biochim. Biophys. Acta, 2013, 1834, 443–453.
6. A. J. Lusis, Nature, 2000, 407, 233–241.
7. G. S. Getz and C. A. Reardon, Curr. Opin. Lipidol., 2004, 15, 261–267.
8. C. J. Ng, D. M. Shih, S. Y. Hama, N. Villa, M. Navab and S. T. Reddy, Free Radical Biol. Med., 2005, 38, 153–163.
9. M. Ben-David, G. Wieczorek, M. Elias, I. Silman, J. L. Sussman and D. S. Tawfik, J. Mol. Biol., 2013, 425, 1028–1038.
10. M. Ben-David, M. Elias, J.-J. Filippi, E. Duñach, I. Silman, J. L. Sussman and D. S. Tawfik, J. Mol. Biol., 2012, 418, 181–196.
11. F. Duarte, B. A. Amrein and S. C. L. Kamerlin, Phys. Chem. Chem. Phys., 2013, 15, 11160–11177.
12. J. E. Chambers, H. W. Chambers, E. C. Meek, K. E. Funck, M. H. Bhavaraju, S. R. Gwaltney and R. B. Pringle, Toxicol. Sci., 2015, 143, 46–53.
13. M. W. Peterson, S. Z. Fairchild, T. C. Otto, M. Mohtashemi, D. M. Cerasoli and W. E. Chang, PLoS One, 2011, 6, e20335.
14. J. Debord, J. C. Bollinger, M. Harel and T. Dantoine, Biochimie, 2014, 97, 72–77.
15. M. Harel, B. Brumshtein, R. Meged, H. Dvir, R. Ravelli, A. McCarthy, L. Toker, I. Silman and J. Sussman, Arch. Ind. Hyg. Toxicol., 2007, 58, 347–353.
16. O. Khersonsky and D. S. Tawfik, J. Biol. Chem., 2006, 281, 7649–7656.
17. Q. A. T. Le, S. Kim, R. Chang and Y. H. Kim, J. Phys. Chem. B, 2015, 119, 9571–9585.
18. D. T. Yeung, D. E. Lenz and D. M. Cerasoli, FEBS J., 2005, 272, 2225–2230.
19. T. Dresselhaus, N. D. Weikart, H. D. Mootz and M. P. Waller, RSC Adv., 2013, 3, 16122–16129.
20. J. Šebera, L. Trantírek, Y. Tanaka, R. Nencka, J. Fukal and V. Sychrovsky, RSC Adv., 2014, 4, 44043–44051.
21. R. Wu, W. Gong, T. Liu, Y. Zhang and Z. Cao, J. Phys. Chem. B, 2012, 116, 1984–1991.
22. H. Hu and W. T. Yang, Annu. Rev. Phys. Chem., 2008, 59, 573.
23. R. A. Friesner and V. Guallar, Annu. Rev. Phys. Chem., 2005, 56, 389–427.
24. O. Khersonsky and D. S. Tawfik, Biochemistry, 2005, 44, 6371–6382.
25. S. Billecke, D. Draganov, R. Counsell, P. Stetson, C. Watson, C. Hsu and B. N. La Du, Drug Metab. Dispos., 2000, 28, 1335–1342.
26. M. Harel, A. Aharoni, L. Gaidukov, B. Brumshtein, O. Khersonsky, R. Meged, H. Dvir, R. B. Ravelli, A. McCarthy, L. Toker, I. Silman, J. L. Sussman and D. S. Tawfik, Nat. Struct. Mol. Biol., 2004, 11, 412–419.
27. A. Aharoni, L. Gaidukov, O. Khersonsky, S. M. Gould, C. Roodveldt and D. S. Tawfik, Nat. Genet., 2005, 37, 73–76.
28. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. A. Montgomery Jr, T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, 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, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzalez and J. A. Pople, GAUSSIAN 03 (Revision D.01), Gaussian, Inc., Wallingford, CT, 2004.
29. H. Li, A. D. Robertson and J. H. Jensen, Proteins, 2005, 61, 704–721.
30. C. R. Sondergaard, M. H. M. Olsson, M. Rostkowski and J. H. Jensen, J. Chem. Theory Comput., 2011, 7, 2284–2295.
31. T. J. Dolinsky, P. Czodrowski, H. Li, J. E. Nielsen, J. H. Jensen, G. Klebe and N. A. Baker, Nucleic Acids Res., 2007, 35, W522–W525.
32. W. Humphrey, A. Dalke and K. J. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
33. W. Versees, K. Decanniere, R. Pelle, J. Depoorter, E. Brosens, D. W. Parkin and J. Steyaert, J. Mol. Biol., 2001, 307, 1363–1379.
34. K. Karaveg, A. Siriwardena, W. Tempel, Z. J. Liu, J. Glushka, B. C. Wang and K. W. Moremen, J. Biol. Chem., 2005, 280, 16197–16207.
35. Y. Tanaka, K. Morikawa, Y. Ohki, M. Yao, K. Tsumoto, N. Watanabe, T. Ohta and I. Tanaka, J. Biol. Chem., 2007, 282, 5770–5780.
36. B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan and M. Karplus, J. Comput. Chem., 1983, 4, 187–217.
37. J. H. Chen, W. P. Im and C. L. Brooks, J. Am. Chem. Soc., 2006, 128, 372–3736.
38. A. D. MacKerell Jr, D. Bashford, M. Bellott, R. L. Dunbrack Jr, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiorkiewicz-Kuczera, D. Yin and M. Karplus, J. Phys. Chem. B, 1998, 102, 3586–3616.
39. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
40. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
41. P. J. Stephens, F. J. Devlin, C. F. Chablowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627.
42. R. Ahlrichs, M. Bar, M. Haser, H. Horn and C. Kolmel, Chem. Phys. Lett., 1989, 162, 165–169.
43. W. Smith and T. R. Forester, J. Mol. Graphics, 1996, 14, 136–141.
44. A. H. de Vries, P. Sherwood, S. J. Collins, A. M. Rigby, M. Rigutto and G. J. Kramer, J. Phys. Chem. B, 1999, 103, 6133–6141.
45. S. R. Billeter, A. J. Turner and W. Thiel, Phys. Chem. Chem. Phys., 2000, 2, 2177–2186.
46. J. Nocedal, Math. Comput., 1980, 35, 773–782.
47. D. C. Liu and J. Nocedal, Math. Program., 1989, 45, 503–528.
48. A. Banerjee, N. Adams, J. Simons and R. Shepard, J. Phys. Chem., 1985, 89, 52–57.
49. A. E. Reed and F. Weinhold, J. Chem. Phys., 1985, 83, 1736–1740.
50. A. E. Reed, R. B. Weinstock and F. Weinhold, J. Chem. Phys., 1985, 83, 735–746.
51. A. E. Reed, L. A. Curtiss and F. Weinhold, Chem. Rev., 1988, 88, 899–926.
52. S. Metz, J. Kastner, A. A. Sokol, T. W. Keal and P. Sherwood, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 101–110.
53. P. Sherwood, A. H. de Vries, M. F. Guest, G. Schreckenbach, C. R. A. Catlow, S. A. French, A. A. Sokol, S. T. Bromley, W. Thiel, A. J. Turner, S. Billeter, F. Terstegen, S. Thiel, J. Kendrick, S. C. Rogers, J. Casci, M. Watson, F. King, E. Karlsen, M. Sjøvoll, A. Fahmi, A. Schäfer and C. Lennartz, J. Mol. Struct.: THEOCHEM, 2003, 632, 1–28.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra09735a

---