Cavity formation energy drives the accumulation of amphiphiles at the air–water interface

Abstract

Amphiphilicity is a key property that governs self-assembly, drug delivery and membrane permeability. We introduce a highly efficient approach based on polarisable continuum models to quantify this important property. Our modelling suggests that the accumulation of amphiphiles at the air–water interface is driven by the release of cavity formation energy in water.

---

Amphiphilic molecules possess both hydrophilic and hydrophobic groups, which enables them to mitigate unfavourable interactions at the interface between two chemical environments of distinct polarity. This property is highly significant in applications ranging from self-assembly (e.g. micelle formation),^1,2 to cleaning (e.g. detergents),³ enzyme modulators,⁴ and drug efficacy and delivery.⁵

In view of its importance, it is surprising that there is a lack of general and quantitative measures of amphiphilicity. Existing measures are limited to specific charged species⁶ and/or motifs such as α-helices.⁷ Yet, amphiphiles span a wide range of compounds including many small drug molecules and organic solvents, in addition to conventional surfactants and phospholipids.⁸ Notably, some sugars which are considered amongst the strongest hydrophiles have been demonstrated to display amphiphilic properties capable of solubilising hydrophobic drug-like molecules such as curcumin.⁹

Very recently, we introduced a general approach to quantify and predict amphiphilicity based on the depth of the free energy minimum at the interface between water and a model lipid membrane (labelled ΔG_int; see Fig. 1).¹⁰ This free energy well is a characteristic of amphiphiles observed in many molecular dynamics (MD) studies^11–16 though it has never been proposed as a measure of amphiphilicity until recently.

Fig. 1 Classical MD simulations indicate that amphiphiles display a characteristic free energy minimum at the water–membrane interface that can be used to quantify amphiphilicity.

A novelty of our approach was the construction of a fatty alcohol (dodecanol) bilayer model in lieu of a more complicated phospholipid membrane model, which would entail extensive configurational sampling and considerable computational resources to extract the desired free energy difference. Whilst fatty alcohols are not known to form bilayers in water, our simulations indicated that this metastable state is sufficiently long-lived (>20 ns), such that it is possible to extract permeation free energy profiles for a variety of small organic solutes.

These (umbrella sampling) free energy calculations are estimated to be ∼40-fold faster than a corresponding phospholipid bilayer simulation, and yielded amphiphilicities that correlated strongly with available experimentally measured surface tension data.¹⁰ Thus, the free energy well at the aqueous-fatty alcohol membrane interface represents a general approach to quantifying amphiphilicity of small drug candidates that will support efforts in drug and materials development.

Despite the substantial speed-up afforded by fatty alcohols over phospholipid models, the generation of the free energy profile in Fig. 1 would typically consume several thousand CPU hours. This communication introduces strategies that can significantly accelerate the calculation of amphiphilicity so that they may be used in routine studies and/or high throughput screening. At the same time, we exploit the simplicity of the new approaches to gain a deeper understanding of the forces that drive the affinity of amphiphiles to these interfaces.

What our fatty alcohol model membrane has demonstrated is a general approach for quantifying amphiphilicity from the free energy at the interface between two media of distinct dielectric constants. In this context, the use of a fatty alcohol bilayer is somewhat arbitrary and serves as a mimic of more complicated membrane models. Accordingly, the simplest interface would be the air–water interface, which could accelerate the simulations by several-fold.

However, before we skip over the aqueous-fatty alcohol membrane model, it is worth pointing out a very straightforward strategy to speed up the calculation of ΔG_int. Specifically, ΔG_int corresponds to the change in solvation free energy of the solute moving from bulk water to the interface, i.e. ΔG_int = ΔG_solv(interface) − ΔG_solv(water). This energy change can be straightforwardly obtained using, for example, a dual-topology free energy perturbation (FEP) simulation^17,18 that alchemically decouples the non-bonded interactions between the solute and bulk water and simultaneously switching on these interactions at the interface. Thus, there is no need to generate the entire membrane permeation free energy profile, though this could be important for applications such as membrane permeation predictions. Table S1 (see SI) shows that the ΔG_int obtained from umbrella sampling and FEP simulations are in excellent agreement (mean absolute deviation = 0.12 kcal mol⁻¹ for 12 solutes).

