Results and discussion

Dataset characteristics

The dataset comprised 62 structurally diverse compounds with reported BSEP inhibition values expressed as pBSEP. While the dataset size is modest, it is consistent with experimentally derived descriptor studies, in which measurement fidelity and mechanistic interpretability are prioritized over dataset size. The response variable exhibited a broad distribution, ranging from 3.00 to 5.60, with a mean value of 4.42 and a standard deviation of 0.68, indicating substantial variability suitable for model development.

The dataset was divided into a training set (n = 50) and an external test set (n = 12) using a response-based stratified sampling approach. The distribution of pBSEP values in both subsets was comparable, with mean values of 4.44 and 4.29 and standard deviations of 0.68 and 0.70 for the training and test sets, respectively. As illustrated in the boxplot (Figure 1), both subsets covered a similar range of inhibitory potency, confirming the absence of systematic bias and supporting the representativeness of the dataset split.

Figure 1.: Boxplot representation of pBSEP values for the training (n = 50) and external test (n = 12) sets. The similar distributions, central tendencies, and spreads indicate a balanced dataset split and adequate coverage of the response space

The biomimetic chromatographic descriptors also exhibited significant variability across the dataset. The CHI IAM values ranged from −10.28 to 62.02 (mean 35.77 ± 14.39), reflecting a wide spectrum of phospholipid affinity, from weak to strong membrane interactions. Similarly, log _k_HSA values spanned from −1.48 to 1.88 (mean 0.97 ± 0.64), indicating substantial differences in plasma protein binding among the compounds.

Inspection of the descriptor limits further supported the dataset's chemical plausibility. The lowest CHI IAM values were observed for isoniazid and salicylic acid, consistent with their small size, polarity, and/or ionizable character and reduced membrane-partitioning tendency. Conversely, the highest CHI IAM value was observed for amiodarone, followed by imatinib, atorvastatin, fluoxetine, and miconazole, which are characterized by greater lipophilicity, aromaticity, and structural complexity [14,16,17]. For log _k_HSA, the lowest values were observed for isoniazid, caffeine, sulpiride, erythromycin, and pravastatin, reflecting comparatively weaker albumin-binding propensity. The highest log _k_HSA value was observed for benzbromarone, followed by ibuprofen, diclofenac, amiodarone, miconazole, and sulfonylureas, including glibenclamide/glyburide and glimepiride. These compounds are generally aromatic, lipophilic, and/or acidic/anionizable, consistent with stronger interactions with HSA [15,20]. These trends indicate that the boundary values of the biomimetic descriptors correspond to chemically interpretable structural features rather than arbitrary outliers.

The broad distribution of both the response variable and the selected descriptors ensure adequate coverage of chemical and physicochemical space, which is essential for developing robust, predictive QSAR models. Overall, the dataset was considered well-suited for modelling BSEP inhibition. No extreme skewness or abnormal distribution patterns were observed that would compromise model stability.

Quantitative structure-activity relationship model development and comparison

A series of QSAR models was developed to evaluate the predictive performance of biomimetic chromatographic descriptors in comparison with conventional physicochemical parameters. The results are summarized in Table 1.

Model | Equation | R | R 2 | R 2 adj | s | F
M1 [CHI IAM] | pBSEP = 2.947 + 0.040·CHI IAM | 0.838 | 0.702 | 0.696 | 0.372 | 113
M2 [CHI IAM + log kHSA] | pBSEP = 2.902 + 0.034·CHI IAM + 0.277·log kHSA | 0.865 | 0.748 | 0.737 | 0.346 | 69.7
M3 [CHI IAM + log kHSA + MW] | pBSEP = 2.646 + 0.029·CHI IAM ++ 0.297·log kHSA + 0.001·MW | 0.898 | 0.806 | 0.794 | 0.307 | 63.8
M4 [log P] | pBSEP = 3.524 + 0.250·log P | 0.599 | 0.359 | 0.346 | 0.546 | 26.9
M5 [log P + MW] | pBSEP = 3.112 + 0.195·log P + 0.001·MW | 0.692 | 0.479 | 0.457 | 0.498 | 21.6
M6 [log D7.4] | pBSEP = 3.971 + 0.209·log D7.4 | 0.583 | 0.340 | 0.326 | 0.554 | 24.7
M7 [log D7.4 + MW] | pBSEP = 3.452 + 0.161·log D7.4 + 0.002·MW | 0.682 | 0.465 | 0.442 | 0.504 | 20.4

