Estimation of Biliary Excretion of Foreign Compounds Using Properties of Molecular Structure

Dataset

The biliary excretion dataset collated by Yang et al. (19), available at http://www.acsu.buffalo.edu/~memorris/Suppl%20Data_Rats_052909_pdb.pdf, plus some additional data from the original literature were used in this study (see supplementary material I for the dataset and the reference list). It consists of in vivo biliary excretion expressed as the percentage of dose excreted as the parent compound intact through the bile (BE%) for 217 compounds in rat after iv or intraperitoneal administration of the compound. The compounds are from different chemical classes such as bile acids, statins, dyes, penicillins and cephalosporins, macrolide antibiotics, quinolone antibiotics, NSAIDs, thrombin inhibitors, analgesics, anti-cancer drugs such as doxorubicin, folates, peptides, anti-HIV agents, quaternary ammoniums, sulphanilamide and arylaminosulphonic acids.

Molecular Descriptors

Molecular descriptors were calculated for the 217 compounds using ACD Labs/logD suite version 12, TSAR 3D version 3.3 (Accelryl Inc), molecular operation environment (MOE) version 2011.10 (Chemical Computing Group) and Symyx QSAR software (Accelryl Inc). The fractions of ionised compounds at pH 7.4 as acid (fiA), as base (fiB) or for zwitterionic compounds ionised as acid and base (fiAB) and fraction unionised (fU) were calculated using Eqs. 1 to 4 and used as additional molecular descriptors as follows:

1

2

3

4

In Eqs. 1 and 2, pKa was the most acidic and the most basic, respectively, which were obtained from ACD Labs pKa database and, in case the experimental pKa was not available, it was calculated by the software.

In MOE, following the wash procedure which removed salt forms, energy minimization was carried in order to calculate atomic coordinates corresponding to the local minimum (the low-energy conformation). Thereafter, self-consistent field calculations were performed, and this latter energy-minimised structure was used for the calculation of all the molecular descriptors. Before building the models, the molecular descriptors were checked to find and discard those columns containing more than 98% constant values or more than 28 (out of 217) missing values. The total number of descriptors used in all statistical analyses was 387.

Model Development and Validation

Several linear and non-linear methods were used for the QSAR model development. The compounds were divided into an external validation set and a training data. Models were developed using training set compounds and assessed using internal and external validation sets. To divide the compounds, they were ordered according to BE%, and from every set of five compounds, four were allocated into the training and one into the external validation set randomly. In this way, training data consisted of 168 compounds and the external validation set consisted of 40 compounds. For the analytical methods that required parameter optimization, a fraction of training set compounds were randomly assigned into internal validation set, or alternatively, cross validation was used if the option was available in the statistical software. STATISTICA Data Miner was the software used for statistical analysis. For the internal validation set, where applicable, the risk estimate and standard error were calculated in STATISTICA software and used as the performance indicators. Risk estimate is calculated as the proportion of residual variance incorrectly estimated by the model. Standard error measures the error of the prediction.

Mean absolute error was used to assess the accuracy of prediction of biliary excretion using the selected models (Eq. 5).

5

In Eq. 5, ‘observed’ refers to the log percentage of intact dose excreted into the bile from in vivo studies (LogBE%), ‘predicted’ is the predicted LogBE% by the QSAR models and N is the number of compounds. In all statistical analyses, the logarithm of percentage dose excreted (LogBE%) was used in the analysis instead of BE%. This was due to the normal distribution of LogBE% as indicated by the skewness comparison with BE%.

Regression Models

Stepwise regression analysis using in vivo rat biliary excretion data as the dependant variable resulted in the MLR model below in which the number of molecular descriptors is limited to eight. The statistical terms of the equation are N the number of compounds, R² the correlation coefficient, S the standard deviation and F Fisher’s statistics and the P value. Observed versus calculated LogBE% by this equation has been plotted (Figure S2 in supplementary material II), with training and validation sets identified in the plot.

MLR model

6

