Empowering Machines to Think Like Chemists: Unveiling Molecular Structure-Polarity Relationships with Hierarchical Symbolic Regression

Rather than employing molecular fingerprints or conventional physicochemical descriptors like common machine learning does, our approach deliberately labels functional groups, just like human chemists typically do. In particular, this study encompasses an array of solute molecules, including those containing carbonyl, hydroxyl, amino, halogen, nitro, and cyano functionalities. This enables a tailored analysis of how these functional groups influence the compound’s $R_{f}$ value across various solvents in terms of type, number, and allocation. For notational simplicity, we utilize the abbreviations from Table 1 for the features employed in our study. This table also provides the descriptions and statistics of these features. Our method merely uses $19$ distinct features for describing solute molecules, and thereby demonstrates an effective reduction of dimensionality compared with MACCS key molecular fingerprints ($167$ dimensions). Furthermore, in contrast to traditional physicochemical descriptors, which are commonly low dimensional, our features offer a more localized and interpretable representation of the target molecule.

To evaluate the effectiveness of the selected features, we conducted a comparative analysis of their representation power (see Materials and Method for more details). As shown in Fig. 2, our proposed features demonstrate comparable fitting performance to MACCS, and slightly surpass conventional physicochemical descriptors. The results emphasize the substantial representational power of our features.

In this study, the stationary phase on the TLC plate is silica gel. The determined $R_{f}$ values epitomize a dynamic equilibrium. It reflects the competitive interplay between the highly polar stationary phase and the solute molecules, which are mediated through interactive forces by changing the polarities of the mobile phase during the experiment. This interplay exhibits two primary patterns: the solvent facilitates the upward transit of solute via capillary action; concurrently, solvent molecules intervene by dislodging solute molecules from the stationary phase’s surface, attenuating their interactions (as illustrated in Fig. 3).

The two polarity indices distinctly capture the interaction dynamics between the solvent compound and the silica gel, as well as between the solute molecule and the silica gel. Importantly, our approach enables a quantitative assessment of the aforementioned interactions. In the following paragraphs, we will separately analyze these two indices.

Considering that the dataset contains $387$ organic compounds, previous research normally involves high dimensional features and physicochemical descriptors to encode the solute molecule into numerical values. However, the former are often challenging for human understanding, while the latter lack structural information.

To carry out this study, we carefully selected $16$ representative functional groups and used their counts as features for the target molecules. Additionally, we characterized the distribution of these functional groups within the molecule by using three features: count of benzene (NBen), molecular skeleton descriptor (MSD), and dipole moment (DM). Although we were able to reduce the feature numbers to $19$, establishing a explainable connection between the polarity of the molecule and these input features remains a major challenge.

We distinguished two types of solute features: functional group (FG) counts and FG distribution, as they may have distinct effects on the polarity of a given molecule (Fig. 3A). For instance, a molecule with a higher count of polarity-functional groups but with a symmetric structure may exhibit lower polarity than molecules with fewer polarity-functional groups. Therefore, we further extracted two indices, namely FG distribution polarity index $\alpha$ and FG polarity index $\beta$, to isolate these effects.

Despite the polarity of functional groups being of great interest, there is a lack of quantitative and scalable methods to compare the polarities of different functional groups; in current practice, chemists often rely mainly on qualitative analysis and their experience. Our FG polarity index provides the first empirical quantification of the impact of different functional groups on molecular polarity. Considering the $16$ functional groups involved in this study, the input features for the sub-model used to obtain the FG polarity index $\beta$ consist of $16$-dimensional vectors, with each dimension representing the number of the respective functional group. To quantify the polarity of a specific functional group FG_i, we configured the input feature vector such that the $i$-th position was set to $1$ and the others to $0$. As illustrated in Fig. 3B, the sub-model’s outputs reveal a distinct polarity order among the functional groups, i.e., amides ($3.36$) $>$ carboxylic acid ($2.90$) $>$ amine ($2.81$) $>$ alcoholic hydroxyl ($2.21$) $>$ phenolic hydroxyl ($2.05$) $>$ aldehyde ($0.46$) $>$ nitro ($-0.05$) $>$ ketone ($-0.31$) $>$ ester ($-0.72$) $>$ fluorine ($-1.25$)$>$ ether ($-1.59$) $>$ methyl ($-1.74$) $>$ bromine ($-2.04$) $>$ cyano ($-2.20$) $>$ chlorine ($-2.22$) $>$ iodine ($-2.94$).

Another important interaction is between the solvent compound and the silica gel. The solvent polarity index $\Psi$ provides an analytical lens through which we are able to assess the influence of disparate solvents on the $R_{f}$ values.

Methanol (MeOH) possesses one hydrogen bond acceptor (HBA) and one hydrogen bond donor (HBD), capable of forming stable six- and seven-membered chelates with the silica gel. It also boasts an HBA that can engage in intermolecular hydrogen bonding with the siloxane (Si-O-Si) bridges of the silica gel. Ethyl acetate (EA) harbors two HBAs that can establish an eight-membered ring complex with the silica gel. The oxygen atoms of the ester moiety can partake in two distinct types of intermolecular hydrogen bonds with the gel. Diethyl ether (Et₂O) contains one HBA and is thus limited to a singular mode of intermolecular hydrogen bonding with the silica gel. Dichloromethane (DCM) is characterized by two chlorine atoms, which can form a relatively less stable eight-membered ring structure with the gel, and induce weak intermolecular hydrogen bonding. The interaction between n-hexane (Hex) molecules and the silica gel is considerably feeble, with negligible capacity to displace solute molecules. Consequently, under the paradigm of “like dissolves like”, we infer a polarity hierarchy as follows:

