Revisiting the Internal Conformational Dynamics and Solvation Properties of Cyclodextrins

Simulações de dinâmica molecular são utilizadas para investigar a dinâmica conformacional interna e propriedades de solvatação das três ciclodextrinas naturais, α-, βe γ-ciclodextrina em água e temperatura ambiente. Estes oligossacarídeos derivados da glicose apresentam uma estrutura molecular peculiar que lhes confere a habilidade de formar complexos de inclusão e modificar as propriedades físico-químicas das moléculas complexadas. As características estruturais das ciclodextrinas em solução são determinantes para a complexação. As análises das trajetórias obtidas mostram que ligações de hidrogênio secundárias entre resíduos de glicose adjacentes estão presentes em solução e apresentam um caráter altamente dinâmico, onde ligações de hidrogênio alternativas podem ser formadas com moléculas de água. Apesar do baixo caráter hidrofílico, as cavidades apresentam-se solvatadas e o número de moléculas de água dentro da cavidade é aproximadamente duplicada por unidade de glicose adicionada ao macrociclo. Os tempos de residência para as moléculas de água dentro das cavidades são inversamente proporcionais ao seu respectivo tamanho.


Introduction
Cyclodextrins (CDs) have been extensively studied since their discovery in the late nineteenth century. 1 Their peculiar properties and behavior have interested scientists in many different fields. 2These cyclic oligosaccharides, obtained from starch by glucanotransferase degradation, are formed by a certain number of α-D-glucopyranose units linked by α(1→4) glycosidic bonds.2][3] The structural features of these saccharides allow them to form inclusion complexes with a great variety of compounds. 3Due to their ability to form host-guest complexes, resulting in significant changes in

Molecular dynamics simulations
All the simulations were performed using the Amber96 force field 54 in conjunction with a modified (optimized charge set, as described below) version of the GLYCAM98 55 force field, implemented in the TINKER program. 56Newton's equations of motion were integrated based on the Verlet algorithm 57 using a 1-fs time step.The RATTLE 58 method was applied to constrain all bond lengths involving hydrogen atoms.The simulations were carried out in the isothermal-isobaric (NPT) ensemble.The temperature was kept at 298 K by coupling the system to a heat bath 59 with a relaxation time of 0.1 ps.The pressure was kept at 1 atm by coupling the system to a pressure bath via isotropic coordinate scaling 59 with a relaxation time of 0.5 ps.A cutoff of 1.0 nm was applied to evaluate non-bonded interactions.The Particle Mesh Ewald method (PME) 60 was employed to account for long range electrostatic interactions.At least 2 ns MD simulations were performed for each CD at 298 K and 1 atm after an initial 100 ps period of heating and equilibration at the same temperature and constant volume.The atomic coordinates were saved every 0.1 ps and the analyses were carried out over the entire production runs, unless otherwise specified.

