5268
2008, 112, 5268-5271 Published on Web 04/08/2008
Spatial Arrangement of r-Cyclodextrins in a Rotaxane. Insights from Free-Energy Calculations Yanmin Yu,† Wensheng Cai,*,‡ Christophe Chipot,§ Tingting Sun,‡ and Xueguang Shao‡ Institute of Green Chemistry and Fine Chemicals, Beijing UniVersity of Technology, Beijing, 100022, P.R. China, Department of Chemistry, Nankai UniVersity, Tianjin, 300071, P.R. China, and Equipe de dynamique des assemblages membranaires, UMR CNRS/UHP 7565, Nancy UniVersite´ , BP 239, 54506 VandœuVre-le` s-Nancy cedex, France ReceiVed: December 3, 2007; In Final Form: February 28, 2008
The rotaxane formed by R-cyclodextrins (R-CDs) threaded onto a poly(ethylene glycol) (PEG) chain was investigated in the gas phase and in an aqueous solution by means of molecular dynamics simulations. The free-energy difference between the three possible spatial arrangements of consecutive R-CDsviz.. head-tohead (HH), head-to-tail (HT), and tail-to-tail, was determined using free-energy perturbation calculations. These simulations reveal that the interaction of R-CD with the PEG chain is very similar in the two surroundings, whereas the mutual interaction of the macrocycles is stronger in the gas phase than in the aqueous solution. Moreover, difference in the overall interaction appears to stem primarily from changes in the electrostatic contribution. Analysis of intermolecular hydrogen bonds indicates that hydrogen bonds created between R-CD and water molecules weaken the hydrogen-bonding interaction of adjacent R-CDs. Comparison of the freeenergy differences characteristic of the three possible spatial arrangements of contiguous R-CDs reveals that the HH motif of the rotaxane is the most stable in the gas phase due to the hydrogen-bond formation between the secondary hydroxyl groups of the two R-CDs, and the slight preference for the HT motif in aqueous solution, which can be related to the directionality of the dipole moment borne by the macrocycles in these two motifs.
Rotaxanes are mechanically interlocked molecular architectures consisting of mobile cyclic molecules threaded onto linear, usually polymeric chains, dissociation of the macrocycles from the latter being hindered by bulky groups, or stoppers, located at both ends of the chain.1,2 In recent years, rotaxanes have been the object of an increased interest on account of their potential to serve as molecular devices, molecular machines, and functional materials.3-7 Among the suitable macrocycle candidates, cyclodextrins (CDs) have been used successfully as hosts to form rotaxanes with linear, polymeric chains.1,2,8 For instance, rotaxanes resulting from the assembly of R-CD with poly(ethylene glycol) (PEG) have recently attracted considerable attention,9-13 Although PEG-R-CD rotaxanes were examined by employing a host of experimental methods, spatial arrangement of the macrocycles threaded onto the wirelike polymeric chain remains in general poorly understood.14,15 To a large extent, this can be explained by the fact that experiment only provides macroscopic and averaged information about the organization into supramolecular assemblies. Investigating the structural properties of rotaxanes featuring different spatial arrangements of R-CDs constitutes the main objective of the present contribution. Since in the course of * To whom correspondence should be addressed. E-mail: wscai@ nankai.edu.cn. Fax: +86-22-23502458. † Beijing University of Technology. ‡ Nankai University. § Nancy Universite ´.
10.1021/jp711413a CCC: $40.75
rotaxane preparation, the supramolecular assembly experiences a variety of environments,11 two distinct media, viz.. the gas phase and an aqueous solution, were considered here. The basic three spatial arrangements of R-CDs forming the rotaxane, viz.. head-to-head (HH), head-to-tail (HT), and tail-to-tail (TT), were investigated in these two surroundings using a combination of molecular dynamics (MD) simulations with the free energy perturbation (FEP) method.16 Detailed analysis of interaction energies and hydrogen-bond formation illuminates key structural differences in the three spatial arrangements examined. These differences are quantified by means of FEP calculations.16 The initial Cartesian coordinates of R-CD were taken from previously published crystal structures.17 The PEG chain featuring two stoppers was constructed using the Insight II modeling platform.18 Because the length of two ethylene glycol units roughly corresponds to the height of the R-CD cavity,9 and the preferred ratio in the rotaxane is 2.6 ( 0.1 oxyethylene units per R-CD,13 the rotaxane examined here consists of two R-CDs and a PEG chain formed by five repetitive oxyethylene units. In addition, since R-CD possesses two different hydroxyl groups at the two mouths of its cavity, it can be threaded onto a PEG chain following three possible spatial arrangementssviz.. HH, HT, and TT. Figure 1 provides a schematic representation of these motifs. In the course of the simulations, the carbon atoms pertaining to the stoppers of the PEG chain were coerced to their initial position using soft harmonic restraints. The initial distance separating the two R-CDs was set to 8 Å. In the case © 2008 American Chemical Society
Letters
Figure 1. Schematic representation of the rotaxane formed by a PEG chain and two adjacent R-CDs organized in three possible spatial arrangements: head-to-head (HH), head-to-tail (HT), and tail-to-tail (TT).
of the aqueous solution, a simulation cell was considered and replicated using periodic boundary conditions. The initial size of this cell, prior to equilibration, was 28.8 × 29.0 × 44.6 Å3 and contained up to 977 water molecules. All MD simulations were performed using the NAMD19 program with the CHARMM force field20 and the TIP3P water model.21 The isolated rotaxane was modeled in the canonical ensemble, whereas the aqueous solution was described in the isobaric-isothermal ensemble. The temperature and the pressure were maintained, respectively, at 298 K and 1 atm, employing Langevin dynamics and the Langevin piston method.22 The equations of motion were integrated with a 1-fs time step in the aqueous solution and a 0.5-fs time step in the gas phase. Lennard-Jones interactions were truncated beyond 14 Å. Longrange electrostatic forces were taken into account by means of the particle-mesh Ewald approach.23 The total simulation time for each spatial arrangement of adjacent R-CDs amounted to 2.5 ns in the gas phase and 5 ns in the aqueous solution. The preferential spatial arrangement of the contiguous R-CDs in the rotaxane was determined following the thermodynamic cycle depicted in Figure 2. The free-energy change associated with each leg of this thermodynamic cycle, wherein one macrocycle is annihilated, was estimated using the FEP methodology16 with the dual-topology paradigm.24 The initial structures utilized for the different FEP calculations were obtained from the preliminary MD simulations described previously. The pathway connecting the reference state of the transformation to the target state was broken down into 106 windows of uneven widths. On account of the singularity at the end point of the transformation, viz.. when λ tends toward 1.0, in addition to the first 99 windows of width ∆λ ) 0.01, 7 extra windows of decreasing width, i.e., ∆λ of 9 × 10-3-10-8, were utilized. In each window, 10 ps of equilibration were followed by 30 ps of data collection, corresponding to a total of 4.24 ns/run. Statistical accuracy of the free-energy calculations, in contrast with precision,25 was assessed by repeating each simulation three times, starting from different initial conditions. In the rotaxane, the R-CDs are threaded onto the PEG chain via two main interactions: the interaction of R-CD with PEG and that of the two R-CDs. The average interaction energies monitored in the course of the MD simulations for the three spatial arrangements, both in the gas phase and in the aqueous solution, are given in the Supporting Information. Comparing for a given spatial arrangement the interaction of R-CD with PEG in the two media, the effect of the surroundings appears
J. Phys. Chem. B, Vol. 112, No. 17, 2008 5269 to be marginal. On the contrary, the interaction of the two R-CDs is much stronger in the gas phase than in the aqueous solution. Moreover, difference in this interaction appears to stem primarily from changes in the electrostatic contribution. The latter contribution follows the ranking HH > HT > TT, identical to that of mutual interaction of R-CDs. The number of hydrogen bonds formed in the course of the MD simulations was monitored, and the corresponding results are provided in the Supporting Information. It can be seen that virtually no hydrogen bonds are formed between PEG and the other molecules of the system. This result can be largely explained by the fact that PEG is essentially engulfed in its entirety in the R-CD cavity. Investigation of the hydrogen bonds formed between the two R-CDs not unexpectedly underscores the agreement between the number of hydrogen bonds being formed with the electrostatic interaction energies. This indicates that, in the gas phase, the hydrogen-bonding interaction of the two R-CDs is prevalent. In the aqueous solution, however, neighboring water molecules weaken this interaction, presumably on account of a competition in the formation of hydrogen bonds between the macrocycles and between the latter and the aqueous solvent. Figure 3 depicts the structure of the rotaxane in an aqueous solution after MD simulation. The HH arrangement of the two R-CDs is highly symmetrical owing to the number of participating hydrogen bonds and the stronger mutual interaction of the macrocycles. A glimpse at the HT motif suggests that the two R-CDs tend to adopt a staggered arrangement. This trend is even more pronounced in the TT arrangement, where the distance separating the two R-CDs is the largest and the marked asymmetry of the assembly weakens the hydrogen-bonding interaction. Noteworthily, these structural features are similar in the gas phase and in the aqueous solution. Based on the point charges utilized in the present work, the dipole moment for the crystal structure of R-CD17 amounts to 9.6 D, which can be compared with the corresponding ab initio value.26 It is collinear to the longitudinal axis of the macrocycle and points from the head (wide mouth) of the latter to its tail (narrow mouth). In the HT motif, the individual dipole moments are essentially parallel, pointing in the same direction, which is conducive to enhanced hydration. In sharp contrast, the direction of the dipole moments in the TT and the HH arrangements is virtually antiparallel, which is anticipated to be in favor of the low-pressure gaseous state. To assess whether the system has reached equilibrium in the data collection part of the FEP calculations, evolution of the free energy with respect to the general extent parameter, λ, was investigated during the annihilation processes (see Figure 1S in the Supporting Information). The behavior of the free energy as a function of λ and the overlap of the thermodynamic ensembles embodied in their density of states suggest that the necessary conditions for an appropriate convergence of the FEP calculations are met. In order to probe the results of the FEP calculations, two distinct routes were selected to perform the transformation between the three possible spatial arrangements. One route connects the TT to the HH motifs via the HT intermediate, while the other route involves the TH intermediate. Each free-energy estimate corresponds to an average value over the three independent runs. The free-energy differences and the standard deviations are gathered in Figure 2. Closure of thermodynamic cycles constitutes a common criterion for appraising the convergence properties of free-energy calculations. The thermodynamic cycle of Figure 2 closes with a residual free energy
5270 J. Phys. Chem. B, Vol. 112, No. 17, 2008
Letters
Figure 2. Thermodynamic cycles used to evaluate the free-energy differences involved in the transformation of HH into TT via either the HT or the TH intermediate. In each leg of the peripheral cycles, an R-CD unit is annihilated and the corresponding free-energy difference, ∆Gij, where j * 1, is estimated using the FEP methodology. All values of ∆Gij in the gas phase and in the aqueous solution are displayed in standard and in boldface characters, respectively.
of +0.2 kcal/mol in the aqueous solution and of -0.3 cal/mol in the gas phase, suggestive that the FEP calculations have reached a reasonable convergence. It may be contended, however, that small residual free energies may conceal larger statistical errors associated with the individual legs of the cycle, likely to cancel out fortuitously. It is apparent from Figure 2 that the standard deviations of the free-energy changes ∆Gij, where j * 1, are generally small compared to the values from which they have been inferreds viz.. typically less than 1%. The same unfortunately is no longer true when opposed to the small net free-energy changes ∆Gi1, which correspond to a difference of two large numbers necessarily affected by a statistical error. In spite of repeated independent runs, reducing this error rapidly proved to be a computational challenge that emphasizes the limitations of FEP calculations when large free-energy changes are involved and reproducible sub-kBT accuracy, i.e., less than 0.6 kcal/mol at 300 K, is sought. In the gas phase, the free energy incurred in the conversion of TT into HT amounts to +1.2 kcal/mol, whereas that required for transforming HT into HH is equal to -3.6 kcal/mol. From this result, the stability of the three spatial arrangements in the rotaxane can be ranked as HH > TT > HT, where HH is the thermodynamically most favorable arrangement in vacuo. These estimates are supported by additional free-energy calculations following the alternate route involving the TH intermediate. From the above analysis and the ranking of the number of hydrogen bonds occurring in vacuo, it can be concluded that the HH preference is mainly due to the contribution of the hydrogen bonds formed between the secondary hydroxyl groups, which was also observed in experiments.14 In the aqueous solution, the free energy required for transforming TT into HT is -2.1 kcal/mol, whereas for HT into HH, it amounts to +1.8 kcal/mol. The stability of the three spatial arrangements in the rotaxane can then be ranked as HT > HH > TT, where HT is the thermodynamically most favorable arrangement in the aqueous solution, in agreement with previous results reported for the R-CD dimer.27 This preferred arrangement stems from the parallel alignment of the
individual dipole moments, which increases the total dipole moment of the assembly, conducive to an enhanced stability in the aqueous solution. The HT arrangement, however, also decreases the number of hydrogen bonds formed between the contiguous R-CDs. Preference for a given arrangement may, therefore, be viewed as a competition between the enhancement of the total dipole moment and the formation of hydrogen bonds between the adjacent CDs. In the aqueous solution, on account of the weakening of the hydrogen-bond interaction by neighboring water molecules, the former is likely to become an important factor in the stability of the arrangement, yielding a slight preference for the HT motif. Yet, the differences between all three motifs remain generally small. The experimental results8 for rotaxanes formed by two R-CDs and a linear aliphatic amino acids show that, in the crystal state, the contiguous R-CDs are arranged in an HH motif, which is in agreement with our prediction in the gas phase, but in the aqueous solution, the other two motifs TT and HT cannot be ruled out from NMR measurements. A molecular necklace formed by ∼22 R-CD units and a PEG chain was studied experimentally by Harada’s group.10,11 The NMR experiments performed, however, only provided macroscopic and averaged information about the spatial arrangement. Further study conducted by this group on the molecular conformation of an R-CD-PEG necklace dried in vacuo by STM method14 revealed that both HH and HT arrangements could coexist, and the corresponding formation ratio was determined to be 2:1. Since different environments, gas phase and aqueous solution, are involved in the preparation of the compound, the ratio of HH to HT may be qualitatively explained by a combination of the preference of HH in vacuo and the slight preference of HT in water. It should be noted that, in a rotaxane consisting of more than two R-CDs, HH and TT arrangements cannot exist independently. The stability of HH and TT motifs, therefore, cannot be compared directly. In summary, intermolecular interaction and hydrogen-bonding patterns have been examined to appreciate the role played by both the spatial arrangement of the macrocycles and their environment on the properties of the rotaxane. Modifying the
Letters
J. Phys. Chem. B, Vol. 112, No. 17, 2008 5271 Supporting Information Available: The average interaction energy and the average number of hydrogen bonds for the three spatial arrangements in the gas phase and in the aqueous solution. The evolution of the free energy change in different environments with respect to the general extent parameter, λ. This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes
Figure 3. Structures of rotaxane in the aqueous solution after MD simulation for the three spatial arrangements: (a) HH, (b) HT, and (c) TT.
macroscopic permittivity of the surroundings affects the mutual interaction of the contiguous R-CDs, primarily at the level of the electrostatic contribution. Hydrogen bonding is observed between the two R-CDs, and between R-CDs and water molecules. No hydrogen bond appears to form between PEG and neighboring molecules. Free-energy calculations reveal that, among the three possible arrangements of adjacent R-CDs, HH is preferred in the gas phase owing to hydrogen-bonding interaction in the secondary hydroxyl side, whereas in an aqueous solution, HT represents the most stable motif. This result can be rationalized by the directionality of the dipole moment borne by the macrocycles and, hence, the differential solvation properties of the supramolecular assemblies. Acknowledgment. This study is supported by the Natural Science Foundation of China (20325517 and 20573102). C.C. acknowledges the French Embassy in Beijing for travel support.
(1) Nepogodiev, S. A.; Stoddart, J. F. Chem. ReV. 1998, 98, 19591976. (2) Wenz, G.; Han, B. H.; Mu¨ller, A. Chem. ReV. 2006, 106, 782817. (3) Harada, A. Acc. Chem. Res. 2001, 34, 456-464. (4) Shigekawa, H.; Miyake, K.; Sumaoka, J.; Harada, A.; Komiyama, M. J. Am. Chem. Soc. 2000, 122, 5411-5412. (5) Pease, A. R.; Jeppesen, J. O.; Stoddart, J. F.; Luo, Y.; Collier, C. K.; Heath, J. R. Acc. Chem. Res. 2001, 34, 433-444. (6) Ballardini, R.; Balzani, V.; Credi, A.; Gandolfi, M. T.; Venturi, M. Acc. Chem. Res. 2001, 34, 445-455. (7) Schalley, C. A.; Beizai, K.; Vo¨gtle, F. Acc. Chem. Res. 2001, 34, 465-476. (8) Eliadou, K; Yannakopoulou, K.; Rontoyianni, A.; Mavridis I. M. J. Org. Chem. 1999, 64, 6217-6226. (9) Harada, A.; Kamachi, M. Macromolecules 1990, 23, 2821-2823. (10) Harada, A.; Li, J.; Kamach, M. Nature 1992, 356, 325-327. (11) Harada, A.; Li, J.; Kamachi, M. J. Am. Chem. Soc. 1994, 116, 3192-3196. (12) Peet, J.; Rusa, C. C.; Hunt, M. A.; Tonelli, A. E.; Balik, C. M. Macromolecules 2005, 38, 537-541. (13) Pozuelo, J.; Mendicuti, F.; Mattice, W. L. Macromolecules 1997, 30, 3685-3690. (14) Miyake, K.; Yasuda, S.; Harada, A.; Sumaoka, J.; Komiyama, M.; Shigekawa, H. J. Am. Chem. Soc. 2003, 125, 5080-5085. (15) Udachin, K. A.; Wilson, L. D.; Ripmeester, J. A. J. Am. Chem. Soc. 2000, 122, 12375-12376. (16) Zwanzig, R. W. J. Chem. Phys. 1954, 22, 1420-1426. (17) Manor, P. C.; Saenger, W. J. Am. Chem. Soc. 1974, 96, 36303639. (18) Accelrys Inc., Insight II User Guide; San Diego: Accelrys Software Inc., 2005. (19) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot C.; Skeel, R. D.; Kale L.; Schulten, K. J. Comput. Chem. 2005, 26, 1781-1802. (20) MacKerell, A. D., Jr.; Bashford, D.; Bellott, R. L.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E., III; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B. 1998, 102, 3586-3616. (21) Jorgensen, W. L.; Chandrsekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (22) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. J. Chem. Phys. 1995, 103, 4613-4621. (23) Darden, T. A.; York, D. M.; Pedersen, L. G. J. Chem. Phys. 1993, 98, 10089-10092. (24) Gao, J.; Kuczera, K.; Tidor, B.; Karplus, M. Science 1989, 244, 1069-1072. (25) Chipot, C., Pohorille, A., Eds. Free Energy Calculations: Theory and Applications in Chemistry and Biology; Springer: Berlin, 2007. (26) Li, X. S.; Liu, L.; Mu, T. W.; Guo, Q. X. Monatsh. Chem. 2000, 131, 849-855. (27) Bonnet, P.; Jaime, C.; Morin-Allory, L. J. Org. Chem. 2002, 67, 8602-8609.
