Theoretical Study on the Alkaline Hydrolysis of Methyl Thioacetate in

ARTICLE pubs.acs.org/JPCA

Theoretical Study on the Alkaline Hydrolysis of Methyl Thioacetate in Aqueous Solution Chein-Wei Fu and Thy-Hou Lin* Institute of Molecular Medicine & Department of Life Science, National Tsing Hua University, Hsinchu, 30013 Taiwan, ROC ABSTRACT: A base-catalyzed hydrolysis reaction of thiolester has been studied in both gas and solution phases using two ab initio quantum mechanics calculations such as Gaussian09 and CPMD. The free-energy surface along the reaction path is also constructed using a configuration sampling technique, namely, the metadynamics method. While there are two different reaction paths obtained for the potential profile of the base-catalyzed hydrolysis reaction for thiolester in the gas phase, a triple-well reaction path is computed for the reaction in the solution phase by two quantum mechanics calculations. Unlike the SN2 mechanism (a concerted mechanism) found for the gas-phase reaction, a nucleophilic attack from the hydroxide ion on the carbonyl carbon to yield a tetrahedral intermediate (a stepwise mechanism) is observed for the solutionphase reaction. Moreover, the energy profiles computed by these two theoretical calculations are found to be very comparable with those determined experimentally.

’ INTRODUCTION Thiolesters are common substances involved in many synthetic and biosynthetic reactions such as formation and degradation of fatty acids and acetyl-CoA. Thiolesters also play a crucial role in tagging proteins with ubiquitin for degradation. The carbonyl carbon on thiolesters and other substitutions on this atom are known to be important in both organic and biological processes because of their participations in a wide variety of chemical and biochemical reactions.13 There have been numerous mechanism studies for the carbonyl groups conducted over the years. Among these reactions, the nucleophilic displacement of substitution exchange is the most widely investigated one because of its crucial involvement in the enzymatic catalysis of carboxylic acid and their derivatives. Moreover, the reaction potential surface for the substituent exchange of carbonyl carbon in the gas phase has been found to be a single- or double-well reaction where some tetrahedral structures are characterized as the transition state (TS).46 If the nucleophile is a hydroxide or alkoxide, the potential surface for the substitution or hydrolysis in the gas phase is a single-well reaction, and the tetrahedral structure characterized by a binding energy of 20 kcal/mol is located at the bottom.7 However, if the nucleophiles involved in the reactions are high-electron-affinity atoms such as Cl and Br, the potential surface is rather characterized by a double-well reaction with a tetrahedral TS. On the contrary, a complex well energy profile is obtained for the potential surface if these carbonyl substitution reactions are conducted in the solution phase.8,9 The gas-phase reactions may differ significantly from the solution-phase ones because of lacking treatment of solvation or desolvation effects in the former. The difference between these r 2011 American Chemical Society

two types of reactions has been extensively studied through comparing the reactivity in gas phase with that including one or more solvent molecules to account for the solvation effect. The energy profiles of chemical reactions have been conveniently studied using the ab initio calculation for the electronic structures. More advanced techniques such as those implementing the electronic structure computation into the molecular dynamics simulation have also been devised. Using density function theory (DFT) and nonlocalized plane wave,10 the CarParrinello molecular dynamics (CPMD)11 proposed can effectively separate the slow dynamics of nuclei from that of the fast-moving electrons. The CPMD program has been widely used in studying the chemical and biological reactions. Metadynamics (MTD) is an another advanced method that uses the history-dependent bias potentials for sampling and reconstructing the free-energy surface.12,13 The method is able to fill the reactant and product wells on the free-energy surface and makes a free-energy surface with a higher probability of barrier crossing than conventional molecular dynamics or Monte Carlo simulation methods. The method has also been extensively used in studying the synthesis of organic compounds,1417 peptides,18,19 and the mechanism of enzymatic reactions in biological systems.20,21 The mechanisms of base-catalyzed ester hydrolysis in aqueous solution have been extensively investigated by several theoretical methods. For example, Haeffner et al.22 have studied the basecatalyzed hydrolysis for methyl acetate at the MP2 (Møller Plesset perturbation theory of the second order) level and used Received: May 19, 2011 Revised: October 4, 2011 Published: October 06, 2011 13523

dx.doi.org/10.1021/jp204658w | J. Phys. Chem. A 2011, 115, 13523–13533

The Journal of Physical Chemistry A

ARTICLE

Scheme 1. Three Possible Mechanisms for the Base-Catalyzed Hydrolysis of Thiolestera

a

In the stepwise mechanism, it is similar to the BAC2 (bimolecular base-catalyzed acyl-oxygen cleavage) mechanism, whereas this type of mechanism is often referred to as the reaction of oxygenester. Here, the stepwise mechanism is used instead.

