ADMET evaluation in drug discovery: 21. Application and industrial validation of machine learning algorithms for Caco-2 permeability prediction

Abstract

The Caco-2 cell model has been widely used to assess the intestinal permeability of drug candidates in vitro, owing to its morphological and functional similarity to human enterocytes. While Caco-2 cell assay is considered safe and cost-effective, it is also characterized by being time-consuming. Therefore, computational models that achieve high accuracies in predicting Caco-2 permeability are crucial for enhancing the efficiency of oral drug development. In this study, we conducted an in-depth analysis of the characteristics of an augmented Caco-2 permeability dataset, and evaluated a diverse range of machine learning algorithms in combination with different molecular representations. The results indicated that XGBoost generally provided better predictions than comparable models for the test sets. In addition, we investigated the transferability of machine learning models trained on publicly available data to internal pharmaceutical industry datasets. Our findings, based on the Shanghai Qilu’s in-house dataset, showed that the boosting models retained a degree of predictive efficacy when applied to industry data. Furthermore, Y-randomization test and applicability domain analysis were employed to assess the robustness and generalizability of these models. Matched Molecular Pair Analysis (MMPA) was utilized to extract chemical transformation rules. We believe that the model developed in this study could represent a reliable tool for assessing Caco-2 permeability during early-stage drug discovery and the chemical transformation rules derived here could provide insights for optimizing Caco-2 permeability.

Scientific contribution

A comprehensive validation of various machine learning algorithms combined with diverse molecular representations on a large dataset for predicting Caco-2 permeability was reported. The transferability of machine learning models trained on publicly available data to internal pharmaceutical industry datasets was also investigated. Matched molecular pair analysis was carried out to provide reasonable suggestions for researchers to improve the Caco-2 permeability of compounds.

Graphical Abstract

Introduction

The oral route is often preferred for drug administration because of its convenience, cost-effectiveness, and high patient adherence [1]. However, about 10% of drugs fail in development due to poor pharmacokinetic properties [2]. One of the foremost challenges faced by orally administered drugs is the traversal across the intestinal epithelial barrier, which determines the rate and extent of absorption in humans, thereby critically influencing their bioavailability [3]. Furthermore, the permeability process itself is intricate. Transport across the intestinal epithelial cell barrier may occur through one or more of four different routes: paracellular and transcellular passive diffusion, the carrier-mediated route, and transcytosis [4]. Therefore, it is paramount to accurately reconstitute the differentiated human epithelial cell monolayer in vitro, as this enables precise prediction of oral drug absorption in humans.

Several cells culture models have been developed to replicate the pertinent characteristics of in vivo absorption, such as Parallel Artificial Membrane Permeability Assay (PAMPA) [5], the human colon adenocarcinoma (Caco-2) cell lines [6], the Madin–Darby canine kidney (MDCK) cell [7], and porcine kidney epithelial cell lines (LLC-PK1) [8]. Among these models, Caco-2 cell monolayers have emerged as the “gold standard” for drug permeability due to their ability to closely mimic the human intestinal epithelium [9]. Indeed, this in vitro model has been endorsed by the US Food and Drug Administration (FDA) for assessing the permeability of compounds categorized under the Biopharmaceutics Classification System (BCS) [10].

Nevertheless, high-throughput screening (HTS) with the traditional Caco-2 cell model poses challenges due to its extended culturing period (7–21 days), which is necessary for full differentiation into an enterocyte-like phenotype [11, 12]. The extended cultivation not only increases the risk of contamination but also imposes significant costs, making it less desirable in the context of drug discovery. Furthermore, there is a growing need to apply selection filters prior to designing and acquiring large compound libraries, thereby shifting the focus towards in silico predictions of ADME-related properties.

