Experimental

Phospholipidosis data collection

The phospholipidosis dataset was compiled from previously published in vitro studies that employed fluorescence-based and high-content screening assays to detect drug-induced phospholipid accumulation [9-11,25-27]. The dataset included compounds evaluated in mammalian cell systems commonly used for PLD screening, including hepatocyte-derived HepG2 cells and macrophage models such as RAW264.7, I-13.35, and primary mouse macrophages [25-27]. These models are relevant because of their sensitivity to lysosomal phospholipid accumulation and their established use for identifying phospholipidosis hazards. PLD induction was assessed using fluorescent phospholipid probes, including NBD-PE and LipidTOX™, which accumulate in lysosomal compartments and serve as experimental surrogates for intracellular phospholipid accumulation [10,11,25,26]. Quantitative endpoints were derived from fluorescence intensity measurements obtained by automated imaging or plate-based detection systems and reflected the extent of intracellular phospholipid accumulation. For each compound, concentration-response data were used to derive potency values expressed as EC₅₀. When multiple experimental values were available for the same compound, a predefined selection hierarchy was applied. First, priority was given to quantitative EC₅₀ values derived from concentration-response experiments that directly measured intracellular phospholipid accumulation using fluorescence-based or high-content PLD assays. Second, values were preferred when the concentration-response relationship was complete, and the reported potency value was clearly derived from a fitted dose-response curve rather than from a single-concentration response or qualitative classification. Third, among otherwise comparable values, preference was given to measurements obtained under assay conditions most consistent with the rest of the curated dataset, including comparable mammalian cell-based PLD endpoints and clearly reported exposure-response information. Values were not retained when the endpoint was not directly related to intracellular phospholipid accumulation, when the concentration-response information was incomplete, or when interpretation was clearly confounded by cytotoxicity or insufficient methodological detail. To enable direct comparison with physicochemical and chromatographic descriptors, EC₅₀ values reported in micromolar units were converted to molar units and transformed to a negative logarithmic scale: pEC₅₀ = -log₁₀ (EC₅₀ / M). Higher pEC₅₀ values, therefore, indicate greater PLD-inducing potency. The complete list of compounds, together with their corresponding potency values and PLD classifications, is provided in Supplementary material, Table S1.

Phospholipid binding (immobilized artificial membrane chromatography)

Immobilized artificial membrane chromatography was used to characterize the membrane affinity of the studied compounds. IAM stationary phases contain phosphatidylcholine analogs covalently immobilized on silica and provide a biomimetic environment that captures key aspects of drug-membrane interactions, including hydrophobic, polar and electrostatic contributions [21,23,24]. Chromatographic analyses were performed using an IAM.PC.DD2 column under gradient elution conditions, with ammonium acetate buffer at pH 7.4 as the aqueous phase and acetonitrile as the organic modifier. Retention times were converted to chromatographic hydrophobicity index (CHI) values, expressed as CHI IAM, using calibration with reference compounds in accordance with established biomimetic chromatography methodology [21,23,24]. The calibration showed excellent linearity, with _R_2 > 0.99. The resulting CHI IAM values were used as experimentally derived descriptors of membrane affinity.

Human serum albumin and α1-acid glycoprotein chromatography

Protein-binding-related distributional properties were characterized using biomimetic chromatography on immobilized human serum albumin (HSA) and α1-acid glycoprotein (AGP) stationary phases. These chromatographic descriptors were included to evaluate whether protein-binding-related retention provides complementary information to membrane-affinity measurements in PLD prediction. Chromatographic separations were performed using gradient elution with ammonium acetate buffer at pH 7.4 and isopropanol as the organic modifier. Retention data were converted to logarithmic retention factors (log _k_HSA and log _k_AGP) using calibration curves derived from reference compounds, following established biomimetic chromatography procedures [22-24]. The calibration curves showed excellent linearity, with _R_2 > 0.99. The resulting descriptors provide experimentally derived measures of drug-protein interactions related to systemic distribution and plasma protein-binding behaviour.

The experimentally determined biomimetic chromatographic descriptors, including CHI IAM, log _k_HSA, and log _k_AGP, are provided in Supplementary material, Table S2.

Physicochemical descriptors

