The alkaline hydrolysis of sulfonate esters: challenges in interpreting experimental and theoretical data.

Sulfonate ester hydrolysis has been the subject of recent debate, with experimental evidence interpreted in terms of both stepwise and concerted mechanisms. In particular, a recent study of the alkaline hydrolysis of a series of benzene arylsulfonates (Babtie et al., Org. Biomol. Chem. 10, 2012, 8095) presented a nonlinear Brønsted plot, which was explained in terms of a change from a stepwise mechanism involving a pentavalent intermediate for poorer leaving groups to a fully concerted mechanism for good leaving groups and supported by a theoretical study. In the present work, we have performed a detailed computational study of the hydrolysis of these compounds and find no computational evidence for a thermodynamically stable intermediate for any of these compounds. Additionally, we have extended the experimental data to include pyridine-3-yl benzene sulfonate and its N-oxide and N-methylpyridinium derivatives. Inclusion of these compounds converts the Brønsted plot to a moderately scattered but linear correlation and gives a very good Hammett correlation. These data suggest a concerted pathway for this reaction that proceeds via an early transition state with little bond cleavage to the leaving group, highlighting the care that needs to be taken with the interpretation of experimental and especially theoretical data.

Introduction

Compared to the large body of work on phosphoryl transfer reactions,[1,2] sulfuryl transfer has received rather less attention. However, recent years have seen a revival of interest in sulfate[3,4] and other related group-transfer reactions.[5,6] Sulfate hydrolysis has directly important biological roles in, for instance, cellular signaling[7] and detoxification,[8] and sulfonate esters can be used by bacteria as sulfur sources in environments with low sulfate concentrations.[9,10] Furthermore, there has been an increased awareness of the prevalence of catalytic promiscuity in enzymes catalyzing phosphoryl and sulfuryl transfer reactions,[11,12] a phenomenon that has been suggested to play an important role both in natural protein evolution[13] and artificial enzyme design.[14]

Recently, a phosphonate monoester hydrolase (PMH) from was found to be one of the most promiscuous hydrolases characterized to date,[15] catalyzing at least five classes of hydrolytic reaction in addition to its native phosphonate monoesterase activity. One of these activities is the hydrolysis of xenobiotic sulfonate monoesters, and this PMH is the only known enzyme capable of catalyzing the hydrolysis of these compounds through direct S-OR cleavage (as opposed to C–O cleavage).[15] Elucidating the precise molecular basis for this promiscuity will enhance our understanding of the evolution of protein function, and an important first step in this direction is obtaining detailed insight into the chemical similarities and differences between the different substrates at the reactant level.

Although less studied than their phosphate counterparts, there has been substantial research effort invested into deciphering the mechanism and nature of the transition states for the hydrolyses of sulfate[3,4,7,16−19] and neutral sulfonate[6,20−29] monoesters. The hydrolysis of a sulfonate monoester can, in principle, proceed through multiple mechanisms, and arguments have been made in favor of both stepwise addition–elimination mechanisms[6,26−28] (involving a pentacoordinate intermediate), as well as concerted mechanisms[24,29] (involving a single transition state) (Figure 1).

Figure 1

Potential stepwise (A) and concerted (B) pathways for the alkaline hydrolysis of aryl benzenesulfonates.

Although not inconceivable, a dissociative DN + AN mechanism has been generally ruled out.[23] In particular, Williams and co-workers[24] examined the reaction of a range of oxyanions with basicities spanning a range of over 8 pKb units with 4-nitrophenyl 4-nitrobenzenesulfonate. They also examined the reaction of phenoxide with substituted phenyl benzenesulfonic acids. In both cases, they found good linear correlations spanning the pKa of the incoming nucleophile, causing them to argue that a single transition state provides the best explanation for the process. Experimental studies on the alkaline hydrolysis of diaryl sulfate diesters have supported a similar mechanistic pathway.[30]

Very recent work[6] has suggested a more complex scenario, presenting a Brønsted plot for the alkaline hydrolysis of aryl benzenesulfonates that breaks at pKa,lg ∼ 8.5, with βlg values of −0.27 and −0.97 for leaving groups with pKas lower and higher than 8.5, respectively (Figure 2A; dashed lines). A curve based on a stepwise reaction involving a change in the rate-determining step[31] also fits the data very convincingly (Figure 2A, solid line; see also discussion in ref (6)). However, this interpretation would suggest that at the breakpoint, the loss of either hydroxide or phenolate occur at equal rates, which is implausible given the very different nucleofugalities of these two leaving groups. For a stepwise process, the rate-limiting step is expected to be the attack of the hydroxide nucleophile for all the compounds studied. Thus, it was proposed that the break in the Brønsted plot corresponds to a switch from a two-step mechanism involving a pentavalent intermediate for poorer leaving groups (Figure 1, pathway A), to a concerted mechanism with good leaving groups (Figure 1, pathway B). It was suggested that for substrates with a sufficiently good leaving group (i.e., with a pKa lower than the breakpoint in the plot), the attack of hydroxide leads to an intermediate that is too unstable to have a significant lifetime, and so the reaction has to become concerted.[6] Although a change in mechanism would usually show curvature that is concave upward in a Brønsted plot, a similar concave-downward Brønsted plot is observed for the reactions of phosphate monoester monoanions and has been rationalized in related terms. However, in this analysis, the rate-limiting step is the breakdown of a tautomer, which is assumed to form in a fast pre-equilibrium step. When this tautomer is so reactive that it is predicted to have a lifetime shorter than a bond vibration, the reaction has to become concerted, with bond cleavage accompanying the proton transfer. The difference with the sulfonate ester hydrolysis reactions discussed here is that the initial attack is always assumed to be the rate-limiting step, even when an intermediate is proposed.[32]

Figure 2

(A) Brønsted correlation for the alkaline hydrolysis of aryl benzenesulfonates. Linear least squares fits to the low (blue) and high (green) pKa data are shown as dotted lines, as reported in ref (6). The solid line is the least squares fit for a stepwise mechanism involving a change in rate-limiting step. (B) Hammett correlation for the same data. The solid is the linear least squares fit, giving ρ = 1.83 ± 0.10 and R2 = 0.981.