Numerous machine learning models have been developed for both classification and regression analysis of Caco-2 permeability based on a variety of molecular representations, such as fingerprints, descriptors and embeddings extracted from neural network models [13,14,15,16,17,18,19]. Notably, there are fewer classification models than regression models for predicting Caco-2 cell permeability, primarily due to the lack of a consensus on the optimal cut-off for modeling and practical application. Here, this study also focused on the regression task. Prior research on Caco-2 permeability prediction has exhibited considerable variation in dataset size, descriptor types, and modeling methods. For instance, Wang et al. constructed different models based on the MOE 2D/3D descriptors to predict P_app values using a structurally diverse data set of 1272 compounds [13]. The comparison between Boosting and other methods revealed that the Boosting model achieved better results, with an R² value of 0.81 and an RMSE of 0.31 for the test set. Additionally, Wang and colleagues proposed the MESN deep learning model, which uses a curated Caco-2 dataset with 4464 compounds to predict the features associated to oral bioavailability [15]. Three molecular embedding approaches were employed in the MESN model, namely, Morgan fingerprints, SMILES-based embedding, and molecular-graph-based embedding. For membrane permeability prediction, the MESN model obtained an MAE score of 0.410 and an RMSE of 0.545. More recently, Gabriela et al. introduced a QSPR approach developed on the KNIME analytical platform, based on a structurally diverse dataset of over 4900 molecules [16]. They employed random forest supervised recursive algorithms for data cleaning and feature selection. The modeling algorithm consisted of a conditional consensus model composed of individual regression random forests, achieving RMSE values ranging from 0.43 to 0.51 and R² values between 0.57–0.61 for all validation sets.

Although progress has been made in predicting Caco-2 permeability, several issues still required discussion. Firstly, Caco-2 permeability is a highly complex process that can occur through various nonlinear routes. The characteristics of the Caco-2 permeability data were not discussed and the scarcity of high-quality Caco-2 permeability data has impeded the development of accurate models with a wide applicability domain. Secondly, there remains a question regarding the performance of deep learning in Caco-2 permeability data. Thirdly, the effectiveness of machine learning models trained on publicly accessible databases in industrial settings remains unclear.

Taking into account the mentioned concerns, we thereby reported a comprehensive validation of various machine learning algorithms combined with diverse molecular representations for predicting Caco-2 permeability. Initially, a relatively larger Caco-2 permeability dataset was collected, and the structural diversity of the compounds was analyzed. Four machine learning methods (XGBoost, RF, GBM, and SVM) and two deep learning models (DMPNN and CombinedNet) were used to build prediction models. For each algorithm, the optimal model was selected for further analysis. Next, we conducted a direct comparison of different in silico predictors through several model validation methods, including the Y-randomization test and application domain (AD) analysis. Additionally, we carried out a performance assessment of different models trained on public data using the Shanghai Qilu’s in-house dataset. Lastly, based on the results of matched molecular pair analysis (MMPA), we provided reasonable suggestions for researchers to improve the Caco-2 permeability of compounds. We anticipate that our findings will shed light on future method development and offer insights for optimizing Caco-2 permeability.

Methods

Data collection and preparation

Experimental values of Caco-2 permeability were obtained from three publicly available datasets [13,14,15]. The first dataset, reported by Wang et al., comprised 1272 compounds used for developing predictive models and an additional 298 compounds for validation [13]. Wang and Cheng collected the second dataset of 1827 compounds to develop QSPR models using neural networks and other machine learning methods [14]. The third dataset, consisting of 4464 compounds, was utilized by Wang et al. for modeling Caco-2 permeability with neural networks [15]. These datasets were combined into an initial dataset of 7861 compounds for subsequent analysis. To ensure data consistency and minimize uncertainty, the following procedures were applied: (1) Permeability measurements were converted to cm/s × 10^–6 and transformed logarithmically (base 10) for modeling. (2) Entries with missing permeability values were excluded. (3) Mean values and standard deviations were calculated for duplicate entries. Only entries with a standard deviation ≤ 0.3 were retained, and the mean values of these compounds were used as standard values for model training. (4) The RDKit module MolStandardize was employed for molecular standardization to achieve consistent tautomer canonical states and final neutral forms, while preserving stereochemistry. After curation, an exhaustive dataset comprising 5654 non-redundant Caco-2 permeability records was compiled. Subsequently, these records were randomly divided into the training, validation, and test sets in an 8:1:1 ratio, ensuring an identical distribution across the datasets. To enhance the robustness of model evaluation against data partitioning variability, the experimental dataset underwent 10 splits using different random seeds. The model was then assessed based on the average performance across these ten independent runs. According to the OECD principles, both internal and external validations are essential to assess model reliability and predictive capability [20, 21]. An additional set of 67 compounds from the Shanghai Qilu’s in-house collection was included as an external validation set to test the prediction performance of the model trained on the public data set on the private data set.

