drug induced liver 
 injury ( dili ) is a major issue worldwide , both 
 for patients and health providers .
it is one of 
 the primary causes for attrition during clinical and preclinical studies 
 and the main reason for drug withdrawal from the market .
dili is divided into types , ( i ) hepatocellular , ( ii ) cholestatic , 
 or ( iii ) mixed ( hepatocellular and cholestatic ) , according to the 
 type of liver damage and the clinical chemistry biomarker alterations .
the cholestatic and mixed hepatocellular and 
 cholestatic type are the two most severe manifestations of dili and 
 yield almost half of the recorded cases of dili .
cholestatic liver injury , or more simply cholestasis , is the disruption 
 of the bile flow , which might be either due to biliary tract obstruction 
 or to complications in bile acid uptake . while the mechanistic basis 
 for hepatocellular dili is still a mystery for the majority of the 
 cases , more knowledge exists for cholestatic dili .
there is growing 
 evidence for a vast amount of cholestasis cases pinpointing the important 
 role of hepatic transporters .
basolateral transporters are responsible for the uptake of drugs 
 and other endobiotics and xenobiotics from the blood , influencing 
 the exposure of the hepatocyte to potential damage .
canalicular transporters 
 regulate the hepatic clearance , as well as the secretion of bile salts 
 and bile conjugates into bile .
any disturbance of the transporters physiological function 
 may result in the accumulation of potentially harmful bile products 
 that can finally cause cholestasis.figure 1 provides an overview 
 on the respective location of hepatocyte transporters .
blue symbols represent 
 mainly the canalicular transporters , and red symbols , the basolateral 
 ones .
mrp16 
 multidrug resistance - associated proteins 16 , ost/ost 
 organic solute transporter , bsep bile salt export pump , bcrp breast 
 cancer resistance protein , mate1 multidrug and toxin extrusion transporter 
 1 , abcg5/g8 atp - binding cassette subfamily g member 5/8 , mdr3 multidrug 
 resistance protein 3 , p - gp p - glycoprotein , atp8b1 atpase - aminophospholipid 
 transporter , oatp organic anion transporting polypeptide , ntcp sodium 
 ( na ) taurocolate cotransporting polypeptide , oct organic 
 cation transporter 1 , oat organic anion transporter .
the most important one , due to its pivotal role 
 in bile salts clearance , is the bile salt export pump ( bsep ) . apart from bsep
, there is evidence for the implication of other 
 canalicular efflux transporters such as the multidrug resistance - associated 
 protein 2 ( mrp2 ) , breast cancer resistance 
 protein ( bcrp ) , multidrug resistance protein 3 ( mdr3 ) , and p - glycoprotein ( p - gp ) .
mdr3 functions as an atp - dependent 
 phospholipid flippase , translocating phosphatidylcholine from the 
 inner to the outer canalicular membrane .
canalicular phospholipids 
 are then solubilized by canalicular bile salts to form mixed micelles , 
 protecting cholangiocytes from the detergent properties of bile salts . 

while p - gp is also not transporting bile salts , it is implicated in 
 cholestasis because of its large amount of substrates and inhibitors 
 which cause drug  drug interactions that disrupt the smooth 
 function of the hepatocyte .
the basolateral 
 transporters play also an important role , both the uptake transporters , 
 such as organic anion transporting polypeptides 1b1 , 1b3 , and 2b1 
 ( oatp1b1 , 1b3 , and 2b1 ) and sodium ( na ) taurocolate cotransporter 
 ( ntcp ) , and the efflux transporters , like multidrug resistance - associated 
 protein 3 and 4 ( mrp3 and mrp4 ) . in 
 particular , in cases of cholestasis , the basolateral uptake transporters 
 ntcp and oatp1b1 have been found down - regulated .
however , in this case , oatp1b3 is up - regulated as a compensatory 
 mechanism for the elimination of xenobiotics from sinusoidal blood . on the contrary , in cases of 
 cholestasis , mrp3 and mrp4
are up - regulated to facilitate the efflux 
 of the toxic bile salts out of the hepatocyte .
thus , simultaneous inhibition of several of these transporters 
 could induce drug toxicity due to inadequate elimination from the 
 blood or increase the cholestatic effect due to accumulation of bile 
 salts in the hepatocyte .
