Unsupervised Hierarchical Symbolic Regression for Interpretable Property Modeling in Complex Multi-Variable Systems

中文 / Bilingual · English version

Siyu Lou ${}^{1,2}$ D | Chengchun Liu ${}^{3}$ D | Dongxiao Zhang ${}^{2}$ | Yuntian Chen ${}^{2}$ D | Fanyang Mo ${}^{3,4,5,6}$

${}^{1}$ School of computer science, Shanghai Jiao Tong University, Shanghai, P.R. China ${}^{1}{}^{2}$ Ningbo Key Laboratory of Advanced Manufacturing Simulation, Eastern Institute of Technology, Ningbo, P.R. China I ${}^{3}$ School of Materials Science and Engineering, Peking University, Beijing, P.R. China I ${}^{4}$ School of Advanced Materials, Peking University Shenzhen Graduate School, Shenzhen, P.R. China I ${}^{5}$ AI for Science (AI4S)-Preferred Program, Peking University Shenzhen Graduate School, Shenzhen, P.R. China I ${}^{6}$ Guangdong Provincial Key Laboratory of Nano-Micro Materials Research, Peking University Shenzhen Graduate School, Shenzhen, P.R. China

Correspondence: Yuntian Chen () | Fanyang Mo ()

Received: 23 October 2025 | Revised: 15 December 2025 | Accepted: 19 December 2025

Keywords: explainable AI | molecular polarity | molecular structure | symbolic regression | TLC

ABSTRACT

Recent advances in artificial intelligence have produced powerful predictive models for chemical analysis, but model interpretability remains a challenge. Here, we introduce Unsupervised Hierarchical Symbolic Regression (UHSR), providing an explainable solution while maintaining competitive predictive performance. With a focus on thin-layer chromatography (TLC), a crucial technique in molecular polarity analysis, UHSR automatically distills chemically-intuitive retention indices and discovers explainable equations that link molecular structures to chromatographic behavior. Experiments have shown UHSR’s capability to derive concise and accurate governing equations linking polarity to molecular structures from the TLC dataset. A survey of 100 expert chemists demonstrates that our UHSR model gains more trust from chemists compared to traditional models. Additionally, we also show its adaptability to other property prediction tasks beyond molecular polarity.

1 | Introduction

Scientific discoveries aim to explore and understand the world through the observation of phenomena, the systematic collection of data, the identification of patterns, and the elucidation of underlying principles that enable the prediction of future events [1, 2]. The derivation of empirical equations and parameters through data-driven methods has been pivotal in elucidating chemical phenomena and predicting experimental outcomes. For example, Henry’s Law $\left( {C = {k}_{H} \cdot P}\right)$ for solubility predictions [3], the Arrhenius equation $\left( {k = A \cdot {e}^{-\frac{{E}_{a}}{RT}}}\right)$ for reaction rate dependence on temperature [4], and the van’t Hoff equation

This is an open access article under the terms of the Creative Commons Attribution License, which cited. $\left( {\ln K = - \frac{\Delta {H}^{ \circ }}{RT} + \frac{\Delta {S}^{ \circ }}{R}}\right)$ for equilibrium constants related to temperature [5]. These empirical equations and parameters, derived through extensive experimental work and large datasets, remain foundational in contemporary chemistry.

Polarity, reflecting the uniformity and symmetry of charge distributions in molecules, is crucial for understanding molecular interactions and characterizing molecules $\left\lbrack {6,7}\right\rbrack$ . Chromatography serves as an essential tool for studying polarity by providing insights into how molecules interact under various chemical environments. Over the years, various chromatographic models have been developed to describe retention behavior and polarity-driven interactions. For example, linear free energy relationship models [8, 9] have focused on analyte-stationary phase interactions, capturing how polarity influences adsorption and desorption processes. Similarly, solvent strength models [10, 11] have been widely used in reversed-phase liquid chromatography. These models are explicitly designed for binary aqueous-organic systems, where polarity differences between water and organic solvents are key drivers of retention behavior.

Siyu Lou and ChengChun Liu contributed equally to this work.

Thin-layer chromatography (TLC) provides critical insights into the behavior of organic molecules in various organic solvent environments. Polarity is a key factor in TLC, as the interactions between solutes, solvents, and the stationary phase are largely governed by differences in polar and nonpolar characteristics. However, TLC is usually labor-intensive and involves many repetitive trials [12, 13]. Despite its importance in determining the retardation factor $\left( {R}_{f}\right)$ , previous prediction models for structure-retardation factors have been limited by the quantity and standardization of TLC data [14, 15]. To address this challenge, our prior work introduced an automatic high-throughput platform for TLC analysis, providing a more extensive and standardized dataset for model training [16, 17]. This enhanced dataset not only improves the performance of predictive models, enabling more accurate and robust predictions of structure-retardation factors in TLC experiments, but also facilitates a better analysis of the liquid-liquid interactions between solutes and solvents.

In recent years, significant progress has been made in artificial intelligence (AI)-assisted chemical analysis, especially in the development of quantitative structure-activity relationships (QSAR) prediction models [18-22] and quantitative structure-property relationships (QSPR) prediction models [23-25]. While these models have achieved remarkable prediction accuracy, the AI models and feature representations they employed often lack interpretability, hindering a deeper understanding of the underlying relationships within the dataset.

The lack of transparency not only restricts the scientific impact of these models, but also undermines researchers’ ability to fully trust and understand their predictions [26-28]. To address this concern, there is an increasing need for (1) transparent models and (2) interpretable feature representations that not only achieve high performance but also align with human reasoning and provide meaningful insights.

One strategy to build transparent models is to articulate explicit formulas between the ${R}_{f}$ values and the solute molecular structures as well as the eluent solvents in TLC analysis. Symbolic regression (SR) has recently emerged as an effective approach for unveiling complex relationships within datasets [29-33], and has significantly contributed to many scientific discoveries [34- 40]. However, applying symbolic regression to datasets with a large number of variables remains challenging, as increasing the number of variables often leads to decreased performance and overly complex formulas, which undermines interpretability. For example, as shown in Supporting Information, when symbolic regression is directly applied to commonly used feature sets, such as physicochemical descriptors, the resulting models have limited predictive performance and reduced interpretability due to the high dimensionality. In terms of feature representations, molecular fingerprints, such as MACCS key molecular fingerprints [41], have been often used in AI models. However, these representations are typically high-dimensional and lack interpretability, making it difficult to understand how specific molecular structures influence the model prediction.

In this work, we present unsupervised hierarchical symbolic regression (UHSR) (see Figure 1), a new SR approach guided by modular neural networks. UHSR introduces novel retention indices, e.g., solvent retention index and solute retention index, which are learned through a modular neural network. Guided by domain knowledge, this model comprises multiple sub-models that govern the mapping relationships between a set of input variables, e.g., solvent compositions, and a specific retention index, e.g., solvent retention index, as well as the mapping between retention indices and the output variable ${R}_{f}$ . It is noteworthy that this process aligns with the chemists’ thinking processes. The human brain struggles with high-dimensional information, opting not to directly comprehend intricate mappings but rather to decompose and analyze through individual sub-models. Then, we can apply the SR method on the output ${R}_{f}$ and the retention indices, effectively reducing the dimension of the input variables of the SR algorithm. Our experiments on the TLC benchmark dataset [17] demonstrate that UHSR achieves competitive prediction performance comparable to DNNs while providing concise formulas linking the ${R}_{f}$ value (indicative of the molecular polarity) with molecular structure features. Importantly, the retention indices and governing equations are evaluated by expert chemists, highlighting the alignment between the decision-making process of UHSR and that of chemists. Additionally, we extend the scope of our study to other chemical property prediction tasks, such as spectra-property prediction [42], where our framework also demonstrates high predictive accuracy and interpretability.

2 | Results

2.1 | Data-Driven Explainable Retention Indices