Molecular descriptors of this equation are not intercorrelated (R² < 0.4). Table I gives a brief description of molecular descriptors used in this model. Vsurf_HB4 is the first molecular descriptor selected by the analysis and it indicates that compounds with high H-bond donor capacity have higher biliary excretion level. AM1_dipole (dipole moment) is the other polarity descriptor which has a positive effect. On the other hand, the equation shows that drugs with greater relative positive partial charge (PEOE_RPC+) have lower biliary excretion. The value of this descriptor is large for small acidic molecules such as benzoic acid and salicylic acid, and therefore the small size of such compounds may be the reason for the reduced biliary excretion. In this equation, fU with a negative coefficient indicates that compounds with higher unionised fraction at pH 7.4 have lower biliary excretion. In other words, although according to fU, acidity and basicity (dissociation in general) increase the biliary excretion of compounds, this is true only for large dissociated molecules. The positive effect of polarity and dissociation on biliary excretion is in agreement with the literature, where for example polar surface area (18) and an acidity indicator (17,23) have been included in linear QSAR models. Also according to this equation, compounds containing many methyl groups (SsCH3) and those that are highly branched containing >CH− groups (SsssCH) have lower biliary excretion. In this equation, SlogP_VSA0 shows the negative impact of the presence of atoms with LogP(o/w) contribution of less than or equal to −0.4. SddssS_acnt indicates the direct effect of sulphate or sulphonamide groups. Sulphate and sulphonamide groups are found in sulphonamide drugs such as succinylsulphathiazole, dyes such as methyl orange and sulphate conjugates such as estrone 3-sulphate which may be substrates of MRP2 or BCRP (24).

Regression Tree Models Using CART

Several RTs were generated using a combination of molecular descriptors while cross-validation was applied (Table II). As seen in Table II, in RT (1), molecular descriptors were selected by C&RT analysis, while in RT (2), the molecular weight and in RT (3), the number of carboxylic acid groups were manually imposed as the first split descriptor using interactive C&RT routine in STATISTICA. These models were developed using the training set while the validation set remained external. The RTs resulting from these trials have been presented in Figs. 2, 3 and 4. In the regression trees, N is the number of compounds, Mu is the average and Var is the variance of LogBE% in each node. The molecular descriptors employed in the trees have been explained in Table I. Table III provides the statistical parameters of the regression trees. Observed versus calculated LogBE% plots are in supplementary material II (Figures S3 and S5) with training and validation sets identified in the plots.

According to RT (1), biliary excretion is much higher for compounds with large hydrophilic volume (vsurf_W3), especially if they are ionised with fU ≤ 0.001 (negligible unionised fractions at pH 7.4). Within the hydrophilic drugs of higher fU values (node 7), those with higher separation of lipophilic interaction sites from the centre of mass (vsurf_ID7 > 0.760) have higher biliary excretion. Surfactant molecules and glucuronide conjugates are examples of such molecules with high VolSurf integy moment (vsurf_ID7) and high biliary excretion. This branch follows to partition the molecules further according to GCUT_SLOGP_1 with compounds of lower hydrophobicity (node 18), and large hydrophobic interaction energy minima (vsurf_EDmin3 > −2.60) showing high biliary excretion (node 23). According to RT (1), the less hydrophilic drugs with vsurf_W3 values below 417.56 can be excreted heavily through the bile if they are highly dipolar (AM1-dipole >4.336) with high ratio of lipophilic to total surface area (FASA_H >0.50), especially if they are predominantly in the ionised form at pH 7.4 (fU ≤ 0.052). On the other hand, compounds with low dipole moment have low biliary excretion specially if they are lipophilic with LogD(6.5) >2.51 (node 9) or otherwise if they contain a high ratio of nitrogen atoms in the molecular structure (node 15). N ratio is low for larger alkaloids such as morphine or non-basic compounds, such as estrone 3-sulphate, which will have moderate biliary excretion especially if they are hydrophilic (PEOE_VSA-0 ≤ 94.24).

RT (2) was a result of molecular weight being employed in the first split using the interactive C&RT analysis in STATISTICA (Fig. 3). The statistically selected molecular weight threshold was 347.9 Da, with the compounds below this weight showing lower LogBE% values than the larger compounds. The tree shows that large (MW > 347.9) hydrophilic compounds (vsurf_CW4 >0.540) have higher biliary excretion, particularly those with large total negative van der Waals surface area (PEOE_VSA_NEG) and low surface area corresponding to highly polarisable groups (SMR_VSA7), especially if they are highly branched (SssssC > −1.812). Within this group of compounds, larger molecules with KierA1 >21.135 will have even higher biliary excretion. Other parameters of RT (2) indicate that high hydrophilic integy moment (vsurf_IW2) (node 13) and fractional negative charge weighted surface area (FCASA−) (node 11) would result in high LogBE% value.