We applied the same FEP strategy to directly quantify the ΔG_int at the air–water interface and examine the correlation with values obtained from the fatty alcohol membrane model. Note that ΔG_int is also the free energy of adsorption at the air–water interface. Remarkably, for a test set of neutral organic amphiphiles of various functionalities (alcohols, amides, esters, ketones), we observed a very strong correlation (r² = 0.93; see Table S2). This validates our supposition that amphiphilicity can be quantified from the free energy at the interface between any two distinct dielectric media. These results are also consistent with our earlier work, which showed that ΔG_int correlated very strongly with experimental surface tension data;¹⁰ a measure of affinity of amphiphilic molecules to the air–water interface.^19,20

Alas, the above strategies still require relatively demanding MD simulations. In a typical MD simulation, one needs to perform thousands of sequential computations, each of which involves the calculation of thousands of not only intermolecular interactions between the solute and the solvent molecules, but also between solvent molecules, and in some cases also intramolecular interactions that contribute negligibly to solvation. Coarse-grained models^21,22 could alleviate the cost, but the MD simulations will still rely on force field parameters that may not be available or require validation.

In this regard, continuum solvent models represent an efficient alternative as the environment around the solute is modelled as a structureless dielectric medium.^23–25 This approach entails a small number of calculations on a single molecule of the solute typically carried out using modest DFT methods. The anisotropic solvent models,^26,27 heterogeneous dielectric generalised born (HDGB)²⁸ and COSMOmic²⁹ (COSMO-RS for micelles) are examples of more advanced continuum solvent models that can predict the free energy profiles across complex interfaces. For example, the COSMOmic approach has been demonstrated to predict the free energy well at the phospholipid DOPC-water interface for a set of 25 small drug-like molecules that are in good correlation (r² = 0.73) with data obtained from MD simulations.³⁰

Recently, Herbert and co-workers implemented and tested a heterogeneous PCM (HetPCM) model that is a relatively straightforward generalisation of standard homogeneous PCM models. In HetPCM, the effective dielectric at each point on the molecular boundary surface is allowed to vary depending on the external environment that each atom is deemed to be in contact with, and a switching function is used to ensure a smooth connection between the different dielectric regions. The authors showed that this is an efficient alternative to generalised Poisson equation solvers.³¹

In principle, details concerning the configurations of the molecule at the interface can be obtained from MD simulations, which we like to avoid for cost reasons. We employed the HetPCM model and considered a major simplification by assigning a dielectric constant (ε) of 78 (water) to all polar functional groups and electronegative atoms; all other atoms are assigned ε = 1 (dielectric constant of air/vacuum). The solvation energy in bulk water is obtained through a corresponding homogeneous PCM calculation. Non-electrostatic contributions to the solvation energy were computed using the SMD model.³²Fig. 2 displays the correlation between the ΔG_int at the air–water interface predicted by MD and HetPCM for a test set of 29 amphiphilic molecules (chemical structures provided in Fig. S1). This test set of neutral molecules includes alcohols, amines, esters, ketones and biomolecules (e.g. cholesterol) which carry a hydrophobic end and a hydrophilic end.

Fig. 2 Plots of amphiphilicities (kcal mol⁻¹) obtained with MD simulations versus those calculated by HetPCM electrostatics plus a non-electrostatic component (top), or versus those for just the non-electrostatic term (bottom).

Remarkably, the HetPCM calculations yielded a reasonably strong correlation (r² = 0.80) while requiring only <0.1% of the time needed for the MD simulation. The HetPCM solvation energies ΔG_int can further be broken down into the electrostatic (ΔG_es) and non-electrostatic cavity formation free energy (ΔG_cav) and dispersion–repulsion (ΔG_dr) contributions. Interestingly, the correlation appears to be dominated by the non-electrostatic component of the solvation energy change (Fig. 2b), where the r² is 0.77 (c.f. r² = 0.01 when only electrostatic contributions are considered; see Fig. S2 SI).

These observations can be rationalised as follows. First, for a typical amphiphilic molecule with a (small) polar head and a (large) non-polar body, the electrostatic component of the solvation free energy is dominated by the polar head group whether the molecule is positioned at the interface or bulk water. For instance, the molecule with the largest difference between the aqueous and interfacial electrostatic solvation energies, 5,7-diisopropyl-8-methylnapthanlen-2-ol, has only an OH head and a 38-atom body. Yet, even in this case, the interfacial value is 47% of the aqueous value. Thus, the similar electrostatic characteristics in bulk water and at the interface cannot account for the distinct differences in the amphiphilicities.

