Theoretical Study of HOCl-Catalyzed Keto-Enol Tautomerization of \xce\xb2-Cyclopentanedione in an Explicit Water Environment
The mechanism of acid-catalyzed keto-enol tautomerization of \xce\xb2-cyclopentanedione (CPD) in solution is studied computationally. Reaction profiles are first calculated for a limited solvation environment using ab initio and density functional methods. Barrier heights for systems including up to one hydration shell of explicit water molecules depend strongly on the number of waters involved in proton transfer, and to a lesser but significant extent on the number of waters forming hydrogen bonds with waters in the proton transfer chain (each such water reduces the barrier by 4.4 kcal/mol on average). Barriers of 8-13 kcal/mol were obtained when a full or nearly full hydration shell was present, consistent with calculations for non-acid-catalyzed keto-enol tautomerization of related molecules. The presence of HOCl reduced the barrier by 4.5 kcal/mol viz-a-viz the gas phase, consistent with the well-known principle that keto-enol tautomerization can be acid or base catalyzed. Reaction was also modeled beginning with snapshots of reactant conformations taken from a 300 K molecular dynamics simulation of CPD, HOCl and 324 explicit waters. Reaction profiles were calculated at a QM/MM level with waters in the first hydration shell either fixed or energy minimized at each step along the reaction coordinate. A substantial variation in barrier height was observed in both cases, depending primarily on electrostatic interactions (hydrogen bonding) with first-hydration-shell waters and to a lesser extent on electrostatic interactions with more distant waters and geometric distortion effects. For the lowest barriers, the extent of barrier reduction by waters involved in proton transfer is consistent with the limited-solvation results, but further barrier reduction due to hydrogen-bonding to waters involved in proton transfer is not observed. It is postulated that this is because highly flexible structures such as extensive hydrogen bonding networks optimal for reaction are entropically disfavored and so may not contribute significantly to the observed reaction rate.
INTRODUCTION<!>COMPUTATIONAL METHODS<!>Reaction in limited solvation environment<!>Reaction in full solution<!>MD simulation and selection of frames<!>Determination of minimum energy path<!>Analysis of barriers<!>CONCLUSION
<p>Keto–enol tautomerization plays a fundamental role in many organic and biochemical reactions, as the activated C-C double bond of the enol form is a useful target for nucleophilic attack. The thermodynamics of reaction has been studied extensively with both experiment1-10 and theory,8,10-31 which are generally in agreement; either of the two tautomers may be favored depending on the compound and solvation conditions (more extensive reference lists can be found in Refs. 23 and 28). The role of water molecules and acid catalysts in the tautomerization mechanism has been studied less extensively. Theoretical studies have shown that the reaction barrier varies substantially with the number and placement of participating waters.21,23-29 However, despite significant advances,21,23-27 the level of solvation required to realistically model keto-enol tautomerization has not been fully established. It is not known whether extensive hydrogen bonding water networks that can significantly reduce the barrier in small water clusters are highly sampled in solution at thermal equilibrium. This paper addresses these aspects of keto-enol tautomerization. Previous work on limited solvation evaluated the enthalpic contribution to the tautomerization reaction barrier for reactions involving proton transfer chains of up to four waters.21,23-27 We focus on proton transfer chains up to two waters in length and compare the influence of solvation in a limited solvation environment and in a fully solvated environment representative of thermal equilibrium at 300 K. We expect the results to be relevant for proton transfer chains of various lengths.</p><p>There is special interest in β-diketones because of their extensive use as organic reagents. Combining this interest with our research group's interest in chloroperoxidase (CPO),32-34 we have focused on β-cyclopentanedione (CPD), a substrate for the chlorination process catalyzed by CPO. We have modeled the acid-catalyzed tautomerization process using hypochlorous acid (HOCl) as the acid catalyst in order to mimic enzymatic reaction conditions; one postulated mechanism for the enzymatic reaction involves catalytic formation of HOCl at the active site of CPO followed by release into solution, where it halogenates CPD.35 Of the two acid species present, hydronium (strong acid, lower concentration) and HOCl (weak acid, higher concentration), we chose HOCl for initial study to avoid the technical and interpretive complications of using a charged species in the calculations, which include small hydrated clusters. Future work will evaluate the use of hydronium for comparison.</p><p>Yamabe et al.23 developed a useful conceptual framework for understanding the role of water in keto-enol tautomerization, which we outline here. As a preliminary note, we will use the terms reactive and solvating to indicate water molecules explicitly involved in proton transfer, and stabilizing but not involved in proton transfer, respectively. Yamabe et al. calculated the barrier for keto-enol tautomerization of small, linear β-diketones, including up to five explicit water molecules. They found that reactive waters reduce the reaction barrier for malonaldehyde from 60 to as little as 13 kcal/mol, and that solvating waters reduce the barrier further to a little as 7 kcal/mol. The determining factor for maximal barrier reduction was the formation of hydrogen bond networks with near-optimum geometry (OH…O distance less that 2 Å and OHO angle near 180°), which was achieved with networks of two or three reactive waters. The energetically most costly step is breaking the C-H bond during nucleophilic attack by a proximal water oxygen (Scheme 1). Yamabe et al. attributed the reduction in barrier primarily to an increase in the nucleophilicity of the attacking water oxygen resulting from the water hydrogen bond network. Stabilization of developing charge on the attacking water by the water hydrogen bond network as the transition is formed provides a complementary explanation.</p><p>Other work on β-diketones and related molecules is largely consistent with the observations of Yamabe et al.26,27 In general, inclusion of a single reactive water lowers the barrier by 30-40 kcal/mol for typical molecules, a second reactive water lowers the barrier a further 8-18 kcal/mol, and each of the first two solvating waters can lower the barrier by up to ~7 kcal/mol; diminishing effects are observed for further incremental solvent waters.21,23,26 Cucinotta et al. applied an ab initio metadynamics treatment to keto-enol tautomerization of acetone in an aqueous solution consisting of 28 water molecules.24,25 They observed a Grotthuss mechanism involving four water molecules in a proton transfer chain and a free energy barrier of 38.5 kcal/mol, as compared to 57.7 kcal/mol for intramolecular proton transfer in the gas phase.</p><p>With the exception of the metadynamics work, previous studies focused on optimized structures with up to five waters distributed between reactive and solvating roles.23,26,27 It has yet to be established whether this degree of hydration is sufficient to reproduce solution-phase barrier heights or reaction mechanisms, or indeed whether it is possible to do so within a minimum energy path (MEP) conceptual framework. Fully optimized water hydrogen bond networks might represent events so extremely rare that they are not predictive of the observed kinetics in solution at ambient temperatures. We address these issues with two studies on CPD. In the first, we systematically increase the number of explicit waters in fully optimized structures and determine the barrier and the reaction energy as functions of the number of waters. We find that the barrier to reaction decreases as the number of waters participating in the proton transfer chain increases, by 36 kcal/mol for the first such water and by 12 kcal/mol for the second. We also find that the barrier decreases as the number of waters forming hydrogen bonds with waters involved in the proton transfer chain increases, by approximately 4.4 kcal/mol for each such water. We find that the presence of the weak acid HOCl decreases the barrier modestly, by 4.5 kcal/mol. In the second study, we use classical molecular dynamics (MD) to generate solution environments representative of thermal equilibrium at 300 K. We then determine reaction profiles at a QM/MM level beginning with snapshot structures of reactants obtained from the MD simulation. We find that changes of the instantaneous water environment result in a significant variation of the reaction barriers. Barriers ranging from 25 to 44 kcal/mol are obtained for ten snapshot structures, with the variation depending primarily on electrostatic interactions (hydrogen bonding) with first-hydration-shell waters and to a lesser extent on electrostatic interactions with more distant waters and other influences mediated by geometric distortion of the reactant CPD. Little stabilization of transition state structures by hydrogen bonding to solvent waters is observed. Within the limits of the QM/MM interaction representation, this suggests that the extensive, optimized hydrogen bonding networks formed in minimum-energy-path studies on limited-solvation systems may have such a high entropic cost that they do not contribute significantly to observed reaction rates. A quantitative assessment will be necessary to address this issue definitively but is beyond the scope of this paper.</p><!><p>Reaction enthalpies in the gas phase and in explicit water clusters were calculated using a modified G3(MP2)//B3LYP method.36 The G3 method37 was introduced to enable calculations of high accuracy while keeping the requirement for computational resources to a minimum. The G3(MP2)//B3LYP method reduces the computation time further while maintaining very similar accuracy over a data base of 299 reactions. In this method, optimized geometries and zero point energies (ZPE) are determined at the B3LYP/6-31+G(d,p) level.38,39 Single point energies at the QCISD/6-31(d) level are modified with a correction scheme to adjust for basis set deficiency and ZPE using the formula:
E[G3(MP2)]=E[QCISD∕6−31G(d)]+ΔEMP2+ZPE
where ΔEMP2 = E[MP2/G3MP2large] – E[MP2/6-31G(d)] is the basis set correction. The standard G3(MP2)//B3LYP scheme also includes a spin-orbit correction and a higher-level correction, but these were not included because they are expected to be very small for our closed-shell, molecular systems. The calculations were performed with the program Gaussian 09.40 Hereafter G3 will be used to refer to this G3(MP2)//B3LYP variant for conciseness. To obtain transition states, an appropriate coordinate, generally a C-H or O-H distance involved in proton transfer, was first scanned using constrained energy minimization. Transition-state optimization methods were applied to the geometry of the energy maximum. Particularly with higher levels of explicit solvation, it was tricky to locate the region of appropriate curvature, and manual adjustment of several atoms was sometimes necessary. Ultimately, all transition states were verified to be first-order saddle points by checking that precisely one harmonic vibrational frequency was imaginary. Single point energies for stationary points of the gas-phase and water-cluster systems were also calculated with the polarized continuum model (PCM)41,42 as implemented in Gaussian, using a dielectric constant representative of water (ε = 78.4). This method represents the extended solvent environment implicitly as a dielectric medium. Geometry optimizations were performed with PCM in a few sample cases; these gave barriers and geometries similar to those obtained with the single-point calculations.</p><p>To study reaction in a fully solvating environment of explicit waters at room temperature, a two-step approach was used. First, classical molecular dynamics was carried out with the program CHARMM43 and the CHARMM22 force field44 on a system containing one molecule each of CPD and HOCl in a truncated octahedral box of waters of side 20.78 Å with periodic boundary conditions. The water molecules were represented with the TIP3P model45 as modified for the CHARMM force field.44,46 Molecular mechanics parameters were developed for CPD and HOCl using the protocol outlined by Mackrell et al.47 and are provided in the Supporting Information (hereafter Supp. Info.). The effective concentration of HOCl was larger than in experiment48-50 (0.17 M vs 0.002 M) in order to efficiently sample conformations in which HOCl and CPD are close to each other. The higher concentration is not expected to significantly affect the water structure in the vicinity of CPD and thus should not affect the calculated reaction profiles. The long range Coulomb interactions were calculated using the particle-mesh Ewald method,51 with a cutoff of 10 Å for real-space interactions and a 1 Å grid with sixth-order B-spline interpolation for reciprocal-space interactions. Lennard-Jones forces were treated with the force switch method52 using a switching range of 8-10 Å. The lengths of all bonds involving hydrogens were constrained using the SHAKE algorithm.53 The simulation was carried out in the NPT ensemble at a constant temperature of 300 K and constant atmospheric pressure of 1 atm using an extended-system algorithm54 implemented in CHARMM (keyword CPT with option PCON) and a 1 fs time step. The simulation was carried out for 10 ns, of which the first 100 ps were not used for analysis to account for equilibration.</p><p>From the MD trajectory, frames were selected for further analysis using a protocol described in the results section below. The protocol identifies structures in which the CPD, the HOCl and one water molecule are all oriented favorably for reaction to form the keto-enol product. Reaction profiles were determined for these frames, with the energy calculated at a QM/MM level using CHARMM interfaced with Q-Chem.55 The B3LYP/6-31+G(d,p) level of theory was used for the QM portion of all QM/MM calculations. In most of the calculations, the CPD, the HOCl and one water, which was involved in proton transfer, were allowed to be flexible and treated at the QM level while the rest of the water molecules were treated at the MM level with their positions kept fixed. We also investigated the influence of allowing waters molecules within 2.6 Å of the QM atoms to be flexible. These waters were allowed to relax in response to the reaction process but were still treated at the MM level.</p><p>The algorithms developed for identifying transition states along the reaction path for small systems are often not well suited to large systems such as proteins or full, explicit solvation.56 Calculating and manipulating an exact or approximate Hessian, which is required in many such methods, may be exhorbitantly expensive. Furthermore, when identifying a transition state and MEP in complex systems, discontinuities are often encountered due to structural relaxation to a conformationally distinct reaction path, resulting in hysteresis effects and making location of a consistent and meaningful MEP difficult.56 A common approach to identifying an MEP in complex systems is the chain-of-replicas class of methods.57-59 In such methods, a chain of conformations, or replicas, representing stages of reaction between reactants and products is created. Simultaneous energy minimization of all replicas with a constraint governing the distance between replicas yields a reaction path, where the distance between pairs of conformations is defined using a root mean square deviation (RMSD) metric.59 Our procedure comprised two steps beginning with the replica path (RP) method59,60 implemented in CHARMM interfaced with Q-Chem.61 The RP method enforces approximately equal spacing of replicas, resulting in MEPs that are accurate but do not locate the transition state with precision. To obtain the transition state, the reaction path was refined using the climbing image (CI) variant62 of the nudged elastic band (NEB) method63 (CI-NEB). The CI option causes a path point near the transition state to climb in energy toward the transition state, while all other path points are energy minimized in directions orthogonal to the MEP as in normal NEB. Details of our usage of the RP and CI-NEB calculations are given in the Supp. Info. This approach was applied to determine the MEP for the selected frames chosen from the MD trajectory.</p><!><p>Table 1 gives reaction barriers (ΔE‡) and reaction energies (ΔE) for the keto-enol tautomerization of CPD in the gas phase and in small clusters of up to nine water molecules. Both G3(MP2)//B3LYP/6-31+G(d,p) and B3LYP/6-31+G(d,p) values are given; they generally agree to within ±2 kcal/mol (average unsigned difference is 1.4 kcal/mol). The barrier for intra-molecular proton transfer in the gas phase is quite high, 62.9 kcal/mol. Inclusion of one reactive water molecule decreases the barrier to 27.3 kcal/mol. Reactant, transition state and product structures are shown for both reactions in Fig. 1. Single-point calculations using the PCM method to model solvation implicitly alter the barrier heights only slightly, by +0.3 and −1.5 kcal/mol without and with a reactive water, respectively. The keto-enol product is slightly more stable than the diketo reactant in both the uncatalyzed and the single-water-catalyzed reactions; PCM-modeled solvation increases the relative degree of product stabilization.</p><p>The effect of additional hydration was evaluated systematically by adding 1-2 waters to the CPD/HOCl complex at a time. Because systems with multiple waters possess a large number of degrees of freedom and great flexibility, interpretable results require careful positioning of waters. Incremental waters were placed in positions highly favorable for forming hydrogen bonds, until the hydrogen bonding positions around the CDP oxygens, the HOCl and waters participating directly in proton transfer processes were saturated. The new waters were added to a previously optimized transition-state structure, and the system was re-optimized to obtain the transition state for the more fully hydrated system. Reactants and products were obtained by performing intrinsic reaction coordinate (IRC) calculations64 in the forward and reverse directions, beginning from the transition-state geometry, followed by full energy minimization. This ensures that reactants, transition state and products pertain to the same MEP. The arrows in Figs. 1-3 indicate the paths of proton transfer as revealed by the IRC calculations. Trends in barrier heights and reaction energies exhibited a systematic relation to level of hydration, as described below. An alternative procedure of simultaneously adding waters to all stationary points followed by optimization yielded unsystematic results, probably indicative of a rough potential energy surface (PES) (details in Supp. Info.).</p><p>The reaction barrier is critically dependent on the positioning of incremental waters. Placing two waters so as to hydrate the HOCl (three waters total) and thereby increase its acidity had little effect on the reaction mechanism or barrier. The process remained concerted, with proton transfer from the alpha carbon to the reactive water and from the reactive water to the keto oxygen occurring simultaneously, as demonstrated by the transition-state Cα-H and O-H distances along the proton transfer path in Fig. 2. When two additional waters (five total) were placed so as to hydrate the reactive water, the barrier decreased dramatically, from 28.9 to 20.1 kcal/mol, an average of 4.4 kcal/mol for each hydrating water. The barrier reduction can be attributed to charge stabilization of the hydronium-like moiety formed in the transition state (Fig. 2); note that proton transfer from the alpha carbon to the reactive water is nearly complete, while proton transfer from the reactive water to the keto oxygen has hardly begun (transition-state Oketo-to-H distance of 2.0 Å vs 1.4 Å for the three-water model). As a result, abstraction of the alpha proton and protonation of the keto oxygen occur consecutively, rather than simultaneously as at the lower levels of hydration, via structures with reduced strain (Figs. 1 and 2). The importance of the water solvating the hydronium-like moiety is also apparent in the normal modes corresponding to the imaginary frequencies (Supp. Info.).</p><p>With the addition of a sixth water, initially placed so as to hydrate the reactive water, the mechanism changes more substantially. After energy minimization, two waters are positioned ready to participate directly in the reaction. When the reaction profile is calculated, a chain of three proton transfers occurs as the enol form is created (Fig. 3). The reaction barrier is reduced from 20.1 to 7.9 kcal/mol. This is the only case in which the keto oxygen interacting with HOCl was the one transformed to a hydroxyl during tautomerization; this is due to an additional hydrating water that preferentially stabilizes proton transfer in this direction (Fig. 3).</p><p>Addition of up to three more waters placed so as to hydrate waters participating directly in the reaction does not alter the mechanism substantially; two waters continue to be involved in the proton-transfer chain. However, the reaction barrier increases slightly, to 9.4 kcal/mol when two more waters are added (eight total), and more substantially, to 13.4 kcal/mol, when a third water is added (nine total). We have not identified an obvious rationalization for the 1.5 kcal/mol increase (7.9 to 9.4 kcal/mol), but the increase of 4.0 kcal/mol (9.4 to 13.4 kcal/mol) can be traced to the placement of the water added to the 8-water system to form the 9-water system. This water is a hydrogen bond donor to the water protonated by the alpha hydrogen during the first step of reaction. The hydrogen bond is weakened as the transition state is formed because the water accepting the alpha hydrogen becomes hydronium-like and thus positively charged; the hydrogen bond H…O distance increases from 1.76 to 1.99 Å as reactants transform to transition state. It is possible that lower-barrier reaction paths for 8- and 9-water models could be identified with effort. Small variations in the barrier height as the number of waters is increased further can be expected due to the complexity of the PES.</p><p>The reaction energy, ΔE, tends to decrease as the level of hydration is increased, with only one exception (three waters). This is probably because the keto-enol species is more polarizable than the diketo species, as chemical intuition suggests on the basis of the charge-separated resonance structure shown in Scheme 1(a) and as calculations confirm [isotropic polarizabilities calculated at the B3LYP/6-31+G(d,p) level are 9.65 (keto-enol) and 9.00 (diketo) in units of 10−24 cm−3]. The solvated keto-enol systems generally have tighter hydrogen bond networks, characterized by shorter hydrogen bonds, than the solvated diketo systems (see example 8-water structures in Supp. Info.). A count of the hydrogen bonds in the diketo and keto-enol structures for 6, 8, and 9 waters underscores the relationship between ΔE and the hydrogen bond network. Using distance and angle criteria of rOH ≤ 2.2 Å and θOHO ≥ 150°, the keto-enol products have zero, one, and two more hydrogen bonds than the diketo reactants for the 6-, 8-, and 9- water systems, respectively. This may account for the reduction of ΔE by 2.9 kcal/mol and 3.2 kcal/mol in going from 6 to 8 and 8 to 9 waters, respectively.</p><p>We establish above that using single-point PCM energies to evaluate the influence of implicit bulk water on the explicit-water models, equivalent to embedding the explicit water models in a dielectric continuum, gives reaction energies and barriers only slightly different from those calculated without PCM (same numbers of explicit waters, same geometries). In principle, the implicit representation of bulk water may affect the optimized geometries as well as the energies of given geometries. Geometry optimization with PCM can be tricky, but for reference we tested its influence on the reaction profiles for the 0-, 1-, and 3-water systems. The reaction energies and barriers obtained are given in parentheses in Table 1. Geometry optimization with PCM changes the values by only 0.5 kcal/mol on average, less than the difference between the in vacuo and single-point PCM values. The structures were only slightly altered by geometry optimization with PCM, the main difference being increased charge localization for the hydronium-like moiety in the transition states for 1 and 3 waters (shorter O-H distances; see Fig. S4 in Supp. Info.).</p><p>We investigated the influence of the acidity of the hydrogen bond donor to the CPD keto oxygen not directly involved in the keto-enol tautomerization (hereafter nonreacting keto oxygen). This was done by removing HOCl from the model altogether or by replacing HOCl with H2O. Table 2 shows a comparison of reaction energies and barriers (a) without any such hydrogen bond donor, (b) with a water (pKa=15.7), and (c) with HOCl (pKa=7.5). One reactive water was also present in all three cases, but there were no additional solvating waters. Compared to the system having no hydrogen bond donor to the nonreacting keto oxygen, hydrogen bond donation by a water or an HOCl reduces the barrier by 3.3 or 4.5 kcal/mol, respectively. Thus interaction of the nonreacting keto oxygen with even a weak acid reduces the reaction barrier substantially, consistent with general organic chemistry principles. In none of our calculations was a proton transferred from the HOCl to a keto oxygen of CPD at any point during reaction, as is normally depicted in mechanisms of acid-catalyzed keto-enol tautomerization. This is presumably because HOCl is a weak acid.</p><p>The work on HOCl-catalyzed reaction of CPD in water clusters is summarized as follows. A single reactive water (a one-water circuit) participating in the reaction lowers the barrier from 62.9 to 27.3 kcal/mol. Hydration of the reactive water can lower the barrier substantially, by as much as 9 kcal/mol, as the five-water system demonstrates. Further hydration makes longer proton-transfer circuits accessible. We have identified three two-catalytic-water circuits with barriers ranging from 7.9 to 13.4 kcal/mol. Higher levels of hydration may lead to yet longer proton-transfer circuits, such as were constructed for malonaldehyde26 and observed in a metadynamics simulation of acetone.25 The contribution of the weak acid HOCl to barrier reduction is approximately 4.5 kcal/mol.</p><p>Explicit water molecules positioned in the manner described stabilize transition states relative to reactants, via a combination of charge stabilization of the hydronium-like moiety in the transition state and structural optimization of water-mediated proton transfer networks. Normally reaction would be predicted to occur primarily via the path with the lowest barrier. However, caution should be exercised in extrapolating from the water cluster results to conclusions about the rate of reaction in solution at ambient temperature. The entropic cost of forming such extended, highly optimized hydrogen bond networks at 300 K may make them minor contributors to the reaction rate in solution at 300 K. Indeed this is demonstrated by our results below for fully solvated systems at 300 K.</p><!><p>This portion of the work was undertaken to evaluate whether the fully optimized structures described above are consistent with solvent structures encountered in full solution at room temperature. Molecular dynamics simulation of one molecule each of CPD and HOCl with 324 explicit water molecules was performed at 300 K and analyzed to obtain snapshots of reaction precursor structures. The MEP was determined for each precursor structure using the RP and CI-NEB methods in the context of a QM/MM model.</p><!><p>Reactant molecules must be oriented favorably with respect to each other for reaction to occur. Such positioning, a 'rare event' occurrence, is unlikely to occur frequently on the nanosecond timescale. Therefore frames representing such 'rare event' occurrences were chosen from the 10-ns MD simulation. We focused on rare events favorable for single-water proton transfer chain mechanisms, in order to compare with water cluster results for the same type of mechanism. The analysis could in principle be extended to multiple-water proton transfer chain mechanisms, but the complexity of the PES renders the calculations more difficult for the reasons described in the Computational Methods section. The frames were selected such that the HOCl was near enough to CPD to form a hydrogen bond with one of the keto oxygens, and a water molecule was positioned ready to accept a proton from the alpha carbon of CPD and to donate a proton to the other keto oxygen. This was accomplished using the following criteria: The distance between one of the keto oxygens of CPD and the hydrogen of HOCl is less than or equal to 2.4 Å (1). There exists a water molecule whose oxygen is within 3.2 Å of the alpha hydrogen of CPD (2) and one of whose hydrogens is within 3.2 Å of the keto oxygen not used in the first criterion (3). We initially used a smaller distance for criterion (3), consistent with a weak hydrogen bond, but that resulted in no frames being selected.</p><p>Statistics regarding the trajectory frames satisfying these criteria are shown in Fig. 4, and a representative selected structure is shown in Fig. 5. It is visually apparent from Fig. 4 that there were ~20 unique encounters during the 10-ns MD simulation, and the encounters had durations of 0.1-1.0 ps and were separated by more than 10 ps. Unique encounter refers to a brief period during which the CPD and HOCl remain close to each other, possibly due to trapping by a solvent cage, such that the three criteria defining a rare event occurrence remain continuously or nearly continuously satisfied. For reference, these encounter durations are on the order of measured65,66 and simulated67-70 first-hydration-shell residence times (0.3-2.5 ps), reorientation times (~2 ps) and water-water hydrogen bonds (0.2-2 ps) for waters in bulk solution or near protein surfaces at temperatures near 300 K. For each unique encounter, the best frame was selected for further consideration. The best frame was defined as the one with the smallest value for criterion (1) above; in the case of a tie, the smallest value for criterion (2) was determining; and if a time still remained, the smallest value for criterion (3) was decisive.</p><!><p>Reaction barriers for the selected frames were determined using a QM/MM model with the RP and CI-NEB chain-of-replicas methods as described earlier. The CPD, the HOCl and a single water molecule participating in the reaction explicitly (i.e. a reactive water) were treated at the QM level, and the rest of the water molecules were treated at the MM level. All the MM waters were kept fixed while the QM atoms were flexible. Thus the barriers calculated are fully optimized with respect to the flexible (QM) degrees of freedom, in a fixed solvent environment (MM) representative of a liquid solution at 300 K. Physically, this represents a reaction that takes place much faster than solvent response. Although the fixed-MM-water model is a simplification of the actual reaction process, it provides a well-defined limiting behavior for analysis. Each frame represents a different solvent environment consistent with a temperature of 300 K. In principle, Boltzmann weighted averaging of barriers across a large number of frames will provide an effective activation barrier for the limiting physical process considered. However, adequate sampling is not practical given the expense of the QM/MM calculation of the MEPs. We have analyzed ten frames to qualitatively assess the variability in barrier height. A Boltzmann-weighted average of these ten frames represents out best estimate of the effective activation barrier. The Boltzmann-weighted average barrier is dominated by the smallest barrier. Such an approach is commonly used in QM/MM calculations of reactions in protein environments.56 Because the MM waters are fixed, their internal energies and the interactions between them are identical for the reactants, transition state and products for a particular frame and cancel for the barrier height. This is convenient, as otherwise thermal fluctuations in these quantities could mask the barrier height dependencies of interest. Details regarding the RP and CI-NEB reaction path calculations, including figures establishing convergence, are given in the Supp. Info.</p><!><p>Table 3 lists the barrier height, ΔE‡,full, for each of the selected frames. The table also gives the component to the barrier due to the QM region only, ΔE‡,QM; a correction, ΔΔE‡,shell, due to the innermost shell of MM waters; and a correction, ΔΔE‡,bulk, due to the rest of the MM waters. The innermost shell was defined as those MM waters containing at least one atom within 2.6 Å of a QM atom (6-11 waters depending on the frame). This definition was chosen with reference to pair distribution functions for water surrounding keto oxygens and for pure water. At ambient temperature, the trough following the first peak in the O-H pair distribution function is in the range 2.6-2.7 Å for water surrounding keto oxygens71,72 and 2.4 Å for pure water modeled with a modified TIP3P force field.73</p><p>The barrier component ΔE‡,QM was obtained from QM single-point calculations on the atoms of the QM region by itself. The geometries used were the stationary points from the QM/MM optimization of the full system, but the MM waters were deleted. The QM and QM/MM calculations used the same level of QM theory, B3LYP/6-31+G(d,p). The correction ΔΔE‡,shell was obtained as the difference between ΔE‡,QM and the barrier for a system consisting of the QM region plus the innermost shell of MM waters (ΔE‡,shell), based on the stationary point structures for the full system and deleting the appropriate atoms. The correction ΔΔE‡,bulk is the difference between ΔE‡,full and ΔE‡,shell (we use bulk to refer to all waters beyond the first hydration shell). For comparison, the barriers were recalculated with the innermost shell of MM waters unconstrained. The change in the barrier height due to this additional flexibility is given in Table 3 as ΔΔE‡,flexshell, and the resulting barrier itself is ΔE‡,full,flexshell.</p><p>A complementary analysis of solvent effects is provided by considering stepwise modification of the subsystem treated quantum mechanically (CPD + HOCl + reactive water molecule). In the first step, the subsystem is distorted from the optimized gas-phase geometry (i.e., optimized without solvating MM waters) to the optimized solution-phase geometry; the difference in energy is Egeomdist. In the second step, the solvating MM waters are included, and the electrostatic (EelecQM/MM) and Van der Waals (EVDWQM/MM) components to the interaction between the QM and MM regions are determined. These three energy components averaged over the ten frames are given for reactants and the transition state in Table 4 (values for individual frames are given in the Supp. Info.). Table 5 presents the contribution to the reaction barrier due to each of the components. These barrier contributions are defined as the difference between the corresponding components of the transition-state and reactant energies, e.g.
ΔE‡,elecQM∕MM=E‡,elecQM∕MM−−ER,elecQM∕MM
Other methods for analyzing contributions to reaction barriers in solution have been devised, including calculation of components such as electronic polarization of solvent; analysis of these factors is beyond the scope of this work.74-77</p><p>There is a substantial variability of 25-44 kcal/mol in the barrier height, ΔE‡,full. The smallest barriers (25 and 27 kcal/mol) are close to the water-cluster barrier with one reactive water and no solvating waters [27.6 and 26.1 kcal/mol in gas phase and in continuum dielectric, respectively, at the B3LYP/6-31+G(d,p) level (Table 1)], and significantly higher than the water cluster barrier with one reactive water and four solvating waters [18.6 and 17.5 kcal/mol in gas phase and in continuum dielectric, respectively (Table 1)]. Because Boltzmann averaging preferentially weights the smaller values, a barrier of 25-27 kcal/mol represents our best estimate, given the limited sampling, of the 'real' barrier that would be obtained by barrier-crossing beginning with a thermodynamic ensemble of reactant conformations.56 Thus a model with slow response of the nearest waters yields barriers close to the cluster simulations without solvating waters.</p><p>On average, the innermost shell of waters raises the barrier by 5 kcal/mol, while the more distant waters have very little effect. The influence of the innermost shell of waters can be traced to a greater degree of hydrogen bonding between solvent waters and the reactive water in reactant than in transition state structures. This is apparent from visual analysis of the frames, which reveals hydrogen bonding structures, together with comparison of the QM/MM electrostatic energy components in Tables 4 and 5. Hydrogen bonding is largely an electrostatic phenomenon, and the average QM/MM electrostatic energy (ΔE‡,elecQM/MM) favors reactants by 6 kcal/mol, which is nearly the same as the average contribution to the barrier height attributed to the innermost solvent shell (ΔE‡,shell, 5 kcal/mol). The QM/MM electrostatic contribution to the barrier is positive for all frames except frame 1, and it is responsible for some of the highest barriers for individual frames (see ΔΔE‡,shell and ΔE‡,elecQM/MM for frames 6 and 8).</p><p>Analysis of individual frames reveals that waters beyond the first hydration shell can influence barriers significantly even though their influence averages to near zero over all barriers. Bulk waters raise or lower the reaction barrier by up to 7 kcal/mol (ΔΔE‡,bulk) via electrostatic interaction with the reactive system (Table 3). This result implies that rapid fluctuations in solvent structure at least as far away as the second hydration shell can be important for reaction. Thus continuum dielectric solvent models such as the PCM model, used to solvate structures listed in Table 1, may correctly capture average effects and yet miss important short-time effects on the reaction barrier. We also note that the TIP3P water model underrepresents structure in the second hydration shell of bulk water; the oxygen-oxygen radial distribution function is too flat in this region.45,78 Since we find that solvent structure tends to increase the barrier, we speculate that a more structured second hydration shell may on average slightly increase barriers.</p><p>On average, the Van der Waals influences on reactants and the transition state are equal and therefore cancel in their effect on the barrier height (Table 4). Even in individual frames, Van der Waals interactions have only a small effect on the barrier height (Table 5). One should, though, be cautious not to overinterpret this result, given the limitations of the force field model via which the Van der Waals interactions are calculated. Geometric distortion contributions also have a negligible average influence on the barrier (Table 4), but they can be significant for individual frames and in fact are the largest component of the QM/MM contribution to the highest barrier (frame 5).</p><p>The predominant solvent effect revealed by this analysis, then, is to raise the barrier via preferential electrostatic stabilization of reactants that reflects hydrogen bonding interactions between innermost-shell solvent waters and the reactive water. Organization of bulk water dipoles can further raise or lower the barrier by up to 7 kcal/mol, although the average influence is zero. Geometric distortion influences can also be important for individual barriers but average to near zero. However, even though solvation of the reactive water tends to raise the barrier, the Boltzmann-weighted average barrier in solution will be dominated by the lowest barriers, which are very close to the gas-phase barrier with one reactive water present.</p><p>This analysis reveals that, for a model of thermally randomized but fixed solvating MM waters, the stationary points on the PES have less than the optimal hydrogen-bonding stabilization identified above in the water cluster analysis. As discussed earlier, this model represents the limit of infinitely slow solvent response to the course of reaction. An alternative limit representing infinitely fast local solvent response was investigated via calculation of ΔΔE‡,flexshell. First-hydration-shell reorganization reduced the average barrier by 3 kcal/mol and reduced but did not eliminate the variation in barrier heights. Thus flexible water molecules have a tendency to reduce the barriers. Depending on the time scale of water response, some level of hydrogen bond stabilization is likely to be observed in aqueous solution. However, regardless of the time scale of response, the level of stabilization is significantly smaller than in the fully optimized hydrogen bond networks obtained in the limited-solvation study.</p><!><p>The influence of solvation on the HOCl-catalyzed keto-enol tautomerization of CPD has been studied with two theoretical approaches: QM-based calculation of reaction paths using fully optimized structures with limited levels of solvation, and QM/MM-based calculation of reaction paths with full solvation. The limited solvation results are consistent with calculations for keto-enol tautomerization of related molecules. We determine that the barrier with inclusion of up to one solvation shell of explicit water molecules depends strongly on the number of waters directly involved in proton transfer (barrier reduction of 36 kcal/mol for the first such water and 12 kcal/mol for second), and to a lesser but significant extent on the number of waters forming hydrogen bonds with waters in the proton transfer chain (barrier reduction of ~4.4 kcal/mol for each of the first two). In accord with the understanding that acids catalyze keto-enol tautomerization, the interaction of HOCl with CPD reduced the barrier by 4.5 kcal/mol; however, proton transfer from HOCl to CPD during formation of an intermediate or transition state was not observed, presumably because HOCl is a weak acid. Barriers of 8-13 kcal/mol were obtained when HOCl and a full or nearly full solvation shell were present.</p><p>In the full solvation work, reaction barriers were determined for rare-event snapshots of reactants positioned favorably for reaction, chosen from a 300 K molecular dynamics simulation of CPD, HOCl and 324 explicit waters. The barriers were determined with a chain-of-replicas method employing a QM/MM energy representation. A substantial variation in barrier heights was observed, depending primarily on electrostatic interactions (hydrogen bonding) with first-solvation-shell waters and to a lesser extent on electrostatic interactions with bulk waters and other influences mediated by geometric distortion of the reactant CPD. For the rate-determining barriers, the extent of barrier reduction by waters involved in proton transfer is consistent with the limited-solvation results, but further barrier reduction due to hydrogen-bonding to waters involved in proton transfer is not seen. Our results show that it is difficult to observe highly flexible structures such as extensive hydrogen bonding networks that are optimal for reaction in solvated simulations at 300 K. Even transient formation of such optimal networks may be so infrequent that they do not contribute significantly to the observed reaction rate in solution. This implies that barriers calculated beginning with fully energy-minimized structures at limited levels of solvation, the basis of a common theoretical approach, may be deceptively small if they involve intricate hydrogen bonding networks and are applied to understand the reaction in solution at ambient temperatures. A quantitative assessment will be necessary to address this issue definitively but is beyond the scope of this paper. On the other hand, study of limited solvation conditions may be revealing of mechanisms available in special environments, such as the interior of enzymes, where the local structure may favor precise positioning of waters otherwise entropically disfavored. A stable water network near a protected active site or access channel in an enzyme could facilitate the diketo—keto-enol tautomerization of CPD.</p>