Development of a mechanistic biokinetic model for hepatic bile acid handling to predict possible cholestatic effects of drugs

Abstract Drug‐induced liver injury (DILI) is a common reason for drug withdrawal from the market. An important cause of DILI is drug‐induced cholestasis. One of the major players involved in drug‐induced cholestasis is the bile salt efflux pump (BSEP; ABCB11). Inhibition of BSEP by drugs potentially leads to cholestasis due to increased (toxic) intrahepatic concentrations of bile acids with subsequent cell injury. In order to investigate the possibilities for in silico prediction of cholestatic effects of drugs, we developed a mechanistic biokinetic model for human liver bile acid handling populated with human in vitro data. For this purpose we considered nine groups of bile acids in the human bile acid pool, i.e. chenodeoxycholic acid, deoxycholic acid, the remaining unconjugated bile acids and the glycine and taurine conjugates of each of the three groups. Michaelis‐Menten kinetics of the human uptake transporter Na+‐taurocholate cotransporting polypeptide (NTCP; SLC10A1) and BSEP were measured using NTCP‐transduced HEK293 cells and membrane vesicles from BSEP‐overexpressing HEK293 cells. For in vitro‐in vivo scaling, transporter abundance was determined by LC‐MS/MS in these HEK293 cells and vesicles as well as in human liver tissue. Other relevant human kinetic parameters were collected from literature, such as portal bile acid levels and composition, bile acid synthesis and amidation rate. Additional empirical scaling was applied by increasing the excretion rate with a factor 2.4 to reach near physiological steady‐state intracellular bile acid concentrations (80 &mgr;M) after exposure to portal vein bile acid levels. Simulations showed that intracellular bile acid concentrations increase 1.7 fold in the presence of the BSEP inhibitors and cholestatic drugs cyclosporin A or glibenclamide, at intrahepatic concentrations of 6.6 and 20 &mgr;M, respectively. This simplified model provides a tool for a first indication whether drugs at therapeutic concentrations might cause cholestasis by inhibiting BSEP. Graphical Abstract Figure. No caption available.