It is essential for chemists to accurately characterize the overall polarity of molecules, as polarity influences molecular interaction, solubility, and chromatographic behavior. For this purpose, there are many known descriptors, such as the molecular polarity index (MPI), which is based on theoretical calculations. While effective in some cases, these descriptors have some limitations. For example, they are usually computationally expensive, especially for large or complex molecules. Also, they typically provide only global polarity information, failing to capture the localized effects of individual functional groups or their specific contributions to intermolecular interactions. Another way to estimate the molecule’s polarity is based on TLC experiments. This technique directly measures the ${R}_{f}$ value, which reflects how a molecule interacts with both the stationary and the mobile phase under specific experimental conditions. TLC experiments are relatively simple, low-cost, and widely accessible compared to computational methods, which often require advanced software, high-performance computing resources, and expert knowledge. By varying the solute or solvent composition and evaluating the resulting ${R}_{f}$ values, chemists can systematically study polarity in a way that is directly applicable to experimental workflows. Building on the strengths of TLC as a practical and experimentally grounded approach to studying molecular polarity, we identified several limitations in traditional methods of analyzing TLC data, which motivated this study. While TLC experiments provide valuable ${R}_{f}$ values that reflect molecular polarity, these raw experimental outputs are often difficult to interpret quantitatively and lack the ability to systematically separate solvent and solute contributions to retention behavior. Therefore, we developed a data-driven method to further refine the analysis. Using the TLC dataset, we extract two retention indices, namely $\Psi$ and $\xi$ , acting as empirical descriptors providing information about the polarity of the solvent system and the solute molecule, respectively. The clear chemical differentiation between solutes and solvents, alongside their near-zero Spearman correlation, compellingly validates our statistical methodology, as detailed in the Figure S1.

FIGURE 1 | Overview of Unsupervised Hierarchical Symbolic Regression (UHSR). (A) Illustration of TLC experiment and the calculation of the retardation factor $\left( {R}_{f}\right)$ . (B) Feature engineering, involving five solvent features based on volume percentages and the decomposition of target molecules into functional groups (FGs). The molecular structure is treated as a composite formed by stacking various functional group modules. (C) UHSR framework with three main stages: chemist-guided feature grouping, modular neural network for latent variable extraction (e.g., solute retention index $\xi$ ), and symbolic regression for discovering explicit equations between the target value and latent variables.

Importantly, compared to traditional physicochemical descriptors, our proposed indices are directly rooted in the experimental TLC dataset, which ensures that these indices are reflective of real-world molecular interactions. Furthermore, $\xi$ provides a localized and detailed perspective on molecular polarity, e.g., how specific functional groups impact the polarity through their specific contributions to solute-solvent and solute-stationary phase interactions. This detailed resolution allows $\xi$ to capture subtle variations in polarity that are often overlooked or averaged out by traditional global descriptors.

TABLE 1 | Description and statistical analysis of the features used in this study. The abbreviations, descriptions, mean values, along with the standard deviations (Std.), are presented. Additionally, Pearson’s correlation coefficients (r) are computed to quantify the linear relationship between each feature and the observed ${R}_{f}$ value.

2.1.1 | Chemically-Intuitive Feature Engineering

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. This makes our features inherently interpretable and accessible to experimental chemists. Specifically, 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 achieves a significant dimensionality reduction, requiring only 18 distinct features for describing solute molecules, compared to MACCS key molecular fingerprints (167 dimensions). While traditional physicochemical descriptors are often low-dimensional, they typically lack the localized and explainable representation that functional group-based features provide. Furthermore, functional groups are well-understood chemical entities, making our results interpretable and actionable for researchers, including entry-level chemists.

It is worth noting that an equation could be derived from various types of descriptors, and alternative descriptors beyond functional group counts, physicochemical descriptors, could also potentially yield high accuracy. However, these alternative descriptors often require specialized knowledge or computational resources, making them less accessible and practical for experimental researchers who need actionable insights. In contrast, our features are straightforward and directly relevant to experimental chemists, therefore, they bridges the gap between computational chemistry and experimental practice.

FIGURE 2 | Illustration of the retention indices and their impact on chromatographic behavior. (A) The input features correspond to different retention indices. Here, FG stands for functional group. (B) TLC experiment can be understood as a process where the stationary phase (silica gel) and the mobile phase (solvent) compete for the solute molecules. This competition’s outcome is reflected in the ${R}_{f}$ value. The factors influencing the ${R}_{f}$ value can be categorized into three types: interaction between the stationary phase and the solvent (I), interaction between the stationary phase and the solute (II), and interaction between the solvent and the solute (III). The two retention indices $\left( {\Psi \text{ and }\beta }\right)$ characterize how solvent and solute impact the chromatographic behavior separately. (C) Illustration of interactions between the stationary phase and different functional groups.

2.1.2 | Retention Indices

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 Figure 2). The two retention 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.

2.1.2.1 | Solute Retention Index $\xi$ . Considering that the dataset contains 387 organic compounds, previous research normally involves high-dimensional features (such as MACCS keys with 167 dimensions) and physicochemical descriptors to encode the solute molecule into numerical values. However, the former are often challenging for human understanding, while the latter fails to meet the needs of experimental chemists for rapid polarity evaluation of new compounds. We also conducted a comparison analysis between the solute retention index $\xi$ and several existing physicochemical descriptors, demonstrating the superiority of the solute retention index $\xi$ in the context of TLC experiments.

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 the number of benzene rings (NBen) and dipole moment (DM). Notably, by further refining the allocation of functional groups within the benzene ring, we were able to improve prediction accuracies; the corresponding results are presented in Supporting Information.

Although we were able to reduce the feature numbers to 18, establishing an explainable connection between the polarity of the molecule and these input features remain a major challenge.

