Abstract
The influence of the size, composition, and atomic distribution of linear guests on -cyclodextrin inclusion complex formation is clarified by means of a molecular dynamics simulation at constant temperature. The intermolecular energy is modelled by a Lennard-Jones potential, where the molecular composition is represented by various parameters and by a continuum description of the guest and cavity walls. It is concluded that the parameters related to the atomic size require minimum values for the confinement of linear molecules inside the cavity. The isomer with optimal affinity for -cyclodextrin as predicted by the free energy presents an asymmetrical molecular structure, and the position probability density shows that the isomer tends to insert the portion with largest atoms into the cavity, although the preferential binding site of the guest is not always located in regions of the host with maximum discriminatory power.
1. Introduction
Cyclodextrins (CDs) are macrocyclic molecules composed of glucose units (6 for -CD, 7 for β-CD, 8 for γ-CD, etc.) forming truncated cone-shaped compounds, giving rise to cavities of different internal diameters capable of including molecules of different structure, size, and composition. The capacity of CDs for catalysis and chiral recognition is due mainly to the formation of inclusion complexes, when some lipophilic part of a molecule enters the hydrophobic cyclodextrin cavity [1–3]. A well-known experimental finding is the existence of an appropriate size of the guest to reach maximum binding affinity in each CD [4–7], but there is no theoretical justification for this result and the factors influencing this affinity are not known. CDs and their inclusion complexes have been theoretically studied using different computational methods: molecular mechanics (MM) [8, 9], molecular dynamics (MD) [6, 10], and Monte Carlo simulations (MC) [11, 12], where all the atoms of both CDs and guest molecules have been described.
The present study is part of a broader work devoted to β-cyclodextrin inclusion complex formation, based on a continuum description of the guest and cavity walls. In this work, different types of guests have been considered: atoms and cyclic, spherical, and linear molecules. Moreover, two simulation methods (MM and MD) have been applied to analyse the potential energy surfaces and mobility of the guest inside and around the CD. We started by studying the interaction energy between β-CD and atoms and cyclic or spherical guest molecules, proposing an analytical model for this energy when the guest’s centre of mass is located along the cavity axis [13]. This model was extended for points away from the axis and used to examine the mobility of these types of guests by means of a molecular dynamics simulation [14]. The interaction between β-cyclodextrin and symmetrical linear molecules with different lengths was carried out by MM and MD [15, 16], where the composition of the guest was characterized by the parameters , . To study the interaction between β-cyclodextrin and asymmetrical linear molecules, new parameters were introduced to represent the molecular composition . These parameters allow different atomic distributions to be represented in the linear molecule which result in positional isomers. The potential surfaces, energy and configuration of inclusion complexes, and their dependence on molecular length, composition, and atomic distribution have already been analysed [17–19]. The curve joining the minimum potential energy for every plane = constant defines the penetration potential, which describes the variation in interaction energy when its path through the cavity is nonaxial. The penetration potential resembles a well potential or two minima separated by a potential barrier, because the interaction energy is deeper inside than outside the cavity, which represents an attractive force to include the guest into a complex. However the shape of this curve is a consequence of the atomic distribution, particularly the position of the larger atoms in the linear guest. The molecular composition is related to the magnitude of the interaction energy, and there are no molecular parameters for which the penetration potential presents special characteristics that can justify the affinity of some types of molecules for β-CD.
Theoretically, there are several necessary conditions for the isomers to be separated by CD: they must be able to enter the cavity (where the complexes formed are more stable), and the differences in retention time of these complexes have to be enough to allow experimental detection. In the mentioned simulation methods, the elution order in an isomeric separation is usually determined from the differences in free energy, where lower interaction energy means a more tightly bound analyte and this is linked with longer retention time. However, MD also allows us to determine the time spent by the guest inside the cavity (residence time) and provides a huge amount of statistical information to also calculate where and how selective binding takes place [20, 21]. The aim of this study was to determine the differences in residence time and free energy for the isomers of a linear molecule, so as to predict the preferential molecular parameters of this type of guests and thus facilitate separation by CD.
2. Simulation Method
The interaction energy between β-CD and a linear molecule is represented by the van der Waals contribution, modelled by a Lennard-Jones (6, 12) potential. It depends on two parameters , where is the depth of the well and is the position where the repulsive branch crosses zero [22–24].
In previous work, we analysed the dependence of the interaction energy between β-CD and linear guest molecules on the length, composition, and atomic distribution of the latter [19]. A minimum value of the molecular length was obtained for the differences in the interaction energy between isomers to be appreciable and thus their capacity to be separated. Therefore the present article deals with linear molecules with length Å. The interaction energy between the atoms of the system is represented by the parameters . With the ensemble , we try to generalize atoms of different size or type without the obscuring complications of too many material parameters, and represents the ratio of each pair , . The positions of different atoms on the linear molecule are represented by the parameters [19]. The same atomic distributions (7 cases) are used as previously to assess the influence of molecular stereochemistry on the interaction energy ( being ):(A)The atoms whose interaction with β-CD is characterized by being distributed together, from the top to the centre of the linear molecule(B)Distributions , , , with at (C), , , with around the centre of the guest(D), , , , (E)As in (B) but with the composition , , (F)As in (C) but with , , (G), , , , [19]The linear molecule can be symmetrical or asymmetrical depending on the atomic distribution and therefore on parameters . For instance, to represent a symmetrical molecule type (C), two parameters and are used, being that and (Figure 1). In contrast, to represent the asymmetrical atomic distribution of linear molecules type (A), only one parameter is needed whose value is .
The interaction energy between β-CD and a linear molecule with a composition characterized by and any atomic distribution can be calculated aswithwhere the guest molecule-CD interaction is represented by taking an average of the uniformly distributed atoms in the CD and in the linear molecule. is the superficial density of atoms in the CD cavity, is the linear density of atoms in the guest molecule, is the differential of length (on the linear molecule) located at , and is the differential of surface (on the cavity) located at (Figure 1). The interaction energy depends on the configuration of the guest molecule (orientation and centre of mass position) represented by . The CD is considered as a truncated cone in which the radius of the larger cone base is b ( Å for β-CD), the radius of the smaller top is a ( Å for β-CD) and the axial length is h ( Å).
The origin of the reference system is located at the centre of mass of the CD and the space-fixed frame is over the principal axis of the β-CD, where -axis is colinear with the cone axis (thus the plane is parallel to the cone base). The configuration of the linear guest is given by the coordinates of its centre of mass and the molecular orientation, defined by the polar angles formed with respect to the absolute frame . The polar angle is related to the orientation of the guest with respect to the cavity axis; for instance, in asymmetrical linear molecules type (A) if , the end of the molecule with values is pointing towards the cone top, and if this end is pointing towards the cone base.
In rigid-body dynamics, the molecular motion can be decomposed into two completely independent parts governed by a basic law of classical mechanics. The translational motion of the molecule depends on the total force acting on the body, whereas the rotation about the centre of mass depends on the total applied torque, this torque, and the angular velocity being perpendicular to the molecular axis at all times [25–27]. The orientational equations of motion are posed in terms of quaternions and they are solved using the same program applied to other types of molecules [16, 28]. To integrate the equations of motion, it is necessary to establish the initial conditions of the guest molecule: centre of mass position, orientation, and velocities (translational and rotational). The magnitude of the initial velocities depends on the temperature of the process (293 K), but the direction of the translational velocity in each trajectory and the orientation and initial centre of mass position are determined randomly. Since we are interested in studying the dependence of inclusion complex formation on molecular parameters, the same initial conditions are considered for the trajectories of the linear molecule, minimising if possible the influence of different initial conditions. Basically there are four relative positions for the molecules: with the centre of mass of the guest near either rim of the CD and with either end of the linear molecule pointing towards the cavity. Twelve trajectories (3 in each relative position) are calculated, the simulation time for each process is 1 ns with a step of 1 fs and the configuration and energies (kinetic and potential) are registered every 100 steps. However, the molecule can only enter the cavity and form an inclusion complex at initial positions of its centre of mass near the rims of the CD, as occurs with cyclic, spherical, and symmetrical linear molecules. When the initial centre of mass position is located outside the cavity walls, the molecule does not reach a stable configuration; it always stays around the cavity moving continually or even tending to move away. We use an in-house program written in Fortran, and the equations of motion to perform constant temperature molecular dynamics are integrated numerically using a variant of the leap-frog scheme (proposed by Brown and Clarke) [29], constraining the rotational and translational kinetic energies separately [30]. The data obtained from this molecular dynamics simulation are as follows: residence time, free energy, position probability, and preferential configuration of the inclusion complex and also how they depend on the molecular parameters (length, composition, and atomic distribution).
In order to analyse the ability of different linear molecules to form β-CD inclusion complexes, the preferential binding site and orientation of the guest are determined by means of the position probability density. We define a grid in which the distance between two consecutive points is 0.5 Å, and the number of guest positions in each volume element is the resulting number density for each trajectory and for the simulation [20, 21]. The position probability density is calculated by dividing the number density in a volume element by the total number of centre of mass positions of the guest. The most probable orientation is that with greatest number density.
The simulation is carried out for molecules with different lengths ( and 10 Å), composition, and atomic distributions (Table 1).
3. Results and Discussion
The initial conditions determine the integration of the equations of motion and also the movements of the guest molecule in each trajectory during the simulation time. However, they affected the simulation in different ways. Whereas the velocities hardly influence the number densities and the mean energy of the process, the greatest differences in these values are due to the initial configuration of the guest, which in turn determines the behaviour of the linear molecule in each process. During trajectories with an initial molecular configuration near the rims of the β-CD, the guest tended to enter the cavity, remain inside forming an inclusion complex, or exit from it after a period (residence time ), depending on the length, composition, and atomic distribution of the guest. To study the dependence of the residence time and free energy on molecular parameters, we averaged the values of these magnitudes obtained for the trajectories ( and , resp.). For a linear molecule with length and composition , the residence time varied from several ps to the total simulation time (1 ns), depending on the initial conditions of the trajectory. Therefore, in some cases, there can be a great difference in the mean values of , so it is represented asymptotically in the figures.
To analyse the influence of linear molecule size and composition on their capacity to stay confined in the β-CD, and therefore on their affinity for this host, we compared the results obtained with linear molecules of different size and composition with atomic distribution type (A). Figure 2(a) represents the mean value of against for molecules with length Å and different potential parameters . The results for linear molecules with L = 9 and 10 Å are similar and the following features can be pointed out:(a)The residence time was always greater for molecules whose composition is represented by two different potential parameters than for only one.(b)For a linear molecule with length L, the residence time grew with until it reaches a maximum of , after which it decreased. The ratio ranged from 0.3 to 0.5 depending on the values of .(c)For linear molecules with composition , the value of was independent of L (Table 2).Figure 2(b) shows the mean value of against (each with the corresponding , Table 1) for linear molecules with Å, kcal/mol, , and different lengths. The following features can be pointed out:(a)For a linear molecule with length L, the residence time grew slowly with until Å, for which rose rapidly, because starting from several initial configurations the molecule spends the simulation time inside the cavity.(b)The residence time for linear molecules with the same composition increased with L, as also seen in Table 2.The residence time depends on and , but there are some similarities in when = constant. Figure 2(c) shows the mean value of against (each with the corresponding , Table 1) for linear molecules with length Å, , and different values of . Figure 2(c) also shows the dependence of on for the same values of . The following features can be pointed out:(a)For = constant, the residence time became longer with and the greater is, the greater the increase in is.(b)For linear molecules whose composition has the same value of , the residence time became longer with .(c)The residence time increased with the difference , t being equal to the simulation time in several trajectories when and Å.
(a)
(b)
(c)
In previous work, the interaction energy was calculated between β-CD and linear guest molecules whose composition was represented by two pairs of potential parameters [17–19]. It was found that the magnitude of the interaction energy depends mainly on the molecular composition, and the possibility of confinement or temporary residence inside the cavity depends on the potential barrier which is a consequence of the atomic distribution. However, after studying the dependence of residence time on guest molecule composition, it can be concluded that the parameter requires a minimum value ( Å for Å, Å for Å) for the confinement of linear molecules inside the cavity. The remaining parameters representing the molecular composition were also restricted: the potential parameter must range from 2.5 to 2.7 Å and the ratio from 0.3 to 0.5, in agreement with the idea of establishing an appropriate size of guest molecule for each CD (according to the cavity size) to obtain optimum results in catalysis and chiral recognition.
One of the theoretical necessary conditions for the isomers to be separated by CD is that they must be able to form inclusion complexes; therefore we only considered the parameters representing the molecular composition that permits the confinement of the guest inside the cavity. Furthermore, the differences in retention time of the complexes formed with the isomers have to be enough to allow experimental detection. Therefore, the aim of this work was to determine the differences in residence time and free energy for the isomers. Greater differences in these magnitudes predict the preferential atomic distribution of the linear guest, so as to be better separated by CD. The mean value of t () in the simulation for a linear molecule with Å, Å, kcal/mol, Å, kcal/mol, , and different atomic distributions is included in Table 3. From the point of view of time differences, it can be concluded that the isomers (C), (D), (F), and (G) were more quickly eluted than (A), (B), or (E), although these are only average values and is the residence time spent by the guest in one, two, or more trajectories. The configurations (C), (D), (F), and (G) are characterized by a symmetrical distribution of atoms in the linear molecule, and therefore molecular asymmetry influences the capacity of the linear guest to remain confined inside the CD cavity. Linear molecules with Å, kcal/mol, Å, kcal/mol, , and length Å showed similar behaviour to molecules with the same composition and Å, although the residence times were shorter even for an asymmetrical distribution of atoms like (E). Thus, the isomers with better affinity for β-CD have an atomic distribution type (A) or (B), and the molecular length contributes less to the differences in residence time among isomers than their composition.
The free energy of the trajectory was calculated over the residence time instead of the simulation time, because in some cases the molecule tends to move away from the cavity and this contribution to the interaction energy is not important for discriminating between isomers. The values of the binding free energies obtained for the trajectories were averaged for a linear molecule with Å, Å, kcal/mol, Å, kcal/mol, , and different atomic distributions (Table 3). There was no proportional relation between the residence time and the free energy of the trajectory; lower interaction energies occur when a larger part of the linear molecule is inside the cavity. This is also the reason for the wide differences between the mean free energies of the processes, which are a consequence of the stable configuration of the guest in the complex and therefore related to the composition and atomic distribution of the isomer. The greatest difference in was found in the configurations (A) and (G), which also had the greatest difference in mean residence times. From the free energy point of view, configuration (A) shows the longest retention time but the results predicted for the other isomers differ from those based on the residence time, since the asymmetrical configuration (B) would be more quickly eluted than the symmetrical configurations (C) or (F). The mean energy for linear molecules with Å, kcal/mol, and length Å was different from Å (Table 3) because it depends on the guest configuration in the complex, principally on which those parts of linear molecules are located inside the cavity. This last factor depends in turn on the atomic distribution. From the energy point of view, configuration A had the longest retention time, as for molecules with Å and the same composition.
To determine the preferential binding site and orientation of the guest molecule, the position probability density was calculated [20, 21]. If the linear molecule tends to exit from the host after the residence time, it does not pass through the cavity in every process, and it can enter and exit from the same rim of the CD, giving rise to a position probability density of about 5%, distributed over the whole cavity (Figure 3(a)). When the guest remains inside the CD for the simulation time , the position probability density increases around the stable configuration of the linear molecule inside the cavity, where it can amount to about 85% (Figure 3(b)). Table 4 shows the centre of mass position and orientation of the guest at the preferential binding site for linear molecules with Å, Å, kcal/mol, Å, kcal/mol, and . It can be seen that the isomers (A), (B), or (E) form stable inclusion complexes with β-CD, where the guest is located inside the cavity and thus has better affinity for the host. The asymmetrical molecules tend to insert that part with larger-sized atoms inside the cavity, forming an angle of about 30° with the cavity axis when the guest centre of mass is near the narrow rim or about 10° near the wider rim. The preferential configuration of the isomer (E) was in accordance with that of the minimum interaction energy and was also found in regions of the host with maximum discriminatory power, that is, those with greatest energy differences (near the cavity walls) [19]. In contrast, the preferential configurations of isomers (A) and (B) were very different from those of minimum energy and were located near the cavity centre, where the potential energy is similar for all isomers and thus less selective. The preferential configuration of the isomers with greater affinity for β-CD was nearly the same for linear molecules with the same composition and atomic distribution, independently of their length (Table 4).
(a)
(b)
The suitability of the model presented in this study is confirmed by comparing it with experimental results such as ultrasonic absorption measurements of some alkylammonium ions with β-CD [31, 32], or a high-performance liquid chromatography (HPLC) procedure for the determination of the positional isomeric impurity in 2-[4-(1-hydroxy-4-[4-(hydroxydiphenylmethyl)-1-piperidinyl]-butyl)-phenyl]-2-methylpropionic acid with a β-CD bonded-phase column [33]. This compound has two isomers with molecular structure type A for the para-isomer and between (E) and (F) for the meta-isomer [19]. From the results obtained in the present work, the meta-isomer would elute first in this chromatographic system, as found experimentally. On the other hand, the authors concluded that propylammonium chloride does not form the inclusion complex, while butyl-, pentyl- and hexylammonium chloride do form the inclusion complex, incorporating the hydrophobic portion of the guest into the CD cavity. These molecules can be represented by the same potential parameters and different length and ratio , as occurs with acrylic esters [17–19]. Moreover, the aforementioned experimental results are compatible with the results obtained in the present work: molecules with this composition and length Å have a mean residence time of about 2 ps, whereas, for Å, the guest can remain confined in the cavity for the simulation time.
Lipkowitz et al. also studied where and how different analytes bind to permethylated β-CD in gas chromatography, using different computational methods (standard molecular dynamics and Monte Carlo simulations) [20, 21]. The structure of some considered analytes (acetyl-alkenes) is essentially linear and in this case the van der Waals forces exert the greatest intermolecular contribution to holding the host-guest complex together. They conclude that the preferred binding site for these analytes is the interior of the macrocycle, in accordance with the results obtained using the model presented in this study.
Other examples of isomers with molecular structure type (A) retained longer in β-CD are diethyl fumarate against diethyl maleate (between (B) and (D)) [19, 34].
4. Conclusions
We have studied the residence time, energy, and configuration differences in the complexes formed by β-CD and the isomers of a linear molecule. The influence of molecular parameters such as length, composition, and atomic distribution on β-CD inclusion complex formation was determined by means of a molecular dynamics simulation at constant temperature. The intermolecular energy was modelled by a Lennard-Jones potential where the molecular composition is represented by different parameters and a continuum description of the guest and cavity walls. From the study of the dependence of residence time on guest molecule composition, we conclude that the theoretical characteristics of molecules with optimal affinity for β-CD are as follows: a molecular composition involving potential parameters with rigid restrictions to allow the confinement of the guest inside the cavity and asymmetrical distributions of atoms in the linear molecule. The maximum discriminatory power of β-CD predicted by both free energy and residence time occurs between an asymmetrical molecular structure type (A) and the other isomers. Type (A) refers to the fact that those atoms whose interaction with β-CD is characterized by are distributed together, from the top to the centre of the linear molecule. The elution order predicted in these cases is similar when the guest forms temporary inclusion complexes with the host.
The position probability density shows that the guest tends to insert the portion with larger atoms into the cavity, the preferential configurations of the resultant inclusion complex coinciding with those of the energy minima in the penetration potential. The configuration of the guest at the preferential binding site is similar to linear molecules with the same composition and atomic distribution, independently of their length, although this preferential configuration is not always located in regions of the host with maximum discriminatory power.
Conflicts of Interest
The author declares that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The author is grateful to the Ministerio de Economía y Competitividad (FIS2013-41532-P and FIS2016-79596-P, AEI/FEDER, UE) for their generous financial support.