Molecular systems
A central problem dealing with carbohydrate simulations consists on the parameterization of a force field to correctly describe conformational preferences and intermolecular interactions simultaneously.][63][64][65][66][67][68] Modeling of CDs, just like other oligosaccharides, involves a particular problem regarding the electrostatic representation.Charge sets in carbohydrate force fields have been often derived based either on small fragments 62,65 or on the glucose molecule. 55However, CDs are oligomers of glucose linked by α(1→4) glycosidic bonds.The groups involved in the glycosidic linkage present a different chemical environment that is expected to alter its electrostatic properties.Aiming to remedy this problem, atomic partial charges were re-optimized to account for the specific molecular characteristics of the CDs.Hartree-Fock calculations at the 6-31G* level on an α(1→4) glucose trimer and a restrained electrostatic potential (RESP) 69 fit based on the optimized geometry were performed using the NWChem4.1 70 program to compute the partial charges.Values of 0.005 au for the harmonic and 0.001 au for the hyperbolic restraints, respectively, were used for all RESP fittings.The final charge set derived for the constituent atoms of the CDs molecules is compared to the corresponding charges of the GLYCAM98 55 force field in Table 1.The numbering used to represent the atoms are displayed in Figure 1, which shows a schematic drawing of the α-CD molecule.
The two charge sets are relatively similar.As expected, the most significant differences are found in the atoms near the glycosidic bond.Charge derivation by means of RESP 69 fitting tends to maintain the dipole moment of a molecule close to its quantum-mechanical value and should lead to a realistic description of intermolecular interactions during the simulation.The charge-derivation procedure was performed according to the Amber96/ GLYCAM98 54,55 force fields philosophy in order to preserve compatibility.
The optimized geometry of a modeled α-CD using the re-optimized force field presented excellent structural agreement with an ab initio optimized structure via HF/ 3-21G* basis set (root mean square atomic position deviation (RMSD) of 0.047 nm between the two structures).The empirical energy minimization was performed using the steepest-descent method implemented in the Tinker package 56 and the ab initio HF/3-21G* calculation was performed using the Gaussian98 package. 71In both calculations an energy gradient of 10 -4 kcal mol -1 Å -1 was used as the convergence criterion.
The descriptive models for α-, βand γ-CD have been created using the derived point charge set along with the Amber96/GLYCAM98 54,55 force field parameters.The systems consisted of one CD molecule in a cubic box containing previously equilibrated TIP3P 72 water molecules.The initial coordinates of the CDs were obtained from the Cambridge Crystallographic Database (http://www.ccdc.cam.ac.uk/).In the simulations of αand β-CD the initial dimensions of the box containing 1450 water molecules were 3.53 × 3.53 × 3.53 nm 3 .The γ-CD simulation involved 1600 water molecules resulting in an initial box of 3.70 × 3.70 × 3.70 nm 3 .The energy of the initial configuration was further minimized in order to relax unfavorable energetic interactions.All atoms of the systems were represented explicitly, CD molecules were considered fully flexible, with constraints only in the bonds involving hydrogen atoms.

Analysis
Time-dependent RMSD was evaluated for the three CD molecules during the simulations.Two different atom sets were considered: i) all solute atoms and ii) ring backbone, excluding hydrogen atoms and lateral hydroxymethyl and hydroxyl groups.The coordinates from the geometry-optimized structures were used as reference structures for the fitting.The corresponding distributions for all atoms were also calculated.The time series of the two glycosidic torsion angles C1 n -O4 n+1 -C4 n+1 -C3 n+1 [ψ], and O5 n -C1 n -O4 n+1 -C4 n+1 [φ] (Figure 1) were evaluated (n refers to any given glucose residue in the macrocycle) for all glucose residues of each CD molecule.The internal hydrogen bonds (H-bonds) between secondary hydroxyl groups (OH2, OH3) in the CD molecules (referred throughout the text as secondary Hbonds) were evaluated considering the time series of the distances between the atoms O2 n …O3 n+1 , O2 n …HO3 n+1 and O3 n+1 …HO2 n of neighbor glucose residues during the simulation (n refers to any given glucose residue in the macrocycle).A 0.2 and 0.3 nm cutoff distances were used  as H-bond definition criteria for hydrogen-oxygen and oxygen-oxygen pairs, respectively.Radial distribution functions, g(r), were evaluated for water around all the oxygen atoms of the CD molecules.The solvation properties of the cavities have been characterized by a number of analyses.The center of the cavity was defined evaluating the geometric center of all the atoms comprising the CD molecule.Radial distribution functions, g(r), and the number of water molecules (coordination number) were calculated as a function of the distance from the center of the cavity.Additionally, the probability of finding water molecules as a function of the distance from the center of the CDs was calculated.
From these analyses an estimation of the average radius of the cavity for each CD molecule has been drawn.Residence time or correlation time for the water molecules inside the cavity was calculated during the simulation as a function of an arbitrary radial cutoff (in the range 0.05 to 1.0 nm from the center of the cavities).The distance of some selected water molecules from the cavity center as well as the number of water molecules inside the region defined by the cavity were calculated as a function of the simulation time.