Therefore, we distinguish between two essential solute characteristics: functional group (FG) counts and their spatial distribution (Figure 2A). Both aspects are critical for elucidating molecular properties, particularly polarity, which profoundly influences molecular behavior in various solvents and reaction environments. FG counts involve 16 predefined structural motifs, offering insights into potential reactivity and solubility. While these counts provide an initial understanding of chemical behavior, molecular properties are not solely additive; linear contributions from groups alone may not fully capture the complexity of chemical phenomena. Thus, the spatial arrangement of these groups is pivotal, as it impacts the molecule’s geometry, interactions, and stereochemical properties, leading to a more comprehensive characterization of molecular polarity. Therefore, we further extracted two indices, namely the FG retention index $\beta$ and the FG distribution retention index $\alpha$ , to isolate these two 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 retention 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 retention index $\beta$ consist of ${16}\mathrm{D}$ vectors, with each dimension representing the number of the respective functional group. To quantify the polarity of a specific functional group ${\mathrm{{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 Figure 2B, the sub-model’s outputs reveal a distinct polarity order among the functional groups, i.e., amides (1.43) > carboxylic acid (0.88) > amine (0.80) > alcoholic hydroxyl (0.75) > phenolic hydroxyl $\left( {0.43}\right) >$ aldehyde $\left( {-{0.20}}\right) >$ ketone $\left( {-{0.20}}\right) >$ ester $\left( {-{0.48}}\right) >$ nitro $\left( {-{0.50}}\right) >$ cyano $\left( {-{1.03}}\right) >$ fluorine $\left( {-{1.03}}\right) >$ chlorine $\left( {-{1.15}}\right) >$ ether (-1.22) > methyl (-1.28) > bromine (-1.55) > iodine (-2.97).

To further validate our proposed solute retention index $\xi$ , we conducted a Spearman correlation analysis with established molecular properties. The solute retention index exhibited a moderate positive correlation with $n$ -octanolwater partition coefficient (LogP, ${r}_{s} = {0.50}$ ), indicating that as the LogP values increase, there is a general trend of increasing the solute retention index $\xi$ , suggesting that our index partially captures hydrophobic characteristics. Conversely, the solute retention index $\xi$ showed a moderate negative correlation with topological polar surface area (TPSA, ${r}_{s} = - {0.61}$ ), highlighting that as the TPSA increases, solute retention index $\xi$ tends to decrease, effectively capturing aspects of molecular polarity related to surface area.

For a comprehensive comparison, we employed our solute retention index $\xi$ , MPI, LogP, TPSA, and DM as 1D solute characteristics, using five different solvent volume percentages as 5D solvent characteristics in a model established through the XGBoost model. This analysis is entirely independent of our UHSR model and was conducted to evaluate the predictive capability of $\xi$ in a separate machine learning context. The XGBoost model served as a robust benchmark for evaluating the predictive capability of different solute characteristics.

Specifically, we compared the predictive performance of the XGBoost model, as shown in Table 2a. The solute retention index $\xi$ demonstrated superior predictive performance, outperforming other solute descriptors. We also performed tree-based feature importance analyses within the XGBoost model, as detailed in Figure S1. These analyses highlighted the predominant influence of the solute retention index $\xi$ over other traditional descriptors.

2.1.2.2 | Solvent Retention Index Ψ. Chromatographic retention on silica gel results from the synergistic interplay of electrostatic (dipole-dipole and ion-dipole), hydrogen-bonding, and van der Waals dispersion forces; molecular polarity alone represents only one component of this adsorption mechanism. Accordingly, we define the solvent retention index $\Psi$ as an empirical measure of the net adsorption strength of each pure solvent on silica.

Methanol (MeOH), with one hydrogen-bond donor (HBD) and one acceptor (HBA), forms both six- and eight-membered chelates with the surface silanol and siloxane groups (see Figure 2C). Ethyl acetate (EA) features two HBAs that engage in an eight-membered ring complex via its ester oxygens. Diethyl ether $\left( {{\mathrm{{Et}}}_{2}\mathrm{O}}\right)$ has a single HBA, allowing only one hydrogen-bonding mode. Dichloromethane (DCM), bearing two chlorides, interacts primarily through weaker dipole-dipole and dispersion forces. In contrast, $n$ -hexane (Hex) adsorbs only via weak dispersion interactions and cannot displace more strongly adsorbing solvents. Under this framework, the adsorption hierarchy is

The input features of the sub-model used to obtain $\Psi$ consist of 5D 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 volume percentage that corresponds to the target solvent compound, while setting the others to 0 . This approach yields the following values (see Figure 2B):

Comparison with experimental miscibility indices (LogP and aqueous solubility [43]) confirms that $\Psi$ captures the combined electrostatic, hydrogen-bonding, and dispersion contributions to solvent-silica adsorption (see Table 2b).

TABLE 2 | Comparison of (a) solute retention index $\xi$ and (b) solvent retention index $\Psi$ with established physicochemical descriptors.

(b) Comparison of solvent retention index $\Psi$ with estabilished solvent properties

2.2 | The ${R}_{f}$ Governing Equation

To comprehensively understand the relationship between the ${R}_{f}$ value and the two retention indices, we used symbolic regression to formulate empirical mathematical expressions connecting the ${R}_{f}$ value and the two retention indices (Figure 3A). The observed ${R}_{f}$ values were obtained directly from our automatic high-throughput platform [17], while the computed ${R}_{f}$ values were calculated from the ${R}_{f}$ governing equation (see Equation 1 below).

To ensure that the computed ${R}_{f}$ values lie within the range from 0 to 1, a sigmoidal function was employed (see Figure 3B). By using symbolic regression, one obtains several candidates for the ${R}_{f}$ governing equation (see Table S1). Since a simple formula is desired, and considering the prescribed fitting accuracy $\left( {{R}^{2} = }\right. \; {0.908},\mathrm{{RMSE}} = {0.104})$ , we selected

In order to separately analyze the individual influence of solvent and solute on the ${R}_{f}$ value, we further decomposed the formula for ${R}_{f}$ in Equation 1 into $h\left( \Psi \right)$ and $g\left( \xi \right)$ . For instance, as illustrated in Figure 3C, when considering a specific solute molecule (i.e., fixing a value for the solute retention index $\xi$ ), we can reformulate Equation 1 as ${R}_{f} = h\left( \Psi \right) = \sigma \left( {{5.15\Psi } + {C}_{1}}\right)$ , where ${C}_{1}$ depends on the chosen solute molecule. The pattern observed in Figure 3C shows that ${R}_{f}$ increases with the solvent retention index $\Psi$ . Since a larger value of solvent retention index $\Psi$ indicates greater polarity, this finding emphasizes that the solvent retention index $\Psi$ correlates with the polarity of the solvent. Furthermore, we observe that the shape of $h\left( \Psi \right)$ ’s graph varies significantly for different solute molecules. In particular, Figure 3C indicates that, within the context of our dataset, changing the solvent’s composition only has a subtle effect on ${R}_{f}$ for solute molecules with extreme polarities. Similarly, when the solvent is fixed, Equation 1 can be reformulated as ${R}_{f} = g\left( \xi \right) = \sigma \left( {{5.15\xi } + {C}_{2}}\right)$ , where ${C}_{2}$ depends on the chosen value for the solvent retention index $\Psi$ . By analysis similar to the previous case (fixed solute), the ${R}_{f}$ value increases with the solute retention index $\xi$ . But now, a larger value of the solute retention index $\xi$ indicates a lower polarity. Moreover, we observe that the shape of $g\left( \xi \right)$ ’s graph is of the same $S$ -shaped type across different solvent systems (see Figure 3C).

Due to the limitations of the data, different modular neural networks can achieve similar accuracy. However, we observed that while the retention indices change, the discovered ${R}_{f}$ equations maintain a similar form with varying parameters, as illustrated in Tables S10-S12. In TLC experiments, the absolute values of the retention indices are less crucial than their relative values. The relative values of the retention indices are what guide the separation process in TLC. Given that the underlying equations maintain a similar form, variations in the retention indices can be viewed as differences in units rather than changes in fundamental properties. For example, given a target molecule and two solvent systems (e.g., the polarity of the first solvent system is smaller than the second solvent system), we may obtain two sets of solvent retention indices $\left( {{\Psi }_{1},{\Psi }_{2}}\right)$ and $\left( {{\Psi }_{1}^{\prime },{\Psi }_{2}^{\prime }}\right)$ , the relative relationship maintain the same, i.e., ${\Psi }_{1} < {\Psi }_{2},{\Psi }_{1}^{\prime } < {\Psi }_{2}^{\prime }$ .

2.3 | Quantification of Solvent Polarity

In addition to characterizing the polarity of the eluent solvent, the UHSR framework offers an explainable solution for understanding how different compounds influence the solvent retention index $\Psi$ . In this study, we consider five distinct solvent compounds: Hex, EA, DCM, MeOH, and ${\mathrm{{Et}}}_{2}\mathrm{O}$ . To achieve the designed solvent polarity, chemists usually adjust the ratios of these compounds in the mobile phase system. By directly applying the SR algorithm, we obtained the empirical mathematical equation between solvent retention index $\Psi$ and its corresponding input variables, i.e., the volume ratio of $\mathrm{{Hex}},\mathrm{{DCM}},{\mathrm{{Et}}}_{2}\mathrm{O},\mathrm{{EA}}$ , and $\mathrm{{MeOH}}$ in the mobile phase system. The following formula achieved a high ${R}^{2}$ of 0.956 and a low RMSE of 0.069 with a simple form,

(2)

FIGURE 3 | Visualization of the latent variables and the decomposition of the retrieved formula. (A) Visualization of the observed and calculated ${R}_{f}$ values with two retention indices solvent retention index $\Psi$ and solute retention index $\xi$ . (B) Fitting observed ${R}_{f}$ values with a Sigmoid function. (C) Decomposition of Equation 1 into $h\left( \Psi \right)$ and $g\left( \xi \right)$ . (D) Visualization the flow of how input variables contribute to the ${R}_{f}$ value, here $C$ represents a constant. (E) Importance of functional groups in determining the functional group retention index ( $\beta$ ). (F) Importance of solvent components in determining the solvent retention index $\left( \Psi \right)$ .

Specifically, there are three mobile phase systems involved in the dataset, which are Hex/EA, Hex/Et, O, and MeOH/DCM. Equation 2 quantifies the impact of the two compounds on solvent polarity within each system. For instance, in the presence of Hex and EA, the solvent retention index $\Psi$ is calculated as $\Psi = - {1.37}\mathrm{{Hex}} - {0.20}\mathrm{{EA}}$ . As the volume of EA increases and the volume of Hex decreases, the solvent polarity exhibits a corresponding increase. Moreover, the constant term of 4.65 with MeOH indicates that even a slight increase/decrease in the volume of MeOH will lead to a significant increase/decrease in solvent polarity.

2.4 | Quantification of Solute Polarity

As discussed earlier, we introduce two additional functional group-related indices, FG distribution retention index $\alpha$ and FG retention index $\beta$ , to enhance our understanding of the relationship between the molecular polarity and the molecular structure. With the help of symbolic regression, we obtained an explicit mathematical relationship between $\xi$ and $\left( {\alpha ,\beta }\right)$ as follows:

Breaking down Equation 3, the molecular polarity can be understood as a combination of the effect of functional groups’ distribution $\left( {{f}_{1}\left( \alpha \right) = {0.37\alpha }}\right)$ , the contribution of the functional groups within the molecule $\left( {{f}_{2}\left( \beta \right) = - {0.35\beta }}\right)$ and the interaction between these two factors $\left( {{f}_{3}\left( {\alpha ,\beta }\right) = - {\alpha \beta }}\right)$ .

To perform a finer-grain analysis of the individual impacts of $\alpha$ and $\beta$ on molecular polarity, we further derived specific equations that decompose the contributions of each feature, as follows,

$$

• {0.64}\mathrm{{NBen}} + {1.66}) - {0.74}\text{ , } \tag{4a} $$

where

(5)

The above equations are indeed more complex and less directly interpretable than the top-level ${R}_{f}$ equation (Equation 1). This complexity arises from the intricate nature of molecular polarity, which involves nonlinear interactions between functional groups and structural features. While the top-level ${R}_{f}$ equation is designed to provide a concise, human-readable relationship between solute polarity, solvent polarity, and the ${R}_{f}$ value, serving as a conceptual tool for understanding how polarity indices influence ${R}_{f}$ values, the more detailed Equations (4) and (5) serve a different purpose. These equations function as computational tools to calculate latent chemical properties, such as $\alpha$ and $\beta$ , from molecular structures. By incorporating terms like ${\gamma }_{1}$ to ${\gamma }_{5}$ , they enable precise quantification of the contributions of specific functional groups and their interactions, which would otherwise be difficult to disentangle.

Although these detailed equations are not intended for direct human interpretation, their value lies in supporting downstream analyses, such as feature importance assessments and visualizations. These analyses, introduced in the following subsection, help bridge the gap between computational complexity and practical insights, making the results more accessible and actionable for researchers.

2.5 | Visualization of the Discovered Equations

To understand how the discovered equations explain the relationship between input features and the final ${R}_{f}$ value prediction, we utilized two complementary visualization approaches: feature importance visualization and direct value contribution visualization. These methods provide different but synergistic perspectives on how features influence the final ${R}_{f}$ values.

Feature importance quantifies the sensitivity of the model output to variations in individual input features. Through the derived equations, we can directly identify which features have the most significant impact on the ${R}_{f}$ prediction. For example, as shown in Figure 3E, the count of halogen atoms (e.g., iodine) and the carboxylic acid group $\left( {\mathrm{{Ct}}}_{{\mathrm{{CO}}}_{2}\mathrm{H}}\right)$ are identified as features with high importance. As shown in Figure 3F, methanol (MeOH) plays a significant role in influencing the ${R}_{f}$ values. This approach is particularly useful for identifying the key drivers of variability in the model’s predictions across the dataset. In addition to understanding feature importance, we also visualized the direct numerical contribution of each input variable to a specific prediction. This approach traces how the “value” of each feature propagates through the discovered equations to influence the final ${R}_{f}$ value. As shown in Figure 3D, this step-by-step visualization reveals how input variables contribute to intermediate indices (e.g., ${\gamma }_{1},\beta ,\xi$ ) and ultimately to ${R}_{f}$ .

For instance, if the count of iodine is set to 1 , we can directly observe how this value influences ${\gamma }_{5}$ , which then propagates to $\beta$ , and finally contributes to the ${R}_{f}$ value. Unlike feature importance, which focuses on sensitivity to changes, this approach provides a numerical breakdown of how much each input feature contributes to the final prediction for a given input configuration.

While feature importance highlights which features are most critical across the dataset, direct value contribution visualization, provides a more detailed view of how the equations process specific inputs. Together, these visualizations offer a comprehensive understanding of the discovered equations, allowing us to both identify key features and trace their influence on the ${R}_{f}$ value.

2.6 | Chemists’ Evaluation on UHSR

2.6.1 | Survey Evaluation

To thoroughly assess whether the UHSR model thinks like chemists, we grounded our evaluation using the judgments of chemists. We conducted a comprehensive survey involving 100 chemists with an average of 10 years of experience in the field of chemistry. Please see the complete survey in the Supporting Information.

According to the survey results, an impressive 97.0% of the chemists indicated that the retention indices derived from the UHSR model align well with their expectations based on their professional experience (see Figure 4A). This high level of alignment shows the model’s ability to capture key aspects of molecular polarity that are consistent with established chemical principles. Similarly, 99.0% of the chemists agreed that the ${R}_{f}$ governing equation meets their expectations. This suggests that the model’s predictions are not only accurate but also resonate with the intuitive understanding of chemists. 99.0% of the respondents felt that the retention indices and the governing equation provided new insights and knowledge. These results indicate that the qualitative concepts chemists had accumulated over time align with the quantitative equations we discovered. Our quantitative equations provide a structured and solidified description of the intuitive conclusions that chemists have developed through their experiences. This solidification of qualitative concepts into quantitative forms offers a valuable bridge between empirical knowledge and rigorous scientific modeling.

To better understand chemists’ preference between the UHSR model and traditional machine learning models (such as artificial neural networks), we categorized the chemists into three levels based on their levels of expertise in AI knowledge [44]: Level 1 indicates minimal AI knowledge, Level 2 indicates moderate AI knowledge, and Level 3 indicates extensive AI knowledge. As shown in Figure 4A, among chemists with Level 3 AI knowledge, 88.2% of chemists expressed a higher level of trust in the UHSR model compared to traditional machine learning models, and the remaining showed no preference between the two types of models. Among chemists with Level 2 AI knowledge, 44.6% trusted the UHSR model more, 43.1% had no preference, and only 12.3% trusted the black box model more. For chemists with minimal AI knowledge, 55.6% of the chemists showed no preference, while 22.2% trusted the UHSR model more. This finding suggests that the transparency and interpretability of the UHSR model are particularly appreciated by those with a deeper understanding of AI.

2.6.2 | Experimental Evaluation

Furthermore, it is crucial to recommend solvent ratios for column chromatography by examining the positioning of compounds in TLC experiments under various solvent ratios. Ideally, the solvent ratio should maximize the difference in ${R}_{f}$ values between the components to be separated, or place the target molecule’s ${R}_{f}$ value between 0.2 and 0.3 [45]. Based on this approach, we explored the consistency of the optimal solvent ratio between the UHSR model (an AI model trained from standardized data) and experimental data under non-standardized experimental conditions (such as different TLC plates, temperature, humidity, and time). Here, we examined the product o-hydroxybenzaldehyde 2 and $p$ -hydroxybenzaldehyde 3 of the Reimer-Tiemann reaction [46] from phenol 1 were detected by thin-layer chromatography (Figure 4B). At different solvent ratios, there was a high correlation between ${R}_{f}$ values of the manual TLC experiment (Exp.) and predicted values of UHSR. By analyzing the ${R}_{f}$ differences between component pairs across 19 different solvent ratios for five solvents, we found that 7 out of 9 optimal solvent ratios recommended by UHSR predictions matched the experimental data, demonstrating the robustness of our model under a variety of experimental conditions. Additionally, the quantification of solvent polarity for the recommended optimal solvent ratios was determined using the proposed solvent retention index $\Psi$ . It was observed that some of these optimal solvent ratios exhibit similar polarity.

2.7 | Transferability and Broader Applications of the UHSR Framework

2.7.1 | Extension to Column Chromatography

Column chromatography (CC) can be viewed as a natural preparative-scale extension of the structure-polarity principles that govern TLC behavior. To evaluate the effectiveness of our UHSR framework in CC, we employ the same chemically interpretable molecular features listed in Table 1 (e.g. functional group counts) alongside solvent composition parameters (volume percentages of mobile phase components), utilizing experimental data from an automated CC platform [47].

In CC optimization, two critical volume metrics are defined: ${v}_{1}$ , the eluent volume consumed when the target compound begins eluting (corresponding to retention factor $k$ ), and ${v}_{2}$ , the total eluent volume required for complete compound elution (incorporating peak broadening effects). These volumes are predicted through polarity-governing equations:

$$

• {1.94}\frac{\mu }{\theta } - {25.5}\frac{1}{\theta } + {25.1} $$

FIGURE 4 | Analysis of the survey and model-experimental data correlation. (A) Survey analysis from chemists with varying levels of AI expertise. (B) Assessment of the consistency between the UHSR model (AI model trained on standardized data) and manual experimental (Exp.) data (nonstandardized experimental conditions, such as different TLC plates, temperatures, humidity levels, and time) on the optimal solvent ratio.

where $\theta$ represents the learned solvent polarity index (mobile phase) and $\mu$ represents the learned solute polarity index (analyte).

The equations exhibit competitive predictive performance compared to “black-box” models, such as artificial neural networks (ANN), as reported in [47]. Specifically, ${v}_{1}$ achieves an ${R}^{2} = {0.841}$ and ${v}_{2}$ achieves an ${R}^{2} = {0.834}$ .

2.7.2 | Beyond Polarity Anlysis

Although in this work, we focused on ${R}_{f}$ prediction in TLC as a case study, the UHSR framework is highly transferable and can be applied to a wide range of property prediction tasks beyond polarity analysis. To demonstrate this, we extended the framework to spectra-property prediction [42].

The spectra-property dataset focused on predicting thDCMe adsorption energy $\left( {\Delta {E}_{ads}}\right)$ and charge transfer $\left( {\Delta e}\right)$ of CO-adsorbed Cu-based MOF systems from their infrared and Raman spectra. Specifically, we selected CO@CuBTC dataset, where each sample consists 18 features, including six frequency values $\left( {{f}_{1},\ldots ,{f}_{6}}\right)$ , their corresponding IR values $\left( {{I}_{1},\ldots ,{I}_{6}}\right)$ , and Raman values $\left( {{R}_{1},\ldots ,{R}_{6}}\right)$ .

Using the modular neural network within the UHSR framework, these features were categorized into three indices: frequency index $\left( \widehat{f}\right)$ , IR index $\left( \widehat{I}\right)$ , and Raman index $\left( \widehat{R}\right)$ . The resulting governing equation for adsorption energy $\left( {\Delta {E}_{ads}}\right)$ and charge transfer $\left( {\Delta e}\right)$ are as follows,

These equations achieve impressive performance, with ${R}^{2} = \; {0.862},\mathrm{{RMSE}} = {0.045}$ for adsorption energy, and ${R}^{2} = {0.990}$ , RMSE $= {0.005}$ for charge transfer, respectively. Notably, these results surpass the baseline performance reported in [42], where ${R}^{2} = {0.834},\mathrm{{RMSE}} = {0.055}$ , and ${R}^{2} = {0.938},\mathrm{{RMSE}} = {0.013}$ .

This case study demonstrates the ability of UHSR to adapt seamlessly to different tasks with high performance and interpretability.

3 | Methods

3.1 | TLC Dataset

Several variables contribute to the variation in ${R}_{f}$ values, including factors such as the nature of the stationary phase, temperature, and humidity. Traditional methods of obtaining TLC data manually suffer from a lack of standardization, making scalability challenging. To address this issue, our prior work [17] introduced a robotic platform designed for high-throughput collection of TLC data (Figure 3A), which not only ensured efficiency in data acquisition but also established a standardized and scalable methodology.

The TLC dataset comprised in total 4944 measurements of the ${R}_{f}$ values, involving 387 organic compounds and three mobile phase systems. These mobile phase systems were $n$ -Hexane/Ethyl acetate, Diethyl ether/ $n$ -Hexane, and Methanol/Dichloromethane systems, with a total of 17 different solvent compositions.

3.2 | Framework of UHSR

To address the inherent complexity imbalance between solvent and solute compounds, our framework, unsupervised hierarchical symbolic regression (UHSR), was strategically organized hierarchically in a top-to-down manner. It was motivated by the rationale grounded in established practices observed in TLC experiments. As illustrated in Figure 5, each stage focused on specific aspects of the feature space, allowing the model to learn different latent variables, which discern patterns at different levels of granularity.

For instance, the top stage provided a broad overview of the interaction between solvent and solute, capturing their global polarity characteristics. Subsequent stages then delved into more detailed and specific characteristics of the solute molecular structure, offering intricate insights regarding the contributions of different molecular features to the overall molecular polarity. As illustrated in Figure 1C, the framework contained three key parts: feature engineering, modular neural network and symbolic regressions. In this subsection, we briefly introduce the algorithms and methods, the detailed information is in the Supporting Information.

3.2.1 | Feature Engineering

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 (Figure 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 and dipole moments. The reasons for choosing the other three descriptors are as follows:

FIGURE 5 | Hierarchical structure of learning latent variables. Three stages are organized from top to bottom. (A) At the first stage, two latent variables, $\Psi$ and $\xi$ , are generated to encapsulate the overall polarity of the solvent and the solute, respectively. (B) At the second stage, two additional latent variables, $\alpha$ and $\beta$ , are learned to assess the solute molecule’s polarity from two distinct perspectives. $\alpha$ is linked to the distribution of functional groups, while $\beta$ pertains to the quantity of individual functional groups. (C) The third stage characterizes five latent variables, each representing the impact of specific groups of functional groups.

• DM (Dipole Moment): The dipole moment is a critical descriptor as it reflects the overall polarity of a molecule, indicating how electron density is distributed. However, the dipole moment alone might not fully capture the complexity of molecular interactions, especially in varied solvent environments. Its exclusion does affect the model, but not to the extent that it could singularly determine all polarity-related interactions. We found that combining the dipole moment with other descriptors provided a more robust prediction.

• NBen (Number of Benzene Ring): The presence and number of aromatic rings significantly influence molecular properties and interactions, particularly in aromatic-aromatic stacking and van der Waals interactions. These factors are crucial in many solute-solvent interactions that the dipole moment alone may not account for. In terms of features related to the solvent, we considered the volume ratios $\left( {V\% }\right)$ of five different solvent compounds, i.e., $n$ -hexane, ethyl acetate, dichloromethane, diethyl ether, and methanol.

We conducted an input feature control experiment to evaluate the representation power of various feature sets. Two baseline feature sets were compared: one-hot encoded categorical features and random features. One-hot encoded features represent input features as binary vectors, where each postion indicates the presence or absence of a specific molecule. Random features, on the other hand, are randomly generated numerical values assigned to each molecule, serving as a non-informative baseline. In addition, we tested another two types of features commonly used in statistical models: 167D MACCS keys and 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).

A reaction

component

chemical-intuitive features [total length d=23]

solute (n=387) solvent (n=5)

e.g. FG counts,

[1 .000, 1.408, 0.000, … 0.980, 0.0196, 0.000]

dipole moment

solute molecule features solvent features

random features (k=18) (k=5)

[total length d =23]

k random features

[0.374, 0.950 0.731 … -0.411, 1.021, 0.000]

for each unique molecule

mol. 1 mol.2 mol. 3 sol.3 sol.4 sol.5

one-hot encoding 1 0 0 … 1 0 0 ]

[total length d =392] 0 1 0 … 0 1 0 ]

