Selection and classification of drugs

The criteria for the inclusion of a drug in our panel were (a) causally linked to the development of DILI by the U.S. Food and Drug Administration (FDA) (Chen et al. 2016), (b) oral administration in clinical practice, (c) able to inhibit bile acid efflux with the half inhibitory concentration (IC50) for bile acid efflux available from an assay with human suspension-cultured hepatocytes and the IC50 being < 100 µM (Zhang et al. 2016), (d) the reported DILI being not immune-mediated, and (e) two or more in vivo pharmacokinetic studies available in the literature to validate the drug PBK model-based predictions for plasma concentrations. The maximal concentration used for the IC50 determination by Zhang et al. (2016) was 100 µM. This resulted in a final inclusion of 18 drugs, i.e., atorvastatin, bicalutamide, bosentan, chlorpromazine, cyclosporine, deferasirox, fluoxetine, flutamide, glimepiride, haloperidol, lovastatin, ketoconazole, pioglitazone, ritonavir, rosiglitazone, saquinavir, trazodone, and troglitazone. Flutamide, lovastatin, and saquinavir were excluded from further simulations based on the results from the generic PBK model, as described in “Glycochenodeoxycholic acid PBK model” section. The drugs resulted in different types of DILI. The LiverTox® database classified chlorpromazine, cyclosporine, and ritonavir as cholestatic DILI, and the remaining drugs were classified as hepatocellular/mixed DILI (http://LiverTox.nih.gov; last accessed on the 24th of August 2023). For atorvastatin, the LiverTox® database reported two cases of mixed, one of cholestatic and one of hepatocellular DILI. Glimepiride and haloperidol were not present in the LiverTox® database. The cholestasis incidence of the 15 remaining drugs was evaluated based on the European database of suspected adverse drug reaction (ADR) reports (www.adrreports.eu; last accessed on the 24th of August 2023). The following adverse reactions were considered cholestatic: cholestasis, cholestatic liver injury, cholestatic jaundice, and cholestatic hepatitis. The incidence of cholestasis was classified by us, using the European ADR database as follows: common (> 0.5% of ADR cholestatic), rare (0.3–0.5% of ADR cholestatic), and no (< 0.3% of ADR cholestatic). The threshold for common incidence was set to 0.5% of ADR to ensure that the drugs chlorpromazine, cyclosporine, and ritonavir, which were identified as cholestatic in the LiverTox® database, were classified as common causes of cholestasis. The 0.3% threshold was set artificially to account for background cholestasis incidence. Additional drugs classified as common for induction of cholestasis according to our classification were: bosentan, ketoconazole, atorvastatin, and glimepiride. Riede et al. (2017) reviewed cholestasis incidence of several drugs based on cohort and retrospective studies. In line with our classification, Riede et al. (2017) classified cyclosporine cholestasis incidence as common and rosiglitazone as not cholestatic, but ketoconazole and atorvastatin cholestasis incidences were classified as rare in contrast to our classification system. No explanation was found for the discrepancy between atorvastatin and ketoconazole cholestasis incidence classification  by these authors and the ADR database, hence, we used our classification system based on the European ADR database. Bosentan was classified as rare or common depending on the dose in the review of Riede et al. (2017). At therapeutic dose level, i.e., 125 mg twice a day, bosentan was considered a rare cause of cholestasis. Since the drugs were evaluated at their therapeutic dose level in the current study, bosentan-induced cholestasis incidence was classified as rare. Troglitazone was banned from the market and therefore not in the European ADR database, and also not in the review by Riede et al. (2017). A review of cohort studies indicated that troglitazone causes hepatocellular liver injury, with rare instances of mixed or cholestatic liver injury (Chojkier 2005). Troglitazone cholestasis incidence was therefore classified as rare. Maximal prescribed daily dosage was used for simulations and obtained from the supplier’s prescription information.

Generic PBK models for drugs

A generic PBK model was used to predict the hepatic concentrations of the selected drugs at therapeutic dose level and above. These concentrations were subsequently used to predict the inhibitory effect on hepatic bile acid efflux and resulting bile acid accumulation using a coupled bile acid PBK model (see “Glycochenodeoxycholic acid PBK model” section).

The generic drug PBK models were adapted from Punt et al. (2022). Briefly, the PBK models consisted of compartments for lung, adipose, bone, brain, heart, intestine, liver, kidney, muscle, skin, spleen, and arterial and venous blood. Different compared to Punt et al. (2022), a blood:plasma ratio of 0.55 was used for acidic compounds (1-hematocrit), and 1 for neutral or basic compounds (Cubitt et al. 2009). Physicochemical properties (pKa, logP, logD, topological surface area, and molecular weight) of the drugs were predicted using Chemicalize, https://chemicalize.com/ developed by ChemAxon (http://www.chemaxon.com). The physicochemical properties were subsequently used to predict tissue:plasma partition coefficients (Berezhkovskiy 2004; Rodgers and Rowland 2006), absorption rate constants, and fractions absorbed (Hou et al. 2004).

As part of this study, several in vitro and in silico methods were evaluated to derive the tissue:plasma partition coefficients, hepatic intrinsic clearance, and fraction unbound in plasma (Fup) (Table 1). In more detail, for the tissue:plasma partition coefficients, values derived using the in silico methods of Rodgers and Rowland (2006) and Berezhkovskiy (2004) were compared. For the fraction unbound (Fup), both the in silico method of Lobell and Sivarajah (2003) and in vitro rapid equilibrium dialysis data using human plasma were evaluated. Clearance data were derived from in vitro hepatocyte studies or the pkCSM in silico tool (Pires et al. 2015). It should be noted that in vitro the hepatic intrinsic clearance was measured, while pkCSM predicted the total clearance, i.e., a combination of hepatic and renal clearance. Where possible, in vitro intrinsic hepatic clearance (CLint) and Fup data were obtained from the high-throughput toxicokinetic (httk) database (Pearce et al. 2017); alternatively, a literature search in Scopus was conducted to obtain the in vitro CLint (Supplementary Table S1). In vitro CLint was determined using fresh or cryopreserved primary human hepatocytes (PHH) cultured in a monolayer or suspension. The CLint was determined based on substrate depletion or metabolite formation. The incubation medium did not contain any serum and was optimized for hepatocyte function, maintaining physiological temperature and pH.

Table 1 Input parameters for the generic PBK model

Corrections for non-specific binding of the compounds to the hepatocytes in vitro were applied based on the calculation method of Kilford et al. (2008). The in vitro and in silico clearance data were scaled to the in vivo situation based on a hepatocellularity of 117.5 × 106 hepatocytes per gram liver (Barter et al. 2007) and a liver weight of 1470 g, or 70 kg body weight, respectively, see Eqs. 1 and 2

$$}_} ,}\;}}} = }_} ,}\;}}} \times } \times } \times 60 \times 10^ ,$$