Structure of the CDs
After the initial equilibration period of 100 ps no significant drifts were observed in the energy (potential, kinetic and total energy) as well as in the density of the systems (data not shown).The RMSD analyses, performed for α-, βand γ-CD during the simulations, indicated no significant difference among the structural fluctuations of the three molecules in solution, as can be seen from the RMSD time evolution (Figure 2, left panel).The Gaussian-like distributions of the corresponding RMSDs for all atoms (Figure 2, right panel) show that all three simulations have reached structural equilibrium.Both all-atoms and backbone RMSDs present the same profile with a small shift in amplitude indicating that the structural changes in CDs come mainly from the fluctuations in the ring backbone.Positional variations of the hydroxyl and hydroxymethyl groups, normally highly-mobile groups in free glucose, here seem to play a less central role for the overall structural fluctuations of these molecules.The phenomenon of conformational flexibility in CDs is still nowadays matter of debate. 73For a long while CD molecules have been considered conformationally quite rigid, but this view has undergone revisions. 73,74Supporting the latter is the ability of CDs to accommodate themselves to the spatial and electronic character of a wide variety of guests.It is hard to conceive such structural adaptation in naturally highly-rigid structures.3][4][5][6]73 Most of the discussions about CDs flexibility have been based on crystallographic information obtained from solid state studies or gas-phase analyses, as discussed in the review by Dodziuk. 73The presence of H-bonds between the neighbor secondary hydroxyl groups OH2 and OH3 is generally considered a factor that enhance the rigidity of these cyclic molecules. 10,75In studies concerning CDs conformational flexibility, the inter-glucose H-bonds are normally considered to be the dominant factor for the macrocycle rigidity, having influence even on the conformational space defined by ψ and φ. 75 In order to probe these assumptions, we have performed analyses of both glycosidic torsion angles (ψ and φ) (Table 2) and inter-glucose H-bonds as a function of time (Figure 3).
The glycosidic torsion angles (Table 2) are clearly more flexible than the ring torsion angles, whose values do not deviate too much from the mean values characterizing the stable chair conformation ( 4 C 1 ) of the glucose residues (data not shown for conciseness).The RMSD analyses have shown that structural fluctuations in the CD molecules were defined basically by the fluctuations in the ring backbone (Figure 2).Since deviations from the average values of the internal torsion  angles in the glucose residues are very small, it can be concluded that the structural flexibility in CD molecules comes mainly from fluctuations of the glycosidic torsion angles.The average values for ψ and φ of the three CD molecules are shown in Table 2 and are in agreement with the calculated mean values obtained from a statistical analysis of these angles in CD crystal structures. 76Interestingly, average standard deviations for β-CD arguably display a smaller internal motion amplitude than its counter parts.While such observation relies on rather small differences over the analyzed final 2 ns, a recent study combining pulse-field-gradient-spin NMR with MD simulations has reached the same conclusion. 43istances between the atoms O2 n …HO3 n+1 and O3 n+1 …HO2 n in each glucose unit of α-, βand γ-CD, as a function time, are illustrated in Figure 3 (a, b and c, respectively).Inter-glucose H-bonds between the secondary hydroxyl groups (OH2 and OH3) are present most of the time in the different glucose residues of the three CD molecules.However, it is important to notice the dynamical behavior of these H-bonds, with apparently no preferential direction (H-bonds are formed between O2…HO3 as frequently as between O3…HO2).This dynamical behavior, so called flip-flop phenomenon, has been already observed experimentally in the crystalline state of β-CD 77,78 and reported in previous MD simulations. 32Our findings are in agreement with references 32 and 77-78 and suggest that the flip-flop phenomenon may be an intrinsic property of these cyclic molecules.It is worth to point that the dynamics of these interactions are probably influenced by the nature of the solvent.In fact, solvation effects on the conformation of carbohydrates are a well-known phenomenon (widely discussed in references 67 and 68).Very likely, when an internal H-bond is broken almost simultaneously there is the formation of an H-bond with water, hence restoring the energetic equilibrium of the system.Therefore, H-bonds between water molecules and polar groups of the CDs would act contributing to increase its conformational diversity.Furthering the assumption, solvents with little or no possibility for intermolecular H-bonds, might hinder such degrees of freedom to certain preferred configurations.

