A new simple and efficient molecular descriptor for the fast and accurate prediction of log P

Dataset

In 2018, Mansouri et al. developed the OPEn structure-activity/property relationship app (OPERA) models based on their dataset^[16]. We used their data as the initial training and testing sets, referring to it as M-dataset. The M-dataset was derived from the available PHYSPOROP physicochemical property and environmental fate datasets, with extensive curation conducted to ensure data quality. It provided log P values and chemical structure formats for 14,050 molecules, along with identifiers such as simplified molecular input line entry specification (SMILES) and structure data file (SDF). In this work, SMILES of 13,963 molecules and their experimental log P values were selected to create the descriptors. Finally, descriptors of 13,952 molecules were calculated using homemade code based on the RDKit library.

The datasets for prediction were obtained from the SAMPL6 challenge and the SAMPL9 challenge. The SAMPL6 challenge provided a set of 11 drug-like molecules in SMILES string format and their experimental values, and the SAMPL9 challenge provided 16 molecules [Figure 1], which had experimental log P values ranging from -1.37 log P units to 4.92 log P units. With the SMILES provided, the descriptors used for prediction were created in the same way as the ones used for training and testing.

Figure 1. (A) Eleven molecules in the SAMPL6 log P challenge; (B) Sixteen molecules in the SAMPL9 log P challenge.

Descriptors

The concept of 3D molecular representation of structures based on electron diffraction (3D-MoRSE) descriptors was initially introduced by Schuur, Selzer, and Gasteiger in 1996^[17,18]. These descriptors are expressed as $$ I(s)=\sum_{i=2}^{N}\sum_{j=1}^{i-1}A_iA_jf(s,r_{ij}) $$, where $$ f(s,r_{ij})=\frac{\mathrm{sin}sr_{ij}}{sr_{ij}} $$ serves as the core function, where s is the scattering parameter ranging from 0 to 31 Å^-1, r_ij is the distance between i th and j th atoms, N is the total number of atoms, and A_i and A_j are the atomic weights, which can be unweighted or represent various atomic properties. The atomic weights used in this work are shown in Table 1.

Our previous works have shown that the ML approach with intermolecular 3D-MoRSE descriptors could realize very high electronic coupling prediction between molecules^[19] and the opt3DM descriptors could act as efficient descriptors for developing interface modifiers in perovskite solar cells^[20]. Based on these works, we here optimized 3D-MoRSE descriptors for the log P prediction. A s_L was added to adjust f(s, r_ij), which is defined as $$ f(s,r_{ij})=\frac{\mathrm{sin}s\times s_L\times r_{ij}}{s\times s_L\times r_{ij}} $$. By meticulously adjusting the coefficients (s × s_L) and the dimension (N_s, range of s), this updated descriptor set has demonstrated the potential to surpass all currently available descriptors in terms of prediction accuracy. As found, the best prediction results were calculated when s_L= 0.5 and N_s = 500.

ML model

The ML algorithms used in this work were implemented with the scikit-learn library. Our ML models consisted of a feature selector and a fitting regressor. The SelectFromModel, a feature selector from the scikit-learn library, served as the selector in the ML model. Various algorithms, including automatic relevance determination (ARD) regression, Bayesian Ridge, and Ridge, were employed in both the feature selection and regression stages. These algorithms were trained on the M-dataset to identify the most suitable ones for the ML model, as illustrated in the subsequent sections. The selected ML model was used for the prediction of the SAMPL6 and the SAMPL9 datasets. The mean absolute error (MAE), RMSE, and coefficients of determination (R²) serve as evaluation metrics, which are defined as

Where n represents the number of molecules, and y_i and $$ \hat{y}_i $$ are the experimental values and the predicted values of log P, respectively.

Prediction on the M-dataset

The M-dataset, created by Mansouri et al., is a dataset providing information for predictions on physicochemical properties and other fields and it is available on GitHub^[21]. The log P distribution is shown in Supplementary Figure 1. Mansouri et al. utilized genetic algorithms (GA) to select features and employed a k-nearest neighbors (KNN) approach for model fitting and prediction^[16]. Their GA-KNN model achieved an RMSE of 0.78 and an R² of 0.86 in the test set.

To achieve the best prediction on the OPERA dataset, selecting descriptors and algorithms is crucial. By introducing the s_L, the range and sensitivity of the coefficients can be adjusted. The s_L serves as a parameter that adjusts the granularity at which the descriptors operate. A lower s_L increases the sensitivity of the descriptors, allowing them to detect finer variations in the molecular stackings, while a higher s_L might simplify the descriptor, potentially overlooking critical details. Following our previous work^[20], we have implemented algorithms referred to in the Methods section and tested the s_L from 0.03 to 3.0 and the N_s up to 500 on the M-dataset with 90% for training and 10% for testing [Supplementary Tables 1 and 2]. Figure 2A shows the s_L-dependent RMSEs of various algorithms with descriptors calculated with 100 dimensions (N_s = 100). The Ridge, Bayesian Ridge, and ARD regression achieved lower RMSEs among the algorithms, in which we found that the optimal s_L is between 0.1-1.0. As shown in Figure 2B, the result showcased that feature selection did not affect the selection of s_L. Second, we increased the dimension to 500 and screened s_L in the optimal range with 10 tests [Table 2]. The best prediction result (with the smallest average RMSE for independent ten predictions) was achieved when s_L = 0.50 and N_s = 500.

