Validating a Strategy for Molecular Dynamics Simulations of

9110

J. Phys. Chem. B 2009, 113, 9110–9122

Validating a Strategy for Molecular Dynamics Simulations of Cyclodextrin Inclusion Complexes through Single-Crystal X-ray and NMR Experimental Data: A Case Study Giuseppina Raffaini,*,† Fabio Ganazzoli,† Luciana Malpezzi,† Claudio Fuganti,† Giovanni Fronza,‡ Walter Panzeri,‡ and Andrea Mele*,† Dipartimento di Chimica, Materiali e Ingegneria Chimica “G. Natta”, Politecnico di Milano, Via L. Mancinelli 7, 20131 Milano, Italy, CNRsIstituto di Chimica del Riconoscimento Molecolare, Via L. Mancinelli 7, 20131 Milano, Italy ReceiVed: February 20, 2009; ReVised Manuscript ReceiVed: May 13, 2009

A theoretical and experimental study about the formation and structure of the inclusion complex (-)-menthylO-β-D-glucopyranoside 1 with β-cyclodextrin (β-CD) 2 is presented as paradigmatic case study to test the results of molecular dynamics (MD) simulations. The customary methodological approachsthe use of experimental geometrical parameters as restraints for MD runssis logically reversed and the calculated structures are a posteriori compared with those obtained from NMR spectroscopy in D2O solution and single crystal X-ray diffraction so as to validate the simulation procedure. The guest molecule 1 allows for a broad repertoire of intermolecular interactions (dipolar, hydrophobic, hydrogen bonds) concurring to stabilize the host-guest complex, thus providing the general applicability of the simulation procedure to cyclodextrin physical chemistry. Many starting geometries of the host-guest association were chosen, not assuming any a priori inclusion. The simulation protocol, involving energy minimization and MD runs in explicit water, yielded four possible inclusion geometries, ruling out higher-energy outer adducts. By analysis of the average energy at room temperature, the most stable geometry in solution was eventually obtained, while the kinetics of formation showed that it is also kinetically favored. The reliability of such geometry was thoroughly checked against the NOE distances via the pair distribution functions, that is, the statistical distribution of intermolecular distances among selected diagnostic atoms calculated from the MD trajectories at room temperature. An analogous procedure was adopted both with implicit solvent and in Vacuo. The most stable geometry matched that found with explicit solvent but major differences were observed in the relative stability of the metastable complexes as a consequence of the lack of hydration on the polar moiety of the guest. Finally, a control set of geometrical parameters of the thermodynamically favored complex matched the corresponding one obtained from the X-ray structure, while local conformational differences were indicative of packing effects. 1. Introduction Cyclodextrins (CD), such as cyclomaltohexaose, cyclomaltoheptaose, and cyclomaltooctaose (R-CD, β-CD, and γ-CD, respectively), are cyclic oligosaccharides with a truncated-cone shape, having an inner hydrophobic cavity and an external hydrophilic surface. The narrow and the wide rims of the macrocycle comprise primary and secondary hydroxyl groups, respectively. Because of these features, apolar organic molecules or polar molecules with an apolar region can be accommodated into the cavity of CDs to form energetically stable host-guest inclusion complexes.1 The formation of such supramolecular complexes is a well assessed phenomenon with applications covering a broad range of fields, including chiral separations,2 drug delivery and targeting,3 food chemistry,4 sensors.5 The driving forces for the formation of host-guest complexes of cyclodextrins are the energetically favorable nonbonding interactions, such as hydrogen bonds (H-bonds), hydrophobic effects, and van der Waals interactions (dispersion forces and dipolar interactions) of the guest included into the host CD * To whom correspondence should be addressed. (G.R.) Fax: 0039 02 23993180. Tel.: 0039 02 23993024. E-mail: [email protected]. (A.M) Fax: 0039 02 23993180. Tel.: 0039 02 23993006. E-mail: [email protected]. † Politecnico di Milano. ‡ CNRsIstituto di Chimica del Riconoscimento Molecolare.

cavity. Additionally, in water the release of relatively isolated small clusters of water molecules from the hydrophobic cavity entropically favors guest inclusion, largely compensating the entropy loss due to inclusion. Clearly, the extent of each contribution to the stabilization of the complex depends on the chemical nature of the guest and host molecules.6 Computational methods based on molecular mechanics (MM) and molecular dynamics (MD) have been extensively used in parallel to experimental structural methods such as NMR, calorimetry, and X-ray diffraction to obtain a comprehensive model of the complexation and molecular recognition phenomena. Typically, a first-approximation model deduced from the experimental data is refined by simple energy minimization (MM methods), or, more recently, through MD simulations in Vacuo or increasingly more often in explicit water, simulating the dynamic behavior on nanoseconds time-scale. Whenever available, the atom coordinates obtained from single crystal X-ray diffraction structure of a given complex are used as starting geometry for MD simulations in water aimed at simulating or predicting the structure of the complex in solution or in an environment mimicking physiological conditions. The availability of suitable force fields and the decreased costs of computing prompted many investigators to massive use of simulations, thus raising the issue of control protocols for the reliability of calculated geometry and conformational properties.

10.1021/jp901581e CCC: $40.75  2009 American Chemical Society Published on Web 06/15/2009

Validating MD Simulations with NMR and X-ray Data

J. Phys. Chem. B, Vol. 113, No. 27, 2009 9111

SCHEME 1: Molecular Structures and Atom Numbering Scheme for Compounds 1 and 2a

a Hydrogen atoms of 1 and 2 are identified by using lower case letters and atom numbers as shown: g ) glucose of menthyl glucoside; m ) menthyl residue; c ) glucose units of cyclodextrin (e.g., g3 refers to the H atom in position 3 of the glucose unit of compound 1, while c3 indicates the homologous proton of compound 2). Symmetry equivalent glucose units of 2 are labelled A-G according to Breslow.15

In the present work a paradigmatic case study is reported and discussed where the results of MM and MD simulations carried out in Vacuo, with implicit solvent and in explicit water on a specifically tailored β-cyclodextrin inclusion complex are validated against solution NMR data and compared with solidstate X-ray results. To this end, the customary approach mentioned abovesfrom experiments to calculationssis logically reversed, going from calculations to experiments. For this study, the guest molecule, (-)-menthyl-O-β-D-glucopyranoside 1 (Scheme 1) was encapsulated into β-cyclodextrin 2 (β-CD, see Scheme 1), and this process was theoretically modeled through computer simulations. The inclusion complex of 1 and 2 was characterized in aqueous solution by NMR spectroscopy and in the solid state by single crystal X-ray diffraction, and the experimental data were thoroughly compared with the simulation results for validating the simulation procedure. The model guest molecule 1, was chosen on the basis of two criteria: (i) Structural suitability for the formation of an inclusion complex with CD exhibiting the broadest spectrum of interactions (hydrophobic, hydrogen bond, dipolar); indeed, 1 fulfills the requirements proposed by Saenger for the ideal guest, namely the presence, in the molecular frame, of “hydrophobic body and hydrophilic “ends” which can form hydrogen bonds to the hydroxyl groups on both sides of the host molecule to stabilize the inclusion complex.”.7 (ii) The possibility of achieving a thorough structural characterization of the inclusion complex both in solution and in the solid state. The quality of the structure determination in solution, in turn, depends primarily on the number of intermolecular NOEs that can be collected, providing a sufficiently wide set of host-guest distances for the formulation of a first-approximation geometry of the complex. As for the solid-state structure, it should be considered that most cyclodextrin inclusion complexes do not afford single crystals suitable for X-ray diffraction studies. Therefore, the choice of the present model system was the result of a large number of attempts of crystallizing inclusion complexes of β-CD with small-medium sized glycoconjugates.8 The plan of the paper is the following: we first report the experimental and computational details, then we discuss the results of the NMR experiments in solution and of the solidstate structure obtained by X-ray single-crystal analysis, and afterward we report the simulation results in solution (by MD runs at room temperature) and in vacuo (by MD runs and energy minimizations) to validate the simulation procedure in comparison with the experimental data.

