2312
J. Phys. Chem. 1995, 99, 2312-2323
Computer Modeling of Phenyl Acetate Hydrolysis in Water and in Reaction.with p-Cyclodextrin: Molecular Orbital Calculations with the Semiempirical AM1 Method and the Langevin Dipole Solvent Model Victor B. Luzhkov’ and Carol A. Venanzi* Department of Chemical Engineering, Chemistq, and Environmental Science, New Jersey Institute of Technology, Newark, New Jersey 07102 Received: September I , 1994; In Final Form: November 29, 1994@
The Langevin dipole (LD) solvent model is used with the AM1 semiempirical quantum mechanical method to compare reaction pathways in the gas phase and in polar solution for phenyl acetate cleavage by hydroxide ion and by P-cyclodextrin. The gas-phase results show that nucleophilic substitution at the acyl carbon atom of phenyl acetate by the hydroxide ion exhibits a downhill potential profile. The introduction of solvent effects using the LD model results in correct reproduction of the experimental values for the activation bamer and the heat of hydrolysis. The reaction of the alkoxide ion of cyclodextrin with bound phenyl acetate in the gas phase is found to have a positive energy of activation, which is further increased in the presence of polar solvent. The results show that acylation at the 3’-hydroxyl of cyclodextrin is favored over the 2’-position by about 15 kcaYmol due to less structural reorganization of the macrocycle during hydrolysis at the 3’-site. The results also show that attack at the 3’-hydroxyl lowers the activation energy barrier by about 10 kcaumol compared to attack by hydroxide ion in solution. The results are discussed in terms of the known physicochemical properties of cyclodextrins and their inclusion complexes.
Monte Carlo simulations of the reaction of H2CO with OHshow that the nucleophilic attack by OH- at the acyl carbon Recently there has been considerable interest in the computer proceeds along a downhill potential in the gas phase, while the modeling of carb0hydrates.l Cyclodextrins are of particular barrier that is observed experimentally in polar solvent is interest due to their extensive applications in both research and introduced by the changes in solvation energy along the reaction path,29c,29f.30.3 1 industry.2-6 Cyclodextrins act as hosts for a variety of included guests7-14 and have been used as a structural template for the The cyclodextrins are known to accelerate the rate of acylation design of artificial enzymes which mimic ~hymotrypsin’~-’~ and of a number of different substrates. Experimental observations other enzymes. l8 Inter- and intramolecular hydrogen-bonding show that cyclodextrins react with phenolic esters via an patterns of cyclodextrin in water have been studied by X-ray alkoxide ion derived from the secondary hydroxyl groups with and neutron diffra~tion,’~ as well as by computer simulation.20 the formation of a covalent intermediate and a subsequent release Molecular mechanics calculations have been used to study of corresponding phenol^.^-^ However, it is not clear whether inclusion complexes of cyclodextrins and various g ~ e s t s . ~ l - ~ ~the reaction occurs at the 2’-hydroxyl or the 3’-hydroxyl. A Previous work from this laboratory has involved gas-phase bound substrate has been observed to donate a nonequilibrating molecular mechanics and dynamics studies of P-cyclodextrin tosyl group to the 2’-hydr0xyl.~~On the other hand, direct (cycloheptaamylose) in order to clarify the steric and electrotosylation of the 3’-hydroxyl has also been noted.34 In addition, static features of the macrocycle that are important to its function a number of metallocene substrates have been synthesized by as an artificial enzyme26 and to investigate the conformational Breslow and co-workers16bscto show the importance of optimal potential energy surface of the molecule.27 In the present work, steric alignment between the cyclodextrin and a substrate in we extend our previous studies considerably by using molecular order to achieve higher reaction acceleration. orbital theory and the Langevin dipole solvation model to study In order to provide further insight into the reactivity of the distinguishing feature of P-cyclodextrin as an artificial cyclodextrins, we present the results of semiempirical molecular enzyme: its ability to increase the rate of acylation of ester orbital (MO) calculations of the reaction path of ester hydrolysis substrates. by OH- and the alkoxide ions of P-cyclodextrin in gas-phase The hydrolysis of phenolic esters is assumed to proceed by and polar environments. The major questions considered in this the mechanism of nucleophilic substitution both in aqueous study are as follows: (1) What is the role of the microenvisolution and in cyclodextrin c o m p l e x e ~ .A~ number ~ ~ ~ of simple ronment provided by the macrocycle cavity in the hydrolysis reactions of nucleophilic substitution at a carbonyl group have of phenyl acetate? (2) Is there a difference in the reactivity of been extensively studied in recent years both in gas-phase and the 2‘- and 3’-secondary hydroxyl oxygen atoms of p-cycloaqueous media using quantum chemical and statistical physics dextrin? To our knowledge, no quantum chemical calculations appro ache^.^^-^* The general feature of different reactions of of ester hydrolysis either in solution or by a cycloamylose have this class is the formation of a tetrahedral intermediate that been previously performed. undergoes further disruption of the C - 0 bond. Ab initio and Methods
Introduction
* Author to whom correspondence should be addressed. Permanent address: Institute of Chemical Physics. Chemogolovka, Moscow Region, Russia 142432. Abstract published in Advance ACS Abstracts, February 1, 1995. @
0022-365419512099-2312$09.00/0
The MO calculations were performed with the semiempirical AMI method35 in the MOPAC-93 p r ~ g r a musing ~ ~ , an ~ ~IBM RSl6000 workstation. While rigorous ab initio techniques are
0 1995 American Chemical Society
Computer Modeling of Phenyl Acetate Hydrolysis preferable for studies of small molecules, semiempirical methods present a reasonable alternative for the study of medium- and large-sized systems. The choice of the AM1 method was determined by its computational efficiency, well-documented accuracy (comparable with ab initio 6-3 lG* calculation^^^) for molecules containing H, C, and 0 atoms, and several advantages compared to the other semiempirical methods, MNDO and PM3. For example, due to an improved description of the hydrogenrepulsion terms, Ah41 reproduces the steric interactions and gasphase transition-state energies better than MND0.39-4* In addition, in contrast to the MNDO method, the AM1 method has the ability to reproduce correctly the energies of hydrogen bonds,36,37,42-44 which are abundantly present in the cyclodextrin structure. While being less accurate than the PM3 method in reproducing the geometries of intermolecular hydrogen-bonded complexes,29athe AM1 method more adequately reproduces the intramolecular hydrogen bonds45and the 0-H bond energies.& For all molecules, full optimization of the molecular structure without any symmetry conditions was carried out using the MOPAC keywords PRECISE and GNORM=O.l (for phenyl acetate hydrolysis in water) and GNORM=l .O (for cyclodextrin and its complexes). The gas-phase-optimized structures were used subsequently in the calculation of solvation energies. Since the systems studied do not belong to a very polarizable class of molecule, the decoupling of the quantum mechanical and the solvation calculations is not expected to have a large effect on the calculated geometries and energies along the reaction path. The study of these molecular systems involved extensive use of molecular graphics within the program CHEM-X!7 The calculations of the solvent effects on the solute electronic structure and of the corresponding hydration enthalpies were performed with the microscopic Langevin dipole (LD) mode1.31,48s49In this approach, the solvent system is divided into two regions by a large sphere around the center of the solute. The solvent molecules within this sphere are treated as polarizable point dipoles, while the solvent in the outer region is treated as a continuum using the Bom reaction field formulae. The average polarization of the solvent dipoles is related to the corresponding local electric field in a self-consistent iterative way using the Langevin-type equation. The solvent term was added to the AM1 Hamiltonian in the MOPAC-93 program by one of us (V.B.L.) in the same way as it was previously incorporated into the AMPAC program.49 The electrostatic solvation energy is given by
where Epem,Eind, and EB~,,,are the contributions to the LD electrostatic energy from the interactions with permanent and induced dipoles and bulk solvent, respectively. The term AHgas-solgives the change in the solute intemal energy (the heat of formation calculated with AM1 method) as a result of its polarization by the solvent molecules. The solvation energy calculated with this approach explicitly takes into account the mutual solute-solvent polarization and the solvent reorganization en erg^.^^,^^ The stability of the reactive species in solution, AHso,,is characterized by the sum of the gas-phase enthalpy of formation and the solvation energy
The method of construction of the solvent dipole grid around the solute followed that described in ref 49. The positions of the solvent dipoles in the inner solvation region were assigned to a cubic grid with step 3.1 A. These grid points were confined by a sphere of radius 19.0 8, around the center of mass of the
J. Phys. Chem., Vol. 99, No. 8, 1995 2313 solute. The dipoles in the first solvation shell were placed at points over the corresponding solvation envelope (surface accessible to the solvent) around the solute. This surface was defined as a superposition of spheres around the solute atoms with the radii equal to the sum of the water radius and the corresponding atomic van der Waals' radius.49 For all molecules, the solvent dipoles were specifically placed at the lone pair sites of the oxygen atoms. The remaining part of the solvation surface was equally covered with the solvent dipoles with the interdipole spacing equal to the step of the cubic grid. The energy of solvation was taken as an average over six different configurations of the dipole distribution over the molecular solvation surface and the corresponding outer cubic grid. The cutoff radius for dipole-dipole interactions was taken to be equal to the radius of the solvation sphere, 19 A. Calculation of the solvation energy was performed using two sets of AM1 atomic charges: conventional (Coulson) charges, Qc, and electrostatic potential-derived charges, QPD. The latter charges were obtained using the procedure of Besler et al.50 given in the MOPAC-93 program, where they were derived by fitting the charges to reproduce the molecular electrostatic potential calculated at solvation envelopes around the molecule. These charges more accurately reproduce the spatial distribution of the electric field and thus might be expected to give better results than the conventional charges when used in the LD model. Both sets of charges were scaled with the corresponding coefficient^^^ in order to reproduce high-level ab initio 6-31G* potential-derived charges. The necessity of using this more realistic intermolecular charge polarization in calculations with the LD model was previously shown in ref 49. the value of the solvent dipole moment in the LD model was adjusted from 1.85 to 1.63 D in order to reproduce the free energies of solvation of ions and polar molecules. Typically, for ions the enthalpy contribution differs only by 5% from the corresponding free-energy We have recently shown that the electrostatic contributions to the relative hydration free energy and relative hydration enthalpy differ by 4% or less for the free base and protonated conformers of the diuretic drug, amil~ride.~'For nonpolar, flexible, neutral molecules, the free-energy term contains also the important entropy parts associated with the hydrophobic effect and with translational and rotational degrees of freedom of the solute. Thus, despite the use in this work of the scaled water dipole moment, 1.63 D, the calculated values of electrostatic solvation energies should be considered to be an estimate of changes in hydration enthalpies along the reaction curve, rather than of hydration free energies. The electrostatic part of the solvation enthalpy is the major contributor to the relative energy changes in the process of ionic reactions of hydrolysis in polar solution.30,31Additional evidence for the prime role of the enthalpy factor in ester hydrolysis is also given in ref 32a,b. 1. Computer Modeling of Phenyl Acetate Hydrolysis in Water. The mechanism of base-catalyzed phenyl acetate hydrolysis is given in Scheme 1. The hydrolysis rate is determined by the reaction of phenyl acetate, 1, with the hydroxide ion, 2 (step 1) via the transition state, 3, giving an unstable tetrahedral intermediate, 4. Step 1 is followed by fast cleavage of the acyl carbon-oxygen bond in 4 (step 2), leading to the reaction products phenoxide ion, 5, and acetic acid, 6, and finally to phenol, 7, and acetate ion, 8, in step 3. Fully-optimized equilibrium molecular structures were determined for all reactive species given in Scheme l. In exploring the reaction path for step 1, the CS-011 bond length was taken as a reaction coordinate. At each value of the reaction coordinate, complete optimization was carried out for all
Luzhkov and Venanzi
2314 J. Pliys. Chem., Vol. 99, No. 8, 1995 SCHEME 1 H.
H
CH3 I
'0H H
4
3
2
1
Step 2
Step 1
HO
\
-$ + '0
\
//c-cH3
0
5
6
7
8
Step 3 geometrical variables, except that the symmetry of the acylCH3 group was kept fixed in order to avoid crossing the potential surface of the hydrogen-abstraction reaction. The transitionstate geometry was recalculated starting from that point by minimization of the gradient norm. The Hessian of this minimized structure was calculated and determined to have one negative frequency. Then, the energies of the reactants and the products given in Scheme 1 were calculated in a polar environment using the LD solvent model and the fixed, gas-phase-optimized geometries. Since there exists extensive experimental data on solvation properties for most of the compounds given in Scheme 1, the calculations provide an additional test of the reliability of solvation energies calculated with the LD model. 2. Computer Modeling of Phenyl Acetate Cleavage in the Reaction with P-Cyclodextrin. The mechanism of phenyl acetate hydrolysis in the reaction with the alkoxide ion of /3-cyclodextrin is given in Scheme 2. This includes the binding of the substrate 1 to /3-cyclodextrin, 9 (step l), ionization of one of the secondary hydroxyl groups of the inclusion complex 10 (step 2), formation of the tetrahedral adduct 13 in the intermolecular reaction (step 3) of the substrate with an alkoxide ion of the secondary hydroxyl group (via the transition state 12), cleavage of the C - 0 bond in 13 with the formation of acylated cyclodextrin, 14, and the phenoxide ion, 15, in step 4. Here, we have not considered the formation of 11 via complexation of 1 with the /3-cyclodextrin ion, since, due to pKa = 12.1,15aP-cyclodextrin exists in solution mainly in the neutral form. As above, the geometries of the reactants and products were first optimized in the gas phase. Then, the energies in the polar solvent were calculated with the LD method using these geometries. Geometry Optimization of p-Cyclodextrin and Its Alkoxide Ions. The starting geometry in the AM1 optimization of /3-cyclodextrin was taken from the neutron diffraction structure
of P-cyclodextrin undecahydrate. 19a In this structure, the hydroxyl groups of the cyclodextrin do not follow a homodromic pattern (with cyclic directionality, i.e., all hydroxyls pointing in one direction) due to the formation of extensive hydrogenbonding networks with the solvent. In geometry optimization of the isolated, neutral molecule 9, all the hydroxyl groups were aligned in homodromic fashion in the starting geometry. This is shown schematically in Figure 1, which also gives the numbering of the cyclodextrin glucose residues and the crystallographic notation for the atoms of each residue. Two conformations of 9 were investigated: 9a corresponds to the fully-optimized geometry of the neutron diffraction structure; 9b corresponds to the "symmetric" conformation obtained by imposing symmetry conditions on each of the seven sets of interand intraglucose torsional angles and on the primary and secondary hydroxyl torsional angles during the optimization. In other words, angles of each type were set equal to each other for all seven residues and were varied during the optimization. The optimized geometry of the &oxide ion of 9 was obtained for the anion at the 0 2 and 03 secondary hydroxyl oxygen position of residue 1. Due to ionization of the hydroxyl group, the homodromic pattern of secondary hydroxyl orientations is broken in the anions of 9 since in this case one pair of adjacent oxygen atoms is not connected by a hydrogen bond. In the f i s t step, the most stable orientations of the secondary hydroxyls of the anions of 9 were determined. The conformations with different mutual orientations of secondary hydroxyl groups were considered, and partial optimization of only the torsional angles involving the secondary hydroxyls was performed. (The structure of the macrocycle was taken from the structure of the 9a conformation of the neutral form.) Thus, the conformation designated as 9a'(02-) corresponds to structure 9a with the anion formed at the 0 2 position of residue 1, which gives the most stable pattern of the mutual orientations of the secondary hydroxyls. Conformation 9b'(03-) is similarly defined. The
J. Phys. Chem., Vol. 99, No. 8, 1995 2315
Computer Modeling of Phenyl Acetate Hydrolysis
SCHEME 2 CH3
I
+
9
1
10
Step 2
Step 1
11
13
12
Step 3
CH3-C
P \
0
+ 14
8 5
Step 4 subsequent conformational study of the anions was performed in two ways. 9a(02-) was obtained from a full geometry optimization of the anion using 9a’(02-) as a starting structure. 9b(02-) was obtained in the same way as 9a(02-) except that 9b’(o2-) was used as the starting structure. Similar calculations with the anion located at the 0 3 position resulted in 9a’(03-), 9a(03-), and 9b(03-). Calculation of the Inclusion Complexes of Phenyl Acetate with B-Cyclodextrin and Its Ions. In the study of inclusion complexes 10 and 11, we considered the orientation of the substrate having the acyl group pointed out of the cavity of 9. This is the type of structure generally observed in the inclusion complex of /3-cyclodextrin with benzenes with hydrophilic substituent^.^^ The starting geometry for the optimization of 10 and 11 was the optimized geometry of 9a, Le., 9a(02-) and 9a(03-), respectively. In the calculation of the structure and stability of the alkoxide ion inclusion complexes 11, we considered all 14 possible ionization sites due to the 7 sets of secondary 0 2 and 0 3 hydroxyls in the optimized structure of 10, since for these ions, the 7-fold symmetry of the structure of
the neutral inclusion complex is lost. A full geometry optimization was performed for each complex. Calculation of the Tetrahedral Adducts of Phenyl Acetate with Alkoxide Ions of B-Cyclodextrin. Two types of tetrahedral adduct, 13, were investigated in which a covalent bond was formed between atom cs8 of the substrate, where the superscript “s” indicates a substrate atom, and either 0 2 - or 03- of the cyclodextrin. Only configurations that have the benzene ring of the substrate pointed inside the cavity were considered since these structures should be the closest to the transition state along the reaction path from the ionized complex 11 to the product 13. For the same reason, as an additional restriction, only the conformers of 13 with the orientation of the benzene ring plane perpendicular to the cavity macrocycle were considered. In the calculation of the structure of 13, the starting geometry of the macrocycle was taken from the 9a(02-) and 9a(O3-) conformations, respectively. The initial positions of the substrate with respect to the macrocycle cavity were generated using the “modify geometry” facility of CHEM-X. Namely, the cs8-
Luzhkov and Venanzi
2316 J. Phys. Chem., Vol. 99, No. 8, 1995
TABLE 1: Selected Geometrical Parameters of the Optimized Structures of 1, 3, and 4" dihedral angle' compd bond lengthbor bond angle' 1
3
4a
7-8 8-9 1-7-8 7-8-9 7-8 8-9 8-11 1-7-8 9-8-1 1 7-8 8-9 8-11 1-7-8 9-8-11
1.377 1.228 120.3 119.8 1.396 1.230 2.680 118.0 97.0 1.519 1.283 1.421 122.8 112.0
1-7-8-9 6-1-7-8 7-1-6-5 7-9-8-10 1-7-8-9 1-7-8-11 6-1-7-8 7-1-6-5 9-8-11-H(11) 1-7-8-9 1-7-8-11 5-6-1-7 6-1-7-8 9-8-11-H(11)
1.o -47.7 174.5 180.0 39.7 -64.8 49.6 174.6 23.4 47.9 -70.7 179.4 -6.7 -23.2
Atom numbering is given in Scheme 1. In angstroms. In degrees.
4
0-
Figure 1. Schematic representation of /3-cyclodextrin with primary and secondary hydroxyls arranged in homodromic fashion. Glucose residues are numbered from 1 to 7. Crystallographic notation is used for the atom labels. 0 2 (CS8-O3) distance was initially set equal to 1.50 8, and the orientations of the chemical bonds of the tertiary cs8carbon atom were set to give approximately an sp3 hybridization. Calculation of the Reaction Path of Nucleophilic Attack in Cyclodextrin Complexes. The simulations of the reaction path of the formation of the tetrahedral adduct 13 were performed by a straightforward scanning of the distance between the acyl carbon, cs8,of the substrate and the oxygen atom of the alkoxide ion with the relaxation of the geometry of the reactive complex. Our attempts to directly follow the reaction path by minimizing the energy of the complex along the reaction coordinate (starting from the structure of the inclusion complex 11 and moving to smaller O-Cs8 distances) gave an unphysically high potential barrier due to large steric interactions between the substrate and the host. Thus, we considered a more clearlydefined reverse process-an abstraction of 1 from the tetrahedral intermediate 13. A full optimization of the product along the reaction curve was performed except that the bond lengths of the six residues in 13 that are not directly involved in chemical bonding with the substrate were kept fixed at their values in 13. The value of the reaction coordinate was changed in increments of 0.1 8, from the equilibrium value of the O-Cs8 distance in 13 up to 3.1 8,. The activation energy, Ea, was estimated from the maximum energy value on the reaction curve. Unlike the modeling of Scheme 1 with OH-, here rigorous determination of the transition-state geometry by gradient minimization was not performe$ due to computational limitations. It was found that at 3.0 A, the substrate already adopts its regular structure and the negative charge is completely localized at the oxygen atom of cyclodextrin. The rest of the curve (corresponding to an initial stage of the attack of substrate on the ionic center) was not considered due to uncertainties involved in the gas-phase modeling of this process. This gap on the reaction curve might be of some minor importance since the reaction barrier is determined by the interaction at a 2-3-A separation of cs8and 0 atoms.
Results and Discussion 1. Computer Modeling of Phenyl Acetate Hydrolysis in Water. Selected geometrical parameters from the optimized
1 3 4 Figure 2. Optimized geometry of 1, 3, and 4a. Selected optimized parameters are given in Table 1.
structures of 1, 3, and 4a are given in Table 1 and in Figure 2. Table 1 shows that the equilibrium conformation of 1 is characterized by almost planar structures of the benzene ring = 174.5') and acyl group ( O ~ - C ~ - O ~ - C ~ O (07-C1-C6-C5 = 180.0'), with the out-of-plane dihedral angle C6-Ci-07c8 equal to -47.7'. Three stable conformers of the product 4 were found. Each had c8 in the plane of the aromatic ring with either (210, 0 9 , or 0 1 1 in the trans position with respect to C1. These conformations are denoted by 4a, 4b, and 4c, respectively. For the tetrahedral intermediate in its most stable conformation, 4a, atom c8 lies in the plane of the aromatic ring and atom Clo of the methyl group is in the trans position to atom C1. The conformations that are obtained by rotation around the 07-Cg bond with atoms O9 or 0 1 1 in the trans position to atom C1 (4b and 4c structures, respectively) have energies higher than the global minimum by 0.6 and 2.4 kcaymol, respectively. The calculated geometrical and thermodynamical characteristics of the isolated acetic acid, phenol, phenyl acetate, and their ions are very similar to the previously reported results of the AM1 calculations of these compound^.^^^^^* Table 2A gives the gas-phase heat of formation and the electrostatic solvation energy for compounds 1-8 in Scheme 1. These values are used along with eq 2 to calculate the activation energy and enthalpy of the reaction for steps 1-3 in Scheme 1. These results are given in Table 2B. The data in Table 2B along with other points along the reaction curve for Scheme 1 are plotted in Figure 3. Generally speaking, the gasphase curve (open boxes) can be considered as a downhill descent from the separated reactants 1 and 2 to the product 4. A shallow local minimum in this curve is observed at the distance of 3.5 A, which can be attributed to theoformation of the ion-dipole complex. The small peak at 2.65 A corresponds to the transition state of the reaction. The fi!al value of the transition-state reaction coordinate equals 2.68 A (Table 1). The
J. Phys. Chem., Vol. 99, No. 8, 1995 2317
Computer Modeling of Phenyl Acetate Hydrolysis
TABLE 2: Energies and Enthalpies for Scheme 1
4
3
5,6
798
2A. Gas-Phase Heats of Formation" and Electrostatic Solvation Energies" of 1-8 compd
m H e a s
1 2 3
-55.3 -14.1 -92.4 - 142.9 -41.1 -103.1 -22.3 -115.5
4 5
6 7 8
A&odQC) -26.1 -107.0 -84.1 -76.6 -74.6 -19.2 -13.4 -87.1
AESOI(@~) -15.9 - 107.1 -82.2
-75.7 -74.5 -11.0 -11.1 -86.5
2B. Calculated and Experimental Activation Energies'sb and Enthalpiesovcof Steps 1-3 in Scheme 1 reaction AM1 AMl/LD(QC) AM1/LD(QPD) expt step 1 Ea -23.0 26.0 17.8 11.9d @ -73.5 -17.0 -26.2 step 2 @ -1.3 -18.5 -11.1 step 3 M 6.4 -0.3 -5.7 -8.3' "In kcal/mol. *Activation energies, Ea, were calculated as the difference in energy between 3 and 1 plus 2 using the data in 2A. The AM1 column is from the AHgasdata. The AMl/LD(QC)and AMl/ LD(QpD)columns are from eq 2 using the AHgasdata and either the AEs,l(Qc) or the AEsol(QpD) data. The enthalpies of reaction, @, were calculated as the difference between the enthalpies of the products and reactants in a similar fashion to E, above. From ref 56. e Estimated using the data for the relative enthalpies of solvation for the (molecule/ ion) pairs 7,5,68.7 kcal/moP3.", and 6,8,77.0 kcaYm01.*~The relative solvation enthalpies for 6, 8 were taken from the work of Wilson et al.53 Those for 7,swere estimated using the approximation of Wilson et al.53 and the data from the work of Pearsod4 and Christensen et a1.55
transition state on the potential surface of step 1 giving the 4b conformer coincides with that of 4a. Figure 3 also gives the AMl/LD(QC) (open circle) and AM1/ LD(QPD) (open diamond) results. The figure shows a major change in the reaction profile due to the influence of the solvent. The gas-phase calculations give a potential curve with no barrier for reaction of 1 with 2,an overestimated enthalpy of the total reaction, and incorrect relative energies for the pairs 5 , 6 and 7, 8 (when compared to experiment in s o l ~ t i o n ~ ~ -The ~~). calculations of the solvation energies for the stationary points on the reaction curve show that the total AEsol for 1 and 2 decreases from -133.1 (-123.0) kcal/mol forreactants to -84.1 (-82.2) kcdmol for the transition state and -76.6 (-75.7) kcaV mol for the tetrahedral intermediate for the Qc (ePD)charge set. The absolute value of the solvation energy for the separated reactants 1 and 2 is 39.3 (37.5) kcal/mol larger than the corresponding energy of the separated products 5 and 6 in the Qc (ePD) charge set. The steady decrease in the solvation energy along the reaction path for step 1 might be explained by the increase in the size of the corresponding ions. The greater stability of the pair of products 7,8 in water relative to the pair 5 , 6 could be also largely attributed to the difference in sizes of ions 5 and 8. Calculation of the activation energy for the formation of tetrahedral intermediate 3 gives 26.0 and 17.8 kcal/mol using the Qc and QpDcharge sets, respectively. The latter value is closer to the experimental value of the enthalpy of activation, 11.9 kcaVm01.~~The calculated heat of reaction for the sum of steps 1 and 2 is -35.5 (-37.3) k c d m o l for the Qc charge sets. In general, calculation of the solvation energy with the potential-derived charge set more closely approximates the experimental results than does use of the conventional AM1 charge set. This is best shown by the activation energy, Ea.and the enthalpy of step 3 (Table 2B). The results show that, for both the activation energy, Ea,and
(ePD)
I
4
LL!c
-60'o -80.0
-7.5
-5.0
-2.5
0.0
2.5
R (C-0) --O- GAS
..._._ +.- LwpD ....0.... IJyc Figure 3. Reaction path for attack by hydroxide ion in the absence of cyclodexth (Scheme 1). Points along the reaction curve are identified by the numbered structures at the top of the figure. Energies are given AHg,(2). Gas-phase AM1 calculation, open relative to Lw,,(1) boxes; Langevin dipole/AMl potential-derived charge set calculation, open diamonds;Langevin dipole/AMl Coulomb charge set calculation, open circles.
+
relative enthalpy of reaction, M,the potential-derived charge set more closely approximates the experimental results: 17.8 k c d m o l for the AM1/LD(QPD) calculation compared to the experimental value of 11.9 kcal/moP for the activation energy of step 1 and -5.7 kcaUmol (AMl/LD(QPD)) vs -8.3 kcal/ mo153-55(experiment) for the heat of reaction of step 3. 2. Computer Modeling of Phenyl Acetate Cleavage in the Reaction with /3-Cyclodextrin. Conformations and Solvation Energies of /l-Cyclodextrin and Its Ions. The optimized structural parameters (inter- and intraglucose torsional and valence angles) of 9a, 9b, and the alkoxide ions of 9a at the 0 2 and 03 positions of residue 1,9a(02-) and 9b(02-), are given in the supplementary material. The projections of the optimized geometry of 9a and 9b are given in Figure 4. The optimized structural parameters show generally similar trends in the values of the torsion angles of the experimental and optimized geometries. The values of the valence C404C1' interglucose angles, where the prime symbol indicates an atom from the (i+l)st residue, are generally 2-3" lower for the calculated geometries compared to the neutron diffraction structure. Structure 9a, the optimized geometry of the neutron diffraction structure, shows the same kind of minor disruption of the regular pattern of glucose orientations that is observed in crystal data. This disruption occurs with residue 6, which has more than an average difference in the interglucose torsional angles that connect it to residues 5 and 7. Comparison of the optimized structural parameters to the neutron diffraction structurelgashows that C3C404Cl' for residue 6 is 140.3" for the neutron diffraction structure and 141.3' for 9a, whereas the same angle is considerably smaller for the other residues in both the neutron diffraction structure and 9a. The same type of behavior is seen for angle C404C1'C; in both structures. This can be seen for 9a in Figure 4a, which shows a slight tilting of residue 6 inside
2318 J. Phys. Chem., Vol. 99, No. 8, 1995
Luzhkov and Venanzi
7
3
TABLE 3: Gas-Phase Enthalpies of Formation and Electrostatic Solvation Energies of 9, 10, and of the Alkoxide Ions of 9”P expt‘ molec AH,,, AEs,l(Qc) AEsOl(QPD) water DMF -1657.8 -64.7 -40.5 -18.0 -30.0 -1656.2 -65.3 -40.6 -1668.5 -82.4 -76.7 -1685.7 -75.3 -70.7 -1685.0 -77.2 -72.9 -1665.7 -81.2 -75.4 -75.0 -70.0 - 1684.4 -1683.8 -74.0 -69.3 -64.6 -40.8 -1713.2 a In kcaymol. See text for a description of the molecule notation. The assignment of the 02- and 0 3 - positions is given in Figure 1. From ref 57.
5
Figure 4. (A, top) Optimized geometry of 9a. (B, bottom) Optimized geometry of 9b. Residue numbering scheme is the same as for 9a. the cavity. (A full analysis of inter- and intraglucose structural features has been given by Lipkowitz et aL21b) The optimized structures of the anions 9a(02-) and 9a(03-) differ significantly in the interglucose torsional angles. Both structures differ from 9a by about 5-10’ for C3C404Cl’ and C404Cl’C; for most of the residues. For 9b, the symmetric conformation (Figure 4b), no such distinguishing feature is noted for residue 6. Rather, the C3C404Cl’ and C40&1’C2’ angles are around 128” and - 133”, respectively, for all the residues. All the glucose rings in 9b are oriented in a similar fashion within a regular annular chain that forms the cavity, and the difference in the other angles is generally of the same order as for 9a. The introduction of the symmetry conditions makes the total energy higher, which shows that the global conformational minimum of 9 might not correspond to a highly symmetric mutual orientation of glucose rings. Extra symmetry conditions in the system seem to create additional strain, which makes the gas-phase enthalpy of formation higher. This effect has been previously noted by Lipkowitz.21d In optimizing the gas-phase geometry of P-cyclodextrin anions, the major effect is associated with the orientation of the 13 secondary hydroxyl groups. Our calculations show that the most stable conformations of 9a are obtained when there are six hydroxyl groups with the clockwise orientation and seven with the counterclockwise for the 0 2 - anion or seven hydroxyl groups with the clockwise orientation and six with the coun-
terclockwise for the 03- anion. Note that in the neutral form (Figure 3), all the hydroxyl groups are oriented in one direction. Comparison of the optimized structural parameters shows that the optimized structures of the anions 9a(02-) and 9b(03-) differs significantly from that of the neutral form. Besides the OH-H distances and orientations of the secondary hydroxyl groups, this difference is observed in the interglucose torsional angles C3C404Cl’ and C404Cl’C; for most of the residues. The corresponding structural changes are explained by strong attractive interactions in the gas phase between the anionic oxygen atom and the secondary hydroxyls. The gas-phase and solvation energies of the neutral and ionized forms of 9 are presented in Table 3. The table shows that the symmetric conformation 9b has a total energy higher than the nonsymmetric conformation 9a by 1.6 kcal/mol. This result shows that the global energy minimum of 9 may not correspond to a highly symmetric orientation of glucose rings. Our subsequent calculations of more symmetric structures of 9 give additional evidence that extra symmetry conditions on the system create additional strain which makes the gas-phase enthalpy of formation higher. For the anions, the gas-phase heats of formation have the lowest values for 9a(02-) and 9a(03-) at - 1685.7 and -1684.4 kcaYmol, respectively. However, 9b(02-) and 9b(03-) have energies within 1 kcaV mol of the above values. For both 9a and 9b, the solvation energy, AEsol(Qc),is much lower than AESo~(QPD). The latter are closer to the absolute experimental enthalpies of solvation of 9 in polar solvents (-18.0 and -30.0 kcaVmo1 in water and DMF, re~pectively).~~ For the anions of 9, the solvation energies lie in the range of -74.0 to -82.4 (-69.3 to -76.7) kcaYmol for the Qc (ePD) charges, respectively, depending on the conformation. From eq 2, the difference of solvation energies for the 9a(02-) and 9a(03-) alkoxide ions equals 1.6 (2.0) kcaY mol, corresponding to 1.2 (1.5) pK, units, for the Qc (ePD) charge sets, respectively. These data agree with the experimental evidence for the higher acidity of the secondary 0 2 hydroxyls. 15as8 Calculation of the Inclusion Complexes of Phenyl Acetate with /?-Cyclodextrin and Its Anions. Comparison of the interglucose torsion values for 9a and 10 in the list of optimized structural parameters (supplementary material) shows that the benzene ring brings about essentially no disruption in the geometry of the host and is easily accommodated inside the cavity. The optimized structure of the neutral complex 10 is given in Figure 5. For the alkoxide ions 11, the calculated geometries show that in all cases, the benzene ring changes its orientation very little inside the cavity. For each of the 14 ion complexes 11, the stable orientations of the secondary hydroxyl groups with respect to the ionic center correspond to the pattems
J. Phys. Chem., Vol. 99, No. 8, 1995 2319
Computer Modeling of Phenyl Acetate Hydrolysis
TABLE 4: Gas-Phase Enthalpies“$and Solvation Enereie9-Jof 11
1
1 2 3 4
-3.3 -4.9 -1.6
5
-0.1
6 7
-0.4 -4.2
0.0
0.0
0.0
1.5
3.3 2.1 6.4 1.3 7.4 5.4
1.6
7.1 0.6 6.3 5.0
1.1 -2.8
-3.3 -4.9 0.9 -5.5
-4.6
2.5
6.1 5.1 4.9 7.0 6.4 7.6
3.3 6.3 6.4 5.4
6.4 7.5 7.5
In kcdmol. The values of Amg,and AAE,,I are given relative to the 02-alkoxide ion of residue 1 in 11. The absolute values of these quantities for this ion are AHg, = -1743.3, MsOl(Qc) = -74.6, and hEsol(QpD)= -70.1 kcdmol, respectively. The assignment of the residue numbers for 11 is given in Figure 1.
Figure 5. Optimized geometry of the neutral inclusion complex 10. H
N
Figure 6. Schematic diagram of the adducts at the 0 2 and 0 3 positions of 13. found for the anions 9a(02-) and 9b(03-). The major structural changes in the ions 11 with respect to the neutral complex 10 occur with the interglucose torsional angles and the torsion angle
C6-C1-07-C8 of the substrate, where the atom numbering scheme is given in structure 3 of Scheme 1. The structures of the most stable gas-phase complexes 11(02-) and 11(03-) are characterized by close contact distances between the oxygen anions ( 0 2 or 0 3 ) and one of the hydrogens of the methyl group of the substrate. Comparison of the enthalpies of formation of 9a, 10 (Table 3), and 1 (Table 2 ) predicts that the inclusion complex 10 has a small stabilization energy in the gas phase of 0.1 kcal/mol with respect to the substrate 1 and the host macrocycle 9a. In solution calculations, 10 is destabilized with respect to 9a and 1 by 26.2 and 15.5 kcaymol, respectively, for the Qc and QPD charge sets. Experimentally, one might expect to obtain here a negative value in the range of several k c a l / m ~ l . One ~ ~ reason for such a discrepancy between the experiment and calculations could be associated with a very small stabilization of the complex 10 in the gas phase, as calculated by the Ah41 method. The other reason may be related to a large overestimation of the solvation energy of the neutral cyclodextrin, 9. The possible explanation of the latter result is that the number of solvent dipoles taken into account by the program may be insufficient to adequately represent several solvation shells for a molecule as large as 9. However, the major focus of this work is the calculation of the energy changes in the course of ionic reactions. Previously, similar LD calculations for large ionic systems such as proteins48and hydrophobic ions49were shown to give results comparable in accuracy to molecular dynamics simulations using all-atom solvent models. The gas-phase and solution enthalpies of the 0 2 - and 03ions of 11 are given in Table 4. The numbering of the glucose residues in 11 corresponds to that of the neutral complex 10 given in Figure 5 . The table shows that, in the gas phase, the most stable 0 2 - alkoxide ion is formed at residue 3, whereas the most stable 0 3 - alkoxide ion is formed at residue 6 . In solution, the most stable 0 2 - and 0 3 - alkoxide ions are formed at residue 1. This result is independent of the atomic charge set used. Although the gas-phase calculations give 11 to be the most stable for alkoxide ions formed at positions other than residue 1, Table 4 shows that for those ions, the increase in gas-phase stability is counterbalanced by the decrease in solvation energy. For 11, the large stability of anions formed at residue 1 compared to anions formed at the other residues might be explained by the fact that, in the former case, the dipole moment of the substrate is aligned in parallel with the total dipole moment of the secondary hydroxyl groups, thus giving the largest solvation energy. Since for 11, the 0 2 - and 03anions formed at residue 1 turn out to be most stable for a given orientation of the substrate (taken from the structure of lo), these two anions were used as the “reference state” of 11 in Figure
2320 J. Phys. Chem., Vol. 99, No. 8, 1995 7a and B, respectively. Correspondingly, the tetrahedral adduct 13 and the transition-state structure 12 were also calculated for the 0 2 - and 0 3 - anions at residue 1. Table S shows that, independent of charge set, the complex formed between the substrate and the cyclodextrin with the anion at the 0 2 position of residue 1 is more stable than that with the anion at the 0 3 position of residue 1. In addition, the energy difference between these complexes is larger than that for the corresponding anions of 9 (Table 3). Calculation of the Tetrahedral Adducts of Phenyl Acetate with the Alkoxide Ions of /3-Cyclodextrin. Selected structural parameters and the gas-phase enthalpies and solvation energies of the tetrahedral adduct 13 are given in Table 5. A schematic diagram of the two possible types of adducts involving the formation of a covalent bond between atom cs8of the substrate and 0 2 - or 0 3 - of the cyclodextrin is given in Figure 6. The structure of the tetrahedral adduct 13 is determined first by the orientation of the newly-formed O-Cs8 bond with respect to the bonds of the corresponding carbon atom of the glucose ring and second by the orientation of the bonds at the tertiary carbon atom c8 (C8-07, CS-09, and Cg-Cio bonds). For both adducts in the most stable orientation, the O-Cs8 bond is pointed inside the cavity with the trans position of the H-C2-02-CS8 and H-C3-03-Cs8 atoms. For the 0 2 - adduct, there is only one direction for the C8-07 bond pointing inside the cavity. For this structure, there are two possible stereoisomers of the tertiary carbon, cs8, that differ by the orientation of Osg- and CS1oin 13. These are listed as 13a(02) and 13b(02) in Table SA. The table shows that the 13a(02) isomer is more stable than the 13b(02) isomer both in the gas phase and in solution. For the adduct at the 0 3 - ionic center, there are two orientations of the bonds at atom Cg that point in the direction of the cavity interior and one bond orientation that is pointed strictly outside the cavity. For the former two cases, we have conformations 13a(03), 13b(03), and 13c(o3) in Table 5. Conformations 13a(03) and 13b(03) have approximately equivalent space orientations of the bonds of atom cs8but differ by the angle of rotation of the benzene ring of the substrate. In the 13a(03) conformation, the plane of the benzene ring is almost parallel to the plane of the cavity, while in 13b(03) they are perpendicular. The equilibrium orientations of the chemical bonds in both adducts largely correspond to the cross-orientation of the substituents on the tertiary carbon c8 and the lone pairs and C-0 bond of atom 07. The exact structures of the adducts are determined by the steric interactions between the atoms of the substrate and the atoms of the glucose residue that are involved in the formation of the C'8-O bond. In the optimized structures of the 0 2 and 03 tetrahedral adducts, the aromatic moiety of the substrate is moved out of the cavity. This indicates that step 3 proceeds via a substantial displacement of the substrate out of its equilibrium position in the inclusion complex 10. Analysis of the structures of the tetrahedral adducts 13 shows that the substrate should be exposed to solvent to a greater extent in the case of an attack at the 0 3 center than at the 0 2 center. In the latter case, the substrate partially remains inside the cavity. This difference in the structure of the adducts is due to the fact that the C2-02 bond is more pointed into the cavity than the C3-03 bond. Calculation of the Reaction Path of Ester C - 0 Bond Cleavage in Cyclodextrin Complexes. The study of the potential energy profile for step 3 was performed using the structures 13a(O2) and 13b(03) as models of the 0 2 and 0 3 tetrahedral adducts, respectively. For reaction at the 0 3 ionic site, the 13b(03) (rather than 13a(03) or 13c(o3)) conformation was selected since it had a perpendicular orientation of the
Luzhkov and Venanzi 11
!
-10.0
12
13
Y 1
-4.0
-3.5
1
-3.0
I
-2.5
1
I
I
-2.0
-1.5
-1.0
R (C-0)
GAS
........0........ LD/PD ....0.... LDIC 11 20.0
12
13
1
10.0
3
$a "Y,
0.0
6
8 z W
-10.0
B
-
-20.0 I
-4.0
-3.5
I
-3.0
1
-2.5
I
-2.0
I
-1.5
1
-1.0
R (C-0) GAS
...-.,.0........ LD/PD
....0....
LD/C
Figure 7. Reaction path for attack by a secondary hydroxyl of cyclodextrin (step 3, Scheme 2). Points along the reaction curve are identified by the numbered structures at the top of the figure. Energies are given relative to A M g a s ( l l ) . Symbols are the same as in Figure 3. (A, top) Reaction path for attack by the 2'-hydroxyl of cyclodextrin. Uses data for 13a(Oz) (Table 5). (B, bottom) Reaction path for attack by the 3'-hydroxyl of cyclodextrin. Uses data for 13a(03) (Table 5).
J. Phys. Chem., Vol. 99, No. 8, 1995 2321
Computer Modeling of Phenyl Acetate Hydrolysis
TABLE 5: Energies and Enthalpies for Scheme 2
1
5A. Gas-Phase Heats of Formation' and Electrostatic Solvation Energies' of 13 compd
msodQ')
12(02) 12(03) 13a(Od 13b(02) 13a(03) 13b(O3) 13c(Oa)
-68.5 -61.8 -73.0 -67.9 -65.3 -65.3 -66.9
-1719.9 -1739.7 -1752.5 -1745.1 - 1760.3 - 1759.0 - 1760.0
msodQPD) -65.3 -59.8 -67.6 -64.3 -61.6 -61.4 -63.0
5B. Activation Energiesalband EnthalpiesaVcof Step 3 in Scheme 2
reaction 13a(Oz) 13b(02) 13a(03) 13b(03) 13C(03)
Ah41
AMl/LD(QC)
Ah4l/LD(pD)
29.6 -7.6
28.3 -6.7
Ea AIP Ea
23.4 -9.2
AIP
-1.8
4.9
4.0
Ea AIP Ea @ Ea
-18.2 2.5 -16.9
- 10.2 13.9 -8.9
-11.9 10.6 - 10.4
-19.9
-13.5
-15.0
AIP
5C. Selected Geometrical Parametersd 13(oz) CZ802-
a b
CZCl -125.1 -90.1
OS9CS802c2
os7Cs802c2
-37.2 135.5
82.8 20.4
c210cs8- CS10h7cs802 oZc2
-168.7 -95.0
-122.3 -179.2
cs2cs1os7Cs8
-42.8 -21.8
13(03) cs803-
a b c
C3cZ -126.6 -126.4 -102.4
os9cs803c3
-52.6 -49.6 36.7
OS7Cs803c3
64.2 69.1 -79.9
cs10cs8- cs10s7CS803 03c3
174.9 178.4 169.4
-157.7 -145.6 164.8
cs2cs1os7Cs8
-11.9 48.3 150.2
"In kcal/mol. bActivation energies, Ea, were calculated as the difference in energy between 12 and 11 using the data in 5A. The AM1 column is from the AHga data. The AMl/LD(Qc) and AM1/ LD(QPD)columns are from eq 2 using the AHg, data and either the A E s o l ( ~or ) the AEso~(QPD) data. Enthalpies of reaction, AEP, were calculated as the difference between the enthalpies of the product (13) and the reactant (11) in a similar fashion to E, above. Superscript , s indicates an atom of the substrate. Selected torsional angles (in degrees) indicate the main differences in the conformations of the tetrahedral intermediates at the 0 2 and 0 3 positions. 'I
9,
benzene ring with respect to the cavity. The points with the highest gas-phase energy along the reaction curves for stretching the Csg-02 and Csg-03 bonds are denoted in Table 5A as 12(02) and l2(O3), respectively. The corresponding potential curves for the C - 0 bond dissociation in the 13a(02) and 13b(O~)conformations are given in Figure 7A,B. The potential curves show that the reaction at the O3 center (Figure 7B) is more thermodynamically and kinetically favorable than at the 0 2 center. The potential barrier obtained for the gas-phase reaction at the 0 3 center is 2.5 kcdmol, which is lower than the calculated barrier of over 20 k c d m o l for reaction at the 0 2 center. The effect of solvent does not significantly change the difference between the heat of reaction and the potential barrier for the two reaction channels. The large difference between the reactivity at the 0 2 and 0 3 sites could be attributed to (1) unfavorable steric interactions with the cavity atoms and (2) adoption of an unfavorable substrate conformation along the reaction curve in the reaction with 0 2 - . Thus, despite the lower ~ ~despite ~ , ~the~ fact that the pKa value of the 0 2 ~ i t e and substrate more snugly fits into the host cavity for the 0 2
3 6
Figure 8. Optimized geometry of the complex of P-cyclodextrin with m-(tert-buty1)phenylacetate. tetrahedral intermediate, a more kinetically and thermodynamically favorable path is obtained for a reaction at another site-the 0 3 atom. , The barrier for step 3 of Scheme 2 in solution turns out to be higher than in the gas phase for the reaction at both the 0 2 and 0 3 centers (Figure 7). As in the case of the reaction of 1 with 2, this is connected with an unfavorable change in the solvation energy of the reactants as the reaction proceeds from the inclusion complex 11 via the transition state 12 to the tetrahedral intermediate 13. But for reaction with cyclodextrin, the extent of this desolvation is much less than for the deacylation reaction in water. The balance between the higher gas-phase reaction barrier and lower desolvation energy gives the overall potential barrier for phenyl acetate hydrolysis by cyclodextrin to be lower (10.6-13.9 kcaVmol for the reaction with the 0 3 hydroxyl, Table 5 ) than in the reaction in water (17.8-26.0 kcdmol, Table 2). Despite a relatively rigid spatial orientation of the substrate inside the macrocyclic cavity in 11, the modeling of the reaction coordinate for step 3 is not a straightforward procedure. One might expect to obtain the best results by using the molecular dynamics simulations with an umbrella sampling approach.60 In that method, the system would be forced to move along the reaction coordinate (C-0 distance) with complete structural relaxation and statistical averaging over the configuration space. The dynamical structural relaxation is especially important in the case of cyclodextrin studies, where the structure is characterized by a broad set of close minima on a potential surface. This type of approach will be considered in future work. Calculation of the Structure of the m-(fert-Buty1)phenyl Acetate Inclusion Complex with /3-Cyclodextrin. In addition to the study of phenyl acetate hydrolysis, we have calculated the structure of the cyclodextrin inclusion complex with m-(tertbuty1)phenyl acetate, 15. The hydrolysis of 15 by /3-cyclodextrin is known to be accelerated by 250 times with respect to the The starting geometry of the corresponding hydrolysis of l.15a inclusion complex was taken from the geometry of complex 10. The optimized structure of the inclusion complex of 9 with 15, i.e., 16, is given in Figure 8. In 16, the benzene ring adopts almost the same position with respect to the macrocycle as in 10. But the orientation of the substituents on the benzene ring in 16 is largely determined by the steric interactions by the bulky tert-butyl group and the primary face of the cyclodextrin cavity. This puts the acetate group in 16 much closer to the secondary hydroxyls of cyclodextrin than in 10. The distance between Csg and the closest secondary hydroxyl is 5.2 8, in 10 but only
2322 J. Phys. Chem., Vol. 99, No. 8, 1995 4.5 A in 16. This kind of preliminary alignment of the reactive groups in 16 might render the gas-phase contribution to the reaction barrier lower than in the case of complex 10 with the unsubstituted phenyl acetate. Conclusions The calculation of the reaction paths of nucleophilic attack at the acyl carbon atom of phenyl ester by hydroxide ion in the gas phase and in water demonstrates the crucial role of solvation in the determination of the activation barrier. The barrier for reaction in water is composed of the negative activation energy of the gas-phase reaction (-23.0 kcdmol) and the large positive solvation contribution (49.0 and 40.8 kcal/mol for the Qc and QPDcharge sets, respectively.) The same reaction carried out in the complex of the substrate with cyclodextrin is expected to proceed in a manner different from the two previous cases. Here, the macrocyclic cavity provides a complicated type of microheterogeneous environment for the reaction center. The reaction proceeds on the border of an aqueous medium and a nonpolar hydrocarbon host molecule. As a gross approximation, the structure of the macrocycle might be represented as a cylindrically-shaped cone with a dielectric constant of 2.0, as was done in model calculations of the decarboxylation reaction in /3-cyclodextrin.61 This kind of macroscopic treatment is not expected to be appropriate in the case of hydrolysis by cyclodextrin since here the macrocycle is directly involved in the reaction as one of the reactants. Our study provides a rigorous microscopic quantum chemical treatment of the effect of cyclodextrin on this reaction. The calculated barrier for hydrolysis by cyclodextrin at the 0 3 center consists of a small positive contribution for the gas-phase reaction (2.5 kcavmol) and a positive solvation contribution (11.5 and 8.1 kcal/mol for the Qc and QPDcharge sets, respectively). The major problem in the correct reproduction of the energetics of the reaction in cyclodextrin turns out to be associated with the calculation of the gas-phase reaction path for this conformationally flexible system. Our quantum mechanical calculations correctly reproduce a lower barrier for hydrolysis in cyclodextrin relative to the aqueous medium and provide a qualitative understanding of the reaction acceleration in cyclodextrin catalysis. Further improvement in the accuracy of the treatment of this reaction might be obtained by using molecular dynamics simulation coupled with an adequate quantum mechanical presentation of the electronic structure of the system in the transition state. Acknowledgment. This work was funded by a grant to C.A.V. from the National Academy of Sciences Cooperation in Applied Science and Technology Program. Supplementary Material Available: Table showing the optimized structural parameters of 9a, 9b, and the alkoxide ions of 9a at the 0 2 and 0 3 positions of residue 1, 9a(02-) and 9b(02-) (3 pages). Ordering information is given on any current masthead page. References and Notes (1) Computer Modeling of Carbohydrates; French, A. D., Brady, J. W., Eds.; ACS Symposium Series 430; American Chemical Society: Washington, DC, 1990. (2) Bender, M. L.; Komiyama, M. Cyclodextrin Chemistry; Springer-Verlag: Berlin, 1978. (3) Breslow, R. Adv. Chem. Ser. 1980, 191, 1. (4) Bender, M. L.; Bergeron, R. J.; Komiyama, M. The Bioorganic Chemist? of Enzymatic Catalysis; John Wiley & Sons: New York, 1984; pp 265-305.
Luzhkov and Venanzi ( 5 ) (a) Saenger, W. In Inclusion Compounds; Atwood, J. L., Davies, J. E. D., MacNicol, D. D., Eds.; Academic: New York, 1984; Vol. 2, pp 231-260. (b) Saenger, W. Angew. Chem., Znt. Ed. Engl. 1980, 19, 344. (6) Szejtli, J. Cyclodextrin Technology;Kluwer Academic: Dodrecht, 1988. (7) (a) Inoue, Y.; Hoshi, H.; Sakurai, M.; Chojjb, R. J . Am. Chem. SOC. 1985, 107, 2319. (b) Inoue, Y.; Kuan, F.-H.; ChOj8, R. Bull. Chem. SOC. Jpn. 1987, 60, 2539. (c) Inoue, I.; Okuda, T.; Miyata, Y.; Chfij8, R. Carbohydrate Res. 1984, 125, 65. (d) Inoue, Y.; Okuda, T.; ChOjb, R. Carbohydrate Res. 1983, 116, C5. (e) Inoue, Y.; KatBna, Y.; Chfij6, R. Bull. Chem. SOC.Jpn. 1979, 52, 1692. (8) Eftink, M. R.; Andy, M. L.; Perlmutter, H. D.; Kristol, D. S. J . Am. Chem. SOC.1989, 111, 6765. (9) Hamilton, J. A.; Chen, L. J. Am. Chem. SOC. 1988, 110, 5833. (10) (a) Takahashi, S.; Suzuki, E.; Nagashima, N. Bull. Chem. SOC.Jpn. 1986,59, 1129. (b) Takahashi, S.; Suzuki, E.; Amino, Y.; Nagashima, N.; Nishimura, Y.; Tsuboi, M. Bull. Chem. SOC. Jpn. 1986, 59, 93. (1 1) Maheswaran, M. M.; Divakar, S. Indian J . Chem. 1991, 30A, 30. (12) Wood, D. J.; Hruska, F. E.; Saenger, W. J . Am. Chem. SOC.1977, 99, 1735. (13) Kano, K.; Mori, K.; Uno, B.; Goto, M.; Kubota, T. J . Am. Chem. SOC. 1990, 112, 8645. (14) Lukovits, I. THEOCHEM 1988, 170, 249. (15) (a) VanEtten, R. L.; Sebastian, J. F.; Clowes, G. A,; Bender, M. L. J . Am. Chem. SOC.1967,89, 3242; (b) Ibid. 1967,89, 3253. (c) VanderJagt, D. L.; Killian, F. L.; Bender, M. L. J . Am. Chem. SOC.1970, 92, 1016. (d) Straub, T. S.; Bender, M. L. J. Am. Chem. SOC. 1972,92, 8875. (e) Bender, M. L. J . Inclusion Phenom. 1984, 2, 433. (f) D’Souza, V. T.; Bender, M. L. Arc. Chem. Res. 1987, 20, 146. (16) (a) Breslow, R.; Czamiecki, M. F.; Emert, J.; Hamaguchi, H. J . Am. Chem. SOC. 1980, 102, 762. (b) Trainor. G. L.; Breslow, R. J . Am. Chem. SOC. 1981, 103, 154. (c) Breslow, R.; Trainor, G.; Ueno, A. J . Am. Chem. SOC. 1983, 105, 2739. (17) Rama, Rao, K.; Srinivasan, T. N.; Bhanumathi, N.; Sattur, P. B. J . Chem. Soc., Chem. Commun. 1990, 10. (18) (a) Breslow, R.; Anslyn, E.; Huang, D.-L. Tetrahedron 1991, 47, 2365. (b) Tong, W.; Ye. H.; Rong, D.; D’Souza, V. T. J . Comput. Chem. 1992, 13, 614. (19) (a) Betzel, C.; Saenger, W.; Hingerty, B. E.; Brown, G. M. J . Am. Chem. SOC. 1984, 106, 7545. (b) Zabel, V.; Saenger, W.; Mason, S. A. J. Am. Chem. SOC.1986, 108, 3664. (c) Steiner, T.; Mason, S. A,; Saenger, W. J . Am. Chem. SOC. 1991, 113, 5676. (20) (a) Koehler, J. E. H.; Saenger, W.; Van Gusteren, W. F. Eur. Biophys. J . 1987, 15, 197. (b) Koehler, J. E. H.; Saenger, W.; Van Gunsteren, W. F. Eur. Biophys. J . 1987, 15, 211. (c) Koehler, J. E. H.; Saenger, W.; Van Gunsteren, W. F. Eur. Biophys. J . 1988, 16, 153. (d) Koehler, J. E. H. In Molecular Dynamics: Applications in Molecular Biology; Goodfellow, J. M. D., Ed.; CRC: Boca Raton, FL, 1990; pp 69106. (e) Koehler, J. E. H.; Saenger, W.; van Gunsteren, W. F. J. Mol. Biol. 1988, 203, 241. (f) Koehler, J. E. H.; Saenger, W.; van Gunsteren, W. J . Biomolec. Struct. Dynam. 1988, 6, 182. (g) Linert, W.; Margl, P.; Renz, F. Chem. Phys. 1992, 161, 327. (21) (a) Lipkowitz, K. B.; Raghothama, S.; Yang, J. J . Am. Chem. SOC. 1992, 114, 1554. (b) Lipkowitz, K. B.; Green, K.; Yang, J. A. Chirality 1992, 4, 205. (c) Lipkowitz, K. B.; Green, K. M.; Yang, J. A,; Pearl, G.; Peterson, M. A. Chirality 1993, 5, 51. (d) Lipkowitz, K. B. J . Org. Chem. 1991, 56, 6357. (22) Thiem, H.-J.; Brandl, M.; Breslow, R. J . Am. Chem. SOC. 1988, 110, 8612. (23) Menger, F. M.; Sherrod, M. J. J . Am. Chem. SOC.1988,110,8606. (24) Linert, W.; Han, L.-F.; Lukovits, I. Chem. Phys. 1989, 139, 441. (25) Ohashi, M.; Kasatani, K.; Shinohara, H.; Sato, H. J . Am. Chem. SOC. 1990, 112, 5824. (26) Venanzi, C. A.; Canzius, P. M.; Zhang, 2.; Bunce, J . D. J . Comput. Chem. 1989, 10, 1038. (27) Wertz, D. A.; Shi, C. X.; Venanzi, C. A. J. Comput. Chem. 1992. 13, 41. (28) (a) Ingold, C. K. Structure and Mechanism in Organic Chemistry; Cornell University: Ithaca, NY, 1969; pp 1128-1191. (b) Kirby, A. J. In Comprehensive Chemical Kinetics; Bamford, C. H., Tipper, C. F. H., Eds.; Elsevier: Amsterdam, 1972; Vol. 10, pp 57-207. (29) (a) Burgi, H. B.; Lehn, J. M.; Wipff, G. J . Am. Chem. SOC.1974, 96, 1956. (b) Alagona, G.; Scrocco, E.; Tomasi, J. J . Am. Chem. SOC.1975, 97, 9876. (c) Miertfis, S.; Kysel, 0.;Krajcik, K. Chem. Zvesti 1981, 35, 3. (d) Williams, I.; Spangler, D.; Femec, D.; Maggiora, G.; Schowen, R. J. Am. Chem. SOC. 1983, 105, 31. (e) Dewar, M. J.; Storch, D. M. J. Chem. Soc., Chem. Commun. 1985, 94. (f) Weiner, S . J.; Singh, U. C.; Kollman, P. A. J. Am. Chem. SOC.1985, 107, 2219. (30) Madura, J. F.; Jorgensen, W. L. J. Am. Chem. SOC.1986,108,2517. (3 1) Warshel, A. Computer Modelling of Chemical Reactions in Enzymes and Solutions; John Wiley & Sons: New York, 1991. (32) (a) DeTar, D. F.; Tenpas, C. J. J . Am. Chem. SOC. 1976, 98,7903. (b) Pop, E.; Huang, M. J.; Brewster, E. M.: Bodor. N. Int. J . Quant. Chem.: Quant. Biol. Symp. 1992, 19, 77. (33) Ueno, A.; Breslow, R. Tetrahedron Lett. 1982, 23, 281.
J. Phys. Chem., Vol. 99, No. 8,1995 2323
Computer Modeling of Phenyl Acetate Hydrolysis (34) Fujita, K.; Tahara, T.; Imoto, T.; Koga, T. J . Am. Chem. SOC.1986, 108, 2030. (35) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J . Am. Chem. SOC. 1985, 107, 3902. (36) Stewart, J. J. P. MOPAC-93.00 Manual; Fujitsu Limited: Tokyo, Japan, 1993. (37) Stewart, J. J. P. J . Camp.-Aided Mol. Des. 1990, 4 , 1. (38) Dewar, M. J. S.; Storch, D. M. J . Am. Chem. SOC. 1985,107,3898. (39) Lluch, J. M.; Bertran, J.; Dannenberg, J. J. Tetrahedron 1988, 44, 7621. (40) Saba, G.; Scano, P.; Sedda, P. A,; Thomson, C. THEOCHEM 1990, 204, 379. (41) Nakayama, A.; Richards, W. D. Quant. Struct.-Act. Relat. 1987, 6, 153. (42) Marcos, E. S.; Maraver, J. J.; Chiara, J. L.; Gomez-Sanchez, A. J . Chem. SOC., Perkin Trans. 2 1988, 2059. (43) Dannenberg, J. J. J. Phys. Chem. 1988, 92, 6869. (44) Galera, S.; Lluch, J. M.; Oliva, A.; Bertran, A. THEOCHEM 1988, 40, 101. (45) Khalil, M.; Woods, R. J.; Weaver, D. F.; Smith, V. H., Jr. J. Compur. Chem. 1991, 12, 584. (46) Tupper, K. J.; Gajewski, J. J.; Counts, R. W. THEOCHEM 1991, 236, 21 1. (47) CHEM-X, developed and distributed by Chemical Design, Ltd., Oxford, England. (48) (a) Warshel, A.; Russel, S. Q.Rev. Biophys. 1984, 17, 283. (b) Lee, F. S.; Chu, Z. T.; Warshel, A. J. Comput. Chem. 1993, 14, 161. (49) Luzhkov, V.; Warshel, A. J. Comput. Chem. 1992, 13, 199.
(50) Besler. B. H.: Merz. K. M., Jr.: Kollman, P. A. J . Comout. Chem. 1960,11, 431.' (51) Buono. R. A.: Venanzi. T. J.: Zauhar, R. J.: Luzhkov, V. B.: Vena&, C. A. J. Am.'Chem. SOC. 1994, 116, 1502. (52) Schroder, S.; Daggett, V.; Kollman, P. J . Am. Chem. SOC. 1991, 113, 8922. (53) Wilson, B.; Georgiadis, R.; Bartmess, J. E. J. Am. Chem. SOC.1991, 113, 1762. (54) Pearson, R. C. J . Am. Chem. SOC. 1986, 108, 6109. (55) Christensen, J. J.; Hansen, D. L.; Izatt, R. M. Handbook ofProton Ionization Heats and Related Thermodynamic Quantities; John Wiley & Sons: New York, 1976. (56) Bunton, C. A.; O'Connor, C.; Turney, T. A. Chem. Ind. (London) 1967, 1835. (57) Danil de Namor, A. F.; Traboulssi, R.; Lewis, D. F. V. J . Am. Chem. SOC. 1990, 112, 8442. (58) Casu, B.; Reggiani, M.; Gallo, G.G.;Vigevani, A. Tetrahedron 1968, 23, 281. (59) Inoue, Y.; Liu, Y.; Tong, L.-H.; Shen, B.4.; Jin, D . 4 . J.Am. Chem. SOC. 1993, 115, 10637. (60) Beveridge, D. L.; DiCapua, F. M. In Computer Simulations o j Biomolecular Systems. Theoretical and Experimental Applications; van Gunsteren, W. F., Weiner, P. K., Eds.; Escom: Leiden, The Netherlands, 1989; pp 1-26. (61) Furuki, T.; Hosokawa, F.; Sakurai, M.; Inoue, Y.; Chiij8, R. J . Am. Chem. SOC. 1993, 115, 2903.
JP942367U