Recent studies by Yang et al. (19) show that presence of carboxylic acid group(s) may indicate a trend towards increased biliary excretion. Therefore, the impact of presence of carboxylic acid group was examined using the interactive C&RT analysis with COOH used as the first partitioning molecular descriptor (Fig. 4). According to RT (3), compounds containing at least one carboxylic acid group have higher biliary excretion levels. Furthermore, RT (3) indicated that compounds with lower total negative partial charge (PEOE_PC−) have much higher biliary excretion (node 6). These are large hydrophilic compounds with many negatively charged atoms. Non-acidic compounds in node 2 will have high biliary excretion if the negative-charge weighted surface area for these molecules is high (node 5). CASA− has an element of size as well as indicating the presence of negatively charged groups such as sulphates.

Boosted Trees

Boosted tree module computes a sequence of simple trees, where each successive tree is built for the prediction of the residuals of the preceding trees. The analysis using various combination of model parameters resulted in two best models selected based on the error level for the internal test set (Table III). In models BT (1) and BT (2), the optimal numbers of trees were 145 and 147, with the learning rate of 0.10 and subsample proportions of 0.55 and 0.60, respectively.

It is possible to elucidate the influential descriptors in boosted trees analysis using variable importance calculation. Variable importance in STATISTICA is calculated as the relative (scaled) average value of the predictor statistic over all trees and nodes; hence these values reflect on the strength of the relationship between the predictors and the dependent variable of interest, over the successive boosting steps (STATISTICA help file, 2009). Included in Table I are the top ten most important molecular descriptors of BT (1) and BT (2) models. Some of the descriptors used by BT models are those already observed in RT and MLR models. For example, LogD (6.5) is present in two RT models and it is amongst the top ten most significant descriptors of both BT models. Other descriptors selected by these models are topological/size descriptors (KierA3, Kier2, Kier3) and other lipophilicity descriptors such as LogD at different pH values and vsurf descriptors. Table III shows the statistical significance of these models. Graphs of average squared error against number of trees for training and cross-validated test sets (Figures S6 and S7 for BT (1) and BT (2)) and plots of observed versus calculated LogBE% using BT (1) and BT (2) (Figures S8 and S9) can be found in supplementary material II.

Random Forest

In RF, the number of trees specifies the number of simple regression trees to be computed in successive forest building steps. The model development used the default values of the software with the number of trees set at 100. The graph of average squared error against number of trees for training and cross-validated test sets indicates that the test and training set errors reach a plateau at around 10–15 trees (see supplementary material II, Figure S10). The best model was achieved with a subsample proportion of 0.60, random test data proportion of 0.2 and number of trees of 100. Table I includes a description of the ten most significant descriptors employed in this model. Table III gives a summary of the statistical parameters of the RF model. Observed versus calculated LogBE% using the RF model can be found in supplementary material II (Figures S11).

Validation of the Models

All models were validated by the same external validation set which had been set aside and not used at any stage of model development. Table IV shows the prediction accuracy of the QSAR models using external validation in terms of the mean absolute error and the number of outliers. In addition, an average estimate of log BE% using all regression trees (RT (1)–RT (3)) was calculated and compared with the observed values to investigate any possible improvements in prediction accuracy. Table IV gives the performance of this estimation method (consensus RTs).

Comparison of the Models

Linear regression equation is one of the simplest and the most common QSAR techniques. This method has the benefit of easy interpretation and it can provide mechanistic insight into the process under investigation. However, it has been argued that many biological processes have more complex relationships with the molecular attributes of the compounds and hence linear regression models may fail to capture these (25). RT offers a suitable alternative to MLR method with the advantage of being flexibly non-linear while retaining the interpretability (26). Ensemble methods such as random forest (22) provide consensus predictions which may have improved accuracy. But this is often accompanied by a loss of interpretability, as the ensemble of many models is often used as a ‘black box’ prediction tool. In this investigation, STATISTICA variable importance analysis was used to find the most significant molecular descriptors in the boosted trees and random forest models.