Let us now turn to the physical basis of the non-electrostatic contribution, which is typically composed of dispersion–repulsion and the cavity formation energy (i.e. energy required to create a cavity to accommodate the solute) in a PCM model.²³ The dispersion–repulsion (ΔG_dr) term can be expected to be attractive between the solute and the solvent under normal conditions. On the other hand, the cavitation term (ΔG_cav) represents the penalty of disrupting the strongly hydrogen-bonded water network when a partly/mostly non-polar molecule is forced from the interface into the bulk water.

The tendency of amphiphiles to accumulate at the air–water surface rather than in bulk water indicates that the repulsive cavitation term dominates the non-electrostatic contributions. The vaporization energy, molar mass, and density of water can be used to estimate the energy of the intermolecular interactions per unit volume, which is 583 kcal L⁻¹. To put this into perspective, the value for benzene is just 83 kcal L⁻¹. Accordingly, it takes many times the energy to displace the same volume of water in comparison with that for a typical organic solvent or vacuum.

Fig. 3 illustrates how cavity formation energy is reduced at the air–water interface using ethanol as an example. In this case, the overall favourable solvation free energy of ethanol in water (ΔG_solv = –5.01 kcal mol⁻¹)³² is presumably due to the electrostatic interactions with water, which is dominated by the polar –OH group (∼85% from our aqueous PCM versus interfacial HetPCM analysis). This outweighs the unfavourable cavity formation energy (Fig. 3; left). A significant part of that cavity formation energy is relieved when the amphiphile resides at the interface (Fig. 3; right), positioning the hydrophobic alkyl tail into the air. Whilst the role of cavitation energy at interfaces has been commented on in other anisotropic continuum solvent model studies,^33,34 the correlation/connection with affinity to the air–water interface has never been examined.

Fig. 3 Schematic representation of solvent cavity formation as a significant driver of the amphiphilic phenomenon.

The importance of the electrostatics provided by the head group is to anchor the molecule at the interface, as one can intuitively expect. As noted earlier, for our set of small surfactants (Fig. S1), the electrostatic interaction of the head group dominates that of the entire molecule. In this case, the average electrostatic binding energy of just the head group in water is –4.2 kcal mol⁻¹. Without this significant binding, e.g., CH₃CH₃ (ΔG_solv = 1.83 kcal mol⁻¹)³² instead of CH₃CH₂OH, we can envisage the molecule being fully expelled to the low-dielectric medium.

In addition to providing insight on the physical forces that drive the affinity of amphiphiles to the air–water interface, the results from Fig. 2b suggest a highly efficient way of predicting and quantifying amphiphilicity. The cavity formation free energy (ΔG_cav), in its current formulation within the PCM framework, does not rely on an underlying electronic structure calculation, but only requires the molecular geometry that goes into the classical equations. Notably, there are analytic forms of the free energy of cavity formation (e.g. proposed by Pierotti based on scaled particle theory)³⁵ used in PCM that require only the density of the solvent, radii of the solute and solvent, and solvent accessible surface area (SASA).

In summary, quick and reliable estimation of amphiphilicity can aid efforts to design more effective drugs and materials. This communication demonstrates how this important property may be quantified from the free energy at the interface of two distinct dielectric environments. Using the simplest air–water interface, this work showed that these free energies can be obtained from atomistic classical MD simulations, as well as from polarisable continuum solvent model calculations. Importantly, the data suggest that much of the driving force for the affinity of typical amphiphilic molecules to the air–water interface is driven by the reduction in cavity formation free energy (ΔG_cav).

The calculation of ΔG_cav is basically instantaneous compared with more advanced continuum solvent models such as COSMO-RS/FLATSURF³⁶ that entail orientation averaging and many DFT computations. The multiple-orders-of-magnitude saving in cost further opens the possibility for high-throughput computations in prospective applications in areas where amphiphilicity matters. This demonstrative preliminary study uses, for the most part, relatively straightforward systems with distinct hydrophilic and hydrophobic ends. Work is underway to extend this approach to more complicated amphiphiles, such as those that contain ionisable groups and/or multiple polar groups that can adopt different orientations at the interface.

JH acknowledges funding from the Australian Research Council (DP210102698 and FT240100584). Computational resources were provided to JH by Australian National Computational Infrastructure via National Computational Merit and UNSW Resource Allocation Schemes (cw7; DOI: 10.26190/PMN5-7J50), and to BC by RIKEN Information Systems Division (RB230026) and Institute for Molecular Science Japan (25-IMS-C070).

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: additional results noted in the main text, computational methods and details, analysis and DFT optimised coordinates. See DOI: https://doi.org/10.1039/d5cc04959h.

Notes and references