The status of the potential intermediate was investigated using hybrid effective fragment potential (EFP)[33−35] and conductor-like polarizable continuum model[36] (C-PCM) calculations (QM(HF)/EFP/PCM approach) using a number of explicit water molecules surrounding the sulfonate. These showed that a pentavalent intermediate was only stable with the poorer leaving groups, but not for the better leaving groups, supporting the proposed mechanistic change. One concern with these calculations is that it was necessary to introduce at least 8 discrete water molecules into the calculations to obtain the intermediate shown in Figure 2 of ref (6). We have discussed the problems such a mixed solvation model introduces in terms of quantitative accuracy of the calculations elsewhere[37] and guide interested readers to this work for further information. Additionally, these calculations were performed at a relatively low level of theory (HF), raising the possibility that this intermediate is actually a simulation artifact (note also that only a single stationary point with no corresponding transition state(s) is presented in ref (6), and that all the water molecules are placed where they selectively stabilize the nucleophile). Thus it is possible that the calculations do not provide unambiguous support for the proposed mechanism.

Overall, there is sufficient doubt about the interpretation and analysis of these compounds to test whether an intermediate is theoretically viable,[6,24,26−29] so we have revisited the LFER presented in ref (6), generating detailed free energy surfaces for the hydrolysis of each compound, following our previous work on the hydrolysis of phosphate and sulfate esters.[3,38−41] We also discuss the broader challenges of obtaining reliable quantitative results for reactions involving both anionic nucleophiles and polarizable atoms such as sulfur using conventional density functional (DFT) approaches.

Finally, we note that while the experimental data appear to generate a smooth, nonlinear Brønsted plot (Figure 2A), a Hammett plot of the same data yields a simple linear correlation (Figure 2B), raising the question of whether correlating the reactivity with the leaving group pKa is appropriate. This has been noted in previous work concerning sulfonate esters: Buncel et al. concluded that σ (rather than σ–, analogous to pKa) was a more appropriate parameter to use for the ethoxide-promoted cleavage of sulfonate esters in anhydrous ethanol,[26] and Um et al. reported a similar result.[28] These analyses led to the suggestion that the reaction proceeds through a stepwise reaction, where the rate-limiting step is the formation of a pentavalent intermediate. However, in subsequent work, Um et al. concluded that a significant contribution from both σ and σ– (derived from a Yukawa–Tsuno analysis) provided the best fit under largely aqueous conditions.[29] These data include substituents in the 2 position, which might be expected to generate scatter due to steric interactions, thus confusing an interpretation based solely on electronic variation. Most recently, Um and Buncel suggest that this is the best analysis to use for the reaction of ethoxide with aryl benzenesulfonate esters, and that the best description of the reaction is concerted, with significant bond formation to the nucleophile but with relatively small bond cleavage to the leaving group in the transition state.[42] Williams et al. also considered these factors in the formation of a cyclic sulfonate with aryloxy leaving groups, and showed that σ– (effectively a Brønsted plot) rather than σ provided the most appropriate analysis.[23] Thus, there is ambiguity about the best way to represent the kinetic data to reveal a reliable correlation.

Hence, to complement these theoretical studies, we also studied the hydroxide-promoted hydrolysis of pyridin-3-yl benzenesulfonate (1), along with its N-oxide (2) and N-methylpyridinium (3) derivatives (Figure 3). The lowest pKa leaving groups in Figure 2A all have substituents that provide strong resonance stabilization of the oxyanion product; if this effect is not significant in the transition state, this may account for why they are less reactive than predicted by the line defined by the poorer (nonresonance stabilized) leaving groups. Therefore, these additional compounds were selected as they also have very good leaving groups, but due to inductive and field effects rather than direct conjugation of the anion with the π-system of a substituent. If the Brønsted plot really does reveal limiting behavior dependent on leaving group ability, these compounds will fall close to previously reported data; if they do not, and the data can be described by a smooth Hammett correlation, then the experimental evidence for discontinuous behavior will no longer be apparent.

Figure 3

Structures of (1) pyridine-3-yl benzene sulfonate, as well as (2) its N-oxide and (3) N-methylpyridinium derivatives.

On the basis of the results of our combined theoretical and experimental studies presented in this report, we conclude that a concerted reaction mechanism involving a single transition state provides a better interpretation of the combined data presented here and in ref (6). Importantly, we believe that this system provides an excellent demonstration of the pitfalls in the interpretation of experimental and especially computational data.

Computational Methods

Figure 4 shows the substituted aryl benzenesulfonates studied in this work and originally presented in ref (6). This set of compounds was then supplemented with compounds 1–3 (Figure 3), as discussed in the Introduction. In our previous computational studies of the LFER for related compounds,[39,41] we have generated 2-D potential energy surfaces for only one compound in the series, and we used this as a starting point to obtain transition states with perturbed leaving groups for the remainder of the series. However, in the present work, as there is potential for two different mechanisms depending on leaving group, we have generated individual potential energy surfaces for all compounds of interest. These were obtained in the space defined by S–O distances to the departing leaving group (x-axis) and the incoming nucleophile (y-axis), spanning a range from 1.5 to 2.8 Å on the x-axis, and 1.5 to 3.3 Å on the y-axis.

Figure 4

Aryl benzenesulfonates examined in this work.

We originally started with surfaces corresponding to three compounds at either extremes of the series and at the break point (i.e., the 3-F-4-NO2, 4-Cl and 3,4-dimethyl compounds), for which we mapped the surface using a finer grid of 0.1 Å increments in the bond distances on each axis. For the remaining compounds, we used a slightly coarser grid of 0.15 Å increments. At each point on the plot, the two degrees of freedom corresponding to the reaction coordinate were kept frozen, while all other degrees of freedom were allowed to freely optimize, and the surface was obtained by carefully pushing the reaction coordinate in all relevant directions until the complete surface was obtained. All points on the surface were obtained using Truhlar’s M06-2X functional,[43] which is a dispersion corrected hybrid metaexchange-correlation functional that has been rigorously parametrized for organic compounds. Solvation was simulated implicitly, using Cramer and Truhlar’s SMD solvation model.[44] In order to save computational cost, initial geometry optimizations were performed using the smaller 6-31+G* basis set, followed by a single point energetic correction to the obtained geometries using the larger 6-311+G** basis set. In each case, the final surface was used to locate an approximate transition state geometry, which was then optimized to a saddle point using an unconstrained transition state optimization. The resulting structures were characterized by frequency calculations, and the minimum energy path connecting reactants to products through this transition state was evaluated by calculating the intrinsic reaction coordinate[45,46] (IRC = ξ) in both the forward and reverse directions in order to determine whether an intermediate could be found.