(1)

$$}_},}\;}}} = }_},}\;}}} \times } \times 60 \times 10^ ,$$

(2)

where in Eq. 1, CLint,in vivo is the intrinsic clearance in vivo in L min−1 entire liver−1, CLint,in vitro the hepatic intrinsic clearance in vitro in µL min−1 10–6 hepatocytes, Hep the hepatocellularity in 106 hepatocytes/g liver, and Vli the weight of the liver in grams. A factor of 10–6 L µL−1 and 60 min h−1 was applied to convert the CLint to units applicable to the PBK model. Equation 2 describes the in silico-to-in vivo scaling of the total clearance (CLtot). CLtot,in vivo is the total clearance in vivo L min−1 entire liver−1, CLtot,in silico the total clearance as predicted by pkCSM in mL min−1 kg body weight−1, and BW the body weight in kg. Here, factors of 10–3 L mL−1 and 60 min h−1 were applied to convert the CLtot to units applicable to the PBK model.

A literature search was performed to compile a dataset on human in vivo drug plasma peak concentrations (Cmax) after a single oral dose of the selected drugs. The following keywords were used in Scopus: ((TITLE (“drug name”) AND ALL (bioavailability OR pharmacokinetics OR kinetics)) AND (( human OR man OR volunteer OR subject)) AND (Cmax OR “c max” OR “maximal concentration” OR “maximum concentration” OR “peak concentration”)). The studies that were identified for each drug were subsequently filtered to exclude (1) results obtained for specific patient groups like patients with renal impairment or gastric by-pass, (2) studies with children, and (3) studies using slow or extended-release formulations. At least two Cmax values per drug were collated from different doses and/or studies. The dose (normalized to bodyweight), Cmax and references were gathered in the file “Invivo.xlsx” in the Github repository.Footnote 1 The final dataset included 115 studies and 179 sets of dose and corresponding Cmax values for 18 drugs. A default bodyweight of 70 kg was assumed if the bodyweight was not reported in the study. Cmax values were averaged if experiments were performed under both fasted and fed conditions. Most studies reported peak concentrations in plasma, but if blood concentrations were reported, they were converted to plasma concentrations based on the blood:plasma ratio (see Table 1). The studies included a variable number of human volunteers per group with an arithmetic mean of 20 (range: min. 5–max. 89).

