|Home | About | Journals | Submit | Contact Us | Français|
The aims were (1) to evaluate the molecular weight (MW) dependence of biliary excretion and (2) to develop quantitative structure–pharmacokinetic relationships (QSPKR) to predict biliary clearance (CLb) and percentage of administered dose excreted in bile as parent drug (PDb) in rats and humans. CLb and PDb data were collected from the literature for rats and humans. Receiver operating characteristic curve analysis was utilized to determine whether a MW threshold exists for PDb. Stepwise multiple linear regression (MLR) was used to derive QSPKR models. The predictive performance of the models was evaluated by internal validation using the leave-one-out method and external test groups. A MW threshold of 400 Da was determined for PDb for anions in rats, while 475 Da was the cutoff for anions in humans. MW thresholds were not present for cations or cations/neutral compounds in either rats or humans. The QSPKR model for human CLb showed a significant correlation (R2=0.819) with good prediction performance (Q2=0.722). The model was further assessed using a test group, yielding a geometric mean fold-error of 2.68. QSPKR models with significant correlation and good predictability were also developed for CLb in rats and PDb data for anions or cation/neutral compounds in rats and humans. Both CLb and PDb data were further evaluated for subsets of MRP2 or P-glycoprotein substrates, and significant relationships were derived. QSPKR models were successfully developed for biliary excretion of non-congeneric compounds in rats and humans, providing a quantitative prediction of biliary clearance of compounds.
Biliary excretion is an important route for the elimination of drugs and/or their metabolites (1). However, it is difficult to experimentally measure the biliary excretion of drugs in humans (2) due to the anatomy of the hepatobiliary system and the incomplete and intermittent elimination of bile into the intestinal tract. Most biliary excretion studies in humans have been performed by collecting bile samples in patients requiring bile duct cannulation, duodenal fluid collection in healthy volunteers, or analyzing fecal samples. There is a need to accurately estimate the biliary excretion of drugs in humans early in drug development. Such information would be valuable for determining the importance of biliary excretion in total clearance, predicting potential drug–drug interactions and assessing the risk of toxicity.
To undergo biliary excretion, drug present in the portal venous or hepatic portal blood must first cross the sinusoidal membrane and enter hepatocytes. Inside the hepatocytes, drug can be metabolized, bind nonspecifically, efflux back into blood across the sinusoidal membrane or cross the canalicular membrane to be excreted into bile. Drug in bile is usually more concentrated than in plasma, and passive transport is thought to be negligible. As a result, biliary excretion represents an active process involving transporters embedded in the canalicular membrane of hepatocytes. These transporters belong to the ATP binding cassette (ABC) superfamily of transporters. The major transporters involved in the biliary excretion of drugs include ABCB1 (P-glycoprotein, MDR1, P-gp), ABCC2 (multidrug resistance-associated protein 2 (MRP2)), and ABCG2 (breast cancer resistance protein (BCRP)). ABCB4 (MDR3) and ABCB11 (bile salt export pump (BSEP)) are two other ABC transporters present at the canalicular membrane but they play a minor role in transporting drugs.
The influence of factors including molecular weight (MW), lipophilicity, polarity, and molecular structure of the drug, as well as species and gender differences, on biliary excretion in animals and humans are not well understood (3). These determinants are likely important for binding to transport proteins on the hepatocyte basolateral or canalicular membrane. A MW cutoff for biliary excretion was originally reported in the 1970s, with a MW threshold of 325±50 for rats and 500±50 for humans (4,5). Based on the above threshold, it is expected that there is a higher chance for a compound with MW over 325 to undergo substantial biliary excretion (defined as >10% of the dose) in rats, while a higher MW is necessary in humans. Importantly, the above thresholds were reported for anions in both rats and humans (4,5). In rats, a MW threshold was reported to be 200±50 for monoquaternary ammonium and 500–600 for diquaternary ammonium compounds (3,6).
Some general comments are warranted regarding the derivation of these MW threshold values for biliary excretion. (1) The MW threshold values determined in rats were mainly based on several series of compounds and their derivatives or metabolites, including benzoic acids, sulfonamides, derivatives of biphenyl, stilbestrol or phenolphthalein, and some dyes, many of which are not drug-like chemicals (5,7,8). It is not known if the threshold will remain the same for drugs and for compounds with broader structural diversity. (2) In cases where the administered compound was mainly present in bile as a polar conjugate (>70–90% of the biliary material), the MW of the formed conjugate was used instead of that of the compound itself (7). In these cases, it is not known whether the administered (pre-formed) conjugate exhibits similar biliary excretion. (3) In most studies, the compounds were administered intraperitoneally rather than intravenously (5,7,8). Since bioavailability was not considered, the results may not be accurate. (4) The MW cutoff value suggested for humans was based on studies with a very limited number of compounds (4, 9–11). (5) There was no statistical basis provided for the derivation of the threshold values. Therefore, one of our goals was to evaluate the evidence supporting a MW cutoff in rats and humans based on a dataset of biliary excretion or biliary clearance data including a diverse range of drug molecules that was assembled from the literature.
In silico models have been explored and applied to predict properties of compounds for many years. One of the well-established methods is quantitative structure–activity relationships (QSAR) (12). The primary goal of QSAR is to construct a relationship between the pharmacological activities of compounds and their physicochemical properties that can be easily calculated using computer software. The relationship is usually described by a mathematical expression and can be used to predict the pharmacological activities of untested compounds. The concept can also be applied to pharmacokinetics by simply replacing the activities with pharmacokinetic parameters. Such quantitative structure–pharmacokinetic relationship (QSPKR) models have been successfully applied in several studies to predict clearance, volume of distribution, bioavailability and other pharmacokinetic parameters (13–21).
Our objectives for this study were: (1) to compile our own dataset of literature data providing accurate and reliable estimates of biliary clearance (CLb) and percentage of dose eliminated in the bile (PDb) in rats and humans; (2) to examine relationships between PDb and MW for the dataset, and for subsets of the data including anions and cations/neutral compounds; and (3) to develop QSPKR models to predict biliary excretion from data published for rats and humans. Since data in the literature are available in two forms, namely CLb and PDb, QSPKR models were constructed for both types of data.
Biliary clearance (CLb) and percentage of dose excreted into bile (PDb) data were collected from studies in rats. Compounds were administered intravenously (bolus or infusion) in all of the studies. Both CLb and PDb values for parent drug were extracted. Most of the biliary clearance values reported in the literature were calculated using one of the following methods:
Xbile is the amount of drug excreted into bile during a certain period, equal to the product of bile concentration of drug and bile volume. AUCplasma is the area under the curve of plasma concentration versus time during the same period. υbile is the rate of biliary excretion, which can be calculated by dividing Xbile with the corresponding time interval. Cplasma,ss is the steady state plasma concentration of drug. Cplasma,midpoint is the plasma concentration at the mid-time point of the corresponding collection period of bile.
In some studies, biliary clearance data were not directly reported but could be calculated using Eq. 4 based on the information available in those studies:
CLsys is the systemic plasma clearance of parent drug after intravenous administration and (fraction of dose)bile is the fraction of the administered dose excreted into bile as parent drug.
CLb and PDb values for parent drug in humans were either reported in the literature or calculated in the same way as for rats. The reports in which compounds were administered intravenously were preferably selected to extract the data. Biliary excretion studies in humans were usually conducted by collecting bile (T-tube in most cases) in patients undergoing certain types of surgery or duodenal fluid in healthy volunteers after administration of drugs. In general, the data from studies conducted in healthy subjects were chosen. Since there were a limited number of biliary excretion studies in humans following the intravenous administration of drugs, some studies with oral administration were also utilized. In these cases, PDb data were corrected by bioavailability (F). Biliary clearance could be calculated from fraction of the dose excreted in the bile and oral clearance (CLpo).
One consideration in using data after the oral administration of drugs is the possibility that first-pass extraction may be due to hepatic biliary elimination. For all drugs except one in this study, where data after oral administration were used, biliary elimination was low (<5%), or oral bioavailability high, indicating that first-pass extraction due to biliary elimination was negligible and should not influence the biliary excretion determination. There was only one compound in our dataset that may exhibit some first pass biliary excretion and that was rosuvastatin. The corrected PDb following oral administration of rosuvastatin was 59%; this value was consistent with the renal excretion of rosuvastatin being 28% and its metabolism being low in humans (22,23).
Finally, the units for CLb in rats were converted to milliliters per minute per kilogram using the body weight values provided in the corresponding studies. CLb values in humans were usually reported with the units of milliliter per minute. If the CLb values were not presented in this manner, they were adjusted to milliliter per minute using a body weight of 70 kg or body surface area of 1.73 m2.
Besides CLb and PDb data, information on administration route, collection period of samples, analytical method, MW of the compounds, charge status of the compounds under physiological conditions (pH=7.4), and transporters known to mediate efflux processes of the drugs were also collected whenever available. For the human studies, population studied (healthy volunteers or patients), number of subjects, and sample collection method (T-tube or duodenal perfusion) were also listed.
Overall, four separate datasets were created: PDb in rats (164 compounds), CLb in rats (55 compounds), PDb in humans (97 compounds), and CLb in humans (56 compounds). For the QSPKR analysis, CLb and PDb datasets in rats and humans were randomly divided into training sets and test sets. The numbers of the compounds in the training sets were two- to threefold greater than the number of compounds in the test sets.
We have used a statistical technique known as receiver operating characteristic (ROC) curve analysis as a diagnostic test to determine whether a MW cutoff for biliary elimination is present and, if present, what the MW cutoff value is in rats and humans. The ROC curve depends on the calculation of sensitivity which is defined as the ratio of positives selected by the test to the true positives and specificity which is defined as the ratio of negatives selected by the test to the true negatives (24). The AUC of a graph of sensitivity versus (1 - specificity) can be used as a diagnostic of the performance of the test. An area of 1 would mean that the test is 100% specific and sensitive, while an area of 0.5 or below would mean that the test arrived at results randomly. Although there is no “gold standard” for AUC, a value larger than or equal to 0.8 is usually regarded as acceptable (25–27).
In our case, true positive is defined as greater than or equal to 10% of administered dose excreted into bile as parent drug and true negatives as less than 10% of dose excreted in bile as parent drug. For each MW cutoff selected, a pair of sensitivity and specificity values can be calculated. ROC curves are then created by plotting sensitivity versus (1 - specificity) and an AUC and its standard error (SE) calculated (28). A MW threshold exists if the AUC is at least 0.8 and the relationship between PDb and MW is significant as defined by a Chi-square test (p<0.05 as criteria for significance). In theory, the point corresponding to the cutoff value should fall at the point on the ROC curve that has the shortest distance to the point (0,1) on the y axis, i.e., the point having both sensitivity and specificity of 1. Based on this distance, the candidates for MW cutoff were located, and the cutoff was accepted if both sensitivity and specificity were greater than or equal to 80%, meaning that false positive and false negative probabilities are less than 20%.
ROC curves were calculated for PDb and CLb in rats and humans to determine whether a relationship exists between MW and biliary excretion and if a relationship exists, what the MW cutoff could be. We also generated ROC curves to determine if there is a relationship between PDb or CLb and log P or log D. The calculation of log D takes into consideration the ionized species along with the unionized species. The log D at pH7.4 was calculated using Chemaxon software (Chemaxon Ltd, Budapest, Hungary). Molecular volume (MOVO) was calculated using the molecular modeling software SYBYL, as described below, and additional analyses examined the relationships between MW and MOVO. ROC curves were used to examine if a relationship between MOVO and biliary excretion exists.
The structures of compounds were obtained from SciFinder scholar database (American Chemical Society, Washington, DC) and schemed using ChemDraw Ultra 8.0 software (CambridgeSoft, Cambridge, MA), through which the SMILES string, a unique code representing each compound, was created. We initially derived QSPKR models from structural descriptors calculated using the molecular modeling software SYBYL. Compounds can be loaded into the software in their ionized form, an added advantage of SYBYL. We calculated 11 descriptors using SYBYL including molecular volume and Connolly surface area. clogP and MW were calculated using the online version of the KowWin program developed by Syracuse Research Corporation (http://www.syrres.com/esc/est_kowdemo.htm) and log D at pH of 7.4 was calculated using Chemaxon software. The QSPKR models were derived using multiple linear regression and partial least squares regression but did not satisfactorily predict the CLb and PDb for both rats and humans. We concluded that this limited set of molecular descriptors was inadequate to explain the complicated process of biliary elimination.
We then proceeded to calculate a total of 116 structural descriptors for each compound based on the SMILES string. cLog P (of unionized form) and MW values were calculated using the online version of KowWin program, and the other 114 structural descriptors were calculated using the QSARis© version 1.2 software (SciVision-Academic Press, San Diego, CA). These 2D and 3D structural descriptors include molecular connectivity Chi indices, Kappa-shape indices, electrotopological state indices, information indices, subgraph count indices, molecular polarizability, weight, and volume.
The selection of the independent variables (structural descriptors) in the QSPKR model was accomplished by stepwise MLR. First, forward selection was performed using QSARis© software to select the most significant descriptor (p<0.05), and this step was repeated until no further significant descriptor could be added. During this process, the number of the descriptors chosen was controlled to be no more than one fifth of the number of compounds to avoid overparameterization (29). Secondly, backward elimination was performed using SAS© statistical software (version 9.0, SAS Institute Inc., Cary, NC) to remove any redundant descriptors using a p value less than 0.05 as the criterion of significance. Lastly, correlation analysis was conducted to ensure no significant colinearity among the descriptors in the model. The final model was generated in the form of:
where Y is the predicted value, xi are the individual structural descriptors, ai are coefficients, and b is a constant.
The performance of the constructed model was then evaluated by cross-validation using the “leave-one-out” method, a commonly used internal validation method (18,19). In this method, the model is reconstructed over the dataset with one compound excluded and the target property of the removed compound is then predicted by the reconstructed model, resulting in a Q2 value.
The prediction performance of the model developed for CLb and PDb in rats and humans, based on the training set, was further assessed with the test set to provide external validation.
The geometric mean fold error was also calculated for CLb data since it gives equal weight to the under- and overestimated values and has been used to evaluate error for clearance data in QSPKR models (19).
A model that predicts the observed data perfectly or without any bias will have a geometric mean fold error of 1. A mean fold error ≤2 is generally considered successful, indicating that most of the predicted values fall within the twofold range of the observed ones (between 0.5 and 2.0) (19).
For the PDb data, the mean error and the root mean squared error (RMSE) were used for assessment of the bias and accuracy of predictions, respectively, as previously used to evaluate error for % bioavailability data (30).
Figure 1a shows a plot of PDb against MW for all 164 compounds. Though the relationship between PDb and MW is significant as defined by the Chi-square test (p<0.01), the area under the ROC curve for all compounds was less than 0.8 (0.761±0.037), suggesting that a MW threshold may not exist with 80% certainty (Table I). The point corresponding to a MW cutoff of 400, which had the shortest distance from the ROC curve to the point (0,1) on the y axis, had a specificity of 72% and sensitivity of 82%. We further evaluated the relationship between PDb and MW for subsets categorized based on charge status of the compound. A significant relationship was derived for anions (p<0.01) with a higher AUC (0.839±0.043) than the entire dataset, indicating that there may be a MW cutoff based on our set criteria (Fig. 1b). By comparing the sensitivity and specificity at a MW cutoff from 350 to 425, 400 Da was determined to be the threshold, while 375 Da almost met the requirements for sensitivity and specificity (Table I). The AUC calculated for the cation subset was below 0.8, and the specificity was only 65% at a MW cutoff of 475. This suggests that a MW threshold is absent for cations in our dataset. The situation is similar for the combined dataset of cations/neutral compounds. Due to the limited number of compounds in the neutral and zwitterion compound groups, a relationship between PDb and MW was not explored for these subsets. Overall, a MW threshold was observed for anions in rats (400 Da) but not for cations alone or cations and neutral compounds.
A relationship between PDb and MW for all 97 compounds was significant using a Chi-square test (p<0.01) and the AUC for the ROC curve was high (0.838±0.049; Fig. 1c). Examination of the distance from the curve to the point of (0,1) revealed a MW cutoff of 475 Da with a sensitivity of 80% and specificity of 79% (just below our requirement of 80%). The presence of a MW threshold was established for anions with an even higher AUC (0.939±0.045) and significant relationship (p<0.01). A MW cutoff of 475 Da met all the requirements and, thus, was determined to be the threshold for anions in humans (Fig. 1d). In contrast, there was no significant relationship between PDb and MW for cations (p=0.068) or cations and neutral compounds (p=0.073) with lower AUC values of 0.697±0.105 and 0.721±0.093, respectively. Both sensitivity and specificity were below 80% for a MW cutoff of 500 for cations and cations/neutral compounds. Hence, a MW threshold is unlikely to be present for cations or cations/neutral compounds in humans. A relationship between PDb and MW was not assessed for neutral or zwitterion compounds, since there was only a limited number of compounds available. In conclusion, a MW threshold of 475 for anions in humans was calculated, and this was higher than the one in rats. A MW threshold was absent for cations or cations/neutral compounds (Table II).
We also conducted ROC curve analysis to evaluate the relationship between CLb and MW in rats and humans using 3 or 10 ml/min/kg as a cutoff for CLb in rats and 3 or 10 ml/min for CLb in humans. These clearance values were chosen based on a review of CLb values in rats and humans in our datasets. No threshold was found in either species. We further assessed the relationship between PDb or CLb and LogP or Log D and did not find any statistically significant relationships (data not shown).
MOVO values were highly correlated with MW based on our rat dataset with a relationship of MW=1.06*MOVO+72.1 (R2=0.868, p<0.001). A ROC curve was generated for PDb and MOVO for the 82 anions in rats. The relationship was significant as indicated by a Chi square test (p<0.01), and AUC was high (0.839±0.048). However, the sensitivity (74%) and specificity (76%) were less than 80% at a MOVO cutoff of 300 Å which had the shortest distance to the point of (0,1) on the y axis of the ROC curve. Therefore, a significant MOVO threshold could not be defined for this dataset. It was noted that a MOVO of 300 Å corresponded to about 390 Da MW based on our correlation between MW and MOVO. This value is very close to the MW threshold determined for anions in rats. No other relationships for MOVO and biliary excretion (total compounds, cations, or cations/neutral compounds) were significant.
There was also a significant correlation between MOVO values and MW in our human dataset, and the relationship was very similar with that in rats (MW=1.07*MOVO+70.7, R2=0.806, p<0.001). A ROC curve was generated for MOVO and PDb for all the 97 compounds in the human dataset, and the relationship was found to be significant with an AUC of 0.837±0.050 (p<0.01, Chi square test). A MOVO cutoff of 375 Å resulted in a sensitivity of 83% and specificity of 79% and almost met our criteria; these results were similar to our results for the analysis for MW threshold. A significant relationship was established for the anion subset in humans with a higher AUC of 0.938±0.045 (p<0.01, Chi square test). A sensitivity of 87% and specificity of 93% were obtained at a MOVO cutoff of 375 Å. In addition, a MOVO of 350 Å also met the requirements for a threshold value, although the specificity (86%) was lower than that at 375 Å. In conclusion, a MOVO cutoff of 350–375 Å for anions in humans was determined. Based on our dataset for humans, a MOVO of 375 Å corresponded to a MW of about 475 Da, which was the MW threshold derived for the anion subset in humans. Significant molecular volume thresholds for cations or cations/neutral compounds were not found.
CLb data were collected for 55 compounds in rats, ranging from 0.012 ml/min/kg (thyroxine) to 37.3 ml/min/kg (BQ-485). Lipophilicity reflected by cLogP varied from −2.83 (ouabain) to 7.28 (indocyanine green), with a median value of 1.85. Univariate analysis showed no significant relationship between CLb and cLogP, indicating the necessity of a more sophisticated analysis. Thus, stepwise MLR analysis based on 116 physicochemical descriptors was conducted to construct a QSPKR model. Before developing the model, the dataset was randomly separated into a training set (37 compounds) and a test set (18 compounds). There was no difference between the two groups in terms of distribution patterns of CLb, lipophilicity, and molecular weight. CLb of the test group varied from 1.41 ml/min/kg (cefixime) to 28.2 ml/min/kg (BQ-123) with cLogP values covering a range from −2.74 (cefotetan, cephamycin) to 6.97 (Merck compound A (31)).
A simple QSPKR model with four structural descriptors was constructed for the training group as follows:
The selected descriptors included the indices of molecular moment (Q), electrotopological state (Hmin), molecular connectivity (xch9), and subgraph count (SsssN_acnt). The predicted CLb values correlated fairly well with the observed values with R2 of 0.689 (Fig. 2). Prediction performance of the model was assessed by cross-validation producing a Q2 of 0.618, suggesting that the model could be used for prediction purposes. Further examination of the training set revealed that 19 out of 37 compounds were MRP2 substrates and 13 compounds were transported by P-gp, based on in vitro or in vivo studies in the literature. Thus, QSPKR models were also developed for CLb of these two subgroups from the training set, MRP2 substrates and P-gp substrates (Fig. 3):
Since there were few known substrates of the ABC transporter BCRP in the training set, a QSPKR model was not developed for that subgroup.
The prediction performance of the general model constructed based on the whole training set was further evaluated with the test group consisting of 18 compounds. The predicted values of nine compounds fell within the twofold error range of the actual values with the prediction of another four compounds close to this zone (within threefold error range of the observed values), resulting in a R2 of 0.262 (Fig. 4a). The geometric mean fold error was calculated as 2.98. Interestingly, some outliers were substrates of MRP2 (CP-671,305) or P-gp (ulifloxacin (UFX)) or both (grepafloxacin). Therefore, the original predictions for the substrates of MRP2 or P-gp in the test group were then replaced with the predicted values derived from the MRP2 or P-gp substrate submodels that are provided above (Table III). In the case where a compound was a dual substrate of MRP2 and P-gp, the final prediction was the average of the values from the two submodels. The geometric mean fold error was further decreased to 1.87 in this way (Table III, Fig. 4b). The correlation between predicted and observed CLb was also improved with an R2 of 0.401.
CLb data were collected for 56 compounds in humans, ranging from 0.16 ml/min (ampicillin) to 400 ml/min (dibromosulfophthalein). Molecular weight varied from 151.2 (paracetamol, acetaminophen) to 1,202 (cyclosporine A). cLogP values covered a range from −2.74 (cefotetan, cephamycin) to 9.5 (metocurine) with a median of 1.98. Similarly with the case in rats, no relationship was found between CLb and cLogP. Since the range of CLb values in humans is quite large, the original data were first transformed by taking the natural logarithm. A QSPKR model was then developed for 38 compounds in the training group as following:
The model was based on seven structural descriptors related with molecular properties (nelem), electrotopological state (SddssS, SHsNH2, and SsCH3), molecular moment (Dz and Qyy), and general 3D property (ovality) of the molecules. Predicted natural log-transformed values from the model correlated very well with the observed LnCLb with a R2 of 0.819 and evenly distributed around the unity line (Fig. 5a). Internal validation (“leave-one-out” method) was performed producing a Q2 of 0.722. The model was further applied to a test set with 18 compounds. Prediction for six out of 18 compounds fell within twofold error range of the observed CLb values, with another seven compounds predicted within a two- to threefold error range (Fig. 5b). Geometric mean fold error was calculated as 2.68, comparable to that for CLb in rats. Overall, these results suggested that the model could be used to predict human CLb of unknown compounds.
In addition, QSPKR models were also successfully developed for 14 MRP2 substrates and 14 P-gp substrates present in the 38 compounds in the dataset (Fig. 6):
Among the 18 compounds in the test set, vindesine has been shown to be a P-gp and MRP2 substrate. The general model shown above underpredicted its CLb (3.62 ml/min) compared with the experimental value (29 ml/min). With the adjustment from both P-gp and MRP2 submodels, the prediction for vindesine was 28.6 ml/min, very close to the observed value. The geometric mean fold error for the test set was decreased to 2.39 with this adjustment. Other compounds in the test set have not been reported to be MRP2 or P-gp substrates, so no further corrections were made.
PDb data were collected for 164 compounds in rats, ranging from 0% (butoprozine, diazepam, felodipine) to 99.7% (J-104132). Molecular weight of these compounds varied from 117 (triethylmethylammonium, TEMA) to 1,255 (actinomycin D). cLogP values covered a wide range from −2.83 (ouabain) to 16.5 (cosalane) with a median value of 2.16. No significant relationship was found between PDb and cLogP or molecular weight. A QSPKR model with a R2 of 0.497 between the predicted and observed values was constructed. However, the Q2 value produced from the “leave-one-out” cross validation was less than 0.5 (Q2=0.437), suggesting that the model was unstable and could not be used for prediction purposes.
Further analysis of the dataset revealed that these compounds could be further subgrouped as 82 anions, 47 cations, 23 neutral compounds, and 12 zwitterions. This classification was based on calculated pKa values from the SciFinder Scholar Database, with 10% ionized at physiological pH (7.4) used as a cutoff. In our datasets, we had a greater range of values of PDb values for anions than we had for cations or neutral compounds. PDb for the 82 anions ranged from 0% (prostacyclin and QMPB) to 99.7% (J-104132), with about a third of the compounds having greater than 50% of the dose excreted in the bile. On the other hand, the values for the 70 cation/neutral compounds varied from 0% (butoprozine, diazepam, felodipine) to 58.3% (ouabain), with only three compounds with biliary excretion equal to or larger than 50%. In addition, it was found that among these compounds, there were 39 MRP2 substrates and 35 P-gp substrates, based on in vitro or in vivo studies in the literature. Almost two thirds of the MRP2 substrates (25 compounds) were anions. In contrast, almost all P-gp substrates (33 compounds) were cations or neutral compounds. Therefore, we evaluated QSPKR models for subsets of anions and cations/neutral compounds separately and also for subsets of MRP2 and P-gp substrates.
A model consisting of five structural descriptors was constructed for the training set of anions (61 compounds):
The PDb predicted by the model correlated very well with the observed values with R2 equal to 0.784 (Fig. 7a). The ‘‘leave-one-out’’ cross-validation generated a Q2 of 0.737, indicating that the constructed model performed stably and could be used for prediction purposes. Similarly, a QSPKR model with R2 of 0.828 and Q2 of 0.731 was established for the training set of cations and neutral compounds (46 compounds; Fig. 7b):
The structural descriptors in the equation represented different categories of molecular properties of compounds, such as hydrogen bond interaction (SHBint9_Acnt), molecular connectivity (xvch5 and xvch7), charge status (MaxNeg), and electrotopological state (SsOH). The predictability of the models was further evaluated with two test sets, i.e., 21 anions and 24 cations or neutral compounds. RMSE was calculated as 11.2 for cation/neutral compounds and 14.0 for anions, and mean error was −3.37 for anions and 4.39 for cations and neutral compounds.
QSPKR models were also successfully developed for subsets of 36 MRP2 and 28 P-gp substrates as following:
There were three MRP2 substrates in the test sets of anions (CP-671,305 and ceftriaxone) and cation/neutral compounds (belotecan). The prediction was greatly improved by the MRP2 submodel compared with the general model. Improvement of prediction using the P-gp submodel was also observed. There were seven P-gp substrates in the test set of cations and neutral compounds. With the adjustment using P-gp and MRP2 submodels, the RMSE and mean error of the prediction for the test set of cations and neutral compounds were greatly decreased to 4.88 and 0.65, respectively.
Overall, PDb data were gathered for 97 compounds in humans, ranging from 0% (diflunisal, leukotriene C4) to 88% (valsartan). Molecular weights of these compounds varied over a wide range from 151 (paracetamol, acetaminophen) to 1,202 (cyclosporine A). cLogP values varied from −3.36 (aztreonam) to 9.5 (metocurine) with a median of 2.1. There was no relationship between PDb and molecular weight or cLogP. A QSPKR model with R2 of 0.569 was constructed for all of the compounds, but the Q2 value was only 0.289, indicating that the model was not stable and, thus, could not be used to predict PDb. The compounds were then separated into two groups based on charge status: anions with 43 compounds and cation or neutral compounds with 44 compounds. The biliary excretion of the compounds in the human dataset was lower, on average, than that in the rat database: only five anions and one cation were excreted into bile in humans at an amount greater than 30% of the dose. The highest values in the anion group were 80% (indocyanine green) and 88% (valsartan), while the largest biliary excretion in the cation or neutral compound group was only 37.5% (hexafluorenium).
QSPKR models were constructed for the training sets of these two groups of compounds separately. A model composed of six structural descriptors was developed for anions (33 compounds):
(Fig. 7c). A Q2 of 0.672 produced by cross-validation suggested that the model could be used for prediction. Likewise, a QSPKR model with R2 of 0.895 and Q2 of 0.838 was established for 30 cations and neutral compounds (Fig. 7d):
SHBint3 and SHBint6_acnt reflect hydrogen bond interaction. MaxNeg reflects charge status, which was also selected in the QSPKR model for PDb of cation or neutral compounds in rats. Px and Q represent molecular moment properties. Gmin is related to electrotopological state. The predictability of the models was further evaluated with two test sets, i.e., 10 anions and 14 cations or neutral compounds. RMSE was calculated as 4.19 for cation/neutral compounds and 11.9 for anions, indicating that our models were accurate in prediction. In addition, there was only modest bias for the prediction indicated by small mean error (1.10 for cation/neutral compounds and 4.67 for anions).
QSPKR models were also evaluated for subsets of 14 MRP2 substrates and 14 P-gp substrates. All the P-gp substrates were cations or neutral compounds. Among the MRP2 substrates, half of them were anions. The submodels were developed as following:
Mean error and RMSE were calculated using predicted values obtained from the leave-one-out method, since there were not enough compounds available to evaluate test sets. Both models had small mean errors (1.98 for MRP2 submodel and 0.12 for P-gp submodel) and comparable RMSE with the models for the anions and the cation/neutral compounds (17.1 and 4.12 for MRP2 and P-gp submodels, respectively). Additionally, improvement of predictability by utilizing MRP2 and P-gp submodels was observed for their substrates in the test sets. A more accurate prediction for cefpriamide (observed value: 28.1%) was obtained from the MRP2 submodel (28.2%) compared with a value of 8.23% from the general model developed for anions; the P-gp submodel offered better prediction (13.2%) for vecuronium (observed value: 10.5%) which was underpredicted by the general model for cations and neutral compounds (6.15%).
Biliary excretion is an important route for the elimination of drugs and/or their metabolites and may contribute to the interindividual differences in response to drugs (32). However, it is difficult to measure biliary excretion in humans due to the invasive nature of studies and the limited techniques available to obtain good clinical data (2). Predicting biliary elimination in silico can add valuable information in the drug discovery and development process.
Historically, MW or lipophilicity has been used to predict biliary elimination of drugs (5). A MW cutoff of 500±50 in humans and 325±50 in rats has been generally accepted as reliable (4,5). Most of the research supporting these MW threshold values were performed in the laboratory of R.T. Williams using rats, rabbits, and guinea pigs. These researchers administered simple organic compounds or glucuronide conjugates, and many of the compounds were not drugs. Most of the compounds were administered intraperitoneally and at different doses to avoid toxicity in the animals. Only bile fluid was sampled, recovery of radioactivity was usually incomplete, and in many studies metabolites were not analyzed separately from parent compound. MW cutoffs for rats, guinea pigs, and rabbits were estimated to be 325±50, 400±50, and 475±50, respectively (33). These cutoffs were determined for organic anions, and do not apply to cations, zwitterions, or neutral compounds, a fact that is not widely recognized. A different cutoff of 200±50 was determined for monoquaternary cations and 500±50 for diquaternary cations in rats (3,6). The number of compounds used to determine the cutoff for cations was limited, and, therefore, less confidence can be placed in these MW threshold values.
In our evaluation of literature data obtained from rat studies, we found a MW threshold for anions of 400 Da based on ROC analysis of the data. ROC analysis has been widely used for diagnostic tests and/or decision making for a binary-classified system in medicine, radiology, and psychology (28). No significant MW threshold was determined for organic cations or organic cations/neutral compounds. The cutoff value for anions in rats was around 400 Da, higher than that reported in the literature. The difference in the MW threshold determined in rats may be attributed to several reasons: (1) our dataset was larger (N=164) and included compounds with a broader structural diversity than the compounds evaluated by the R.T. Williams group; (2) our data were extracted or calculated from literature only for parent drug, thus, avoiding the interference by metabolites and more reflecting the behavior of the compounds themselves; (3) our data were strictly collected from those studies in which the compounds were administered intravenously, therefore, circumventing any variability in the data due to absorption or first-pass biliary extraction; (4) our conclusion was based on several criteria including significance by Chi-square test, high AUC for the ROC curve, and generally accepted values for both sensitivity and specificity in ROC analysis. In contrast, there was no description concerning the derivation of the 325±50 threshold in this earlier analysis (5). In addition, since the majority of the compounds in our dataset are drugs, the threshold derived from our study may provide better predictions for biliary excretion in drug discovery and development.
There have been no studies dedicated to determining a MW cutoff for humans. The MW threshold for humans that is usually cited is 500±50 Da, a value that was originally suggested by Dr. Peter Millburn (4); however, the cutoff value was based on a very limited number of compounds eliminated in bile in humans and involved anions only. This cutoff value has been generally applied to all compounds, not just anions. We calculated a MW threshold in humans using ROC curve analysis, based on a much larger dataset (N=97) with compounds with MW values ranging from 199 to 1,127. Interestingly, our finding of a MW threshold of 475 Da for PDb for all compounds in humans (almost significant with a sensitivity of 80% and specificity of 79%) is very similar to the suggested value of 500 Da. A MW threshold of 475 Da was significant for anions but no threshold MW was identified for cations or cations/neutral compounds based on our statistical criteria.
Although molecular weight has been traditionally used as a determinant of biliary excretion of drugs, MOVO may be a more appropriate descriptor. MOVO has been previously used in other QSPKR analysis (34) and has been suggested as a determinant in binding to and transport by P-gp (35,36). A MOVO cutoff may represent a lower limit for the molecular size of a compound that can fit in the substrate binding site of an ABC or other transporter. Using ROC analysis, we found a MOVO threshold of 375 Å for the anion subset in humans and a potential value of 300 Å for anions in rats. MOVO was highly correlated with MW in our datasets in both rats and humans. MOVO values of 300 and 375 Å correlate to MW values of 400 and 475 Da, respectively, similar to the MW threshold values. This suggests that the molecular weight threshold for biliary excretion may reflect, at least in part, the molecular size of a compound.
It has been suggested in the literature that lipophilicity represents an important determinant of biliary elimination (5). We performed univariate analysis to determine if any relationship exists between PDb and Log P or Log D, and failed to find a relationship. We also generated ROC curves to determine if there is any relationship between CLb and Log P or Log D and between PDb and Log P or Log D. We found no relationship suggesting that Log P and Log D are not strongly correlated with the biliary elimination of drugs in rats and humans. Additionally, we were not able to identify a MW threshold for CLb. This may reflect the fact that two compounds (for example, pitvastatin and pravastatin) may have similar MWs and PDb values, but very different CLb values due to differences in the total clearances. PDb provides an estimate of the CLb, normalized for total CL. Further study is needed to understand these relationships.
MW is just one of the many factors involved in the biliary elimination of drugs. The volume of biliary fluid is very small, resulting in very high drug concentrations and bile/plasma concentration ratios greater than 1, but most commonly between 10 and 1000 (37). Thus, passive transport of compounds is negligible and hepatobiliary elimination requires efflux transporters to move compounds from hepatocytes into the canalicular space. ABC transporters that are important for the biliary elimination of drugs include P-gp, MRP2, and BCRP. Elucidation of the structural features necessary for binding to and transport by these ABC transporters may allow a better prediction of the biliary elimination of drugs. It is well accepted that P-gp transports mostly cation and neutral compounds while MRP2 transports predominantly anions, although they have substrates in common. Structure–activity relationships (SAR) and QSAR to uncover structural similarities in substrates for P-gp have been explored (38,39), but it has been challenging to describe the wide variety of P-gp substrates with a few descriptors. MRP2 and BCRP structure activity studies have been explored as well (40,41), although most studies have evaluated transport inhibitors. Additionally, most studies have considered datasets of congeneric series of compounds and the results cannot be extended to non-congeneric series of compounds.
We have applied an in silico method, QSPKR, to predict biliary excretion of non-congeneric compounds in rats and humans. In the past, QSPKR has been successfully employed to predict a number of pharmacokinetic properties, including clearance, volume of distribution, bioavailability, and protein binding (13–21). To our knowledge, this study is the first to apply such a method to predict biliary excretion. When first used, QSPKR models were predominantly used to examine relationships between pharmacokinetic properties and a single conventional physicochemical descriptor. The majority of the studies focused on the octanol-water partition coefficient, logP or logD, which was shown to be related with renal clearance or non-renal clearance of several series of compounds including beta-adrenergic blockers (42,43). Generally, renal clearance decreased with an increase in lipophilicity, while non-renal clearance showed a positive correlation with LogP or LogD. Thus, we also performed a univariate analysis trying to correlate biliary clearance with cLogP. However, a relationship was not observed in either rats or humans. Comparably, a recent analysis of pharmacokinetic parameters of 670 drugs in humans suggested that systemic clearance could not be quantitatively predicted from any single molecular property examined in that study (lipophilicity reflected by cLogP, polar surface area, sum of hydrogen acceptors, and donors) (44). One of the reasons for the lack of such a relationship may be the complexity associated with biliary excretion. To be excreted into bile, a compound has to be taken up into hepatocytes from blood across the sinusoidal membrane. Once entering the hepatocytes, parent compound may undergo biliary elimination, metabolism, or efflux back into blood. Another potential reason for the lack of relationship with a single molecular property is that the compounds evaluated in our study are not from a congeneric series. The diversity of the structures makes it difficult for one single descriptor to explain the pharmacokinetic properties of all the compounds. Our initial studies with a limited number of molecular descriptors failed to predict the biliary elimination of drugs in both rats and humans.
Considering the complexity of biliary excretion and the diversity of the compounds, we developed QSPKR models based on a large number of structural descriptors calculated using the principle of quantum mechanics. These 2D and 3D structural descriptors include molecular connectivity Chi indices, Kappa-shape indices, electrotopological state indices, information indices, subgraph count indices, molecular polarizability, weight, and volume. Previously, we have used the software (QSARis©) to determine structural descriptors for a group of flavonoids and developed a QSAR model to predict the inhibitory potency of flavonoids on BCRP (40). A QSPKR model was developed for biliary clearance for 37 compounds with data collected from studies in rats. A model with four parameters predicted CLb fairly well with a R2 of 0.689. Its prediction power was also demonstrated by internal validation generating a Q2 of 0.618, which is comparable to the one from a QSPKR study using PLS to develop a model predicting systemic clearance of 12 adenosine A1 receptor agonists in rats (20). More importantly, the predictability of the model was further verified by external validation when the model was tested with 18 compounds not present in the training set. The model predicted CLb of the test group reasonably well, producing a geometric mean fold error of 2.98, indicating that most of the predicted values fell between threefold and one third of the observed values. The removal of one outlier with a predicted negative value for CLb (UFX) resulted in a drop in the geometric mean fold error to 2.32, which was comparable to that reported from a physiologically based model used to predict total clearance and volume of distribution of 68 compounds in rats (30). In that study, the PBPK model was applied to six chemical classes and generated mean fold errors from 1.5 to 2.4 for total clearance with an average of 1.8. It should be noted that this PBPK model relied on the measurement of a series of parameters from in vitro experiments, such as solubility, permeability, and especially intrinsic clearance which had to be determined using hepatocytes, while our present QSPKR model is a pure in silico method. We found that QSPKR submodels built with MRP2 or P-gp substrates in the training set provided more accurate predictions of CLb than the general model, when evaluating MRP2 or P-gp substrates in the test group. We speculated that the improvement may be partially due to the similarities of the molecular properties shared among substrates of each transporter. With the modifications from these two submodels, the geometric mean fold error of the final prediction reached 1.87, indicating that the QSPKR models could provide very good quantitative predictions.
A QSPKR model was also successfully developed for human CLb data, based on a training set of 38 compounds, most of which were approved drugs including β-lactam antibiotics, quinolone antibiotics, non-steroid anti-inflammatory drugs, chemotherapeutic agents, and histamine H2-receptor antagonists. Predicted values for 29 compounds (76% of 38 compounds) fell within the threefold error range of observed values. Prediction performance of the model was also supported by the results of the leave-one-out method, with a high Q2 value of 0.722 suggesting that the model could be used for prediction purposes. Geometric mean fold error was calculated as 2.68 for another 18 compounds that represented the test set. This error value was comparable to those (1.84 to 2.55) from a QSPKR study using stepwise regression to construct models predicting apparent volume of distribution of a heterogenous series of drugs in humans (14). In addition, we also obtained good correlations when examining subsets of MRP2 or P-gp substrates; however, our data were limited. BCRP substrates could not be evaluated due to a lack of identified BCRP substrates in our dataset. Knowledge of ABC transporters involved in the biliary excretion of compounds may further improve the predictability of the general model for human CLb, based on our experience with CLb data in rats.
The majority of the biliary elimination data in the literature is provided as percentage of dose excreted into bile; in some cases, these data can be utilized to calculate biliary clearance based on total clearance data in the same study or more commonly using clearance data obtained from other investigations. Therefore, although our intent was to determine QSPKR models for CLb, an intrinsic property of a compound, we also looked at the possibility of obtaining QSPKR models using data on PDb. In rats and humans, significant correlation coefficients were obtained using all compounds; however, the predictive ability of the models (Q2=0.437 in rats and 0.289 in humans) was low.
QSPKR models were then evaluated by dividing the compounds into anion or cation/neutral groups. Our reasons for such subdivision were based on the following: first, it has been indicated by early published studies and our current investigation that there may be differences in the molecular weight cutoff for biliary excretion (for >10% of the dose eliminated in bile) between anions and cations, suggesting differences in the relationships between their physicochemical properties and PDb (3,6,33); secondly, there may be differences in the characteristics of MRP2 substrates (mainly anions) and P-gp substrates (largely cations and neutral compounds), although such separation based on charge status is not a perfect division of substrates of these transporters. When subdivided in this manner, better QSPKR models were obtained with greatly improved predictability, as suggested by good Q2 values generated from internal validation. Zwitterions were evaluated separately and demonstrated a significant QSPKR, although the data were limited (data not shown). In addition, QSPKR models were also successfully developed for the subsets of P-gp and MRP2 substrates. When the models developed for anions or cations/neutral compounds in rats and humans were applied to a number of compounds not included in the training set (i.e., test set), they provided good predictions, as suggested by small mean error and reasonable RMSE values. Improvement of predictability by utilizing MRP2 and P-gp submodels will be further confirmed when clinical data become available for other substrates of these transporters.
Importantly, it is necessary to emphasize that the relationships were established for parent drugs only. In some cases, metabolites, such as glucuronide conjugates, undergo extensive biliary excretion. However, we did not evaluate metabolite biliary excretion unless a preformed metabolite (e.g., UFX) was administered. Additionally, there are a number of limitations for our current study. First, although significant relationships were derived without considering whether a compound is a substrate for any transporter in our general models, we suggested that relationships might be improved by considering whether the drug was substrate for MRP2 or P-gp. We did not examine the potential influence of hepatic uptake transporters, although hepatic uptake may be rate-limiting for some substrates (e.g., pravastatin). Second, although protein binding is an important factor in determining the pharmacokinetics of drugs, we did not consider the effects of protein binding on the biliary elimination of drugs. To fully consider the effects of protein binding, it may be necessary to use intrinsic clearance values for biliary elimination in the derivation of QSPKR models; however, this information is not readily available in the literature. It is interesting to note, though, that QSPKR models for PDb should not be influenced by protein binding. Third, a linear relationship was assumed between pharmacokinetic parameters and independent descriptors, an assumption of many other QSPKR studies, since linear regression is often adequate to predict a response reasonably well. Nonlinearity can be taken into account by adding quadratic or cross-product terms of independent variables in the linear regression models, and this approach was demonstrated to increase model fitting and predictability for some parameters significantly in a QSPKR study predicting distribution parameters of 17 drugs (17). Lastly, since the descriptors were obtained based on principles of quantum mechanics, it is difficult to interpret the selected descriptors in terms of physicochemical and structural parameters. Nevertheless, the molecular property categories to which these descriptors belong suggest that some chemical features such as molecular size, charges, and hydrogen bond formation may be important for biliary excretion, which is consistent with the previous dogma that high MW and the presence of polar groups favor biliary excretion (45). For example, xc3, xch5, and xvch7 were present in the models for PDb in rats. These molecular connectivity indices represent mainly the size of molecules. MaxNeg, which is the largest negative charge over the atoms in molecule, was shared by models for PDb data of cations and neutral compounds groups in rats and humans. Both models also contained descriptors associated with hydrogen bond interaction (SHBint9_acnt in the model for rats, SHBint6_acnt, and SHBint3 for humans). In addition, parameters relevant to dipole or quadrupole moments were included in most models for PDb and CLb, indicating the contribution of van der Waals force. However, ambiguity in the physicochemical meanings of these descriptors makes it difficult to directly guide specific chemical modifications of existing molecules. Nonetheless, such a detailed interpretation may be unnecessary for a QSPKR model to be used mainly for prediction purposes instead of drug design.
In conclusion, we compiled literature data on MW and biliary elimination of drugs and derived ROC curves to establish the presence of MW thresholds for the biliary elimination of drugs for anions in rats and humans. MW threshold values for organic anions in rats of 400 Da and in humans of 475 Da were determined; cations or cations/neutral compounds did not demonstrate statistically significant MW threshold values for biliary excretion. QSPKR models were successfully developed for biliary clearance and percentage of dose excreted into bile using 2D and 3D quantum chemical descriptors. Internal and external validation suggested that the derived models could provide quantitative prediction for other compounds. Importantly, a QSPKR model with good predictability for the biliary clearance of drugs in humans has been successfully developed. Many QSAR/QSPKR studies have been conducted with a series of analogs and, thus, may have limitations with predicting the activities of structurally unrelated compounds, while our models were constructed based on non-congeneric compounds with diversified structures and hence may have broader applicability.
This work was supported in part by Pfizer Inc. We thank Dr. Lisa J. Benincosa from the Groton/New London Laboratories at Pfizer Inc. for her support and suggestions regarding this research. We thank Dr. Xueya Cai from Division of Biostatistics of Indiana University for her help with receiver operating characteristic curves analysis.
The complete dataset of compounds with biliary clearance and percentage of dose eliminated in bile data, in both rats and humans, along with information on administration route, collection period of samples, analytical method, MW of the compounds, charge status of the compounds under physiological conditions (pH=7.4) and transporters known to mediate efflux processes of the drugs, are available as Excel files. A table containing the QSARis© descriptors is also present. This information can be found at http://www.buffalo.edu/~memorris