the polarizable continuum model (PCM) to represent the solvent molecules. Both the stepwise and concerted mechanisms are investigated, and two water molecules are included in their calculation.22 Their results show that the methanol elimination step is catalyzed by water molecules. Tantillo and Houk23 have reported a theoretical study for the base-catalyzed hydrolysis of phenyl acetate at the MP2 level by including the solvent effect with the PCM and SCI-PCM (self consistent isodensity surfacePCM) methods. However, only the stepwise mechanism is explored, and the second step of elimination of the corresponding alcohol is found to proceed almost without a barrier. Furthermore, for the hydrolysis of p-nitrophenyl acetate, there is no barrier found in the second step of the hydrolysis reaction, and the reaction progresses as a concerted process. More recent theoretical studies on the base-catalyzed hydrolysis of methyl acetate and methyl formate in aqueous solution have been reported by Zhan et al.24,25 A continuum solvent model and several explicit water molecules are utilized to obtain the activation barriers for both the stepwise and concerted mechanisms.25,26 Their results show that stepwise is the main hydrolysis mechanism

for these RCOOR esters.25 Moreover, several different PCM solvent models are subsequently used to investigate both the stepwise and concerted mechanisms for some alkyl esters.26 Many experimental studies conducted have helped us to understand the mechanism of hydrolysis of esters in solution.8,2739 Some mechanistic experimental works have demonstrated that the hydroxide ion reacts with esters of carboxylic acids through a nucleophilic attack at the carbonyl carbon to form a tetrahedral intermediate (INT) (Scheme 1). This reaction results in losing an alkoxide ion and then is followed by a proton transfer to yield the carboxylate anion and alcohol as the final products.8,2730 Using some isotopes of heavy atoms to investigate the mechanism, Merlier33 has shown that the rate-determining step in the basecatalyzed hydrolysis reaction for thiolester is the formation of a tetrahedral INT structure. However, depending on pH or whether OR is a strong leaving group, the aforementioned mechanism can be switched from a stepwise to a concerted (the SN2 mechanism) one. In the latter case, the tetrahedral INT structure converts to a transition state and the reaction eventually proceeds as an SN2 mechanism.3539 It has also been proposed that the nucleophile 13524

dx.doi.org/10.1021/jp204658w |J. Phys. Chem. A 2011, 115, 13523–13533

The Journal of Physical Chemistry A attack in alkaline solution may be assisted by water molecules, and the reaction progresses as a general base-catalyzed mechanism, as that depicted in Scheme 1. In a general base-catalyzed reaction, a water molecule within the first shell of the hydroxide ion will act as a nucleophile in the reaction and attack the carbonyl carbon by transferring a proton to the hydroxide ion. Moreover, a kinetic study for hydrolyses of formate esters34 has shown that the esters are indeed hydrolyzed by a general basecatalyzed mechanism. In this report, we have employed two theoretical calculations, Gaussian09 and CPMD, to study the base-catalyzed hydrolysis of s-methyl thioacetate in both the gas and solution phases. We treated water molecules explicitly in the quantum mechanics calculation to account for the fact that they may play some important roles as hydrogen bond (HB) acceptors or donors. As shown in Scheme 1, the approaching hydroxide ion and its subsequent nucleophilic attack on the carbonyl carbon of the thiolester are the rate-limiting step of the reaction. This ratelimiting step is studied in both the gas and solution phases using Gaussian09. Further, the free-energy surface of the entire hydrolysis process is studied and constructed using the CPMD program. The solvent effect of water on the hydrolysis process is also studied via a static electronic structure calculation using the 14 water clusters in Gaussian09. Then, the CPMD simulation accompanied by the MTD technique is used to study the hydrolysis process in the bulk water system. The energy barrier of the hydrolysis process thus computed is also compared with those determined experimentally.40,41

’ METHODS Gaussian09. The static electronic calculations for the basecatalyzed hydrolysis reaction for thiolester in both gas and solution phases were performed using the Gaussian09 package. The optimized geometries required for calculation in both the gas and solution phases were obtained from DFT at the B3LYP/ 6-31++G(d,p) level. The energy minima and saddle points searched were confirmed by the frequency calculation, and only one imaginary vibrational mode was allowed for each TS identified. Each TS identified was further confirmed by the intrinsic reaction coordinate (IRC) calculation. The solvation effect in the presence of water on the system was studied using the PCM implemented in the package. CPMD with MTD. A classical molecular dynamics simulation program, namely, NAMD42 was used to prepare the initial system for ab initio molecular dynamics calculations. The initial system prepared consisted of 1 molecule of s-methyl thioacetate (CH3COSCH3), 1 hydroxide ion, 37 water molecules, and a sodium ion added to neutralize the charge. The system was simulated in a NPT ensemble with periodic boundary conditions employed and at temperature 300 K for 1 ns. The dimension of the cubic simulation box created was 10.9 10.9 10.9 Å3, and the average box volume estimated throughout the simulation was 1295 Å3. For further CPMD calculations, the dimension of the cubic simulation box was adjusted to 11 11 11 Å3 to render the water density estimated as 0.993 g/cm3. The ab initio molecular dynamics simulations were performed using the CPMD simulation package. In these simulations, the system was modeled by a periodically repeated cubic cell of length 11 Å, where 1 s-methyl thioacetate, 1 hydroxide ion, 37 water molecules, and a sodium ion were added. All of the hydrogen atoms were substituted by deuterium ones to maintain