A set of conventional physicochemical descriptors was collected for all compounds to enable comparison with biomimetic chromatographic parameters and to support modelling analyses. These descriptors included lipophilicity parameters (log P and log D at pH 7.4, log _D_7.4), MW, TPSA, and charge-related descriptors, namely the fractions of positively charged species, negatively charged species, and zwitterionic forms at physiological pH. Hydrogen-bonding properties were described using both structural and solvation-based descriptors. Specifically, HBD and HBA counts were included as simple structural descriptors, while Abraham’s solvation parameters for hydrogen bond acidity (A) and basicity (B) were also considered, providing a more quantitative representation of intermolecular hydrogen-bonding interactions. The combined use of these descriptors enables capture of both the presence and strength of hydrogen-bonding capacity, which is known to influence membrane partitioning, protein binding, and intracellular distribution. All physicochemical descriptors were calculated using the Percepta platform (ACD/Labs, Advanced Chemistry Development, Toronto, Canada) [28], which provides standardized computational estimates of molecular properties relevant to drug disposition. A complete list of calculated descriptors for all compounds is provided in Supplementary Material Table S3. These descriptors were included to evaluate the performance of traditional property-based models against biomimetic chromatographic descriptors and to assess their contribution to predicting phospholipidosis potential.

Dataset splitting and characterization

The dataset, comprising 65 compounds, was divided into a training set (n = 52) and an independent external test set (n = 13) using response-stratified sampling to ensure approximately 80:20 coverage of the pEC₅₀ range. This approach was selected to preserve the distribution of biological activity across both sets, thereby avoiding bias toward specific activity regions and supporting the development of predictive and generalizable models [29,30]. Descriptive statistics (mean, median, standard deviation, minimum, maximum, and range) were calculated to characterize the distribution of pEC₅₀ values, while compounds were split ordinally and their activity class distribution was assessed to ensure balanced representation between the training and test sets (Table S1 and Table S4). The distribution of pEC₅₀ values was examined using descriptive statistics and frequency histograms (Figure S1) to assess the consistency of the activity range and variability between the two sets. In parallel, the distribution of ionization states was assessed by classifying compounds into charge categories (base, weak base, neutral, acid, zwitterion) (Figure S2), as charge-related properties are known to influence membrane interactions and transporter behaviour. To evaluate the representativeness of the training and test sets in chemical space, PCA was performed on all biomimetic and physicochemical descriptors. This analysis was employed to ensure that the external test set falls within the descriptor space defined by the training set, thereby supporting reliable external validation and minimizing extrapolation.

Multiple linear regression

MLR models were developed to investigate the relationships among biomimetic chromatographic descriptors, physicochemical properties and phospholipidosis potential, expressed as pEC₅₀ [29,30]. These models used the training and external test sets described above and were fitted by ordinary least squares regression. Descriptor selection was based on mechanistic relevance, coefficient-level statistical significance, and avoidance of multicollinearity. Individual regression coefficients were evaluated using t-tests, and descriptors were retained in the presented MLR equations only when statistically significant at p < 0.05. Multicollinearity among descriptors was assessed using variance inflation factors (VIFs). For each descriptor, VIF was calculated by Equation (1):

where _R_j2 is obtained by regressing that descriptor against the remaining descriptors in the model.

Thus, VIF quantifies how much the variance of a regression coefficient is inflated by correlation with other predictors. Models with VIF values below 5 were considered acceptable and were retained. Model performance was evaluated using the coefficient of determination (_R_2), adjusted _R_2 (_R_2adj), and the F-statistic. Internal validation was performed using 5-fold cross-validation, yielding the cross-validated coefficient (_Q_2cv) and the RMSEcv. External predictive performance was assessed on the test set by calculating _Q_2ext and RMSEP. Definitions and formulas for the validation metrics used in this study, including _Q_2cv, _Q_2ext, RMSEcv and RMSEP, follow established QSAR validation recommendations [29,30].

Partial least squares

PLS regression was applied as a complementary multivariate method to assess the potential impact of descriptor interdependence on model performance. PLS is inherently robust to collinearity, as it projects the original variables onto latent components that maximize covariance with the response variable [31]. Models with suspected descriptor interdependence were developed using the same training and external test sets as those used in the MLR analysis. The optimal number of latent components was determined based on cross-validation (5-fold), selecting the model that maximized predictive performance (_Q_2cv) while avoiding overfitting. Model performance was evaluated using R2 for the training set, _Q_2cv for internal validation, and _Q_2ext for external validation, along with the corresponding RMSE values. The PLS analysis was used solely to confirm the robustness of the MLR findings with respect to descriptor interdependence.

Applicability domain

The applicability domain of the developed MLR models was evaluated using the leverage approach (Williams plot) [30]. Leverage values (_h_i) were calculated from the hat matrix, and a warning leverage (h) was defined as h = 3(p + 1)/n, where p is the number of model parameters, and n is the number of training compounds. Standardized residuals were calculated based on the standard deviation of the training residuals. Compounds with _h_i > h* were considered structurally influential, while those with standardized residuals outside ±3 were identified as response outliers.

Ordinal regression modelling