2. Materials, Experimental Methods, and Computational Details Cyclomaltoheptaose (β-cyclodextrin, 2) was purchased from Sigma-Aldrich and used without any further purification. (-)Menthyl-O-β-D-glucopyranoside (1) was obtained from San Giorgio Flavours (Torino, Italy) and used as such. Deuterium oxide (99.97%) was obtained from Cortec (France). 2.1. Molecular Modeling and Data Analysis. The simulations were performed with InsightII/Discover 2000,9 using the consistent valence force field CVFF.10 This force field describes nonbonded interactions through van der Waals and Coulombic terms only, with no extra term for H-bonds. CVFF, originally designed to model proteins, was later augmented to include additional functional groups, accounting also inter alia for the acetal moiety related with the anomeric effect in carbohydrates.11 Therefore, CVFF can be satisfactorily used for these molecules,12 even though this effect is not conformationally dominant in CDs because of the geometric constraint imposed by the macrocycle. Extensive tests carried out in comparison with NMR data do support the general accuracy of CVFF for oligosaccharides.13 The geometries of β-CD and of 1 were generated with the available templates of InsightII, then subjected to an MD run in vacuo at room temperature and finally optimized up to an energy gradient lower than 4 × 10-3 kJ mol-1 Å-1. 1 was then placed close to β-CD but outside the hydrophobic cavity, selecting 12 unbiased trial geometries (Figure 1) for the initial optimizations. The arrangements were chosen so that all the main different sides of 1 could approach the two rims and the outer surface of β-CD. No inclusion complex was assumed in this stage, unlike what it is often done.14 The whole system was then hydrated by adding more than one thousand water molecules at the density of 1 g cm-3 in a cubic cell with edges of 33 Å and periodic boundary conditions (constant-volume conditions), and after full geometry optimization all the resulting adducts were subjected to independent MD simulations. Further simulations were carried out for the same starting geometries in vacuo and with implicit solvent in an effective dielectric medium to check whether these much faster procedures produce satisfactory results. The dynamic equations were integrated using the Verlet algorithm with a time step of 1 fs at a temperature of 300 K, controlled through the Berendsen thermostat, and the instantaneous coordinates were periodically saved for further analysis. The system equilibration was monitored by the time change of the potential energy and of its components, and of

9112

J. Phys. Chem. B, Vol. 113, No. 27, 2009

Raffaini et al. r from atoms i (or from a specific point such as their center of mass), and is defined here in non-normalized form as

gij(r) ) d〈Nij(r)〉/dV(r)

(2)

where d〈Nij(r)〉 is the average number of times the j atoms are comprised in a spherical shell of thickness dr at a distance r from atom i within an MD run. Here gij(r) is not normalized, because we consider only finite systems (the host and the guest), so that gij(r) f 0 for r f ∞, and upon volume integration gij(r) yields the total number of the i,j atom pairs, 〈Nij〉:

〈Nij〉 ) 4π

Figure 1. The starting geometries adopted to simulate the formation of the host-guest complex of (1) with β-CD close to the secondary CD rim (first column), to the primary rim (second column) and to the outer surface (third column). The oxygen atoms are shown in darker gray, and the hydrogen atoms were omitted for clarity.

relevant intra- and intermolecular distances, such as those between opposite pairs of glycosidic or hemiacetalic oxygen atoms, and those between the centers of mass of each ring of 1 and of the β-CD macrocycle. A further check of the system equilibration was also carried out by the similarity maps of the host, of the guest and of the whole complex (see next paragraph). On the basis of these equilibration criteria, the MD runs were carried out for 0.5 ns in explicit water and for 2 ns in Vacuo and with implicit solvent. In all cases, the main changes only took place in the initial part of the runs. Note that the system thermalization was significantly faster in water than in the other cases due not only to the thermostat, but also to the random collisions with the solvent molecules. Simulation Data Analysis: Similarity Maps (rmsd) and Pair Distribution Function (PDF). For an analysis of the unlike conformations found in the MD runs, the similarity maps, or rmsd (i.e., root-mean-square distance) maps, are constructed by calculating for a given set of n instantaneous conformations, or frames, the n × n root-mean-square distances among sets of atoms (for instance the guest and host molecule), and plotting them as a function of the frame indices as a bidimensional map with an appropriate color coding. The rmsd between two frames is given by the minimum of the function

{

N

∑

1 (ra - rib)2 N i)1 i

}

1/2

(1)

where the index i runs over the N selected atoms, the a, b superscripts indicate two different frames, and r is the atomic vector position, while the minimum is chosen to remove trivial effects due to rigid-body translations or rotations of the whole system. Thus, similar conformations belonging to the same family of conformers have a small rmsd, and unlike conformations have a large rmsd. The equilibrium geometries sampled during an MD run, including the distances between selected sets of atoms, are best analyzed through the pair distribution function gij(r) (or simply PDF) calculated from the trajectories of the simulations. The PDF gives the probability density of finding atoms j at a distance

∫0∞ r2gij(r) dr

(3)

2.2. NMR Spectroscopy. The samples for NMR spectroscopy were prepared by dissolving a suitable amount of crystallized 1:1 complex of 1 and 2 in 750 µL of D2O in order to obtain 3 mM solutions. The NMR spectra were carried out on an Avance 500 spectrometer (Bruker, Karlsruhe, Germany) at 305.0 ( 0.1 K. All chemical shifts were referenced to internal methanol (3.310 ( 0.001 ppm) according to Funasaki et al.16 Phase-sensitive two-dimensional rotating frame nuclear Overhauser enhancement correlation experiments (2D-ROESY) were carried out by using the pulse sequence proposed by Desvaux et al.17 for enhancing sensitivity and reducing artifacts due to TOCSY type magnetization transfer. Data were collected over 2K points in F2, 512 increments zero filled to 1K in the processing stage. Spin lock radiofrequency power and offset were calibrated for θ ) 40°, with θ ) arctan [(γB1/2π)/∆ν] and ∆ν ) frequency difference (Hz) between hard pulse and spin lock. Spin lock duration was 300 ms. Control experiments in identical conditions were also carried out on 1 alone in order to spot intramolecular interactions involving protons of the menthyl and the glucopyranose ring of 1 and to differentiate them from the intermolecular ones observed between the menthyl group of 1 and the glucopyranose rings of 2. The stoichiometry of the complex was assessed via the continuous variation method (Job plot):18 stock solutions of 1 (solution A) and 2 (solution B) in D2O, 15 mM each, were prepared and used for chemical shift titration. The NMR solutions were prepared by mixing increasing volumes VB of solution B and keeping constant VA + VB ) 1 mL for each sample. The complexation induced chemical shift variation of c3 of 2 was monitored for the assessment of the stoichiometry of the complex. All the experiments were repeated in triplicate for the sake of reproducibility. 2.3. X-ray Crystallography. The inclusion complex of 1 with 2 was crystallized by slow evaporation from an aqueous solution of the components. Suitable crystals for X-ray diffraction analysis were obtained as very thin plates. The selected crystal, 0.6 × 0.3 × 0.04 mm in size, was glued to a glass fiber. The 1:1 complex crystallizes in the orthorhombic space group P212121 (No. 19), a ) 15.342(3), b ) 21.983(3), c ) 22.960(4) Å, Z ) 4, Dc )1.383 g cm-3, µ ) 1.07 mm-1, F(000) ) 3421, room temperature. The X-ray diffraction data were collected on a Siemens P4 diffractometer with graphitemonochromated Cu KR radiation (1.54179 Å), with the ϑ/2ϑ scan technique. Unit cell parameters were determined at room temperature by using 56 reflections in the range 16 e 2ϑ e 50°. A total of 8606 reflections (7940 unique, Rint ) 0.0316) were collected up to 136° in 2ϑ; 5957 of them showed I g 2σ(I) and were considered. The stability of the crystal was