n features for each 0 0 1 …

0 0 1 ]

molecule

FIGURE 6 | Comparison of input feature control experiments. (A) Schematic of our features vs. random features and one-hot encoded categorical features. (B) Comparison of coefficient of determination $\left( {R}^{2}\right)$ and RMSE values between input features.

The five sets of input features were tested on six different machine learning models, including linear regression (Lin. Reg.), partial least squares (PLS), k-nearest neighbors (kNN), a simple neural network (NN), random forest (RF) and XGBoost. The results are presented in Figure 6.

3.2.2 | Architecture Design of Modular Neural Network

The architecture of our network was motivated by the rationale grounded in established practices observed in TLC experiments. TLC experiments conventionally entailed the execution of controlled experiments designed to meticulously dissect the influence of solute polarity and solvent polarity on ${R}_{f}$ . To meet the specific demands of this task, we modified conventional Multi-Layer Perceptron (MLP) architecture into a modular neural network, designed to enhance interpretability and align with domain knowledge.

Feature Grouping. Let us consider a specific modular neural network. Given an input sample $\mathbf{x} \in {\mathbf{R}}^{d}$ , the initial step involved dividing the input features into $k$ feature groups, denoted as ${\mathbf{x}}_{1},\ldots ,{\mathbf{x}}_{k}$ . Each feature group ${\mathbf{x}}_{i} \in {\mathbf{R}}^{{d}_{i}}$ was non-overlapping, and their combined dimensionality satisfied $\mathop{\sum }\limits_{{i = 1}}^{k}{d}_{i} = d$ . These feature groups were determined using a chemist-guided strategy, and the specific grouping applied at each stage is illustrated in Figure 5. Specifically, in the first stage, we categorized features into two distinct groups: solvent features $\left( {{\mathbf{x}}_{\text{ solvent }} \in {\mathbf{R}}^{5}}\right)$ and solute features $\left( {{\mathbf{x}}_{\text{ solute }} \in {\mathbf{R}}^{18}}\right)$ . In the second stage, the solute features were further categorized into two groups: FG features $\left( {{\mathbf{x}}_{\mathrm{{FG}}} \in {\mathbf{R}}^{16}}\right)$ and FG distribution features $\left( {{\mathbf{x}}_{\text{ FGdist }} \in {\mathbf{R}}^{2}}\right)$ . In the third stage, the FG features were further categorized into five distinct groups: ${\mathbf{x}}_{\mathrm{{FG}}1} \in {\mathbf{R}}^{2},{\mathbf{x}}_{\mathrm{{FG}}2} \in {\mathbf{R}}^{3},{\mathbf{x}}_{\mathrm{{FG}}3} \in {\mathbf{R}}^{2},{\mathbf{x}}_{\mathrm{{FG}}4} \in {\mathbf{R}}^{5}$ and ${\mathbf{x}}_{\mathrm{{FG}}5} \in {\mathbf{R}}^{4}$ .

