From a generated PES, one can determine the relative energies of species involved, the sequence in which they occur, and the activation barrier(s) associated with individual steps or the overall mechanism. Furthermore, they can provide more insights than a simple indication of a path of sequential mechanistic structures and their energetic relationships. The investigation into the activation of O2 by alpha-ketoglutarate-dependent dioxygenase (AlkB) clearly shows the opportunity for spin inversion, where one can see that the lowest energy product may be formed via several possible routes. In the investigation of uroporphyrinogen decarboxylase III (UROD), the use of QM/MM methods allowed for the inclusion of the anisotropic protein environment providing greater insight into the rate-limiting barrier. Lastly, the mechanism of 6-phospho-α-glucosidase (GlvA) was discussed using different active site models. In particular, a continuum model PES was compared to the gas-phase PES.

1. Introduction

For a chemical reaction, enzymatic or nonenzymatic, the reactants, intermediates, and products all exist on a multidimensional surface. With this surface, a reaction is perfectly described by the statistical average of all possible paths from reactants to products via all possible intermediates [1]. However, a system with N atoms would require computing a (3N-6)-dimensional surface. Hence, if X number of points are to be computed for each of the (3N-6) dimensions, then calculations must be done [1]. Thus, such an undertaking is usually only computationally feasible for chemical models that consist of a few atoms. In contrast, studies on enzymatic mechanisms often necessarily require large chemical models consisting of important active site functional groups and cofactors and the substrate. Indeed, the cluster-based approach for investigating biochemical reactions typically use chemical models containing 200 atoms or more [2]. Thus, for such systems, it is impossible to determine the complete PES. Instead, a “slice” of the surface is typically constructed that involves only two coordinates, energy and reaction coordinate, and is commonly referred to as the PES [3].

The use and applicability of such “reduced-dimensionality” PESs reflects the fact that chemists and biochemists are usually only interested in key, mechanistically relevant structures including, for example, the (i) isolated reactants, (ii) reactive complex, (iii) transition structures (TS), (iv) intermediate(s), (v) product complex, and (vi) the separated products (Figure 1) [1]. It is noted that in general, computational studies on enzymatic reactions are investigated with the enzyme-substrate complex already formed. This is due in part to the inherent difficulties associated with modeling substrate binding and product release such as, for instance, consideration of clathrate waters. Furthermore, as previously stated [4], in some cases substrate binding is “seldom the most interesting part of the potential surface”.

A “reduced dimensionality” PES can provide a considerable amount of information and insight into a chemical system. For example, a common goal of computational chemistry is the accurate determination of the thermochemistry of a particular system. From a PES, one can determine the relative energies of species involved, the sequence in which they occur, and the activation barrier(s) associated with individual steps or the overall mechanism. Consequently, one can determine if a reaction step or pathway is feasible under the given conditions. With the appropriate corrections, these potential energies can be converted to enthalpies or free energies as has been previously discussed in detail [3]. For enzyme-catalysed reactions, free energy barriers between 55–75 kJ mol−1 are typical with the upper thermodynamic limit generally held to be 84–105 kJ mol−1 [4]. For barriers greater than 40 kJ mol−1, transition state theory is sufficiently accurate to enable determination of the associated rate constant, albeit perhaps to within several orders of magnitude (1) [4]

The reliability and accuracy of any PES is dependent on the computational model used; that is, the choice of computational method plus chemical model. A common approach, in particular for investigations of enzymatic processes, is the use of “cluster chemical models” in combination with a reliable quantum mechanical (QM) or density functional theory (DFT) methods. Typically, in such cases, only those regions of the active site-substrate complex involved in bond making and breaking processes and active site residues or functional groups that play a direct role are included [5]. For the study of biocatalytic processes, DFT methods, in particular, B3LYP, are presently the most widely used [68]. This is due to the fact that they include electron correlation effects yet are computationally less expensive than more conventional wave function-based electron correlation methods. Hence, they can be applied to comparatively large chemical models [4, 9]. Furthermore, they have been shown to often be able to provide accurate and reliable structures and thermochemistry. Indeed, B3LYP has been shown to often be highly accurate [5] with relative energy errors for 1st and 2nd row atom and transition metal-containing systems of just ~13.0 and ~21.0 kJ mol−1, respectively [4]. The surrounding protein environment, more specifically its general polarity, is then included via use of a polarizable continuum model (e.g., IEF-PCM) [10, 11].