J. Phys. Chem. B, Vol. 113, No. 27, 2009 9113 The full list of data collection, structure determination and refinement is reported in the Supporting Information, Table S1.

Figure 2. Continuous variation plot (Job plot) for the formation of the inclusion complex of 1 and 2. The quantity ∆δ is the complexation induced chemical shift variation of proton c3. [H] and [G] indicate the molar concentration of host 2 and guest 1, respectively.

periodically monitored during data collection and no crystal decay was observed. The structure was solved by direct methods (SIR97) providing almost all the non-H atoms of the β-CD complex. Subsequent difference Fourier syntheses and leastsquares refinement (SHELX-97 program)19 converged to a final R(F2) ) 0.0550 for the observed reflections and wR(F2) ) 0.1580 for all data. The hydrogen atoms were generated in idealized position by the XP module of SHELX-97 and refined using a riding model. Five fully positioned water molecules were found on the electron density map, while minor peaks of electron density suggested partially occupied water sites: only those with occupation factor refined to values g0.20 were retained in the final refinement. The H atoms of the water molecules were not introduced. Because of the large number of parameters, the final cycle of refinement was carried out in two blocks: one for β-CD alone (702 refined parameters) and one for guest and water molecules (366 refined parameters). The maximum values for positive and negative residual electron density were 0.31 and -0.30 e Å-3, respectively. The goodness of fit was S ) 1.040.

3. Results and Discussion 3.1. The Solution Structure of the Complex: NMR. The 1 H NMR spectrum of a mixture of 1 and 2 in D2O showed the typical complexation-induced chemical shift variations of selected protons: deshielding of the signals of the guest molecule protons and shielding of c3 and c5 of cyclodextrin (see Scheme 1). These findings provided the first evidence of the formation of a true inclusion complex in solution.20 The continuous variation plot (Job plot) reported in Figure 2 shows a sharp maximum at r ) [H]/([H] + [G]) ) 0.5, consistent with 1:1 stoichiometry of the inclusion complex. A thorough characterization of the solution geometry of the complex was obtained by the analysis of the intermolecular nuclear Overhauser enhancements in the rotating frame (ROEs). An expansion of the 2D-ROESY experiment on the inclusion complex is shown in Figure 3. The one-dimensional projection on the F1 dimension is related to signals assigned to the protons of the menthyl moiety, while the F2 projection displays the spectral region of the nonanomeric protons of CD. Intense cross-peaks relating protons c3 of cyclodextrin and m6eq, m3eq, and m2 of the menthyl residue of 1 indicate that the latter protons are in spatial proximity to the secondary OH rim of cyclodextrin. Interestingly, m6eq, m3eq, and m2 do not give rise to cross peaks with c5. In a complementary way, the cross peaks corresponding to the pairs c5/m4eq and c5/m5ax point out that the molecular frame consisting of the carbon atoms in position 4 and 5 (including also the methyl group 7) of the menthyl moiety is close in space to the primary OH rim of cyclodextrin. Once again, m4eq and m5ax do not show dipolar contact with c3. Furthermore and quite surprisingly, also proton m8 shows no ROE with either c3 or c5. These findings are consistent with the menthyl residue of 1

Figure 3. Detail of the contour plot of 2D-ROESY of the inclusion complex. Horizontal (F2) and vertical (F1) projection are referred to the menthyl moiety of 1 and to the nonanomeric protons of 2, respectively. The assignment of relevant signals is given according to the atom labeling used in Scheme 1. Cross-peaks arising from intermolecular ROEs are framed in rectangular boxes (for ROEs of methyl groups see Figure 4). Significant intramolecular ROEs of proton m1 at 3.63 ppm (m1/m6eq, m1/m5ax, m1/m3ax, and m1/m10) are also marked with an asterisk.

9114


Raffaini et al.

Figure 4. Detail of the contour plot of 2D-ROESY of the inclusion complex. The cross-peaks related to the methyl groups of 1 and the inner protons of 2 are showed. The assignment of relevant signals is given according to the atom labeling used in Scheme 1. Cross-peaks arising from intermolecular ROEs are framed in rectangular boxes. No contribution of intramolecular ROEs of 1 is detectable.

entering the hydrophobic cavity of 2 from the secondary OH rim and fully inserted into the cavity. The intermolecular ROEs due to the methyl groups of 1 with the inner protons c3 and c5 of 2 are displayed in Figure 4. The observed pattern is fully consistent with the complexation geometry inferred from the data reported above. The two diastereotopic methyl groups m9 and m10 of the isopropyl substituent in position 2 of the menthyl residue show strong correlations with c3, and no correlation at all with c5. The isopropyl substituent is therefore close to the large rim of the CD cavity. Conversely, the methyl group in position 7, namely far away from the glycosidic junction, shows selective ROE with c5 and no correlation with c3. More interestingly, a weak but significant ROE between m7 and both c6 protons is clearly detectable, suggesting that such methyl group is protruding out of the cavity of CD and close in space to the CH2OH tails of CD. The discussion so far assumed that inter- and intramolecular ROEs could be unambiguously distinguished. Any possible misinterpretations due to overlap of inter- and intramolecular ROEs was double checked by control experiments carried out on the guest 1 alone. ROESY experiments on 1 showed that m1 and g1, falling in the glucopyranose region of the spectrum (3.2-4.5 ppm), are the only protons giving rise to strong intramolecular ROEs with the menthyl protons falling in the high field region of spectrum (0.7-2.2 ppm). The signals of m1 (3.63 ppm) and g1 (4.49 ppm) are well separated from those of c3 (3.91 ppm) and c5 (3.74 ppm), and therefore the intramolecular cross-peaks do not overlap the intermolecular ones. The general solution geometry of the complex can be summarized as follows: (i) formation of a 1:1 inclusion complex;