ARTICLE

Scheme 2. CV1 and CV2 Distances Defined

the adiabaticity.43 The change of hydrogen to deuterium mass would not affect the calculation results presented. Moreover, all of the hydrogen atoms in the diagrams and figures subsequently presented were not changed to the deuterium ones. The TroullierMartins pseudopotentials44 were used to describe the core electrons, while the Perdew, Burke, and Ernzerhof (PBE) functional was used to describe both the chemically active valence electrons and the electron correlation energy.45 An energy cutoff of 80 Ry was used for the plane-wave basis sets, which was adequate enough to carry out these calculations. The CPMD simulations were conducted in a NVT ensemble at a temperature of 300 K with a NoseHoover chain thermostat employed.46 A fictitious electron mass of 400 amu was used to separate the motion of electrons from that of slowly moving nuclei and also for maintaining the electron adiabaticity. A time step of 4 au (0.1 fs) was used throughout the CPMD simulations. The CPMD simulations lasted for 2 ps to equilibrate the system and also to proceed with the MTD calculation. MTD is a configuration sampling method designed for crossing an energy barrier that is higher than the thermal energy (kBT) for chemical reactions.12,13,47,48 In this method, the free energy is assumed to depend on n collective variables (CVs). The i-th CV collected is denoted as Si. For each CV collected, an auxiliary particle si with fictitious mass m and coupling force constant ki with Si is also assigned. The basic principle of the method lies in the extended Lagrangian described as follows L ¼ L0 þ

∑i mi s2i ∑i 2 ki ðSi si Þ2 V ðsÞ 1

where L0 is the original Lagrangian and is determined from the CPMD calculation. The second term is the kinetic energy of the auxiliary particle si, and the third term is the harmonic potential between Si and si. The last term is the history-dependent bias potentials, which are dependent on the dynamics of si. In a MTD simulation, the system is originally placed at the bottom of the free-energy surface. Then, the bias potential is added to si as a time-dependent force component. The well is gradually filled by the bias potential until the system has enough energy to cross the barrier. The CVs assigned for a typical chemical reaction may contain a reactant well, a TS, and a product well. At the beginning, the system is placed at the bottom of the reactant well. As the MTD simulation proceeds, the reactant well is gradually filled by the bias potential to assist the system in crossing over the TS toward the product well. Then, the product well is filled by the bias potentials again until the system recrosses the energy barrier of the TS. After the simulation, the negative sum of the bias potentials added is used to construct the freeenergy surface. The 2D Gaussian potentials were used in our MTD calculation to construct the free-energy surface for the hydrolysis reaction of thioester in the solution phase. The 2D Gaussian function used was characterized by a width of 0.4 and height of 0.002 au. Both the width and height were scaled by a factor during the simulation 13525



ARTICLE

Scheme 3. TE (trans-Ester) Path of Base-Catalyzed Hydrolysis of Thiolester in the Gas Phasea

In the TE path, the hydroxide and methyl thioacetate (CH3COSCH3) form the first minimum and RC, respectively, and the hydrolysis reaction proceeds with a SN2 mechanism. a

Scheme 4. PE (Proton Exchange) Path of Base-Catalyzed Hydrolysis of Thiolester in the Gas Phasea

In the PE path, the RC structure formed is comprised of an enolate (CH2COSCH3) and H2O. The hydrolysis reaction of the PE path undergoes a SN2 mechanism with a H2O instead of hydroxide as the nucleophile. a

for exploring the narrow wells of possible intermediate states. The total MTD simulation time used was over 10 ps. We chose the distance between the approaching oxygen atom of the hydroxide ion (O1 of the hydroxide ion) and the carbon atom of the thiolester (C1 of s-methyl thioacetate) as our first CV, CV1, because they were both involved in the first step of the basecatalyzed hydrolysis reaction for the thiolester (Scheme 1). The second CV, CV2, was chosen as the distance between C1 of s-methyl thioacetate and S1 of s-methyl thioacetate. Both CV1 and CV2 chosen are depicted in Scheme 2. A force constant k of 0.8 au and a fictitious mass m of 100 au were also chosen for both CVs. The minimum and maximum numbers of CPMD steps executed between the two Gaussian hills were set as 200 and 600, respectively. Moreover, some distance constraints were imposed on the HO bond of the surrounding four water molecules to

prevent proton transfer from the water hydrogen to the oxygen atom of the hydroxide ion.

’ RESULTS AND DISCUSSION Gas-Phase Computation by Gaussian09. The hydrolysis reaction of the thiolester has been less widely studied than that of the ester. Here, we employed several quantum mechanics calculations to study the reaction in both the gas and solution phases for a thorough understanding the mechanism and also to construct a complete potential energy surface for the reaction. We studied the hydrolysis process of the thiolester in the gas phase using the static electronic calculation provided by the Gaussian09 package. These calculations were performed at the B3LYP/6-31++G(d,p) level of theory, and each corresponding 13526