More recently, QM/MM methods have been increasingly applied to the elucidation of enzymatic reaction PESs. This is due in part to the fact that (i) they are able to model a greater portion of the enzyme through their combined use of QM and molecular mechanics (MM) methods, and thus (ii) are able to model the polar and steric nonhomogeneity of the protein environment surrounding the active site. Furthermore, the computed energetics of the system have been shown to often converge faster with increasing chemical model size than for the above cluster/QM-based approach [12]. In addition, the role of the residues and mechanism is also less sensitive to increasing or varying model size [13]. However, QM/MM-based approaches suffer some of the same limitations as the alternative cluster/DFT approach such as the fact that the dynamic behaviour of the enzyme is not fully taken into account [13, 14]. The application of these methods to the elucidation of enzymatic pathways has been previously reviewed; see, for example, Senn and Thiel [15] and Llano and Gauld [16].

There are computational techniques available to determine a statistical average of possible alternative pathways of an enzymatic reaction [1, 3, 12, 17]. One common approach is to use a semiempirical (SE) method to describe the QM layer, while an MM method is used to describe the outer layer. The use of SE methods is due to the fact that extensive sampling of the protein, substrate, and solvent must be performed in order to calculate the free energies. This sampling requires large numbers of calculations that also include the inner layer, hence the need for highly efficient methods that describe this region. Alternatively, the CPMD method typically describes the inner layer using DFT methods. However, due to the added computational costs, simulation timescales are limited to 0.1 ns, significantly shorter than those noted above using SE methods [18]. In addition, QM/MM MD can also be used to determine free energies. Similar to CPMD, QM/MM MD uses DFT methods to describe the inner layer, where the complete QM/MM model is simulated for a period of time [19]. For electron transfer (ET) processes, a combined QM/MM and MD approach can be used to obtain free energies [20]. More specifically, structures of the ET reactant and product are obtained using QM/MM. An MD simulation is then run on these structures, where the equilibrium averages of the MD trajectories are used to calculate the free energies. It is noted that Zhang et al. [13] investigated the effects of structural fluctuations on the reaction energy barrier in the catalytic mechanism of acetylcholinesterase. Using eight different starting conformations obtained from a 1 ns MD simulation, they found that while the enzyme-substrate structural fluctuations led to differences in the calculated barrier of approximately ±8 kJ mol−1, the mechanistic details remained very consistent.

Potential energy surfaces, however, can provide more insights than a simple indication of a path of sequential mechanistic structures and their energetic relationships. For example, species found along an enzymatic mechanism may have more than one electronic state that lie close in energy and within the boundaries imposed by the “enzyme thermodynamic limit” (see above), particularly in the case of metalloenzymes. Often, such species contain one or more unpaired electrons. Importantly, however, it is possible that as a reaction progresses, the energy differences between states may vary, and in fact, may even change in their energetic ordering. Thus, a mechanism may proceed via “switching” between states, that is, undergoing a spin inversion (SI). Such reactions are termed multistate reactions (MSR). Hence, in the case of MSR pathways, a fundamental property is that the PES for each state may cross or may remain in close proximity due to the near degeneracy of the d-orbitals of metals within the complex [22]. The ability in changing of states may have a significant impact on the reaction. For example, for a given transition metal complex, a particular reaction may have a significant barrier, and thus may not feasibly occur. However, for another possible state of the complex, this reaction may have a markedly lower barrier such that it may now be feasible. While SI is typically forbidden, it has been suggested that it may be common in transition metal chemistry [23] and by extension, metalloenzymes. The points at which SI may occur on a full-dimensionality PES are challenging to describe accurately due to the possibility of strong spin-orbit coupling between the two states. However, these crossing points are not generally at stationary points [23]. Fortunately, modern computational tools are good at locating stationary points [23]. Thus, one can map a “reduced-dimensionality” PES for each state, and where they cross indicates a possible region in which SI may occur [23].

In this present paper, some of the ways in which PESs may be exploited to provide insights into enzymatic mechanisms and some factors that may influence their accuracy and reliability are discussed. These are illustrated through a review of several enzymatic systems that we have recently examined. More specifically, this article highlights three major topics in the application of PES in the elucidation of enzymatic mechanisms through examples: (i) the use of multistate reactivity PESs in elucidating possible pathways of O2 activation in the nonheme iron DNA repair metalloenzyme alpha-ketoglutarate-dependent dioxygenase (AlkB) and the movement of electrons within the active site during the reactions progression [21], (ii) how differences in model and method (QM versus QM/MM) choice influence the mechanistic PES of Uroporphyrinogen Decarboxylase III (UROD), a key “protein-only” enzyme in porphyrin biosynthesis, [24, 25], and (iii) the influence of solvation and the use of “broken” PESs in investigations on the mechanism of the NAD+ and Mn2+ dependent metalloenzyme 6-phospho-α-glucosidase (GlvA) [26].

2. Multistate Reactivity Potential Energy Surfaces: Activation of O in AlkB