Figure S1A (Supporting Information) shows an overlay of the calculated IRC from each transition state for each compound (normalized to reactants), using the M06-2X functional, and Figure S1B (Supporting Information) shows the change in S–Onuc and S–Olg distances along the reaction coordinate for a representative compound (4h). From these figures, it can be seen that, apart from the energetics of the product state, the energy profiles for the hydrolysis are very similar, with early reactant-like transition states. The final structure from following the IRC in the reactant direction was then subjected to an unconstrained geometry optimization in order to obtain the geometry of the reactant complex, which was then used to evaluate the activation barrier for the reaction. For comparison, the transition states obtained using M06-2X were then reoptimized using a number of other DFT functionals: the popular B3LYP[47−49] functional, the dispersion corrected ω-B97X-D,[50] and finally CAM-B3LYP[51] (which was designed to correct delocalization error). New IRC calculations were run from these TS in order to verify the identity of the reoptimized transition states. We have also explored the effect of including an increasing number of discrete water molecules (2, 8, and 10) in the calculations on the nature of the transition state using the Hartree–Fock (HF) method, for direct comparison with previous work.[6] Finally, free energies in the absence of explicit water molecules were obtained by evaluating the zero point energies and entropies from calculated vibrational frequencies using the 6-311+G** basis set, SMD solvation model and relevant DFT functional. Note that we did not include the free energy corrections when examining systems with extra discrete water molecules, as the frequency calculations become increasingly unreliable as the total number of degrees of freedom is increased. All calculations presented in this work have been performed using the Gaussian 09 quantum chemical suite of programs.[52]

Results and Discussion

Free Energy Surfaces and Initial Characterization of Transition States

To probe the potential mechanism(s) for the hydrolysis of the compounds listed in Figure 4, our starting point involved the generation of 2-D potential energy surfaces for hydroxide attack on each of these compounds, as outlined in the Computational Methods. Here, the first question involves identifying the relevant conformation of the two aromatic rings with respect to each other for our optimization, as they could take either of two conformations as illustrated in Figure 5. In the first of these, the two aromatic rings have a weak π-stacking interaction, whereas in the second, the rings avoid each other. The latter conformation is the one used in the calculations of ref (6); however, we were only able to obtain optimized transition states in implicit solvent using the conformation in which the two rings stack together. We have performed an energy scan of the O–S–Olg–C dihedral angle in the parent compound (Figure 5) to characterize the effect of adjusting the position of the two aryl rings with respect to each other.

Figure 5

Scan of the O–S–Olg–C dihedral angle in the ground state of the parent compound.

From this figure, it can be seen that there are two minima along the reaction coordinate, corresponding to each of the conformations outlined above, and that the conformation with the weak π-stacking interaction appears to be the favored conformation in the ground state. Therefore, this is the conformation we have used for all our subsequent calculations in implicit solvent. It is important to bear in mind that this could be an artifact of the implicit solvent model. However, as discussed in Section 3.2, we would like to emphasize that performing the calculations using the alternate conformation (corresponding to that used in ref (6)) in the presence of a number of explicit water molecules does not change the qualitative nature of our results, and either conformation would in principle form a viable starting point to similar effect.

Potential energy surfaces (PES) for hydroxide attack on 3-fluoro-4-nitrophenyl (pKa 6.11), 3-nitrophenyl (pKa 8.35), 4-chlorophenyl (pKa 8.61), and 3,4-dimethyl benzene sulfonate (pKa 10.36) are shown in Figure 6A–D, with the PES for the remaining compounds being presented in Figure S2A–D (Supporting Information). The corresponding transition state geometries, obtained using optimization with M06-2X/6-311+G**, are shown in Figure 7A–D and Figure S3A–D (Supporting Information), respectively. The compounds highlighted in Figures 6 and 7 correspond to the two pKa extremes of the LFER presented in Figure 1 in ref (6), as well as the two compounds that lie on the possible break-point of this plot. As can be seen from these figures, in all cases the potential energy surface indicates the presence of a single, concerted transition state, with no intermediate. The approximate transition state from each surface was then subjected to a full geometry optimization and characterized as described in the Computational Methods, and the resulting energetics (including a breakdown of different contributions to the calculated activation barrier) as well as S–O distances to the incoming nucleophile and departing leaving group at the transition state are highlighted in Tables 1 and 2, respectively. The fact that these are concerted transition states is further verified by IRC calculations in both forward and reverse directions, which all lead to either reactant or product complexes. Note that all energetics presented in Table 1 are relative to the corresponding reactant complex obtained by following the IRC as far as possible in the reactant direction, and concluding with a final geometry optimization to obtain the stationary point. In terms of geometric parameters (Table 2), while there appears to be negligible change in the S–Olg distance moving across the series, there is a gradual increase in S–Onuc distance, from 2.35 Å for the compound with the poorest leaving group (3,4-dimethyl benzene sulfonate) to 2.47 Å for the compound with the best leaving group (3-fluoro-4-nitro benzene sulfonate), suggesting that the better the leaving group, the earlier there is commitment to the reaction, resulting in a gradually more expansive transition state. The shifting transition states observed across these series of compounds follow a trend we previously observed for both methyl arylphosphate diesters[41] and fluorophosphates[39] (although these transition states are very slightly more expansive than the ones obtained here).

Figure 6

Free-energy surfaces for the alkaline hydrolysis of aryl benzenesulfonate. Shown here are 3-fluoro-4-nitrophenyl, 3-nitrophenyl, 4-chlorophenyl, and 3,4-dimethylphenyl benzenesulfonates. The approximate positions of the relevant transition states are indicated by ‡, and the actual optimized structures are shown in Figure 7

Figure 7

Geometries of optimized transition states for alkaline hydrolysis of (A) 3-fluoro-4-nitrophenyl, (B) 3-nitrophenyl, (C) 4-chlorophenyl, and (D) 3,4-dimehtylphenyl benzenesulfonates. These structures were obtained by optimization of the approximate geometries highlighted on the surfaces shown in Figure 6. The corresponding transition states for the remaining compounds shown in Figure 3 are presented in Figure S3 (Supporting Information).