According to Table IV, the most predictive model with the lowest estimation error for the external validation set is RT (1) followed by BT (1) and BT (2) and then RT (2). In other words, increasing the complexity of the models by allowing non-linear relationships and an ensemble of such models has been able to improve the prediction accuracy in comparison with a simple linear regression model (MLR). Table IV shows the number of outliers from each of the models. According to this table, RT (1) followed by BT (2) and BT (1) and then RT (2) are the best externally validated models with the lowest numbers of outliers in the validation set. The advantage of RT is the obvious simplicity and interpretability which can make it more popular with the end users in drug discovery disciplines. For example, when using the tree for a new compound, the molecular descriptors used in the tree will need to be calculated for the compound and then the terminal node (leaf) where the compound falls according to the molecular descriptor values should be identified. The average LogBE% of the terminal node (Mu) is the estimate of the tree for this compound. Despite that RT provides discrete predictions of a continuous observation which is not ideal, this is a much more straightforward procedure than using BT or RF for the estimation of BE%. These models are ensemble of many trees, and therefore the prediction has to be performed by the computer rather than manually.

An interesting observation was made as MW and COOH were not significant in MLR equation when forced into stepwise regression analysis (P > 0.05). Despite this, incorporation of these two parameters was statistically significant in CART analysis resulting in RT (2) and RT (3). This indicates the non-linear nature of the impact of these two parameters on biliary excretion. Average prediction by the three RT models was also considered and found to be of similar accuracy to RT (1) (Table IV).

In this work, the MLR model based on the training set of 168 compounds had the second poorest prediction accuracy after RF. Studies by Yang et al. (19) and Chen et al. (23) report MLR models based on training sets of 37 and 46 compounds, respectively. The proposed model by Yang et al. incorporated molecular connectivity indexes and atom-type electrotopological indexes which have also been used in this study. The model proposed by Chen et al. also incorporated similar molecular descriptors to our study, with the addition of Abraham descriptors representing polarisability and hydrogen bond acceptor capacity. Although we have not used Abraham’s descriptors, there are other molecular descriptors in our set of 386 descriptors that measure the same properties. Examples are the number of hydrogen bonding acceptor atoms and atomic charge on the most negatively charged atom in the molecule which may represent hydrogen bond acceptor ability (27) and molar refractivity descriptors which may indicate molecular polarisability (28).

In another study, Luo et al. (17) used 50 proprietary compounds from Bristol-Myers Squibb Co. for model development. They also developed a multiple linear regression model, but in addition to more common molecular descriptors, they employed free energy of aqueous solvation calculated from a self-consistent reaction field method. In analysing this model, Gandhi and Morris (18) found that the model failed to generalise further to the new set of compounds and specifically free energy of aqueous solvation was not statistically significant. They argued that a complex process such as hepatobiliary excretion cannot be captured by simple physicochemical properties when examining chemically dissimilar compounds. Indeed, such extrapolations to external compounds will fail when the compounds are outside the domain of applicability of the QSAR models. Incorporation of a larger dataset in this work may provide the opportunity for capturing an extended chemical space. This will be discussed further when analysing the outliers in the next two sections.

Structural Features of Compounds for Biliary Excretion

Table I gives a brief description of the significant molecular descriptors used in the models. For the sake of this discussion, the descriptors in this work can be classified roughly into five categories as follows: lipophilicity, ionisation, molecular size and topological and constitutional descriptors.

It can be seen in Table I that lipophilicity descriptors such as log D at different pH levels and surface area of hydrophilic molecules (SlogP_VSA0) are present in all models. In all interpretable models (except for the linear regression equation), lipophilicity descriptors show a negative effect on the biliary excretion of compounds. This may relate to the fact that highly lipophilic compounds are known to be highly extracted and metabolised in the liver (29) rather than being excreted unchanged through bile or kidney. For example, metabolism by cytochrome P450 enzymes (30) and (UDP)-glucuronosyltransferase (31) is mainly controlled by lipophilicity and increased for more lipophilic compounds. There have been inconsistent findings in the literature regarding the effect of lipophilicity on the biliary excretion of xenobiotics. Proost et al. found no significant correlation between lipophilicity and biliary excretion of a series of bulky organic cations despite it being the predominant factor for the degree of plasma protein binding and hepatic uptake rate (29). Similar observations have been made for other compilations of biliary excretion data (19). Other studies indicate negative effect of lipophilicity on the biliary excretion within the range of compounds studied (17,32). Lipophilicity has been associated with many models of ADME properties (33). It is a well-established fact that compounds with higher logP have poor aqueous solubility and are more likely to pass through lipid bilayer of biological membranes (34). The general trend in the literature with regards to the role of lipophilicity in pharmacokinetic processes indicates that more lipophilic compounds have higher oral absorption, plasma protein binding, and volume of distribution (35–37) and are more prone to P450 metabolism (30,35). This may lead to the reduced chance of excretion through bile as the intact drug.