An organisms “blueprint” is encoded in its DNA, in particular the sequence of its nucleobases. Thus, the fidelity and integrity of these DNA components is essential to its proper functioning and genetic transmission. However, both internal (e.g., metabolic byproducts) and external (e.g., radiation) factors can damage nucleobases via depurination, oxidation or deamination processes [2734]. Alternatively, they may be alkylated at their oxygen or nitrogen centres [35]. Consequently, cells have developed several approaches to mediate or repair such alkylation damage [3640]. In particular, the AlkB family of enzymes have evolved a novel mechanism by which they oxidatively dealkylate several such damaged nucleobases, specifically 1-meA, 3-meC, 3-meT, and 1-meG (Scheme 1).


The AlkB proteins belong to the α-ketoglutarate-Fe(II)-dependent dioxygenase superfamily but are the only members to catalyze oxidative dealkylation [41, 42]. More specifically, they are nonheme iron containing proteins in which the Fe(II) ion is ligated by a facial triad of residues consisting of two histidines (His131, and His187) and an aspartate (Asp133) [43, 44]. In addition, they also require α-ketoglutarate (α-KG) as a cosubstrate [41, 42] that is converted to succinate during the mechanistic activation of O2.

The overall mechanism of AlkB is thought to occur in three stages: [21] (i) O2 activation to give an ferryl-oxo (FeIV= O) moiety, (ii) reorientation of the FeIV = O oxygen from the axial to an equatorial position, and (iii) oxidative dealkylation of the alkylated nucleobase [45]. It is noted that a reorientation of the FeIV = O has been suggested to occur in other enzymes of this superfamily such as the clavaminate synthase [46], while the dealkylation step has been previously computationally studied [21].

The overall proposed mechanism by which the O2 moiety is activated with concomitant oxidation of α-KG to succinate and CO2 is shown in Scheme 2. In particular, following binding of α-KG and O2, a self-redox process results in formation of a ferric-superoxide FeIII- species [47, 48]. This then nucleophilically attacks the α-keto group of α-KG, resulting in its decarboxylation (i.e., loss of CO2) and the formation of succinate with concomitant generation of an “axial-positioned” FeIV = O intermediate. However, the exact details by which the electrons released from the α-KG cosubstrate are used in the O2 activation process, as well as the influence of the systems electronic state upon the choice of preferred pathway, remain unclear.


We examined [21] this key stage using a cluster/DFT-based computational model. In particular, optimized geometries were obtained at the PB-SCRF-B3LYP/LACVP(d) level of theory with a dielectric constant of 4.0. Frequency calculations were performed on these structures in order to characterize them as minima or transition states. Relative total energies were calculated via single-point energy calculations at the PB-SCRF/B3LYP/LACV3P+(d,p) level of theory based on the above geometries. All calculations were performed using the Gaussian 03 [49] and Jaguar 5.5 [50] software packages. For the chemical model, the coordinates of the key active site residues and groups were obtained from a crystal structure of AlkB (PDB: 2FD8) and is shown in Scheme 3 [45]. In order to maintain integrity of the model, a minimum number of atoms remote from the reactive regions were held fixed at their crystal structure positions. Complete computational details are provided in the article by Liu et al. [21].


The potential energy surfaces obtained for activation of O2 in AlkB to give the highly reactive oxo- and oxyl-type intermediates (IC4) are shown in Figure 2. It should be noted that the oxidation state of the Fe centre in each complex was in part determined by calculating its spin density [21]. In addition, in order to provide further insights into electron transfer processes, a stylized electronic representation of the bond making/breaking region is also shown for the minimum energy path. Assuming the commonality of spin changes in transition metal chemistry, this minimum energy path involves spin inversion from the septet to quintet surface [23].

The preferred mode for O2 binding to the Fe(II) centre is end-on (). Four overall spin combinations are possible for the resulting dioxygen-bound α-KG−Fe(II) active-site complex (RC): singlet, triplet, quintet, and septet. The overall spin-triplet complex, , was chosen to be the reference energy level throughout this work.

In both the singlet and triplet states, O2 binding has occurred without electron transfer from the Fe(II) to O2 moiety. The spin-singlet state, , lies significantly higher in energy than by 57.1 kJ mol−1 (not shown). This is likely due to less favourable antiferromagnetic alignment of the unpaired electrons of 1O2 with the low-spin Fe(II) centre in , compared to 3O2 with the high-spin Fe(II) centre in [21]. In contrast, in the quintet and septet spin states, O2 binding occurs with charge transfer. That is, a high-spin sextet Fe(III) ion (S = 5/2) is formed with concomitant reduction of 3O2 (valence electron configuration in a fragment orbital picture) to . Importantly, the spin-septet 7[FeIII-] complex lies 12.5 kJ mol−1  lower in energy than , while the spin-quintet complex is just 5.7 kJ mol−1 higher in energy than .

While was chosen as the “reference relative energy” for all complexes, the mechanistic details of the spin-triplet PES are not discussed herein due to its considerably higher activation barriers compared to the spin-quintet and spin-septet surfaces (Figure 2).