Latent Representation and Modularity. A key innovation of our architecture lay in the formulation of the target latent representation layer $j$ , where ${z}_{i}^{\left( j\right) } = {g}_{i}\left( {\mathbf{x}}_{i}\right)$ . Each function ${g}_{i}$ operated independently within the modular neural network, meaning there was no parameter sharing across ${g}_{1},\ldots ,{g}_{k}$ . Thus, we could consider ${g}_{i}$ as a sub-model. This design ensured the preservation of feature-specific information. The final output, $o = h\left( {\mathbf{z}}^{\left( j\right) }\right)$ , encapsulated the information from the target latent representation layer.

This design restricted the receptive field of the latent variables, therefore ensuring that each latent variable depended solely on a specific group of input features. For instance, in the first stage, the latent variables learned from solvent features and solute features were interpreted as the solvent retention index $\left( \Psi \right)$ and the solute retention index $\left( \xi \right)$ , respectively. These indices captured distinct contributions from the input features, reflecting their individual effects on ${R}_{f}$ .

Training Procedure. The training procedure consisted of three stages, with each stage involving a modular neural network to extract retention indices that reflected different levels of granularity, as illustrated in Figure 5. In the first stage, the input features were divided into two feature groups: solvent features and solute features. These were processed independently by the modular neural network, resulting in two latent variables: the solvent retention index $\left( \Psi \right)$ and the solute retention index $\left( \xi \right)$ . These indices represented the overall contributions of solvent and solute features to the ${R}_{f}$ value.