All models presented in this work indicate the significant role of ionisation and polarity through molecular descriptors such as COOH, fU, FCASA− and SddssS_acnt. Acids are able to ionise into anions which are substrates of several transporters (generally organic anion transporters). Compounds that carry positive as well as negative charge or partial charges can use both the ‘organic anion’ and the ‘organic cation’ transport systems (38). For example, OAT3 accepts various kinds of bulky hydrophobic anions, while OAT1 can transport relatively hydrophilic small molecules, such as nucleoside analogues (39). Besides, monocarboxylate transporters (MCT1 to MCT14) constitute a family of proton-linked plasma membrane transporters that carry molecules having one carboxylate group. MCT1 is expressed nearly all over in every tissue in the human body and also in rat and calves hepatocytes (40). MCT2 is abundant on the surface of human, rat and hamster hepatocytes (41). MCT5 and MCT8 are also known to play transporting role in rat hepatocytes (41). Studies of biliary excretion of exogenous compounds have indicated the relation between polarity and biliary excretion stating that possession of a strongly polar anionic group was important factor in appreciable biliary excretion (17,42). In all the interpretable models reported here, polarity descriptors show a positive impact on biliary excretion. Examples are the positive coefficients of dipole moment (AM1_dipole) in the linear regression equation and higher percent of compounds with lower unionised fractions at pH 7.4 (fU) in RT (1).

Molecular size is the other important factor in biliary excretion represented in the models by molecular descriptors such as kappa shape indexes, hydrophobic volumes (vsurf_W1 and vsurf_W3) and surface areas of atoms with specific charge or lipophilicity ranges (e.g. PEOE_VSA_NEG and PEOE_VSA_HYD). These molecular descriptors show positive effect on biliary excretion level in all models. This is in line with the common understanding that a molecular weight threshold may apply to biliary excretion of compounds, and that high molecular weight compounds may be predominantly excreted through bile (19,32,42). Yang et al. (19) suggested a molecular weight threshold value of 400 Da for biliary excretion of anionic drugs in rats using 164 drugs. In this study, regression tree analysis found the threshold value for molecular weight to be at 347.9 Da for biliary excretion in rat (RT (2)). Incidentally, this regression tree had the second highest prediction accuracy for the external validation set amongst the RT models. This was despite the fact that molecular weight was not the descriptor of choice by C&RT analysis.

The incorporation of some structural fragments in the models gave interesting information regarding molecular requirements for biliary excretion. Examples include SddssS_acnt and SsssCH which indicate higher biliary excretion of compounds containing sulphate groups and non-branched structure (MLR). Compounds containing carboxylic acid groups are also more likely candidates for biliary excretion according to RT (3). Up to half of compounds in our dataset contain –COOH groups (103 compounds out of 217). Sixty-five out of 103 COOH containing compounds had biliary excretion of >20%. Varma et al. (32) have analysed the interconnection between physicochemical requirements of OATP substrates and the biliary excretion rates. It was then suggested that substrate specificity of OATPs including acidity may primarily indicate the elimination through bile (32).

Analysis of the Outliers

There are a number of compounds that are outliers from majority of the models. Analysis of outliers may provide interesting information regarding the applicability of the models. Within the BE% range, it could be observed that compounds with low biliary excretion show a higher average error in general (Table V). For example, the average error by all seven models was the highest for the six compounds with the extremely low biliary excretion (BE% <0.23), followed by the compounds with 0.23 < BE% < 1.23 (−0.64 < LogBE% < 0.09). A closer inspection of the data reveals that despite the high average error for the six compounds with low biliary excretion, the estimation may still be acceptable as all these compounds have been estimated to have a BE% value <4% (average of all models) and below 0.6% by RT1 model with only one exception (benzoic acid). A hypothesis here could be that these compounds may have suitable properties for higher biliary excretion, but other routes of elimination are predominating. For example, it has been shown for benzoic acid that when clearance by the kidney is prevented, biliary excretion increases by 10% (12). Out of 217 compounds in the dataset, the predominant routes of elimination are biliary excretion for 115 compounds, renal excretion for 65 compounds and metabolism for 37 compounds. However, the outlier compounds do not belong to any single groups above in terms of the predominant routes of elimination.