1. E. W. Kaler, A. K. Murthy, B. E. Rodriguez and J. A. N. Zasadzinski, Science, 1989, 245, 1371–1374.
2. T. Maki, R. Yoshisaki, S. Akama and M. Yamanaka, Polym. J., 2020, 52, 931–938.
3. D. Lombardo, M. A. Kiselev, S. Magazù and P. Calandra, Adv. Condens. Matter Phys., 2015, 2015, 151683.
4. G. Savelli, N. Spreti and P. Di Profio, Curr. Opin. Colloid Interface Sci., 2000, 5, 111–117.
5. J. Xi and H. Liu, Adv. Ther., 2020, 3, 1900107.
6. H. Fischer, M. Kansy and D. Bur, Chimia, 2000, 54, 640.
7. D. Eisenberg, R. M. Weiss and T. C. Terwilliger, Nature, 1982, 299, 371–374.
8. Water in Biological and Chemical Processes: From Structure and Dynamics to Function, ed B. Bagchi, Cambridge University Press, Cambridge, 2013, pp. 243–260.
9. S. Wong, J. Zhao, C. Cao, C. K. Wong, R. P. Kuchel, S. De Luca, J. M. Hook, C. J. Garvey, S. Smith, J. Ho and M. H. Stenzel, Nat. Commun., 2019, 10, 582.
10. N. A. Ahmad and J. Ho, J. Chem. Inf. Model., 2025, 65, 417–426.
11. W. Shinoda, Biochim. Biophys. Acta, Biomembr., 2016, 1858, 2254–2265.
12. R. W. Tejwani, M. E. Davis, B. D. Anderson and T. R. Stouch, Mol. Pharm., 2011, 8, 2204–2215.
13. A. L. Lomize and I. D. Pogozheva, J. Chem. Inf. Model., 2019, 59, 3198–3213.
14. C. H. Tse, J. Comer, S. K. Sang Chu, Y. Wang and C. Chipot, J. Chem. Theory Comput., 2019, 15, 2913–2924.
15. M. Sugita, S. Sugiyama, T. Fujie, Y. Yoshikawa, K. Yanagisawa, M. Ohue and Y. Akiyama, J. Chem. Inf. Model., 2021, 61, 3681–3695.
16. O. Engin, A. Villa, M. Sayar and B. Hess, J. Phys. Chem. B, 2010, 114, 11093–11101.
17. P. H. Axelsen and D. Li, J. Comput. Chem., 1998, 19, 1278–1283.
18. D. A. Pearlman, J. Phys. Chem., 1994, 98, 1487–1493.
19. C. M. Phan, ACS Omega, 2023, 8, 47928–47937.
20. M. Peng, T. T. Duignan, C. V. Nguyen and A. V. Nguyen, Langmuir, 2021, 37, 2237–2255.
21. J. Jin, Y. Han, A. J. Pak and G. A. Voth, J. Chem. Phys., 2021, 154, 044104.
22. S. K. Singh, A. Noroozi and A. Soldera, J. Chem. Phys., 2025, 162, 144501.
23. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3093.
24. C. J. Cramer and D. G. Truhlar, Chem. Rev., 1999, 99, 2161–2200.
25. J. M. Herbert, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2021, 11, e1519.
26. L. Frediani, R. Cammi, S. Corni and J. Tomasi, J. Chem. Phys., 2004, 120, 3893–3907.
27. M. F. Iozzi, M. Cossi, R. Improta, N. Rega and V. Barone, J. Chem. Phys., 2006, 124, 184103.
28. S. Tanizaki and M. Feig, J. Chem. Phys., 2005, 122, 124706.
29. A. Klamt, U. Huniar, S. Spycher and J. Keldenich, J. Phys. Chem. B, 2008, 112, 12148–12157.
30. M. Paloncýová, R. DeVane, B. Murch, K. Berka and M. Otyepka, J. Phys. Chem. B, 2014, 118, 1030–1039.
31. P. E. Bowling, M. Gray, S. K. Paul and J. M. Herbert, J. Chem. Theory Comput., 2025, 21, 1722–1738.
32. A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2009, 113, 6378–6396.
33. J.-B. Wang, J.-Y. Ma and X.-Y. Li, Phys. Chem. Chem. Phys., 2010, 12, 207–214.
34. K. Mozgawa and L. Frediani, J. Phys. Chem. C, 2016, 120, 17501–17513.
35. R. A. Pierotti, Chem. Rev., 1976, 76, 717–726.
36. A. Klamt, Fluid Phase Equilib., 2016, 407, 152–158.

---