Table 1

Energy Decomposition for ΔG‡calc at the M06-2X/6311+G** Level of Theorya

substituted phenol	pKa	ΔE‡gas	ΔΔGsolv	ΔZPEb	–TΔS‡b	ΔG‡calc	ΔG‡expc
3F-4-NO2	6.11	–21.9	33.4	–0.1	3.2	14.5	17.3
4-NO2	7.14	–22.9	34.8	–0.1	3.1	15.0	18.0
4-CN	7.95	–22.2	34.6	–0.2	3.0	15.2	18.1
3-NO2	8.35	–21.8	34.1	0.1	3.2	15.5	18.2
3-CN	8.61	–21.8	34.4	0.1	3.3	15.9	18.3
4-Cl	9.38	–20.9	34.0	0.0	3.8	16.8	19.3
H	9.95	–19.5	33.1	0.2	3.6	17.3	20.2
3,4-dimethyl	10.36	–19.9	33.9	–0.5	2.2	15.8	20.9

All energies are in kcal/mol, relative to the reactant complex.

Zero point energies and entropies were obtained from frequency calculations at 323.15 K.

All experimental values have been corrected for the entropic cost of bringing the reacting fragments into the reacting cage (K = 0.017 M–1), following ref (53). This is important for consistency, as our reference state in the calculations is a reactant complex (obtained by following the minimum energy path from the relevant transition state) and not the individual fragments at infinite separation from each other.

Table 2

Comparison of S–O Distances to the Nucleophile (S–Onuc) and Leaving Group (S–Olg) at the Reactant and Transition States for the Alkaline Hydrolysis of the Benzene Sulfonates Shown in Figure 4a

		RSb		TSb		ΔΔRS → TSb
substrate	pKa	S–Onuc	S–Olg	S–Onuc	S–Olg	S–Onuc	S–Olg
3F-4NO2	6.11	3.62	1.65	2.47	1.72	–1.15	0.07
4-NO2	7.14	3.63	1.64	2.41	1.72	–1.21	0.08
4-CN	7.95	3.59	1.64	2.41	1.72	–1.19	0.08
3-NO2	8.35	3.63	1.64	2.42	1.71	–1.21	0.07
3-CN	8.61	3.61	1.64	2.40	1.71	–1.21	0.08
4-Cl	9.38	3.67	1.63	2.38	1.71	–1.29	0.08
H	9.95	3.63	1.63	2.38	1.71	–1.25	0.08
3,4-CH3	10.36	3.64	1.63	2.35	1.71	–1.29	0.08

All geometry optimizations were performed using the M06-2X functional, and all distances are shown in Å.

RS and TS denote reactant and transition states, respectively. ΔΔRS → TS denotes the change in distance upon moving from RS to TS respectively.

As an aside, an important issue to take into account in our calculations is the challenges of obtaining quantiative accuracy in calculations involving hydroxide as a nucleophile (see also discussion and other examples in the literature[39,54−57]), particularly when using an implicit solvent model. That is, including a charged nucleophile creates a major underestimate in the calculated energetics; similar problems have been seen before, both in our simulations of hydroxide attack on the phenyl phosphate dianion,[39] and in independent studies of the hydrolysis of p-nitrophenyl phosphate[55] and acetate[56] (to name a few examples). We believe that this underestimate arises from the combination of a number of factors. The first issue is simply due to the known problems with correctly solvating anionic species and the hydroxide ion in particular using dielectric continuum models (see discussion in the literature[57]). That is, while the precise value of the solvation free energy of the hydroxide ion remains controversial[58,59] due in part to the challenges evaluating ΔsG*(H+) and misunderstandings about corrections for gas/solution phase standard states,[60] dielectric continuum models tend to significantly understimate hydroxide solvation, which would in turn lead to an underestimate of the calculated activation barrier. In combination with problems with incorrect solvation of the hydroxide ion, DFT approaches tend to underestimate barrier heights,[61−64] although this problem is somewhat mitigated by the M06-series of functionals,[65] which was part of our rationale for using this functional in this work. One source of this problem is the delocalization error,[62,66] which refers to the tendency of the approach to artificially spread the electron density out too much. Tying in with this is the static correlation error.[61] These issues can be magnified in large systems, because the delocalization error increases with system size.[62] Additionally, one would expect them to be compounded when more challenging atoms such as sulfur,[67,68] which is polarizable and has low-lying vacant d-orbitals, are introduced into the calculation, and where the quantitative reliability of DFT approaches becomes extremely unpredictable and functional/basis set dependent.[68−70] However, it should be noted that much more reasonable energetics have been obtained for sulfate hydrolysis using a neutral nucleophile,[3] making it likely that the presence of the charged nucleophile is a significant part of the problem. Finally, implicit solvent models are known to underestimate activation barriers due to the fact that they neglect nonequilibrium solvation effects.[71]

Clearly, the majority of these issues are extremely challenging, and resolving them is out of the scope of the present work. However, one trivial source of error that can be addressed is the solvation free energy of the hydroxide ion. That is, as there is no bond formation to the nucleophile in the reactant complex, one can simply adjust the calculated value of the solvation free energy of the reactant complex by the error in the calculated solvation free energy of hydroxide ion compared to the experimental value, thus introducing a constant correction to all calculated solvation free energies (see also related discussion[72]). The only caveat with this approach is the continued uncertainty as to what this solvation free energy actually is. Values in the range of −90.6 to −110.0 kcal/mol[58−60,73−75] have been reported for the solvation free energy of the hydroxide ion in the literature, although at present most sources appear to converge on −104.5 kcal/mol[58,60] as being the most reliable estimate. The value obtained using 6-311+G**/M06-2X/SMD is −97.3 kcal/mol, which is clearly far too low; therefore, we have decided to follow the advice of the literature[58,60] and use −104.5 kcal/mol as the “experimental” solvation free energy for hydroxide, based on the related discussion about the accuracy of values for the free energy of hydration of H+. This results in a correction of 7.2 kcal/mol for the solvation free energy of OH–, which we have added to the ΔΔGsolv values shown in Table 2 (see Table S1 (Supporting Information) for the error and corresponding correction we obtain using other functionals). As can be seen from Table 1, once the calculated solvation free energies are adjusted for this correction, we obtain better quantitative accuracy with experiment, with calculated activation barriers that lie within ∼4 kcal/mol of the experimental value. It should be pointed out that the inclusion of this correction is based on the assumption that the error in the solvation free energy of hydroxide is absent (or reduced by a constant value) at the transition state along this series of compounds. As the distance of the hydroxide ion from the reacting center increases with increasing acidity of the leaving group (see below), there is a risk that this is not completely true. Nevertheless, while obtaining absolute quantitative accuracy is challenging in these cases, the relative trends should still be informative.