The predicted plasma Cmax was compared with the observed plasma Cmax as obtained from the meta-analysis. The ratio predicted:observed Cmax was calculated for each study and/or dose. This resulted in a number of ratios per compound. The median of this ratio was calculated per compound, and for further simulations, the combination of input parameters that gave a median ratio predicted:observed Cmax closest to 1 was selected. Only the drugs of with a median ratio predicted:observed Cmax within tenfold were used for further analysis (n = 15). An overview of the metabolizing enzymes or transporters involved in the kinetics of the drugs was made to find explanations for over- or underpredictions. The information about the involved enzymes and transporters was obtained from literature (Wishart et al. 2018; Elsby et al. 2012; Cockshott 2004; Hebert 1997; Klatt et al. 2011; Treiber et al. 2007).

Glycochenodeoxycholic acid PBK model

A PBK model describing synthesis, circulation and excretion of bile acids in healthy individuals was based on our previous work (de Bruijn et al. 2022, 2023). The conceptual model is presented in Fig. 1. The enterohepatic circulation was modeled as a circulation of GCDCA between the liver (extracellular and intracellular), gallbladder, and intestine. The intestinal uptake and the hepatic uptake and efflux were described using carrier-mediated transport processes, i.e., ASBT, NTCP, or BSEP-mediated, respectively. The NTCP-mediated hepatic uptake of GCDCA was modeled permeability-limited as described in our previous work (de Bruijn et al. 2023). The kinetic parameters for ASBT-mediated transport were obtained using Caco-2 cells cultured on permeable cell culture inserts and scaled from the in vitro to in vivo situation as described in our previous work (de Bruijn et al. 2023). GCDCA de novo synthesis in the liver was set equal to its excretion via the feces. GCDCA was actively transported from the liver to the common bile duct by BSEP following Michaelis–Menten kinetics. The BSEP-mediated efflux of GCDCA was described by Eq. 3

$$E = \frac}}} \times \left[ }} \right]}}}}} + \left[ }} \right]}},$$

(3)

where E is the BSEP-mediated efflux in µmol/h, Vmax is the maximum efflux rate of GCDCA in blood in µmol/entire liver/h, [Cliveriw] the free concentration of bile acids in intracellular water in liver in µmol/L and Km,BSEP the Michaelis–Menten constant in µmol/L for BSEP-mediated GCDCA efflux. The Vmax and Km for BSEP-mediated transport of GCDCA were obtained from a vesicular transport assay in a baculovirus-infected Sf9 system (Kis et al. 2009).

Fig. 1

Conceptual model for the PBK modeling of bile acid homeostasis and the influence on this homeostasis by drugs. Conceptual model for GCDCA PBK model was taken from de Bruijn et al. (2023), and conceptual model for generic drugs PBK models was taken from Punt et al. (2022). CLint: intrinsic clearance; Fup: fraction unbound plasma; GCDCA: glycochenodeoxycholic acid; GFR: glomerular filtration rate

The differential model equations were encoded and solved using the deSolve package version 1.32 in R version 4.1.0 (Soetaert and Petzoldt 2010; R Core Team 2022). The model code can be found in the Github repository (see Footnote 1).

Inhibitory effect of drugs on hepatic bile acid efflux

The hepatic-free concentration of the within tenfold predicted drugs at the maximal prescribed daily therapeutic dosage was used to simulate their inhibitory effect on hepatic GCDCA efflux. The conceptual PBK model for the drugs was combined with the PBK model for GCDCA as displayed in Fig. 1. The concentrations required to reduce hepatic bile acid efflux by 50% (IC50) were derived from a study using suspension-cultured primary human hepatocytes (PHHs) (Zhang et al. 2016) and corrected for in vitro non-specific binding (Kilford et al. 2008). The IC50 values obtained using PHHs were for all drugs except pioglitazone lower than the IC50 values obtained using BSEP-transfected membrane vesicles (Supplementary Table S2). Therefore, the results obtained from suspension-cultured PHHs were used for further simulations as a worst-case estimate. Membrane vesicles and suspension-cultured PHHs provide different insights in hepatic bile acid efflux. PHHs are known to have a physiologically relevant expression of several transporters and metabolizing enzymes, rather than the exclusive BSEP expression in membrane vesicles. Hence, in contrast to membrane vesicles, the drug-induced inhibitory effect on bile acid efflux from primary suspension-cultured human hepatocytes is not necessarily caused by an exclusive inhibition of BSEP-mediated transport, but could also be caused by inhibition of the uptake, basolateral efflux by MRP3/4 or the conjugation process. Nevertheless, the net drug-induced bile acid efflux inhibition was incorporated in the equation for BSEP-mediated efflux. The PHHs were commercially obtained and pooled from 10 donors, including 5 males and 5 males aged 17–65. Cells were cultured in suspension at a density of 0.25 × 106 cells/mL in the presence of 10 µM cholic acid and various drug concentrations. Upon incubation, the PHHs conjugated the cholic acid with a glycine group resulting in glycocholic acid. After incubation, intra- and extracellular bile acid concentrations were quantified using LC–MS/MS. In the absence of specific data on the effects of all the selected drugs on bile acid efflux for GCDCA, glycocholic acid efflux was used as a surrogate for all drugs (Zhang et al. 2016). Glycocholic acid showed greater inhibition (lower IC50) in assays with PHHs than GCDCA and was considered as a worst-case estimate of GCDCA efflux (Chothe et al. 2021; Yucha et al. 2017). Competitive inhibition was assumed, as this is the typical mode of drug-transporter-inhibition (Kenna et al. 2018). The following formula to calculate Ki applies for competitive inhibitors like the drugs considered here (Eq. 4) (Yung-Chi and Prusoff 1973):