The model based solely on the phospholipid affinity descriptor CHI IAM (M1) already demonstrated strong predictive ability (_R_2 = 0.702), indicating that membrane partitioning plays a significant role in BSEP inhibition. The positive coefficient of CHI IAM suggests that compounds with higher phospholipid affinity exhibit increased inhibitory potency, consistent with enhanced accumulation in membrane environments where BSEP is located.

Incorporating plasma protein binding via log _k_HSA (M2) led to a notable improvement in model performance (_R_2 = 0.748), highlighting the additional contribution of protein binding to BSEP inhibition. The positive contribution of log _k_HSA suggests that compounds with higher affinity for serum albumin tend to exhibit increased inhibitory potential, likely reflecting the interplay between binding, distribution, and effective exposure at the hepatocellular membrane. The log _k_HSA should not be interpreted as a causal determinant of BSEP inhibition, but rather as an experimentally derived surrogate descriptor that captures systemic distribution characteristics that covary with membrane exposure and transporter interactions.

Further inclusion of molecular weight (MW) in the hybrid model (M3) yielded the highest predictive performance (_R_2 = 0.806), suggesting that size-related effects may also influence BSEP inhibition, potentially through steric or permeability-related mechanisms. However, the relatively modest improvement compared to M2 suggests that the primary determinants of BSEP inhibition are captured by the biomimetic descriptors.

In contrast, models based on conventional lipophilicity descriptors showed substantially lower predictive performance. The log P-based model (M4) yielded _R_2 = 0.359, while the log _D_7.4 model (M6) showed similar performance (_R_2 = 0.340), indicating that traditional bulk lipophilicity descriptors alone are insufficient to describe BSEP inhibition. The inclusion of molecular weight slightly improved these models (M5 and M7), but their performance remained significantly inferior to that of the biomimetic models.

The weaker performance of the conventional lipophilicity-based models may partly reflect the phenomenological limitations of log P and log _D_7.4. Although log P is widely used in QSAR modelling, it describes the partitioning behaviour of the neutral form of a compound and does not account for ionization at physiological pH. Since many drugs are partially or fully ionized under biological conditions, log P may not accurately reflect their in vivo distribution properties. Log _D_7.4 partially addresses this limitation by incorporating the apparent distribution of ionized and unionized species at physiological pH; however, it remains a simplified bulk descriptor. It does not explicitly capture microenvironmental interactions such as phospholipid affinity, membrane-specific partitioning, or plasma protein binding. These interactions may be influenced by ionization state, charge distribution, aromaticity, and molecular structure [14-16]. In contrast, biomimetic chromatographic descriptors inherently reflect experimentally observable interactions with membrane-like and protein-binding environments, contributing to their improved predictive performance and mechanistic relevance in the present study [14-17].

These findings highlight a key limitation of conventional descriptors, which primarily reflect bulk partitioning properties but fail to capture microenvironmental distribution phenomena, such as membrane affinity and protein binding [14-16]. In contrast, CHI IAM and log _k_HSA provide experimentally derived proxies for membrane partitioning and plasma protein binding, respectively, enabling a more mechanistically relevant description of compound behaviour in biological systems.

The superior performance of biomimetic models demonstrates that experimentally derived chromatographic descriptors can significantly improve predictions of BSEP inhibition. The results support the hypothesis that local exposure at the membrane interface, rather than bulk physicochemical properties alone, is a critical determinant of transporter inhibition. This is consistent with the broader view that transporter-mediated toxicity is governed by the interplay between intrinsic potency and effective exposure at the site of interaction [1,4].

Model validation and predictive performance

The predictive performance of the developed QSAR models was evaluated using both internal (Q2cv, RMSECV) and external (Q2ext, RMSEP) validation metrics (Table 2).