Hydration and interaction with water
Radial distribution functions, g(r), for water around the oxygen atoms of α-, βand γ-CD molecules are shown in Figure 4 (a, b and c, respectively).g(r) for atoms O2, O3, O5 and O6 present structural features characteristic of h-bonding to water molecules (correlation peaks at 0.18 nm and 0.28 nm for hydrogen and oxygen atoms of water, respectively).These findings confirm our previous statement that energetically stable arrangements in these systems are also possible when H-bonds are formed between water and CD.Nonetheless, the correlation peaks for the water hydrogen around the oxygen atoms of the CDs sometimes present intensities smaller than one, showing that the backbone hinder water molecules from freely approaching the hydrophilic sites.This hindrance becomes so effective for the O4 atoms that no significant correlation was observed between these sites and the surrounding water molecules.The effect is also evident for O5 (located inside the ring structure) and O2 and O3, reflecting the presence of H-bonds between OH2 and OH3 groups.
The g(r) functions calculated for water around the cavity center of the three CDs (Figure 5a) show two regions of correlation, one closer to the center, indicating the existence of water molecules inside the cavity, and another from 0.85 to 1.1 nm, corresponding to the hydration shell molecules outside the cavity.To fully describe the correlation between water and CDs, coordination number and radial probabilities have also been evaluated.The radial probability of finding water around the geometric center of the CD molecules increases until this distance approaches the inner surface of the cavity (Figure 5b, inset graph), since the CD backbone may be pictured as an excluded volume to the water molecules.Accordingly, the probability must increase again after the distance becomes larger than the outer radius of the CD (Figure 5b).Therefore, the first minimum in the radial probability of finding a water molecule around the geometric center of the CD molecules may be considered as a parameter to evaluate the size of their cavities.The radii found in this way were 0.25, 0.29 and 0.37 nm for α-, βand γ-CD, respectively (Figure 5b, inset graph), in good agreement with experimental data obtained from hydrated crystals (0.24-0.27 nm, 0.30-0.33nm, 0.38-0.42nm for α-, βand γ-CD, respectively 75 ).Note that the predicted radii are towards the lower end value measured in the crystal.Partial solvation and intrinsic flexibility restraints of the crystal are expected to favor intra-molecular interactions in the CDs, by fully-solvating these molecules an increase in the entropy of the hydroxyl groups would probably lead to an average slightly smaller cavity radius.Coordination numbers show that the larger the cavity radius, the more water molecules are found inside it (Figure 5c).
The number of water molecules inside the cavity as a function of the time in the last 1 ns of each simulation (Figure 6, left panels) shows the correlation between the cavity size and the number of guest water molecules.In α-CD, most of the time there are 2 water molecules inside the cavity while in β-CD this number increases to an average of 4. In γ-CD a further increase to an average of 8 water molecules inside the cavity is observed.
Given the geometric center of two molecules, it is possible to evaluate the autocorrelation function that determines the probability of a pair of two such molecules to still exist after a time lag t, given that it existed at time t 0 . 79The autocorrelation function integral is a characteristic time, which may be regarded as a residence time of the molecular dimer.Residence times for water as a function of an arbitrarily-defined radial cutoff from the center of the cavity of the CD molecules are shown in Figure 7.As can be seen, α-CD presents a different behavior when compared to βand γ-CD.Water molecules present a residence time of approximately 98 ps in the spherical region defining the cavity of α-CD (ca.0.25 nm).This residence time is effectively longer than the residence time observed for water in the cavity of β-CD (27 ps) and γ-CD (25 ps) defined approximately by the spherical radii 0.29  nm and 0.37 nm, respectively.These results indicate that the water molecules entering in the cavity of α-CD spend, on average, more time than they would do in the cavities of βand γ-CD.This finding may be related with the size of the cavities, i.e., bigger cavities permit a higher flux of water molecules, as shown in Figures 5c and 6. α-CD cavity holds a smaller number of water molecules that can eventually remain trapped inside the cavity, giving rise to the elevated residence time.It is worth to note that the trapped molecules continue to exhibit relatively high rotational and translational motion inside of the cavity.Probably due to such mobility, we were unable to identify any specific interactions responsible for this observation.The small likelihood of γ-CD to trap water molecules within the cavity is also evident and agree with a recently reported experimental study. 80The dynamics of the water molecules through the CDs cavities is better appreciated in Figure 6, right panels.The distance between the four water molecules presenting the longest residence times inside the cavity and the geometric center of the CDs is plotted as a function of time during the final 1 ns of each simulation (only the four most stable dimers are shown, i.e., the number of water molecules in the cavities at any given time is actually larger than the  one shown in the plot).α-CD presents stable dimmers most of the time, while βand γ-CD have some periods in which more transient dimers were found (Figure 6, right panels).