(ii) the menthyl group of 1 is located inside the cavity of 2; (iii) the glycosidic linkage of 1 is on the secondary OH rim of CD; (iv) the glucose moiety is outside the cavity, exposed to water solvent and likely to establish hydrogen bonds with OH groups of CD. The latter conclusion can only be indirectly inferred from points i-iii as no direct ROEs between hydrogen atoms of the glycan of 1 with H atoms of 2 could be observed due to spectral overlap. 3.2. The Solid State Structure of the Complex: X-ray Diffraction. β-CD and (-)-menthyl-O-β-D-glucopyranoside form an host-guest inclusion complex in the 1:1 ratio. The solid-state structure of the complex is shown in Figure 5, while some geometrical parameters of β-CD are reported in Table 1, in comparison with analogous quantities obtained by the computer simulations (see later). The β-CD macrocyclic ring, characterized by the usual truncated-cone shape, shows a pseudo 7-fold symmetry, apart from some tilt of the sugar ring, as indicated by the geometry of the almost planar heptagon defined by the O4 atoms (see Table 1): the sides of the heptagon range from 4.315(6) to 4.512(6) Å (mean value 4.372 Å) and O4(n - 1) · · · O4(n) · · · O4(n + 1) angles range from 123.7(1) to 132.4(2)° (mean value 128.2° to be compared with the ideal value of 128.6° for regular heptagon). The perimeter of the heptagon defined by the O4 atoms is 30.60 Å. For the sake of comparison, the values of the geometrical parameters reported above are 4.36 ( 0.12 Å and 128.3 ( 3.7° for uncomplexed β-CD.21 Furthermore, the O4 atoms show relatively small deviations from planarity, as reported in Table 1 through the distances from their mean plane, the largest deviations being +0.240(3) Å and -0.282(3) Å.


Figure 5. The structure of the inclusion complex viewed along two directions in order to illustrate the geometry of inclusion and the intermolecular hydrogen bonds. For the sake of clarity, all H atoms not involved in hydrogen bonds have been omitted.

A good descriptor of the conical form of the β-CD molecule is the average value of the dihedral angles between the O4 atoms plane and the least-squares plane of the C2, C3, C5, and O5 atoms of each glucose residue: as shown in Table 1, these values range from 67.9(1)° to 87.2(2)°, with a mean value of 78.9(2)°. The degree of inclination of each residue relative to the molecular axis is well measured by the tilt angle, defined as the dihedral angles between the O4(n) plane and the planes defined by the O4(n + 1)-C1(n) · · · C4(n)-O(4)n atoms; it ranges from 7.0° to 29.5°. The cyclic rigid shape of the β-CD molecule is mainly stabilized, as expected, by the intramolecular H-bond between the O2-H group of a glucopyranoside residue and the O3-H group of the adjacent one, (see Table 2). Within each glucosyl residue, bond angles and distances show no significant differences. The conformation of the glucosyl residues, all in the 4C1 conformation, is clarified by values reported in Table 1. The most relevant differences are in the orientation of the C6-O6 bonds. In fact, apart the O6G-H group, all the primary hydroxyl groups point outside the cavity and show the preferred gauche-gauche conformation with mean torsion angles C4-C5-C6-O6 and O5-C5-C6-O6 of 56.2(6) and -65.0(6)°, respectively, as shown in Table 1. The C6G-O6G bond is rotated inward and shows a gauche orientation with respect to the O5G-C5G bond and a trans orientation with respect to the C4G-C5G bond. The detailed values of the torsion angles involving the orientation of the C6-O6 bonds are affected in general by crystal packing interactions (see later).

J. Phys. Chem. B, Vol. 113, No. 27, 2009 9115 Accordingly, they show significant deviations from the theoretical values obtained for the isolated complex by the simulations in vacuo, a feature that shall be discussed in some detail in a later section. In the complex, the guest is partially located inside the β-CD cavity showing the menthyl residue fully inserted into the hydrophobic cavity (Figure 5). The glucose residue of the guest 1 protrudes out of the cavity and lays in the intermolecular space in proximity of the secondary OH rim of cyclodextrin. The stability of the inclusion complex is well ensured by two different factors: the hydrophobic interaction of the menthyl moiety with CD cavity, and the host-guest hydrogen bonds. In the present complex, an unusually high number of hydrogen bonds has been found. Indeed, each host is connected to its guest by one O-H · · · O and two C-H · · · O short interactions, as indicated in Figure 5. In the crystal assembly of cyclodextrins inclusion complexes, the hydrogen bond network stabilizes both the complex and the crystal packing. We have demonstrated8a that H-bonds between host (H ) and guest (G ), either within the complex (Hi-Gi ) or across two adjacent complexes (Hi-Gj ), are specific stabilizing factors of the complex. Similarly, H-bonds between two separate host molecules close in the crystalline assembly (Hi-Hj ) or between guest molecules close in space but belonging to distinct complexes (Gi-Gj ), mainly contribute to the stability of the crystal packing. Further consolidation of the packing structure is also achieved by a number of H-bonded water molecules. The crystal packing of the title complex is shown in Figure 6 and the list of the observed H-bonds and their classification in the above-mentioned types is summarized in Table 2. A sketch of the position of a given complex with the closest in-plane neighbors is also reported in Scheme 2. As can be seen from Table 2, each inclusion complex interacts with other 10 neighboring complexes. In the crystal, the β-CD complexes are arranged as monomeric entities in an herringbone fashion, assuming a cage-type structure and forming infinite chains parallel to the a axis. In the intermolecular space there are 10 water molecules per unit cell; 5 of them are affected by disorder and are distributed over 13 positions. No water molecules are located within the macrocycle cavity. An extended hydrogen bonding cluster among the water molecules and the primary and secondary hydroxyl groups of the β-CD contributes to the stabilization of the crystal packing. 3.3. Computer Simulations. Computer simulations show that in vacuo isolated β-CD has a truncated conical shape with little distortion from C7 symmetry, while minor distortions are present in water with a slight tilt of a few glucoside rings which does not greatly affect the symmetry of the cavity, unlike what happens for R-, γ-, or δ-CD.22 Compound 1 has an elongated shape, and the MD runs in vacuo show that at room temperature it has a limited flexibility, localized at the torsion angles around the bonds adjacent to the glycosidic oxygen. These angles show large fluctuations, with full relative rotations of the two ring systems. In the optimized geometry, 1 can approach β-CD with its different sides at either rim, or at the outer surface of β-CD as shown in Figure 1. No host-guest complex was assumed a priori, and guest inclusion was eventually obtained only through MD runs. All the geometries shown in Figure 1 were used as starting points for the initial energy minimizations in explicit water, in vacuo, and with implicit solvent (see the Experimental Section). In water, these minimizations produced high-energy outer adducts, with roughly the same arrangement as in Figure 1. The

9116


Raffaini et al.

TABLE 1: Some Relevant Geometrical Parameters for the Host-Guest Complex of 1 with β-CD Obtained by the Single Crystal X-ray Analysis (Bold-Face Values) and by the Computer Simulations for the Most Stable Arrangement Obtained in Vacuo, Very Similar to IA in Figure 7 (Values in Italic) glucosyl residue of β-CD distances (Å) O4n-O4(n + 1) a

anglesb (deg) O4n-O4(n+1)-O4(n + 2) distancesc (Å) O4n-plane Id dihedral angles (deg) between planes I and IId,e dihedral (or tilt) angles (deg) between planes I and IIId,f torsion angles (deg) O5n-C5n-C6n-O6n torsion angles (deg) C4n-C5n-C6n-O6n

A

B

C

D

E

F

G

4.410 4.585 123.7 129.1 0.240 -0.064 77.6 87.6 12.9 2.9 -76.5 172.5 45.7 -65.5

4.314 4.539 130.2 130.4 0.017 0.123 83.2 91.4 7.1 3.2 -71.2 169.3 48.4 -68.5

4.520 4.473 130.4 128.3 -0.082 0.011 81.4 83.9 9.0 6.4 -52.0 171.3 68.0 -66.6

4.348 4.507 126.4 125.9 -0.112 -0.136 83.3 96.5 14.4 4.5 -65.9 169.6 57.0 -68.4

4.321 4.605 125.8 130.2 0.215 0.061 67.9 86.0 23.2 7.6 -62.7 168.0 58.1 -70.1

4.370 4.479 132.5 130.8 0.004 0.097 87.2 93.0 7.0 4.3 -61.8 172.7 60.4 -65.5

4.328 4.480 128.6 124.9 -0.282 -0.092 71.8 98.0 29.5 6.8 59.4 65.6 178.5 -171.8

a Average values with standard deviations: 4.37 ( 0.07 Å from the X-ray results and 4.52 ( 0.05 Å from the computer simulations. Average values with standard deviations: 128.2 ( 3.1° from the X-ray results and 128.5 ( 2.3° from the computer simulations. c Standard deviations of ( 0.003 Å from the X-ray results. d Plane I is the average plane defined by the O4 atoms of β-CD. e Plane II is the plane defined by atoms C2, C3, C5, and O5. f Plane III is the plane defined by atoms O4′, C1, C4, O4.

b

TABLE 2: List of Hydrogen Bonds Involving a Complex Moiety with All Its First Neighboring Symmetry-Related Complexes (See Text for a Discussion of Types Interactions and Scheme 2 for a Graphical View). The Symbols H and G Indicate the Host and Guest Molecules, Respectively D-H · · · A

D-H (Å)

H · · · A (Å)

D · · · A (Å)

∠(DHA) (deg)

code for interaction

O2A-H2OA · · · O3B O2A-H2OA · · · O4B O3B-H3OB · · · O2A O3B-H3OB · · · O4B O3C-H3OC · · · O2B O2D-H2OD · · · O3D O3D-H3OD · · · O2C O3E-H3OE · · · O2D O2G-H2OG · · · O3A O2G-H2OG · · · O4A

0.75 0.75 0.70 0.70 0.85 0.85 0.85 0.82 0.89 0.89

Type of Interaction: Host-Host (H-H ), Intramolecular 2.13 2.866(5) 165.1 2.32 2.732(5) 115.8 2.17 2.866(5) 177.1 2.50 2.835(5) 111.6 1.90 2.742(5) 174.4 2.47 2.899(5) 112.5 2.04 2.885(6) 174.3 1.97 2.790(5) 173.9 2.03 2.906(5) 168.0 2.30 2.753(4) 111.4

H 1-H 1 H 1-H 1 H 1-H 1 H 1-H 1 H 1-H 1 H 1-H 1 H 1-H 1 H 1-H 1 H 1-H 1 H 1-H 1

O3A-H3OA · · · O2 C3-H3C · · · O2A C5-H5C · · · O3B

Type of Interaction: Host-Guest (H-G ) within the Same Complex 0.85 1.94 2.775(4) 166.5 1.14 2.58 3.688(5) 163.9 0.94 2.57 3.478(5) 162.0

H 1-G 1 G 1-H 1 G 1-H 1

O6C-H6OC · · · O3a O2D-H2OD · · · O4b O6E-H6OE · · · O2c C1E-H1CE · · · O3c O2-H2O · · · O5Ed O2-H2O · · · O6Ed O3-H3O · · · O6Fd O4-H4O · · · O2De O6-H6O · · · O3Ee

Type of Interaction: Host-Guest (H-G ) of Adjacent Complexes 0.88 1.81 2.692(4) 171.4 0.85 2.49 2.823(4) 104.2 0.95 2.40 3.340(5) 170.8 1.07 2.35 3.297(5) 146.4 1.01 2.03 2.951(3) 150.7 1.01 2.63 3.340(5) 127.2 0.85 1.84 2.671(4) 163.7 0.82 2.01 2.823(4) 171.6 0.85 1.78 2.635(4) 179.6

H 1-G4 H 1-G 5 H 1-G 7 H 1-G 7 G 1-H 8 G 1-H 8 G 1-H 8 G 1-H 10 G 1-H 10

O6A-H6OA · · · O3Gf O6B-H6OB · · · O6Fg O6B-H6OB · · · O5Fg O6D-H6OD · · · O6Ah C1C-H1CC · · · O6Gh O3G-H3OG · · · O6Ed O6G-H6OG · · · O6Di C1F-H1CF · · · O6Aj

0.63 0.96 0.96 0.85 0.97 0.95 0.85 1.10

Type of Interaction: Host-Host (H-H ), Intermolecular 2.08 2.676(5) 156.5 2.23 3.131(6) 156 2.42 3.104(6) 128 1.99 2.839(5) 179.5 2.55 3.481(7) 161.0 1.71 2.606(5) 155.1 1.95 2.799(6) 179.5 2.31 3.399(6) 167.6

H 1-H 2 H 1-H 3 H 1-H 3 H 1-H 6 H 1-H 6 H 1-H 8 H 1-H 9 H 1-H 11

Symmetry codes: a 1/2 - x, -y, 1/2 + z. b 1/2 + x, -1/2 - y, 1 - z. c 1/2 - x, -y, 1/2 + z. d 1/2 -x, -y, z - 1/2. e x - 1/2, -y - 1/2, 1 - z. f x - 1/2, 1/2 - y, 1 - z. g x - 1, y, z. h -x, y - 1/2, 1/2 - z + 1. i -x, 1/2 + y, 1/2 -z + 1. j 1/2 + x, 1/2 - y, 1 - z.

MD runs at 300 K of these adducts allowed the systems to largely explore the configurational space, and to optimize the geometry of the inclusion complex in some local, deep energy minimum. To find the most stable arrangement at room

temperature, some MD runs were repeated two or three times with unlike random initial velocities to ensure a better sampling. When 1 was initially placed close to the outer surface of β-CD, it could simply “diffuse” on the outer surface, or fly away,


Figure 6. A view of a two-dimensional section of the crystal packing, projected along the b axis. Hydrogen bonds are drawn as dashed lines. For the sake of clarity, all water molecules have been omitted.

SCHEME 2: Sketch of the Position of a Given Complex in the Crystal Structure with the Closest in-Plane Neighbors. H ) host, G ) guest. The Numbers Refer to an Arbitrary Numbering of 1:1 Complexes

J. Phys. Chem. B, Vol. 113, No. 27, 2009 9117 the thermalization of the system through random collisions, slowing down however the relaxation because of the solvent “viscosity”. The lowest-energy complex shows inclusion of the menthyl group in the hydrophobic cavity from the secondary rim: minor distortions are present, and the glucoside ring is well hydrated by the solvent (not shown in Figure 7) and forms an H-bond with the secondary hydroxyls. However, two different arrangements are still possible, due to the possible variation of the glycosidic torsion angles of 1. In water, the most stable arrangement shows the isopropyl group of the menthyl moiety close to the β-CD secondary rim (IA in Figure 7), while the methyl group buried in the cavity approaches the primary rim, as also found experimentally in solution and in the solid state. A further geometry (IB in Figure 7) shows a rotated menthyl moiety, so that the isopropyl and the methyl groups swap their position. In water, at room temperature such arrangement is less stable by about 12.5 kJ/mol than the former one: this value favors the most stable geometry by a factor of about 150, as obtained through the Boltzmann statistical weights. The arrangement IA is robust, and the depth of the guest into the cavity remains almost constant through the major part of the MD run (see also later). Conversely, in higher-energy geometries the fluctuations are much larger, almost leading to the guest expulsion from the cavity. For instance, in some MD runs the menthyl group entered the cavity with an unfavorable orientation, and its center of mass showed large fluctuations with respect to the center of mass of the macrocycle in the distance range 2 j d j 5 Å, suggesting a dynamic equilibrium between inner and outer positions with weak attractive interactions. Similar large fluctuations were also observed within the MD runs for the inclusion complex formed in water by benzyl alcohol in β-CD, where the small size of the guest compared to the cavity ensured a significant conformational freedom.23 A few snapshots of the inclusion kinetics of 1 in β-CD in water to the most stable arrangement IA are reported in Figure 8. In the figure, we display the six water molecules initially clustered inside the cavity (the MD simulations indicate an average value of 6.3 in isolated β-CD22), while the bulk solvent is not shown. These molecules are eventually set free upon guest inclusion, increasing the system entropy and favoring the complex formation. The kinetics of inclusion is summarized in Figure 9, where the potential energy Epot and the distance d between the centers of mass of the menthyl group and of the macrocycle are plotted as a function of time. In either case, as well as in all the other MD runs, the inclusion kinetics follows a simple exponential decay described by

y ) y∞ + A exp(-t/τ)

maximizing the overall configurational entropy. Otherwise, the MD runs stabilized the system forming a small number of inclusion complexes with a few unlike arrangements. The final optimizations eventually led to four basic inclusion types, labeled as I to IV, with two minor variants indicated with the letters A or B, for a total of six geometries shown in Figure 7. We note that the kinetics of inclusion is strongly affected by the solvent: the water molecules act as a thermal bath and hasten

(4)

where y, y∞, can be Epot, Epot,∞, or d, d∞, the ∞ subscript denotes the equilibrium values achieved at very long time, and τ is the corresponding characteristic time (τE or τd). The relative stability of the various geometries in water at 300 K was established through the (average) values Epot,∞. The best-fit parameters are reported in the caption of Figure 9. The figure shows that the apolar moiety quickly reaches its final position in the cavity, with a characteristic time τd ) 11.5 ( 0.4 ps. On the other hand, since τE ) 25.2 ( 1.4 ps, Epot settles to a constant value at a somewhat longer time. We conclude that according to the MD simulations the guest quickly reaches its final, most favorable position, whereas minor, somewhat slower rearrangements of β-CD are required to optimize the host-guest interactions. Interestingly, all the other runs producing higher-

9118


Raffaini et al.

Figure 7. The geometries of inclusion obtained after the MD runs in water. The labels refer to the various inclusion arrangements as described in the text. The most stable inclusion complex obtained in water corresponds to arrangement IA. In both panels, the guest is shown in a ball-and-stick representation, and the β-CD by a line drawing for clarity. The oxygen atoms are in a darker gray color, and the hydrogen atoms were omitted for simplicity, while the water molecules are not shown for clarity.

Figure 8. Snapshots of the inclusion of the guest 1 (shown with balls and sticks) in the β-CD cavity obtained in the MD trajectory in water and leading to the most stable geometry IA. The concomitant expulsion of the water molecules initially present in the hydrophobic cavity (shown as balls and sticks for clarity) is also shown. The bulk water outside the β-CD cavity is not shown for clarity. The simulation times of the snapshots are: (a) 0 ps; (b) 1 ps; (c) 2.5 ps; (d) 5 ps; (e) 10 ps; (f) 12.5 ps; (g) 17.5 ps; (h) 25 ps; (i) 40 ps.

energy systems did always show characteristic times in the range 22-25 ps both for τE and for τd, indicating that the thermodynamically most stable state is also kinetically favored. Other metastable inclusion complexes are possible. One of them (complex II in Figure 7) involves an incomplete inclusion of the glucose ring of 1 in the cavity of 2 through the secondary

rim, with an (average) potential energy larger by 25.1 kJ/mol than the most-stable state, and a distance between the centers of mass of the sugar ring of 1 and of the macrocycle that is about twice as large. Thus, the glucoside moiety of 1 is effectively close to the secondary hydroxyls of β-CD, while the menthyl group in water suffers a poor hydration, due to the


Figure 9. The time evolution of the potential energy Epot (upper plot) and of the distance d between the centers of mass of the macrocycle and of the menthyl group of compound (1) (lower plot), calculated from the MD runs in water starting from the uppermost arrangement in the left column of Figure 1. The continuous gray lines are the bestfit curve obtained from eq 4 with the fitted values Epot,∞ ) -40480 ( 7 kJ/mol, AE ) 1523 ( 59 kJ/mol, and τE ) 25.2 ( 1.4 ps for Epot, and d∞ ) 0.87 ( 0.01 Å, Ad ) 6.0 ( 0.1 Å, and τd ) 11.5 ( 0.4 ps for d.

lack of H-bonds. Interestingly, this inclusion geometry was also obtained starting from the same initial arrangement (uppermost arrangement in the left column of Figure 1) yielding also the most stable final inclusion complex in a separate MD run. In the present case, the inclusion process, driven by a different set of initial velocities, showed at first a fleeting, shallow inclusion of the menthyl group, followed by its expulsion, an overall rotation of 1 with respect to β-CD, and a two-step inclusion of the sugar ring. Such kinetics is shown in Figure 10 through the time dependence of the distances between the centers of mass of the two ring systems of 1 and of the macrocycle, and through the similarity map (rmsd map, see section 2.1) of the whole system. The latter plot clearly shows the distance between the conformations in all the frame pairs saved in the MD trajectory, allowing detection of the different families of “conformers” clustered close to the main diagonal of the rmsd map with a small rms distance. Such families effectively correspond to the initial arrangement, the fleeting inclusion and then the expulsion of the menthyl group, the rotation of 1, and the final, two-step inclusion of the sugar ring in a robust metastable geometry. Further metastable inclusion geometries may involve menthyl or glucoside insertion of the guest through the primary rim (geometries III and IV in Figure 7). This process is however less favorable than the previous ones, due in part to the smaller size of the opening. Actually, the size difference does not affect the kinetics of inclusion, which shows again characteristic times τE in the 22-25 ps range, but it does affect the complex energy due to the macrocycle strain: thus, insertion of the menthyl group leads to a complex less stable than the best one by 21.0 kJ/ mol, while a larger relative destabilization (amounting to 75.3 kJ/mol) is found upon glucoside inclusion. Moreover, energy minimizations of instantaneous geometries taken during the MD trajectories often lead to a partial expulsion of the included

J. Phys. Chem. B, Vol. 113, No. 27, 2009 9119

Figure 10. (i) Upper plot: the distance between the centers of mass of the ring systems of 1 (indicated in the legend) and of the β-CD macrocycle plotted as a function of time. (ii) Lower plot: the rmsd map of the whole system calculated from the same MD run: a frame was saved every 0.5 ps, and the pixel color shows the distance (shown in the left scale in Å) between the frames plotted along the two axes. Darker pixels correspond to similar conformations with a small rms distance [as defined in eq 1], white pixels correspond to highly different conformations separated by an rms distance lager than 2 Å.

glucoside, favored by H-bonds and dipolar interactions with the hydroxyl groups at the β-CD rim, and by hydration. In conclusion, the simulations in explicit water show that the arrangement IA is thermodynamically the most favorable one, and that its formation is also kinetically favored. A further MD run in water at room temperature lasting for additional 500 ps was then carried out for this arrangement to analyze the distribution of the intermolecular distances among diagnostic hydrogen atoms for a comparison with the results of the NOE experiments discussed before. The usual threshold of 4 Å for a vanishing NOE is considered, since the intermolecular crossrelaxation occurs within the long-lived inclusion complex (vide supra).24 The distribution of the intermolecular distances related to NOE contacts is best analyzed through the PDF plots (see section 2.1), yielding the probability density of finding selected sets of atoms as a function of their separation. In this way, thermal fluctuations and molecular motions at equilibrium are taken into account in the calculation of statistical averages to be compared with the experimental data. Selected PDF obtained for geometry IA are reported in Figure 11, showing results for the diagnostic β-CD hydrogens c3 and c5 (see Scheme 1 for the nomenclature), and for the main hydrogens of the menthyl moiety of compound 1. Figures 11a and 11b report the PDF between hydrogens m9 and m10 and hydrogens c3 and c5. The first peak at r = 2.4 Å for the distances among c3 and m9, m10 shows the close proximity of these atoms, indicating that the isopropyl methyl groups are close to the secondary rim, as also found through the NOE data. Conversely, the much larger distance of protons c5 from the same atoms (their PDF is nonvanishing only for r > 4 Å) as well as the large separation of m8 from c3 and c5 shown in Figure 11c (m8 is closer to c3 than to c5, but the distances are anyway larger than 4 Å) are again consistent with the lack of NOE among these atoms. Also, Figure 11d shows the closer proximity of hydrogens m7 to c5 than to c3, again in keeping with the observed intense NOE. The same panel also shows

9120


Raffaini et al.

Figure 11. The PDF of selected sets of diagnostic atoms (see text) shown in the legends of each panel, plotted as a function of their separation r (in Å). These plots were calculated from the MD runs in explicit water after full equilibration of the inclusion complex with the most stable geometry IA shown in Figure 7.

the broad distribution of distances between m7 and c6: the low value of the PDF at small r values (2 Å < r < 4 Å) entails a small probability of achieving a short distance between these sets of atoms, hence a very weak NOE. The close separation between the menthyl axial proton m5 and proton c5 should also be noted: the corresponding PDF (see Figure 11e) is much larger than that involving protons m5 and c3, also in agreement with the NOE pattern. We further report in Figure 11f the PDF between the equatorial proton m4 and protons c3 and c5, showing the closer proximity of m4 with the latter β-CD atom than with the former one, a pattern opposite to that shown by m2 (see Figure 11g), once more in agreement with the NOE results. Other PDF plots for additional host-guest proton pairs are reported in the Supporting Information (Figure S1).

We should mention here that the formation of an inclusion complex formed by benzyl alcohol and β-CD was already studied through MD runs in water starting from an outer arranngement of the guest.23 In this case, the main interest was on the conformational properties of the host after guest inclusion, for instance at the glycosidic bonds of the macrocycle, including also the distribution of the tilt angles. The inclusion geometries were clustered in families of conformers based on the depth of inclusion, with minimal energy diferrences, however, which suggest large fluctuations favored by the small size of the guest. The analysis of the inclusion geometry was relatively limited, since the interest of the paper was mainly devoted to alternative methodologies to sample the configurational space of the system, and therefore, for instance, the kinetics of the process was largely

Validating MD Simulations with NMR and X-ray Data ignored. Therefore, no meaninglful comparison can be made with our metastable geometries or our calculated distribution of distances. We now briefly comment on the simulations carried out in vacuo and with implicit solvent starting again with the arrangements of Figure 1. In both cases, the initial energy minimizations already produced an inclusion complex, or sometimes a loosely bound outer adduct. After the MD runs, the latter adducts could also form an inclusion complex. Interestingly, with implicit solvent if the guest was initially placed close to the side surface of β-CD no inclusion took eventually place. This pattern is the same as in water, but it differs from what found in vacuo, where inclusion could also be obtained from this initial arrangement. The MD runs carried out with implicit solvent yielded the same inclusion geometries and relative stability as in explicit water. In fact, the most stable geometry is the same as type IA reported before, which in turn fully agrees with the experimental results, both in solution and in the solid state. The main difference in the host-guest arrangement compared to the simulations in water consists in a slightly deeper inclusion within the cavity: the average distance between the centers of mass of the menthyl ring and of the macrocycle amounts to 0.314 ( 0.002 Å, less than in water, where the H-bonds formed with the glucose ring favor a slightly shallower inclusion. The relative stability of the higher-energy metastable states, calculated by optimization of many instantaneous conformations of the MD trajectory, is again the same as in water. In implicit solvent, the inclusion complex having the menthyl group entering the cavity through the primary rim is less stable by 9.6 kJ/mol, whereas the complexes with an included glucoside ring have a still higher energy, amounting to 10.9 and 13.0 kJ/mol above the lowest energy minimum if inclusion takes place from the secondary or the primary rim, respectively. These energy differences should be compared with those calculated in water, where the included glucose ring yields a smaller stability than with implicit solvent due to lack of H-bonds with water. The simulations in vacuo, after the MD runs and the geometry optimizations of many instantaneous geometries, yielded the same arrangements discussed before. The most stable inclusion complex basically shows the same geometry as the most stable one found in explicit water. On the other hand, noticeable differences are found for the relative stabilities of the metastable inclusion complexes. In fact, and quite surprisingly, inclusion of the glucose ring from the primary rim appears to be quite stable in vacuo, its energy being only larger by 1.3 kJ/mol than the most stable one. Inclusion of the glucose ring from the secondary rim yields an energy that is larger than the most stable one by 5.9 kJ/mol, and inclusion of the menthyl group from the primary rim produces a complex with a still higher energy by 7.5 kJ/mol. This implies that in vacuo the relative populations of these states are in the order 1:0.61:0.09:0.05, so that no unique state would really dominate. The relatively large stability of the host-guest complexes with an included glucose moiety is clearly related to the lack of hydration for the whole system, and with the strong H-bonds and dipolar interactions that the sugar can show at the two CD rims with a shallow inclusion. This result confirms also that the simulations in vacuo should be treated with great care when a guest with a polar moiety is included in water. The most stable geometry obtained in vacuo can be compared with the solid-state structure, if crystallization is under thermodynamic control, as it is likely in the present case (slow evaporation of the solvent). This geometry is basically the same as found in solution, and compares reasonably well with the

J. Phys. Chem. B, Vol. 113, No. 27, 2009 9121 solid-state structure determined through the X-ray crystal analysis reported in section 3.2. Differences between the two geometries are indicative of crystal packing effects. Some geometrical parameters of the most stable geometry calculated in vacuo are reported in Table 1 in comparison with the corresponding X-ray values. The crystal packing effects are most evident in the details of the host conformation concerning in particular the torsion angles of the primary hydroxyl groups (last two lines in Table 1). In fact, the simulations in vacuo indicate that these hydroxyls mainly form intramolecular H-bonds, whereas in the solid state structure they form intermolecular H-bonds involving both the neighboring complexes and the water molecules trapped in the crystal, leading to large differences in the torsion angles (see Table 1) between calculated and observed values. For the same reason, the distances between the hydroxyl groups of adjacent rings within the host are significantly different, and also the tilt angles show smaller values and a narrower distribution in the simulation results in vacuo compared to the X-ray data (see Table 1) to optimize the intramolecular H-bond at both rims. Interestingly, a difference is also found for the torsion angle around the C5-C6 bond of the glucose moiety of the guest molecule: such difference is again due to an intermolecular H-bond of the primary hydroxyl in the crystal structure, which becomes a host-guest H-bond involving a secondary hydroxyl of β-CD in the simulation geometry. 4. Conclusions This paper reports a theoretical and experimental study of a host-guest inclusion complex formed by a menthyl glycoconjugate with a macrocyclic oligosaccharide, β-cyclodextrin (βCD). The theoretical study was carried out with computer simulations in explicit water and did not assume any prior guest inclusion within the host molecule, but rather considered the possible outer arrangements, and inclusion was only obtained through molecular dynamics (MD) runs at room temperature. The possible inclusion geometries were ranked on the basis of the average energy at equilibrium achieved within the MD runs. The most favorable inclusion complex obtained in water by this systematic ab initio procedure were fully validated by comparison of the distribution of the intermolecular distances with those inferred from the experimental NOE data in solution. Moreover, the analogous simulations performed in vacuo yielded a similar arrangement, with a good agreement with the solidstate structure shown by single-crystal X-ray analysis. Local conformational differences with the latter geometry provided a clear indication of the intermolecular interactions due to crystal packing, most notably through H-bonds with neighboring complexes and/or water molecules. The simulations carried out in this paper also suggest that a computationally efficient strategy can be adopted to establish the preferred host-guest geometry and its stability compared to other possible inclusion arrangements. Such strategy involves simulations with implicit solvent through energy minimizations and MD runs, and comparison of the resulting arrangements ranking them in terms of relative stability through their potential energy. However, the present results suggest that additional, shorter runs in explicit water should be carried out to check for the relative stabilities and obtain the final energies with greater accuracy. In this way, hydration effects due to the solvent, and in particular H-bonds with the water molecules, would be easily determined. For similar reasons, the in vacuo simulations should be treated and interpreted with great care. Nevertheless, the most stable geometry eventually obtained in vacuo essentially yields

9122


the solid-state structure, and we suggest that it can provide a benchmark for assessing specific conformational effects due to the intermolecular packing interactions. Supporting Information Available: Figure S1 with further plots of the Pair Distribution Functions for selected atom pairs calculated from the MD runs. Table S1 with a summary of the X-ray data collection procedure. This material is available free of charge via the Internet at at http://pubs.acs.org. References and Notes (1) (a) Saenger, W. Angew. Chem., Int. Ed. 1980, 19, 344–362. (b) Jozwiakowski, M. J.; Connors, K. A. Carbohydr. Res. 1985, 143, 51–59. (c) Szejtli, J. In ComprehensiVe Supramolecular Chemistry; Atwood, J. L., Davies, J. E. D., Macnicol, G. D., Vo¨gtle, F., Eds.; Pergamon: Oxford, 1996; Chapter 5, p 189-203. (2) Mitchell, C. R., Armstrong, D. W. In Chiral Separations; Gu¨bitz, G., Schmid, M. G., Eds.; Springer: Berlin, 2004; Chapter 3, p 61-112. (3) (a) Hirayama, F.; Uekama, K. AdV. Drug DeliVery ReV. 1999, 36, 125–141. (b) Uekama, K.; Hirayama, F.; Irie, T. Chem. ReV. 1998, 98, 2045– 2076. (4) Szejtli J. Cyclodextrin Technology; Kluwer Academic: Dordrecht, The Netherlands, 1988; p 307-328. (5) Liu, Y.; Shi, J.; Guo, D.-S. J. Org. Chem. 2007, 72, 8227–8234. (6) Rekharsky, M. V.; Inoue, Y. Chem. ReV. 1998, 98, 1875–1917. (7) An˜ibarro, M.; Geβler, K.; Usoń, I.; Sheldrick, G. M.; Saenger, W. Carbohydr. Res. 2001, 333, 251–256. (8) (a) Malpezzi, L.; Fronza, G.; Fuganti, C.; Mele, A.; Bru¨ckner, S. Carbohydr. Res. 2004, 339, 2117–2125. (b) Fronza, G.; Fuganti, C.; Genesio, E.; Mele, A. J. Incl. Phenom. Macrocyc. Chem. 2002, 44, 225–228. (9) Accelrys Inc. InsightII 2000; San Diego, CA. See also the URL: http://www.accelrys.com/.

Raffaini et al. (10) Dauber-Osguthorpe, P.; Roberts, V. A.; Osguthorpe, D. J.; Wolff, J.; Genest, M.; Hagler, A. T. Proteins: Struct., Funct., Genet. 1988, 4, 31– 47. (11) Bush, C. A.; Martin-Pastor, M.; Imberty, A. Annu. ReV. Biophys. Biomol. Struct. 1999, 28, 269–293. (12) von der Lieth, C.-W.; Kozar, T. J. Mol. Struct. (Theochem) 1996, 368, 213–222. (13) Asensio, J. L.; Martin-Pastor, M.; Jimenez-Barbero, J. Int. J. Biol. Macromol. 1995, 137–248. (14) (a) Franchi, P.; Lucarini, M.; Mezzina, E.; Pedulli, G. F. J. Am. Chem. Soc. 2004, 126, 4343–4354. (b) Dodziuk, H.; Lukin, O. Chem. Phys. Lett. 2000, 327, 18–22. (c) Dodziuk, H.; Lukin, O.; Nowin´ski, K. S. J. Mol. Struct. (Theochem) 2000, 503, 221–230. (d) Zubiaur, M.; De Federico, M.; Burusco, K. K.; Bea´, I.; Virgili, A.; Sanchez-Ferrando, F.; Jaime, C. J. Incl. Phenom. Macrocycl. Chem. 2005, 51, 241–247. (e) Nuń˜ez-Aguëro, C.-J.; Escobar-Llanos, C.-M.; Dıáz, D.; Jaime, C.; Gardunõ-Jua´rez, R. Tetrahedron 2006, 62, 4162–4172. (15) Breslow, R.; Doberty, J. B.; Guillot, G.; Hersh, C. L. J. Am. Chem. Soc. 1978, 100, 3227–3229. (16) Funasaki, N.; Nomura, M.; Yamaguchi, H.; Ishikawa, S.; Neya, S. Bull. Chem. Soc. Jpn. 2000, 73, 2727–2728. (17) Desvaux, H.; Berthault, P.; Birlirakis, N.; Goldman, M.; Piotto, M. J. Magn. Reson. A 1995, 113, 47–52. (18) Connors, K. A. Binding Constants. The Measurement of Molecular Complex Stability; John Wiley & Sons: New York, 1987, p 24-28. (19) Sheldrick, G. M.; Schneider, T. R. In Methods in Enzymology; Carter, C. W., Jr., Sweet, R. M., Eds.; Academic: San Diego, CA 1997; Vol. 277, pp 319-343. (20) Inoue, Y. Annu. Rep. NMR Spectrosc. 1993, 27, 59–101. (21) Lindner, K.; Saenger, W. Carbohydr. Res. 1982, 99, 103–115. (22) Raffaini, G.; Ganazzoli, F. Chem. Phys. 2007, 333, 128–134. (23) Varady, J.; Wu, X.; Wang, S. J. Phys. Chem. B 2002, 106, 4863– 4872. (24) Bagno, A.; Rastrelli, F.; Saielli, G. Progress NMR Spectrosc. 2005, 47, 41–93.

JP901581E

Validating a Strategy for Molecular Dynamics Simulations of

Recommend Documents