A GNN-Based QSPR Model for Surfactant Properties

Abstract

Surfactants are among the most versatile molecules in the chemical industry because they can self-assemble in bulk solutions and at interfaces. Predicting the properties of surfactant solutions, such as their critical micelle concentration (CMC), limiting surface tension (${\gamma }_{cmc}$), and maximal packing density (${\Gamma }_{max}$) at water–air interfaces, is essential to their rational design. However, the relationship between surfactant structure and these properties is complex and difficult to predict theoretically. Here, we develop a graph neural network (GNN)-based quantitative structure–property relationship (QSPR) model to predict the CMC, ${\gamma }_{cmc}$, and ${\Gamma }_{max}$. Ninety-two surfactant data points, encompassing all types of surfactants—anionic, cationic, zwitterionic, and nonionic—are fed into the model, covering a temperature range of [20–30 °C], which contributes to its generalization across all surfactant types. We show that our models have high accuracy (R² = 0.87 on average in tests) in predicting the three parameters across all types of surfactants. The effectiveness of the QSPR model in capturing the variation of CMC, ${\gamma }_{cmc}$, and ${\Gamma }_{max}$ with molecular design parameters are carefully assessed. The curated dataset, developed model, and critical assessment of the developed model will contribute to the development of improved surfactants QSPR models and facilitate their rational design for diverse applications.

1. Introduction

Surfactants are amphiphilic molecules featuring hydrophilic (water-attracting) and hydrophobic (water-repelling) motifs simultaneously. Their amphiphilic nature facilitates the interaction between immiscible phases such as oil and water, thereby reducing their interfacial tension [1]. The latter underpins their extensive usage in numerous consumer and industrial applications. For example, surfactants are essential in detergent formulation, where they play a critical role in allowing aqueous solution to penetrate deeply into soiled fabrics and emulsify, eventually removing hydrophobic contaminants [2]. In the petroleum industry, surfactants facilitate the release of oil droplets trapped on oil-wet surfaces, thereby enhancing oil recovery. Similarly, surfactant usage is integral in pharmaceutical, food, and agriculture industries [3,4,5,6].

Most of the functionalities of surfactants originate from their ability to self-assemble in bulk solutions and at fluid–fluid interfaces. The self-assembly of surfactants and their impact on fluid–fluid interfaces can be characterized by many properties; the most important include ${\Gamma }_{max}$, CMC, and ${\gamma }_{cmc}$. The maximum surface excess concentration, ${\Gamma }_{max}$, measures the maximum density of surfactant molecules adsorbed at an interface. The critical micelle concentration, CMC, denotes the concentration at which surfactant molecules in solution begin to self-assemble and transition from monomeric to aggregated forms. The surface tension at CMC, ${\gamma }_{cmc}$, represents the surface tension of a solution when micelles first begin to form, marking the point at which the addition of more surfactant no longer notably lowers the surface tension. These three parameters are crucial metrics in the practical application of surfactants. For example, a lower CMC indicates a more efficient surfactant, as it requires a lower surfactant concentration to achieve micellization, which is critical in applications ranging from detergency [7] to drug delivery systems [8]. Likewise, a lower ${\gamma }_{cmc}$ corresponds to a more potent surfactant.

Designing surfactants with tailored ${\Gamma }_{max}$, CMC, and ${\gamma }_{cmc}$ is a central problem in their applications. Such design, when practiced by trial-and-error, can be time-consuming and costly due to the vast chemical space of surfactant molecules. Quantitative structure–property relationship (QSPR) models that predict these properties from surfactant molecules’ chemical structure can significantly accelerate their design. QSPR models for ${\Gamma }_{max}$, CMC, and ${\gamma }_{cmc}$ of surfactants can be developed using theoretical, simulation, and data-driven approaches. Models built by theoretical approach often fall short in accuracy, as they over-simplify the complex surfactant–surfactant and surfactant–fluid interactions at play [9]. For example, among polyethylene oxide surfactants (PEO), ${\Gamma }_{max}$ is a function of the number of ethylene oxide units only. The pure simulation approach powered by molecular dynamics (MD) modeling, while potentially more accurate, is computationally expensive and time-consuming. This is because ${\Gamma }_{max}$, CMC, and ${\gamma }_{cmc}$, all related to surfactant self-assembly, are governed by intricate collective, emergent behaviors of surfactant molecules that can only be accurately captured through long simulations of large MD systems. The experimental approach also has drawbacks to accurately measuring surfactant properties. The most common way of calculating ${\Gamma }_{max}$ is to utilize the Gibbs adsorption equation $\mathsf{\Gamma }~$ d$\gamma$/d ln(c), but this method requires a very large number of measurements and long equilibration time to be accurate, especially in low concentration ranges [10]. Similarly, inconsistent CMC values because of different determination techniques has been reported [11]. Given these limitations, there is a pressing need for a data-driven approach that can leverage datasets curated from experimental and simulation studies of surfactants [12,13]. QSPR models developed using this approach could potentially offer a cost-effective and reliable approach to designing surfactants for specific applications.