ARTICLE

Figure 1. The optimized structures of paths TE and PE in the gas phase are drawn. RC stands for the reactant complex, TS stands for the transition state, and PC stands for the product complex.

energy minimum and saddle point found were confirmed by the frequency calculation. The TS structures found from these calculations were confirmed by the IRC calculation. We found that there are two different reaction paths, namely, TE (transester) and PE (proton exchange), for the hydrolysis reaction (Schemes 3 and 4) that take place in the gas phase, as presented in both Figures 1 and 2. In the TE path, the hydroxide and methyl thioacetate (CH3COSCH3) form the first minimum, which is the reactant complex (RC), and the hydrolysis reaction proceeds with an SN2 mechanism. In the PE path, the hydroxide does not form a reactant complex with methyl thioacetate directly but rather receives a hydrogen from the methyl group of methyl thioacetate due to the acidity of the α-hydrogen of the methyl group. The RC structure formed in the PE path is comprised of an enolate (CH2COSCH3) and H2O, which is different from that formed in the TE path. Therefore, the hydrolysis reaction of the PE path proceeds with an SN2 mechanism where a nucleophile of H2O is involved instead of a hydroxide. The energy profile of path TE computed has revealed that there are two minima and one TS existing in the path (Figure 2). The first minimum (RC) found is characterized by a structure (Figure 1) where the O1C1 (CV1) and C1S1 (CV2) distances (see also Scheme 2) determined are 3.7 and 1.8 Å, respectively. This is the reactant well where the hydroxide ion is moving closer to the carbonyl carbon of the thiolester. The O1C1 and C1S1 distances determined for the structure of the second minimum (product complex, PC) identified are 1.3 and 3.9 Å, respectively (Figure 1). This corresponds to the product well where an iondipole complex is formed between the methanethiol (CH3SH) and acetic acid (CH3COOH) (Figure 1). The highest point in the energy profile (Figure 2)

Figure 2. The energy profile of paths TE and PE in the gas phase obtained by the DFT B3LYP/6-31++G(d,p) level of theory are compared. In the Reactant, the s-methyl thioacetate and hydroxide ion are located at an infinte distance. RC stands for the reactant complex, TS stands for the transition state, and PC stands for the product complex. In the Product, the methanethiol and acteic acid are also located at an infinite distance. To compare the energy difference between these states, the Reactant is represented as 0 kcal/mol.

computed is identified as the TS of the path and is characterized by a tetrahedral structure where the distances determined for O1C1 and C1S1 are 2.76 and 1.83 Å, respectively (Figure 1). Besides the fact that there is a difference in the energy profile computed (Figure 2), there is also a significant difference in structures for RC and PC identified between path TE and PE (Figure 1, Schemes 3 and 4). For example, in path PE, the 13527



ARTICLE

Figure 3. The optimized structures of model SOL3 computed by the DFT B3LYP/6-31++G(d,p) level of theory with PCM are drawn. RC stands for the reactant complex, TS1 stands for the transition state 1, INT stands for the intermedium, TS2 stands for the transition state 2, and PC stands for the product complex.

Figure 4. The energy profiles of model SOL2, SOL3, and SOL4 computed by the DFT B3LYP/6-31++G(d,p) level of theory with PCM are compared. In the Reactant, the s-methyl thioacetate and hydroxide ion are located at an infinte distance. RC stands for the reactant complex, TS1 stands for the transition state 1, INT stands for the intermedium, TS2 stands for the transition state 2, and PC stands for the product complex. In the Product, the methanethiol and acteic acid are located at an infinite distance. To compare the energy difference between these states, the Reactant is represented as 0 kcal/mol.

O1C1 and C1S1 distances computed for RC are 3.4 and 2.0, while those computed for PC are 1.3 and 3.9 Å, respectively. Being identified in the middle of the energy profiles, the TS

structure of path PE also differs significantly from that of path TE (Figure 1). The O1C1 and C1S1 distances computed for the TS of path PE are 2.7 and 1.8 Å, respectively (Figure 1). Moreover, the RC of path TE is bound via an iondipole interaction and is stabilized by the CH 3 3 3 O interaction from the alkyl and acetyl sides. Therefore, to attack the carbonyl carbon, the incoming hydroxide ion must push away these two alkyl groups. On the contrary, the RC of path PE is stabilized by the C 3 3 3 HO interaction generated by a H2O formed by PE between the hydroxide ion and acetyl group. While the TSs of both paths are all stabilized by the CH 3 3 3 O iondipole interaction, the energy barriers determined between them differ significantly and are 8.4 and 26.9 kcal/mol in paths TE and PE, respectively. A spontaneous PE between the hydroxide ion and methyl group of the thiolester may take place in path PE because an isolated hydroxide ion is not stable. A similar PE process has been observed in the kinetic study of hydrolysis of the oxygen ester system.1 However, the tetrahedral structure of hydrolysis of the thiolester presented here in the gas phase is a maximum, while that of hydrolysis of the oxygenester system is at the minimum of the energy profile. Moreover, for hydrolysis of the thiolester in the gas phase, the reaction proceeds with an SN2 mechanism (a concerted mechanism), and no stable tetrahedral INT is observed along the reaction path. Solution-Phase Computation by Gaussian09. The static electronic calculations for hydrolysis of the thiolester in the solution phase were also conducted with the PCM solvent model accompanied by 14 water clusters implemented in the Gaussian09 package. These calculations were performed at the B3LYP/6-31++G(d,p) level of theory, and each corresponding 13528



