Experimental

Dataset and endpoint definition

A dataset comprising 62 structurally diverse compounds with reported BSEP inhibition data (IC50 values) was compiled from published in vitro vesicle-based transporter assays and curated datasets, primarily from Morgan et al. [6], Pedersen et al. [7], Montanari et al. [8], and Dawson et al. [18]. When multiple values were available, preference was given to data obtained under comparable experimental conditions, and inconsistent entries were excluded to ensure dataset reliability. The endpoint was expressed as pBSEP, defined as the negative logarithm of the experimental BSEP inhibition metric (IC50), standardized to a consistent molar concentration unit (M). Prior to analysis, chemical structures were curated to ensure consistency in representation.

The dataset was divided into a training set (n = 50) and an external test set (n = 12) using response-based stratified sampling. Specifically, compounds were first sorted by their pBSEP values, and samples were then randomly assigned to the two subsets to ensure that both subsets adequately covered the full range of the response variable. Summary statistics confirmed comparable distributions between the training and test sets in terms of range, mean and variability. The training set was used for model development and internal validation, while the test set was used exclusively for external validation. Literature pBSEP data for each set are included in Table S1 in the Supplementary material.

Biomimetic chromatographic descriptors

Biomimetic chromatographic descriptors related to phospholipid affinity and plasma protein binding were determined using established HPLC-based methodologies described in the literature [14-16]. These descriptors were selected for their mechanistic relevance to drug-membrane interactions and protein binding, which are key factors influencing BSEP inhibition. The derived chromatographic descriptors, CHI IAM and log _k_HSA, were used as proxies for membrane partitioning and plasma protein binding, respectively, in subsequent QSAR modelling of BSEP inhibition. All compounds were obtained from commercial sources (Sigma-Aldrich, Merck) and prepared as 10 mM stock solutions in dimethyl sulfoxide (DMSO). Prior to analysis, samples were appropriately diluted and injected into an Agilent 1100 HPLC system equipped with a diode array detector (DAD) (Agilent Technologies, Santa Clara, CA, USA).

Phospholipid binding (immobilized artificial membrane chromatography)

Phospholipid affinity was evaluated using an immobilized artificial membrane (IAM) column (IAM.PC.DD2, 100×4.6 mm), which mimics the phospholipid environment of biological membranes [14,16]. Chromatographic measurements were performed using a 50 mM ammonium acetate buffer (pH 7.4) as the aqueous phase and acetonitrile as the organic modifier. A gradient was applied to reach 90 % acetonitrile within 4.75 min, followed by a short hold and re-equilibration, resulting in a total run time of 6 min. Retention times were calibrated using a set of reference compounds with established CHI IAM values reported in the literature. The calibration exhibited excellent linearity between gradient retention times and CHI IAM values (_R_2 > 0.99), confirming the method's robustness. Measurement reproducibility was within ±0.005 min. The Chromatographic Hydrophobicity Index for IAM (CHI IAM) approximates the effective organic modifier concentration at elution and serves as a descriptor of phospholipid binding [14,16]. CHI IAM descriptors for this study’s dataset are included in Table S2 in the Supplementary material.

Plasma protein binding (human serum albumin chromatography)

Plasma protein binding was assessed using a Chiralpak HSA column (50×3 mm, 5 μm). The mobile phase consisted of 50 mM ammonium acetate buffer (pH 7.4), with isopropanol as the organic modifier, and the gradient reached 35 % within 3.5 min [14,15]. The total analysis time was 6 min, including re-equilibration. Retention times were calibrated using reference compounds with known plasma protein-binding values reported in the literature. The calibration demonstrated strong linear correlation (_R_2 > 0.99), supporting the reliability of the derived descriptors. The reproducibility of retention time measurements was within ±0.01 min. The chromatographic retention data were converted to logarithmic retention factors (log _k_HSA), which provide an estimate of a compound's affinity for human serum albumin [15]. Chromatographic HSA affinity descriptors are summarized in Table S3 in the Supplementary material.

In silico descriptors