For , the lowest energy reactant complex, the first step is nucleophilic attack of the ferric-superoxide (7[FeIII-]) group at the pyruvate’s C2-carbon. It is noted that the ability of FeIII-  to act as a strong nucleophile has been previously investigated by both experimental [5154] and theoretical [5558] studies. This step proceeds via with a barrier of only 19.1 kJ mol−1 to give the intermediate complex lying just 1.8 kJ mol−1 lower in energy than (Figure 2). In , the O2 moiety forms a peroxy bridge between Fe(III) and the pyruvate cosubstrates C2-centre [21]. In the subsequent step, the cosubstrate’s C2–COO  bond is cleaved with the electrons formally moving into the Fe-peroxo moiety. As a result, the iron is reduced to Fe(II), while the peroxy moiety now has a charge of −2, and consequently, its O–O bond lengthens considerably. This step proceeds via with a barrier of 34.4 kJ mol−1 to give the ferric-peroxide intermediate lying 15.9 kJ mol−1 higher than . The peroxide O–O bond is then cleaved with concomitant two-electron oxidation of the Fe(II) centre to Fe(IV) via at a cost of just 5.7 kJ mol−1 with respect to . The resulting Fe(IV)-oxo-type product complex is considerably lower in energy than by 168.6 kJ mol−1. Thus, overall, this pathway is enzymatically feasible.

On the quintet-state PES, the first step is also found to be nucleophilic attack of the 5[FeIII-] group at the cosubstrate’s C2 centre. However, this now occurs with decarboxylation of the pyruvate and concomitant reduction of both Fe(III) and to Fe(II) and , respectively. Importantly, this reaction proceeds via with a barrier of only 32.4 kJ mol−1 to give the energetically very low lying intermediate (−160.8 kJ mol−1). Based on their structural similarity, it is suggested that could be regarded as the reaction-coordinate equivalent of [21]. Subsequently, can undergo a stepwise two-electron transfer from Fe(II) to the O–O unit [21]. The first electron transfer occurs via without a barrier. It is noted that optimized structures are obtained by minimization of the electronic energy at 0 K. As a result, on flat PESs the inclusion of energy corrections via, for example, single-point calculations, enthalpy, or Gibb’s free energy corrections can result in TSs having lower relative energies than the reactant or product which it interconnects. This is typically taken to indicate that the particular reaction essentially occurs without a barrier. In the resulting intermediate , the electron has transferred from the Fe(II) into the orbital. As a result, the peroxo moiety, in a fragment MO picture, has the electronic configuration . The second electron transfer proceeds via at a cost of 7.1 kJ mol−1 to give the ferryl-oxo product complex. is the ferryl-oxo compound and is found to lie 25.3 kJ mol−1 lower in energy than . This one-electron oxidation of Fe(III) again promotes an electron into to the orbital, thus giving the peroxo moiety the highly unstable electron configuration and resulting in scission of the O–O linkage. The complex lies significantly lower in energy than by 220.1 kJ mol−1. Importantly, however, it is also markedly lower in energy than both and by at least 40.8 kJ mol−1 (Figure 2).

Thus, by elucidating the pathways for the various possible spin states one can see that the lowest energy product may in fact be formed via several possible low energy routes. In particular, it may occur either at the beginning of the mechanism via an initial spin inversion from to , and thus providing a pathway with an overall barrier of 50.6 kJ mol−1. Alternatively, spin inversion may occur along the pathway, for example, as reacts and proceeds towards and thus to directly give , allowing for cleavage of the peroxo O–O bond with a possible pathway barrier of less than 50.6 kJ mol−1. Importantly, spin inversion would enable the mechanism to access and proceed via considerably more exoergic intermediates, for example, lies 181.3 kJ mol−1 lower in energy than [21]. The probability of SI occurring at a crossing-point is determined by the strength of spin-orbit coupling between the two states. Fortunately, as noted in the introduction, such crossing-points are not generally stationary points [23]. In addition, it is rarely found that the spin-orbit coupling is so strong that it affects the adiabatic energies of the PESs of the two states [23].

3. Effects of Explicitly Modeling the Protein Environment on the Catalytic Mechanism of Uroporphyrinogen Decarboxylase III

Porphyrin is an important biomolecule for all organisms. It plays a range of diverse key roles in, for example, proteins and enzymes involved in ligand transport, electron transfer, light harvesting, and redox mechanisms [6265]. Within cells, its biosynthesis occurs via a multistage multienzymatic process in which Uroporphyrinogen Decarboxylase III (UROD) catalyses the first branching point, the fifth step. Specifically, it catalyses the sequential nonsymmetric decarboxylation of the four acetates of uroporphyrinogen III (URO-III) to give coproporphyrinogen III (CP-III) [6264, 66].