consequently , drug - induced liver injury 
 and cholestasis are important 
 toxicity alerts to be considered in drug development .
interestingly , 
 there are only a few computational studies for the prediction of cholestasis 
 reported in literature . with respect to the 
 involvement of hepatic transporter ,
there are some in vitro studies 
 correlating cholestasis with transporter inhibition , such as bsep , mrp3 , mrp4 , and ntcp .
also several in 
 silico studies for the identification of potentially cholestatic compounds 
 via modeling of transporters and then associating them with the cholestatic 
 effect of their inhibitors have been conducted .
a characteristic example 
 is the study by greupink et al . in 2012 , who developed a pharmacophore 
 approach for ntcp in order to identify 
 potentially ntcp inhibitors . under the same principles , in 2014 ritschel 
 and colleagues performed a 3d ligand - based 
 pharmacophore model for bsep inhibition .
however , in most of these 
 cases the amount of validated drugs is small and what is basically 
 described is the association between transporter inhibition and cholestasis . 

thus , as the respective is associated with cholestasis , it is assumed 
 that an inhibitor is causing cholestasis .
most recently , muller et 
 al . , in order to model dili , also modeled 
 some more hepatotoxicity end points , including cholestasis .
moreover , 
 mulliner et al . presented a multilevel 
 modeling approach for dili , where cholestasis was also included as 
 a morphological hepatobiliary finding .
however , examining the liver 
 transporters contribution was not within the scope of their work . 

dilisym ( www.dilisym.com ) is a mechanistic mathematical model of dili , that has been used 
 to investigate the effects of bsep inhibition on drug - induced liver 
 injury , as well as on bile acid - mediated 
 dili .
developed 
 a multiscale , liver - centric in silico modeling framework for acetaminophen 
 pharmacology and metabolism that can be extended in predicting hepatotoxicity 
 due to acetaminophen overdosing . in this study
we present a classification scheme in order to predict 
 cholestasis from a public data set , using physicochemical descriptors 
 as well as predicted transporter inhibition profiles . for the latter 
 we used our in house classification models for bsep , bcrp , p - gp , oatp1b1 , and oatp1b3 .