Model | R 2 | Q 2 cv | Q 2 ext | RMSECV | RMSEP
Phospholipid Binding (M1: CHI IAM) | 0.702 | 0.671 | 0.664 | 0.384 | 0.398
Biomimetic (M2: CHI IAM + log kHSA) | 0.748 | 0.696 | 0.667 | 0.369 | 0.370
Hybrid (M3: CHI IAM + log kHSA + MW) | 0.806 | 0.742 | 0.578 | 0.340 | 0.435
Conventional (M5: log P + MW) | 0.479 | 0.327 | 0.588 | 0.549 | 0.440
Conventional (M7: log D7.4 + MW) | 0.465 | 0.278 | 0.500 | 0.568 | 0.488

The biomimetic model based on CHI IAM and log _k_HSA demonstrated strong and well-balanced performance, with _R_2 = 0.748, _Q_2cv = 0.696, and _Q_2ext = 0.667. The close agreement between internal and external validation metrics indicates good model stability and generalizability, suggesting that the model is not overfitted and retains predictive power for unseen compounds.

To further evaluate model parsimony, the predictive performance of the IAM-only model (M1) was compared with that of the two-descriptor biomimetic model (M2). The CHI IAM descriptor alone provided strong predictive capability, with _R_2 = 0.702, _Q_2cv = 0.671, _Q_2ext = 0.664, RMSECV = 0.384, and RMSEP = 0.398, indicating that membrane affinity represents the dominant biomimetic determinant of BSEP inhibition in the present dataset. Addition of log _k_HSA modestly improved the training and internal validation statistics, increasing _R_2 to 0.748 and _Q_2cv to 0.696, with a slightly lower RMSECV of 0.369. However, this improvement did not translate into a substantial gain in external predictive performance, as M2 showed a very similar _Q_2ext value of 0.667 and RMSEP of 0.370. These results suggest that the IAM-only model may represent the most practical first-tier screening approach, particularly when minimizing experimental time and cost is important. The combined IAM and HSA model remains mechanistically informative because it integrates membrane affinity with albumin-binding propensity, but its use should be justified by the need for broader biomimetic characterization rather than by a clear improvement in external prediction. This comparison highlights the trade-off between experimental simplicity and mechanistic coverage in chromatography-based modelling of BSEP inhibition.

The hybrid model, which additionally incorporates molecular weight, achieved the highest goodness-of-fit (_R_2 = 0.806) and internal predictive performance (_Q_2cv = 0.742). However, its external predictive ability was lower (_Q_2ext = 0.578, RMSEP = 0.435), indicating reduced generalizability relative to the simpler biomimetic model. This behaviour suggests that including molecular weight, although it improves model fit, may introduce overfitting or capture dataset-specific variance that does not generalize to external compounds.

In contrast, models based on conventional physicochemical descriptors showed weaker internal predictive performance, with _Q_2cv values of 0.327 for the log P-based model and 0.278 for the log _D_7.4 model. Although the log P-based model exhibited moderate external predictivity (_Q_2ext = 0.588), the discrepancy between internal and external validation metrics indicates limited robustness and reduced reliability. The log _D_7.4-based model showed consistently lower performance across both internal and external validation metrics.

Among the multi-descriptor models, the biomimetic model provided the best balance among goodness-of-fit, robustness, interpretability, and external predictive performance. At the same time, the IAM-only model showed comparable external predictive performance, supporting its potential as a simpler first-tier screening tool. These findings highlight the advantage of experimentally derived chromatographic descriptors in capturing mechanistically relevant aspects of compound distribution, such as membrane affinity and protein binding, which are not adequately described by conventional lipophilicity parameters [14]. Most published models for BSEP inhibition are based on classification approaches and report metrics such as AUC and accuracy rather than continuous predictive performance [9-12]. In this context, the present model provides a quantitative prediction of inhibitory potency (_Q_2ext = 0.667) and offers enhanced mechanistic interpretability through biomimetic descriptors.

Importantly, the results support the concept that BSEP inhibition is not solely governed by intrinsic potency or bulk physicochemical properties but is strongly driven by membrane-proximal exposure [1,4]. By incorporating experimentally derived descriptors that reflect these processes, the biomimetic model offers improved predictive performance and mechanistic interpretability compared to traditional QSAR approaches [14-16].

To further assess the robustness and reliability of the selected biomimetic model, graphical validation techniques and applicability domain analysis were performed.

The observed-versus-predicted plot for the biomimetic model (Figure 2) demonstrates strong agreement between experimental and predicted pBSEP values for both the training and external test sets. Most data points are clustered near the identity line, indicating that the model accurately captures the relationship between the selected descriptors and BSEP inhibition.