Figure 2. The RMSEs of predictions for the test set of M-dataset with various algorithms and optimized 3D-MoRSE descriptors calculated by various s_L. (A) The s_L-dependent RMSE without descriptor selection; (B) The s_L-dependent RMSE with descriptor selection. The testing ratio was 0.1. The descriptors were selected using ARD regression. RMSEs: Root mean square errors; 3D-MoRSE: 3D molecular representation of structures based on electron diffraction; ARD: automatic relevance determination.

To abstract as much information as possible from molecules, a large N_s is required. However, the large N_s increases the dimension, and also includes some reductant information. Therefore, a feature selection was necessary. We employed the SelectFromModel module (Meta-transformer for selecting features based on importance weights, which is a module in the scikit-learn library) with ARD regression as the estimator for feature selection. In this process, we used 90% of the dataset to select features and to train the ML model, as 10% of the dataset was set for testing. The dataset was split randomly with a random state ranging from 0 to 20 in order to obtain robust and reliable results. Moreover, a threshold was set for finding the optimal dimension and we found the optimal threshold was 0. 1 [Supplementary Tables 3-8]. However, noticing that the reduction of feature dimension may enhance the generalized capability, features selected at a threshold of 0.05 will be considered. To test the average performance of these descriptors, we employed these descriptors to fit the model and predict the test set based on those divided datasets. The average performance of descriptors selected in various thresholds and the dimensions of these descriptors are shown in Figure 3.

Figure 3. The RMSEs and feature dimensions of models using descriptors selected with various thresholds. (A) The RMSEs of models fitted with various algorithms. The optimal threshold was around 0.001 and 0.05; (B) The feature dimension decreased rapidly with the increased threshold. RMSEs: Root mean square errors.

To have a direct comparison with the work of Mansouri et al., we used the same dataset and training and testing splitting ratio of 0.75:0.25^[16]. To ensure robustness, the dataset was partitioned into training and test subsets 20 times with different random states, resulting in 20 distinct tests [Supplementary Table 9]. This iterative process mitigates the risk of coincidental outcomes and substantiates the model’s dependability. The best and the worst predictions have not much difference in prediction accuracy, as shown in Figure 4, demonstrating that our models are very stable and robust. Table 3 lists the performance of the OPERA model and our models. The OPERA model predicts on test set for one time and achieves an RMSE of 0.78, while our models demonstrate an average RMSE of 0.68. Collectively, these results indicate that our model demonstrated superior performance compared to the OPERA model.

Figure 4. The experimental log P values and the ML predicted log P values. (A) The ARD regression model with the best (A) and the worst (B) performance in 20 simulations with different random states; The Bayesian Ridge regression model with the best (C) and the worst (D) performance in 20 simulations with different random states. The Ridge regression model with the best (E) and the worst (F) performance in 20 simulations with different random states. All models used the coefficient set s_L = 0.5 and N_s = 500. The test set was 0.25. ML: machine learning; ARD: automatic relevance determination.

Predictions on the SAMPL6 dataset and the SAMPL9 dataset

The SAMPL6 and the SAMPL9 challenges are the two most famous log P challenges, having received over 100 submissions from various groups. To verify the reliability of our models, we utilized the selected descriptors and models trained on the M-dataset to predict the log P in the SAMPL6 and SAMPL9 challenges.

We employed the ARD regression, Bayesian Ridge Regression and Ridge Regression to train our models, which achieved average RMSEs of 0.309, 0.326 and 0.329 [Supplementary Tables 10 and 11], respectively. Notably, the minimum values achieved was an RMSE of 0.29 [Figure 5A]. The comparative performance of our models and other state-of-the-art models on the SAMPL6 dataset is detailed in Supplementary Table 10 and Figure 6A. The other models encompassed various methodologies, including deep learning, QC calculations, MD simulations, and other complex approaches. Our models outperformed most of the challengers. Among them, only the deep neural network with tautomer (DNNtaut) and the deep neural network with original SMILES (DNNmono) demonstrated similar performance compared to ours. The DNNmono, a deep learning model developed by Ulrich et al., reported the lowest RMSE of 0.31 and typical RMSE of 0.33 in their other models^[12]. Our models show the lowest RMSE on the SAMPL6 challenge among all models, outperform QC- and MD-based models, and also surpass deep neural network models. They offer advantages in terms of complexity, indicating faster prediction times. Additionally, these models exhibited commendable stability, which is reflected in the compact distribution of MAE, RMSE, and R² across 20 stochastic parameter configurations. Overall, our models provided the most accurate prediction on the SAMPL6 dataset among more than 100 submissions, emphasizing the superior efficiency of our developed opt3DM descriptors.