Molecular representations

To incorporate comprehensive chemical information at both global and local levels, three types of molecular representation methods were employed to depict the structural features of molecules. Specifically, the molecular representations utilized were as follows:

• 1. Morgan fingerprints with a radius of 2 and 1024 bits: The RDKit implementation of Morgan fingerprints was utilized [22].

• 2. RDKit2D descriptors: Normalized descriptors from descriptastorus (https://github.com/bp-kelley/descriptastorus), which wraps the RDKit implementation and normalizes values using a cumulative density function from Novartis’ compound catalog.

• 3. Molecular graphs: For the message-passing neural network, molecular graphs $G=(V,E)$ served as the foundational representation, where $V$ represents atoms (nodes) and $E$ represents bonds (edges). This approach was implemented using the open-source ChemProp package to enhance model performance and effectively capture nuanced molecular features [23].

Furthermore, a combination of Morgan fingerprints and RDKit2D normalized descriptors was employed to train all methods, with the exception of DMPNN and CombinedNet, which utilized molecular graphs. Specifically, CombinedNet employed a hybrid approach that combined Morgan fingerprints and molecular graphs: the former providing information on substructure existence and the latter conveying connectivity knowledge.

Model construction and evaluation

In this study, we conducted a comprehensive comparison of multiple machine learning and deep learning algorithms for quantitative predictions of Caco-2 permeability. The evaluated algorithms included Random Forest (RF) [24], extreme gradient boosting (XGBoost) [25], Support Vector Machine (SVM) [26], Gradient Boosting Machine (GBM) [27], Directed Message Passing Neural Network (DMPNN) [28], and CombinedNet (Figure S1) [29]. The first four methods were implemented using the scikit-learn package, while the latter two deep learning techniques were implemented using the ChemProp package [23, 30]. Hyperparameters for all algorithms, except DMPNN and CombinedNet, were optimized through fivefold cross-validation and grid search. The DMPNN algorithm employed the Adam optimizer for gradient descent optimization, while CombinedNet utilized the LAMB optimizer with a mean squared error (MSE) loss function [31]. Hyperparameters in neural networks were fine-tuned using Bayesian optimization algorithms, and to prevent overfitting, an early stopping strategy was applied. Further details about each algorithm can be found in the Table S1.

To ensure robust generalization ability of the QSPR models in predicting Caco-2 permeability for new chemical entities, both fivefold cross-validation and independent test sets were employed. The performance of different machine learning models on the Caco-2 permeability datasets was evaluated utilizing the following metrics: Coefficient of determination (R²), Mean absolute error (MAE), Root mean squared error (RMSE), Pearson correlation coefficient (Pearson’s r), and Spearman’s rank correlation coefficient (Spearman’s rho).

Definition of applicability domain

Due to the limitations inherent by the chemical space represented by the training set, each machine learning model demonstrated a predisposition towards predicting certain types of compounds [32]. This suggests that predictions may be more reliable when the compounds being evaluated fall within the specific applicability domain of each machine learning model. In this study, we adopted a methodology named the Euclidean distance-based method for applicability domain analysis [33]. This method introduces a distance cutoff value ${D}_{c}$ that defines a similarity threshold for external compounds. Morgan fingerprints was utilized to calculate the similarity between molecules. The detailed formula is as follows:

$$D_{c} = \overline{y} + Z\sigma$$

where $\overline{y }$ represents the average Euclidean distance of the k nearest neighbors of each compound in the training set, $\sigma$ is the corresponding standard deviation of Euclidean distances, and $Z$ is an optional parameter to control the significance level. If the distance from an external compound to any of its nearest neighbors in the training set is above the distance cutoff ${D}_{c}$, the compound is considered to be outside of domain; otherwise, it is considered to be within the domain.

Matched molecular pairs (MMP) analysis

Matched molecular pair (MMP) analysis elucidates the impact of specific chemical transformations on activity or property changes within a molecule [34]. This approach has become widely adopted for optimizing both biological activity and ADMET properties [35, 36]. To gain deeper insights into how minor chemical modifications influence Caco-2 permeability properties and to guide ligand optimization, we constructed an MMP database. This database incorporates transformation rules linked to shifts in logP_app derived from our curated Caco-2 permeability datasets. For the creation of this MMP database, we utilized the open-source mmpdb package [37]. The process involved two main stages: fragmentation and indexing. During the fragmentation phase, we permitted up to three cuts per molecule, with constraints on the maximum number of heavy atoms (100) and rotatable bonds (10). Chirality was preserved when cutting bonds. In the subsequent Indexing phase, variable fragments were defined with a minimum of 2 and a maximum of 15 heavy atoms. The maximum radius for indexing was set to 3. The portion of heavy atoms in changing fragments was no more than half of the molecule. Default values were applied for other parameters as per the mmpdb package settings. We included only the transforms that contained at least ten matched pairs and showed a significant property shift with a p-value < 0.05. This rigorous approach ensures that the results comprise meaningful and actionable insights into how specific chemical alterations affect Caco-2 permeability.

Results and discussion

Data set analysis

We collected data from public databases and previous studies, resulting in 7861 compounds being included in this research. After data preparation, a total of 5654 compounds with structural diversity were obtained. The distribution of the logP_app values was shown in Fig. 1a. The logP_app values ranged from − 8.778 to − 3.510, with a mean of − 5.340. The permeabilities were then categorized as follows: compounds with a permeability less than 1 × 10^–6 cm/s are considered low permeability, those up to 7 × 10^–6 cm/s are classified as medium permeability, and compounds with a permeability greater than 7 × 10^–6 cm/s are regarded as highly permeable. It can be observed that the curated dataset is highly unbalanced and contains a higher number of highly permeable compounds (n = 2746) compared to medium (n = 1747) and low (n = 1161) permeable compounds. Figure 1b shows that most structural clusters are quite mixed in terms of permeability classes. However, these clusters are not equally representative of all permeability classes, as many are predominantly composed of highly permeable (green) compounds. The presence of a biased chemical space presents a substantial complication in the accurate prediction of Caco-2 permeability. To further explore the chemical space of the Caco-2 data set, the Tanimoto similarity based on the ECFP4 fingerprints was calculated and the Murcko scaffolds were analyzed. The Murcko scaffolds of the total data set were extracted by removing side chain substituents but retaining the linkers and ring systems with the RDKit package. The overall color of the Tanimoto similarity heat map was light blue with an average similarity of 0.5 (Fig. 1c), indicating the structural diversity of the data set. Additionally, we detected 2979 different Murcko scaffolds from the data set, suggesting that each Murcko scaffold shared an average of ~ 1.9 molecules. The most common scaffolds in the data set were shown in Fig. 1d along with their frequency. Apparently, the six-membered ring was the most frequent scaffold in many drug compounds and molecules with polycyclic scaffolds may be the focus of Caco-2 permeability studies. In a word, the above analysis demonstrated the structural diversity and biased distribution of the current Caco-2 permeability data set.

Fig. 1

Chemical space and scaffold analysis of the curated Caco-2 permeability data set. A The experimental logP_app value distribution showing that most of the compounds belong to the highly permeable class (n = 2746) while the least number of structures are categorized as low permeable. B t-SNE plot of curated Caco-2 permeability data set (green: high permeable compounds; blue: medium permeable; orange: low permeable). C Heat map of Tanimoto similarity based on with the ECFP4 fingerprints of the total data set. D Frequency of the Murcko scaffolds in the data set

Performance of Caco-2 permeability prediction models

First, we evaluated the ability of machine learning approaches to predict Caco-2 permeability. To improve the representation of the chemical and structural features of molecules for classical machine learning and DNN models, various combinations of molecular descriptors and fingerprints were considered. The detailed statistical results for all models across the training, validation, and test sets were summarized in Table S2. The performance was quantified using the R² and RMSE metrics on the test set molecules. Notably, large differences in predictive performance were observed among the data sets, with the R² values ranging from ~ 0.4 to ~ 0.65 log units (Fig. 2a). This trend was also evident in the overall performance of the test set molecules, with the RMSE values ranging from ~ 0.5 to ~ 0.6 log units (Fig. 2b). The performance variations were primarily attributed to the chosen molecular descriptor, with rdkit2d descriptors yielding the relatively low average prediction error. Binary fingerprints generally performed worse than nonbinary descriptors, with higher variation across data sets. Additionally, it was observed that a combined molecular representation would benefit machine learning algorithms.

Fig. 2

Performance of different machine learning methods. The summary of A R² and B RMSE for different machine learning models on test set (indicated by colors). C Global ranking of all methods using PCA, scaled between best and worst performance. Every point represents a different combination of the machine learning method and the descriptor it relied on. “Worst” and “Best” indicated the worst and best performance obtained across all data sets, respectively

To provide a comprehensive assessment of methodological performance across the analyzed data sets, we conducted principal component analysis (PCA) on the obtained RMSE values (Fig. 2c). To enhance the PCA results, we added additional rows representing the best and worst RMSE values for each dataset along the direction defined by both the best and worst outcomes, as observed in previous studies [38, 39]. The proximity of a method’s position to the “best” point in the PCA plot signified its superior overall performance across the considered data sets. Furthermore, the magnitude of orthogonal deviation from the best–worst line indicated the extent of variability in a method’s performance based on the data set. Compared to machine learning models, deep learning models were obviously more sensitive to the specific data used during training. Additionally, the results confirmed the higher impact of molecular descriptors than the chosen machine learning algorithm on the model performance. On average, XGB coupled with a combined molecular representation proved to be the best method for predicting Caco-2 permeability. Machine learning models coupled with ECFPs outperformed DMPNN based on graphs, which was somewhat surprising given that ECFPs were derived from molecular graphs. The finding highlights a current gap in efficiently learning features from “raw” molecular representations in the small-data regimes typical of drug discovery. However, the deep learning model coupled with ECFPs showed comparable results to machine learning models coupled with ECFPs. These results once again emphasize the importance of including both generic molecular representations and local chemical information in MPNN models, which is in alignment with the prior research findings [40].

Understanding of the generalization of Caco-2 permeability prediction models

The model with the most superior performance throughout 10 independent runs for all machine learning algorithms was selected to carry out further analysis. As shown in Fig. 3, a similar tendency was discerned, wherein the majority of combinations of models and molecular representations attained favorable performance in predicting highly permeable compounds. However, medium and, in particular, low permeable compounds exhibited a propensity to be overestimated. We supposed that the suboptimal prediction performance regarding low-permeable compounds is not attributable to the improper choice of algorithms and descriptors. Instead, the model’s bias towards highly-permeable compounds might give rise to an overestimation of their permeability (Fig. 1a and b).

Fig. 3

Comparison of Caco-2 permeabilities predictions to experimentally measured values for test set (blue dash line: two-fold error; orange part of each subplot: low permeable; blue part: medium permeable; green part: high permeable)

To ensure the robustness of our QSPR model against random generation, we employed the Y-randomization test to validate its reliability. This method involved random shuffling of the logP_app values to disrupt their original order [41]. It was expected that the resulting QSPR models would exhibit poor performance when applied to the validation dataset. This randomization process was repeated 100 times, and the collective results were illustrated in Figure S2. It was obvious that the performance of the randomly shuffled models was quite different from the real model, with R² values of approximately − 0.5 to 0. The limited predictive ability of the randomized models demonstrated that our predictive models have adeptly captured the genuine relationship between molecular properties and logP_app, rather than spurious correlations.

Defining the domain of applicability for any predictive model is critical in order to gain an understanding of the model generalization and the origin of estimated errors [32]. Herein, we calculated the similarity between a test compound and its nearest neighbor in the training set to quantify its distance to the model’s domain of applicability and to understand further the potential link with the uncertainty in prediction. More specifically, the structural similarity between any two samples was measured using the Tanimoto coefficient based on the ECFP4 fingerprints with a folding size of 1024 bits. The resulting similarity score was binned in intervals of 0.1 units. Figure 4 highlighted the correlation between the prediction error and the binned average structural similarity for all compounds in the Caco-2 data set. A general trend was observed where the model performance, regardless of machine learning algorithms, increased as the similarity for the test compound against the training set increased, which was in agreement with the conclusions from previous studies [40]. In this case, we believed that compounds with similarities higher than 0.70 to the training set tended to obtain more reliable prediction results.

Fig. 4

Bar chart plots showing the R², RMSE and MAE in prediction as a function of binned structural similarities for different models

The Euclidean distance-based method was then used to evaluate the domain of application. According to equation, it can be presumed that the elevation of the Z value and reduction in k led to a corresponding augmentation of the distance threshold and a consistent diminution of compounds external to the AD. We calculated the average Euclidean distance of the compounds in the training set of 9.203 and the standard deviation of the Euclidean distance of 0.920. For each model, the statistical results before and after considering the applicability domain were collected in Table 1 with the R² values for model fit and RMSE for the test set. As anticipated, the improvement in model performance correlated positively with the reduction in the distance threshold. It makes sense because decreasing the distance threshold causes a greater number of compounds with lower similarity scores to fall outside AD (Figure S3 showed the examples of out-of-domain compound structures). To minimize the loss of valuable chemical space, the value of k and Z were set to 6 and 0.9, respectively. Accordingly, the final distance threshold was 8.375. As a result, the XGBoost model and SVM model performed best. For instance, XGB_combined and SVM_combined model yielded the similar performance with R² value of 0.663, 0.667 and RMSE of 0.479, 0.475, respectively.

Table 1 The performance of different models considering the applicability domain on the test set

Model validation on the Shanghai Qilu’s in-house dataset

We next assessed how well machine learning models trained on public data can be transferred to pharmaceutical industry data by applying them to the Shanghai Qilu’s in-house dataset. We also compared the established best-performing model with existing Caco-2 prediction model in the ADMET prediction platform ADMETlab 2.0 [42]. The regression performance statistics for prospective validation of aforementioned prediction models on the Shanghai Qilu’s in-house dataset were given by Table 2 and further visualized by Fig. 5. In general, the performance of these models on the external test set was significantly lower than that on the test set in terms of R² value. There may be two reasons that lead to such a poor performance and negative R² values. First, the similarity of molecules between the Shanghai Qilu’s internal dataset and training set is low. As illustrated in Figure S4, the highest Tanimoto similarity was less than 0.45. Second, most compounds in the Shanghai Qilu’s internal data set belong to the low/medium permeability class (Figure S5), while the high permeability bias in the curated data set tends to overpredict the low and medium permeable compounds. However, particular models showed a high correlation between their prediction results and experimental values with Pearson’s r of at least 0.6 and Spearman’s rho of at least 0.6. More specifically, two boosting algorithms (XGBoost, GB) coupled with ECFPs outperformed RF in all cases in terms of Pearson’s r, Spearman’s rho, MAE and RMSE. SVM seemed to work well with all representations but did not outperform the boosting models. The optimal model XGB_fp and GB_fp achieved the Pearson’s r of 0.702, 0.719, Spearman’s rho of 0.704, 0.725, MAE of 0.629, 0.659 and RMSE of 0.771, 0.797, respectively. Similarly, the deep learning approach showed improvements across metrics of Pearson’s r and Spearman’s rho compared to the other machine learning models but did not outperform XGBoost and GB. Although in the last few years, it is a general belief that the increase of the dimension of descriptors gives models with better predictive performances [43]. Based on our findings, the extra addition of 2D descriptors conversely increases the redundancy of features and makes the model become more complicated, thus inversely reducing the prediction performance. Noticeable exceptions from this conclusion seemed to be the Project3, where GB_combined model was the best method. Additionally, the CombinedNet model which used a feature combination of the molecular graph and the molecular fingerprints achieved acceptable performance with the Pearson’s r of 0.677, Spearman’s rho of 0.677, MAE of 0.730 and RMSE of 0.873. However, it did not further enhance prediction performance when compared to boosting algorithms. We also investigated for each project, which of the machine learning methods worked best (see Table S3). Here only projects with more than 8 molecules were included. It can be noted that boosting models were advantageous for the majority of projects in general. Subsequently, we used these models to predict the in-domain compounds in the in-house data set and the performance was provided by Table S4 and Figure S6. It can be found that the GB_fp model had the best overall predictive performance, but it may not necessarily apply to specific project. For example, GB_combined model achieved the best performance in Project 3 (Table S5). However, it is worth to mention that the overall performance of all models has decreased when compared to Table 2. The results suggested that defining a strict similarity threshold correlating with decreased generalization error is challenging and potentially highly dependent on the explored chemical space and in vitro ADME property.

Table 2 The performance of different models on the Shanghai Qilu’s internal data set

Fig. 5

Correlation of logP_app measurements (x-axis) against machine learning model predictions (y-axis) across 7 Shanghai Qilu’s internal projects (orange part of each subplot: low permeable; blue: medium permeable; green: high permeable)

Together, these results suggested that (1) the boosting model presents a favorable option for the prediction of Caco-2 permeability, (2) the public model transfer is observed with a certain performance decrease but is still valuable for prospective predictions, (3) deep learning models outperforms few methods investigated in this study on industrial data but may not be the best solutions, and (4) it is worth keeping in mind that the feature combination becomes effective only when the added features indeed complement the lacked information and sometimes less can be more.

Representative MMP transformations of Caco-2 permeability

The diffusion of a compound is a complex process influenced by various types of interactions [4]. In an attempt to understand how different functional groups affect in Caco-2 permeability, we applied the open-source package mmpdb version 2 to generate a MMP knowledge base tool using all the 5654 compounds in the curated dataset. The three-cut method was used for the analysis, resulting in 65,440 matched molecular pairs and a total of 415,716 rules. A representative set of MMP rules reflecting common medicinal chemistry transformations was listed in Fig. 6 for the optimization of Caco-2 endpoints. For this analysis, we included only transformations that contained at least 10 matched pairs and showed a significant property shift as measured with a p-value < 0.05. In the Fig. 6, “MMP transformations” refers to the transform rules of molecular pairs, “ΔCaco-2 ± std” denotes the average difference and standard deviation in the value of logP_app of molecular pairs, and “nPairs” is the number of molecular pairs. Entries in Fig. 6 can broadly be classified in two groups. The first group involves the introduction of hydrocarbon groups, which tends to increase the logP_app values. Because of the presence of phospholipid bilayer, high lipophilicity of a compound is favorable for efficient permeability. The second group includes transformations of polar groups. As we know, high polarity is detrimental to permeability [44]. Removing or masking hydrogen-bond donors seems to be a good tactic to improve permeability. The introduction of heteroatoms (e.g., O, N) tends to reduce logP_app. Additionally, the halogen substitution in the molecules may play an important role in increasing the polarization and thus decreasing the water solubility and absorption of the molecule. In conclusion, the above rules are generally consistent with our chemical intuition that Caco-2 permeability is related to key physicochemical parameters of molecules, such as lipophilicity, ionization, hydrogen bonding, and molecular size [45]. We believe that these rules will be an important assistive tool for Caco-2 prediction and provide a reference for medicinal chemists to optimize lead compounds.

Fig. 6

Representative medicinal chemistry transforms with corresponding changes of Caco-2 permeability

Conclusions

Understanding Caco-2 permeability is crucial for predicting drug behavior within the gastrointestinal tract and guiding decisions in the early stages of drug discovery and development. This study systematically evaluated Caco-2 permeability using machine learning models trained on public data. Regression models for Caco-2 permeability were developed using four machine learning methods (RF, SVM, GBM, and XGBoost) and two deep learning methods (DMPNN and CombinedNet). Various molecular representations, including ECFP4, RDKit2D descriptors, and combinations thereof (ECFP4 + RDKit2D), were explored to assess their impact on model performance. Our results reinforced previous evidence suggesting that deep learning methods currently do not consistently outperform simple machine learning methods in drug discovery applications [46]. Boosting models demonstrated superior overall performance compared to deep learning models across both training and test datasets. Approaches leveraging human-engineered molecular descriptors surpassed graph-based deep learning methods in predictive accuracy, while the deep learning model coupled with ECFPs yielded comparable results to machine learning methods. Furthermore, the models trained on public data were transferred and evaluated on Shanghai Qilu’s in-house dataset, revealing that particular models such as XGBoost_fp, GBM_fp, and CombinedNet maintained robust predictive quality to some extent with Pearson’s r and Spearman’s rho of about 0.65 to 0.70, despite potential differences in data distributions between public and industry databases. Additionally, the contributions of common substituents to logP_app were elucidated using Matched Molecular Pair Analysis (MMPA), offering insights into structural optimization strategies.