ARTICLE

Scheme 5. Pathway of Base-Catalyzed Hydrolysis of Thiolester in the Solution Phase

Table 1. Detail Energy Profile of Each State of Models TE, SOL3, and CPMD Computed by the CPMD-MTD Simulationsa Reactant

RC

TS1

INT

TS2

PC

Product

TE Gas (kcal/mol)

0

24.84

16.36

30.93

45.50

60.07

47.39

SOL3 (kcal/mol) CPMD (kcal/mol)

0 0*

4.45 4.45

8.21 9.35

6.04 11.02

7.67 6.09

22.82 53.06

19.60 19.60*

ΔGsolv (kcal/mol)

55.52

43.13

40.96

39.15

39.19

39.55

47.38

In the TE, SOL3, and ΔGsolv models, the structure of each state is computed at the DFT B3LYP/6-31++G(d,p) level of theory. The “Reactant” consists of the following two molecules, OH(H2O)3 and CH3COSCH3, and the “Product” consists of the following two molecules, CH3COO(H2O)3 and CH3SH. Note that in the Reactant and Product computed by the CPMD model*, it is difficult to estimate OH(H2O)3 and CH3COSCH3 at longer distances or at a distance greater than 7 Å due to the limitation of box size used. Therefore, the energy of Reactant and Product of the SOL3 model is used instead. a

energy minimum and saddle point found was confirmed by the frequency calculation. The TS structures found from these calculations were confirmed by the IRC calculation. The four models with 14 water molecules generated by the PCM solvent model are designated as SOL1, SOL2, SOL3, and SOL4. As shown in Figures 3 and 4 for the energy profiles of these hydrolysis reactions of the SOL2, SOL3, and SOL4 models built, there are three wells (minima) detected with a tetrahedral INT identified in the middle. However, there are only two wells and one saddle point detected in the energy profile computed using the SOL1 model, which is similar to that computed in gas phase (data not shown). For the SOL3 model used and presented in Figure 3, the O1C1 and C1S1 distances computed for RC are 5.2 and 1.78 Å, respectively. This corresponds to the reactant well where the hydroxide ion stays close to the carbonyl carbon of the thiolester (Figure 3). The O1C1 and C1S1 distances computed for the second minimum, which is the INT, are 1.53 and 1.98 Å, respectively (Figure 3). A tetrahedral structure is also observed for the INT, as shown in Figure 3. As characterized by

distances of 1.33 and 7 Å computed, respectively, for O1C1 and C1S1, the third minimum identified is the product well, where a PC is formed by methanethiol and acteic acid with water molecules. Between RC and INT and also between INT and PC, there are two saddle points detected corresponding, respectively, to the two TSs, TS1 and TS2 (Figure 4). The energy barrier determined between RC and TS1 is 13 kcal/mol, and the corresponding O1C1 distance computed for TS1 is 1.96 Å. However, a smaller energy barrier of 2 kcal/mol is computed between INT and PC, and also, a shorter distance of 1.5 Å is computed for the O1C1 distance of TS2. Therefore, the energy of all of the states identified follows the trend TS1 > TS2 > INT > RC > PC (see also Figure 4). Moreover, the static calculation results for the hydrolysis of the thiolester in solution clearly demonstrate the important role played by the water molecules. The tetrahedral INT may be stabilized through strong hydrogen bonding with the surrounding water molecules. Unlike the gasphase reaction between OH and CH3COSCH3 presented in this study, the energy profile computed for hydrolysis of the 13529



ARTICLE

Figure 5. (a) The free-energy surface of hydrolysis of the thiolester in bulk water computed at 300 K using the CPMD-MTD simulation technique. (b) The top view of the free-energy surface computed in (a).

thiolester progresses through a stepwise mechanism, and a tetrahedral INT is observed along the reaction path (Scheme 5). However, the TS that corresponds to a general base-catalyzed mechanism has not been observed in the SOL1, SOL2, SOL3, and SOL4 models used. The solvation energy computed for different steps in model SOL3 is presented in Table 1. Apparently, the solvation energy computed for both the TS and INT structures is similar (Table 1). The hydroxide ion appears to be stabilized more by water molecules, and this may give rise to the higher activation energy for the reaction to take place in aqueous solution compared to that in the gas phase (Table 1). Additionally, there is no significant barrier observed between CH3COO H 3 3 3 SCH3 and CH3COO 3 3 3 HSCH3 in the product well of the gas and solution phases obtained by the Gaussion09 computation. We chose CH3COO 3 3 3 HSCH3 as the structure and energy of PC in Figures 14 because its energy computed is lower than that of CH3COOH 3 3 3 SCH3. Solution-Phase Computation by CPMD with MTD. The free-energy surface for hydrolysis of the thiolester in bulk water is constructed through the CPMD simulation accompanied by the MTD. Approaching the thiolester carbon atom by the oxygen

