Thermodynamics of Supramolecular Associations with Macrocyclic Water-Soluble Hosts.

The thermodynamic study of the complexation of the β-cyclodextrins and p-sulfonatocalix[n]arenes (CnS) with the 4-aminoazobenzene was reported and was carried out by molecular dynamics simulations. We determined the whole thermodynamic properties (K, Δr G°, Δr H°, and TΔr S°) using the potential of mean force (PMF) technique and more precisely the adaptive biasing force method. Depending on both the nature of the host molecule and the pH of the solution, the PMF profiles present different shapes and energy minima. Considering the complexity of these PMF profiles, we are also interested in the structural properties of these associations. Hence, we calculated hydrogen bonds, Lennard-Jones and electrostatic energies, the number of atoms of the guest molecule inserted inside the cagelike host molecule, and the number of water molecules expelled from the cavity.

Introduction

Host–guest systems consist generally of two parts, the molecular receptors (usually small molecules with a hydrophobic cavity such as cyclodextrins (CD),[1] calixarenes,[2] or cucurbiturils[3−5]) and the guests that are able to bind with both high affinity and/or selectivity. For technological applications in the field of photoswitchable materials,[6−11] azobenzene derivatives have been widely investigated and cyclodextrins have been largely used as hosts for these derivatives. In the literature, some stability constants are available revealing some discrepancies.[12−18] In fact, supramolecular assemblies typically provide weak force associations that are quite difficult to detect and quantify with both precision and accuracy. Moreover, the experimental conditions are usually so different that the values can hardly be compared.

We previously reported experimental and computational investigations on the association of water-soluble macrocyclic hosts (β-cyclodextrin and p-sulfonatocalixarenes) with 4-aminoazobenzene (4AA)[19] in water (Figure ). UV–visible spectroscopy was used to provide the whole thermodynamic characterizations of the association (ΔrG, ΔrH, and TΔrS). Molecular simulations were also performed to interpret experimental results and to obtain structural data for these complexes. The pH effect was studied. Whereas the thermodynamic parameters found for β-cyclodextrin were almost identical at the two pH values studied (pH 1 and 7), these parameters differed significantly for p-sulfonatocalixarenes. Stable inclusion complexes with calixarenes were found at acidic pH, but our investigations were unable to confirm the inclusion of the guest at neutral pH. Knowing that the binding abilities of such calixarenes are dramatically pH-controlled, a weaker association was indeed expected.[20] However, a clarification is required to provide a definitive answer regarding the possible binding of the 4AA to calixarenes at neutral pH.

Figure 1

Representation of (a) β-Cd, (b) 4AA, and (c) p-sulfonatocalix[4]arene.

Molecular description combined with thermodynamic characterization using numerical experiments is an efficient approach for the study of the binding process.[21−26] In the present paper, we have re-examined the association process using the free energy approach to obtain a thermodynamic characterization using numerical experiments. This approach corresponds to the calculation of the free energy along a reaction coordinate. The calculation of the potential of mean force (PMF) along complex chemical reactions was originally introduced by Kirkwood in 1935.[27] The comparison between the binding properties previously obtained from UV–visible spectroscopy for β-cyclodextrin and p-sulfonatocalixarenes with 4-aminoazobenzene in water at acidic pH allowed the validation of the computational procedures. As a result, we have given a quantitative answer to the possible association between p-sulfonatocalixarenes and 4-aminoazobenzene in water at neutral pH.

Section discusses the main results, and our conclusions are given in Section . Section describes the calculation of the potential of mean force using the adaptive biasing force (ABF) and of the corresponding thermodynamic properties of association. The computational details of the molecular dynamics are also provided.

Results and Discussion

β-Cyclodextrin/4AA

The potential of mean force characterizing the association process of the 4AA with the β-cyclodextrin is shown in Figure at pH = 7. No significant changes are observed by changing the pH in concordance with our experimental results.[19] The Gibbs free energy profile in Figure shows some local minima characterized by separation distances that clearly implies the insertion of the 4AA into the CD cavity in line with our previous results. Because the azobenzene moiety is longer than the cyclodextrin cavity, the 4AA is included longitudinally in the β-cyclodextrin cavity in perfect agreement with some reported NMR results.[15] Interestingly, three quite different local minima were found; the first at the distance of −1.3 Å is the deeper free energy minimum. It corresponds to the configuration in which the first phenyl ring of the azobenzene moiety is inside the cavity with the hydrophilic head, the amino group, exposed to water. The second minimum at 0.75 Å corresponds to the configuration in which the two phenyls rings of azobenzene are exposed to the bulk water. Finally, in the third minimum at 3.95 Å, it is the second phenyl ring that is included in the cavity of β-cyclodextrin.