Conclusions
MD simulations were used to investigate the internal conformational dynamics and solvation of the three natural CDs, α-, βand γ-CD, in water at room temperature.A modified version of Amber96/GLYCAM98 force field was developed by the re-optimization of the partial charges to account for some of the specific molecular characteristics of the CDs.The results obtained have shown that these molecules present high conformational freedom in solution.Overall, flexibility seems to be defined by fluctuations in the glycosidic torsion angles (φ and ψ) and interaction with water.The ring torsion angles do not scatter too much from the average values defining the stable 4 C 1 chair conformation for the glucose residues.Secondary H-bonds are formed in solution and present a fast dynamic flip-flop character.This observation agrees with other studies (experimental and theoretical) on the structure of CDs either in solution or in the solid state.The molecule of γ-CD presents in solution a more symmetric structure when compared to α-, β-CD, as shown by the persistent chain of secondary H-bonds present in this molecule.We speculate that due to solvation effects, the secondary H-bonds do not confer as much rigidity to the CD molecules as in the solid state.If that would be the case, interactions between the CDs and high-polar solvents, such as water, may increase molecular flexibility allowing dynamical rearrangements of the structure with minimum energetic costs.Additional simulations of the CDs in low-dielectric solvents would be required to asses such assumption.The CDs cavities are known to be relatively hydrophobic, however the present simulation results show that the cavities are able to host water molecules, whose number and residence times depend on the size of the cavity.In summary, conformational variability in the molecular structure of CDs is assisted by interactions with water molecules suggesting that flexibility is strongly dependent on the chemical environment.The formation of H-bonds with neighbor glucose residues can significantly affect their conformational changes, and certainly the intrinsic flexibility of the CD molecules has a strong influence on the complexation mechanisms.Additionally, the findings presented here offer ground basis information for the novel design of carbohydrate-derived cyclic compounds and insights into their complexation mechanisms.

Figure 2 .
Figure 2. RMSD time series (left panel) to α-, βand γ-CD during the simulations performed.Black lines show the RMSD for all atoms of the CD molecule, gray lines show the RMSD for the ring backbone atoms.Normalized statistical distributions for the all-atoms RMSD values are shown in the right panel.

Figure 3 .
Figure 3.Time evolution from 0 to 2 ns, at a 0.1-ps interval, of the distances between the O2 n …HO3 n+1 (in black) and O3 n+1 …HO2 n (in gray) atoms of α-(A), β-(B), and γ-CD (C) during the simulation.The number of plots per panel corresponds to the number of glucose residues in each macrocycle.

Figure 5 .
Figure 5. Radial distribution function, g(r) (top), radial probability (RP) (mid) and coordination number (CN) (bottom panel) for water from the center of the cavity of the CD molecules.The two inset graphs present the radial probability and coordination number in a scale from 0 to 0.4 nm to amplify the differences.

Figure 4 .
Figure 4. Radial distribution functions, g(r), for water around the oxygen atoms of α-(A), β-(B) and γ-CD (C).The number of plots per panel corresponds to the number of oxygen atoms in each glucose residue.

Figure 6 .
Figure 6.Number of water molecules inside the CD cavity as a function of time, for the final 1 ns of each simulation (left panel), and the distance between the four water molecules presenting the longest residence time inside the cavity and the geometric center of the CDs as a function of time during the final 1 ns of each simulation (right panel, only the four most stable dimmers are shown, arbitrarily represented by four different symbols).

Figure 7 .
Figure 7. Residence time for water molecules as a function of an arbitrarily-defined radial cutoff (0.05 to 1.0 nm) from the geometrical center of the cavity of the CD molecules.

Table 1 .
Partial atomic point charges derived for the glucose atoms in the CD molecules and the corresponding point charges for the glucose molecule in the GLYCAM98 force field