In the second stage, we focused on solute features only. The solute features were further divided into two subgroups: functional group (FG) counts and the distribution of functional groups within the molecule. The modular neural network processed these subgroups to extract two additional indices: the FG retention index $\left( \beta \right)$ , which captured the contribution of different FGs to ${R}_{f}$ , and the FG distribution retention index $\left( \alpha \right)$ , which represented the influence of how functional groups were distributed across the molecule. These indices provided a more detailed understanding of how functional groups impact molecular polarity and retention behavior.

The third stage involved a fine-grained analysis of functional group contributions. Here, the FG features were further subdivided into five specific FG categories. Using the modular neural network, five fine-grained FG retention indices $\left( {{\gamma }_{1},\ldots ,{\gamma }_{5}}\right)$ were extracted. These FG retention indices captured the unique contributions of specific types of functional groups to ${R}_{f}$ .

This hierarchical decomposition, progressing from broad solvent and solute features to detailed FG-specific contributions, ensured that the model disentangles complex relationships between molecular features and the ${R}_{f}$ value.

3.2.3 | Symbolic Regression

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 explainable models that directly represent the underlying patterns in the data. SR is usually implemented by evolutionary algorithms, such as genetic programming (GP) $\left\lbrack {{29},{30},{48}}\right\rbrack$ .

The most common way to visualize a symbolic expression is a tree-structure with nodes and branches, as shown in Figure 3A. 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 ${}^{1}$ [49]. The hyperparameter configurations for our implementation are detailed in the Tables S8 and S9. Moreover, when deriving the ${R}_{f}$ governing equations, we imposed a constraint requiring all parameters to be positive. If any parameter turned out to be negative, we added a minus sign before the learned retention index to standardize the results.

It is worth mentioning that our method is not limited to GP algorithms. In fact, the UHSR framework is designed to accommodate various symbolic regression methods; comparable results obtained from deep reinforcement learning based SR methods (such as DISCOVER ${}^{2}$ [50]) are included in the Table S5.

3.3 | Expert Chemist Evaluation

The expert chemist evaluation consists of a quantitative survey assessment in which chemists were provided with a QR code or website link to complete the survey (please see Supporting Information for the complete survey). The study was submitted for ethical review and received approval from the Peking University Institutional Review Board (No. IRB00001052-2408). Participation was not incentivized with compensation, as it was considered to contribute to the general advancement of research. The study was anonymous, and only the lead researcher had access to the assignment of the experimental IDs. Participants were chemical professionals, including faculty members, industry chemists, and graduate students, representing a wide range of experiences and roles in chemistry-related fields (Please see Figure S3 for geographic statistics of the participants). This survey aimed to evaluate the chemists’ opinions on the interpretability of the UHSR model vs. traditional machine learning models. Participants provided their feedback based on their professional experience, and demographic data such as years of experience, educational background, and geographical region were collected to ensure a comprehensive analysis. A total of 101 chemists agreed to participate. One response was excluded due to an invalid answer, leaving 100 valid responses.

3.4 | Experimental Settings and Training Details

When constructing the modular neural network, at each stage, we set the number of neurons in the hidden layers of each sub-model to 50. Furthermore, two hidden layers were employed. To facilitate the discovery of concise equations, we leaned toward employing a relatively simple network structure, which helps to prevent overfitting to experimental noise. The activation function chosen for each sub-model was LeakyReLU. In the first stage, the activation function of the output layer was set to a sigmoid, which ensures a standardized output within the range of 0 and 1 (Figure 3B). The batch size was set to 2048, and the Adam optimizer implemented in PyTorch was applied. The learning rate was 0.01 .

${}^{1}$ https://github.com/MilesCranmer/PySR

${}^{1}$ https://github.com/MilesCranmer/PySR

${}^{2}$ https://github.com/menggedu/DISCOVER

${}^{2}$ https://github.com/menggedu/DISCOVER

To evaluate the performance and generalizability of UHSR, we employed 10-fold cross-validation. The full dataset was randomly partitioned into 10 equal subsets. For each fold, 9 folds were used to train the model, and the remaining fold was used for testing. This process was repeated 10 times, ensuring that each data point was used for both training and testing. Each modular neural network was trained for 1000 epochs, then the best validation checkpoint was selected for testing. We conducted the evaluation in two parts: the modular neural networks to produce retention indices, and the derived ${R}_{f}$ equations.

The performance of the modular neural networks was evaluated in each fold using the $R$ -squared coefficient $\left( {R}^{2}\right)$ and the root mean square error (RMSE), which are defined as

where $N$ is the total number of samples, ${y}_{i}$ and ${\widehat{y}}_{i}$ represent the label and the prediction value for the $i$ -th sample, and $\bar{y}$ denotes the mean of the ground truth values.

The results, presented in Table S2, demonstrate consistently high ${R}^{2}$ values and low RMSE values across all folds, indicating the robustness of the modular neural networks. We also derived the ${R}_{f}$ governing equations of each fold, as shown in Table S3. Most equations share a similar structure. Although the coefficients for $\Psi$ and $\xi$ vary slightly across folds, but generally remain within a consistent range.

4 | Discussion and Conclusion

Chellenges in current approaches to explainability. Recent efforts to establish quantitative structure-polarity relationships primarily rely on machine learning models and the development of informative feature representations [16, 23]. While these approaches achieve high prediction accuracy, a significant limitation lies in their lack of interpretability. This lack of transparency hinders a deeper understanding of the intrinsic relationships within the data, consequently reducing the scientific impact of these models by limiting researchers’ ability to trust and fully comprehend their predictions [26-28].

To address this challenge, explainable methods have been developed, broadly categorized into intrinsic and post-hoc approaches. Post-hoc methods, such as LIME [51] and Shapley values [52], provide insights into the importance of input variables relative to model outputs. However, these methods are limited in their ability to explain the model’s internal mechanisms, offering only a partial understanding of how features interact or contribute to predictions. They essentially provide a narrow view into the “black box,” revealing feature importance without elucidating the relationships that drive the model’s decisions. Conversely, intrinsic models, such as multiple linear regression or partial least squares, are inherently interpretable but often lack the accuracy required to handle complex datasets (Figure 6B).

The UHSR framework introduced in this study offers an alternative solution. By integrating equation-based modeling with hierarchical insights into feature importance. Unlike traditional methods, UHSR not only identifies key features but also explains how these features interact and propagate through the hierarchical structure to influence the final output. Key to this framework are two retention indices, the solvent retention index $\left( \Psi \right)$ and the solute retention index $\left( \xi \right)$ . These indices provide interpretable insights into the structure-polarity relationship while enabling direct ${R}_{f}$ prediction through a formula-driven approach. By leveraging only a few indices (in this case, two) closely tied to the dataset, UHSR enhances interpretability and reduces complexity compared to high-dimensional feature representations. Moreover, the predictive accuracy of UHSR is comparable to that of existing DNN models, demonstrating that interpretability does not need to come at the cost of performance.