for compiling the dili training data set 
 we searched in pubmed ( http://www.ncbi.nlm.nih.gov/pubmed ) , google , scopus ( https://www.scopus.com/ ) , and the sider database v2 using the search terms :  drug - induced cholestasis  
 or 
the retrieved publications were then investigated manually for human data , i.e. compounds that 
 are positive or negative for drug - induced cholestasis in humans .
those 
 compounds  in principle drugs  were added to those obtained 
 from the sider v2 database .
unfortunately , cholestasis 
 is an end point that is not widely examined in terms of experimental 
 or in silico studies that would potentially guide us to big data sets . 

thus , even though we were able to compile several drugs positive for 
 cholestasis , there was almost no information in terms of the negatives . 

on the other hand , dili in general is studied quite extensively and 
 there are several respective data sets . since choleastasis is a possible 
 manifestation for dili , we can consider safely that any compound negative 
 for dili will definitely be also negative for cholestasis .
thus , the 
 negative compounds for dili that we had compiled and curated in a 
 previous work were also used as negatives for this study .
the data 
 set was carefully curated according to the following rules : all inorganic compounds were removed 
 based on their 
 chemical formula in moe 2014.09.salt parts and compounds containing metals 
 and/or rare 
 or special atoms were removed and the chemical structures were standardized 
 using the standardiser tool created by francis atkinson.duplicates and 
 permanently charged compounds were removed 
 using moe 2014.09 .
here we have to note 
 that stereoisomers , even if biologically can be considered as different 
 compounds , were considered as duplicates , since they give the exactly 
 same vector of descriptors .
if two ( or more ) stereoisomers are of 
 the same class , only one was kept .
if they were of different classes , 
 they were all removed.3d structures 
 were generated using corina ( version 3.4 ) , and their energy was minimized with moe 2014.09 , using default settings , but changing the gradient 
 to 0.05 rms kcal/(mol a ) .
all inorganic compounds were removed 
 based on their 
 chemical formula in moe 2014.09 .
salt parts and compounds containing metals 
 and/or rare 
 or special atoms were removed and the chemical structures were standardized 
 using the standardiser tool created by francis atkinson .
here we have to note 
 that stereoisomers , even if biologically can be considered as different 
 compounds , were considered as duplicates , since they give the exactly 
 same vector of descriptors .
if two ( or more ) stereoisomers are of 
 the same class , only one was kept .
3d structures 
 were generated using corina ( version 3.4 ) , and their energy was minimized with moe 2014.09 , using default settings , but changing the gradient 
 to 0.05 rms kcal/(mol a ) .
in addition , the existing chirality 
 was preserved . after these curation steps 
 152 compounds remained as 
 positives for cholestasis .
the negatives for dili , and subsequently 
 for cholestasis , were 466 compounds . however , when uniting the data , 
 there were compounds with contradictory class assignments .
these compounds 
 were removed from the data set , yielding a data set of in total 578 
 compounds ( 131 positives and 447 negatives ) . the compiled data set 
 is provided in the supporting information . recently  and after having 
 already compiled our training set for cholestasis and developed the 
 respective model  a data set covering multiple levels of hepatotoxicity 
 was published by mulliner and co - workers .
the data are hierarchically clustered by the authors into three 
 levels of hepatotoxicity : level 0 corresponds to general hepatotoxicity , 
 level 1 corresponds to clinical chemistry findings and morphological 
 finding as distinguished parts of general hepatotoxicity , and level 
 2 discriminates both clinical chemistry and morphological findings 
 into hepatocellular and hepatobiliary injury .
we use only clinical 
 data , i.e. the human data , of morphological findings for hepatobiliary 
 injury as an external test set for validating the developed cholestasis 
 model . once more , we performed chemical data curation and removed 
 the compounds overlapping with the training set , which led to 1347 
 compounds ( 230 positives and 1117 negatives ) as an external test set .
we also merged the training with the 
 test set and tried to generate a model based on the united data .
the 
 merged data set comprises 1904 compounds : 355 positives and 1549 negatives 
 for cholestasis . here
we must note that the total number of compounds 
 in the merged data set is not the same as the sum of the compounds 
 of training set and test set .
the reason is that when removing the 
 overlapping compounds between training set and test set , all of the 
 compounds were removed from the test set , since we had already selected 
 the final model for cholestasis . from those overlapping compounds , 

this did not matter so much , since 
 they were totally removed from the test set .
however , when merging 
 the two data sets for modeling , we did not want to decide regarding 
 the class of those compounds . for both data sets , several types 
 of molecular descriptors have been calculated , such as all 192 2d 
 moe descriptors , the 3d volsurf series of descriptors , as well as ecfps ( extended connectivity fingerprints ; 
 ecfp6 ) , using rdkit ( http://www.rdkit.org/ ) .
the list of the final descriptors 
 used for the proposed model is given in the supporting information ( table s1 ) .
in addition to this , predicted hepatic 
 transporter inhibition profiles were also included in the list of 
 descriptors .
the transporters investigated comprise bsep , p - gp , bcrp , 
 oatp1b1 , and oatp1b3 . in particular , for basolateral transporters 
 we calculated the predictions for four in silico classification models 
 built upon padel descriptors for oatp1b1 
 and oatp1b3 inhibition .
for obtaining 
 the predictions we use the models version implemented in etoxlab , an open source modeling framework for implementing 
 predictive models . out of each model
we use the sum of these binary scores , 
 denoted  sum binary score  . the sum binary score can 
 take values between 0 ( if all models predict the compound as negative ) 
 and 4 ( if all models predict the compound as positive ) . for basolateral 
 transporters
a more thorough 
 description of the transporters inhibition models is provided in the supporting information ( table s2 ) , where the 
 size of the training set , the inhibition threshold of the training 
 set , and the algorithm and performance of each individual model based 
 on auc values are provided .
the two - class classification models 
 were built using the software package weka ( version 3.7.12 ) .
we investigated the performance of several base 
 classifiers , such as logistic regression , tree methods ( random forest 
 and j48 tree ) , support vector machines ( smo in weka with polynomial , 
 rbf , and puk kernels ) , nave bayes , and k - nearest 
 neighbors .
moreover , because the data set is slightly imbalanced , 
 in order to equilibrate the effect of the majority class on model 
 performance , we also applied the cost - sensitive meta - classifier metacost .
the cost matrix applied corresponds to the imbalance 
 ratio of the data ( 3:1 ) .
additionally , several meta - classifiers were 
 explored for attribute selection ( attributeselectedclassifier ) , as 
 well as for improving the statistical performance , such as bagging and boosting .
the models were originally validated 
 via 10-fold cross validation , which is considered a quite trustworthy 
 method of validation . the best models  according 
 to 10-fold cross - validation evaluation  were further validated 
 via using the data set by mulliner .
subsequently , 
 for the best obtained models , 50 iterations were performed by changing 
 the cross - validation seed ( for splitting the data within cross validation ) 
 and the respective performance parameters were calculated . in order 
 to compare whether the inclusion of the transporters predictions in 
 the descriptors set improves significantly the model s performance , 
 a two - sample t test was performed in r. the statistics metrics taken into consideration 
 were accuracy , sensitivity , specificity , matthews correlation coefficient 
 ( mcc ) , area under the curve ( auc ) , precision and weighted average 
 precision .
apart from the conventional metrics of accuracy , sensitivity , 
 and specificity , we also use auc , since it is a measure of the models 
 capability to rank the predictions , while it decouples classifiers 
 from class imbalance and error costs .
moreover receiver operating 
 characteristic ( roc ) graphs are very useful for visualization of the 
 models results .
the mcc value , considering 
 its formula , takes into account all values of the confusion matrix : where , 
 tp is true positives , fp is false positives , 
 tn is true negatives , and fn is false negatives .
thus , it is considered 
 more balanced and informative than the column- or linewise metrics .
weighted average precision is the average precision 
 obtained for the two classes but weighted from the total number of 
 instances of the classes .
it is a quite 
 helpful parameter in multiclass classification problems , as well as 
 for imbalanced data sets where the number of negatives is greater 
 than the number of positives .
especially for the latter case , due 
 to the definition of precision [ ppv = tp/(tp + fp ) ] , its value for 
 the positive class would be low , which not necessarily means that 
 the total performance of the model is bad .
of course , since we are 
 dealing with a toxicity classification problem , like cholestasis , 
 the metrics that is of particular interest and that should by no means 
 drop below 0.5 is sensitivity or true positive rate . in order 
 to be confident regarding the validity of the models we used , we investigated 
 the coverage of the transporters models for the cholestasis data . 

additionally , we checked how reliable the predictions of the cholestasis 
 model for the cholestasis test set are .
the applicability domain was 
 checked on knime with the enalos nodes that compute 
 the applicability domain on the basis of the euclidean distances .
the number of compounds within the model s 
 applicability domain for each model and for each cholestasis data 
 set is provided in the supporting information ( table s3 ) .
several 
 combinations of descriptors and classifiers were investigated and 
 the optimal classification model was selected on the basis of the 
 results of 10-fold cross validation . with respect to the classifier , 
 the best results were obtained using as base classifier ibk  the k - nearest neighbors implementation 
 in weka  with k = 5 .
the meta - classifier metacost 
 was also applied , with the application of the cost matrix [ 0.0 , 1.0 ; 
 3.0 , 0.0 ] , i.e. weighting the minority class 3 times more than the 
 majority class , in order to cope with the slightly imbalanced training 
 set .
2d moe descriptors were performing better than fingerprints and/or 
 volsurf descriptors , especially for sensitivity , mcc and auc . combining 
 the volsurf descriptors with 2d moe descriptors
also did not provide 
 any significant improvement of the results . from the whole set of 
 2d moe descriptors we decided to use a subset of 93 interpretable 
 descriptors that give almost the same performance compared to using 
 all 2d moe descriptors .
apart from the 93 2d descriptors , we also 
 included the predicted transporter inhibition profiles . in order to 
 assess the importance and significance of this additional information 
 individually , we used them in different combinations : all transporters , 
 only bsep , all transporters excluding either bsep , or p - gp , or bcrp , 
 or the oatps . this led to in total seven models ( table 1 ) .
it is interesting to mention that 
 we also exchanged the training - test 
 set roles between the two data sets and tried to generate a model 
 for cholestasis based on the bigger data set from the work of mulliner 
 et al .
however , even though the results 
 for 10-fold cross validation were equivalent to the model we propose 
 ( generated on the compiled training set of 578 compounds ) , the results 
 for the external validation of the 578 compounds were rather poor 
 ( results not shown ) .
furthermore , we tried to generate a cholestasis model on the merged 
 data for the subset of 93 2d moe descriptors , with or without the 
 transporters interaction profiles .
interestingly , the k - nearest neighbor settings ( k = 5 ) , which gave quite 
 satisfactory results for 10-fold cross validation while modeling either 
 the training or the test set standalone , did not have the same effect 
 for the united data . for
the merged data set svm ( smo implementation 
 in weka ) using a polynomial kernel , with exponent equal to 2 , performs 
 better .
the use of metacost with a cost matrix of [ 0.0 , 1.0 ; 5.0 , 
 0.0 ] , due to the new imbalance ratio of the data , is also necessary . 

additionally , under these settings , the performance of the model is 
 significantly better when using the transporters predictions as additional 
 descriptors .
the obtained performance of this model , as well as the 
 respective p - values of the performed two - sample paired t test out of 50 iterations , is presented in the supporting information ( table s4 ) .
inspecting 
 the obtained results in table 1 , it becomes obvious that the best settings 
 for the model for 10-fold cross validation are achieved with the inclusion 
 of all transporter inhibition predictions in the list of descriptors . 

nevertheless , this is not the case for the external validation , where 
 including predicted inhibitor profiles for all transporters yields 
 lower accuracy and specificity values , while sensitivity remains almost 
 the same .
interestingly , the use of bsep inhibition prediction stand - alone 
 does not seem to be sufficient .
there is a drop in the statistics  especially 
 for sensitivity  in comparison to the use of the whole set of 
 transporter predictions , both for 10-fold cross - validation and for 
 the external test set . in order to assess if the predicted transporter 
 inhibition profiles indeed statistically significantly improve the 
 models , we performed 50 iterations of 10-fold cross validation followed 
 by a two sample t test on the performance parameters . 

for this we used the models with 2d moe descriptors , 2d moe plus all 
 transporters , 2d moe plus bsep , and 2d moe plus all transporters without 
 bsep ( table 2 ) . analyzing the p - values 
 for the pairwise comparisons 
 (
supporting information table s5 ) the main 
 conclusion is that indeed the use of liver transporter inhibition 
 predictions contributes significantly to the models performance when 
 compared to the use of only 2d physicochemical descriptors .
interestingly , 
 it is not only the bsep inhibition contribution , which matters . on 
 the contrary , when only bsep predictions are used ,
there is only a 
 slight increase in sensitivity and specificity of the model , while 
 sensitivity decreases .
the other way round , if all transporters predictions 
 are used , sensitivity increases , but accuracy and specificity reach 
 their peak only after the inclusion of bsep predictions .
nevertheless , 
 it does not contain all the important information , despite the fact 
 that it is the most widely discussed transporter in literature with 
 respect to cholestasis . a possible explanation for this interesting result could be the 
 different thresholds required for bsep inhibition to cause cholestasis 
 versus the threshold of the bsep inhibition model . for the bsep model , 

every compounds with an ic50 > 50 m was classified 
 as noninhibitors , while ic50 < 10 m was the threshold 
 for characterizing a compound as inhibitor .
however , for the development 
 of cholestasis , an ic50 < 300 m has been reported 
 as sufficient .
the other transporters 
 included in our study are in general not well described in literature 
 via experimental procedures , but they are rather pinpointed due to 
 the fact that they are transporting bile salts or bile conjugates 
 ( with the exception of p - gp , whose role is attributed mainly to drug 
thus , our study gives extra weight to literature indications 
 concerning bcrp , p - gp , oatp1b1 , and 1b3 .
it was quite curious 
 that even though the transporters predictions 
 significantly increase the performance both for the model built on 
 the training set of 578 compounds and for the model trained with the 
 merged data set of 1904 compounds ( but with a different base classifier ) ; 
 this was not the case for the prediction on the test set .
the performance 
 of all statistics metrics decrease , when transporters predictions 
 are used as descriptors .
we can only speculate a different way of class 
 assignment between the two data sets , since they are coming from different 
 sources .
another plausible explanation could be the fact that the 
 external validation set had a contradiction rate of almost 20% regarding 
 the class labels of those compounds shared with the training set ( 71 
 out of 419 shared compounds had contradictory class labels ) .
we assume 
 that this is due to the drawbacks of the toxicity reporting system : 
 under - reporting , voluntarily basis , difficulties to obtain 
 the data , which are often proprietary , as well as the lack of the prerequisite of a causal relationship 
 between drug and adverse event . in any 
 case , despite these contradictions between the training and the test 
 set , the model retained its satisfactory performance . furthermore
, 
 we would like to mention that there is also experimental 
 evidence for the implication of the basolateral efflux transporters 
 mrp3 and mrp4 , as well as for the 
 canalicular efflux transporter mrp2 .
we are aware of this fact , but for these transporters there are 
 currently not sufficient data available to develop high quality models 
 that can be further used for contributing to our cholestasis model . 

additionally , we should pinpoint the fact that we are using predictions 
 for the inhibition of transporters rather than real in vitro data . 

this could potentially accumulate some additional error in the final 
 predictions of the cholestasis model . in any case , despite any deteriorating 
 factors , our final in silico models for cholestasis were extensively 
 validated with 10-fold cross validation and statistical tests .
furthermore , 
 the external validation set was of a significant size being even bigger 
 than the training set . for both cases
moreover , for both training sets ( regular and merged one ) there is 
 a trend indicating the importance of transporter predictions in the 
 development of cholestasis .
the performance of the classification 
 models from which we obtained the predictions provides us with enough 
 confidence to use them in our input matrix .
moreover , the percentage 
 of the cholestasis data that are within the applicability domain of 
 the transporters models ranges between 93.1% and 99.5% ( table s3 , supporting information ) , which is quite satisfactory 
 and further enhances our confidence in using predicted transporter 
 interaction profiles as descriptors .
in 
 this study we present a two - class classification model for the 
 prediction of cholestasis ( or cholestatic dili ) based on a public 
 data set of 578 compounds .
the model is incorporating information 
 both from 2d physicochemical descriptors , as well as predictions of 
 inhibition of the hepatic transporters bsep , bcrp , p - gp , oatp1b1 , 
 and oatp1b3 .
the performance of the resulting model is rather satisfactory 
 and is validated both via 50 iterations of different 10-fold cross 
 validations , as well as an external test set of over 1300 compounds .
our results demonstrate that adding transporter 
 predictions as additional descriptors to the list of 2d physicochemical 
 descriptors is improving model performance .
this is in alignment with 
 evidence from literature which shows that inhibition of selected hepatic 
 transporters contributes to cholestasis .
interestingly , the 
 increase in model performance can not be attributed 
 solely to bsep inhibition , which is the transporter that is most correlated 
 to cholestasis in literature .
the performance increases only when 
 the whole list of transporter inhibition predictions is included . 

this result points toward a rather synergistic effect of several transporters , 
 including the less elucidated role of oatps , bcrp , and p - gp in cholestasis .
our study is the first of its kind regarding combining physicochemical 
 descriptors and predicted transporter information in order to predict 
 cholestasis .
 cholestasis represents 
 one out of three types of drug induced liver 
 injury ( dili ) , which comprises a major challenge in drug development . 
 in this study 
 we applied a two - class classification scheme based on k - nearest neighbors in order to predict cholestasis , using 
 a set of 93 two - dimensional ( 2d ) physicochemical descriptors and predictions 
 of selected hepatic transporters inhibition ( bsep , bcrp , p - gp , 
 oatp1b1 , and oatp1b3 ) . in order to assess the potential contribution 
 of transporter inhibition , we compared whether the inclusion of the 
 transporters inhibition predictions contributes to a significant 
 increase in model performance in comparison to the plain use of the 
 93 2d physicochemical descriptors . 
 our findings were in agreement 
 with literature findings , indicating a contribution not only from 
 bsep inhibition but a rather synergistic effect deriving from the 
 whole set of transporters . 
 the final optimal model was validated via 
 both 10-fold cross validation and external validation . 
 it performs 
 quite satisfactorily resulting in 0.686  0.013 for accuracy 
 and 0.722  0.014 for area under the receiver operating characteristic 
 curve ( auc ) for 10-fold cross - validation ( mean  standard deviation 
 from 50 iterations ) . 
Introduction
Methods
Results and Discussion
Conclusions