Figure 2

Gibbs free energy profile obtained for the inclusion complex (β-Cd, 4AA) as a function of the reaction coordinate. In the inset, the PMF is shown in the whole range of sampled separation distances.

To explain the driving forces of β-cyclodextrin complexation, we report in Figure a the energy contributions and in Figure b the percentage of inserted atoms into the CD cavity at different separation distances. The number of water molecules inside the cavity of CD is also detailed along the coordinate reaction in Figure c. As expected for an insertion complex with the β-CD, the van der Waals interaction is predominant over the electrostatic energy contribution. A strong correlation was also observed between the LJ contributions and the number of inserted atoms. This contribution is more favorable between the distance of −4 and 4 Å of the two centers of mass when the percentage of inserted atoms is maximum. Note that the Lennard-Jones (LJ) energy minimum corresponds to the first deeper Gibbs free energy minimum of the PMF at −1.3 Å (Figure ). Regarding the number of water molecules inserted into the cavity, two minimums are shown at −0.7 and 3.3 Å, respectively (Figure b). These separation distances are in line with the two local minima of the potential of mean force corresponding to the locations of the two phenyl rings of azobenzene within the cavity. Clearly, these conformations involve larger dehydration of the host. We also observe that the cavity recovers some water molecules when the two phenyl groups of 4AAB are exposed to water; this situation corresponds to the second minimum at 0.75 Å in the free energy profile (Figure ).

Figure 3

(a) Energy contributions as a function of the distance between centers of mass for (β-Cd, 4AA). The blue curve corresponds to the Lennard-Jones energy, the green to the electrostatic part, and the red is the sum of both contributions. (b) Percentage of 4AA atoms inserted inside the β-Cd cavity as a function of the distance between centers of mass. (c) Number of water molecules inside the host cavity as a function of the distance between the host and the guest.

To conclude on these aspects, these results indicate that the more negative free energy minimum at −1.3 Å in the free energy profile given in Figure is a result of the combination of both the larger dehydration of the cavity and the more negative LJ energy contributions. Finally, for the same system, we have also reported the number of hydrogen bonds formed by the 4AA with water molecules as a function of the separation distance of the two centers of mass as illustrated in Figure . Obviously, as 4AA penetrates into the cavity of the host, the total number of hydrogen bonds formed by the 4AA decreases. As shown in Figure , this decrease is mainly due to the azo group that is deeply inserted into the cavity, preventing any formation of hydrogen bonds.

Figure 4

Number of hydrogen bonds (H-bonds) formed by the 4AA with water molecules. The colored curves correspond to a given type of H-bond, whereas the black curve corresponds to the sum of all contributions.

Calixarenes/4AA

Now, we focus on the potential of mean force characterizing the association process of the 4AA into p-sulfonatocalixarenes in both neutral and acidic solutions (Figure ). It should be noted that compared to the previous system, the 4AA does not pass through the cavity of the calixarenes. The results obtained for C4S and C6S are rather similar in terms of the shape and well depth of the PMF curve. However, the PMFs are strongly dependent on the pH. Indeed, the PMF obtained in acidic solution shows deeper minimum than in neutral medium and the location of 4AA into the cavity is always deeper at pH = 1 (at 4.1 and 2.5 Å for C4S and C6S, respectively) than that at pH = 7 (at 5.35 and 3.35 Å for C4S and C6S, respectively). Interestingly, for both C4S and C6S at pH = 7, one can observe a small energy barrier (about 5 kJ mol–1) that precedes the local free energy minimum.

Figure 5

Free energy profiles obtained for (a) (C4S, 4AA) and (b) (C6S, 4AA). The plain line corresponds to acidic pH, whereas the dotted line refers to the neutral medium.

We propose here to focus our attention on the C4S. To further understand these PMFs, the partitioning of the free energy profiles into the Lennard-Jones and electrostatic parts of the energy contributions is reported as a function of the separation distance in Figure for (C4S, 4AA) in water. The total energy contributions are also reported in the same figure. At both pH values, the variations of the LJ contributions as a function of the distance are similar. When the 4AA approaches the C4S, these contributions become largely favorable, with the insertion of the 4AA into the cavity starting at a distance of about 8 Å. Figure a shows the percentage of inserted atoms into the C4S cavity, and Figure b shows the number of water molecules inside the host cavity. Interestingly, the more negative LJ contribution at the Gibbs free energy minima in acidic medium corresponds to 40% of the guest insertion, whereas it is only 20% in neutral solution. In both cases, at the Gibbs free energy minima, the dehydration of the cavity represents at least 50% of the total number of the water molecules inside the cavity.