atom of the hydroxide ion is crucial at the beginning of the hydrolysis reaction or the reactant state. The hydroxide ion is placed at a distance of 4 Å from the carbonyl carbon atom in the reactant state. Proton transfer from a neighboring water molecule to the hydroxide ion after 0.5 ps of simulation is found to take place. In the first solvent shell of the hydroxide ion, the integrity of the hydroxide ion is maintained by imposing some distance constraints to avoid proton transfer between the hydroxide ion and surrounding water molecules.49 The free-energy surface is constructed using the Vreco module of the CPMD package after the system has been subjected to 2 ps of equilibration and 10 ps of MTD simulation steps. The free-energy surface computed at this stage consists of triple wells (reactant, intermediate, and product wells) where a tetrahedral complex is lying in the middle plain. At the first stage of simulation, the reactant well is filled with a bias potential at the 190th step of the MTD simulation to let the thiolester system cross the first saddle point and move through the intermediate toward the product well. Then, after 200 steps of MTD simulation, the product well is filled with the bias potential to make the system move over the second TS toward the INT one. At last, the system is brought back to the reactant 13530



ARTICLE

Figure 6. The snapshots taken at simulation stages A, B, C, and D (see Figure 5b for labels) during the CPMD-MTD simulation process are shown. Stage A corresponds to the 36th meta step at the reactant well where RC is located, stage B corresponds to the 231st meta step at the product well where PC is located, stage C corresponds to the 389th meta step at the INT well, and stage D corresponds to the 482nd meta step where TS1 is recrossed.

well from the INT one by continuously filling of the bias potential. The magnitudes of two CVs are monitored throughout the MTD simulation process. The system is considered to return back to the reactant well when the magnitudes of two CVs monitored are restored back to their initial values. The MTD simulation is considered to be completed at this final stage, and the total number of simulation steps used is 530. Moreover, the wells of the reactant, INT, and product are filled with the bias potential. The free-energy surface of the hydrolysis reaction of the thiolester in bulk water computed by the CPMD-MTD simulation is presented in Figure 5. Apparently, there are three minima and two TSs appearing on the free-energy surface (Figures 5 and 6). The first minimum is identified to be the reactant well (RC) where the hydroxide ion is moving toward the carbonyl carbon of the thiolester (Figures 5 and 6). The corresponding O1C1 (CV1) and C1S1 (CV2) distances computed are 2.2 and 1.9 Å, respectively. The second minimum identified is the INT, where a tetrahedral structure is formed, and the corresponding O1C1 and C1S1 distances computed are 1.2 and 1.7 Å, respectively (Figure 5). With computed O1C1 and C1S1 distances of, respectively, 1.1 and 3.2 Å, the third minimum identified is the product well, where a PC is formed by the methanethiol, acetic acid, and surrounding water molecules (Figure 6). There are also two saddle points appearing in the

energy profile that correspond to the two TSs, namely, TS1 and TS2, where one is located between the reactant well and INT and the other one is located between the INT and product well. The TS1 is formed by the approaching hydroxide ion to the carbonyl carbon atom of the thiolester to generate the CO bond. The corresponding O1C1 distance of TS1 computed is 1.8 Å. While the height of the first energy barrier estimated between the reactant well (RC) and TS1 is 14 kcal/mol, the second energy barrier estimated between INT and TS2 is 4 kcal/mol (Figure 5). The C1S1 distance estimated for the TS2 formed is 2.3 Å, where the SC bond between the S1 and C1 atoms of s-methyl thioacetate is breaking. The solvent effect given by bulk water on the hydroxide ion and the thiolester is also analyzed by determining the formation of HBs during the CPMD-MTD simulation. It is found that at the equilibrium state, strong HBs are usually formed between the hydroxide ion and the surrounding three or four water molecules. However, there are only one to two HBs formed between the carbonyl oxygen of thiolester and some bulk water molecules. The corresponding estimated residence time of these water molecules is shorter than that estimated with the hydroxide ion. In the MTD steps, some HBs are formed between the hydroxide ion and the surrounding four water molecules before the tetrahedral INT is formed. From RC to INT, the number of 13531


