Modeling radical scavenging activity using molecular descriptors of unclassified compounds

2.1 Data collection

The dataset included 3,133 records of in vitro antioxidant measurements, compound names, experimental activities, and simplified molecular-input line-entry system (SMILES) data, sourced from an antioxidant database (https://aodb.idruglab.cn/) (Deng et al., 2023, Weininger, 1988). The radical scavenging activity against DPPH was expressed as IC50 (nM). IC50 values were converted to pIC50 values by applying the -log10 method. The pIC50 values ranged from -0.98 to 10.30. SMILES comparisons were used to identify duplicate structures, and rows with duplicate smiles were removed to keep a single representation to avoid data bias. Oversampling and undersampling were not performed due to the small number of outliers (Supplementary Fig. S1).

The dataset’s compounds were divided into training and test sets, where the training set included 2,506 compounds and the test set contained 627 compounds.

2.2 Calculation of molecular descriptors

Molecular descriptors are simple mathematical expressions used to encode the chemical and structural features of a molecule. Molecular descriptors were calculated using the Python open-source RDKit library (Landrum, 2013). The descriptors were calculated using the RDKit module after the conversion of SMILES to Mol.

During model development, highly collinear molecular descriptors were carefully excluded. Testing and removal of descriptors that encode similar molecular information and those with little variation (descriptors with a Pearson correlation coefficient absolute value greater than 0.9) are important because collinear descriptors encode similar information. Subsequently, high-impact descriptors were selected using an RF model. This method involves plotting the correlation-coefficient values obtained for the different cutoffs of the descriptor significance included in each model to analyze the statistical quality improvement of the model. Hence, 82 descriptors that could be used to develop linear and non-linear models were screened.

2.3 QSAR model development

QSAR was used to predict the antioxidant activity of the compounds. The QSAR model was trained using three algorithms: XGBoost, RF, and SVM.

2.4 K-fold-cross validation

Cross-validation is a model validation technique used to test what kind of results an analysis created by a machine learning algorithm would produce on an independent dataset. For this study, we chose a k value of 5 for k-fold cross-validation, which means we divided the dataset into five subsets with similar proportions of each class represented in each subset. Training and testing were performed in five iterations. In each iteration, one subset was used for testing, and the rest (four subsets) were used for training. Each subset was also used for testing once at the same time (Anguita et al., 2012).

2.5 Model optimization

The hyperparameters of the learning algorithm were optimized using a grid search method, which is a decision-theoretic approach that searches for the optimal configuration by sequentially applying fixed hyperparameter values (Yang and Shami, 2020). Select the hyperparameters to tune and determine the search range. Get all possible combinations and optimize them to have the smallest root mean square error (RMSE) value.

2.6 Model evaluation and validation

Model validation is crucial in the development of QSAR models. Various methods and statistical parameters were used for model evaluation. When there are no restrictions on the applicability domain, each model can theoretically forecast the activity of any compound, including those with structures entirely distinct from those in the training set. Consequently, lacking an applicability domain, a critical element of QSAR models, may cause undue extrapolation across chemical space. This can significantly increase the risk of making inaccurate predictions.

Therefore, it is advised to adhere to a specific set of statistical criteria that the prediction model must meet. For continuous QSAR, the recommended guidelines for creating activity/property predictors include the following. The slope (K, K’) between the predicted and actual observed values should be calculated to evaluate the linear relationship between the two variables. Cross-validation determination coefficient (Q²) should be used to assess the generalization ability of the model, calculated through cross-validation of the training set. The correlation-coefficient (R) measures the strength of the correlation between predicted and observed values. The determination coefficient (R²) indicates how well the model predictions explain the actual observed values. R₀² and R’₀² are calculated through linear regression using the regression of observed values for the predicted values and the regression of predicted values for the observed values, respectively. These were used to assess the extent to which the model predictions matched the actual observed values. The validation criteria outlined by Tropsha for a model to be deemed acceptable include the following specifications:

Models that meet these conditions are considered to have reliable predictive power (Tropsha, 2010).

3.1 Descriptor screening

As mentioned previously, the quantitative relationship between the activity and structural characteristics of antioxidant compounds was investigated using linear and non-linear models. The main step in descriptor screening is the selection of the most influential descriptors. During model development, care was taken to avoid including highly linear molecular descriptors. Descriptors with high Pearson correlation coefficients may encode similar molecular information, and the model may learn redundant information, leading to overfitting. Therefore, it was important to test the descriptors and remove those with low variation and those that encode hereditary information (Fatemi and Gholami Rostami, 2013).

The pIC50 values in our study show a large variation from -0.98 to 10.30, which suggests that while covering a wide spectrum from low to high activity increases the versatility of the model, it may also increase the prediction error somewhat in the extreme value range. Therefore, further experimental validation is suggested as future work to improve the accuracy in this range.

We performed a Pearson correlation analysis of the descriptors to remove those with absolute correlation values greater than 0.9 (Fig. 1). Then, 82 descriptors were selected by using an RF algorithm to evaluate attribute importance. Subsequently, descriptors with attribute importance greater than 0.0025 were screened (Fig. 2). These selected descriptors were used for training the QSAR models. In our research, we employed three ML algorithms, XGBoost, RF, and SVM, to train the QSAR models.

Close

Fig. 1.

Heatmap of Pearson correlation coefficients for the selected 82 descriptors.

Export to PPT

Close

Fig. 2.

The relationship between R² values and the count of descriptors, along with their significance in the random forest model.

Export to PPT

3.2 Hyperparameter optimization and QSAR model SHAP analysis

In the QSAR modelling process, we used a grid search method to select the hyperparameters of the model and found the hyperparameters with the lowest RMSE, i.e. C, the optimal hyperparameters. The lower the RMSE, the better the performance of the model. We applied the 5-fold cross-validation method to optimize the hyperparameters using only the training data set and calculated the cross-validation correlation coefficient (Q²). When comparing the RMSE values of the three ML algorithms, the XGBoost model had the lowest value (Hodson, 2022). The best-performing XGBoost model and its hyperparameter exploration coverage have been shown in Table 1. The optimal hyperparameter values were n_estimators = 400, learning_rate = 0.05, max_depth = 7, and colsample_tree = 0.3. Fig. 3 represents the plot of experimental pIC50 values against the calculated pIC50 values using the XGBoost, SVM, and RF models. A good match between the calculated and experimental values indicates the fitness of the developed models. Fig. 4 displays the residuals of the predicted pIC50 values compared to the experimental data. The spread of residuals around the zero line suggests that the XGBoost model developed does not contain systematic errors. According to the results of the SHAP analysis for the XGBoost model, the impact of the descriptors have been depicted in Fig. 5. Table 2 outlines the names and definitions of these descriptors, all of which represent the topological and electronic dimensions of the chemical structures of compounds that influence antioxidant activity. Among these, the most critical descriptor is fragment-Ar-OH, highlighting its effectiveness over others by indicating the presence of phenolic groups within a compound’s structure. Similarly, the fourth most important descriptor, NHOHCount, represents the number of NH and hydroxyl groups in a compound. As previously discussed, the existence and location of hydroxyl groups on the benzene ring markedly affect the compound’s ability to scavenge radicals (Duchowicz et al., 2019; Shi, 2021). The values of this descriptor positively impact the antioxidant activity. This observation can be explained by the hydrogen atom transfer reaction pathway (Mayer, 2011). The second most important descriptor, MinPartialCharge, is also related to the density of electrons within the compound. A lower minimum partial charge within the compound positively impacts antioxidant activity. This observation can be explained according to the literature (Klein et al., 2005).

Close

Fig. 3.

The plot of XGBoost (a), Random Forest (b), Support vector machine (c) predicted versus experimental of values pIC50.

Export to PPT

Close

Fig. 4.

The plot of XGBoost (a), Random forest (b), and Support vector machine (c) residuals versus experimental of values pIC50.

Export to PPT

Close

Fig. 5.

(a) The impact of the descriptor on the output of the model (SHAP value) and (b) average of the absolute values of the impacts.

Export to PPT

3.3 QSAR model validation and evaluation

We adopted the methodology outlined by Tropsha to evaluate the applicability of the QSAR model. A high value of Q² alone is not a sufficient criterion for a QSAR model to have high predictive power, so researchers have described reliable predictive performance not only in terms of the value of Q² but also in terms of fulfilling the conditions of eq (2) to eq (5) (Golbraikh and Tropsha, 2002; Tropsha et al., 2003). Table 3 outlines the statistical metrics for the external validation set used across the three models. The results from the XGBoost model align closely with the previously defined limits, further underscoring the XGBoost model’s strong predictive capability. A comparison of these statistics with those obtained from the RF and SVM models demonstrates the superiority of the DT-based models. Among the DT-based models, those applying the boosting algorithm proved superior in predicting the pIC50 values. We also compare the QSAR model in our study with the reported QSAR models. Table 4 shows that the reported QSAR models used classified compounds as training data, while our QSAR model used unclassified compounds as training data. Although our model did not outperform Shi’s SVM or Duchowicz et al. (2019) MLR model in terms of R2 (0.85-0.93) for specific flavonoid or phenolic acid datasets (23-75 compounds) in terms of accuracy, However, our model includes 3,133 unclassified compounds, which ensures broader structural coverage. Although the R² is 0.81, greater chemical diversity increases real-world applicability, demonstrating the robustness of the model across different compound classes. As a result, the XGBoost model provides a useful alternative for those looking for a ML solution to predict the antioxidant activity of unclassified compounds.