UROD exists as a homodimer, and it has been shown experimentally that each active site is independent from the other [24]. Three electrostatic regions have been identified within each active site: a negative, a polar-positive, and a hydrophobic region [62]. The negatively charged region contains an invariant aspartyl (hUROD: Asp86) that is thought to help orientate the substrate for catalysis as well as stabilize various mechanistic intermediates [62, 67]. Helping in substrate binding and recognition, the polar-positive region contains several residues that interact with the carboxylates of URO-III [67, 68]. Importantly, of these polar-positive residues, it has been found that one or more active-site arginyl residues (hUROD: Arg37, Arg41 and Arg50) are catalytically essential [62, 6871]. While all four acetate decarboxylations are believed to involve the same catalytic residues, unfortunately, the exact mechanistic role of various active site residues remains unclear.

A general mechanism has been proposed for decarboxylation of the pyrrole-acetate and is shown in Scheme 4 [5961]. Specifically, two acids are thought to be involved. In the first step, an acid (HA) protonates the C2 carbon of the pyrrole ring. This destabilizes the carboxylate group, resulting in cleavage of the C3′–C3′′OO bond, that is, decarboxylation. The second acid (HB) then protonates the newly formed methylene C3′ carbon, while the first acid, now in its conjugate base form (A), abstracts a proton from the––group, thus regenerating HA.


Several experimental X-ray crystal structure and kinetic studies, based on this “blueprint”, have proposed possible identities for the two mechanistic acids. Specifically, from enzyme X-ray structures they obtained coupled with manual docking of the substrate, Martins et al. [64] suggested that HA is Asp86, while HB is Tyr164. In contrast, based on a X-ray crystal structure of a URODproduct complex, Phillips et al. [67] suggested that HA is an H2O but were unable to conclusively determine HB. More recently, Lewis and Wolfenden, [69] using kinetics and p measurements, suggested that HA and HB are Asp86 and Arg37, respectively.

Silva and Ramos [25] performed a computational investigation on the first acetate decarboxylation of URO-III as catalysed by UROD. More specifically, they used a cluster/DFT-type approach with a chemical model consisting of the aspartyl (Asp86) and arginyl (Arg37) R-groups and 1,3,4-methyl-2-acetyl pyrrole for the URO-II substrate (Figure 3). This model was derived from a crystal structure of a URODproduct complex [67]. A PCM-solvation approach was used to model the general affects of the surrounding protein environment. Using this computational model, they concluded that decarboxylation could proceed in accordance with the “blueprint” mechanism and that it was thermodynamically and enzymatically feasible with an overall barrier of 89.5 kJ mol−1 [25]. Notably, this energy corresponded to the initial step, that is, proton transfer from Arg37 onto C2 of the substrate pyrrole ring. Thus, the R-group guanidinium of the active-site residue Arg37 was concluded to be both suitably positioned and capable of acting as the general acid HA that protonates the substrate’s C2 centre. They further concluded that the second proton, from HB, is donated by the solvent. However, based on the experimental observation that stereochemistry is retained at the C3′ centre from which the CO2 is lost, it has been concluded that HB must be an active site residue [59, 60, 72].

We reexamined the catalytic mechanism of UROD using an alternate computational approach. In particular, it was investigated [24] through combined QM and MM methods in the ONIOM formalism with mechanical embedding [7381] as implemented in the Gaussian 03 [44] program suite. Such an approach enabled us to consider the role of a second arginyl (Arg50) residue as well to take into account the anisotropic environment surrounding the active site. The resulting QM/MM chemical model included the substrate URO-III and all active site residues immediately surrounding it, that is, first-shell residues (Figure 3). In addition, for those portions of the substrate exposed to solvent, the first solvation shell was retained. A subset of the complete model centered on the reactive region of the active site was then selected for the high-level QM treatment consisting of two arginyls (Arg37 and 50), an aspartyl (Asp86) residue, and the first substrate pyrrole that is decarboxylated [24]. Optimized geometries were obtained at the ONIOM(B3LYP/6-31G(d):AMBER94) level of theory [82]. However, in order to obtain more reliable calculated relative Gibbs free energies, single-point energy calculations on the above optimized structures were performed at the ONIOM(B3LYP/6-311+G(2df,p):AMBER94) level of theory. The free energies were obtained by adding the necessary energy corrections calculated at standard ambient temperature and pressure (SATP) as implemented in Gaussian 03.