Figure 5. The experimental and predictive log P using the ARD regression model with the best random state. (A) This image showcases the prediction performance on the SAMPL6 dataset. The error range was 0.30; (B) This image showcases the prediction performance on the SAMPL9 dataset. The error range was 1.00. The MAE, RMSE, and R² are provided. The model was trained on the M-dataset using the coefficient set s_L = 0.5 and N_s = 500. ARD: Automatic relevance determination; MAE: mean absolute error; RMSE: root mean square error; R²: coefficients of determination.

Figure 6. The RMSEs of various models and our models. (A) The RMSEs of models in SAMPL6 challenges and our models; (B) The lowest RMSEs in one literation of models in SAMPL9 challenges and our models. RMSEs: Root mean square errors.

Our models not only demonstrated robust performance in the prediction of the M-dataset and the SAMPL6 dataset, but also showed good prediction in the SAMPL9 challenge. From 2021 to 2023, submissions to the SAMPL9 challenge showed varied results, with RMSEs ranging from 1.52 to 2.95^[15,22] [Supplementary Table 12]. Later, Nevolianis et al. replicated the success of the COSMO-RS approach in the SAMPL9 challenge achieving an RMSE of 1.23^[23]. Despite their sophistication, these physics-based approaches demand substantial computational resources as their model requires QC calculations. Besides, the submission by Nevolianis et al. to SAMPL9, using the COSMO-RS approach, required 3.7 h to predict the log P of just sixteen molecules, indicating that predictions for a larger number of molecules would be time-consuming^[23]. By using ARD regression, we achieved the lowest average RMSE of 1.01 and attained a minimum RMSE of 0.97 in a single iteration [Figure 5B and Supplementary Table 13]. The best prediction on SAMPL9 was reported by Zamora et al.^[13]. However, duplicate entries with SAMPL9 data were found in their training dataset^[15]. Therefore, we can conclude that our model again gives the most accurate prediction in the SAMPL9 challenge with blind test, as shown in Figure 6B. Additionally, we trained the ML models on the M-dataset consisting of the n-octanol/water log P (log P_o/w), and realized accurate prediction on the SAMLP9 data that are all toluene/water log P (log P_tol/w), demonstrating the generalization ability of our model in predicting log P.

To identify the molecular features with great importance to the log P, we performed Shapley Additive Explanations (SHAP) analysis^[24], and the results are shown in Figure 7A and B. In previous studies, the first half of the sine wave period is pronounced. When s_L = 0.5, the half-sine wave period is about 3.14 Å for MorU-2, MorRC-2, and MorIP-2. They are longer than most of the chemical bonds (about 1.5 Å) in molecules. As the atomic weight is 1.0 for MorU, the descriptor value is solely determined by bond lengths. As the bond length decreases, the MorU-2 value increases. Unsaturated carbon chains, compared to alkanes, generally have stronger hydrophilicity. This is because the double or triple bonds present in unsaturated carbon chains introduce polarity, reducing the overall non-polar nature and enhancing the polarity of the molecule. In unsaturated carbon chains, the presence of double bonds leads to an uneven distribution of electron cloud density, which imparts polarity to the molecule, thereby increasing the interaction with water molecules and making these molecules more hydrophilic. In contrast, alkanes, which only have single bonds, have a more uniform electron cloud distribution and are non-polar overall, thus exhibiting stronger hydrophobicity. The more unsaturated carbon chains a molecule has, the denser the molecule is, the larger the MorU-2 value is, and the more hydrophilic the molecule becomes. It is the same for MorRC-2, which gives carbon a larger weight. Therefore, MorU-2 and MorRC-2 are negatively correlated with log P. The MorIP-2 gives hydrogen a larger weight. A greater MorIP-2 value indicates more alkanes, leading to hydrophobicity. Therefore, MorIP-2 is positively correlated with log P. For the MorE-9, the sine wave period is about 1.40 Å. As shown in Figure 7C, its value increases with the length of bonds in most cases. More alkanes cause a smaller MorE-9 value and lead to hydrophobicity. Therefore, MorIP-2 is negatively correlated with log P. For MorRC-26, the sine wave period is about 0.48 Å. As shown in Figure 7C, the C–F and C–Cl bonds in halogenated hydrocarbons which are hydrophobic give a negative value on the MorRC-26 curve; The C–O and C–N bonds in hydrophilic groups have little impact on the MorRC-26 value. Therefore, MorRC-26 is negatively correlated with log P.

Figure 7. The SHAP values of twenty ARD regression models on (A) M-dataset, (B) SAMPL 6 and SAMPL 9 and (C) The f(s, r_ij) of MorX-2 (X = U, RC, IP), MorRC-26 and MorE-9 and the function values of atoms with various distances in different chemical bonds. SHAP: Shapley Additive Explanations; ARD: automatic relevance determination.

Figure 8 illustrates the RMSE of different ML models applied to the SAMPL6 dataset, along with their complexity and computational expense. Notably, our model, being the simplest, achieved the lowest RMSE. Its minimal complexity also resulted in significantly reduced computational costs compared to other models, highlighting the efficiency of our opt3DM descriptors. These descriptors have significantly advanced rapid and precise log P prediction, facilitating extensive drug and material screening.

Figure 8. The RMSE of predictions, model complexity and categories of various approaches on SAMPL6 dataset. RMSE: Root mean square errors.