Figure 6

Energy contributions as a function of the distance between centers of mass for (C4S, 4AA) at (a) acidic pH and (b) neutral pH. The blue curve corresponds to the Lennard-Jones energy, the green to the electrostatic part, and the red is the sum of all contributions. The dotted lines represent the location of the deep minimum in the PMF profile for (C4S, 4AA).

Figure 7

(a) Percentage of 4AA atoms inserted inside the C4S cavity as a function of the distance between centers of mass. (b) Number of water molecules inside the host cavity as a function of the distance between the host and the guest. The dotted lines represent the location of the deep minimum in the PMF profile for (C4S, 4AA).

According to the PMFs discussed above, it is clear that changing the pH impacts significantly the electrostatic interactions. In acidic medium, these contributions are negative and gradually decrease when the guest approaches the cavity. On the contrary, the electrostatic part is always positive along the host–guest separation distance in the neutral medium. Interestingly, the penetration of 4AA into the cavity (from a distance of 8 Å) is associated with an electrostatic energy cost in both cases. Note that this cost is in the range of 400 kJ mol–1 at pH = 1, whereas it is rather on the order of 40 kJ mol–1 at pH = 7. However, the electrostatic energy remains largely negative and it is partially compensated by the negative Lennard-Jones energy contribution in acidic solution. This is no longer true in the neutral medium as shown in the total energy contributions (Figure ), as the energy balance is clearly unfavorable. This leads to the small energy barrier mainly observed in the PMF curves obtained in the neutral medium.

We propose here to explain these unfavorable electrostatic terms observed at 8 Å at both pH values and discussed above. To do so, we assume that the guest changes its orientation to form the energetically favored inclusion structure for complexes. In acidic solutions, the approach of the guest to the cavity allows strong electrostatic interactions of protonated 4AA with the sulfonate groups, with the amino group of the guest being oriented toward the sulfonate groups of the C4S. This orientation may also promote hydrogen bond formation between the amino group of 4AA and oxygen atoms of sulfonate groups. Upon complexation, the penetration into the host cavity of the phenyl part of the guest rather than the amino group will enhance its hydrophobic interactions with C4S. As a consequence, the orientation of the guest must change. To validate this hypothesis, we will deeply investigate the orientation coefficient α and the hydrogen bonds in the following part.

Further insight into the orientation of the guest as a function of the separation distance may be obtained from the calculation of α value. This parameter depends on the orientation of the guest within the calixarene cavity as explained in Figure S5 (Supporting Information). Three zones are defined, and at larger separation distances of the two centers of mass (zone C), a significantly larger α value (α > 1) is calculated. This result confirms a privileged position of the amino group toward the sulfonate groups, especially true in acidic solution (with the protonated 4AA). In zone B, the α value fluctuates in agreement with the change in the guest orientation upon complexation. Finally, the insertion of the guest inside the cavity leads to a smaller value of α (Zone A) that is thoroughly consistent with the penetration of the phenyl part of 4AA rather than the amino group.

Finally, we also report in Figure the total number of hydrogen bonds as a function of the separation distance for both pH values. We observe that the total number of hydrogen bonds is fewer at pH = 7 than that at pH = 1. From 10 Å, the amino group of 4AA forms hydrogen bonds with the sulfonate group of C4S. Upon complexation, these H-bonds are replaced by H-bonds with the oxygens of the water due to the change in the 4AA orientation in agreement with the fluctuation of the α value. Once the insertion has been made, the amino group is able to recover hydrogen bonds with the sulfonate groups as shown in Figure . This leads to a favorable negative electrostatic energy obtained from 5 Å as shown in Figure . As a conclusion, the formation of H-bonds and the evolution of α validate our hypothesis previously made.

Figure 8

H-bonds formed as a function of the distance between centers of mass for (C4S, 4AA) association at (a) acidic pH and (b) neutral pH. The colored curves correspond to a given type of H-bond, whereas the black curve corresponds to the sum of all H-bonds. The dotted lines represent the location of the deep minimum in the PMF profile for (C4S, 4AA).