Introduction
Drug-induced liver injury (DILI) is a major cause for drug withdrawal from the market. The most common form is cholestatic DILI, which can be distinguished into three phenotypes: cholestatic, hepatocellular and mixed. One of the mechanisms of cholestatic DILI, as proposed by Hofmann (1999) and Wagner et al. (2009), includes the inhibition of canalicular hepatocellular bile acid (BA) excretion, resulting in the accumulation of bile acids in the hepatocytes (Hofmann, 1999;Wagner et al., 2009). If toxic intracellular bile acid concentrations are reached, this can lead to liver injury due to necrotic and/or apoptotic cell death. Several transporters are involved in the hepatocellular transport of bile acids. For uptake, the sodium-dependent taurocholate transporter (NTCP, SLC10A1) and the organic anion transporters (OATPs, SLCO) are assumed to be important players. Multidrug resistance related proteins 3 and 4 (MRP3, ABCC3; MRP4, ABCC4) and the organic solute transporter (OSTα/β; SLC51A/SLC51B) are considered to be involved in sinusoidal efflux (Dawson et al., 2009), which is low under normal physiological conditions. Most importantly, the canalicular efflux transporter bile salt export pump (BSEP) seems to be the key player in the development of cholestatic DILI. This has been summarized by Vinken et al. in an Adverse Outcome Pathway (AOP), and further confirmed by Woodhead et al. by sensitivity analysis of a mechanistic model of cholestatic DILI (Vinken et al., 2013;Woodhead et al., 2014).
To predict cholestatic effects of newly developed drugs, animal models are currently used. However, there is increasing awareness that the current experimental animal models are not predictive for the human clinical situation. This discrepancy can be due to species differences in bile acid pool composition (García-Cañaveras et al., 2012;Setchell et al., 1997), in individual bile acid concentrations and in abundance and functionality of the transporters and enzymes involved in bile acid handling, as well as in their susceptibility for interference by the drug (Ulzurrun et al., 2013). Moreover, the animals used for drug testing are genetically more uniform than the human patient population and this can result in skewed results. Experimental human in vitro or ex vivo models have been developed for testing cholestatic effects of drugs, such as sandwich and 3D spheroid cultures of human hepatocytes (Guo et al., 2016;Hendriks et al., 2016), human liver slices (Vatakuti et al., 2016a(Vatakuti et al., , 2016b and human liver cell lines (HepaRG) (Hendriks et al., 2016;Qiu et al., 2016). Next to these animal and human models, in silico models can be applied to contribute to the preclinical prediction of cholestasis in man. An in silico model was recently described using DILIsym (Woodhead et al., 2014), in which BA handling in rat and human was modelled in the presence of a cholestatic drug (glibenclamide for rats and a theoretical cholestatic agent for humans). However, in this human model, many of the mechanistic parameters on hepatocellular bile acid handling were not experimentally determined, but were optimized by reverse-modeling, taking the clinical BA profile in humans as a starting point. To date, molecular pharmacological data on bile acid handling are not available in the literature. Kinetic parameters for the interaction of human transporters (NTCP, BSEP and OATPs) and bile acids are either not published with appropriate scaling factors or at best are available only for taurocholate, which is a minor BA in humans. In addition, Woodhead et al. showed that BSEP abundance is the most important unknown kinetic parameter, that has a great influence on the outcome of the BA concentration in the liver in the presence of a BSEP inhibitor (Woodhead et al., 2014).
Here, we developed a bottom-up human biokinetic model of hepatocyte bile acid handling. This model is based on existing and new human in vitro kinetic data on hepatocellular bile acid handling, as well as abundance data of key bile acid transporters. We aimed to predict the effect of competitive BSEP inhibitors on bile acid accumulation in the human hepatocyte. As prototypical BSEP inhibiting drugs we used glibenclamide and cyclosporin A. Both drugs have been associated with cholestatic DILI and are inhibitors of BSEP (Byrne et al., 2002;Goodman et al., 1987;Morgan et al., 2010;van Basten et al., 1992;Varma et al., 2014).

Selection of Bile Acids Incorporated in Model
The following bile acids were selected to be included in the biokinetic model: deoxycholic acid (DCA), chenodeoxycholic acid (CDCA), and their glycine-and taurine-amidated forms (GDCA, TDCA for DCA, and GCDCA and TCDCA for CDCA). Together, they account for the majority of the bile acids present within the human bile acid pool and represent all three classes of bile acids of the human body: primary, secondary and conjugated bile acids. Because kinetic data on transporter-mediated handling of these bile acids were lacking, we determined the affinity (Km) and transport capacity (Vmax) and the expression of the transporters, as described below. The remaining bile acids in the human bile acid pool were combined into 3 groups for modeling purposes: unamidated, glycine amidated and taurine amidated.

Development of the Mechanistic Biokinetic Model
In order to calculate the bile acid concentrations in the hepatocyte, passive and active bile acid uptake and efflux were taken into account as well as bile acid synthesis and amidation. Thus the intracellular concentration over time was calculated for unamidated bile acids (Eq. (1)) and amidated bile acids (Eq. (2)), in which C is the concentration, A is the amount and Vc is the volume per million hepatocytes: Unfortunately, human data to describe all of these parameters is far from complete. Therefore, the following assumptions were made to reduce the number of parameters and to include only the most relevant available human data. The first assumption is that only the first periportal hepatocytes that are exposed to the portal blood are modelled. These cells are exposed to the highest bile acid concentration, as it is well-known that this concentration decreases along the sinusoid (Groothuis et al., 1982) and as such will give the 'worst case scenario' for bile acid exposure-related toxicity. Secondly, it is assumed that the portal concentrations of bile acids do not change rapidly in response to the intake of a cholestatic drug, which is acceptable due to the large size of the bile acid pool and that strong effects on the portal bile acid concentration will only occur later on in the cholestatic process. In addition, it is assumed that the synthesis rate of bile acids is constant and that 5% of the bile acid pool is newly synthesized per day (Chiang, 2009). Regulation of bile acid synthesis could not be included in the model as no quantitative data is available for this. The fourth assumption is that the amidation rate is very high and that all bile acids are instantaneously amidated (Falany et al., 1994). Hence, uptake and synthesis rate of unamidated bile acids was assumed to be rate limiting. Moreover it is anticipated that amidation is not saturable, given the high enzymatic conversion capacity present in hepatocytes (Falany et al., 1994). This would result in the absence of unamidated bile acids intracellularly, which was indeed confirmed by analysis of human intracellular bile acid profiles (García-Cañaveras et al., 2012). Finally, it is assumed that because of their hydrophilicity and resulting poor membrane permeability, glycine and taurine amidated bile acids will not undergo uptake or efflux by passive diffusion. This leads to the following simplification of the model in which each amidated bile acid can now be described as: Several transport proteins are involved in the uptake and efflux of bile acids. To simplify the model, we only focused on the major bile acid uptake transporter NTCP and the canalicular efflux transport BSEP. We are aware that other uptake transporters such as OATP1B1, OATP2B1 and OATP1B3 are involved, as well as the canalicular efflux transporter MRP2 and the sinusoidal efflux transporters MRP3, MRP4 and OSTα/β. However, NTCP and BSEP are generally assumed to be the major players and the kinetic parameters of the other transporters are not available at this moment.
Our final assumption includes competitive inhibition of bile acids for both the uptake and efflux transporters using the adapted formula for enzymes (Asante-Appiah and Chan, 1996;Martinez-Irujo et al., 1998;Pocklington and Jeffery, 1969;Schäuble et al., 2013). In addition, the drugs we choose as model drugs, glibenclamide and cyclosporin A, are both competitive inhibitors of BSEP. For both transporters the uptake/efflux rate for each of the individual bile acids is calculated according to Eq. (8).
Vmax x is the maximum transport rate, S x is the substrate concentration, Km x is the concentration at which the transport rate is 50% of Vmax x in the absence of other bile acids and inhibitor, I is the concentration of the inhibitor, Km i is the concentration of the inhibitor at which the transport rate is 50% reduced in the absence of other BAs, and x is one of the bile acids 1-9.
1: CDCA. 2: GCDCA. 3: TCDCA. 4: DCA. 5: GDCA. 6: TDCA. 7: remaining bile acids. 8: remaining glycine amidated bile acids. 9: remaining taurine amidated bile acids. I: competitive inhibitor (in this study glibenclamide or cyclosporin A). The Km and Vmax values of the six bile acids, mentioned above, were derived from uptake and efflux studies in HEK293 cells as described below and were scaled to the human hepatocyte using the abundance of the transporters in overexpressing cells and human liver. However, as mentioned before, kinetic parameters for human transporters (NTCP and BSEP) have not been published with appropriate scaling. In literature and in our own data both the Km and Vmax were rather similar within the group of unamidated, glycine amidated and taurine amidated. Therefore, we took the average of the values obtained in each of these three groups and used these values for the unamidated, glycine amidated and taurine amidated remaining bile acid groups (see Table 3). The transport rates were then calculated using Eq. (4).
As mentioned above it was assumed that 5% of the total bile acid pool is newly synthesized per day. In addition, it was also assumed that both primary and secondary bile acids were synthesized in the hepatocytes and that bile acids were synthesized in the same ratio as they appear in blood (Scherer et al., 2009). In addition, it was assumed that all bile acids were amidated with glycine and taurine instantaneously relative to the ratio of these amidated bile acids present in the hepatocyte, i.e. 70% and 30% respectively (García-Cañaveras et al., 2012). The bile acid synthesis rate was scaled to a million human hepatocytes. Ho et al. described a bile acid pool of 0.136 mmol/kg body weight (Ho et al., 1980) in a 70 kg person with 1.5 kg of liver. This results in a bile acid synthesis rate of 0.22 pmol bile acid/min/g liver, which is 1.8 pmol bile acids/min/10 6 hepatocytes, assuming a hepatocellularity of 122 × 10 6 ·hepatocytes per gram liver (Arias, 1988;Wilson et al., 2003;Sohlenius-Sternbeck, 2006) (Arias, 1988;Sohlenius-Sternbeck, 2006;Wilson et al., 2003). For each bile acid the synthesis rate was calculated proportional to their fractional blood concentration (see Table 1). These values were all included in Eq. (3), and the intracellular concentration of each of the 9 groups of bile acids and their total concentration were simulated using Berkeley-Madonna (version 8.3.18, University of California; http://www.berkeleymadonna.com), assuming a hepatocyte volume of 0.01 ml/10 6 hepatocytes. The details of the script are given in the supplementary data.

Culture and Transduction of HEK293 Cells
HEK293 cells used for transduction with either NTCP or BSEP were cultured under standard conditions at 37°C and 5% CO 2 . Cells were cultured using Dulbecco's Modified Eagle Medium (DMEM) (ThermoFisher scientific) supplemented with 10% Fetal Calf Serum (FCS) (Greiner bio-one). HEK293 cells were seeded in 96-well plates coated with Poly-D-Lysine (Corning) and incubated for 24 h at 37°C and 5% CO 2 . Then, transduction was performed by adding 30 μl of baculovirus-containing solution, encoding NTCP, BSEP, or Enhanced Yellow Fluorescent Protein (EYFP) as a negative control. Following addition of the virus, sodium butyrate was added to a final concentration of 5.0 mM in order to enhance transcription and ultimate transporter protein expression. After another 72 h of incubation at 37°C and 5% CO 2 (uptake experiments were performed with cells expressing NTCP or EYFP and vesicles were prepared from cells transduced with BSEP or EYFP).

Production of BSEP Membrane Vesicles
Vesicles were obtained following transduction of HEK293 cells with baculovirus encoding BSEP or EYFP as a negative control. Following the incubation for 72 h at 37°C and 5% CO 2 as described above, the cells were harvested and were then lysed and fused into membrane vesicles (8) using a LV1 Low Volume High Shear microfluidizer (Microfluidics, Westwood, MA), followed by centrifugation at 4000g at 4°C for 20 min and subsequent centrifugation of the resulting supernatant at 25,000g at 4°C for 90 min in order to obtain the membrane fractions. Preparations, consisting of a 1:1 ratio of inside-out and right-side-out vesicles, were then resuspended in Tris-Sucrose (TS) buffer (10 mM Tris, 250 mM Sucrose, pH 7.4) and stored at − 80°C until further use.

Transporter Kinetics
In order to determine the Km and Vmax of bile acid transport by NTCP, HEK293 cells transduced with NTCP or EYFP (control) were used. Cells (10,000 cells/well) were cultured in poly-D-lysine coated 96well plates. After 24 h, cells and were washed with 37°C HBSS-HEPES pH 7.4 buffer, before incubation with increasing concentrations (0-50 μM) of one of the bile acids of interest in HBSS-HEPES pH 7.4. After 1 min of incubation at 37°C, the plate was inverted and immediately washed with ice-cold HBSS-HEPES pH 7.4 containing 0.5% Bovine Serum Albumin in order to bind excess bile acids. Then, the cells were washed twice with ice-cold HBSS-HEPES pH 7.4 and after removal of the washing buffer, the plates containing the cells were stored at − 20°C until subsequent LC-MS/MS analysis of intracellular bile acid concentrations, as described below.
The Km and Vmax of ATP-dependent BSEP transport of the bile acids of interest were determined using membrane vesicles containing BSEP or EYFP (control). Vesicles were thawed and transferred to a 96well plate, containing 25 μl of a mix containing TS buffer, supplemented with 10 mM MgCl 2 , 4.0 mM ATP pH 7.4 and increasing concentrations (0-50 μM) of the bile acid of interest. Following incubation of the plate at 37°C for 1 min, the plate was placed on ice and 150 μl of ice-cold TS buffer was added to stop the reaction. The samples were then transferred to a 96 well filter plate (Millipore) pre-incubated with TS buffer. The filter membrane was then washed twice with ice-cold TS buffer. Subsequently, the filters were punched out of the plate and used for LC-MS/MS mediated measurement of intravesicular bile acid concentrations as described below. In parallel experiments, vesicles derived from EYFP transduced cells were used to determine the passive transport rate of the various bile acids.

LC-MS/MS Analysis of Bile Acids
Bile acid concentrations were determined with LC-MS/MS, using an Acquity UPLC (Waters, Milford, MA, USA) coupled to a Xevo TQ-S (Waters) triple quadrupole mass spectrometer. The compounds were separated using a Zorbax Eclipse Plus C18 analytical column (Rapid Resolution HD 1.8 μm; 50 × 2.1 mm, Agilent, USA). The elution gradient was as follows: 0 min 70%A/30% B, 2-2.5 min 20%A/80% B, and 2.6-4.0 min 70%A/30% B. Solvent A consisted of 0.1% formic acid in water and Solvent B consisted of 0.1% formic acid in acetonitrile. The column temperature was set at 40°C, and the flow rate was 500 μl/min. The effluent from the UPLC was passed directly into the electrospray ion source. Negative electrospray ionization was achieved using nitrogen as a desolvation gas with the ionization voltage set at 3.0 kV. The source temperature was set at 500°C and argon was used as col-

Simulation of Intrahepatic Glibenclamide and Cyclosporin A Concentrations
Physiologically-based pharmacokinetics (PBPK) models were used to predict the portal vein and intrahepatic concentrations of glibenclamide and cyclosporin A. These concentrations were generated using Simcyp® Simulator version 13.2. For glibenclamide, a full PBPK model published by Varma et al. (2014) was used in a 7-day simulation with an oral dose of 7.5 mg given every 12 h to 50 healthy volunteers. This is the maximum recommended dose for glibenclamide. For cyclosporin A, a minimal PBPK model present within the Simcyp® Simulator was transformed into a full PBPK model by incorporating an intravenous clearance value and subsequent scaling of the volume of distribution to values reported in literature. The simulation was performed for 50 healthy volunteers taking 525 mg of cyclosporin A every 12 h for 7 days. This is the recommended maximum dose of 15 mg/kg for prophylaxis of solid organ transplant rejection during the first 1-2 weeks after transplantation for a 70 kg human being.

BSEP and NTCP Abundance Measurements
Abundance measurements were performed as described previously (Van De Steeg et al., 2013). Protein abundance of BSEP and NTCP in the plasma membrane was measured in transduced HEK293 cells (overexpressing BSEP or NTCP) and in a pooled human liver sample originating from three individual donors. Human liver samples were obtained from patients undergoing partial hepatectomy for liver metastasis of colorectal carcinoma, as described before (de Graaf et al., 2010). The absolute transporter protein abundance was corrected for the abundance of the membrane protein ATP1A1, used as an internal standard. The Vmax values of bile acid transport by BSEP and NTCP were used to calculate bile acid fluxes (according to the Eqs. (3)-(5) and (7) outlined above), which is by definition dependent on the expression of the transporters at the plasma membrane. In order to scale from in vitro to in vivo human liver, the Vmax values from the in vitro systems were expressed per pmol transport protein and subsequently scaled to the amount of transport protein present in the human liver (see Table 3). Similarly the passive diffusion rate constant was expressed per μg plasma membrane protein and subsequently scaled to the amount of plasma membrane protein in the human liver. In addition this rate constant was corrected by 1/6 assuming that only 1/6 of the plasma membrane in the liver will be on the sinusoidal side where the passive transport takes place.

Kinetic Parameters of Bile Acid Transport by NTCP and BSEP
The results of the uptake of bile acids by HEK293 cells overexpressing NTCP are shown in Fig. 1. Km values for these 6 bile acids were rather similar and ranged from 4.9-10 μM. However, there was a clear difference in Vmax between the unconjugated and conjugated bile acids, with the unconjugated showing a lower maximal transport rate (40-50 pmol/mg protein/min) than the conjugated (1200-1700 pmol/ mg protein/min).
For BSEP transport studies, we used vesicles isolated from HEK293 cells over-expressing BSEP and compared the data to those obtained with control vesicles containing EYFP. The results from these studies are presented in Fig. 2. Whereas TDCA, GCDCA and TCDCA showed similar Km values between 4.1 and 7.4 μM, we detected a relatively lower affinity of GDCA (Km = 19 μM). We did not determine Michaelis-Menten kinetics for the two unconjugated bile acids included in the present study, as one of the key assumptions was instant amidation of unconjugated bile acids upon uptake by the hepatocyte, thus excluding their biliary excretion.

Simulated Intrahepatic Concentrations of Glibenclamide and Cyclosporin A
In order to simulate the effect of cholestatic drugs in the model, we selected two known cholestatic compounds (glibenclamide and cyclosporin A) and we used the Simcyp simulator to predict their intrahepatic concentrations upon reaching steady-state levels under commonly applied dosing regimen. In short, the predicted intrahepatic concentration for glibenclamide at a dosing regimen of 7.5 mg twice daily was 0,02 uM. The total predicted steady-state concentration in the liver after cyclosporin A administration (525 mg bid) was 6.6 μM. As no good estimate of intracellular unbound concentrations was available for CsA, we ran subsequent simulations in the biokinetic bile acid handling model for the worst case scenario, i.e. assuming all cyclosporin A within the hepatocyte to be unbound and available for BSEP inhibition.

Abundance Measurement of NTCP and BSEP
In the HEK293 cells transduced with NTCP, the abundance of NTCP was 156 fmol/μg plasma membrane protein. The abundance of transport protein in vesicles of BSEP-overexpressing HEK293 cells was 59.4 fmol/μg plasma membrane protein. In the human liver sample the abundances were 5.9 fmol/μg plasma membrane protein and 6.7 fmol/ μg plasma membrane protein, respectively (Table 2).

Biokinetic Mechanistic Model
The biokinetic model used after simplification as described in the methods section, is represented in Fig. 3.In order to parameterize our biokinetic model with the results from the HEK293 cellular overexpression studies, a correction is needed for the different abundances of the transporters in the cell systems and the human liver. This was achieved by scaling the NTCP and BSEP expression in the different systems to an average expression in human liver (Table 2). For the abundance of NTCP in human liver 5.9 fmol/μg plasma membrane protein and of BSEP 6.7 fmol/μg plasma membrane protein was found (Table 2), resulting in scaling factors of 4.01 × 10 − 7 and 7.93 × 10 − 7 . The scaled abundance data were used to calculate the Vmax in human liver for the individual bile acid groups. For this purpose, Vmax values from the in vitro systems were expressed per pmol transport protein and subsequently scaled to the amount of transport protein present in the human liver. HEK293 cells overexpressing NTCP contained 1.9 μg plasma membrane protein/mg total protein. HEK293 vesicles overexpressing BSEP contained 60 μg crude membrane protein/mg total protein, resulting in 32 μg plasma membrane protein/mg total protein (unpublished data). For the BSEP expressing membrane vesicles it was assumed that 50% were oriented inside-out and 50% rightsideout. This implies that only 50% of the expressed BSEP protein is functionally involved in bile acid uptake by the vesicles. The pooled liver contained 4.2 μg plasma membrane protein per million cells assuming 122 × 10 6 ·hepatocytes per gram liver (Arias, 1988;Wilson et al., 2003;Sohlenius-Sternbeck, 2006). Table 3 shows the Km and Vmax values for human liver used in the biokinetic model (obtained in this study) as well as the portal concentrations of the respective bile acids used for modeling (Scherer et al., 2009;Vatakuti et al., 2016a).
Passive transport was determined by the rate observed in the EYFPtransduced HEK293 cells, assuming that this transport was passive. A slight overestimation of transport rate cannot be excluded due to nonspecific binding. However this will not be of major influence on the model outcome as the passive transport rate is low, between 0.37 and 0.53 pmol/min/10 6 hepatocytes/μM, as described in Table 3.

Biokinetic Mechanistic Model for Bile Acid Handling
All values used to parameterize the biokinetic model are given in Table 3. When the values for Km and Vmax of the nine bile acids as well as the data for bile acid synthesis were incorporated, the model did not reach a physiological steady-state concentration of bile acids in the hepatocyte as depicted in Fig. 4 A, indicating missing or inadequate data. As can be seen in Fig. 4B, the efflux system is saturated rapidly due to the fast increase of bile acid concentration and operates at maximum velocity, leading to a continuous increase in intracellular bile acid concentrations in time. Apparently, the uptake rate of bile acids exceeds the excretion rate, even under Vmax conditions of BSEP. However, by additional empirical adjustment of BSEP abundance, based on the assumption that our model likely reflects an   (Arias, 1988;Wilson et al., 2003;Sohlenius-Sternbeck, 2006).

S. Notenboom et al.
European Journal of Pharmaceutical Sciences 115 (2018) 175-184 underestimation of the BSEP abundance, we found that increasing the abundance by 2.4 times, resulted in near physiological intracellular bile acid concentrations after exposure to physiological portal vein bile acid levels (Fig. 4C, D).

Effect of the Presence of Cholestatic Drugs on the Intrahepatic Bile Acid Concentration
Upon intrahepatic exposure of 6.6 μM cyclosporin A, an increase in intracellular bile acid concentration is predicted. This results in a new steady state with a 1.7 fold increase in hepatocyte bile acid concentration (Fig. 5A, B). However, simulation of exposure to the therapeutic liver concentration of glibenclamide (0.02 μM) did not result in an increase in bile acid concentration (Fig. 5C). We used the model to find a concentration of glibenclamide that would results in a similar increase in bile acid concentration as 6.6 μM cyclosporin A. This appeared to be a concentration of 20 μM, which is 1000-fold higher than the predicted therapeutic glibenclamide concentration (Fig. 5D).

Discussion
We developed a mechanistic biokinetic model of hepatocyte bile acid handling based on human in vitro data (Fig. 3), which simulates intrahepatic bile acid concentrations in the absence and presence of cholestatic drugs. This model can be used for prediction of cholestatic potential of drugs during drug development. We measured and included kinetic data of two selected bile acids and their glycine and taurine conjugates for the two most important human bile acid transporters (NTCP and BSEP) and the remaining bile acids were taken as one group as well as their glycine and taurine conjugates. The abundance of NTCP and BSEP was measured in our in vitro HEK293 cell models and in human liver for scaling purposes, and we used data obtained from literature on bile acid synthesis and amidation. (Table 1).
An essential step for development of the model described here was the generation of kinetic data for the interaction of bile acids with NTCP and BSEP. The only literature data available for BSEP are Km values for TCDCA, TDCA, and GCDCA, which range between 4.2 and 12.5, 34.4, and 1.7-7.5, respectively (Hayashi et al., 2005;Kis et al., 2009;Yamaguchi et al., 2009). The Km values obtained for TCDCA and GCDCA in the present study are similar to these previously reported values. However, the Km we determined for TDCA was lower than published by Kis et al., which could be accounted for by the use of different model systems. For the other three bile acids in our model, we have generated kinetic data on BSEP-mediated transport. To our knowledge, this is the first time that transporter kinetics have been determined for the interaction of such an extensive set of bile acids (including DCA and CDCA) with the two major bile acid transporters, BSEP and NTCP. Interestingly, the four conjugated bile acids studied here show rather similar affinities and maximum transport rates for each of the transporters, but these values are clearly different from those of the two unconjugated bile acids. They appear to be transported by both transporters at a much lower rate than their conjugated counterparts and are mainly taken up by passive diffusion as measured by the uptake rate in mock-transduced HEK293 cells, lacking a bile acid  uptake transporter. This can be explained by their higher lipophilicity (LogP > 2.5). Therefore, we included parameters for passive uptake in our model (Table 3). We assumed that bile acid amidation is very rapid and not saturated under physiological conditions with a Km around 1 mM for taurine conjugation and 6 mM for glycine conjugation (Falany et al., 1994). Both values are much higher than the intracellular bile acid concentration; therefore excretion of unconjugated bile acids from the liver will play only a minor role and can be neglected in our model. The unexpected accumulation of bile acids in the simulations where the experimentally found BSEP abundance was used could be caused by an overestimation of uptake rate or underestimation of excretion rate. Scaling of BSEP abundance was adapted so that it resulted in near physiological intracellular bile acid concentrations under normal exposure to portal vein bile acid levels. The relatively small scaling factor of 2.4 needed to reach a physiological total bile acid concentration in the hepatocyte shows that the model is very sensitive to BSEP abundance, as was also predicted by the DILIsym simulations (Woodhead et al., 2014). As we only measured BSEP abundance in a pooled sample of 3 human livers, no conclusion can be drawn on the variation of its expression in human livers. However, based on the data of Ji et al., it is expected that the expression of BSEP is subject to > 2fold variation. They analyzed 15 human livers and found BSEP expressions between 3 and 7 fmol/ug membrane protein (Li et al., 2009). The expression of 7 fmol/ug membrane protein found in our pooled liver sample is within this range. In addition to variations in abundance, variability in BSEP functional activity could contribute to differences in Fig. 4. The predicted intracellular concentrations (A) and canalicular efflux rates (B) of bile acids in the human hepatocyte following exposure to 60 μM bile acids on the portal side. The black dotted line in 3B represents the uptake rate of total bile acids (TBA) by NTCP, showing that uptake > efflux (B). After fitting the model to intracellular bile acid concentrations within the physiological range as measured by Starokozhko et al. (unpublished data), the model shows the adjusted predicted intracellular concentrations (C) and canalicular efflux rates (D) of bile acids in the human hepatocyte following 60 μM bile acids exposure on the portal side. the population as a result of polymorphisms  or subcellular distribution in membrane vesicles (Kubitz et al., 2004). Further studies with more human liver samples and more sophisticated separation of canalicular membranes from other intracellular membranes/vesicles may provide more insight into this scaling factor. The high sensitivity of intracellular bile acid concentration to BSEP activity as shown in our model and in the DILIsym model can explain the large variation in sensitivity for cholestatic drug effects in the human population. The value found for NTCP of 5.9 fmol/μg plasma membrane protein might have been an overestimation of active transporter as it is higher than found by others (Qiu et al., 2013;Wang et al., 2015) and it has been reported that the transporter can be also found in submembrane vesicles (Sommerfeld et al., 2015;Webster et al., 2002). The role of OATPs in the uptake of bile acids was assumed to be low, but this should be incorporated in a future version of the model. A more likely explanation of the accumulation of bile acids might be the presence of other biliary efflux transporters, such as MRP2 at the canalicular membrane. However, as MRP2 has always been described to transport sulfate and glucuronide conjugates, which constitute only a very small fraction of the bile acid pool, its contribution to biliary bile acid excretion rate might be limited. Nevertheless this fraction includes sulfated LCA, which is supposed to be toxic to the hepatocyte. Also MRP3, MRP4 and OSTα/β on the basolateral membrane might be involved in sinusoidal excretion, especially under cholestatic conditions, as the expression of MRP3 and 4 are low under normal conditions but are upregulated during cholestasis (Gradhand et al., 2008;Vinken et al., 2013). Adding these efflux transporters to the model might improve the simulations, but quantitative data is lacking to date. In addition, our model does not take into account sinusoidal excretion and reuptake of bile acids in the hepatocytes downstream along the sinusoids. As this sinusoidal excretion is low, additional availability of bile acids for uptake is also limited, which results in lower portal bile acid concentrations down the sinusoids. It should also be noticed that the BSEP abundance measured in this study is an average of all hepatocytes, while it was shown that bile acid excretion is lower in perivenous cells than periportal cells if exposed to the same bile acid concentration. (Groothuis et al., 1982).
As already indicated above, bile acid amidation was considered to be immediate and complete, which is supported by the high Km values (1-6 mM) found for this metabolism (Falany et al., 1994). However it cannot be excluded that under cholestatic conditions amidation is saturated or undergoes product inhibition, which should be subject of future investigations.
The synthesis rate of new bile acids was based on the general assumption that synthesis equals bile acid loss, which is estimated to be 5% of the bile acid pool per day. Although it is generally known that proteins involved in bile acid synthesis and transport are regulated via FXR, DILIsym results indicated that this regulatory pathway does not have a substantial effect on the intracellular concentration of bile acids (Woodhead et al., 2014). Therefore, and because experimental quantitative data are lacking, we did not incorporate any regulation in our model.
Subsequently, we simulated the effect of the presence of two known cholestatic drugs, cyclosporin A and glibenclamide, which resulted in an increase in intracellular bile acid concentrations. For cyclosporin A we used the maximal intrahepatic concentration (6.6 μM) after an oral dose of 525 mg twice daily (15 mg/kg/day), as estimated with the Simcyp simulator. This induced a 1.75-fold increase of the bile acid concentration to 140 μM. The question remains how much bile acid concentrations should rise to result in cholestatic liver injury. Recent studies with rat precision cut liver slices (PCLS) (Starokozhko et al., 2017) showed that a 50% decrease in viability (measured by ATP content) was observed at a cyclosporin A concentration that raised the bile acid content in the slices by a factor of 2.5. Studies with human PCLS are ongoing.
For glibenclamide an intrahepatic concentration of 0.02 μM was predicted using the Simcyp simulator for a therapeutic dose of 15 mg/ day. In our biokinetc model this concentration did not result in changes in bile acid concentrations, which can be explained by the high Km of glibenclamide for BSEP (27.5 μM) (Byrne et al., 2002). This is in line with clinical findings that at therapeutic doses of 7.5 mg bid glibenclamide-induced cholestasis is an idiosyncratic phenomenon, occurring only in a minority of patients. At a 1000-fold higher concentration, an increase in bile acid levels in the hepatocyte was predicted, which is well in line with the DILIsym data on rats with very high plasma concentrations of glibenclamide (> 20 μM) that displayed increased plasma concentrations of bile acids (Woodhead et al., 2014). Moreover, a 2.4-fold increase in liver bile acid concentration was found in rat PCLS exposed to 180 μM (=89 μg/ml) glibenclamide, but not at lower concentrations (Starokozhko et al., 2017). These data indicate that glibenclamide-induced cholestasis can be due to BSEP inhibition only in patients with a genetic variant or low abundance of BSEP. Alternatively, it cannot be excluded that glibenclamide induces cholestasis by an as yet unknown mechanism (Starokozhko et al., 2017).
In summary, the mechanistic biokinetic model presented here is able to provide a first indication on whether drugs at therapeutic concentrations might cause cholestasis by inhibiting BSEP. Next to an inhibitory constant of a drug for BSEP, a PBPK model to simulate its intrahepatic concentration will be needed to predict whether this simulated cholestatic concentration will actually occur at the intended therapeutic plasma concentrations. This would enable performing initial risk assessments based on human in vitro data only, thereby serving as a replacement of animal experiments.
In the future our model can be optimized and refined by adding more data on specific bile acids, especially the sulfation and glucuronidation rates of the toxic lithocholic acid. In addition more data are needed on bile acid transport by OATPs, MRP3, MRP4 and OST α/β, saturation of amidation, regulation of transporter expression and synthesis rate, as well as protein binding of the bile acids and drugs. Realizing all these shortcomings, the predictive value of the model for cyclosporine A and glibenclamide appeared remarkably well, indicating that the most important parameters have been included.