The input features of the sub-model used to obtain $\Psi$ consist of five-dimensional vectors, with each dimension representing the volume percentage of a solvent compound. Consequently, variations in the compound ratio within the solvent system directly influence the value of $\Psi$, thereby changing the polarity of the solvent. Moreover, to isolate the effect of each compound, we can assign the value $1$ to the for the volume percentage that corresponds to the target solvent compound, while setting the others to $0$. This approach yields the following values (see Fig. 3B):

As shown in Table 2, we also conducted a comparative analysis for our empirical solvent descriptor $\Psi$ and experimental miscibility indices such as LogP and aqueous solubility (?). The observed numerical consistency and trends align with our preliminary postulates (see Table 2).

Organic molecules are composed of various functional groups arranged along their molecular skeletons. The arrangement of these groups influences the polarity of a given molecule. Organic chemists typically assess the polarity of target molecules by considering the type, number, and distribution of functional groups. Here, we focused on $16$ representative functional groups (Fig. 1C). Besides the type and number of functional groups, their spatial distribution is also crucial for molecular polarity. To account for this, we included additional descriptors such as the number of aromatic rings, dipole moments, and molecular skeleton descriptors (detailed in the Supplementary Materials).

In terms of features related to the solvent, we considered the volume ratios ($V\%$) of five different solvent compounds, i.e., n-Hexane, Ethyl acetate, Dichloromethane, Diethyl ether and Methanol.

To evaluate the representation power of the selected features, we conducted a comparative analysis. This analysis employed the XGBoost model and examined three different sets of features. Hereto, like in previous studies  (?), we primarily focus on solute molecular features. To maintain consistency, the features related to the eluent solvent were kept the same across all experiments. Feature group A consists of $19$-dimensional solute features (Table 1). Feature group B encompasses $167$-dimensional features, containing $167$-dimensional MACCS keys for describing molecular structure. Feature group C contains $8$-dimensional features, focusing on conventional physicochemical descriptors. Specifically, these descriptors include molecular weight (MW), topological polar surface area (TPSA), the logarithm of $n$-octanol–water partition coefficient (LogP), hydrogen bond donor number (HBD), hydrogen bond acceptor number (HBA), molecular polarity index (MPI), polar surface area percentage (PSA) and dipole moment (DM).

The architecture of our network is motivated by the rationale grounded in established practices observed in TLC experiments. TLC experiments conventionally entail the execution of controlled experiments designed to meticulously dissect the influence of solute polarity and solvent polarity on $R_{f}$. In this context, we modified conventional Multi-Layer Perceptron (MLP) architecture to meet the specific demand of our task. Specifically, given an input sample $\mathbf{x}\in\mathbf{R}^{d}$, the initial step involves the classification of all features into $k$ clusters, denoted as $\mathbf{x}_{1},\dots,\mathbf{x}_{k}$. Here, $\mathbf{x}_{i}\in\mathbf{R}^{d_{i}}$ and their dimensions sum up to $d$, i.e., $\sum_{i=1}^{k}d_{i}=d$. To enhance interpretability of these clusters, we incorporated a chemists-guided strategy, and the clusters used at each stage are demonstrated in Fig. 5.

Then, a key innovation lies in the formulation of the target latent representation layer $j$, where $z^{(j)}_{i}=g_{i}(\mathbf{x}_{i})$. Each function $g_{i}$ operates independently in the modified MLP, in the sense that there is no parameter sharing across $g_{1},\dots,g_{k}$. Thus, we can consider $g_{i}$ as a sub-model. This design ensures the preservation of feature-specific information. The final output $o=h(\mathbf{z}^{(j)})$ encapsulates the information from the target latent representation layer.

This design restricts the receptive field of the latent variables, therefore ensuring that each latent variable solely depends on a specific group of input features. For instance, in the first stage, we considered all features and categorized them into two distinct clusters, solvent features and solute features. This categorization aligns with practical considerations where the polarity of the solvent and the solute is often analyzed separately. Consequently, the learned latent variables in this layer can be regarded as the solvent polarity index and the solute polarity index, respectively. The connections between the latent variables and the output layer capture the nonlinear relationship inherent in the data. This approach enables the model to learn informative latent variables in an unsupervised manner, which disentangles and highlights the distinct contributions of different groups of input features.

SR aims to identify mathematical expressions that capture the inherent relationships within a given dataset. Unlike approaches that involve fitting parameters to overly complex general models (such as machine learning models), SR explores a space of simple analytic expressions. The goal is to discover accurate and interpretable models that directly represent the underlying patterns in the data. SR is usually implemented by evolutionary algorithms, such as genetic programming (GP)  (?, ?, ?).

The most common way to visualize a symbolic expression is a tree-structure with nodes and branches, as shown in Fig. 4A. This structure includes primitive functions (e.g., $+,-,\times,\div,\exp)$ and terminal nodes (input features and numeric constants). Through a series of mutation operations, the GP algorithm seeks to determine the optimal number of nodes and terminals that provide the best fit to a given dataset. For our implementation of the widely used GP-SR algorithm, we employed the open-source Python library PySR¹¹1https://github.com/MilesCranmer/PySR  (?). The hyperparameter configurations for our implementation are detailed in the Supplementary Materials.

It is worth mentioning that our method is not limited to GP algorithms. In fact, the UHiSR framework is designed to accommodate various symbolic regression methods; comparable results obtained from deep reinforcement learning based SR methods (such as DISCOVER  (?)) are included in the Supplementary Materials.