Unambiguously, the small energy barriers obtained at 7.75 and 7.5 Å for C4S and C6S, respectively, at pH = 7 are related to the fact that the orientation of the guest is changed upon complexation. As already mentioned in the Introduction part, our previous results[19] tend to show that the 4AA is not significantly complexed by both C4S and C6S in neutral medium. The calculation of the free energy profiles obtained here and the presence of these barriers give a quantitative answer regarding this complexation and confirm that there is no association between the calixarene and the 4AA in water, in perfect agreement with our experimental and simulation results.[19] The same conclusions can be drawn for the (C6S, 4AA) association from the analysis of Figures S2–S4 (see the Supporting Information).

Thermodynamics Properties

The thermodynamic properties of the association calculated from the PMFs are reported in Table . As underlined in the methodological section, the calculation of the thermodynamic properties requires the knowledge of the d parameter. The dependence of the calculation on the value of d is shown for each complex in Figures S7–S9 (see the Supporting Information). In the case of the complex studied here, the use of eqs –9 represents an approximation because of the nonsphericity of the conformations of the guest and the host. However, we can estimate the thermodynamic properties that are listed in Table . For each calculation, we report the value of d, which remains reasonable for an association in terms of separation distances.

Table 1

Thermodynamic Properties of (Host, 4AA) Complexes Obtained from Free Energy Profilesa

host	ΔrG° (kJ mol–1)	ΔrH° (kJ mol–1)	TΔrS° (kJ mol–1)
β-Cdb (pH = 7)	–14	–21	–7
β-Cd (exp)	–19	–8	11
C4Sc (pH = 1)	–14	–25	–11
C4S (exp)	–11	–25	–13
C6Sd (pH = 1)	–10	–20	–10
C6S (exp)	–14	–18	–4

Value of d is given as a function of rmin defined as the position of the Gibbs free energy minimum.

d = rmin.

d = rmin + 2.5 Å.

d = rmin + 1.5 Å.

The calculation of the thermodynamic properties of the association from the PMFs shows a weak association in line with experimental interpretations. The complexes are characterized by the insertion of the guest into the host in line with negative values of ΔrH°. The values of the Gibbs free energy of association are on the same order of magnitude as those obtained from spectroscopy. In the case of the calixarenes, the deviation between simulations and experiments does not exceed 10% on the values of the enthalpy. The largest deviation is obtained with the β-Cd and can be explained by the fact that the reaction path is much more complicated for this association process. Indeed, the host molecule is then accessible from both sides and the guest is longitudinally inserted. These two features make the original assumptions of the spherical shapes of the guest and the host no longer valid. The deviation of the value of the enthalpy explains the change in the sign of the entropic term.

There are two predominant effects that influence the entropy effects: the loss of degrees of freedom and the desolvation of both guest and host molecules. From our simulation results, we saw that the entropic term is negative for all complexes. The loss of degrees of freedom is rather predominant. Nevertheless, for (BCD, 4AA) association, the value of the entropy is higher than the one for (CnS, 4AA) associations. This implies that the loss of degrees of freedom is higher for CnS than for BCD. These results are in perfect agreement with the ones we previously obtained.[19]

Conclusions

Molecular simulations have been performed to calculate the potential of mean force in the association between the 4AA and C4S, C6S, and β-CD in water. The calculation of PMF represents a valuable property to definitely conclude on the association between two molecules. The PMFs have been calculated using the method of the adaptive biasing force. Indeed, these PMFs can help in the interpretation of experiments when the experimental data do not allow a definitive conclusion of the association in neutral medium. The shape of the simulated PMF was rather complex, and additional simulations were performed to investigate the changes in the structure and energy along the reaction coordinate.

In a neutral medium, we have demonstrated that an energetic barrier of about 5 kJ mol–1 prevents any association, as suggested by experiments. This is the main conclusion of this work. We have also explained the reasons for the formation of the energetic barrier from the analysis of the van der Waals and electrostatic energy contributions. In acidic medium, the simulations confirm a weak association and an insertion complex in terms of thermodynamic properties. These results agree with experiments showing that the process is enthalpically favored. We conclude that only a combined approach of molecular simulations (structure and energetics) and experiments is able to elucidate the formation of these weak complexes.

Experimental Section

Methodology

Adaptive Biasing Force