Using this larger and more complete model, the initial step is again proton transfer from the guanidinium of Arg37 onto the C2 centre of the substrate pyrrole. However, this step now occurs via TS1 at a cost of 43.1 kJ mol−1 relative to the initial substrate-bound active site complex RC (Figure 4). This is approximately half the size of the barrier for this initial step obtained by Silva and Ramos [25] using a cluster/DFT approach. Experimentally, human UROD has been measured to have a value of 0.16 s−1 [5, 69, 83], corresponding to an overall barrier of 77.4 kJ mol−1. Thus, the smaller cluster/DFT approach would appear to give better agreement than the above larger and more extensive QM/MM-based approach. However, as summarized by Juárez et al. [84], decarboxylation of URO-III generating the 7-carboxylate intermediate, that is, the first decarboxylation, is most likely not the rate-limiting step. In fact, for several variants of UROD, the rate-limiting step appears to be the decarboxylation of the 7-carboxylate intermediate, that is, the second acetate decarboxylation. It is noted that for UROD from various species, the experimentally determined barriers for the first decarboxylation lie in the range of 8.4–51.5 kJ mol−1 [68, 84]. Thus, the QM/MM calculated barrier is in fact in better agreement with related experimentally reported values for the decarboxylation of ring D.

The resulting C2-protonated intermediate (I1) lies higher in energy than RC by just 6.3 kJ mol−1. Importantly, with the explicit addition of the protein environment (as opposed to the use of a PCM approach), a stabilizing effect occurs. This stabilizing effect is seen to reduce the relative energy of I1 by 41.0 kJ mol−1 in comparison to that obtained using the cluster/DFT approach in which I1 was calculated to lie higher in energy than RC by 47.3 kJ mol−1 [25].

Following protonation of C2 by Arg37, the next step is the decarboxylation of the acetate moiety (Scheme 4). Using a cluster/DFT-based approach, Silva and Ramos [25] obtained a barrier for this step of 82.4 kJ mol−1. Furthermore, based on the optimized structures obtained, they concluded that the release of CO2 was hindered by the hydrogen bond interaction between the carboxylate of Asp86 and the pyrrole ring’s HN moiety. In contrast, using a QM/MM-based approach, we found that this occurs via TS2 with a relative energy barrier lower than that of I1 (again due to inclusion of single-point energy and Gibb’s free energy corrections; see section on Multistate Reactivity PESs) [24]. That is, the loss of CO2 from the acetate essentially occurs without a barrier. Furthermore, a more complete and explicit inclusion of the protein environment enables a lengthening of the Asp86C2-protonated pyrrole interaction. In turn, this leads to a destabilization of the C2-protonated pyrrole and thus enhancement of the rate of decarboxylation [25].

The subsequent and final step is protonation of the newly formed methylene carbon from which the CO2 was lost. Unfortunately, direct comparison of the results obtained using the cluster/DFT-based [25] and QM/MM-based [24] approaches for this step is not possible. Using the former approach, the proton transferred to C′3 was proposed to originate from the solvent. In contrast, for the latter QM/MM-based study, it was found that a second arginyl active site residue (Arg50) could act as the second required mechanistic acid HB (Scheme 4) that protonates the methylene carbon. Furthermore, the barrier for this step was just 3.1 kJ mol−1 [24]. It is noted that the identification of an active site residue as HB is supported by experimental conclusions based on the observed retention of stereochemistry at the pyrrole-CH3 group formed [5961].

4. Elucidating the Catalytic Redox Mechanism and “Driving Force” in 6-Phospho-α-Glucosidase (GlvA)

In aqueous solution under standard conditions, glycosidic bonds are remarkably resistant to hydrolysis [8589]. Glycoside hydrolases (GH) are the family of enzymes that catalyse their hydrolytic cleavage within organisms and, in fact, are amongst some of the most efficient enzymes known [90100]. In general, these enzymes are highly stereospecific; however, the GH4 subfamily is able to catalyse the hydrolysis of either or both α- and β-glycosidic bonds [101, 102]. The enzyme 6-phospho-α-glucosidase (GlvA) is a member of this subfamily that uses NAD+ and a divalent Mn2+ as cofactors in its redox mechanism of glycosidic bond cleavage [97, 103, 104]. In particular, the mechanism has been proposed to proceed via two half-reactions shown in Scheme 5 [99, 105107]. The first is thought to be initiated by proton transfer from the C3-OH group to an Mn2+ bound hydroxyl and hydride transfer from C3–H to the NAD+ moiety, ultimately resulting in cleavage of the glycosidic bond at C1. The cleaved product moiety ROH is then replaced within the active site by a water molecule. The second half-reaction is believed to essentially be the reverse of the first with the H2O moiety in place of the ROH.


In order to investigate the proposed [99, 105107] catalytic mechanism of GlvA, we employed [26] a cluster/DFT-based approach in which a large active site model was used with the surrounding protein simply modeled using a PCM-based solvation method. In particular, optimized geometries were obtained at the B3LYP/6-31G(d) 5D level of theory. Relative energies were obtained at the B3LYP/6-311+G(2df,p)//B3LYP/6-31G(d) 5D level of theory (i.e., gas-phase). Corrections for the surrounding environment were obtained at the same level of theory using the integral equation formalism (IEF-) PCM method with a dielectric constant of ε = 4.0; a value typically used to model the polarity of an internal protein environment [26]. The chemical model used (Scheme 6) was derived from an X-ray crystal structure of GlvA cocrystallized with the substrate analog 6-phospho-β-glucose bound within its active site (PDB ID: 1U8X) [100]. The coordinates of proposed key mechanistic residues were then extracted, and a minimum number of atoms remote from the reaction region were held fixed at these positions in order to maintain integrity of the model during computations. Full computational details are provided in [22].