Figure 2.: Observed versus predicted pBSEP values for the IAM + HSA biomimetic model. The dashed line represents the line of identity, indicating perfect agreement between experimental and predicted values. Training set compounds are shown in grey, while external test set compounds are highlighted in orange. Marker shapes indicate compound charge class

The training-set compounds are tightly clustered around the diagonal, reflecting good model fit and consistency across the calibration dataset. Importantly, the test set compounds also follow the same trend, confirming the model’s ability to generalize to unseen data and supporting its external predictive performance. Although a small number of compounds deviate slightly from the ideal line, no systematic overestimation or underestimation is observed across the range of pBSEP values. The deviations appear randomly distributed, suggesting that prediction errors are not associated with specific regions of the response space. Furthermore, the model maintains good predictive performance across compounds with different ionization states, as indicated by the distribution of acids, bases, and neutral molecules in the plot. This suggests that the biomimetic descriptors effectively capture underlying distributional properties relevant across diverse chemical classes. The absence of systematic deviation at higher or lower pBSEP values further indicates that the model does not suffer from range-dependent bias.

Residual analysis (Figure 3) shows that standardized residuals are randomly distributed around zero across the full range of predicted pBSEP values. No systematic patterns, trends, or curvature are observed, indicating that the assumptions of linearity and homoscedasticity are satisfied.

Figure 3.: Standardized residuals versus predicted pBSEP values for the biomimetic model (CHI IAM + log kHSA) for both training and test sets

The residuals are symmetrically distributed around the zero line, with most values falling within the acceptable range (|_r_ᵢ| ≤ 2), and all values remaining within the critical threshold (|_r_ᵢ| ≤ 3). This suggests the absence of extreme outliers and confirms the model's stability. Importantly, no increase in residual variance is observed for higher or lower predicted values, indicating that the model does not exhibit heteroscedasticity. The consistent spread of residuals across the prediction range supports the model's reliability for both weak and strong inhibitors.

The applicability domain of the biomimetic model was evaluated using the Williams plot (Figure 4), where standardized residuals were plotted against leverage values. Most compounds are located within the defined boundaries (|_r_ᵢ| ≤ 3 and _h_ᵢ ≤ h* = 0.18), indicating that predictions for most compounds are reliable and fall within the model’s domain.

Figure 4.: Williams plot for the biomimetic model (CHI IAM + log kHSA), showing standardized residuals versus leverage values (hᵢ). The horizontal lines represent the ±3 residual limits, and the vertical line indicates the warning leverage threshold (h* = 0.18)

A small number of compounds, namely isoniazid, ibuprofen, and salicylic acid, exhibited leverage values exceeding the warning threshold (_h_ᵢ > h*), indicating that they occupy influential regions of the descriptor space. However, their standardized residuals remained within acceptable limits (|_r_ᵢ| ≤ 3), suggesting that they are not response outliers and that their predicted values remain consistent with the model [19]. The high leverage observed for ibuprofen and salicylic acid can be rationalized based on their structural and physicochemical characteristics. Both compounds are relatively small aromatic carboxylic acids within the acidic NSAID/salicylate-like chemical space.

Their ionizable acidic functionality, combined with aromatic character and appreciable albumin-binding propensity, gives them a distinctive biomimetic profile [15,20]. In particular, these compounds combine comparatively low-to-moderate IAM affinity with relatively high HSA interaction, positioning them at the edge of the descriptor space relative to the broader dataset [14-16]. Other NSAIDs present in the dataset, such as diclofenac, indomethacin, etodolac, and sulindac, exhibit similar acidic/aromatic features but do not exceed the leverage threshold because their IAM and HSA descriptor values lie closer to the model's main chemical space. Therefore, ibuprofen and salicylic acid should be interpreted as edge cases within the acidic NSAID-like region rather than as evidence of systematic class-wide outlier behaviour. Importantly, no compounds exceeded both the leverage and residual thresholds simultaneously, confirming the absence of influential outliers and supporting the model's robustness within its applicability domain [19].

To assess the possibility of chance correlation, Y-randomization tests were performed by randomly permuting the response variable while keeping the descriptor matrix unchanged. The models generated from randomized data exhibited significantly lower _R_2 and _Q_2 values than the original model (Table 3).