To validate our results, we repeated our transition state optimizations and IRC calculations using a number of different functionals with various corrections implemented into them to see how much these would effect our results. Our starting point was just to compare our approach to the standard and popular B3LYP functional,[47−49] which does not include a complete dispersion treatment (in contrast to M06-2X). We then also included CAM-B3LYP,[51] which has been developed to correct the delocalization error, and another dispersion corrected functional, ω-B97X-D,[50] which also includes a correction for long-range effects, as discussed in the Computational Methods. A comparison of calculated total activation energies using these different functionals, as well as the resulting rate constants (obtained using transition state theory), are shown in Table 3 and Figure 8, respectively. From Figure 8, it can be seen that while we have fairly consistent deviations between calculated and experimental values across the series, we are not able to qualitatively reproduce the experimental values and instead obtain very scattered plots. The problem here is that the observed changes in activation barrier across the series are so small (4 kcal/mol between the most reactive and least reactive compounds) and are overshadowed by the error margin in the entropies obtained from the vibrational frequencies. This is supported by the fact that much better qualitative agreement is obtained between theory and experiment if only the trends obtained from ΔE‡ are used, without the addition of any free energy corrections.

Table 3

Calculated Activation Barriers (ΔG‡calc, kcal/mol) and Rate Constants (log kcalc) for the Different Functionals Tested in This Worka

substituted phenol		3-F-4-NO2	4-NO2	4-CN	3-NO2	3-CN	4-Cl	H	3,4-dimethyl
functional	pKa	6.11	7.14	7.95	8.35	8.61	9.38	9.95	10.36
M06-2X	ΔG‡calc	14.5	15.0	15.2	15.5	15.9	16.8	17.3	15.8
	log kcalc	3.00	2.68	2.52	2.32	2.06	1.44	1.11	2.16
B3LYP	ΔG‡calc	19.7	21.4	22.4	21.3	24.3	23.0	23.1	24.6
	log kcalc	–0.52	–1.64	–2.34	–1.55	–3.63	–2.71	–2.81	–3.83
ω-B97X-D	ΔG‡calc	20.2	21.2	21.2	20.8	21.0	22.1	22.3	23.0
	log kcalc	–0.81	–1.50	–1.53	–1.22	–1.36	–2.15	–2.25	–2.73
CAM-B3LYP	ΔG‡calc	19.3	21.6	22.9	21.7	22.3	23.2	23.5	21.7
	log kcalc	–0.25	–1.75	–2.65	–1.82	–2.22	–2.84	–3.10	–1.82
experiment	ΔG‡exp	17.3	18.0	18.1	18.2	18.3	19.3	20.2	20.9
	log kexp	1.12	0.68	0.56	0.51	0.42	–0.25	–0.86	–1.28

The corresponding experimental values (ΔG‡exp and log kexp) are also presented here for comparison. All energies are presented in kcal/mol. Note that the experimental energetics have been corrected for the cost of bringing the reacting fragments from infinite separation into the encounter complex, as discussed in the caption to Table 1.

Figure 8

A comparison of calculated (pink circles) and experimental (blue circles) log k obtained when taking into account just ΔE‡ (A) and also when including zero point energy and entropy corrections (B) using the M06-2X[43] (1), B3LYP[47−49] (2), ω-B97X-D[50] (3), and CAM-B3LYP[51] (4) functionals, respectively.

The corresponding energy breakdowns for each functional are shown in Table S2 (Supporting Information), and a comparison of the changes in S–Onuc/lg distances for each functional are shown in Table S3 (Supporting Information). In all cases, the calculated solvation free energies have been corrected for the experimental solvation free energy of the hydroxide ion, as was done in Table 1 for M06-2X and outlined above. Finally, Figures 9A and B show changes in S–Onuc and S–Olg distances at the transition state for different leaving groups with different DFT functionals, respectively. The first thing that can be observed is that, while there are clearly changes in absolute values upon moving to different functionals, overall trends are not affected by changing the functional. In all cases, we obtain very early transition states with very little bond cleavage to the leaving group (S–Olg distances in a range of 1.70 to 1.80 Å), and some bond formation beginning to occur to the nucleophile (S–O distances in the range of 2.28 to 2.53 Å). Additionally, in all cases, S–Onuc distance increases with increasing acidity of the leaving group, whereas the S–Olg distance remains largely unchanged. This is also reflected in the calculated charges (Figure S5, Supporting Information), which range from −0.5 to −0.8 on the leaving group oxygen atom and −1.08 to −1.32 on the nucleophile oxygen, depending on level of theory and approach used to calculate the charges. Despite the challenges with quantitatively reproducing the experimental activation barriers, the qualitative data these calculations provide all lead to a convergent picture that is in agreement with experimental observables. The crucial factor here is that all functionals tested provide the same trend, making it less likely that this is an accidental observation.

Figure 9

(A) S–Onuc and (B) S–Olg distances at the transition state for different leaving groups with different DFT functionals, respectively.

Exploring the Effect of Explicit Water Molecules on the Transition States

Following on from our calculations using implicit solvation, we have also explored the effect of introducing discrete water molecules into our calculations. Although increasingly popular, there are a number of potential pitfalls associated with such an approach. Most critically, as discussed in ref (37) and references cited therein, it is unclear whether using a mixed implicit/explicit solvation model with a limited number of discrete water molecules reproduces the correct polarization boundary conditions between the explicit solvent and the bulk continuum. If it does not, this would potentially result in overpolarization of the water molecules in the first solvation shell, or, to put it more simply, discrete water molecules embedded into a continuum do not necessarily behave like water molecules embedded into a large, explicit water sphere. In addition to this, including explicit water molecules in the geometry optimization makes it necessary to take into account the entropic cost associated with releasing these water molecules from the artificial QM-optimized positions (which can be quite large[37]). Nevertheless, if treated carefully, as shown by the literature,[55,56] the inclusion of a cluster of explicit water molecules can provide more reasonable energetics when dealing with alkaline nucleophiles, where workers other than us have also observed the continuum model alone to drastically underestimate the calculated activation barrier.