The PESs obtained for the two half-reactions are illustrated in Figure 5. The relative energies given in parenthesis were obtained in the gas phase, that is, without use of the PCM-based approach to model the general effects of the surrounding polar protein environment, and are included for comparison (see below). The catalytic mechanism of GlvA is initiated by an oxidation of the glucopyranose ring at the 3-position. Specifically, the 3′–OH group first transfers its proton to the metal bound hydroxide ligand (Mn2+OH) to give a metal-bound water and a now negatively charged 3′–O moiety on the sugar ring (IC1). This step occurs with a barrier of 22.7 kJ mol−1 when within the protein environment (“solution phase”), while the resulting intermediate IC1 lies 46.3 kJ mol−1 lower in energy than the RC. This is then followed by a hydride transfer from the 3′C–H group to the pro-R face of the C4 centre of the NAD+ [26] cofactor via TS2 at a cost of 80.8 kJ mol−1 with respect to IC1. This results in the formation of the 3-keto intermediate IC2 which lies only slightly higher in energy than IC1 by 7.2 kJ mol−1.

With the formation of IC2, a double bond has now been introduced into the substrate via the formation of a keto group. In the subsequent steps, the double bond is then effectively shifted around the ring via a series of keto-enol tautomerizations coupled with proton transfers. First, a glutamyl-tyrosyl catalytic diad deprotonates the 2′C–H group of IC2 via TS3 to form the enolate (2′C = 3′C–O) containing intermediate IC3 (Figure 5). The barrier for this process is 54.0 kJ mol−1, while IC3 lies 14.6 kJ mol−1 lower in energy than RC. This is subsequently followed by proton abstraction from the C2–OH group by an adjacent aspartyl (Asp172) residue. Simultaneously, however, the Asp172 donates its proton to the leaving CH3O moiety. Importantly, it is in this step that the actual cleavage of the glycosidic bond occurs via TS4 at cost of only 58.9 kJ mol−1 with respect to IC3. NBO and second-order perturbation analyses were used to analyse the “driving force” that leads to this greatly reduced barrier for heterolytic glycosidic bond cleavage and the regioselectivity of GlvA [22]. In particular, it was found that in the first half-reaction, electronic destabilizing effects within the sugar ring are enhanced due to several electron delocalizations that strongly favour the heterolytic cleavage of glycosidic bonds. Furthermore, this effect was greater in the case of axial (α) compared to equatorial (β) glycosidic bonds. It is noted that the first half-reaction is exothermic as the resulting product complex IC4MeOH lies 29.9 kJ mol−1, respectively, lower in energy than RC.

For the second half-reaction, the cleaved MeOH moiety must first be replaced by a water to give IC4H2O. By comparing the energies of [IC4MeOH + H2O] versus [IC4H2O  +  MeOH], it is found that this replacement is slightly exothermic by 7.3 kJ mol−1. Thus, while this exchange requires a “splitting” of the calculated overall PES, one is still able to maintain an energetic relationship between the two surfaces. With the formation of IC4H2O the second half-reaction is effectively the reverse of the first half-reaction [26]. For instance, it is initiated by nucleophilic attack of the water oxygen at the C1 centre of IC4. This occurs via TS5 at a cost of 77.7 kJ mol−1 in which the attacking H2O simultaneously transfers a proton to the Asp172 residue which itself transfers its proton onto the C2–O oxyanion centre. Similar to the analogous complex IC3, the resulting intermediate IC5 lies only slightly lower in energy than RC by 18.4 kJ mol−1. In the subsequent step, the C2 centre is reprotonated by the catalytic tyrosyl-glutamyl diad via TS6 at a cost of 79.5 kJ mol−1 to give the keto intermediate IC6 lying 39.6 kJ mol−1 lower in energy than RC. The last two steps in the overall mechanism begin with a hydride transfer from the NADH moiety onto the C3 centre via TS7, with a barrier of 72.2 kJ mol−1, to give the oxyanionic intermediate IC7. Notably, both IC7 and the preceding intermediate IC6 are essentially thermoneutral with their analogous first half-reaction complex IC2 and IC1, respectively (see Figure 5). The next and final step is a proton transfer from the metal-bound H2O onto the newly formed C3–O group via TS8 at a cost of 98.6 kJ mol−1. The overall mechanism is slightly exothermic, with the product-bound active site complex PC lying just 4.0 kJ mol−1 lower in energy than RC.