The adaptive biasing force (ABF) method is an adaptative method that computes the potential of mean force w(ξ) along a reaction coordinate ξ. It was designed by Darve et al.[28,29] in 2001. Since then, many improvements have been proposed. For example, we could cite the extended-system adaptive biasing force (eABF) designed by Lesage et al.[30] or more recently the meta-eABF designed by Fu et al.[31] In ABF method, the total reaction coordinate ξ is divided in an appropriate number of bins. For each bin, the mean force, ⟨Fξ|ξ*⟩, is evaluated from a running averageBy applying an opposite force −⟨Fξ|ξ*⟩, the total average force acting on this system becomes close to zero after a brief equilibration. This allows avoiding significant energy barriers. The reaction coordinate is then well sampled. Now, let us look at the way we can compute the potential of mean force (PMF) with this method. By considering that n(Nstep) is the number of samples collected in bin k after Nstep, the running average force Fξ acting on ξ in the bin k is expressed aswhere corresponds to the lth force sample when ξ is in the bin k. In the LAMMPS collective variable module, the average force for the ABF method is calculated fromwith α(Nξ) being a scaling factor (going from 0 to 1) depending on the number of samples Nξ and ∇G̃(ξ) the estimate of the gradient of free energy for the point ξ in the reaction coordinate. Finally, the Gibbs free energy along the reaction coordinate from a state a to b, ΔGa→b, can be obtained by integrating the total biased force according toNevertheless, if the sampling is poor for a given bin k (typically at the beginning of the simulation), it is possible to have a nonequilibrium effect due to poor estimation of ⟨Fξ|ξ*⟩ that leads to an incorrect bias. To avoid these effects, we multiply Fξ by a function R(n(Nstep)) = min(1, n(Nstep)/n0). In practice, we took n0 = 200. The running average force becomes

Thermodynamic Properties

The thermodynamic properties of the association are calculated from the profile of the potential of mean force (PMF). The association constant[32]Ka can be calculated using the following expressionwhere NA corresponds to the Avogadro constant, kB to the Boltzmann constant, and w(h) to the PMF profile. rcyl corresponds to the mean radius of the cylinder where the guest can freely evolve for its in-plane movement when it is associated inside the cavity. This mean radius is calculated at each step h of the reaction coordinate. When the guest molecule is out of the cavity, we consider r as a constant equal to the highest value (i.e., when the guest is placed above the bigger rim). The parameter d, which is introduced in the different integrations, defines the upper limit of the association in terms of separation distances. This is not an easy parameter to define especially when the shapes of the guest and hosts are no longer spherical and quite different between different association processes. We took the route of representing different integrals as a function of the d parameter in the Supporting Information (Figures S7–S9).

The derived properties of the association such as the Gibbs free energy, enthalpy, and entropy can be calculated using the following expressions

Computational Details

As shown by experimental results,[19] the stoichiometry of these complexes is 1:1. As a consequence, the simulation boxes were composed of one host, e.g., a calixarene (see Figure c for C4S) or a β-Cd (Figure a); one guest, e.g., a 4AAB (Figure b); and 2000 water molecules. To maintain charge neutrality in the box, we may add, if needed, some Cl– or Na+. Our simulation boxes were cubic with a box length of 42 Å. The periodic boundary conditions were applied in the three directions. A scheme of each molecule is given along with the nomenclature in the Supporting Information (Figures S1–S3).

Concerning the force fields, we chose the TIP4P2005 model for water molecules[33] and the General AMBER Force Field (GAFF) for hosts and guests molecules.[34] However, for the guest molecule, the dihedral angle C–N=N–C was taken from the work of Heinz et al.[35] MD simulations were performed with the LAMMPS package.[36] The simulation parameters were chosen as follows:

The SHAKE algorithm[37] was used to constrain H-based bonds and the HOH angle for water.

Lennard-Jones crossing parameters were calculated using Lorentz–Berthelot rules (i.e., and , where i and j refer to the force centers and ϵ and σ are the energy parameter and diameter of atoms of types i and j, respectively).

Our simulations were performed in the NPT ensemble with a velocity-Verlet integrator and a timestep of 2 fs.

We used a Langevin[38,39] thermostat with a relaxation time of 2 ps and kept the temperature fixed at 300 K.

A Berendsen barostat[40] was used with a relaxation time of 2 ps and a bulk modulus of 1000 atm. The pressure was kept fixed at 1 atm.

We considered both long-range dispersion–repulsion tail correction[41] and long-range electrostatic interactions. Long-range Coulomb interactions were calculated using the PPPM style (particle-particle particle-mesh) with a relative error force of 10–4 which maps atom charge to a three-dimensional mesh.[42,43]

Our simulations were composed of an equilibration phase of 1 ns and an acquisition phase of 140 ns for a PMF going from 0 to 15 Å (280 ns for PMF going from −15 to 15 Å).