DFT calculations for the equilibrium isotope effect for deuterium substitution at the anomeric centre Cα in 2-(p-nitrophenoxy)tetrahydropyran with continuum solvation show significant variation in the range of relative permittivity 2 ≤ ε ≤ 10. One-dimensional scans of potential energy (with implicit solvation by water) or of free energy (from QM/MM potentials of mean force with explicit aqueous solvation with a hybrid AM1/OPLS method) for heterolysis of the bond between Cα and the nucleofuge do not show a transition state. A two-dimensional free-energy surface that considers also the distance between Cα and a nucleophilic water indicates a pre-association DN*ANint‡ mechanism with a transition state involving nucleophilic attack upon an ion-pair intermediate, and this is supported by good agreement between the mean values of the calculated and experimental α-D KIEs. However, the magnitudes of the standard deviations about the mean values for the making and breaking C–O bonds suggest that the transition state is rather plastic, with Cα–Onu≈2 ± 0.4 Å and Cα–Olg≈3 ± 0.5 Å. Not only is nucleophilic solvent assistance necessary, but there is also evidence for electrophilic assistance through specific hydrogen bonding to the nucleofuge.
Kinetic isotope effects (KIEs) are a powerful experimental tool for transition-state (TS) analysis of mechanisms of reactions in chemistry and biochemistry, provided that their values may be meaningfully interpreted. It is often assumed that trends in KIEs may be related to structural changes in TSs, understood in terms of geometries and bond orders, but recent computational studies have suggested that changes in the electrostatic environment of a TS may be significant , . If this is true, then interpretation of experimental KIEs for enzyme-catalysed reactions ought to include consideration of the electrostatic environment of the TS as well as of the geometrical changes accompanying bond making/breaking and rehybridisation around the reaction centre.
Glycosides are ubiquitous in biology , from the “bulk chemicals” aspect of organismal cell wall materials and biomass conversion through to the “fine chemicals” aspect of post-translational modification of proteins critical to specific cellular biological functions, and thus, the enzymes (glycosyltransferases) that control these chemical processes are immensely important . In particular, glycoside hydrolases, which catalyse the transfer of a glycosyl unit to water, thereby breaking polysaccharides into smaller fragments, have both historical and topical relevance. Unlike simple acetals, glycosides – which are acetals derived from carbohydrates – are usually rather unreactive towards hydrolysis under all pH conditions, and require enzymic catalysis, which may occur with retention or inversion of configuration at the anomeric carbon; the mechanisms by which these enzymes operate remains a subject of great interest and importance .
The simple model substrate 2-(p-nitrophenoxy)tetrahydropyran 1 does undergo uncatalysed pH-independent hydrolysis over a wide pH range, by virtue of possessing a good nucleofuge (leaving group) but lacking the strongly deactivating hydroxyl substituents of carbohydrate-derived acetals. There is general consensus for a DN+AN (or SN1) mechanism in non-acidic aqueous solution: the C–O bond between the anomeric carbon and the oxygen of the p-nitrophenoxide leaving group of 1 heterolyses, giving an oxacarbenium intermediate 3, which then undergoes nucleophilic attack by water, yielding the hemiacetal product. Maskill and co-workers  determined the α-deuterium secondary kinetic isotope effect (α-KIE) for the uncatalysed hydrolysis of 1 to be 1.17 in water at 46°C (corresponding to about 1.19 at 25°C). Since this value was appreciably higher those that had been typically reported for enzyme-catalysed glycolysis , , , it was interpreted as indicating that there may be less carbenium-ion character in the TS for rate-limiting step in enzymic glycolysis than had been previously thought.
As the initial stage of a larger programme of research aimed at performing computational simulations of mechanisms for enzymic and non-enzymic glycolysis making reliable predictions of kinetic isotope effects, we now present preliminary results for the uncatalysed hydrolysis of 1. We have studied the dependence of the α-deuterium equilibrium isotope effect (EIE) for heterolysis of aqueous 1 yielding intermediate 3 in a dielectric medium of lower relative polarity, using density-functional theory and the polarisable continuum model (PCM) for solvation. The purpose of this was to check what the extent to which the isotope effect depends on the change in the electrostatic environments of 1 and 3 (or of the TS 2 leading to the intermediate). We have also employed a hybrid quantum-mechanical/molecular-mechanical (QM/MM) method with molecular dynamics (MD) to generate a two-dimensional free-energy surface for the hydrolysis of 1 in water.
Heterolysis of 1 in implicit solvent
Geometry optimization and harmonic frequency calculations were performed with Gaussian09 (revision A.02) . The IEFPCM solvation model was used with the B3LYP density functional , , ,  and the aug-cc-pVDZ basis set , as employed in our previous work , . The default scaled van der Waals surface solute cavity was used with Universal Force Field (UFF) atomic radii  for C, O and H; all other solvent-cavity parameters were set to the default values for water, but the value of the relative permittivity was varied.
The α-D EIE was evaluated as the ratio of isotopic partition function (IPFR) ratios fε (eq. 1) for 1 in water (ε=80) and 3 in a continuum solvent (2≤ε≤80) using the harmonic-oscillator, rigid-rotor approximation with standard expressions for translational, rotational and vibrational partition functions (eq. 2) as previously described , , with unscaled frequencies. The prime (′) in eq. 2 denotes the heavy isotopologue L=D (Scheme 1) as opposed to the unprimed L=H; |I| is the determinant of the moment of inertia tensor, M is the molecular mass, N is the number of atoms, v is a vibrational frequency, u=hv/kBT where h, kB and T are the Planck and Boltzmann constants and the absolute temperature.
The α-D KIE in PCM  water was estimated according to eq. 3 for 1 and an approximate TS, for which the IPFR fTS includes a quantum correction for the imaginary frequency (eq. 4).
A scaling factor equal to 0.978 was applied to the B3LYP/aug-cc-pVDZ harmonic frequencies, as previously reported ; its value was by comparison of 24 vibrational frequencies for four isotopologues of the methyl cation (CH3+, CH2D+, CHD2+ and CD3+) calculated by means of the vibrational configuration interaction CCSD(T*)-F12a/cc-pVQZ-F12 method which includes anharmonicity.
Hydrolysis of 1 in explicit water
QM/MM molecular dynamics, energy minimization and Hessian-evaluation calculations were performed by means of the fDynamo library . The initial structure of 1 was calculated using the B3LYP/aug-cc-pVDZ functional and basis set, which was then added to a 31.4 Å cubic box of 1018 pre-equilibrated water molecules. Any waters with an oxygen atom lying within a radius of 2.8 Å from an atom of acetal 1 were removed. The geometries of the remaining water molecules were then optimized. The QM system for one-dimensional treatment contained acetal 1 and five water molecules (one as a nucleophile and the others making hydrogen bonds to the nucleofuge). The QM system for two-dimensional treatment contained acetal 1 and one water molecule as a nucleophile. The QM region was described by the semi-empirical AM1 method. The MM subsystem comprised the remaining solvent water molecules and was described by the TIP3P potential within the OPLS-AA force field . The whole system was then equilibrated by MD at 300 K in the NVT ensemble with periodic boundary conditions.
Progression of the reaction was described by two reaction coordinates, defined as the difference between the distances between Cα and the oxygen atoms of the nucleofuge and the nucleophile: RC1=d(Cα–Olg) and RC2=d(Cα–Onu). A one-dimensional potential of mean force (1D-PMF) was constructed by averaging the atomic motions of the whole system at 300 K, as governed by the QM/MM potential energy function and forces for QM 1+5 H2O surrounded by 1010 MM water molecules, in each of 64 simulation windows, using umbrella sampling with a harmonic constraint of 2500 kJ mol−1 Å−2 applied to the coordinate in the range 1.4≤RC1≤4.6 Å; using a time-step of 1 fs in Langevin-Verlet MD, 15 ps of equilibration was followed by a 15 ps production trajectory, from which the PMF was obtained by means of the weighted-histogram analysis method (WHAM) . Similarly, a 2D-PMF was constructed with respect to RC1 and RC2 using double umbrella sampling in 36×24 simulation windows in the ranges 1.4≤RC1≤4.0 Å and 1.4≤RC2≤3.0 Å, as previously described ; the QM region comprised 1 and one water molecule as the nucleophile. The resulting surface provides information about the change in Helmholtz free energy as a function of both exocyclic Cα–Olg bond breaking and Cα–Onu bond making; the difference between this and the Gibbs energy surface is negligible, and hereafter we refer simply to relative free energy.
The α-D KIE for hydrolysis of 1 (AM1) in explicit water was obtained as the quotient of average IPFRs for the RS 1 and TS 2. Thirty-two snapshot configurations from a constrained MD simulation (RC1=3.0 Å and RC2=2.0 Å) were subjected to relaxation of the QM region to a local minimum (for 1) or first-order saddle point (for 2) within a frozen MM solvent environment, and a subset hessian was evaluated for the Ns QM atoms in each individual configuration. The IPFR for each configuration was evaluated using the harmonic approximation for all 3Ns subset degrees of freedom, according to eqs. 6 and 7, where the brackets denote averages over thermally accessible configurations as previously described , .
Results and discussion
Influence of dielectric environment of 3 on α-D EIE for heterolysis of aqueous 1
Figure 1 shows a plot of the natural logarithm of the calculated EIE (=KH/KD) for deuterium substitution on the anomeric carbon in 1 at 25°C for heterolysis of aqueous 1 yielding the oxacarbenium ion 3 in a continuum dielectric environment of varying relative permittivity ε. All values of the EIE are normal, as expected for the hybridisation change from sp3 in 1 to sp2 in 3. The function (1–1/ε) on the abscissa is a factor in the Born model for the electrostatic contribution the free energy of solvation of a point charge within a spherical cavity . Since the logarithm of an isotope effect (either a KIE or an EIE) is proportional to a change in free energy, this plot represents an example of a novel type of free-energy relationship. The most important point is that the magnitude of the EIE varies with the relative permittivity of the environment surrounding the cationic species 3; it changes from ln(EIE)=0.211 (KH/KD=1.235) for ε=80 (i.e. 1–1/ε=0.9875) to ln(EIE)=0.230 (KH/KD=1.258) for ε=2 (i.e. 1–1/ε=0.5); the form of the variation in the plot is described very well by a quadratic function. Note that more than 80% of the total variation in ln(EIE) occurs in the range of low polarity, 2≤ε≤10 (i.e. 0.5≤ 1–1/ε≤0.9), usually associated with the effective value of the relative permittivity in an enzyme active site .
The EIE for the heterolysis (1)water→(3)hexane (or 3 in any other low-polarity solvent) is unlikely to be studied experimentally any time soon, but the significance of this computational study is that it provides a simple model for an enzyme-catalysed glycolysis. The TS for this reaction would be expected to closely resemble the oxacarbenium ion, and therefore the KIE should be close in value to the calculated EIE. If an experimental KIE were measured under non-saturating conditions as the isotope effect on kcat/Km, then the reactant state (RS) would be the free glycoside in water whereas the TS would be an oxacarbenium-like species in an enzyme active site characterised by a low value of the effective relative permittivity. TS-stabilising interactions between the oxacarbenium-like species and the enzyme might influence the TS structure, perhaps in accordance with a Hammond-postulate prediction, shifting its position along the reaction coordinate, and thereby affecting the value of the KIE. This type of “solvent effect” upon a KIE has been discussed previously, for example, for the Menshutkin reaction in a range of solvents . Solvent variation may also change a mechanism, for example by altering the identity of a rate-limiting step, as discussed by Shiner for solvolytic substitutions . However, to the best of our knowledge, the only previous recognition that changing the solvent might also influence an isotope effect by affecting isotopically-sensitive vibrational frequencies directly, regardless of TS structure or the identity of a rate-determining step, was by Keller and Yankwich , , , in their largely-overlooked model studies of medium effects on heavy-atom KIEs; otherwise this possibility seems to have been totally neglected by the physical-organic chemistry community until now.
KIEs offer powerful probes for mechanism and TS structure in enzyme-catalysed reactions if their experimentally determined values and variations can be correctly interpreted. It is common for analyses of TS structure, based upon KIEs for multiple isotopic substitutions, to consider force-constant changes only as functions of molecular geometry . However, the present work suggests that consideration of the electrostatic environment might also be necessary, as force-constant changes may depend to some extent on the relative permittivity of the medium. It is conceivable that a reaction involving charge redistribution, separation or neutralization within an enzyme active site could manifest variations in KIEs, as between a wild-type and a mutant form of the enzyme, that originate from changes in the effective dielectric within the active-site environment; if so, there would be important implications for the interpretation of experimental KIEs in mechanistic enzymology.
Potential energy profile and estimated KIE for heterolysis of aqueous 1
Figure 2 shows the potential energy (PE) profile for heterolysis of aqueous 1 as a function of the exocyclic C–O bond length, calculated by means of the B3LYP/aug-cc-pVDZ/PCM method with ε=80. The dashed line indicates the asymptote for separate and independently solvated oxacarbenium cation 3 and p-nitrophenoxide anion. Clearly there is no indication of a maximum in this profile corresponding to a transition structure for formation of an ion pair. This means that a true KIE cannot be calculated between 1 and a genuine TS for heterolysis in PCM water. However, it is possible to consider a range of non-stationary points on the profile in the region of negative curvature and to calculate the IPFR for each (using eqs. 3 and 4) as if each were indeed a TS. These estimated KIEs (Table 1) may then be compared with the experimental value of 1.17 for kH/kα-D at 46°C . Using unscaled frequencies in the IPFR calculations, the best agreement is obtained for an exocyclic C–O bond length equal to about 2.23 Å, which therefore provides an estimate of the extent of bond breaking in the approximate TS. Using scaled frequencies, best agreement is obtained at C–O≈2.25 Å. The ratio of natural logarithms of the estimated KIE and the calculated EIE gives ln(1.171)/ln(1.210)=0.83, which may be considered as a rough measure of the position of the approximate TS along the reaction coordinate for heterolysis. However, these estimates should be regarded with caution in view of the limitations of the PCM method for solvation and the absence of a maximum in the 1D PE profile.
Potential of mean force, free-energy surface and KIE for hydrolysis of aqueous 1
The open squares in Fig. 2 show the AM1/PCM potential-energy profile for heterolysis of 1 in implicit water whereas the apparently solid black line (which is actually made up of closely-spaced dots) shows the AM1/OPLS free-energy profile for this process in explicit water. The two profiles are essentially the same until Cα–Olg≈1.9 Å, but for longer distances the free-energy profile plateaus out whereas the potential-energy profile tends towards the asymptotic limit of infinitely separated cation and anion in the continuum solvent. The free-energy profile has been extended as far as Cα–Olg=9.4 Å but shows no obvious maximum that would indicate a TS for formation of any sort of ion pair as an intermediate in the overall hydrolysis.
Figure 3 shows the 2D AM1/OPLS free-energy surface for the hydrolysis of 1 as a function of the coordinates RC1 and RC2, which describe Cα–Olg bond breaking and Cα–Onu bond making, respectively. The energy well on the left corresponds to 1 and that on the right to the protonated hemiacetal initial product, denoted as A.
A cross-section of the 2D free-energy surface along the front left-hand edge (denoted by a black line), corresponding to elongation of the Cα–Olg bond while the nucleophile remains at about 3 Å, closely resembles the 1D free-energy profile in Fig. 2, with no distinguishable maximum present along this reaction coordinate. However, there is a shallow channel of lower free-energy for structures (B) with Cα–Onu≈2.3 Å, which suggests a degree of nucleophilic solvent-assistance during Cα–Olg heterolysis, corresponding to an SN2 (intermediate) or DN*ANint mechanism , . Previous experimental studies of this reaction have suggested a degree of solvent participation in the TS; however, this was attributed to assistance of the nucleofuge via hydrogen bonding . Inspection of structures on the free-energy surface in the vicinity of the “shoulder” from the reactant valley onto the plateau region (Cα–Onu≈2.5 Å, Cα–Olg≈2.0 Å) does indeed suggest a specific hydrogen bond between a water molecule and the Olg atom. The additional stabilisation that this electrophilic assistance provides to the nucleofuge probably explains the divergence (noted above) between the PCM potential-energy profile with implicit solvation (which cannot describe specific hydrogen bonds) and the QM/MM free-energy profile with explicit solvation.
A free-energy barrier at Cα–Onu≈2.0 Å clearly separates channel B from product valley A. The free energy of structures in the vicinity of C on the free-energy surface is about 110 kJ mol−1 relative to 1. This AM1/OPLS estimate is in good agreement with the experimental value of ΔG‡=100±9 kJ mol−1 reported by Maskill and co-workers . Figure 4 shows the structures of points of interest on the free-energy surface: (a) a protonated hemiacetal+unbound nucleofuge, from region A; (b) a putative SN2 (intermediate) structure, from region B, showing electrophilic assistance for nucleofuge departure; (c) a putative transition structure, from region C. The hybridisation changes at Cα are evident: from sp3 in 1, to sp2 in B and C, and back to sp3 in A.
A sample of 32 snapshot configurations, from constrained MD simulations with the AM1/OPLS method, in the vicinity of region C was used to obtain an average IPFR for deuterium substitution at the anomeric position in the TS which, combined with the average IPFR for the RS, gave the α-D KIE values shown in Table 2 at both 25 and 46°C. The average TS bond length to the nucleofuge was ⟨Cα–Olg⟩=3.02±0.54 Å (cf. 1.44±0.01 Å in RS) and to the nucleophile was ⟨Cα–Onu⟩=1.97±0.38 Å. The relatively large standard deviations about the mean values of these bond distances reflects the plasticity of the TS and its sensitivity to the dynamic solvation environment in water. Furthermore, this gives rise to the relatively large standard deviation about the mean value of the IPFR for the TS at each temperature and consequently to the much larger standard deviation on the calculated KIE than the quoted error for the experimental KIE at 46°C. The experimental KIE at 25°C was, of course, estimated in the usual manner by means of the experimental activation enthalpy, and is included here for ease of comparison with other published α-D KIEs. The agreement between the mean values of the calculated and experimental KIEs is very satisfying, as it lends support to a DN*ANint‡ mechanism with a TS involving nucleophilic attack upon an ion-pair intermediate but comprising an ensemble of individual transition structures spanning a range of geometries with Cα–Onu≈2±0.4 Å and Cα–Olg≈3±0.5 Å.
The magnitude of the B3LYP/aug-cc-pVDZ/PCM calculated EIE for deuterium substitution at Cα during heterolysis of aqueous 1, yielding the oxacarbenium ion 3 in a continuum dielectric environment, varies with the value of the relative permittivity, particularly in the range of low polarity (2≤ε≤10) usually considered to describe an enzyme active site. This suggests that it is necessary to consider not only geometrical changes but also changes in the electrostatic environment in studies of reactions involving charge redistribution, separation or neutralization within an enzyme active site, and that variations in KIEs, as between a wild-type and a mutant form of the enzyme, could originate from changes in the effective dielectric within the active-site environment; if so, there would be important implications for the interpretation of experimental KIEs in mechanistic enzymology.
A 1D PE scan at the B3LYP/aug-cc-pVDZ/PCM level with respect to Cα–Olg bond breaking does not indicate the existence of a TS for heterolysis of aqueous 1, but comparison of the approximate calculated α-D KIE with experiment suggests this might occur at distance of about 2.23–2.25 Å. The AM1/OPLS 1D-PMF profile with explicit aqueous solvation follows the corresponding 1D PE scan with respect to Cα–Olg bond breaking until about 1.9 Å, after which additional stabilisation arises due to specific hydrogen bonding of a water molecule to the Olg atom of the nucleofuge.
The AM1/OPLS 2D-PMF with explicit aqueous solvation shows a broad and rather flat plateau containing a shallow channel with a QM nucleophilic water molecule at a distance Cα–Onu≈2.3 Å for a range of Cα–Olg distances between 2 and 4 Å (at least). This region of pre-associated species is separated from the (protonated) hemiacetal product valley by a ridge characterised by Cα–Onu≈2.0 Å. This free-energy surface suggests neither a DN+AN (or SN1) mechanism nor an ANDN (or SN2) mechanism but rather a pre-association DN*ANint‡ mechanism with a TS involving nucleophilic attack upon an ion-pair intermediate. This is supported by the good agreement between the mean values of the calculated and experimental α-D KIEs. However, the magnitudes of the standard deviations about the mean values for Cα–Onu and Cα–Olg suggest that the TS is rather plastic, with Cα–Onu≈2±0.4 Å and Cα–Olg≈3±0.5 Å. Not only is nucleophilic solvent assistance necessary, but there is also evidence for electrophilic assistance through specific hydrogen bonding to the nucleofuge.
Further work with higher-level corrections to the QM/MM calculated 2D PMF ,  is in progress to determine whether the mechanism suggested by the topography of the AM1/OPLS free-energy surface is confirmed by use of a DFT description of the QM region and by more extensive sampling of solvent configurations. This includes the refining of stationary points of interest using DFT methods; however, this is computationally demanding and beyond the scope of the present study.
A collection of invited papers based on presentations at the 24th IUPAC International Conference on Physical Organic Chemistry (ICPOC 24) held in Faro, Portugal, 1–6 July 2018.
Funding source: EPSRC
Award Identifier / Grant number: EP/L016354/1
Funding statement: This work was supported by the EPSRC Centre (Funder Id:http://dx.doi.org/10.13039/501100000266, Grant Number: EP/L016354/1) for Doctoral Training in Sustainable Chemical Technologies (Studentship to J.H.G.) and made use of the Balena High Performance Computing (HPC) Service at the University of Bath.
 J. Kamerling, Ed. Comprehensive Glycoscience: From Chemistry to Systems Biology, Elsevier, Amsterdam (2007).Search in Google Scholar
 G. Davies, M. L. Sinnott, S. G. Withers. in Comprehensive Biological Catalysis: A Mechanistic Reference, M. Sinnott (Ed.), Academic Press, San Diego (1997).Search in Google Scholar
 I. A. Ahmad, S. L. Birkby, C. A. Bullen, P. D. Groves, T. Lankau, W. H. Lee, H. Maskill, P. C. Miatt, I. D. Menneer, K. Shaw. J. Phys. Org. Chem. 17, 560 (2004).10.1002/poc.803Search in Google Scholar
 M. 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, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, D. J. Fox. Gaussian 09, Revision A.02. Gaussian, Inc., Wallingford, CT (2009).Search in Google Scholar
 I. H. Williams. J. Mol. Struct., Theochem. 94, 275 (1983).Search in Google Scholar
 M. J. Field, M. Albe, C. Bret, F. Proust-De Martin, A. Thomas. J. Comput. Chem. 21, 1088 (2000).10.1002/1096-987X(200009)21:12<1088::AID-JCC5>3.0.CO;2-8Search in Google Scholar
 V. J. Shiner. in Isotope Effects in Chemical Reactions, C. J. Collins, N. S. Bowman (Eds.), Van Nostrand Reinhold, New York (1970).Search in Google Scholar
 P. J. Berti, K. S. E. Tanaka. Adv. Phys. Org. Chem. 37, 239 (2002).Search in Google Scholar
© 2020 IUPAC & De Gruyter, Berlin/Boston