Comparison of the above “solution-phase” relative energies with the gas-phase energies enables one to gain insight the role and effect of the general protein environment. It can be seen from Figure 5 that the effect of the polar environment is not systematic. For example, inclusion of the polar environment stabilizes IC1 and TS7 by 17.3 and 42.3 kJ mol−1, respectively, but destabilizes IC2 and TS4 by 0.5 and 10.0 kJ mol−1. In general, however, larger effects are observed for transition structures; the observed absolute changes in relative energies for TSs range from 3.3–42.3 kJ mol−1 while for the intermediates and product complexes, they range from 0.5–27.3 kJ mol−1.

More importantly, however, the inclusion of the general protein environment can have marked effect on the catalytic mechanism. For example, in the gas phase the overall mechanism is endothermic with the PC lying 11.4 kJ mol−1 higher in energy than RC. However, in the “solution-phase”, it is now slightly exothermic by 4.0 kJ mol−1 (Figure 5). Moreover, with the inclusion of the protein environment the rate-limiting step has changed. Specifically, in the gas phase, the proton transfer from the Mn2+-bound water to the 3′C–O group via TS8 represents the rate-limiting step with a reaction barrier relative to RC of 78.3 kJ mol−1. However, with inclusion of the protein environment (i.e., use of a PCM-based approach), the energy of TS8 is reduced significantly by ~26.2 kJ mol−1 to 52.1 kJ mol−1 respect to RC. As a consequence, TS6 (61.1 kJ mol−1) now has the highest energy relative to RC. Thus, reprotonation of C2 by the catalytic tyrosyl-glutamyl diad has become the rate-limiting step. In agreement with experiment, glycosidic bond cleavage is not the rate-controlling step.

5. Summary

Potential energy surfaces (PESs) are powerful tools for elucidating the catalytic mechanisms of enzymes. Indeed, they can be more than a simple indication of a path of sequential mechanistic structures and their energetic relationships. In this paper the application of the tools and techniques of computational chemistry in combination with PES investigations has been described. Through a review of several of our recently published studies on enzymatic systems, how such an approach may be applied in order to provide deeper insights into the movement of electrons during the course of a reaction, the mechanistic role of active site residues and the surrounding environment, and the utility of multistate PESs in metalloenzymes has been highlighted. In addition, some of the challenges that these studies may encounter, such as appropriate choice of computational model on the reliability and accuracy of a PES, have also been discussed.

For example, the catalytic activation of O2 by AlkB to give the highly reactive oxo- and oxyl-type intermediates was examined for several possible spin states. It was found that at least in the initial steps of the process, multiple states were close in energy with each other and, furthermore, within the enzymatic thermochemical limit. As a result, and considering the likelihood of spin changes in transition metal chemistry, several minimum energy pathways to give the lowest energy product complex were identified, each involving spin inversion (SI) from a septet to quintet PES [23]. Importantly, SI enables the catalytic mechanism of AlkB to access and proceed via highly exoergic intermediates.

For the enzyme UROD, a comparison was made of its catalytic mechanism as elucidated using a combined QM and MM method in the ONIOM formalism with mechanical embedding and a cluster/DFT+PCM approach. [25] Use of a QM/MM approach enables one to more directly and extensively model the anisotropic environment surrounding an enzyme active site. In the case of UROD, with the QM/MM-based approach, it was found that the stabilizing effect of the surrounding environment provided a significant reduction in both the calculated rate-limiting proton transfer as well as the barrier for acetate decarboxylation. In fact, the more accurate inclusion of the environment enabled us [19] to obtain an energy for the rate-limiting barrier in good agreement with related experimentally determined values.

The applicability of cluster/DFT-based approaches to the study of enzymatic systems was examined further using the case of the glycosidic hydrolase GlvA. In particular, such an approach was shown to be able to provide reliable and accurate insights when used in combination with a well-chosen chemical model. Inclusion of the general effects (polarity) of the environment surrounding the active site was shown to be important in choosing a suitable model. Indeed, its inclusion had significant effects on the mechanism obtained such as shifting it from an overall endothermic to exothermic process. More importantly, however, the rate-limiting step changed. Specifically, in the gas phase, a proton transfer from a metal-bound water to the 3′C–O group represented the rate-limiting step. However, with the inclusion of the polar protein environment using a PCM-based approach, the rate-limiting step was instead protonation of the sugar ring’s C2 centre by the catalytic tyrosyl-glutamyl diad.


The authors gratefully acknowledge the Natural Sciences and Engineering Research Council (NSERC, Canada) for funding and CGS & PGSD scholarships (E. A. C. Bushnell). SHARCNET (Canada) is acknowledged for a graduate scholarship (W.-J. Huang) and additional computational resources.