The Journal of Physical Chemistry A HBs formed by the same oxygen atom is reduced from four to two. However, there is no significant difference in the HB formed by the carbonyl oxygen of the thiolester and that from bulk water. Furthermore, a tetrahedral INT is observed in the reaction path, and the reaction appears to follow a stepwise mechanism based on the free-energy surface calculated by the CPMD-MTD. In addition, the TS that corresponds to a general base-catalyzed mechanism is not observed in the results given by MTD. The activation energies for the hydrolysis reaction of the thiolester estimated by the Gaussian09-PCM solvent model and the CPMD-MTD computations are 13 and 14 kcal/mol, respectively. Both results computed here are comparable with those measured experimentally, which range from 12 to 13 kcal/mol.40,41 However, the freeenergy surface of the product well (PC) computed by the CPMDMTD is about 20 kcal/mol lower than that computed by the Gaussian09-PCM solvent model. The difference in free energy may be attributed to the effect that desolvation or solvent reorganization caused by moving methanethiol in bulk water toward the acetic acid molecule during the CPMD-MTD simulation steps. The difference in the desolvation energy computed between the Gaussian09-PCM and CPMD-MTD models may be ascribed to the limited size of the simulation box employed in the latter. The box size and total water molecules of the CPMD system used here represent a compromise between the sufficient solvation of the charge ion and the affordable simulation time required to obtain the results. The CPMD system used is comprised of 37 water molecules, and they are filled within the two hydration shells to surround the charge ion. The two water shells added are generally assumed to be adequate for accounting for the free energy of desolvation. However, different numbers of hydration shells added may influence the desolvation energy of the system computed differently.50 In addition, the movement of counterions may be also influenced by the limited box size used. One can use the implicit background charge to replace the counterions, which may influence substantially the solvation structure of the system in a smaller simulation box.51 Although the free energy computed by the Gaussian09-PCM model agrees better with the experimental value, the same mechanism of the hydrolysis reaction in aqueous solution is obtained by both models. Moreover, the Gaussian09-PCM model used does not give insight into the dynamics of the solvent and the important coupling interaction between solvent and solute.

’ CONCLUSION In this study, we have focused on the hydrolysis of the thiolester in both the gas and solution phases by using both the Gaussian09 and CPMD-MTD simulations. As revealed by the potential profile computed by Gaussian09, the hydrolysis for the thiolester in the gas phase proceeds with an SN2 mechanism with two different reaction paths, namely, TE and PE, identified. However, in the solution phase, the hydrolysis reaction for the thiolester takes a stepwise mechanism along a triple-well reaction path as given by the Gaussain09 computation. We also used the CPMD-MTD simulation to account for the solvent effect of the hydrolysis reaction taking place in the solution phase. There are 37 water molecules explicitly included, and all of the molecules are treated with the quantum mechanical calculation. A triplewell process (a stepwise mechanism) as shown by the corresponding free-energy surface computed by this approach is also obtained. In general, the rate-determining step is the nucleophilic

ARTICLE

attack of the hydroxide ion to the carbonyl carbon atom of the thiolester. The second and third wells correspond, respectively, to the INT and product wells. We have also found that the activation energies computed by these two theoretical calculations are in good accord with those measured experimentally. The base-catalyzed hydrolysis of thiolesters in aqueous solution proceeds primarily by the stepwise mechanism, unlike the aforementioned gas-phase reaction between OH and CH3COSCH3, and the solvation effects may have two major influences, namely, (i) the hydroxide ion is stabilized by water molecules, and this will give a higher activation free energy computed for the reaction in the aqueous compared to the gas phase, and (ii) the solvent molecules added tend to stabilize the tetrahedral INT and the structure of the TS in the stepwise mechanism, and this will render the reaction to favor the stepwise mechanism.

’ AUTHOR INFORMATION Corresponding Author

*Fax: 886-3-571-5934. E-mail: [email protected].