The data-driven approach for QSPR model development, however, is not without its challenges. One significant hurdle is a sparsity of available data [12]. In molecular design, the combinatorial space of possible surfactant structures is virtually infinite. One can narrow down this space by leveraging chemical knowledge. For example, a surfactant molecule can be divided into head- and tail groups. The former can be nonionic, cationic, anionic, or zwitterionic, while the latter can feature hydrocarbon, fluorocarbon, and silicone tails of different lengths. However, even with chemical knowledge, the potential molecular designs remain extraordinarily high-dimensional (e.g., numerous nonionic surfactant head groups are possible). Building a dataset for ${\Gamma }_{max}$, CMC, and ${\gamma }_{cmc}$ that comprehensively covers the high-dimensional chemical space of surfactants is difficult, if not impossible. The sparsity of data, coupled with the complexity and variability of surfactant behavior, complicates the development of robust data-driven models. As a result, there is an ongoing need for innovative approaches that can efficiently explore and model this high-dimensional space to achieve accurate predictions of ${\Gamma }_{max}$, CMC, and ${\gamma }_{cmc}$.

Many methods have been used to develop data-driven QSPR models with data scarcity. One classical method involves building a property function using molecular descriptors through molecular vector embedding. The molecular descriptors are selected either by domain knowledge [14,15,16,17,18,19,20], multi-linear regression [21,22,23], Pearson coefficient [24], support vector machine, [25,26,27] or random forests [28,29,30]. While they can show good performance in predicting surfactant properties, the descriptors from molecular descriptor libraries are not readily interpretable. For example, Seddon et al. reported good performance in predicting ${\Gamma }_{max}$ and CMC, but the descriptors revealed as important input features, e.g., VE3Sign_RG (the logarithmic coefficient sum of the last eigenvector from the reciprocal squared geometrical matrix) are hard to interpret intuitively. One can mitigate the data sparsity issue by narrowing the data scope to certain types of surfactants. For example, work reported by Kania et al. [28] investigated only three types of functional groups in nonionic surfactants. Similarly, work by Wu et al. [14] is limited to only nonionic surfactants, and work by Creton et al. [31] focuses only on perfluoroalkyl substance. While this method can construct good models for certain surfactant types, its relatively narrow chemical space limits its utility in exploring and searching for optimal surfactants for target applications.

With the rise of machine learning (in particular, deep learning), the graph neural network (GNN) has become a popular data-driven method for developing QSPR models. In GNN, a molecule with atoms and bonds can be treated as a graph with nodes and edges. Since GNN takes the constitutional atomic features directly from the molecule structure information, it enables neural network models to learn the complicated relationship without additional embedding vectors. Qin et al. [32] studied the CMC prediction of surfactants covering all main categories using the graph convolutional network (GCN) algorithm. They also built saliency maps to study the positive or negative contribution of atoms on CMC. Moriarty et al. [33] also used the GCN architecture and Gaussian process to predict CMC with uncertainty quantification. Moreover, Brozos et al. studied the predictive models of ${\Gamma }_{max}$ and CMC [34] and investigated the temperature dependence of CMC using a similar architecture [35].

The above recent studies suggest that GNN and its derivatives can be a powerful method for developing QSPR models for surfactants. Nevertheless, these models have not yet been utilized to predict the full set of ${\Gamma }_{max}$, CMC, and ${\gamma }_{cmc}$ when the dataset has a modest size, contains inevitable noise, yet covers a chemical space spanned by a wide spectrum of surfactants. Therefore, the question of whether an effective QSPR model for these properties can be developed under such practical constraints remains open. Furthermore, chemical knowledge of surfactant properties is not directly encoded in data-driven QSPR models (e.g., how ${\gamma }_{cmc}$ varies with the tail length), and overfitting remains a major concern even with advanced ML algorithms. Therefore, the question of whether these models can capture variations in (${\Gamma }_{max}$, CMC, and ${\gamma }_{cmc}$) while surfactant structure features remain open. Indeed, chemists have acquired such knowledge through decades of experimental, theoretical, and computational studies. It is worth examining to what extent the predictions of a data-driven QSPR model align with such knowledge.

In this work, we develop GNN-based QSPR models to predict ${\Gamma }_{max}$, CMC, and ${\gamma }_{cmc}$ using a dataset of 92 surfactants with a wide variety of head- and tail groups. The rest of the manuscript is organized as follows. Section 2 discusses the curation and noise/uncertainty of the dataset, the architecture of the QSPR model, and the selection of the model’s hyperparameters. Section 3 analyzes the performance of the QSPR model and assesses its capability to capture the variation of ${\Gamma }_{max}$, CMC, and ${\gamma }_{cmc}$ as a function of surfactant design parameters, including head- and tail types, and tail length.

2. Materials and Methods

2.1. Dataset Curation

To build a dataset of ${\Gamma }_{max}$, ${\gamma }_{cmc}$, and CMC for surfactant solutions, we collected the surface tension isotherms of 92 surfactants measured experimentally at water-air interfaces (see Table S1 in the Supplementary Materials for a list of the surfactants and their origins). There are 78 surfactants with hydrocarbon tails, 6 with fluorocarbon tails, and 8 with silicone tails. The number of anionic, cationic, zwitterionic, nonionic surfactants is 23, 5, 8, and 56, respectively. The sodium ion is the counterion for anionic surfactants, while the bromide ion is the counterion for cationic surfactants. Therefore, our dataset spans a vast chemical space despite its small size. A sample surface tension isotherm is shown in Figure 1a for the aqueous solution of sodium dodecyl sulfate (SDS), a canonical anionic surfactant. To extract ${\Gamma }_{max}$, ${\gamma }_{cmc}$, and CMC of the surfactant solution, we fit the surface tension isotherms to the Szyszkowski equation at surfactant concentration below CMC and ${\gamma }_{cmc}$ at higher concentrations:

$$\gamma =\left\{\begin{array}{c}{\gamma }_0-RT{\Gamma }_{max}\ln \left(1+K_{L}c\right) if c<CMC, nonionic surfactant\\ {\gamma }_0-2RT{\Gamma }_{max}\ln \left(1+K_{L}c\right) if c<CMC, ionic surfactant\\ r_{cmc} if c>CMC\end{array}\right.$$

(1)

where ${\gamma }_0$ is the surface tension in the absence of surfactants (72 mN/m for the air-water interfaces at 25 °C), R is the ideal gas constant, T is the absolute temperature, $K_L$ is the Langmuir constant, and c is the surfactant concentration. We only gather isotherms at temperatures between 20 and 30 °C, and no additional salts are presented in the surfactant solution. The maximum surfactant density at interfaces ${\Gamma }_{max}$ reflects the slope of the surface tension curve before reaching the CMC, and its value is calculated based on the nonlinear least square optimization function scipy.optimize.curve fit in the scipy library [36], implemented in Python (version 3.12.1). For ${\gamma }_{cmc}$ and CMC values, we take either reported values from the literature or the intersection of the two guidelines before and after CMC (see Figure 1a) suggested by authors of the surface tension isotherm. Values of ${\Gamma }_{max}$, ${\gamma }_{cmc}$, and CMC thus obtained are listed in Table S1 of the Supplementary Materials. Figure 1b–d show the histograms of ${\Gamma }_{max}$, ${\gamma }_{cmc}$, and logCMC of the surfactants curated. The histogram of ${\Gamma }_{max}$ is skewed toward left while that of ${\gamma }_{cmc}$ is skewed toward right. The logCMC histogram is balanced. The ranges of ${\Gamma }_{max}$, ${\gamma }_{cmc}$, and logCMC are [1.6, 11] × 10⁻⁶ mol/m², [14.5, 45.5] mN/m, and [−4.8, −0.9] (mol/L) in logarithmic scale, respectively. Our dataset thus covers a rather broad range of ${\Gamma }_{max}$, ${\gamma }_{cmc}$, and CMC, which is consistent with the wide variety of surfactants it includes.

For many surfactants in our dataset, complete information for their properties is not available: often the surface tension isotherm was provided only at concentrations considerably below CMC and CMC was not reported. For example, the surface tension isotherm of sodium octyl carboxylate (C₈OONa) reported by Neys et al. [38] shows a monotonic decrease in surface tension, showing that micelles do not form at the surfactant concentration studied. On the other hand, the data of sodium octyl sulfate from Feinerman et al. [39] show its full profile that contains CMC information. In total, only 64 of the surfactants have known CMC values and the surface tension at CMC ${\gamma }_{cmc}$. The rest do not have reported CMC values while ${\Gamma }_{max}$ can still be determined with the reported surface tension isotherm.

It is also found that the ${\Gamma }_{max}$ extracted from surface tension isotherms is sometimes sensitive to the surface tension data near CMC. As illustrated in Figure 2, for sodium octyl sulfate (C₈SO₄Na), ${\Gamma }_{max}$ is fitted to be 5.02×10⁻⁶ mol/m² when the seven data points [(logC [M], γ [mN/m]) = (0.002, 71.388), (0.006, 70.214), (0.020, 63.732), (0.038, 55.975), (0.055, 50.566), (0.089, 43.931)] and CMC data (0.114, 40.563) are used in fitting. However, the fitted ${\Gamma }_{max}$ decreases to 4.01 ×10⁻⁶ mol/m² when an additional point near CMC [(logC [M], γ [mN/m]) = (0.128, 40.614)] (colored as a green dot in Figure 2) is used in fitting. Such sensitivity, along with the fact that CMC is not always known with great accuracy [11], implies that some noise and uncertainty inevitably exist in our dataset.

2.2. Architecture of Machine Learning Model

Our ML model seeks to predict ${\Gamma }_{max}$, ${\gamma }_{cmc}$, and CMC of a surfactant using its chemical structure (specifically, its SMILES string) as an input. Figure 3 shows the overview of the architecture of our model. The model consists of three modules: (1) an initial feature encoding module that converts the SMILES string of a surfactant into an adjacency matrix, an atom feature matrix, and a bond feature matrix, (2) a molecular embedding module that learns the molecular embedding vector via message passing networks, and (3) a fully connected neural network that has three end nodes corresponding to ${\Gamma }_{max}$, ${\gamma }_{cmc}$, and logCMC of the surfactant. Below we outline only the first two modules because the third module is straightforward.

In the initial feature encoding module, a SMILES string representing the surfactant structure is used to generate an adjacency matrix, an atom feature matrix, and a bond feature matrix. The adjacency matrix is a square matrix encoding the connectivity between atoms (1 as connected or 0 not connected), and its size is the product of the number of atoms in the surfactant molecule (N). The atom feature matrix stores atomic features. We selected the features listed in

Table 1. Atom features considered for surfactant molecules ^1.

Atom Features	Descriptions	Size
Atom type	Type of atom (e.g., C, N, O), by atomic number	100
Number of bonds	Number of bonds the atom is involved in	6
Formal charge	Integer electric charge assigned to atom	5
Chirality	Unspecified, tetrahedral CW/CCw, or others	4
Hybridization	sp, sp2, sp3, sp3d, or sp3d2	5
Number of hydrogens	Number of bonded hydrogen atoms	5
Aromaticity	Whether this atom is part of an aromatic system	1
Atomic mass	Mass of the atom, divided by 100	1
Partial charge 2	Non-integer electric charge assigned by the OPLS force field	1
$\mathrm{LJ} \mathrm{potential} \sigma$ 2	The distance parameter of the Lennard-Jones potential assigned by the OPLS force field	1
$\mathrm{LJ} \mathrm{potential} ϵ$ 2	The potential well depth parameter of the Lennard-Jones potential assigned by the OPLS force field	1

¹ All features are one-hot encoded except the atomic mass and MD-computed features. ² MD-computable features.

, where the i-th row and j-th column of the matrix corresponds to the j-th feature of i-th atom. The bond feature matrix contains the features of chemical bonds between atoms, and

Table 2. Bond features considered for surfactant molecules ^1.

Bond Features	Descriptions	Size
Bond type	Single, double triple, or aromatic	4
Conjugated	Whether the bond is conjugated	1
In-ring	Whether the bond is part of a ring	1
Stereo	None, any, E/Z, or cis/trans	6

¹ All features are one-hot encoded.

lists the features we selected. In the molecular embedding module, we chose the directed message-passing neural networks (D-MPNN) algorithm developed by McGill and colleagues [40,41] as the message-passing network. This is because D-MPNN has excellent capability in molecular property prediction. For example, its effectiveness has been demonstrated with [21] public datasets including QM9, ESOL, and FreeSolv, as well as datasets collected by other authors. [42,43,44] In D-MPNN, the molecular embedding module starts by introducing the directed bond feature $e_{vw}^d$. This vector is devised to avoid messages passed along any path of the form v₁, v₂, …, v_n (v_i = v_i+2 for any i) [39] and is given by:

$$e_{vw}^d=cat(x_v,e_{vw})$$

(2)

where $x_v$ and $e_{vw}$ are the atom feature vector of atom v and the bond feature vector of the bond between atoms v and w. These vectors are a subset of the atom feature matrix and bond feature matrix obtained from the previous feature encoding module, respectively. It is noted that $x_w$ is excluded from the concatenation because it is directional. The directed edge feature vectors are used to initialize the edge hidden feature $h_{vw}^0$ with a learned matrix $W_i$ and activation function $\tau$, followed by the message-passing process between neighbor edge hidden features.

$$h_{vw}^0=\tau \left(W_i{e}_{vw}^d\right)$$

(3)

$$h_{vw}^{t+1}=\tau (h_{vw}^0+W_j\sum _{k\in (N\left(v\right)\backslash w)}{h}_{kw}^t)$$

(4)

Note that during one cycle of this message-passing layer to update $h_{vw}^t$, the information from the bond itself, the neighbor atoms, and the other bonds of one of the neighbor atoms are shared. This update is iterated across all bonds as one cycle of the message-passing layer. Once this process iterates T times of message-passing layers to obtain final hidden bonds feature vectors $h_{ij}^T$ where T is the number of messages passing layers and a hyperparameter, they are used to build the hidden atom features $h_v$ via

$$q=cat(x_v,\sum _{k\in \left(N\left(v\right)\right)}{h}_{kw}^T)$$

(5)

$$h_v=\tau \left(W_{k}q\right)$$

(6)

The hidden atom features $h_v$ are aggregated into the molecular embedding vector $h_m$ by taking the average of the hidden atom features. $x_w$ is an optional vector of additional molecular descriptors that can concentrate with the molecular embedding vector. Here we add two geometric descriptors that can be computed using MD simulations (Gromacs 2020.3, access date: 1 November 2023–13 November 2024), and they are described in the next section.

$$h_m^{\prime }={\displaystyle \frac{1}{n}}\sum _{v\in V}{h}_v$$

(7)

$$h_m=cat(h_m^{\prime },x_m)$$

(8)

where n is the size of $h_v$. The molecular embedding vector is fed into the fully connected layers which have three nodes at the end targeting the surfactant properties ${\Gamma }_{max}$, ${\gamma }_{cmc}$, and logCMC.

2.3. Additional Features and Descriptors by MD Simulations

Directly predicting surfactant properties such as maximum surface density or CMC through MD simulations is challenging, as it requires long simulation times. However, MD simulations can still provide valuable insights. To explore their potential utility, we investigated the use of computationally inexpensive MD simulations as molecular features and descriptors to enhance the predictive power of our property models. Therefore, in addition to the default list of the features, we introduce additional atom features and descriptors that are computable using MD simulations (see Figure 4). The partial charge and the Lennard-Jones parameters $\sigma$ and $ϵ$ of each atom in a surfactant molecule are essential atom-level parameters that govern the strength of electrostatic and van der Waals interactions between atoms. Therefore, they are introduced and later examined whether they help improve model performance. The MD-computable features in Table 1 are obtained by the OPLS force fields using the LigPargen web server (https://zarbi.chem.yale.edu/ligpargen/, accessed on 13 November 2024) [45] and are concatenated in the initial atom features $x_v$. In Figure 4a, the MD-computable features of a carbon, an oxygen, and a sulfur atom of SDS are displayed. It is noted that the hydrogen atoms connected to the carbon atoms are weakly positively charged, so the tail remains overall neutral as intended.

It is known the geometrical properties of surfactant molecules play a key role in their self-assembly and adsorption at water-air interfaces [1]. Although measuring these geometrical properties is difficult in experiments, they can be computed through MD simulations. Therefore, it is worth computing some of these geometrical properties and investigating whether incorporating them into the ML model can improve its performance. Here we build a 5 × 5 × 5 nm³ size MD system where a single surfactant molecule is immersed in water (see Figure 4b). The system is run for 15 ns in the NVT ensemble, with the last 10 ns taken as the production run. The subgroup atoms (e.g., headgroup or tailgroup) of the surfactant molecule are selected and fit to an ellipsoid by equating the moment of inertia of the selected atoms to that of an ellipsoid (see Supplementary Materials for the full detail of the MD setting and ellipsoid fitting). We then use the surface area of the fitted ellipsoid encapsulating headgroup and tailgroup as descriptors (see

Table 3. The list of additional molecular descriptors.

Additional Molecular Descriptors	Descriptions	Size
Surface area of tailgroups	Surface area of a fitted ellipsoid encapsulating tailgroups	1
Surface area of headgroups	Surface area of a fitted ellipsoid encapsulating headgroups	1

). These descriptors are concatenated in the molecular embedding vector as $x_m$.

2.4. Hyperparameter Setting

We use random data splitting based on 0.8/0.1/0.1 for the train/validation/test.

Table 4. Hyperparameters of the ML model, with the selected values highlighted in bold.

Hyperparameters	Values
Number of message-passing layers	2, 3, 4
$\mathrm{The} \mathrm{size} \mathrm{of} \mathrm{hidden} \mathrm{bond} \mathrm{feature} h_{vw}$	200, 300, 400
Number of fully connected layers	2, 3, 4
Dropout rate	0, 0.1, 0.2

shows the list of hyperparameters tried and their final selection obtained using the grid-search algorithm. The criteria to determine the best set of hyperparameters is to minimize the average of the three surfactant properties’ root mean square error (RMSE). Each RMSE value is averaged from five ensembles initialized with different random weights, which has been demonstrated to help increase model accuracy [41,46].

3. Results

3.1. Model Evaluation

We vary the choice of the initial atom features $x_v$ as well as the molecular descriptors $x_m$ to investigate how different sets of inputs affect the model performance.

Table 5. The list of atomic features and additional molecular descriptor sets (the default set refers to the first eight atom features in Table 1).

Set	Atom Feature Set	Additional Molecular Descriptors
1	Default	None
2	$\mathrm{Default}+\mathrm{partial} \mathrm{charge}+\sigma$$+ϵ$	None
3	Default	Surface areas of tailgroup and headgroup
4	$\mathrm{Default}+\mathrm{partial} \mathrm{charge}+\sigma$$+ϵ$	Surface areas of tailgroup and headgroup

lists the combination of the atom features and the additional molecular descriptors used as inputs. Specifically, we tried four sets of inputs by including or excluding the MD-computable atom features (partial charge, $\sigma$, and $ϵ$) and MD-computed molecular descriptors (the surface area of surfactant tailgroup and headgroup). The default set in Table 5 refers to the first eight atom features listed in Table 1, which do not require MD simulations to obtain.

Figure 5 shows the RMSE values of the target surfactant properties with different input sets. For ${\Gamma }_{max}$, Set 3 shows the smallest RMSE (1.04) in the test set. Set 3 also shows the smallest test RMSE (2.46) for ${\gamma }_{cmc}$. Set 1 shows the smallest test RMSE (0.28) for logCMC. It appears that the concatenation of information regarding the surface area of a surfactant’s headgroup and tailgroup into the molecular embedding vector helps the model better capture its properties. For example, the RMSE is reduced by 7% in ${\gamma }_{cmc}$ with the surface area information, and the RMSE is reduced by 3% when predicting ${\Gamma }_{max}$. However, its effect does not always improve the model’s performance. When predicting logCMC, the RMSE increases by 10%. Overall, introducing MD computed features into our model only marginally affects its performance. We envision that other MD-computed features or other ways of incorporating MD-computed features into the model may improve its performance. A systematic investigation of such practice is beyond the scope of this work. Instead, below we select Set 1 in Table 5 as the baseline input set to further analyze the model performance.

Figure 6 shows the parity plot of the three target properties of all 92 surfactants in our dataset. The values are the averaged values of the five ensembles predicted based on Set 1 input in Table 5. The R² value in the train set for ${\Gamma }_{max}$, ${\gamma }_{cmc}$, and logCMC is 0.94, 0.95, and 0.95, respectively, and that in the test set is 0.82, 0.90, and 0.90, respectively. Though all three surfactant properties are predicted reasonably well by the model in the test set as well as in the training set, it appears that the model is overfitting with respect to ${\Gamma }_{max}$. This is likely due to the inherent uncertainty/noise of ${\Gamma }_{max}$ as illustrated in Figure 2. A similar phenomenon has also been reported by Seddon et al. [28], where the prediction of ${\Gamma }_{max}$ is weaker than that of other properties logKL and logCMC. Nevertheless, the parity plots show that our model performs well in predicting the three surfactant properties considered here.

3.2. Surfactant Parameter Space Exploration

The ultimate goal of a QSPR model is to predict surfactant properties over a wide chemical space. The good performance of our model shown in the previous section for the surfactants in our dataset is encouraging. However, the dataset does not cover the high-dimensional chemical space of surfactants well (e.g., many hydrocarbon surfactants with nonionic headgroups are included, but far fewer fluorocarbon surfactants with similar headgroups exist in the dataset), and overfitting is possible even with advanced ML models. Therefore, it is important to examine the behavior of the model with a wide variety of chemical structures. Because rigorous validation is difficult due to a lack of experiment data, we shall compare the semi-quantitative trend and scaling law predicted by the model with those expected from chemical knowledge. Because our model is trained using very sparse data over the high-dimensional space of surfactant structure and the data likely contain some noise, capturing those trends and scaling law, even though possible for experienced chemists, is not guaranteed a priori for the data-driven model we develop here.

Below, we generate the landscape of surfactant properties as a function of several molecular design parameters based on our models. Specifically, we vary the surfactant structure to investigate the variation of ${\Gamma }_{max}$, ${\gamma }_{cmc}$, and logCMC. While numerous groups of surfactants vary in structure, we select the nonionic head (polyethylene oxide, PEO), and anionic head types (sulfate, sulfonate, and carboxylate) as the representative surfactants, because these functional groups are commonly found in surfactants. Figure 7 shows the schematic of the surfactant head- and tail types used in this study. For PEO, we vary the number of carbons in its tail $N_C$ and the number of ethylene oxide (EO) units $N_{EO}$ in its headgroup. With the PEO head, we limit the choice of tail types only to hydrocarbon. For anionic surfactants, we vary not only the surfactant head types (sulfonate, sulfate, and carboxylate) but also the tail types to hydrocarbon (C_xH_y) and fluorocarbon (C_xF_y) tails. The chain length $N_C$ of the hydrocarbon and fluorocarbon tails is also varied.

3.2.1. Polyethylene Oxide (PEO) Surfactants

Figure 8a depicts the variation of ${\Gamma }_{max}$ as a function of $N_C$ and $N_{EO}$ for 6 ≤ $N_C$ ≤ 20 and 4 ≤ $N_{EO}$ ≤ 16 predicted by our model, and the experimental data available in our dataset are marked by color-coded crosses and triangles (all symbols are color-coded based on $N_{EO}$). Despite the minimal experimental data in the parameter space examined here, the predicted dependence of ${\Gamma }_{max}$ on $N_C$ and $N_{EO}$ is smooth, suggesting that our model is likely not overfitted.

The predicted dependence of ${\Gamma }_{max}$ on $N_C$ and $N_{EO}$ is largely consistent with the trend observed for PEO surfactants and that expected for nonionic hydrocarbon surfactants. Specifically, ${\Gamma }_{max}$ is predicted to increase with decreasing $N_{EO}$, which agrees with that observed for PEO surfactants [1] and also aligns with the expectation that the packing density of nonionic surfactants at fluid interfaces typically increases as the size of their hydrophilic head decreases [1]. The model predicts that ${\Gamma }_{max}$ generally increases with the tail length $N_C$ although an oscillatory behavior is observed at $N_C$ ≳ 12 when $N_{EO}$ is less than about 10. The former is consistent with the observation that ${\Gamma }_{max}$ of PEO surfactants generally increases with the length of their hydrophobic group, which can be attributed to the stronger lateral interactions and thus denser packing of their tails [1]. The reliability of the more complicated dependence of ${\Gamma }_{max}$ at small $N_{EO}$ and large $N_C$ predicted by the model cannot be ascertained yet. However, at small $N_C$, the favorable PEO-water and PEO-PEO interactions and hydrophobic interactions between hydrocarbon tails may be insufficient to enforce a dense packing of long hydrocarbon tails due to the entropic penalty associated with the latter. Therefore, at small $N_{EO}$, the dependence of ${\Gamma }_{max}$ on $N_C$ likely would not follow the monotonic trends found at large $N_{EO}$, and a more complicated relation, as suggested by our model, may emerge.

At a more quantitative level, the model seems robust enough to overcome likely noise in the training dataset and can capture an important scaling law. At ($N_C$, $N_{EO}$) = (12, 6), the predicted ${\Gamma }_{max}$ is 3.16 µmol/m², which is 26% lower than the experimental data; a larger discrepancy of 75% is observed at ($N_C$, $N_{EO}$) = (16,6). In light of the uncertainties of experimentally derived ${\Gamma }_{max}$ highlighted in Section 2 and Figure 2, it is likely that the reported ${\Gamma }_{max}$ has substantial noise and the model has avoided overfitting these data. The minimum area per PEO surfactant molecule, given by Amin = 1/${\Gamma }_{max}$, follows the well-known scaling law $A_{min}$∼ $\sqrt{N_{EO}}$ [47]. Figure 9 shows the relationship between $A_{min}$ and $\sqrt{N_{EO}}$. While the model prediction satisfies the scaling law well, the experimental data do not align with the scaling law. The latter suggests that the experimental data in Figure 8a likely contain significant uncertainty. Despite such uncertainty, the model can capture the scaling law of ${\Gamma }_{max}$ established in the literature. This is likely because the experimental data within our dataset but outside of the ($N_C$, $N_{EO}$) range considered in Figure 8 work as regularizers to prevent the overfit to the noisy ${\Gamma }_{max}$. Moreover, the nature of multi-task prediction utilized in our model framework helps to generate robust high-level representations [48].

Figure 8b shows the predicted variation of ${\gamma }_{cmc}$ for 6 ≤ $N_{C }$ ≤ 20 and 4 ≤ $N_{EO}$ ≤ 16, and the experimental data available in our dataset are marked by color-coded crosses and triangles. The consistency between the model prediction and available experimental data is greatly improved compared to that for ${\Gamma }_{max}$. The largest difference in ${\gamma }_{cmc}$ is 11% and occurs at ($N_C$, $N_{EO}$) = (10, 4), where the model prediction and experimental data are 29.86 and 26.57 mN/m, respectively. In terms of qualitative trend, ${\gamma }_{cmc}$ generally decreases as $N_C$ increases and $N_{EO}$ decreases. These trends are difficult to establish using the experimental data in the ($N_C$, $N_{EO}$) space examined here alone. Nevertheless, the reduction of surface tension by surfactants generally increases with their tail length [49] and their packing density at water-air interfaces [1]. Figure 8a shows that, overall, ${\Gamma }_{max}$ decreases (increases) with an increase of $N_{EO}$ ($N_C$). Therefore, the trends shown in Figure 8b seem reasonable and warrant further experimental validations.

Figure 8c shows the predicted variation of CMC with $N_C$ and $N_{EO}$. The agreement between the predictions and available experimental data is comparable to that for ${\gamma }_{cmc}$. Overall, logCMC decreases as $N_C$ increases for a given $N_{EO}$, as has been observed for a wide variety of surfactants [50]. The predicted logCMC decreases as $N_{EO}$ decreases for a given $N_C$, and the reduction of logCMC with $N_{EO}$ is much weaker than that per methylene unit in the hydrophobic chain. Both trends are consistent with the literature [1], although the validity of the more subtle details of logCMC variation shown in Figure 8c, e.g., the insensitivity of logCMC relative to $N_{EO}$ at $N_C$ ≤ 10, are difficult to ascertain.

3.2.2. Anionic Surfactants

In this section, we investigate the variation of ${\Gamma }_{max}$, ${\gamma }_{cmc}$, and logCMC of anionic surfactants with different types of head- and tail groups, and the length of tail groups. We consider three anionic head groups (sulfate, sulfonate, and carboxylate) and two tail group types (hydrocarbon and fluorocarbon). The number of carbon atoms $N_C$ of hydrocarbon and fluorocarbon tails is varied from 6 to 20.

${\Gamma }_{max}$ of anionic surfactants. Comparison of the data in Figure 10a,d shows that ${\Gamma }_{max}$ is more affected by the tail type than the head type, with ${\Gamma }_{max}$ of fluorocarbon surfactants always higher than that of hydrocarbon surfactants for all three head types. This result aligns well with the fact that the cross-section area of the fluorocarbon chain (about 27 $A^2$) is much higher than that of the hydrocarbon chain (about 20 $A^2$) [51] For both fluorocarbon and hydrocarbon surfactants, ${\Gamma }_{max}$ increases with the tail length for $N_C$ ≤ 12; at larger $N_C$, the effects of tail length on ${\Gamma }_{max}$ becomes marginal. These trends for the hydrocarbon surfactants can be loosely inferred from the limited, and somewhat noisy experimental data in our dataset (see the crosses and triangles in Figure 10a), and agree well with the experimental observations that, for ionic surfactants, their ${\Gamma }_{max}$ increases with $N_C$ for short tails but varies little once the $N_C$ exceeds 10 [1].

The model predicts a complicated dependence of ${\Gamma }_{max}$ on the head type. For hydrocarbon surfactants with the same $N_C$, ${\Gamma }_{max}$ follows the order of carboxylate > sulfonate > sulfate at $N_C$ ≥ 10; at smaller $N_C$, the order changes to sulfonate > carboxylate ≈ sulfate. For fluorocarbon surfactants, ${\Gamma }_{max}$ consistently follows the order of carboxylate > sulfate > sulfonate. These trends do not always align with the trend that may be estimated from the limited experimental data in Figure 10a. Inferring these trends from data for other surfactants is not straightforward given that ${\Gamma }_{max}$ is governed by the complex interplay between tail–tail, head–head, and head–water interactions, which depends on the structure of head and tail in a manner that defies simple extrapolation in the high-dimensional chemical space. While this makes validating the model prediction difficult, it also highlights the potential of ML-based models in suggesting new experimental studies to explore complex structure–property relations.

${\gamma }_{cmc}$ of anionic surfactants. Comparison of the data in Figure 10b and e shows that for surfactants with the same anionic headgroup, the predicted ${\gamma }_{cmc}$ is lower for fluorocarbon tailgroups than for hydrocarbon surfactants. While experimental data supporting this trend is unavailable in the parameter space we examine here, this trend aligns well with the fundamental studies of surface tension reduction by surfactants. Specifically, as reviewed by Czajka, Hazell, and Eastoe [52], surface tension is governed by the additive contributions of dispersion and polar interactions. For surfactants with the same headgroups but different tail groups, their difference in ${\gamma }_{cmc}$ is dominated by the difference in the dispersion interactions. Because of the bulkier size of fluorocarbon tails and the lower polarizability of fluorine atoms than hydrogen atoms, the dispersion contribution of fluorocarbon tails is smaller than that of hydrocarbon tails [52]. Consequently, ${\gamma }_{cmc}$ of fluorocarbon anionic surfactants is lower than their hydrocarbon counterparts.

Our model predicts three main trends for the dependence of ${\gamma }_{cmc}$ on the tail length of sur factants considered. First, ${\gamma }_{cmc}$ of all hydrocarbon surfactants decreases with the tail length $N_C$, which does not contradict the two data points in our dataset. Second, for fluorocarbon surfactants, the decrease of ${\gamma }_{cmc}$ with $N_C$ is observed for sulfate and sulfonate headgroups, in contradiction to the data for sulfonate surfactants (cf. the two blue crosses in Figure 10f), which is likely due to experimental noise. Finally, for fluorocarbon surfactants with carboxylate headgroups, their ${\gamma }_{cmc}$ is independent of tail length $N_C$. The first two predictions are consistent with the general trend that ${\gamma }_{cmc}$ decreases with increasing tail length. The fact that this trend is recovered robustly in the presence of some noise in the training data is encouraging. The third prediction is likely valid. The limiting surface tension of fluorocarbon surfactants is usually in the range of 15–25 mN/m [52]. For fluorocarbon carboxylate surfactants, a ${\gamma }_{cmc}$ of about 18 mN/m is already achieved at $N_C$ = 6, making the reduction of ${\gamma }_{cmc}$ with increasing $N_C$ unlikely.

Our model predicts that the ${\gamma }_{cmc}$ of surfactants with the same tailgroup follows the trend of sulfonate > sulfate > carboxylate, i.e., carboxylate headgroups are most effective in lowering surface tension, followed by sulfate and then sulfonate headgroups. Such a trend is difficult to establish with the experimental data shown in Figure 10b,e. However, in a related study, Liu and colleagues [53]. studied the surface activities of sodium alkyl polyethylene oxide carboxylate (AE₃C), sodium alkyl polyethylene oxide sulfate (AE₃S), and sodium alkyl polyethylene oxide sulfonate (AE₃SO). The terminal motifs of these surfactants’ headgroups are the same as those studied here, and their tailgroups have 12–14 carbon atoms. They found that the ${\gamma }_{cmc}$ of surfactants with sulfonate group is highest, followed by those with sulfate group, and then by those with carboxylate. This study thus lends indirect support to our model prediction, even though a direct prediction of our prediction is not yet available.

logCMC of anionic surfactants. Comparison of the data in Figure 10c,f shows that logCMC of hydrocarbon surfactants decreases sharply and almost linearly with the tail length $N_C$. For fluorocarbon surfactants, the decrease of logCMC with $N_C$ is much slower. The former prediction is consistent with the general trend of logCMC of anionic hydrocarbon surfactants [1] and the data for C_nH2_n+1OSO₃Na (n = 8 to 14) reported by Srinivasan and Blankschtein [54]. The latter prediction, however, is inconsistent with some reported data in the literature. For example, the data presented by Kancharla and colleagues showed that the logCMC of C_n−1F_2n−1COONa decreases by 1.7 as n increases from 5 to 9 [55], which is even more steep compared to that for C_nH_2n+1OSO₃Na [54] The models by Creton [31] also support the steep slope of the graph, showing that the logCMC of C_n−1F_2n−1COONa decreases from −1 to −3 as n increases from 6 to 10, and logCMC of C_n−1F_2n−1SO₃Na decreases from −2 to −3 as n increases from 6 to 9. This inconsistent trend observed in our model compared to the references arises likely due to the limited number of fluorinated surfactant data.

For hydrocarbon surfactants with the same tail length, our model predicts that their logCMC follows the order of sulfate ≈ sulfonate > carboxylate. While there is no logCMC data in our dataset to support these trends, they seem unreliable based on the data reported by Sadeghi and Shahabi: they found that the logCMC of C₁₂H₂₅SO₃Na is higher by about 0.18 than that of C₁₂H₂₅SO₄Na [56]. For fluorocarbon surfactants with the same tail length, our model predicts that their logCMC follows the order of carboxylate > sulfate > sulfonate. This trend seems to be supported by a pair of data points in our dataset (cf. the red triangle and blue cross in Figure 10f), but its robustness is unclear at present.

4. Conclusions

In summary, we develop a GNN-based QSPR model to predict surfactant properties ${\Gamma }_{max}$, ${\gamma }_{cmc}$, and CMC.Ninty-two surface tension isotherms are curated and used to extract the three properties. Atom and bond features of surfactant molecules are incorporated as model inputs; MD-computed atom features and molecular descriptors are also introduced to test whether these physics-inspired features help improve the model. The full feature vectors in the model are updated through the D-MPNN algorithm and aggregated into the molecular embedding vector, which is connected to fully connected neural network layers to predict the three surfactant properties.

Our model predicts the three surfactant properties reasonably well: the average R² is 0.87 and the RMSE is 1.07, 2.64, and 0.28 for ${\Gamma }_{max}$, ${\gamma }_{cmc}$, and CMC, respectively. The model performance is only marginally improved by introducing MD-computed features of surfactant molecules. We also systematically examine the variation of predicted properties for surfactants with different types and sizes of head- and tail groups. With a few exceptions, within the chemical space explored, the trends (and, sometimes, the scaling law) of the three properties predicted by the model often agree well with those established theoretically and experimentally in the literature, despite that little experimental data regarding the explored chemical space is available in the training dataset and that those trends are not encoded directly into the model.

The fact that our QSPR model was trained on a small and noisy dataset in the vast chemical space of surfactants not only describes the dataset reasonably well but also captures the variation of surfactant properties in the chemical space is encouraging. The quality of the model can be improved by expanding the training dataset. In this regard, analysis of the trends predicted by the current model can help researchers gather new data strategically. Indeed, it would be helpful to experimentally explore the chemical space in which the current model predicts unusual trends or trends contradictory to expectations based on conventional ‘chemical’ wisdom or data from related surfactants. Such experiments help explore potentially new surfactant behaviors to expand our chemical knowledge, and their data will improve the reliability of the QSPR model, both of which will benefit the rational design of surfactants for diverse applications.