As a starting point, we have performed a full QM(HF) (rather than QM(HF)/MM with the water molecules treated classically) optimization of the phenyl benzene sulfonate intermediate presented in ref (6) using the Cartesian coordinates provided in the Supporting Information of ref (6) to determine whether this is a true stationary point or a Hartree–Fock artifact.[76−78] While the structure presented in ref (6) is not a true minimum using a full QM description in which the water molecules are treated as part of the quantum system (based on frequency calculations which give an imaginary frequency of −142.7633 cm–1), it was nevertheless possible to use this as a starting point to obtain both an intermediate structure as well as transition states for both addition and elimination steps of the reaction. These were then connected to each other by means of IRC calculations to obtain a consistent profile and relevant stationary points (Figure 10). From Figure 10, it can be seen that while this reaction is tending toward a mechanism that is stepwise in terms of bonding pattern when including explicit water molecules, the apparent “intermediate” from the optimization is in fact not thermodynamically stable. Note that although the profiles in Figure 10 do not include a correction for the configurational entropy, this would be expected to be the same for both transition states and the intermediate, and will therefore be unlikely to change the relative energetics of these species. Additionally, as would be expected, as one moves to a better leaving group (e.g., X = 4-CN, Figure 10), the reaction becomes unquestionably concerted and the shoulder observed with the poorer leaving group vanishes. However, it is important to emphasize that even the “intermediate” obtained with the parent compound is not thermodynamically stable and appears to be just an inflection point along the intrinsic reaction coordinate. Additionally (and more critically), it is striking that all the water molecules in the calculation in ref (6) are placed so as to stabilize the nucleophile, with none stabilizing the leaving group. This could create an artifactual intermediate by artificially stabilizing only one side of the reacting species. To test this, we added 2 extra water molecules to stabilize the leaving group. Upon doing this, we were no longer able to obtain a stable intermediate at the HF level of theory, highlighting the dangers of selective microsolvation using explicit water molecules. Overall, it would appear that despite its simplicity, implicit solvation is sufficient to obtain a reliable description of the reacting system in this case (see also discussion in the Supporting Information of ref (41)).

Figure 10

Energy profile for the hydrolyses of (A) phenyl benzene diol and (B) the corresponding 4-cyano substituted compound, at the HF/6-31+G** level of theory, including eight explicit water molecules and the CPCM continuum solvent model. Here, the two aromatic rings were placed in the same conformation as in ref (6), on the basis of having used the coordinates provided in the Supporting Information of this work as a starting conformation.

Discussion of Kinetic Data

The rates of reaction of compounds 1–3 under the same conditions as the data shown in Figure 2 are added to the previously reported data in Figure 11, and the raw data is shown in Figure S6 (Supporting Information). Inspection of Figure 11A shows that the apparent smooth curvature of the Brønsted plot has disappeared to become a linear, somewhat scattered correlation. Figure 11B shows that the new data contributes to a Hammett plot that retains a much better linear correlation. This correlation shows a slight concave-downward curvature, consistent with Hammond movement of the transition state. The much better correlation with σ rather than pKa suggests that the leaving group has not broken its bond to the sulfur to any significant extent in the transition state, and so mesomeric interactions do not have a large impact. These data, which suggest little bond cleavage to the leaving group in the transition state and a small amount of structural variation with leaving group ability, is also consistent with the theoretical picture. According to the calculations, the main changes that occur lie in the interaction of the nucleophile with the sulfur atom. This will be affected by the electronic character of the leaving group oxygen, which will perturb the electrophilicity of the sulfur atom and so affect the reaction rate. The minimal variation in leaving group bond lengths suggests that bond cleavage is not extensive; a slight lengthening is to be expected as the nature of the sulfur atom in the bond is altered. The better leaving groups have a longer bond to the nucleophile in the transition state, suggesting that the nature of the electrophilic sulfur has changed less, and a correspondingly smaller variation with σ is observed in this region of the plot; i.e., there is slightly concave downward curvature of the Hammett plot.

Figure 11

(A) Brønsted correlation for the alkaline hydrolysis of aryl benzenesulfonates (black points reported in ref (6); red points correspond to compounds 1–3). The red dotted line is the linear least-squares fit to all the data, giving β = −0.67 ± 0.07 and R2 = 0.912. (B) Hammett correlation for the same data. The red dotted line is the linear least-squares fit to all the data, giving ρ = 1.61 ± 0.07 and R2 = 0.983.

The ρ value is 1.61; this can be compared with the ρ value for the overall equilibrium process to provide some insight into the reaction progress as described by the effective charge on the leaving group. Williams has measured the β for the equilibrium as 1.8 (which also leads to an assignment of an effective charge of 0.8[79] on the leaving oxygen atom in the starting substrate). This β value can be converted to an equilibrium value for ρ by using the ρ value for the ionization of phenols (2.1), and so the equilibrium ρ = 2.1 × 1.8 = 3.8. This in turn leads to an estimate of the Leffler index for this reaction as 1.61/3.8 = 0.42, which further reinforces the early nature of the transition state. This corresponds to a change of effective charge on the leaving oxygen atom from +0.8 to about 0 in the transition state (0.8 – 0.42 × 1.8 = 0.04). Williams has thoroughly described this method of analysis,[80] showing that the principal assumptions linking the analysis of type II LFER apply to both Hammett and Brønsted correlations. In molecular terms, as the hydroxide attacks the sulfur, the leaving group is converted from a strongly electron-withdrawing group to an essentially neutral substituent. The ρ value reflects this change in character, rather than any significant cleavage of the sulfur–oxygen bond.

The conclusions from this analysis are broadly consistent with the earlier analysis of Buncel, concerning the reactions of sulfonate esters with ethoxide in ethanol.[26,27] These authors show a clearly improved correlation with σ rather than σ– and concluded that this could be intepreted in either the rate-limiting formation of an intermediate or a transition state that closely resembles this. The large ρ values measured for these reactions led Buncel to suggest that the reaction closely resembles the formation of a pentavalent species in the rate-limiting step, or a transition state that closely resembles a putative pentavalent intermediate. Likewise, the recent work of Um and Buncel suggests a small resonance demand from the leaving group in the transition state in anhydrous ethanol is consistent with this description of the reaction character.[42] The experimental and theoretical analysis here suggests that the reaction in water is rather earlier in character.

