- Email: [email protected]
Molecular dynamics study of the influence of solvents on the chiral discrimination of alanine enantiomers by β-cyclodextrin
Tetrahedron: Asymmetry 24 (2013) 1198–1206
Contents lists available at ScienceDirect
Tetrahedron: Asymmetry journal homepage: www.elsevier.com/locate/tetasy
Molecular dynamics study of the influence of solvents on the chiral discrimination of alanine enantiomers by b-cyclodextrin Elena Alvira ⇑ Departamento de Física Fundamental II, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain
a r t i c l e
i n f o
Article history: Received 20 June 2013 Accepted 13 August 2013
a b s t r a c t The influence of solvents on the separation of alanine enantiomers using b-cyclodextrin as a chiral selector was studied by means of a molecular dynamics simulation at a constant temperature. The potential energy of the interaction is modelled by the AMBER force field, where different polar and non-polar solvents are represented by the dielectric constant e and two configurations for the amino acid derived from its electric charge distribution: the AMBER data base or its zwitterion state. The L enantiomer has more positions inside the cavity of a b-cyclodextrin where it is more stable than the D-enantiomer in vacuo and solution, except for solvents such as hydrocarbons in which most positions of the D-alanine inside and outside the cavity are more stable. In all cases, the greatest differences are located near the cavity walls. Molecular dynamics simulations show that Ala is able to form inclusion complexes with b-cyclodextrin in vacuo and in solvents such as hydrocarbons, benzene, acetone, ethanol or water. The chiral discrimination of Ala by b-cyclodextrin is mainly due to the adaptation of the guest to the host in the presence of non-polar agents, whereas the nonbonded interaction is the driving force for zwitterions. The elution order depends on the type of organic modifiers while a reversal of the enantiomeric elution order can be observed in solvents with higher dielectric constants. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction Cyclodextrins (CDs) are macrocyclic molecules composed of glucose units (6 for a-CD, 7 for b CD and 8 for c-CD, etc.) that form truncated cone-shaped compounds. These have cavities of different internal diameters capable of including molecules of different structure, size and composition. Their capacity for catalysis and chiral recognition is mainly due to the formation of inclusion complexes, when some lipophilic part of a molecule enters the hydrophobic cyclodextrin cavity.1,2 Such capabilities have been applied to various research fields such as solubility enhancement, drug delivery, chemical protection, separation technology and supramolecular chemistry.3,4 High-performance liquid chromatography (HPLC) and capillary electrophoresis (CE) using CDs as chiral selectors are frequently used for the enantioseparation of chiral compounds. Amino acids are biologically important organic compounds, which perform critical roles as neurotransmitters, and in transport and synthesis. In the form of proteins they comprise the second largest component (other than water) of human muscles, cells and other tissues. The L- and D-forms of amino acids differently adsorbed on a chiral surface constitute a large group of substances, ⇑ Tel.: +34 922318258; fax: +34 922318320. E-mail address: [email protected] 0957-4166/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.tetasy.2013.08.006
which are separable by CDs. Most of the experimental and theoretical studies carried out with CDs deal with amino acids such as Val, Leu, Ile, Tyr, Phe, Trp and mixtures of them, and aromatic and other types of amino acid derivatives. The reason for this selection of potential guest molecules is presumably the known complexing affinity of CDs for molecules of appropriate sizes and for aromatic functional groups in aqueous solutions.5,6 In fact, there is experimental evidence of inclusion complex formation by these amino acids with CDs in solution,7,8 as well as in gas phase.9 The amino acid Ala is smaller than the aforementioned amino acids but Ramirez et al. also demonstrated the inclusion complex formation of Ala with CDs in the gas phase.9 Experimental results also reveal the influence of organic modifiers on the enantiomeric elution order for some amino acids using CDs as chiral selectors.10 Theoretical methods such as molecular mechanics11 and molecular dynamics simulations9,12 are used to help explain the enantiodifferentiation mechanism and predict the separation process. In these theoretical studies, the presence of organic modifiers is represented by different values of the dielectric constant in the electrostatic contribution to the interaction energy. However under some pH conditions or with polar solvents, amino acids appear as zwitterions with the NHþ 3 and COO groups instead of NH2 and COOH, and this implies not only a different electric charge distribution, but also a different molecular configuration. Therefore, we can represent polar solvents such as ethanol and non-polar solvents
1199
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
such as hydrocarbons not only by different values of dielectric constant, but also by different configurations of the amino acids. This allows us to study the influence of solvent polarity on the separation of Ala by b-CD. In previous work, we analysed the chiral discrimination of alanine enantiomers with cellulose (crystalline and amorphous) using molecular mechanics simulations.13,14 The amorphous cellulose exhibits spatial helicoidal arrangements of the glucose units, allowing partial inclusion of the analyte molecules. We concluded that chiral discrimination increases when Ala is located inside the cavity, and it is maximum for similar sizes of helix radii and pitch. Cyclodextrins and cellulose derivatives are usually considered different types of chiral stationary phases (CSP)15,16 although they are composed of glucose rings (differently linked). The line of demarcation between them (suggested by Wainer and Lipkowitz) is the extent of guest inclusion into a chiral cavity. Herein our aim was to study theoretically inclusion complex formation and separation process of Ala with b-CD, in the presence of different solvents. We determined the influence of organic modifiers on the inclusion complex formation and elution order of the separation of alanine enantiomers by b-CD using molecular dynamics (MD). The potential energy of the interaction is modelled by the AMBER force field, although in recent years there has been an increasing development in new force fields that are applicable to carbohydrates, such as GLYCAM17 or q4md-CD, built from a combination of the GLYCAM04 and Amber99SB force fields to correctly describe the geometrical, structural, dynamical and hydrogen bonding aspects of heterogeneous cyclodextrin based systems.18 In the molecular simulation methods, the elution order in an isomeric separation is usually determined from the differences in free energy. The smallest value means a more tightly bound analyte and this is linked with the longest retention time. The chromatographic retention time is an experimental result that depends on macroscopic parameters; in contrast, free energy is calculated theoretically from microscopic modelling on the molecular scale; therefore it is only possible to compare the two sets of results qualitatively. However, MD also allows us to determine the time spent by the guest inside the cavity and provides a huge amount of statistical information to calculate where and how selective binding takes place.19,20 These results are not usually provided by the theoretical studies of separations of amino acids by CDs, but they can help to optimise these processes. The interaction potential and simulation method used herein are presented in Section 2. Section 3 is devoted to the analysis of the interaction energies between b-CD and the alanine enantiomers, and a discussion of the main results of molecular dynamics simulations related to inclusion complex formation and elution order.
2. The Model 2.1. Expression of the interaction potential In order to determine the interaction energy between b-CD and Ala, we used the force field proposed by Weiner et al. for the molecular mechanics simulation of nucleic acids and proteins.21,22 The interaction potential Etotal is the sum of the intramolecular Eintra and intermolecular Einter energies:
Etotal ¼ Eintra þ Einter
ð1Þ
The intramolecular energy is modelled by the equation:
Eintra ¼
X
kr ðr req Þ2 þ
bonds
X angles
kh ðh heq Þ2 þ
X Vn ½1 þ cosðn/ cÞ 2 dihedrals ð2Þ
The bond stretching and bending functions are quadratic, which allows an adequate description of the structure and energies for proteins and nucleic acids. The torsional energy is represented by a Fourier series approach. The intermolecular energy includes the nonbonded (Lennard–Jones), electrostatic and H-bond terms:
Einter ¼
" X Aij
Bij
Rij
Rij
12
i
þ 6
# " # X C ij qi qj Dij þ 12 eRij R10 Hbonds Rij ij
ð3Þ
where Rij represents the distance between the ith atom of the guest and the jth atom of b-CD. The atomic coordinates of b-CD are taken from the literature,23 and its net atomic charges and the AMBER force field parameters (Aij, Bij, Cij and Dij) are employed.21,22 We also considered the AMBER data as the molecular configuration and atomic point charges for the amino acid,21 but under some pH conditions or with polar solvents, the amino acid appears as a zwitterion (with the NHþ groups 3 and COO instead of NH2 and COOH). This presents different atomic point charges24 and molecular configurations from the AMBER force field. We represent P as the zwitterion configuration, NP as the non-polar configuration of Ala and the solvents by different values of the dielectric constant (Table 1). We considered all of the atoms in the molecule of alanine because the difference between the two configurations is related to the positions of some H atoms that can contribute decisively to the formation of hydrogen bonding between host and guest, and may be reflected in the interaction energy Etotal. We placed the origin of the reference system at the centre of the mass of the cavity and the space-fixed frame over the principal axis of the CD, in which the inertia tensor is diagonal. The position of Ala is given by the coordinates of its centre of mass, and the molecular orientation is defined by the relation between its principal body-fixed system and the axis system fixed in space. The potential energy for each enantiomer is determined by Eq. (1) for different positions and orientations of the guest centre of the mass, inside and outside the CD. In each plane Z = constant, about 2500 points are explored, the range of variation along the Z axis is 10 Å (with a path of 0.1 Å) and about 23,000 orientations of the guest are considered at each grid point. The results obtained are represented as potential energy surfaces, penetration potentials and the inclusion complex configuration, as in previous work.25,26 The curve joining the minimum intermolecular energy for every plane Z = constant defines the penetration potential W, which describes the variation in Einter when its path through the cavity is non-axial. The position of the guest molecule for which we obtain the absolute minimum potential energy Etotal gives the geometry of the inclusion complex for the different solvents. However the enantiomers do not generally adopt the lowest energy configuration due to the separation process, and can form complexes with configurations different from the minima that lead to a different enantiomeric elution order. We considered the regions with maximum chiral discrimination to be those localised in the positions of potential surfaces with greatest differences in energy, but we assigned to each grid point the average Boltzmann energy corresponding to different guest orientations, instead of the lowest energy.19,20 2.2. Simulation method The molecular dynamics (MD) is based on the resolution of classical equations of motion to determine the trajectories of the particles, depending on the initial conditions of the guest molecule: position, orientation (Euler angles) and velocities (translational and rotational). The magnitude of the initial velocities depends on the temperature of the process (293 K) but the directions of the translational and rotational velocities in each trajectory, as well as the initial centre of mass position are determined randomly. The
1200
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
Table 1 The minimum interaction energy Emin obtained with the AMBER force field for each enantiomer and the different contributions: Lennard–Jones ELJ, electrostatic Eele, hydrogen bonding EH-bond, bond stretching Ebond, angle bending Eangle and torsional energy Etorsion Con.
P P P NP NP
e 1 26 80 2 21
Emin (kcal/mol)
ELJ (kcal/mol)
Eele (kcal/mol)
EH-bond (kcal/mol)
Ebond (kcal/mol)
Eangle (kcal/mol)
Etorsion (kcal/mol)
L
D
L
D
L
D
L
D
L
D
L
D
L
D
52.95 4.67 3.47 5.24 4.00
51.99 4.51 3.28 3.71 3.76
10.11 10.71 10.51 2.91 11.23
8.97 10.67 10.79 3.20 10.63
50.77 1.89 0.55 8.31 0.88
51.00 1.83 0.53 7.26 0.85
0.48 0.49 0.91 0.50 0.00
0.49 0.47 0.49 0.00 0.49
1.72 1.74 1.80 5.09 4.97
1.66 1.66 1.74 5.21 5.01
6.61 6.59 6.62 1.30 1.29
6.61 6.61 6.59 1.34 1.30
0.08 0.08 0.08 0.09 0.09
0.20 0.20 0.20 0.20 0.20
simulation time for each trajectory is 3 ns with a step of 1 fs and the configuration and energies (kinetic and potential) were written every 100 steps. We used an in-house programme written in Fortran and the equations of motion in order to ensure constant temperature molecular dynamics. They are integrated numerically, using a variant of the leap-frog scheme (proposed by Brown and Clarke),27 which constrains the rotational and translational kinetic energies separately.28 This programme has been used to study the inclusion complex formation of different types of guests in b-CD.25,26 We calculated 6 different trajectories with initial configurations of the guest outside the b-CD: 3 near the primary (narrow end) and 3 near the secondary rims (wide end) of the CD. In order to avoid the influence of different initial conditions, in each trajectory the alanine enantiomers have the same initial centre of mass position and molecular orientation, so they differ in their atomic positions. In order to analyse the ability of Ala to form b-CD inclusion complexes with different solvents, we determined the preferential binding site of the guest molecule in the simulation by means of the probability density of presence in a volume element. We defined 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 guest.19,20 The probability density of presence is calculated by dividing the number density in a volume element by the total number of centre of mass positions of the guest. We determined the elution order by calculating DF = FL FD, with F being the binding free energy in the simulation for each enantiomer.29–31 In order to compare the time spent by the guest inside the cavity in different solvents with F, we averaged the values of the time obtained in MD simulations.
3. Results and discussion 3.1. Interaction energy The values of dielectric constant e considered herein are similar to those of agents such as hydrocarbons, benzene, toluene, acetone, ethanol or water. The interaction energy Etotal between b-CD and Ala for these values of e is negative inside the cavity with any of the guest configurations (NP or P), which means that the amino acid is able to form inclusion complexes. However, in vacuo (e = 1) the energy Etotal is positive inside the b-CD cavity when the configuration of Ala is NP, in this case the inclusion complex can only be formed by the configuration P of zwitterions. Ramirez et al. assumed a dielectric of 1.0 in the molecular modelling simulations where they provided theoretical evidence of inclusion complex formation of several amino acids in the gas phase,9 but they did not study alanine enantiomers. Figure 1 represents the penetration potential W, Lennard–Jones (LJ) and electrostatic (ELE) terms for the interaction between b-CD and the enantiomers of Ala, with different solvents. In the case of a polar configuration
with e = 1 (Fig. 1a), e = 26 (Fig. 1b) and e = 80 (Fig. 1c), as e increases the electrostatic term contributes less to the total energy, and we can see that this contribution is nearly the same for both enantiomers. In contrast, the Lennard–Jones term is similar in magnitude for the different solvents but presents greater differences between the enantiomers. The sharpened area on the curve representing the LJ term inside the cavity is due to the onset of hydrogen bonding between the host and the guest (easier to see in e = 80). This term is highly selective in position and direction and is expected to freeze the centre of mass and orientation of the guest molecule, thus contributing to the discrimination between enantiomers.14 The minimum value of the dielectric constant that allows inclusion complex formation for non-polar solvents is 2 (Fig. 1d). The ELE term inside the cavity is positive in this case, but the contribution of LJ and H-bond terms makes the energy Einter attractive inside the cavity. The electrostatic potential in this case is the main contributor to enantiodifferentiation; its greatest value and the minimum LJ potential occur at positions of the guest near the narrow rim. The external positions of the guest are more stable (less energy) than the internal position, the absolute minimum of Einter was located outside the cavity due to the formation of H-bonds between Ala and b-CD. The ELE term is also positive inside the cavity for the NP configuration of Ala and e = 21 (Fig. 1e), but it is similar for both enantiomers and amounts to approximately 10% of the Einter energy. The minimum interaction energy Emin with the different contributions is given in Table 1, and the complex configurations for each case are represented in Figure 2. The bond angle terms are the main contribution to Eintra, while the O–C–CH–N bond of Ala is the only torsion that contributes to the torsional potential of the guest. The intramolecular energy Eintra is greater for zwitterions than the NP structures of Ala in the absolute minimum, and the solvent affects the conformational adaptation of the guest more than the host as a consequence of their different sizes. The values of Eintra calculated by the GLYCAM force field would be nearly the same because the parameters of the bonds and atoms involved in the system under study are similar to the AMBER parameters. There are some small differences in the LJ potential parameters, but the influence of the force field used on the results is mainly due to the H-bond contribution to the interaction energy. The greatest differences in Einter between the enantiomers in the minima are due to the LJ potential except for e = 2 (NP structure) (as we have seen in Figure 1). The enantiomers with a P configuration are located near the cavity centre in these minima, with the NHþ 3 group pointing to the wide rim of the CD. The complex configurations for each enantiomer in vacuo and with e = 26 (polar solvent) are very similar despite the significant differences in the contributions of Einter. The minimum is outside the CD, near the narrow rim for non-polar solvents with e = 2. Figure 3 represents the projections in the XY and XZ planes of the regions with maximum chiral discrimination of the enantiomers with polar and non-polar solvents. The positions of potential energy surfaces in which the L-Ala is more stable are represented
1201
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
(a)
(b)
(c)
(d)
(e)
Figure 1. The penetration potential W, Lennard–Jones (LJ) and electrostatic (ELE) terms for the interaction between b-CD and Ala (a) in vacuo e = 1 (P structure), and in different solvents: (b) e = 26 (P structure) (c) e = 80 (P structure) (d) e = 2 (NP structure) (e) e = 21 (NP structure).
by crosses and those for D-Ala by circles; the more intense the symbol the greater the difference in energy it represents. Figure 3a shows the main features of the regions with maximum chiral discrimination in vacuo. It can be seen that the L-enantiomer has more positions inside the cavity where it is more stable than the D-enantiomer, with the greatest differences being located near the cavity walls. The characteristics of regions with maximum chiral discrimination for polar solvents with e = 26 (Fig. 3b), water and e = 21 (NP structure) are similar to those in vacuo, although the differences between the enantiomers are smaller. This similarity can be justified, despite the influence of the solvents and the guest configurations, since the electrostatic term is very similar for L- and D-Ala and the differences in Etotal are mainly due to the nonbonded term. On the other hand, most of the positions of D-Ala inside and outside the cavity are more stable for e = 2 (NP structure) (Fig. 3c). The greatest difference in the absolute minimum potential energy (Table 1) and penetration potential W (Fig. 1d) between the enantiomers corresponds to this case, although when we consider the average Boltzmann energy instead of the lowest energy at each grid point, these differences decrease considerably (the maximum value being about 3 kcal/mol) and can even change the sign.
3.2. Molecular dynamics The movements of the guest molecule in each trajectory during the simulation time are different because the initial conditions determine the integration of the equations of motion. However, these initial conditions affect the simulation in different ways: while 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 also determines the behaviour of the amino acid in each process. The simulation time ts for each trajectory is 3 ns but Ala spends this time inside the cavity only in vacuo (e = 1), in the remaining cases it tends to enter the cavity, remain inside to form an inclusion complex and exits from it after a period (residence time t). This effect is illustrated in Figure 4, which shows the projections in the XY and XZ planes of the position probability density for the inclusion of Ala in b-CD, and includes a schematic representation of the projections of the cavity in those planes. The preferential binding site for each enantiomer in vacuo is different (Fig. 4a), and during the simulation time the enantiomers reach stable configurations with small variations in the molecular orientation and
1202
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
Figure 2. The complex configurations of minimum interaction energies Emin for e = 1 (P structure), e = 26 (P structure), e = 80 (P structure), e = 2 (NP structure) and e = 21 (NP structure).
centre of mass position around the most probable configuration (approximately 70%). Although these are not the configurations of absolute minimum energy for L- and D-Ala they are located in favourably enantioselective positions, as shown in Figure 3a. If the amino acid stays temporarily confined inside the CD, it does not pass through the cavity in every process. It can enter and exit from the same rim of the CD, which means greater delocalization for the guest and thus, smaller position probability density distributed over the whole cavity. Figure 4b represents the preferential
binding site for alanine enantiomers in a polar solvent with e = 80, which has a maximum value (20%) greater than the rest of the polar and non-polar solvents considered herein. In these cases there are small differences in the probability of the presence for each enantiomer while the guest is moving inside the cavity, modifying its centre of mass position and orientation continuously. In order to compare the residence time spent by the guest in different solvents, we averaged the values of this magnitude obtained for the trajectories (tmean) (Table 2). It can be seen that tmean is consid-
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
1203
(a)
(b)
(c)
Figure 3. The projections in the XY and XZ planes of the regions with maximum chiral discrimination of Ala enantiomers for (a) e = 1 (P structure) (b) e = 26 (P structure) (c) e = 2 (NP structure). The positions of the potential energy surfaces in which the L-Ala is more stable are represented by crosses, those for D-Ala by circles and the more intense the symbol the greater difference in energy it represents. A schematic representation of the projections of b-CD in those planes has been included.
erably greater for P configurations and decreases as the dielectric constant increases, but the effect is opposite in the NP configuration. In order to study the solvent effects on the interaction energy in MD simulation, we averaged the values of Etotal obtained for t during the trajectories (Emean) and their different terms (Table 2). Emean
is greater than Emin because the energy of every position with a probability different from zero contributes to the average energy. This implies positive values of Emean for the NP configuration of Ala since Eintra >|Einter| in the major part of positions inside the cavity. We appreciate that the intramolecular contribution Eintra for NP is smaller than for zwitterions, but the main difference between
1204
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
(a)
(b)
Figure 4. The projections in the XY and XZ planes of the position probability density for the inclusion of Ala in b-CD for (a) e = 1 (P structure) (b) e = 80 (P structure). A schematic representation of the projections of b-CD in those planes has been included.
1205
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
Table 2 The average total energy Emean, the intramolecular Eintra and intermolecular Einter contributions, the average residence time for the trajectories tmean, the binding free energy for each enantiomer F and the elution order obtained in the molecular dynamics simulation Con.
P P P NP NP
e 1 26 80 2 21
Emean (kcal/mol)
Eintra (kcal/mol)
Einter (kcal/mol)
tmean (ps)
F (kcal/mol)
L
D
L
D
L
D
L
D
L
D
47.68 1.36 0.34 5.12 0.03
49.72 1.30 0.13 6.30 0.21
8.42 8.44 8.44 6.53 6.54
8.49 8.54 8.54 6.65 6.63
56.10 9.80 8.78 1.41 6.57
58.21 9.84 8.67 0.35 6.42
ts 1142.42 939.72 5.40 242.53
ts 1611.58 701.02 8.75 123.77
55.29 7.07 6.17 1.71 4.95
56.72 7.31 5.70 1.53 4.44
Elution order
L L D L D
ts is the simulation time for each trajectory (3 ns).
the two types of configurations is that Eintra amounts to between 50% and 80% in Emean for NP and from 14% to 50% for zwitterions. The ELE term depends on e and the guest configuration (P or NP) according to the interaction energy between b-CD and Ala, but this contribution is very similar for both enantiomers except for e = 2 (NP structure). The H-bond term only contributes during the trajectories in those positions of the guest compatible with the restrictions of distance and orientation of this type of potential; its maximum contribution is about 0.5 kcal/mol. The greatest differences in the mean energy Emean between L- and D-Ala in the simulation are due to the different configurations adopted by the guest moving inside the cavity, and so they are directly related to the Lennard–Jones potential. Therefore, we can conclude that the adaptation of the guest to the host is a decisive factor that influences the chiral discrimination of Ala by b-CD in the presence of non-polar agents, whereas the nonbonded interaction is the driving force for zwitterions. We determined the elution order of the separation of Ala with b-CD in the presence of different solvents by calculating DF = FL FD, where F is the binding free energy in the simulation for each enantiomer (Table 2). It can be seen that F depends on the dielectric constant and also on the alanine configuration, because similar values of e (21 and 26) give rise to opposite signs of DF for P or NP configurations. The first eluted enantiomer in vacuo is the L-Ala in accordance with the experimental findings,9 but this is only possible with the P configuration of enantiomers. Ala modifies the sign of the elution order increasing the value of e, but since zwitterions need a polar solvent such as water to make this possible, a value of e = 21 with the NP configuration has opposite sign to e = 2. These results are in accordance with the experimental separations of amino acids by capillary electrophoresis using cyclodextrins as chiral selectors, performed by Wan et al.10 who observed a reversal of the enantiomeric elution order of some amino acids in the presence of organic modifiers. The first eluted enantiomer has a smaller residence time for each solvent, meaning that the more negative F is, the more stable is the complex formed, while the guest tends to spend more time inside the cavity. The small difference in free energy and residence time obtained for the enantiomers can be due to the AMBER force field parameters that underestimate the H-bond contribution to the interaction energy for these type of complexes, as demonstrated by Paton et al.32 The H-bond term obtained for the interaction between alanine and cellulose was considerably greater, despite the less favourable structure for complexation presented by the amorphous cellulose.14 In that study, the H-bond term was represented by the Lippincott–Schroeder instead of the 10,12 Lennard–Jones potential, but the magnitude of this contribution is not as important as its effect on the enantioselectivity, because although achiral, this term is highly selective in the position and direction and is expected to freeze the centre of mass and orientation of the guest molecule, so influencing the rest of the energy potentials and contributing to discrimination between enantiomers. Nevertheless, the discrimina-
tory ratio obtained for the alanine–CD complex in vacuo (2%) is greater than for alanine–amorphous cellulose (1%).
4. Conclusion We have studied the influence of organic modifiers on the separation of alanine enantiomers using b-CD as a chiral selector, by means of a molecular dynamics simulation at a constant temperature. The potential energy of the interaction is modelled by the AMBER force field, where we represent different polar and non-polar solvents by the dielectric constant e and two different configurations for Ala (NP and P). The interaction energy between b-CD and Ala for different values of e is negative within the cavity for any of the guest configurations, however in vacuo this potential energy is attractive only for the zwitterion configuration. The L-enantiomer has more positions inside the cavity where it is more stable than the D-enantiomer, except for solvents with e = 2 (NP structure) in which the majority of positions of D-Ala inside and outside the cavity are more stable. In all cases, the greatest differences are located near the cavity walls. We concluded from the MD simulation that Ala is able to form inclusion complexes with bCD in vacuo and in the presence of solvents such as hydrocarbons, benzene, acetone, ethanol or water. The chiral discrimination of Ala by b-CD is mainly due to the adaptation of the guest to the host in the presence of non-polar agents, whereas the nonbonded interaction is the driving force for zwitterions. The first eluted enantiomer in vacuo is L-Ala, although it is obtained only with the P structure of guests. The elution order depends on the type of organic modifier and a reversal of enantiomeric elution order can be observed in solvents with higher dielectric constants.
References 1. Szejtli, J.; Osa, T. Comprehensive Supramolecular Chemistry; Pergamon Press: Oxford, 1996. 2. Dodziuk, H. In Cyclodextrins and their Complexes, Chemistry, Analytical Methods, Applications; Wiley-VCH: Weinheim, 2006. 3. Stella, V. J.; Rao, V. M.; Zannou, E. A.; Zia, V. Adv. Drug Deliv. Rev. 1999, 36, 3–16. 4. Schneiderman, E.; Stalcup, A. M. J. Chromatogr., B 2000, 745, 83–102. 5. Lipkowitz, K. B. J. Chromatogr., A 1995, 694, 15–37. 6. Lipkowitz, K. B. Chem. Rev. 1998, 98, 1829–1873. 7. Cunniff, J. B.; Vouros, P. J. Am. Soc. Mass Spectrom. 1995, 6, 437–447. 8. Ramanathan, R.; Prokai, L. J. Am. Soc. Mass Spectrom. 1995, 6, 866–871. 9. Ramirez, J.; Ahn, S.; Grigorean, G.; Lebrilla, C. B. J. Am. Chem. Soc. 2000, 122, 6884–6890. 10. Wan, H.; Engström, A.; Blomberg, L. G. J. Chromatogr., A 1996, 731, 283–292. 11. Elbashir, A.; Suliman, F. O. J. Chromatogr., A 2011, 1218, 5344–5351. 12. Choi, Y.; Jung, S. Carbohydr. Res. 2004, 339, 1961–1966. 13. Alvira, E.; Vega, L.; Girardet, C. Chem. Phys. 1987, 118, 233–240. 14. Alvira, E.; Breton, J.; Plata, J.; Girardet, C. Chem. Phys. 1991, 155, 7–18. 15. Wainer, I. W. TrAC, Trends Anal. Chem. 1987, 6, 125–134. 16. Cezard, C.; Trivelli, X.; Aubry, F.; Djedaini-Pilard, F.; Dupradeau, F. Y. Phys. Chem. Chem. Phys. 2011, 13, 15103–15121. 17. Kirschner, K. N.; Yongye, A. B.; Tschampel, S. M.; JORGE González-Outeiriño, J.; Daniels, C. R.; Foley, B. L.; Woods, R. J. J. Comput. Chem. 2008, 29, 622–655. 18. Lipkowitz, K. B. J. Chromatogr., A 2001, 906, 417–442.
1206
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
19. Lipkowitz, K. B.; Pearl, G.; Coner, B.; Peterson, M. A. J. Am. Chem. Soc. 1997, 119, 600–610. 20. Lipkowitz, K. B.; Coner, B.; Peterson, M. A. J. Am. Chem. Soc. 1997, 119, 11269– 11276. 21. Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S., Jr.; Weiner, P. J. Am. Chem. Soc. 1984, 106, 765–784. 22. Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J. Am. Chem. Soc. 1986, 108, 230–240. 23. Klingert, B.; Rihs, G. J. Chem. Soc., Dalton Trans. 1991, 10, 2749–2760. 24. Girardet, C.; Vega, L. Surf. Sci. 1985, 151, 447–458. 25. Alvira, E.; Mayoral, J. A.; Garcia, J. I. J. Inclusion Phenom. Macrocycl. Chem. 2008, 60, 103–113.
26. 27. 28. 29.
Alvira, E. Supramol. Chem. 2010, 22, 156–162. Brown, D.; Clarke, J. H. R. Mol. Phys. 1984, 51, 1243–1248. Fincham, D.; Quirke, N.; Tildesley, D. J. J. Chem. Phys. 1986, 84, 4535–4541. Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: San Diego, 2002. 30. Rapaport, D. C. The Art of Molecular Dynamics Simulation; Cambridge University Press: Cambridge, 1995. 31. Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. 32. Paton, R. S.; Goodman, J. M. J. Chem. Inf. Model. 2009, 49, 944–955.
Contents lists available at ScienceDirect
Tetrahedron: Asymmetry journal homepage: www.elsevier.com/locate/tetasy
Molecular dynamics study of the influence of solvents on the chiral discrimination of alanine enantiomers by b-cyclodextrin Elena Alvira ⇑ Departamento de Física Fundamental II, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain
a r t i c l e
i n f o
Article history: Received 20 June 2013 Accepted 13 August 2013
a b s t r a c t The influence of solvents on the separation of alanine enantiomers using b-cyclodextrin as a chiral selector was studied by means of a molecular dynamics simulation at a constant temperature. The potential energy of the interaction is modelled by the AMBER force field, where different polar and non-polar solvents are represented by the dielectric constant e and two configurations for the amino acid derived from its electric charge distribution: the AMBER data base or its zwitterion state. The L enantiomer has more positions inside the cavity of a b-cyclodextrin where it is more stable than the D-enantiomer in vacuo and solution, except for solvents such as hydrocarbons in which most positions of the D-alanine inside and outside the cavity are more stable. In all cases, the greatest differences are located near the cavity walls. Molecular dynamics simulations show that Ala is able to form inclusion complexes with b-cyclodextrin in vacuo and in solvents such as hydrocarbons, benzene, acetone, ethanol or water. The chiral discrimination of Ala by b-cyclodextrin is mainly due to the adaptation of the guest to the host in the presence of non-polar agents, whereas the nonbonded interaction is the driving force for zwitterions. The elution order depends on the type of organic modifiers while a reversal of the enantiomeric elution order can be observed in solvents with higher dielectric constants. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction Cyclodextrins (CDs) are macrocyclic molecules composed of glucose units (6 for a-CD, 7 for b CD and 8 for c-CD, etc.) that form truncated cone-shaped compounds. These have cavities of different internal diameters capable of including molecules of different structure, size and composition. Their capacity for catalysis and chiral recognition is mainly due to the formation of inclusion complexes, when some lipophilic part of a molecule enters the hydrophobic cyclodextrin cavity.1,2 Such capabilities have been applied to various research fields such as solubility enhancement, drug delivery, chemical protection, separation technology and supramolecular chemistry.3,4 High-performance liquid chromatography (HPLC) and capillary electrophoresis (CE) using CDs as chiral selectors are frequently used for the enantioseparation of chiral compounds. Amino acids are biologically important organic compounds, which perform critical roles as neurotransmitters, and in transport and synthesis. In the form of proteins they comprise the second largest component (other than water) of human muscles, cells and other tissues. The L- and D-forms of amino acids differently adsorbed on a chiral surface constitute a large group of substances, ⇑ Tel.: +34 922318258; fax: +34 922318320. E-mail address: [email protected] 0957-4166/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.tetasy.2013.08.006
which are separable by CDs. Most of the experimental and theoretical studies carried out with CDs deal with amino acids such as Val, Leu, Ile, Tyr, Phe, Trp and mixtures of them, and aromatic and other types of amino acid derivatives. The reason for this selection of potential guest molecules is presumably the known complexing affinity of CDs for molecules of appropriate sizes and for aromatic functional groups in aqueous solutions.5,6 In fact, there is experimental evidence of inclusion complex formation by these amino acids with CDs in solution,7,8 as well as in gas phase.9 The amino acid Ala is smaller than the aforementioned amino acids but Ramirez et al. also demonstrated the inclusion complex formation of Ala with CDs in the gas phase.9 Experimental results also reveal the influence of organic modifiers on the enantiomeric elution order for some amino acids using CDs as chiral selectors.10 Theoretical methods such as molecular mechanics11 and molecular dynamics simulations9,12 are used to help explain the enantiodifferentiation mechanism and predict the separation process. In these theoretical studies, the presence of organic modifiers is represented by different values of the dielectric constant in the electrostatic contribution to the interaction energy. However under some pH conditions or with polar solvents, amino acids appear as zwitterions with the NHþ 3 and COO groups instead of NH2 and COOH, and this implies not only a different electric charge distribution, but also a different molecular configuration. Therefore, we can represent polar solvents such as ethanol and non-polar solvents
1199
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
such as hydrocarbons not only by different values of dielectric constant, but also by different configurations of the amino acids. This allows us to study the influence of solvent polarity on the separation of Ala by b-CD. In previous work, we analysed the chiral discrimination of alanine enantiomers with cellulose (crystalline and amorphous) using molecular mechanics simulations.13,14 The amorphous cellulose exhibits spatial helicoidal arrangements of the glucose units, allowing partial inclusion of the analyte molecules. We concluded that chiral discrimination increases when Ala is located inside the cavity, and it is maximum for similar sizes of helix radii and pitch. Cyclodextrins and cellulose derivatives are usually considered different types of chiral stationary phases (CSP)15,16 although they are composed of glucose rings (differently linked). The line of demarcation between them (suggested by Wainer and Lipkowitz) is the extent of guest inclusion into a chiral cavity. Herein our aim was to study theoretically inclusion complex formation and separation process of Ala with b-CD, in the presence of different solvents. We determined the influence of organic modifiers on the inclusion complex formation and elution order of the separation of alanine enantiomers by b-CD using molecular dynamics (MD). The potential energy of the interaction is modelled by the AMBER force field, although in recent years there has been an increasing development in new force fields that are applicable to carbohydrates, such as GLYCAM17 or q4md-CD, built from a combination of the GLYCAM04 and Amber99SB force fields to correctly describe the geometrical, structural, dynamical and hydrogen bonding aspects of heterogeneous cyclodextrin based systems.18 In the molecular simulation methods, the elution order in an isomeric separation is usually determined from the differences in free energy. The smallest value means a more tightly bound analyte and this is linked with the longest retention time. The chromatographic retention time is an experimental result that depends on macroscopic parameters; in contrast, free energy is calculated theoretically from microscopic modelling on the molecular scale; therefore it is only possible to compare the two sets of results qualitatively. However, MD also allows us to determine the time spent by the guest inside the cavity and provides a huge amount of statistical information to calculate where and how selective binding takes place.19,20 These results are not usually provided by the theoretical studies of separations of amino acids by CDs, but they can help to optimise these processes. The interaction potential and simulation method used herein are presented in Section 2. Section 3 is devoted to the analysis of the interaction energies between b-CD and the alanine enantiomers, and a discussion of the main results of molecular dynamics simulations related to inclusion complex formation and elution order.
2. The Model 2.1. Expression of the interaction potential In order to determine the interaction energy between b-CD and Ala, we used the force field proposed by Weiner et al. for the molecular mechanics simulation of nucleic acids and proteins.21,22 The interaction potential Etotal is the sum of the intramolecular Eintra and intermolecular Einter energies:
Etotal ¼ Eintra þ Einter
ð1Þ
The intramolecular energy is modelled by the equation:
Eintra ¼
X
kr ðr req Þ2 þ
bonds
X angles
kh ðh heq Þ2 þ
X Vn ½1 þ cosðn/ cÞ 2 dihedrals ð2Þ
The bond stretching and bending functions are quadratic, which allows an adequate description of the structure and energies for proteins and nucleic acids. The torsional energy is represented by a Fourier series approach. The intermolecular energy includes the nonbonded (Lennard–Jones), electrostatic and H-bond terms:
Einter ¼
" X Aij
Bij
Rij
Rij
12
i
þ 6
# " # X C ij qi qj Dij þ 12 eRij R10 Hbonds Rij ij
ð3Þ
where Rij represents the distance between the ith atom of the guest and the jth atom of b-CD. The atomic coordinates of b-CD are taken from the literature,23 and its net atomic charges and the AMBER force field parameters (Aij, Bij, Cij and Dij) are employed.21,22 We also considered the AMBER data as the molecular configuration and atomic point charges for the amino acid,21 but under some pH conditions or with polar solvents, the amino acid appears as a zwitterion (with the NHþ groups 3 and COO instead of NH2 and COOH). This presents different atomic point charges24 and molecular configurations from the AMBER force field. We represent P as the zwitterion configuration, NP as the non-polar configuration of Ala and the solvents by different values of the dielectric constant (Table 1). We considered all of the atoms in the molecule of alanine because the difference between the two configurations is related to the positions of some H atoms that can contribute decisively to the formation of hydrogen bonding between host and guest, and may be reflected in the interaction energy Etotal. We placed the origin of the reference system at the centre of the mass of the cavity and the space-fixed frame over the principal axis of the CD, in which the inertia tensor is diagonal. The position of Ala is given by the coordinates of its centre of mass, and the molecular orientation is defined by the relation between its principal body-fixed system and the axis system fixed in space. The potential energy for each enantiomer is determined by Eq. (1) for different positions and orientations of the guest centre of the mass, inside and outside the CD. In each plane Z = constant, about 2500 points are explored, the range of variation along the Z axis is 10 Å (with a path of 0.1 Å) and about 23,000 orientations of the guest are considered at each grid point. The results obtained are represented as potential energy surfaces, penetration potentials and the inclusion complex configuration, as in previous work.25,26 The curve joining the minimum intermolecular energy for every plane Z = constant defines the penetration potential W, which describes the variation in Einter when its path through the cavity is non-axial. The position of the guest molecule for which we obtain the absolute minimum potential energy Etotal gives the geometry of the inclusion complex for the different solvents. However the enantiomers do not generally adopt the lowest energy configuration due to the separation process, and can form complexes with configurations different from the minima that lead to a different enantiomeric elution order. We considered the regions with maximum chiral discrimination to be those localised in the positions of potential surfaces with greatest differences in energy, but we assigned to each grid point the average Boltzmann energy corresponding to different guest orientations, instead of the lowest energy.19,20 2.2. Simulation method The molecular dynamics (MD) is based on the resolution of classical equations of motion to determine the trajectories of the particles, depending on the initial conditions of the guest molecule: position, orientation (Euler angles) and velocities (translational and rotational). The magnitude of the initial velocities depends on the temperature of the process (293 K) but the directions of the translational and rotational velocities in each trajectory, as well as the initial centre of mass position are determined randomly. The
1200
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
Table 1 The minimum interaction energy Emin obtained with the AMBER force field for each enantiomer and the different contributions: Lennard–Jones ELJ, electrostatic Eele, hydrogen bonding EH-bond, bond stretching Ebond, angle bending Eangle and torsional energy Etorsion Con.
P P P NP NP
e 1 26 80 2 21
Emin (kcal/mol)
ELJ (kcal/mol)
Eele (kcal/mol)
EH-bond (kcal/mol)
Ebond (kcal/mol)
Eangle (kcal/mol)
Etorsion (kcal/mol)
L
D
L
D
L
D
L
D
L
D
L
D
L
D
52.95 4.67 3.47 5.24 4.00
51.99 4.51 3.28 3.71 3.76
10.11 10.71 10.51 2.91 11.23
8.97 10.67 10.79 3.20 10.63
50.77 1.89 0.55 8.31 0.88
51.00 1.83 0.53 7.26 0.85
0.48 0.49 0.91 0.50 0.00
0.49 0.47 0.49 0.00 0.49
1.72 1.74 1.80 5.09 4.97
1.66 1.66 1.74 5.21 5.01
6.61 6.59 6.62 1.30 1.29
6.61 6.61 6.59 1.34 1.30
0.08 0.08 0.08 0.09 0.09
0.20 0.20 0.20 0.20 0.20
simulation time for each trajectory is 3 ns with a step of 1 fs and the configuration and energies (kinetic and potential) were written every 100 steps. We used an in-house programme written in Fortran and the equations of motion in order to ensure constant temperature molecular dynamics. They are integrated numerically, using a variant of the leap-frog scheme (proposed by Brown and Clarke),27 which constrains the rotational and translational kinetic energies separately.28 This programme has been used to study the inclusion complex formation of different types of guests in b-CD.25,26 We calculated 6 different trajectories with initial configurations of the guest outside the b-CD: 3 near the primary (narrow end) and 3 near the secondary rims (wide end) of the CD. In order to avoid the influence of different initial conditions, in each trajectory the alanine enantiomers have the same initial centre of mass position and molecular orientation, so they differ in their atomic positions. In order to analyse the ability of Ala to form b-CD inclusion complexes with different solvents, we determined the preferential binding site of the guest molecule in the simulation by means of the probability density of presence in a volume element. We defined 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 guest.19,20 The probability density of presence is calculated by dividing the number density in a volume element by the total number of centre of mass positions of the guest. We determined the elution order by calculating DF = FL FD, with F being the binding free energy in the simulation for each enantiomer.29–31 In order to compare the time spent by the guest inside the cavity in different solvents with F, we averaged the values of the time obtained in MD simulations.
3. Results and discussion 3.1. Interaction energy The values of dielectric constant e considered herein are similar to those of agents such as hydrocarbons, benzene, toluene, acetone, ethanol or water. The interaction energy Etotal between b-CD and Ala for these values of e is negative inside the cavity with any of the guest configurations (NP or P), which means that the amino acid is able to form inclusion complexes. However, in vacuo (e = 1) the energy Etotal is positive inside the b-CD cavity when the configuration of Ala is NP, in this case the inclusion complex can only be formed by the configuration P of zwitterions. Ramirez et al. assumed a dielectric of 1.0 in the molecular modelling simulations where they provided theoretical evidence of inclusion complex formation of several amino acids in the gas phase,9 but they did not study alanine enantiomers. Figure 1 represents the penetration potential W, Lennard–Jones (LJ) and electrostatic (ELE) terms for the interaction between b-CD and the enantiomers of Ala, with different solvents. In the case of a polar configuration
with e = 1 (Fig. 1a), e = 26 (Fig. 1b) and e = 80 (Fig. 1c), as e increases the electrostatic term contributes less to the total energy, and we can see that this contribution is nearly the same for both enantiomers. In contrast, the Lennard–Jones term is similar in magnitude for the different solvents but presents greater differences between the enantiomers. The sharpened area on the curve representing the LJ term inside the cavity is due to the onset of hydrogen bonding between the host and the guest (easier to see in e = 80). This term is highly selective in position and direction and is expected to freeze the centre of mass and orientation of the guest molecule, thus contributing to the discrimination between enantiomers.14 The minimum value of the dielectric constant that allows inclusion complex formation for non-polar solvents is 2 (Fig. 1d). The ELE term inside the cavity is positive in this case, but the contribution of LJ and H-bond terms makes the energy Einter attractive inside the cavity. The electrostatic potential in this case is the main contributor to enantiodifferentiation; its greatest value and the minimum LJ potential occur at positions of the guest near the narrow rim. The external positions of the guest are more stable (less energy) than the internal position, the absolute minimum of Einter was located outside the cavity due to the formation of H-bonds between Ala and b-CD. The ELE term is also positive inside the cavity for the NP configuration of Ala and e = 21 (Fig. 1e), but it is similar for both enantiomers and amounts to approximately 10% of the Einter energy. The minimum interaction energy Emin with the different contributions is given in Table 1, and the complex configurations for each case are represented in Figure 2. The bond angle terms are the main contribution to Eintra, while the O–C–CH–N bond of Ala is the only torsion that contributes to the torsional potential of the guest. The intramolecular energy Eintra is greater for zwitterions than the NP structures of Ala in the absolute minimum, and the solvent affects the conformational adaptation of the guest more than the host as a consequence of their different sizes. The values of Eintra calculated by the GLYCAM force field would be nearly the same because the parameters of the bonds and atoms involved in the system under study are similar to the AMBER parameters. There are some small differences in the LJ potential parameters, but the influence of the force field used on the results is mainly due to the H-bond contribution to the interaction energy. The greatest differences in Einter between the enantiomers in the minima are due to the LJ potential except for e = 2 (NP structure) (as we have seen in Figure 1). The enantiomers with a P configuration are located near the cavity centre in these minima, with the NHþ 3 group pointing to the wide rim of the CD. The complex configurations for each enantiomer in vacuo and with e = 26 (polar solvent) are very similar despite the significant differences in the contributions of Einter. The minimum is outside the CD, near the narrow rim for non-polar solvents with e = 2. Figure 3 represents the projections in the XY and XZ planes of the regions with maximum chiral discrimination of the enantiomers with polar and non-polar solvents. The positions of potential energy surfaces in which the L-Ala is more stable are represented
1201
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
(a)
(b)
(c)
(d)
(e)
Figure 1. The penetration potential W, Lennard–Jones (LJ) and electrostatic (ELE) terms for the interaction between b-CD and Ala (a) in vacuo e = 1 (P structure), and in different solvents: (b) e = 26 (P structure) (c) e = 80 (P structure) (d) e = 2 (NP structure) (e) e = 21 (NP structure).
by crosses and those for D-Ala by circles; the more intense the symbol the greater the difference in energy it represents. Figure 3a shows the main features of the regions with maximum chiral discrimination in vacuo. It can be seen that the L-enantiomer has more positions inside the cavity where it is more stable than the D-enantiomer, with the greatest differences being located near the cavity walls. The characteristics of regions with maximum chiral discrimination for polar solvents with e = 26 (Fig. 3b), water and e = 21 (NP structure) are similar to those in vacuo, although the differences between the enantiomers are smaller. This similarity can be justified, despite the influence of the solvents and the guest configurations, since the electrostatic term is very similar for L- and D-Ala and the differences in Etotal are mainly due to the nonbonded term. On the other hand, most of the positions of D-Ala inside and outside the cavity are more stable for e = 2 (NP structure) (Fig. 3c). The greatest difference in the absolute minimum potential energy (Table 1) and penetration potential W (Fig. 1d) between the enantiomers corresponds to this case, although when we consider the average Boltzmann energy instead of the lowest energy at each grid point, these differences decrease considerably (the maximum value being about 3 kcal/mol) and can even change the sign.
3.2. Molecular dynamics The movements of the guest molecule in each trajectory during the simulation time are different because the initial conditions determine the integration of the equations of motion. However, these initial conditions affect the simulation in different ways: while 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 also determines the behaviour of the amino acid in each process. The simulation time ts for each trajectory is 3 ns but Ala spends this time inside the cavity only in vacuo (e = 1), in the remaining cases it tends to enter the cavity, remain inside to form an inclusion complex and exits from it after a period (residence time t). This effect is illustrated in Figure 4, which shows the projections in the XY and XZ planes of the position probability density for the inclusion of Ala in b-CD, and includes a schematic representation of the projections of the cavity in those planes. The preferential binding site for each enantiomer in vacuo is different (Fig. 4a), and during the simulation time the enantiomers reach stable configurations with small variations in the molecular orientation and
1202
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
Figure 2. The complex configurations of minimum interaction energies Emin for e = 1 (P structure), e = 26 (P structure), e = 80 (P structure), e = 2 (NP structure) and e = 21 (NP structure).
centre of mass position around the most probable configuration (approximately 70%). Although these are not the configurations of absolute minimum energy for L- and D-Ala they are located in favourably enantioselective positions, as shown in Figure 3a. If the amino acid stays temporarily confined inside the CD, it does not pass through the cavity in every process. It can enter and exit from the same rim of the CD, which means greater delocalization for the guest and thus, smaller position probability density distributed over the whole cavity. Figure 4b represents the preferential
binding site for alanine enantiomers in a polar solvent with e = 80, which has a maximum value (20%) greater than the rest of the polar and non-polar solvents considered herein. In these cases there are small differences in the probability of the presence for each enantiomer while the guest is moving inside the cavity, modifying its centre of mass position and orientation continuously. In order to compare the residence time spent by the guest in different solvents, we averaged the values of this magnitude obtained for the trajectories (tmean) (Table 2). It can be seen that tmean is consid-
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
1203
(a)
(b)
(c)
Figure 3. The projections in the XY and XZ planes of the regions with maximum chiral discrimination of Ala enantiomers for (a) e = 1 (P structure) (b) e = 26 (P structure) (c) e = 2 (NP structure). The positions of the potential energy surfaces in which the L-Ala is more stable are represented by crosses, those for D-Ala by circles and the more intense the symbol the greater difference in energy it represents. A schematic representation of the projections of b-CD in those planes has been included.
erably greater for P configurations and decreases as the dielectric constant increases, but the effect is opposite in the NP configuration. In order to study the solvent effects on the interaction energy in MD simulation, we averaged the values of Etotal obtained for t during the trajectories (Emean) and their different terms (Table 2). Emean
is greater than Emin because the energy of every position with a probability different from zero contributes to the average energy. This implies positive values of Emean for the NP configuration of Ala since Eintra >|Einter| in the major part of positions inside the cavity. We appreciate that the intramolecular contribution Eintra for NP is smaller than for zwitterions, but the main difference between
1204
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
(a)
(b)
Figure 4. The projections in the XY and XZ planes of the position probability density for the inclusion of Ala in b-CD for (a) e = 1 (P structure) (b) e = 80 (P structure). A schematic representation of the projections of b-CD in those planes has been included.
1205
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
Table 2 The average total energy Emean, the intramolecular Eintra and intermolecular Einter contributions, the average residence time for the trajectories tmean, the binding free energy for each enantiomer F and the elution order obtained in the molecular dynamics simulation Con.
P P P NP NP
e 1 26 80 2 21
Emean (kcal/mol)
Eintra (kcal/mol)
Einter (kcal/mol)
tmean (ps)
F (kcal/mol)
L
D
L
D
L
D
L
D
L
D
47.68 1.36 0.34 5.12 0.03
49.72 1.30 0.13 6.30 0.21
8.42 8.44 8.44 6.53 6.54
8.49 8.54 8.54 6.65 6.63
56.10 9.80 8.78 1.41 6.57
58.21 9.84 8.67 0.35 6.42
ts 1142.42 939.72 5.40 242.53
ts 1611.58 701.02 8.75 123.77
55.29 7.07 6.17 1.71 4.95
56.72 7.31 5.70 1.53 4.44
Elution order
L L D L D
ts is the simulation time for each trajectory (3 ns).
the two types of configurations is that Eintra amounts to between 50% and 80% in Emean for NP and from 14% to 50% for zwitterions. The ELE term depends on e and the guest configuration (P or NP) according to the interaction energy between b-CD and Ala, but this contribution is very similar for both enantiomers except for e = 2 (NP structure). The H-bond term only contributes during the trajectories in those positions of the guest compatible with the restrictions of distance and orientation of this type of potential; its maximum contribution is about 0.5 kcal/mol. The greatest differences in the mean energy Emean between L- and D-Ala in the simulation are due to the different configurations adopted by the guest moving inside the cavity, and so they are directly related to the Lennard–Jones potential. Therefore, we can conclude that the adaptation of the guest to the host is a decisive factor that influences the chiral discrimination of Ala by b-CD in the presence of non-polar agents, whereas the nonbonded interaction is the driving force for zwitterions. We determined the elution order of the separation of Ala with b-CD in the presence of different solvents by calculating DF = FL FD, where F is the binding free energy in the simulation for each enantiomer (Table 2). It can be seen that F depends on the dielectric constant and also on the alanine configuration, because similar values of e (21 and 26) give rise to opposite signs of DF for P or NP configurations. The first eluted enantiomer in vacuo is the L-Ala in accordance with the experimental findings,9 but this is only possible with the P configuration of enantiomers. Ala modifies the sign of the elution order increasing the value of e, but since zwitterions need a polar solvent such as water to make this possible, a value of e = 21 with the NP configuration has opposite sign to e = 2. These results are in accordance with the experimental separations of amino acids by capillary electrophoresis using cyclodextrins as chiral selectors, performed by Wan et al.10 who observed a reversal of the enantiomeric elution order of some amino acids in the presence of organic modifiers. The first eluted enantiomer has a smaller residence time for each solvent, meaning that the more negative F is, the more stable is the complex formed, while the guest tends to spend more time inside the cavity. The small difference in free energy and residence time obtained for the enantiomers can be due to the AMBER force field parameters that underestimate the H-bond contribution to the interaction energy for these type of complexes, as demonstrated by Paton et al.32 The H-bond term obtained for the interaction between alanine and cellulose was considerably greater, despite the less favourable structure for complexation presented by the amorphous cellulose.14 In that study, the H-bond term was represented by the Lippincott–Schroeder instead of the 10,12 Lennard–Jones potential, but the magnitude of this contribution is not as important as its effect on the enantioselectivity, because although achiral, this term is highly selective in the position and direction and is expected to freeze the centre of mass and orientation of the guest molecule, so influencing the rest of the energy potentials and contributing to discrimination between enantiomers. Nevertheless, the discrimina-
tory ratio obtained for the alanine–CD complex in vacuo (2%) is greater than for alanine–amorphous cellulose (1%).
4. Conclusion We have studied the influence of organic modifiers on the separation of alanine enantiomers using b-CD as a chiral selector, by means of a molecular dynamics simulation at a constant temperature. The potential energy of the interaction is modelled by the AMBER force field, where we represent different polar and non-polar solvents by the dielectric constant e and two different configurations for Ala (NP and P). The interaction energy between b-CD and Ala for different values of e is negative within the cavity for any of the guest configurations, however in vacuo this potential energy is attractive only for the zwitterion configuration. The L-enantiomer has more positions inside the cavity where it is more stable than the D-enantiomer, except for solvents with e = 2 (NP structure) in which the majority of positions of D-Ala inside and outside the cavity are more stable. In all cases, the greatest differences are located near the cavity walls. We concluded from the MD simulation that Ala is able to form inclusion complexes with bCD in vacuo and in the presence of solvents such as hydrocarbons, benzene, acetone, ethanol or water. The chiral discrimination of Ala by b-CD is mainly due to the adaptation of the guest to the host in the presence of non-polar agents, whereas the nonbonded interaction is the driving force for zwitterions. The first eluted enantiomer in vacuo is L-Ala, although it is obtained only with the P structure of guests. The elution order depends on the type of organic modifier and a reversal of enantiomeric elution order can be observed in solvents with higher dielectric constants.
References 1. Szejtli, J.; Osa, T. Comprehensive Supramolecular Chemistry; Pergamon Press: Oxford, 1996. 2. Dodziuk, H. In Cyclodextrins and their Complexes, Chemistry, Analytical Methods, Applications; Wiley-VCH: Weinheim, 2006. 3. Stella, V. J.; Rao, V. M.; Zannou, E. A.; Zia, V. Adv. Drug Deliv. Rev. 1999, 36, 3–16. 4. Schneiderman, E.; Stalcup, A. M. J. Chromatogr., B 2000, 745, 83–102. 5. Lipkowitz, K. B. J. Chromatogr., A 1995, 694, 15–37. 6. Lipkowitz, K. B. Chem. Rev. 1998, 98, 1829–1873. 7. Cunniff, J. B.; Vouros, P. J. Am. Soc. Mass Spectrom. 1995, 6, 437–447. 8. Ramanathan, R.; Prokai, L. J. Am. Soc. Mass Spectrom. 1995, 6, 866–871. 9. Ramirez, J.; Ahn, S.; Grigorean, G.; Lebrilla, C. B. J. Am. Chem. Soc. 2000, 122, 6884–6890. 10. Wan, H.; Engström, A.; Blomberg, L. G. J. Chromatogr., A 1996, 731, 283–292. 11. Elbashir, A.; Suliman, F. O. J. Chromatogr., A 2011, 1218, 5344–5351. 12. Choi, Y.; Jung, S. Carbohydr. Res. 2004, 339, 1961–1966. 13. Alvira, E.; Vega, L.; Girardet, C. Chem. Phys. 1987, 118, 233–240. 14. Alvira, E.; Breton, J.; Plata, J.; Girardet, C. Chem. Phys. 1991, 155, 7–18. 15. Wainer, I. W. TrAC, Trends Anal. Chem. 1987, 6, 125–134. 16. Cezard, C.; Trivelli, X.; Aubry, F.; Djedaini-Pilard, F.; Dupradeau, F. Y. Phys. Chem. Chem. Phys. 2011, 13, 15103–15121. 17. Kirschner, K. N.; Yongye, A. B.; Tschampel, S. M.; JORGE González-Outeiriño, J.; Daniels, C. R.; Foley, B. L.; Woods, R. J. J. Comput. Chem. 2008, 29, 622–655. 18. Lipkowitz, K. B. J. Chromatogr., A 2001, 906, 417–442.
1206
E. Alvira / Tetrahedron: Asymmetry 24 (2013) 1198–1206
19. Lipkowitz, K. B.; Pearl, G.; Coner, B.; Peterson, M. A. J. Am. Chem. Soc. 1997, 119, 600–610. 20. Lipkowitz, K. B.; Coner, B.; Peterson, M. A. J. Am. Chem. Soc. 1997, 119, 11269– 11276. 21. Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S., Jr.; Weiner, P. J. Am. Chem. Soc. 1984, 106, 765–784. 22. Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J. Am. Chem. Soc. 1986, 108, 230–240. 23. Klingert, B.; Rihs, G. J. Chem. Soc., Dalton Trans. 1991, 10, 2749–2760. 24. Girardet, C.; Vega, L. Surf. Sci. 1985, 151, 447–458. 25. Alvira, E.; Mayoral, J. A.; Garcia, J. I. J. Inclusion Phenom. Macrocycl. Chem. 2008, 60, 103–113.
26. 27. 28. 29.
Alvira, E. Supramol. Chem. 2010, 22, 156–162. Brown, D.; Clarke, J. H. R. Mol. Phys. 1984, 51, 1243–1248. Fincham, D.; Quirke, N.; Tildesley, D. J. J. Chem. Phys. 1986, 84, 4535–4541. Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: San Diego, 2002. 30. Rapaport, D. C. The Art of Molecular Dynamics Simulation; Cambridge University Press: Cambridge, 1995. 31. Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. 32. Paton, R. S.; Goodman, J. M. J. Chem. Inf. Model. 2009, 49, 944–955.