Calculated physicochemical descriptors were obtained using the widely used and validated ACD/Percepta Platform (Advanced Chemistry Development, Toronto, Canada) for predicting physicochemical properties. The selected descriptors included lipophilicity (log P), distribution coefficient at physiological pH (log _D_7.4), molecular weight (MW), hydrogen-bond acceptors (HBA), hydrogen-bond donors (HBD), ionization state (net charge at pH 7.4), Abraham solvation parameters (A and B), and topological polar surface area (TPSA). These descriptors were selected for their established relevance and statistical significance in quantitative structure-activity/property/toxicity relationship (QSAR/QSPR/QTAR) studies, particularly in modeling drug distribution, membrane interactions, and transporter-related effects. The conventional lipophilicity descriptors log P and log _D_7.4 were included as reference physicochemical parameters for comparison with the experimentally derived biomimetic descriptors. Log P describes the partitioning behavior of the neutral form of a compound and therefore does not account for ionization at physiological pH. Log _D_7.4 partially addresses this limitation by representing the apparent distribution of all ionization states at pH 7.4. However, both descriptors remain bulk physicochemical parameters and do not explicitly capture microenvironmental interactions such as phospholipid affinity or plasma protein binding. All calculations were performed using standardized molecular structures under default software settings. The complete set of calculated descriptors for all compounds is provided in Table S4 in the Supplementary material.

Model development

Quantitative structure-activity relationship (QSAR) models were developed using multiple linear regression (MLR) implemented in IBM SPSS Statistics v27 (IBM Corp., Armonk, NY, USA). Model development was restricted to the training set. Several descriptor combinations were explored, including biomimetic descriptors (CHI IAM, log _k_HSA), hybrid models incorporating physicochemical data, and baseline models based on conventional lipophilicity descriptors (log P, log _D_7.4). Model selection was based primarily on predictive performance, particularly external validation, rather than solely on goodness-of-fit, in accordance with established QSAR modelling and validation principles [19]. The general pattern of the model is presented by Equation (1):

where ŷ is the predicted pBSEP, _b_₀ is the intercept, _b_ᵢ are the regression coefficients and _x_ᵢ are the selected descriptors.

Cross-validation

Internal validation was performed using 5-fold cross-validation on the training set [19]. The dataset was partitioned into five subsets, and each subset was used once as validation data, with the remaining subsets used for model training. Predictions from all folds were combined to obtain cross-validated predictions.

The cross-validated predictive coefficient (_Q_2cv) was calculated by Equation (2):

where _y_ᵢ represents the observed values, _ŷ_i,cv the predicted values obtained from cross-validation, and _ȳ_train the mean of the observed values in the training set.

The root-mean-square error of cross-validation (RMSECV) was calculated by Equation (3):

where _n_train is the number of compounds in the training set.

External validation

The final model was developed using the entire training set and applied to the external test set.

The external predictive coefficient (_Q_2ext) was calculated by Equation (4):

where _ŷ_ᵢ,test are the predicted values for the external test set.

The root-mean-square error of prediction (RMSEP) was calculated by Equation (5):

where _n_test is the number of compounds in the external test set.

Goodness-of-fit

The coefficient of determination (_R_2) was calculated by Equation (6):

where _y_ᵢ represents the observed values, _ŷ_ᵢ the predicted values from the model, and _ȳ_train the mean of the observed values in the training set.

The adjusted coefficient of determination (_R_2adj) was calculated by Equation (7):

where n is the number of compounds in the training set, and p is the number of descriptors included in the model. Additional statistical parameters, including the correlation coefficient (R), F-statistic, and standard error of estimate (s), were evaluated to assess model robustness and statistical significance.

Residual analysis and Y-randomization

Residual analysis was performed by plotting standardized residuals (_r_ᵢ) against predicted values. Standardized residuals, expressed in units of standard deviation, were expected to be randomly distributed around zero with no systematic patterns. A random scatter of residuals without discernible trends was considered indicative of model adequacy, linearity, and absence of heteroscedasticity.

Y-randomization (response permutation) was performed by randomly shuffling the response variable while keeping the descriptor matrix unchanged [19]. Multiple randomized models were generated and compared with the original model. A significant reduction in _R_2, _Q_2, and other predictive metrics for the randomized models confirmed that the original model was not due to chance correlation. Additionally, the intercept values of the randomized models were examined to ensure the absence of systematic bias.

Applicability domain

The applicability domain of the final model was evaluated using the Williams plot, where standardized residuals (_r_ᵢ) were plotted against leverage values (_h_ᵢ) [19].

The warning leverage threshold (h*) was defined by Equation (8):

where p is the number of descriptors in the model and _n_train is the number of compounds in the training set.

Compounds with standardized residuals |_r_ᵢ| > 3 were considered response outliers, while compounds with leverage values _h_ᵢ > h* were considered structurally influential. The absence of compounds exceeding both thresholds supported the reliability and robustness of the model predictions within its applicability domain. Compounds falling within these boundaries were within the model’s applicability domain and therefore associated with reliable predictions.