To account for the ordered nature of phospholipidosis severity, ordinal regression models were developed using a proportional-odds logistic regression framework [32]. This approach is appropriate when the response variable consists of ordered categories rather than independent nominal classes. In the present study, compounds were categorized according to pEC₅₀ into non-inducers (pEC₅₀ < 4), weak/moderate inducers (4 ≤ pEC₅₀ ≤ 5), and strong inducers (pEC₅₀ > 5). These thresholds correspond to PLD-induction EC₅₀ values of >100 μM, 10 to 100 μM and <10 μM, respectively, where EC₅₀ denotes the concentration producing half-maximal intracellular phospholipid accumulation under the assay conditions. The thresholds were selected as pragmatic, order-of-magnitude potency categories for early PLD risk stratification rather than as regulatory cutoffs. This categorization separates compounds that show no or only low-potency PLD induction at high micromolar concentrations from those that produce phospholipid accumulation at lower micromolar concentrations. The use of concentration-response-derived PLD potency values is consistent with previous fluorescence-based and high-content in vitro phospholipidosis screening studies [9-11,25-27]. The model estimates the cumulative probability that a compound belongs to a given class or a lower class, assuming that the effect of each descriptor is constant across class thresholds.

Model fit was assessed using log-likelihood, AIC, BIC and McFadden’s pseudo-_R_2 [32-35]. Log-likelihood is the logarithm of the probability of the observed classification data under the fitted model, given the estimated model parameters [32,35]. Higher, or less negative, LL values indicate better model fit, although LL was interpreted alongside AIC and BIC, as it is affected by model complexity. AIC and BIC were calculated by Equations (2) and (3) [33,34]:

where LL is the maximized log-likelihood, k is the number of estimated model parameters, and n is the number of observations.

McFadden’s pseudo-_R_2 was calculated by Equation (4) [35]:

where LLmodel and LLnull are the log-likelihoods of the fitted and intercept-only models, respectively.

McFadden’s pseudo-_R_2 expresses the improvement of the fitted model relative to a null model containing only intercept terms; it is not directly equivalent to the _R_2 used in linear regression but provides an indication of relative model fit. Predictive performance was evaluated by assigning each compound to the class with the highest predicted probability. Classification accuracy was calculated as the proportion of correctly classified compounds. Because the three PLD classes were not treated as interchangeable and class balance is important for risk stratification, macro-averaged F1 score (measure of the harmonic mean of precision and recall) was also calculated as the unweighted mean of the class-specific F1 scores. Thus, macro-F1 gives equal weight to non-inducers, weak/moderate inducers, and strong inducers.

Receiver operating characteristic analysis was performed using a one-vs-rest strategy for multiclass classification. In this approach, each PLD class is considered in turn as the positive class, with the remaining two classes combined as the negative class. The corresponding AUC values therefore describe the model's ability to discriminate each class from all others. Macro-AUC was calculated as the unweighted average of the class-specific one-vs-rest AUC values. Confusion matrices were generated for both the training and external test sets to evaluate the distribution of correct and misclassified classifications across ordered PLD categories.

Model comparison strategy

The performance of biomimetic, protein-binding, and conventional physicochemical models was systematically compared across both regression and classification frameworks. Emphasis was placed on evaluating whether biomimetic chromatographic descriptors (CHI IAM, log _k_HSA, log _k_AGP) offer improved predictive power over traditional descriptors such as log P and log D, while considering both statistical robustness and mechanistic interpretability. Model comparison was based on a combination of statistical robustness (_Q_2cv, _Q_2ext), predictive accuracy (RMSEP, classification metrics), and mechanistic interpretability.

For regression models, predefined acceptance criteria were used to support model interpretation. Retained descriptors were required to show coefficient-level statistical significance (p < 0.05), and multicollinearity was considered acceptable when VIF values were below 5. Models interpreted as predictive were expected to show stable internal validation and external predictive performance, as indicated by positive _Q_2cv and _Q_2ext values, preferably above 0.5, together with low RMSEcv and RMSEP relative to the dataset’s pEC₅₀ range. Applicability-domain assessment was also required, with most compounds expected to fall below the warning leverage threshold and within ±3 standardized residuals. Models that did not satisfy these predictive criteria were retained only as comparative or baseline models, not as primary predictive models. For ordinal regression models, acceptable performance required improvement over the null model, assessed by McFadden’s pseudo-_R_2, together with balanced class performance based on macro-F1 and discrimination ability based on macro-averaged one-vs-rest AUC. Ordinal models were considered useful for risk stratification when they exhibited balanced classification performance and avoided systematic confusion between distant PLD classes.

All statistical analyses and data visualization were performed using Python with the statsmodels and scikit-learn libraries.