According to Table V, highly lipophilic compounds (log P > 5.35) and low molecular weight compounds (MW ≤ 280) also show higher error rates, and this may need to be considered when using the models for the prediction of external compounds.

Table VI gives a list of the compounds that are outliers in six or seven models out of the seven models proposed here. In addition, there are four compounds which were outliers in four or five models but had exceptionally high average error from the seven models. These compounds were part of the training or validation sets, but none were omitted from average error calculations.

The outliers in Table VI have been over- or under-predicted by the models. One compound in the table has shown underestimation by some and overestimation by other models; biliary excretion of benzoic acid was overestimated by all models except for BT (1) and BT (2). It can be seen in Table VI that fosmidomycin, nelfinavir and 2-ethyl-5-methyl-3,3-diphenyl-1-pyrroline (EMDP) are over-predicted by five or more models. Benzoic acid is rapidly cleared by the kidney, so it may not have enough time to pass into the bile (12). Abou-El-Makarem and his colleagues examined this possibility by tying up the renal pedicles in rats, so that clearance by the kidney was prevented, and the results indicated that when clearance by the kidney is prevented, biliary excretion increased by 10% (12). Fosmidomycin has a short half-life of 1.7 h and is rapidly cleared by the kidneys (43). It is a small molecular weight polar agent which may not be cleared in high quantities through bile according to the molecular weight threshold hypothesis. Despite the use of molecular size descriptors, this compound still appeared to be overestimated by all seven models, even using RT (2) which has employed MW for the first branching. The problem with RT (2) in relation to this compound is that although this compound falls into node 2 along with 44 other low molecular weight compounds, this node has an average LogBE% of 0.42 which is much lower than node 3 with an average LogBE% of 1.27 but not low enough for this compound. Likewise, other models have indicated low biliary excretion of small-sized compounds, but somehow, estimation is higher than what is actually observed.

Nelfinavir has a half-life of 3.5 to 5 h and is eliminated via metabolism by the cytochrome P450 enzyme system (44). This is a highly lipophilic compound which is poorly excreted through bile, and is predicted as such by the models (predicted BE% below 2% using all models except for RT (3) and RF which predict 13 and 7.6%, respectively).

EMDP is a major metabolite of methadone which has been over-predicted by most models despite a very low biliary excretion. As with nelfinavir, the predicted BE% for this compound by most models is quite low at <4% (MLR is an exception) and the selected model, RT (1), predicts a biliary excretion value of ∼0.3%. Despite this, in comparison with the extremely low observed value of 0.05%, the predicted values are much higher, leading to a numerically large average error, even though qualitatively, the predicted biliary excretion may be reasonably low.

EDDP, PAEB (procaine amide ethobromide) and tolrestat are the under-predicted compounds. All these compounds have high BE% values at 36, 32 and 54%. This is despite the relatively low molecular weights of EDDP and PAEB which are below the defined MW threshold of 347 Da for biliary excretion. The exact mechanism of high biliary excretion of these compounds warrants further investigation to explore the reasons behind such high biliary excretion despite the low molecular weight.

Tolrestat has a relatively high molecular weight suitable for biliary excretion and a COOH group making it a suitable substrate for OATPs (32). Despite this, the hydrophilic volume calculated by the VolSurf descriptor vsurf-W3 is not high enough to put this compound in node 3 rather than node 2 of RT (1) model. In RT (2), the compound falls into node 16, which is due to the lack of non-aromatic branched structure which would place it in node 17 with a higher predicted BE%. Likewise, in RT (3), this compound fails to be placed in node 7 and falls in node 6 instead due to the low total negative charge (> −2.33) as a result of the low number of negatively charged atoms. This indicates a shortcoming in the abovementioned models which lack suitable parameters that can capture the relative polarity in relation to the molecular size.