’ ACKNOWLEDGMENT This work is supported in part by a grant from the National Science Council (NSC99-2627-B007-005 and NSC99-2628B007-001-MY3), Taiwan. The Gaussian09 and CPMD-MTD computations were conducted at the National Center for High Performance Computing, Taiwan. ’ REFERENCES (1) Chen, X.; Brauman, J. I. J. Phys. Chem. A 2005, 109, 8553. (2) Page, M. I. Organic and Bio-organic Mechanisms; Harlow Longman: London, 1997. (3) McMurry, J. Organic Chemistry, 6th ed.; Thomson Brooks/Cole: Belmont, CA, 2004. (4) Asubiojo, O. I.; Brauman, J. I. J. Am. Chem. Soc. 1979, 101, 3715. (5) Baer, S.; Stoutland, P. O.; Brauman, J. I. J. Am. Chem. Soc. 1989, 111, 4097. (6) Vanderwel, H.; Nibbering, N. M. M. Recl. Trav. Chim. Pays Bas 1988, 107, 479. (7) Wilbur, J. L.; Brauman, J. I. J. Am. Chem. Soc. 1994, 116, 9216. (8) Bender, M. L. Chem. Rev. 1960, 60, 53. (9) Lee, K.; Sung, D. D. Curr. Org. Chem. 2004, 8, 557. (10) Car, R.; Parrinello, M. Phys. Rev. Lett. 1985, 55, 2471. (11) CarParrinello molecular dynamics; copyrighted jointly by IBM Corp and by Max-Planck Institute, Stuttgart. (12) Iannuzzi, M.; Laio, A.; Parrinello, M. Phys. Rev. Lett. 2003, 90, 238302. (13) Parrinello, M.; Laio, A. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562. (14) Sagui, C.; Asciutto, E. J. Phys. Chem. A 2005, 109, 7682. (15) Boero, M.; Ikeshoji, T.; Liew, C. C.; Terakura, K.; Parrinello, M. J. Am. Chem. Soc. 2004, 126, 6280. (16) Boero, M.; Tateno, M.; Terakura, K.; Oshiyama, A. J. Chem. Theory Comput. 2005, 1, 925. (17) Ensing, B.; Laio, A.; Gervasio, F. L.; Parrinello, M.; Klein, M. L. J. Am. Chem. Soc. 2004, 126, 9492. (18) Nair, N. N.; Schreiner, E.; Marx, D. J. Am. Chem. Soc. 2008, 130, 14148. (19) Schreiner, E.; Nair, N. N.; Marx, D. J. Am. Chem. Soc. 2008, 130, 2768. (20) Babin, V.; Roland, C.; Darden, T. A.; Sagui, C. J. Chem. Phys. 2006, 125, 204909. (21) Boero, M.; Ikeda, T.; Ito, E.; Terakura, K. J. Am. Chem. Soc. 2006, 128, 16798. 13532



ARTICLE

(22) Brinck, T.; Haeffner, F.; Hu, C. H.; Norin, T. J. Mol. Struct.: THEOCHEM 1999, 459, 85. (23) Tantillo, D. J.; Houk, K. N. J. Org. Chem. 1999, 64, 3066. (24) Zhan, C. G.; Landry, D. W.; Ornstein, R. L. J. Phys. Chem. A 2000, 104, 7672. (25) Zhan, C. G.; Landry, D. W.; Ornstein, R. L. J. Am. Chem. Soc. 2000, 122, 2621. (26) Zhan, C. G.; Landry, D. W. J. Phys. Chem. A 2001, 105, 1296. (27) Bender, M. L.; Dewey, R. S. J. Am. Chem. Soc. 1956, 78, 317. (28) Bender, M. L.; Ginger, R. D.; Unik, J. P. J. Am. Chem. Soc. 1958, 80, 1044. (29) Shain, S. A.; Kirsch, J. F. J. Am. Chem. Soc. 1968, 90, 5848. (30) Johnson, S. L. In Advances in Physical Organic Chemistry; Gold, V., Ed.; Academic Press: New York, 1967; Vol. 5, p 237. (31) Humphreys, H. M.; Hammett, L. P. J. Am. Chem. Soc. 1956, 78, 521. (32) Bruice, T. C.; Holmquis., B J. Am. Chem. Soc. 1968, 90, 7136. (33) Marlier, J. F. J. Am. Chem. Soc. 1993, 115, 5953. (34) Stefanidis, D.; Jencks, W. P. J. Am. Chem. Soc. 1993, 115, 6045. (35) Basaif, S.; Luthra, A. K.; Williams, A. J. Am. Chem. Soc. 1989, 111, 2647. (36) Guthrie, J. P. J. Am. Chem. Soc. 1991, 113, 3941. (37) Hengge, A. C.; Hess, R. A. J. Am. Chem. Soc. 1994, 116, 11256. (38) Fernandez, M. A.; de Rossi, R. H. J. Org. Chem. 1999, 64, 6000. (39) Williams, A. Acc. Chem. Res. 1989, 22, 387. (40) Noda, L. H.; Kuby, S. A.; Lardy, H. A. J. Am. Chem. Soc. 1953, 75, 913. (41) Rylander, P. N.; Tarbell, D. S. J. Am. Chem. Soc. 1950, 72, 3021. (42) Schulten, K.; Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L. J. Comput. Chem. 2005, 26, 1781. (43) Blochl, P. E.; Parrinello, M. Phys. Rev. B 1992, 45, 9413. (44) Troullier, N.; Martins, J. L. Phys. Rev. B 1991, 43, 1993. (45) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865. (46) Nose, S. J. Chem. Phys. 1984, 81, 511. (47) Micheletti, C.; Laio, A.; Parrinello, M. Phys. Rev. Lett. 2004, 92, 170601. (48) Qian, X. H.; Dong, H. T.; Nimlos, M. R.; Himmel, M. E.; Johnson, D. K. J. Phys. Chem. A 2009, 113, 8577. (49) Sprik, M. Faraday Discuss. 1998, 110, 437. (50) Chandler, D.; Geissler, P. L.; Dellago, C. J. Phys. Chem. B 1999, 103, 3706. (51) Blumberger, J.; Ensing, B.; Klein, M. L. Angew. Chem., Int. Ed. 2006, 45, 2893.

13533


Theoretical Study on the Alkaline Hydrolysis of Methyl Thioacetate in

Recommend Documents