We note that it is difficult to find aryloxy leaving groups that have pKas lower than 8 that do not have 2-substituents and/or strong resonance interactions with the substituent. These substituted pyridin-3-yl types of compounds are the only readily available aryloxy leaving groups that we have been able to find that are not substituted in the 2 position and have low pKas due to inductive and field rather than resonance effects (and so do not have dissimilar σ and σ– parameters, in contrast to most phenols with low pKas). For many slow reactions, good leaving groups are required for practical purposes: however, it is possible that these data are skewed by the extensive use of substituents that rely on resonance interactions to perturb the pKa significantly. Herschlag and Zalatan[81] have noted that these types of compounds may be outliers in the plot of methyl aryl phosphate diesters and analyzed a series of reactions to justify treating 4-NO2 and 4-CN substitutents as outliers in Brønsted plots; in ref (6), it was noted that perhaps the explanation advanced in terms of intermediate instability might provide an alternative explanation for that data too. However, the data presented here suggests that this explanation is not required. A reasonable approach to correlating the phosphate diester series is to apply a Yukawa–Tsuno (or Young–Jencks) analysis to the data, which reveals a significant (∼33%) resonance demand to be included alongside the inductive and field effects (Figure S7, Supporting Information). Thus, in these reactions, the degree of bond cleavage at the transition state is greater, and delocalization plays a more significant role in its stabilization. This is supported by theoretical calculations on this series using the B3LYP functional,[41] which show shifting transition states with P–Onuc distances similar to those presented in Table S3 (Supporting Information), but with slightly longer P–Olg distances compared to the arylsulfonates presented in the present work. The greater degree of bond cleavage presumably allows for better synchronization of the charge delocalization with the charge development on the leaving group, although again there is limited variation in P–Olg distances across the series. One can reasonably expect a continuum: as the correlation becomes better with σ–, the Hammett ρ or Brønsted β values will increase. At low dependence, one might expect delocalization to be less important (as a result of either an early transition state or imperfect synchronization of charges between nucleophile and leaving group) and Hammett to provide the more appropriate correlation. Combining these parameters through a Yukawa–Tsuno (or Young–Jencks[82]) analysis probably provides a better approach to correlating the effect of varying phenolic leaving groups on rates of reaction. Recent work[83] has proposed the introduction of an additional saturation term (for electron-releasing substituents) in an extended Yukawa–Tsuno equation to analyze phenolates, but there are seldom enough data in these reactivity-based LFER to justify the inclusion of a third term. Furthermore, because for slow reactions the phenols chosen tend to be more rather than less acidic than the parent and so these leaving groups feature less heavily than those with strongly electron withdrawing groups, in the data presented here, leaving groups with strongly electron-releasing groups do not feature. Lastly, in considering how to correlate reactivity data with the electronic properties of a leaving group, we can suggest that the substituted pyridyl leaving groups provide a useful series that extend the series of leaving groups only affected by inductive and field effects down to low pKa, which is likely to be useful in analyzing leaving group ability, perhaps better than, and certainly a complement to, the widely utilized 4-nitrophenyl substituent.

A Revised Mechanistic Picture Based on Computational and Experimental Evidence

On the basis of the computational and experimental evidence presented in this work, the data are most simply analyzed in terms of a concerted pathway, where bond cleavage to the leaving group is not greatly advanced in the transition state. It is possible that there is a transient intermediate, but it is not required and there is no evidence for it.

The previous explanation[6] suggested that there might be a transition from a marginally stable intermediate to one that has no lifetime. This may be true but would not be expected to make a significant difference to the structure of the rate-limiting transtion state, and hence would not lead to markedly different sensitivity to the properties of the leaving group due to this qualitative transition. The scenario here is somewat different to the case where the rate-limiting step leads to an intermediate that becomes too unstable to exist, at which point the reaction is forced to become concerted. This can lead to downward curvature but is similar to a change in rate-limiting step, rather than a change in mechanism.

Conclusions

We find that in all cases the free energy surface provides a single, concerted reaction pathway, with compact transition states that become gradually more dissociative with increasing acidity of the leaving group. In addition, we have explored the effect of including 2, 6, and 8 explicit water molecules in the calculation and demonstrate that the intermediate obtained using HF/6-31++G* has no appreciable lifetime. This intermediate appears to be a simulation artifact, because it also disappears as one adds extra water molecules to stabilize the leaving group as well as the nucleophile. Overall, these results are consistent with our previous studies on the alkaline hydrolysis of methyl arylphosphate diesters,[41] as well as the corresponding fluorophosphates,[39] and experimental work arguing in favor of a rate-determining transition state with little bond cleavage to the leaving group.[26−29] Additionally, we demonstrate, in alignment with the observations of other workers,[68−70] that performing reliable quantitative studies of hydrolysis of sulfur-containing compounds is extremely challenging, particularly if alkaline nucleophiles are involved. However, through validation using several different theoretical approaches, we are reasonably confident with the qualitative information obtained from our calculations. Specifically, on the basis of the free energy surfaces presented in this work, as well as transition states obtained at several levels of theory both with and without the inclusion of discrete water molecules, we demonstrate that the alkaline hydrolysis of all compounds in this series proceeds via a single-step mechanism with a concerted transition state, in line with previous experimental interpretations,[24] arguing against a change in mechanism or rate-determining step as one moves across the series. Finally, exploring the effect of adding explicit water molecules shows that while seemingly a stationary point, the structure shown in ref (6) does not correspond to a thermodynamically stable intermediate even at the HF level of theory, and particularly if additional water molecules are included to stabilize the leaving group in addition to the nucleophile, any hint even of an intermediate vanishes.