Model | R 2 | Q 2
Original | 0.748 | 0.696
Random 1 | 0.12 | -0.05
Random 2 | 0.08 | -0.12

In particular, the randomized models showed very low _R_2 values (≤ 0.12) and negative or near-zero _Q_2 values, indicating a complete loss of predictive ability. In contrast, the original biomimetic model demonstrated substantially higher performance (_R_2 = 0.748, _Q_2 = 0.696), confirming that the observed correlation is not due to random effects. These results strongly support the statistical robustness of the developed QSAR model and confirm that its predictive performance arises from meaningful relationships between descriptors and response rather than chance correlation.

Mechanistic interpretation

The superior performance of the biomimetic model can be interpreted in the context of BSEP biology and the distributional behaviour of drug molecules. BSEP is an ATP-dependent efflux transporter located at the canalicular membrane of hepatocytes and represents the rate-limiting step in bile acid secretion. Because inhibition occurs at a membrane-embedded target, local exposure at the membrane interface is likely to be more relevant than bulk physicochemical properties alone. The International Transporter Consortium has emphasized that the interpretation of BSEP inhibition should take into account in vivo exposure and local concentrations rather than relying solely on in vitro inhibitory potency [1,4].

The positive contribution of CHI IAM is mechanistically consistent with this view. IAM chromatography is widely used as an experimental surrogate for phospholipid affinity and membrane partitioning, and recent reviews describe IAM retention as informative for permeability, tissue binding, and off-target binding processes [14-16]. Thus, higher CHI IAM values can reasonably be interpreted as indicating a greater tendency of a compound to partition into membrane-like environments, potentially favouring local enrichment near the canalicular membrane, where BSEP resides.

This interpretation is also supported by the BSEP literature, which shows that inhibition is strongly associated with lipophilicity. Pedersen et al. reported that BSEP inhibition correlates strongly with compound lipophilicity, whereas positive molecular charge is associated with reduced inhibition. They further note that transported BSEP substrates are preferably monovalent, negatively charged bile acids, and the few known non-bile acid substrates also carry a negative net charge at physiological pH [7]. Taken together, these observations support the idea that membrane-affine, non-cationic compounds are more likely to occupy the physicochemical space compatible with BSEP interaction.

The inclusion of log _k_HSA further improved model performance, although its role requires careful interpretation. Albumin binding itself is not a direct mechanistic driver of BSEP inhibition; rather, log _k_HSA should be viewed as an experimentally derived descriptor of a broader distributional phenotype [15]. Compounds with higher albumin affinity often share physicochemical features, such as hydrophobicity and anionic character, that are also associated with hepatobiliary distribution and transporter interaction. Human serum albumin is the major plasma carrier protein and binds a wide range of drugs, thereby limiting their free concentrations and modulating systemic transport and distribution [15]. Structural and biochemical studies show that HSA contains major drug-binding pockets, including the Sudlow sites I and II; site II preferentially binds aromatic compounds such as ibuprofen, whereas site I binds ligands such as warfarin [20]. These complementary mechanistic roles of CHI IAM and log _k_HSA are schematically summarized in Figure 5.

Figure 5.: Mechanistic framework linking biomimetic chromatographic descriptors to BSEP inhibition. CHI IAM reflects phospholipid affinity and membrane partitioning, while log kHSA captures protein-associated distribution and albumin binding. Together, these descriptors approximate local exposure at the canalicular membrane, providing a mechanistic basis for predicting BSEP inhibition

More generally, HSA is recognized as an important carrier for many acidic drugs [15,20]. In this sense, log _k_HSA complements CHI IAM by capturing systemic distribution behaviour that covaries with membrane exposure. The combination of these descriptors therefore provides a more complete representation of microenvironmental exposure conditions relevant to BSEP interaction than conventional descriptors such as log P or log _D_7.4 alone [1,4,14].

The present findings support a mechanistically plausible framework in which BSEP inhibition is favoured by compounds that combine appreciable membrane affinity with distribution characteristics typical of albumin-binding. While these descriptors do not replace transporter assays, they provide experimentally grounded proxies for the physicochemical and distributional determinants that govern BSEP interactions, thereby supporting their use in early-stage risk assessment. This framework may therefore support early prioritization of compounds with reduced cholestatic liability during drug discovery.