Advantages of UHSR beyond interpretability. The symbolic regression component of UHSR offers several notable advantages over traditional machine learning models, particularly for deployment in automated platforms and resource-constrained environments. (1) Deployment efficiency and low computational cost. Although the training cost of UHSR is higher than standard machine learning models due to the additional step of extracting interpretable equations, its key advantage lies in deployment efficiency. Once the symbolic equations are discovered, inference becomes nearly instantaneous, requiring only a few arithmetic operations. In contrast, black-box models like DNNs involve evaluating large numbers of parameters during inference, leading to significantly higher computational demands. This makes UHSR particularly advantageous for automated platforms, where real-time predictions with minimal computational overhead are critical. (2) Robustness and transferability. Equation-based models have fewer parameters compared to traditional machine learning models, resulting in improved robustness and transferability. These models are less prone to overfitting and can adapt to new conditions with minimal retraining. Unlike DNNs, which often require extensive fine-tuning and retraining when applied to new scenarios, symbolic equations derived from UHSR offer cost-effective adaptability while retaining interpretability.

Current limitations and future directions. Despite these strengths, the current study has several limitations, which in turn represent exciting future directions. (1) Dataset Coverage. First, the UHSR model is restricted to the solvent systems represented in the TLC dataset. While the five solvents included are widely used in TLC experiments and adequately represent typical solvent systems (as indicated by 82.0% of surveyed chemists), the current framework cannot address more complex eluents with additives (e.g., AcOH or ${\mathrm{{NEt}}}_{3}$ ) or solvent systems beyond those represented in the dataset. Extending the framework to such systems would require incorporating additional experimental data and retraining the model on an expanded dataset. Second, the diversity of analytes in the training set is also limited, as it predominantly consists of substituted arenes and heteroarenes, which are rigid, planar scaffolds. These analytes’ retention behavior can be effectively captured by the chemically intuitive features developed in this study, such as functional group counts and a single orientation descriptor. However, this limits the model’s applicability to more diverse molecular structures. For example, heterocycles without benzene rings (e.g., thiophene, pyridine) or molecules with multiple benzene rings and fused rings (e.g., naphthalene) would require the introduction of additional chemically intuitive features, such as pyridine nitrogen, thiophene sulfur, fused ring descriptors, or steroid core types. While the current framework is not yet equipped to handle these cases, its flexibility allows for adaptation to new molecular classes once the dataset is expanded. Third, the model’s focus on analytes that are predominantly ${\mathrm{{sp}}}^{2}$ - rich. The chemically intuitive features developed in this study may not generalize well to conformationally flexible, ${\mathrm{{sp}}}^{3}$ -rich molecules, even if they share similar functional groups. Future work will address this limitation by enriching the dataset with diverse 3D scaffolds and developing feature sets that incorporate 3D descriptors (e.g., conformer-ensemble statistics, dihedral angle distributions, and 3D autocorrelation functions). These efforts will help maintain the symbolic regression framework’s interpretability while extending its applicability to a broader range of chemical structures.

(2) Non-Uniqueness of Derived Equations. It is important to note that the retention indices and governing equations derived from the UHSR model are not unique solutions. Due to the inherent limitations of the data, different neural networks (and their corresponding equations) can achieve similar levels of accuracy (see Tables S10-S12 for details). This phenomenon is common in scientific discovery, where multiple models can explain observed phenomena equally well. For example, in reaction kinetics, different rate laws derived from experimental data often describe the same behavior. UHSR is designed to provide one interpretable solution to the problem, aiming to enhance understanding and inspire further research within the field.

(3) Reliance on prior knowledge for feature selection and grouping. The current implementation of the UHSR framework relies heavily on expert-defined feature selection and grouping strategies, which were used to construct the solvent retention index $\left( \Psi \right)$ and solute retention index $\left( \xi \right)$ . While this approach ensures that the derived equations are chemically meaningful and interpretable, it introduces a dependency on prior chemical knowledge. This limits the scalability of the framework, particularly when applied to novel chemical spaces or datasets with less defined chemical properties. Future work will explore more automated approaches to feature selection and grouping to reduce reliance on domain expertise. For instance, unsupervised learning methods such as clustering molecular descriptors could be employed to identify functional groups or suggest feature groupings in a data-driven manner.

Broader implications. While this study focuses on applications in chemistry, the development of explainable frameworks has broader relevance. By improving the interpretability of machine learning models, our framework can enhance trust in predictions and foster collaboration across scientific disciplines, contributing to advancements in data-driven discovery.

Acknowledgements

We thank the National Natural Science Foundation of China (Grant No. 12572266) and the New Generation Artificial Intelligence-National Science and Technology Major Project (grant no. 2025ZD0121904). We also thank the High Performance Computing Platform of Peking University and the High Performance Computing Centers at Eastern Institute of Technology, Ningbo, for machine learning model training. We also thank Mengge Du for assisting in the implementation of DISCOVER. [Correction added on 5-February-2026, after first online publication: New funder is added in the acknowledgment section.]

Ethics Statement Heading

The survey conducted in the study was submitted for ethical review and received approval from the Peking University Institutional Review Board, No. IRB00001052-2408.

Conflicts of Interest

The authors declare no conflicts of interest.

Data Availability Statement

The original code for the experiments and data analysis in this work has been deposited at the website: https://github.com/SiyuLou/ UnsupervisedHierarchicalSymbolicRegression

References

1. H. W. de Regt, Understanding Scientific Understanding (Oxford University Press, 2017).

2. A. Potochnik, Idealization and the Aims of Science (University of Chicago Press, 2017).

3. W. Henry and J. Banks, “Iii. Experiments on the Quantity of Gases Absorbed by Water, at Different Temperatures, and under Different Pressures,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 93 (1803): 29-274.

4. S. Arrhenius, “Über Die Reaktionsgeschwindigkeit Bei Der Inversion Von Rohrzucker Durch Säuren,” Zeitschrift für physikalische Chemie 4, no. 1 (1889): 226-248.

5. J. H. van’t Hoff, Études de Dynamique Chimique (F. Muller & Co., Amsterdam, 1884).

6. G. M. Sessler, Physical Principles of Electrets (Springer, Berlin, Heidelberg, 1987), https://doi.org/10.1007/3540173358_10.

7. A. R. Katritzky, D. C. Fara, H. Yang, K. Tämm, T. Tamm, and M. Karelson, “Quantitative Measures of Solvent Polarity,” Chemical Reviews 104, no. 1 (2004): 175-198.

8. M. H. Abraham, A. Ibrahim, and A. M. Zissimos, “Determination of Sets of Solute Descriptors from Chromatographic Measurements,” Journal of Chromatography A 1037, no. 1-2 (2004): 29-47.

9. C. F. Poole, “Reversed-Phase Liquid Chromatography System Constant Database over an Extended Mobile Phase Composition Range for 25 Siloxane-Bonded Silica-Based Columns,” Journal of Chromatography A 1600 (2019): 112-126.

10. C. F. Poole and N. Lenca, “Reversed-Phase Liquid Chromatography,” in Liquid Chromatography (Elsevier, 2017), pp. 91-123.

11. U. D. Neue and H.-J. Kuss, “Improved Reversed-Phase Gradient Retention Modeling,” Journal of Chromatography A 1217, no. 24 (2010): 3794-3803.

12. K. Ciura, S. Dziomba, J. Nowakowska, and M. J. Markuszewski, “Thin Layer Chromatography in Drug Discovery Process,” Journal of Chromatography A 1520 (2017): 9-22.

13. C. F. Poole, “Thin-Layer Chromatography: Challenges and Opportunities,” Journal of Chromatography A 1000, no. 1 (2003): 963-984.

14. Ł. Komsta, “A Functional-Based Approach to the Retention in Thin Layer Chromatographic Screening Systems,” Analytica Chimica Acta 629, no. 1 (2008): 66-72, https://doi.org/10.1016/j.aca.2008.09.036.

15. S. Yousefinejad, F. Honarasa, and N. Saeed, “Quantitative Structure-Retardation Factor Relationship of Protein Amino Acids in Different Solvent Mixtures for Normal-Phase Thin-Layer Chromatography,” Journal of Separation Science 38, no. 10 (2015): 1771-1776.