The results presented here are qualitatively very similar to those obtained for similar calculations on fluorophosphates[39] and methyl arylphosphate diesters[41] and show a shifting transition state that becomes gradually more reactant-like as the leaving group becomes more acidic. The key difference between these compounds is in the solvation effects, with the neutral sulfonate esters being substantially solvent destabilized at the TS (because of better solvation of the hydroxide ion in the ground state), whereas the transition states of the monoanionic phosphate diesters and fluorophosphates appear to be stabilized by solvent. Therefore, while the transition states for these compounds are geometrically quite similar, as in the related comparison of phosphate and sulfate monoesters,[3] there are some fundamental chemical differences between these compounds, putting significant challenges on unusual enzymes like BcPMH, which can nevertheless catalyze both reactions within the same active site.[15] The precise molecular basis for why this happens is still unresolved, but the present work is an important step toward understanding such catalytic promiscuity at the atomic level.

Experimental Methods

Compound 1 was synthesized from phenyl sulfonyl chloride and 3-hydroxypyridine in the presence of pyridine and triethylamine. 1 was converted to 2 using 3-chloroperbenzoic acid in chloroform, and to 3 by reaction with excess methyl iodide in refluxing acetone.

The rates of reaction of compounds 1–3 were monitored by UV spectroscopy at 50 °C in solutions of KOH with the ionic strength maintained at 0.5 M with KCl. 1: 300 nm, [KOH] 0.01–0.1 M. 2: 320 nm, [KOH] 0.004–0.07 M. 3: 320 nm, [KOH] 0.002–0.008 M. Reactions were initated by adding 5–15 μL of a stock solution in DMSO or dioxan (0.05 M) to 3 mL of KOH solution that had equilibrated at 50 °C. All reactions showed excellent first order behavior, and observed rate constants were obtained by fitting the absorbance change to the integrated first order rate equation. Each compound showed good first order dependence on hydroxide when the observed rate constants were plotted against [KOH] (see the Supporting Information). The rate constant calculated for 0.1 M KOH was used in the Brønsted and Hammett plots shown, to allow direct comparison with the data as reported in ref (6). The pKa values of the leaving groups were obtained from the literature. 3-Hydroxypyridine: 8.42.[85] 3-Hydroxypyridine-N-oxide: 6.45.[86] 3-Hydroxypyridine-N-methyl iodide: 4.96.[85] The Hammett sigma values used for these substituents were as follows: 3-hydroxypyridyl, 0.67;[87] 3-hydroxypyridyl-N-oxide, 1.59;[88] 3-hydroxy-N-methylpyridinium, 2.31.[87] All other sigma values were taken from the literature.[89] See the Supporting Information for a discussion of these parameters.

Pyridin-3-yl benzenesulfonate[90] (1)

Benzenesulfonyl chloride (2.56 mL, 20 mmol) and pyridin-3-ol (1.90 g, 20 mmol) were stirred in THF (50 mL) at room temperature, and triethylamine (3.35 mL, 24 mmol) and pyridine (1.93 mL, 24 mmol) were added dropwise. The reaction was stirred at room temperature for 16 h, the precipitate was filtered off, and the solvent was removed in vacuo. The residue was dissolved in dichloromethane (20 mL) and washed with water (10 mL) and then sodium bicarbonate solution (10 mL). The solution was dried over sodium sulfate before removing the solvent and purifying the crude product by column chromatography on silica (67% 40–60 petroleum ether:33% ethyl acetate) to yield 3.54 g of 1 (74%) as a colorless solid: mp 46–47 °C; δH (400 MHz, CDCl3) 7.32 (1 H, dd, 4.6, 8.3), 7.48 (1 H, ddd, 1.3, 2.7, 8.3), 7.58 (1 H, t, 7.5), 7.73 (1 H, t, 7.5), 7.86 (1 H, d, 7.5), 8.18 (1 H, d, 2.7), 8.53 (1 H, dd, 1.3, 4.6); δC (100 MHz, CDCl3) 124.2, 128.5, 129.4, 130.1, 134.69, 134.74 (C–S), 144.0, 146.4 (C–O), 148.3; ESI-MS positive ion mode m/z 236 (MH+, 100%); HRMS (TOF mode) calculated for C11H10NO3S, 236.0381, found 236.0370.

1-Oxidopyridin-3-yl benzenesulfonate (2)

3-Pyridyl benzenesulfonate (0.50 g, 2.1 mmol) and 3-chloroperbenzoic acid (0.55 g, 3.2 mmol) in chloroform (50 mL) were stirred at room tempertaure for 24 h. The solvent was removed in vacuo, and the crude product was purified by column chromatography on silica (95% ethyl acetate:5% methanol) to yield 0.24 g of 2 (45%) as a white solid: mp 97–98 °C; δH (400 MHz, CDCl3), 7.14 (1 H, d, 8.5), 7.26 (1 H, dd, 6.5, 8.5), 7.62 (1 H, t, 7.5), 7.77 (1 H, t, 7.5), 7.89 (1 H, s), 7.90 (1 H, d, 7.5), 8.13 (1 H, d, 6.5); δC (100 MHz, CD3OD) 122.9, 126.6, 128.4, 129.6, 134.0 (C–S), 134.6, 135.3, 138.1, 148.2 (C–O); ESI-MS positive ion mode m/z 252 (MH+, 100%); HRMS (TOF mode) calculated for C11H10NO4S, 252.0331, found 252.0331.

1-Methyl-3-[(phenylsulfonyl)oxy]pyridinium iodide[90] (3)

3-Pyridyl benzenesulfonate (0.50 g, 2.1 mmol) and methyl iodide (1.32 mL, 21 mmol) in acetone (50 mL) were refluxed for 16 h. The solvent was removed in vacuo, and the residue was dissolved in the minimum amount of methanol. The product was precipitated by the dropwise addition of diethyl ether, filtered, and washed with diethyl ether to yield 0.2 g of 3 (25%) as a yellow solid: mp 144.8–145 °C; δH (400 MHz, D2O) 4.24 (3 H, s), 7.57 (1 H, t, 7.5), 7.76 (1 H, t, J = 7.5), 7.81 (1 H, d, 7.5), 7.90 (1 H, dd, 6.1, 8.7), 8.09 (1 H, d, 8.7), 8.68 (1 H, d, 6.1); δC (100 MHz, CD3OD) 48.4 (CH3), 128.6, 128.8, 130.0, 133.5 (C–S), 135.7, 139.0, 141.0, 144.7, 148.3 (C–O); ESI-MS positive ion mode m/z 250 (M+ – I–, 100%); HRMS (TOF mode) calculated for C12H12NO3S, 250.0538, found 250.0527.