$$K_}} = \frac}_ }}}} }}}},$$

(4)

where Ki is the inhibitory constant in µM, IC50 the half maximum inhibitory concentration of the drug in µmol/L, [S] the substrate concentration in µM, and Km the Michaelis–Menten constant in µM.

Subsequently, the inhibitory effect of the drugs on GCDCA efflux was incorporated in the PBK model equation describing BSEP-mediated efflux. In line with competitive inhibition, the Km,BSEP in the corresponding Michaelis–Menten reaction (Eq. 3) was modified to the apparent Km (Km,BSEP,app) as follows (Eq. 5):

$$K_},},}}} = K_},}}} \left( } }} }}} \right),$$

(5)

where Km,BSEP,app is the apparent Michaelis–Menten constant in µM, [I] the unbound hepatic concentration of the inhibitor (= drug) in µM and Ki the inhibitory constant in µM.

The dose-metric used to evaluate intracellular GCDCA accumulation was the area under the curve (AUC), because it has been acknowledged AUC is the most relevant for endpoints that are influenced by total dose over time resulting in an accumulation (Rietjens et al. 2019).

Simulating sensitive individuals

The developed PBK approach was also employed to evaluate drug effects on intrahepatic GCDCA accumulation in sensitive individuals. For these studies, cyclosporine was selected as the model drug, because the maximal prescribed daily therapeutic dose resulted in intrahepatic concentrations equivalent to the Ki for bile acid efflux inhibition, facilitating detection of changes in the bile acid accumulation. In our previous work, we established that an over 1.5-fold increased total bile acid pool size posed an individual at risk for intrahepatic bile acid accumulation as a result of BSEP-inhibition (de Bruijn et al. 2022). Furthermore, a low BSEP abundance was identified as a potential risk factor for the development of cholestasis. Therefore, as an example, the effects of cyclosporine on intrahepatic accumulation were simulated for (a) a reference individual, (b) an individual with a 1.5-fold increased total bile acid pool size compared to the reference individual, (c) an individual with low hepatic BSEP abundance, or (d) an individual with an increased pool size and a low BSEP abundance. The low and reference BSEP abundances were derived from a meta-analysis of hepatic transporter abundances in healthy Caucasians (Burt et al. 2016). The in vitro to in vivo extrapolation of Vmax was based on the BSEP abundance. A lower BSEP abundance thus resulted in a lower in vivo Vmax. The reference individual had a BSEP protein abundance of 0.84 pmol BSEP protein per a million hepatocytes. Low BSEP abundance was set to the reported mean minus three times the standard deviation and amounted to 0.23 pmol BSEP protein per million hepatocytes.

Sensitivity analysis

To assess the influence of the parameters on the model outcome, a sensitivity analysis was performed for the plasma Cmax of the drugs and the intrahepatic GCDCA levels. The drug’s doses were set to 1 mg/kg body weight. For plasma Cmax, all potential combinations of input parameters were evaluated. The normalized sensitivity coefficients (NSC) for hepatic GCDCA levels were calculated using the combination of input parameters giving the drug’s Cmax in closest agreement with the in vivo data (Supplementary Table S3). Based on the method reported by Evans and Andersen (2000), the normalized sensitivity coefficients (NSCs) for the model parameters were calculated as follows:

$$} = \frac - C}} - P}} \times P/C,$$

(6)

where C indicates the initial value of the model output, and C′ indicates the modified value of the model output resulting from an increase in the parameter value. P indicates the initial parameter value and P′ indicates the modified parameter value after a 5% increase of its value, keeping all other parameters at their original value.