16. H. Xu, J. Lin, Q. Liu, et al., “High-Throughput Discovery of Chemical Structure-Polarity Relationships Combining Automation and Machine-Learning Techniques,” Chem 8, no. 12 (2022): 3202-3214, https://doi.org/ 10.1016/j.chempr.2022.08.008.

17. H. Xu, D. Zhang, and F. Mo, “High-Throughput Automated Platform for Thin Layer Chromatography Analysis,” STAR Protocols 3, no. 4 (2022): 101893.

18. A. F. Zahrt, J. J. Henle, B. T. Rose, Y. Wang, W. T. Darrow, and S. E. Denmark, “Prediction of Higher-Selectivity Catalysts by Computer-Driven Workflow and Machine Learning,” Science 363, no. 6424 (2019): eaau5631, https://doi.org/10.1126/science.aau5631.

19. E. N. Muratov, J. Bajorath, R. P. Sheridan, et al., “Qsar Without Borders,” Chemical Society Reviews 49, no. 11 (2020): 3525-3564.

20. C. B. Wahl, M. Aykol, J. H. Swisher, J. H. Montoya, S. K. Suram, and C. A. Mirkin, “Machine Learning-Accelerated Design and Synthesis of Polyelemental Heterostructures,” Science Advances 7, no. 52 (2021): eabj5505, https://doi.org/10.1126/sciadv.abj5505.

21. N. I. Rinehart, R. K. Saunthwal, J. Wellauer, et al., “A Machine-Learning Tool to Predict Substrate-Adaptive Conditions for Pd-Catalyzed C-N Couplings,” Science 381, no. 6661 (2023): 965-972, https://doi.org/10.1126/science.adg2114.

22. J. Lin, H. Zhang, M. Asadi, et al., “Machine Learning-Driven Discovery and Structure-Activity Relationship Analysis of Conductive MetalOrganic Frameworks,” Chemistry of Materials 36, no. 11 (2024): 5436-5445, https://doi.org/10.1021/acs.chemmater.4c00229.

23. Q. Yang, Y. Li, J.-D. Yang, et al., “Holistic Prediction of the Pka in Diverse Solvents Based on a Machine-Learning Approach,” Angewandte Chemie International Edition 59, no. 43 (2020): 19282-19291, https://doi.org/10.1002/anie.202008528.

24. Y. Yao, Q. Dong, A. Brozena, et al., “High-Entropy Nanoparticles: Synthesis-Structure-Property Relationships and Data-Driven Discovery,” Science 376, no. 6589 (2022): eabn3103, https://doi.org/10.1126/science.abn3103.

25. B. A. Koscher, R. B. Canty, M. A. McDonald, et al., “Autonomous, Multiproperty-Driven Molecular Discovery: From Predictions to Measurements and Back,” Science 382, no. 6677 (2023): eadi1407, https://doi.org/10.1126/science.adi1407.

26. C. Rudin, “Stop Explaining Black Box Machine Learning Models for High Stakes Decisions and Use Interpretable Models Instead,” Nature Machine Intelligence 1, no. 5 (2019): 206-215.

27. R. Geirhos, J.-H. Jacobsen, C. Michaelis, et al., “Shortcut Learning in Deep Neural Networks,” Nature Machine Intelligence 2, no. 11 (2020): 665- 673.

28. Z. C. Lipton, “The Mythos of Model Interpretability: In Machine Learning, the Concept of Interpretability Is Both Important and Slippery,” Queue 16, no. 3 (2018): 31-57.

29. L. Billard and E. Diday, “From the Statistics of Data to the Statistics of Knowledge: Symbolic Data Analysis,” Journal of the American Statistical Association 98, no. 462 (2003): 470-487.

30. M. Quade, M. Abel, K. Shafi, R. K. Niven, and B. R. Noack, “Prediction of Dynamical Systems by Symbolic Regression,” Physical Review E 94, no. 1 (2016): 012214.

31. S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering Governing Equations from Data by Sparse Identification of Nonlinear Dynamical Systems,” Proceedings of the National academy of Sciences of the United States of America 113, no. 15 (2016): 3932-3937.

32. B. K. Petersen, M. Landajuela, T. N. Mundhenk, C. P. Santiago, S. Kim, and J. T. Kim, “Deep Symbolic Regression: Recovering Mathematical Expressions from Data via Risk-Seeking Policy Gradients,” in Proceedings of the International Conference on Learning Representations (2021).

33. M. Cranmer, A. Sanchez-Gonzalez, P. Battaglia, et al., “Discovering Symbolic Models from Deep Learning with Inductive Biases,” Proceeding International Conference on Neural Information Processing Systems 33 (2020) 17429-17442.

34. N. J. O’Connor, A. Jonayat, M. J. Janik, and T. P. Senftle, “Interaction Trends between Single Metal Atoms and Oxide Supports Identified with Density Functional Theory and Statistical Learning,” Nature Catalysis 1, no. 7 (2018): 531-539.

35. B. Weng, Z. Song, R. Zhu, et al., “Simple Descriptor Derived from Symbolic Regression Accelerating the Discovery of New Perovskite Catalysts,” Nature communications 11, no. 1 (2020): 3513.

36. S.-M. Udrescu and M. Tegmark, “Ai Feynman: A Physics-Inspired Method for Symbolic Regression,” Science Advances 6, no. 16 (2020): eaay2631.

37. A. A. Kaptanoglu, J. L. Callaham, A. Aravkin, C. J. Hansen, and S. L. Brunton, “Promoting Global Stability in Data-Driven Models of Quadratic Nonlinear Dynamics,” Physical Review Fluids 6, no. 9 (2021): 094401.

38. A. E. Allen and A. Tkatchenko, “Machine Learning of Material Properties: Predictive and Interpretable Multilinear Models,” Science Advances 8, no. 18 (2022): eabm7185.

39. J. Zeng, H. Xu, Y. Chen, and D. Zhang, “Deep Learning Discovery of Macroscopic Governing Equations for Viscous Gravity Currents from Microscopic Simulation Data,” Computers & Geosciences 27 (2023): 987-1000.

40. Z. Li, H. Yuan, W. Han, et al., “Bi-Level Identification of Governing Equations for Nonlinear Physical Systems,” Nature Computational Science 5 (2025): 456-466.

41. J. L. Durant, B. A. Leland, D. R. Henry, and J. G. Nourse, “Reopti-mization of Mdl Keys for Use in Drug Discovery,” Journal of Chemical Information and Modeling 42, no. 6 (2002): 1273-1280.

42. Y. Chong, Y. Huo, S. Jiang, et al., “Machine Learning of Spectra-Property Relationship for Imperfect and Small Chemistry Data,” Proceedings of the National Academy of Sciences of the United States of America 120, no. 20 (2023): e2220789120.

43. I. Smallwood, Handbook of Organic Solvent Properties (Butterworth-Heinemann, Oxford, 2012).

44. Y.-J. Cha, S. Baek, G. Ahn, et al., “Compensating for the Loss of Human Distinctiveness: The Use of Social Creativity under Human-Machine Comparisons,” Computers in Human Behavior 103 (2020): 80-90.

45. J. Sherma and B. Fried, Handbook of Thin-Layer Chromatography (CRC, Eastern, Pennsylvania, 2003).

46. H. Wynberg, “The Reimer-Tiemann Reaction,” Chemical Reviews 60, no. 2 (1960): 169-184.

47. W. Wu, H. Xu, Y. Xu, et al., “Intelligent Column Chromatography Prediction Model Based on Automation and Machine Learning,” Chem 11 (2025): 102598.

48. I. Arnaldo, K. Krawiec, and U.-M. O’Reilly, “Multiple Regression Genetic Programming,” in Proceeding Annual Conference on Genetic and Evolutionary Computation (2014), 879-886, https://doi.org/10.1145/ 2576768.2598291.

49. M. Cranmer, “Interpretable Machine Learning for Science with Pysr and Symbolicregression. Jl,” 2023.

50. M. Du, Y. Chen, and D. Zhang, “Discover: Deep Identification of Symbolically Concise Open-Form Partial Differential Equations via Enhanced Reinforcement Learning,” Physical Review Research 6 (2024): 013182.

51. M. T. Ribeiro, S. Singh, and C. Guestrin, “Why Should I Trust You? Explaining the Predictions of Any Classifier,” in Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (2016): 1135-1144.

52. S. M. Lundberg and S.-I. Lee, “A Unified Approach to Interpreting Model Predictions,” Advances in Neural Information Processing Systems 30 (2017): 4768-4777.

Supporting Information

Additional supporting information can be found online in the Supporting Information section.

Supporting File: advs73626-sup-0001-SuppMat.pdf.