In Silico Identification and in Vitro Validation of Potential Cholestatic

Apr 8, 2014 - ABSTRACT: Drug-induced cholestasis is a frequently observed side effect of drugs and is often caused by an unexpected interaction with t...
0 downloads 0 Views 4MB Size
Download PDF
Article pubs.acs.org/crt

In Silico Identification and in Vitro Validation of Potential Cholestatic Compounds through 3D Ligand-Based Pharmacophore Modeling of BSEP Inhibitors Tina Ritschel,† Susanne M. A. Hermans,†,‡ Marieke Schreurs,‡ Jeroen J. M. W. van den Heuvel,‡ Jan B. Koenderink,‡ Rick Greupink,‡ and Frans G. M. Russel*,‡ †

Computational Discovery and Design (CDD) Group, Centre for Molecular and Biomolecular Informatics (CMBI), Radboud university medical center, P.O. Box 9101, 6500 HB Nijmegen, The Netherlands ‡ Department of Pharmacology and Toxicology, Radboud Institute for Molecular Life Sciences (RIMLS), Radboud university medical center, P.O. Box 9101, 6500 HB Nijmegen, The Netherlands S Supporting Information *

ABSTRACT: Drug-induced cholestasis is a frequently observed side effect of drugs and is often caused by an unexpected interaction with the bile salt export pump (BSEP/ABCB11). BSEP is the key membrane transporter responsible for the transport of bile acids from hepatocytes into bile. Here, we developed a pharmacophore model that describes the molecular features of compounds associated with BSEP inhibitory activity. To generate input and validation data sets, in vitro experiments with membrane vesicles overexpressing human BSEP were used to assess the effect of compounds (50 μM) on BSEP-mediated 3 H-taurocholic acid transport. The model contains two hydrogen bond acceptor/anionic features, two hydrogen bond acceptor vector features, four hydrophobic/aromatic features, and exclusion volumes. The pharmacophore was validated against a set of 59 compounds, including registered drugs. The model recognized 9 out of 12 inhibitors (75%), which could not be identified based on general parameters, such as molecular weight or SlogP, alone. Finally, the model was used to screen a virtual compound database. A number of compounds found via virtual screening were tested and displayed statistically significant BSEP inhibition, ranging from 13 ± 1% to 67 ± 7% of control (P < 0.05). In conclusion, we developed and validated a pharmacophore model that describes molecular features found in BSEP inhibitors. The model may be used as an in silico screening tool to identify potentially harmful drug candidates at an early stage in drug development.

■

INTRODUCTION Bile acid homeostasis is an important task of the liver, which is regulated by the interplay between enzymatic and membrane transport processes. Bile acids are secreted by active transport across the canalicular membrane of the hepatocytes into intrahepatic bile canaliculi. Once excreted into the gut lumen, more than 90% of the bile acids is reabsorbed into the blood and is taken up again into the hepatocytes by influx transporters in the sinusoidal membrane. This enterohepatic circulation ensures the continuous recycling of bile acids, thus maintaining bile acid homeostasis.1,2 The bile salt export pump, BSEP/ABCB11, is located in the canalicular membrane of hepatocytes. It has 12 transmembrane helices and two Walker domains for ATP binding.1 BSEP is the most important transporter mediating the efflux of bile acids, in particular conjugated monovalent bile acids. Mutations in the BSEP gene can result in impaired function of the protein and lead to an obstructed bile flow, cholestasis, and the accumulation of bile acids within hepatocytes can ultimately lead to severe hepatotoxicity. Bile acid retention can be also © 2014 American Chemical Society

caused by compounds that inhibit BSEP activity. Several therapeutic drugs are known to induce cholestasis, such as rifampicin, cyclosporine A, and troglitazone, which was withdrawn from the market because of hepatotoxicity.1,3,4 Because drug-induced cholestatic liver injury remains an important reason for late stage attrition, predictive methods are required to identify compounds that can inhibit BSEP at an early stage in drug development.5,6 This could be achieved by using molecular modeling to predict compounds that can induce cholestatic effects. Previous in silico models for the identification of BSEP inhibitors are based on (physico-) chemical molecular properties or 2D chemical structures.7,8 A 3D structural model that describes the interaction of smallmolecule compounds with this transporter may help to refine the predictions made with such models. However, to the best of our knowledge, a 3D pharmacophore model of BSEP inhibition is currently not available. The absence of 3D models is, in part, Received: February 6, 2014 Published: April 8, 2014 873

dx.doi.org/10.1021/tx5000393 | Chem. Res. Toxicol. 2014, 27, 873−881

Chemical Research in Toxicology

Article

in ice-cold homogenization buffer (0.5 mM Na2PO4 and 0.1 mM EDTA, pH 7.4) supplemented with protease inhibitors (100 μM phenylmethylsulfonyl fluoride, 5 μg/mL aprotinin, 5 μg/mL leupeptin, 1 μM pepstatin, and 1 μM E-64) and shaken for 30 min. Lysed cells were centrifuged at 100,000g for 30 min, and the pellets were homogenized in ice-cold TS buffer (10 mM Tris-HEPES and 250 mM sucrose, pH 7.4) supplemented with the same mixture of protease inhibitors using a tight fitting Dounce homogenizer for 25 strokes. After centrifugation at 4000g for 20 min, the supernatant was centrifuged at 100,000g for 60 min. The resulting pellet was resuspended in TS buffer and passed through a 27-gauge needle 25 times. Protein concentration was determined with a Bio-Rad protein assay kit. Crude membrane vesicles were dispensed in aliquots, frozen in liquid nitrogen, and stored at −80 °C until use. Characterization of BSEP Membrane Vesicles. Uptake studies of 3H-TCA into BSEP membrane vesicles were performed using a rapid filtration technique. TS buffer supplemented with a substrate mixture of 0.1 μM 3H-TCA and 5 mM ATP (or 5 mM AMP in the case of controls), and 12 mM MgCl2 was added to the membrane vesicles in a final volume of 27.5 μL. The reaction was started by incubating the mixture at 37 °C. After 5 min, the reaction was stopped by placing samples on ice. Immediately, 150 μL of ice cold TS buffer was added to each reaction well. Diluted samples were filtered by a Multiscreen HTS-Vacuum Manifold filtration device (Millipore, EttenLeur, The Netherlands) through 0.45-μm-pore 96-wells Millipore filters that were preincubated with TS. After adding 2 mL of OptiFluor scintillation fluid (PerkinElmer, Waltham, MA, USA) to each filter and subsequent liquid scintillation counting, uptake of 3H-TCA into membrane vesicles was determined by measuring radioactivity associated with the filters. In control experiments, ATP was substituted with AMP, and the EYFP vesicles were used in the presence of ATP and AMP. Net ATP-dependent transport was calculated by subtracting average radioactivity values measured in the control experiments from those measured in BSEP vesicles in the presence of ATP. Each experiment was performed in triplicate. Inhibition Experiments. To evaluate the inhibitory effects of initial compounds, validation compounds, and selected compounds, the 3H-TCA transport assay described above was performed in the absence or presence of 50 μM (see Discussion and Conclusions for more details on why this concentration was chosen) of each of the compounds for a period of 5 min. As common for these types of studies, inhibition experiments were conducted at a substrate concentration ≪Km so that observed inhibitory effects are reflective of intrinsic inhibitory potency of the tested compounds. As the Km of 3 H-TCA for BSEP is 4.1 ± 0.5 μM,15 we used a 3H-TCA substrate concentration of 0.1 μM. The known BSEP inhibitor rifampicin (50 μM) was used as the control. When DMSO was used to dissolve the test compounds, the end concentration in each well did not exceed 1%, and appropriate vehicle controls were included in the experiments. Generation of Multiple Conformations. All in silico modeling steps and computational analyses were performed with the molecular software package MOE, version 2013.08.14 Three-dimensional molecular structures were imported from PubChem16 if available; otherwise, structures were drawn manually in MOE. The preparation of the molecules prior to performing modeling tasks was done following the MOE wash protocol. This adjusts the protonation state of the compound to that present at physiological pH (7.4). In the next step, the database energy minimization protocol with force field MMFF94x was used to enforce low energy conformations of the molecules.14 The MOE pharmacophore screening method uses a rigid placing algorithm to fit the compounds in the pharmacophore model. Since molecules are not rigid, multiple conformations need to be generated for each compound prior to the pharmacophore screening. For the initial compounds and the validation compounds, multiple conformations were generated with a MOE conformational search (a low mode search with an energy window = 20 kcal/mol and conformation limit = 250, and RMSD limit = 0.25 Å). The low mode algorithm for conformation generation uses a short (∼1 ps) run of molecular dynamics, followed by energy minimization. Only the conformations with the lowest overall energy are kept, having a

due to the fact that the crystal structure of human BSEP is yet to be discovered. To circumvent this, we now explore a ligandbased approach as the starting point for the development of a pharmacophore model. This is an in silico model describing functional groups common to molecules associated with BSEP inhibition and also fixing the intramolecular distances between these important molecular features. We present a combination of in silico studies with an in vitro testing approach to develop and validate a pharmacophore model of human BSEP inhibition.9−12

■

MATERIALS AND METHODS

Compound Sets. SMILES or two-dimensional-representations along with molecular properties for all compounds are given in the Supporting Information. We used a number of compound sets throughout the article. The initial compounds are drug-like molecules, which have been investigated in prior work for their inhibitory effects against the sodium-taurocholate cotransporting polypeptide (NTCP, SLC10A1), a hepatocyte influx transporter for bile acids (Table S1, Supporting Information).11 As NTCP and BSEP display substrate overlap with respect to bile acid transport, we hypothesized that this data set would be enriched with respect to BSEP inhibitors as well, thus reducing the number of compounds that will have to be screened to identify an adequate number of hits for building a pharmacophore model. The validation compounds are 60 drug-like compounds, many of which are on-the-market drugs (Table S2, Supporting Information). The selected compounds are six compounds (Table S3, Supporting Information), selected from the results of a virtual screening of the CoCoCo (Commercial Compound Collection) database maintained by the BioChemoInformatics Laboratory, University of Bologna, Italy. The CoCoCo database is a publicly available compound database specifically designed for in silico drug design projects.13 For the initial compounds, the validation compounds, and the CoCoCo database, the octanol/water partition coefficient (SlogP), molecular weight, charge, number of rotatable bonds, number of hydrogen bond acceptors, and number of hydrogen bond donors were calculated with the molecular software package MOE, version 2013.08.14 Chemicals. DMSO and all other reagents were of the highest obtainable grade available and were purchased from local suppliers. 3 H-taurocholate (3H-TCA) (specific activity 5.0 Ci/mmol) was obtained from PerkinElmer (Waltham, MA, USA). The initial compounds were obtained from TimTec (Newark, DE, USA), Sigma (Zwijndrecht, The Netherlands), Life Chemicals (Niagara-onthe-Lake, Ontario, Canada), Chembridge (San Diego, CA, USA), TCI (Zwijndrecht, The Netherlands) and Enamine (Kiev, Ukraine). The validation compounds were obtained from Sigma, Sequoia Research Products (Pangbourne, United Kingdom), and local suppliers. From the selected compounds, A and B were obtained from ChemDiv (San Diego, CA, USA) and C, D, E, and F from Enamine (Kiev, Ukraine) (Table S3, Supporting Information). All compounds were dissolved in ultrapure water or DMSO, depending on solubility in either vehicle. In all experiments, the appropriate vehicle controls were included. Cell Culturing and Membrane Vesicle Preparation. Overexpression of human BSEP was achieved in Human Embryonic Kidney 293 (HEK293) cells. BSEP was cloned into the Gateway entry vector, and baculovirus was produced as described in the Bac-to-Bac manual (Life Technologies, Bleiswijk, The Netherlands). As a control, the enhanced yellow fluorescent protein (EYFP) was also cloned into the gateway entry vector (Invitrogen, Bleiswijk, The Netherlands) and expressed. HEK293 cells were grown until 40% confluency in 500 cm2 triple flasks (Sanbio, Uden, The Netherlands) in DMEM supplemented with 10% fetal calf serum at 37 °C under 5% CO2-humidified air. The culture medium was removed, and 10 mL of virus and 25 mL of medium were added. The cells were incubated for 20 min at 25 °C, after which 40 mL of medium containing sodium butyrate (5 mM) to increase transporter expression was added. HEK-BSEP cells were harvested 3 days after transduction by centrifugation at 3000g for 10 min. All steps were performed at 4 °C. The pellets were resuspended 874

dx.doi.org/10.1021/tx5000393 | Chem. Res. Toxicol. 2014, 27, 873−881

Chemical Research in Toxicology

Article

maximal energy difference of 20 kcal/mol (energy window) and an RMSD greater than 0.25 Å (RMSD limit) to the conformation with the lowest energy.17 For the virtual screening of the CoCoCo database, the CoCoCo-MC database as provided by the authors was transformed into a MOE database.13,14,18 Pharmacophore Modeling. Experimentally identified inhibitors were used to deduce structural information for this transporter. A 3D alignment of the compounds 02, 12, 15, 26, and 32 was made using the flexible alignment (default settings) option in MOE to superpose chemically similar groups of the compounds. The final 3D alignment was used to generate pharmacophore features by calculating a pharmacophore consensus model (with tolerance = 1.3 and threshold = 75%) in MOE. This will mark the most prominent chemical properties within these compounds with a pharmacophore feature. The resulting pharmacophore model is an abstract representation of the important chemical properties needed for BSEP inhibition. The model was further optimized by placing exclusion volumes around the pharmacophore model to resemble the shape of the transporter binding site. The shape of the transporter can be estimated with respect to the difference in fitting of the inhibitors and noninhibitors into the pharmacophore model. Therefore, a ligand shape (2 Å) was placed around the top 10 conformations of the inhibitors detected by the pharmacophore model searching the initial compounds. Comparing this shape with the top 10 conformations of the noninhibitors detected by the pharmacophore model allows selecting all atoms of the noninhibitors that were outside the shape. On this selection of noninhibitor atoms, exclusion volumes with a radius of 1 Å were placed. The exclusion volumes were added to the pharmacophore model (Figure S2, Supporting Information).10,14 Validation Pharmacophore Model. The pharmacophore model was validated by performing a virtual pharmacophore screening (pharmacophore search in MOE) using the experimental results obtained for a different compound set (the validation compounds). The pharmacophore search fits a rigid conformation of a compound in the pharmacophore model and attempts to match the chemical groups of the compound to the corresponding pharmacophore features. When at least one of the conformations of a compound matches the pharmacophore, the compound is considered to be a hit (for more details see Results).14 Virtual Screening. A virtual screening of the CoCoCo database was performed using the pharmacophore model. The resulting hits were first filtered, using an adjusted Lipinski filtration (accepting all compounds with −2 < charge < +1; −1.35 < SlogP < 6.0; molecular weight ≤1000 Da, number of rotational bonds ≤10; number of hydrogen bond acceptors ≤10; number of hydrogen bond donors ≤5) to correct for the large molecular weight of the natural substrates of BSEP.19 For every compound, only the best conformation was selected, based on the root-mean-square deviation (RMSD) value, which quantifies the overlap between the pharmacophore features and the matching molecular features of the compounds. Driven by vendor availability, compounds with different scaffolds and charge were subsequently selected from the hit list for experimental validation in the lab. Evaluation of Pharmacophore Screening. TP = number of true positives, FP = number of false positives, TN = number of true negatives, and FN = number of false negatives. Precision = TP/(TP + FP), recall = TP/(TP + FN), and specificity = TN/(TN + FP). Matthews correlation coefficient (MCC) is calculated with eq 1. The MCC ranges from −1 (no correlation) to 1 (full correlation).

MCC =

performed to analyze all tested compounds. For the inhibition experiments, a comparison between every single value and the control was performed, which was analyzed with a Dunnett’s posthoc test. When all groups were compared, Bonferroni’s multicomparison posthoc test was used. Differences were considered statistically significant when P < 0.05. To check for correlation between the effect on BSEP-mediated 3H-TCA uptake and molecular properties of the compounds, several plots were made using MOE and GraphPad (Figure S1, Supporting Information). Regression and subsequent analyses for these plots were also performed in GraphPad.

■

RESULTS Characterization of BSEP Membrane Vesicles. The membrane vesicles were characterized by comparing the amount of 3H-TCA taken up in BSEP vesicles in the presence and absence of ATP, as well as by comparing uptake in BSEP vesicles to uptake by EYFP vesicles. As can be seen in Figure 1,

Figure 1. Characterization of the BSEP-overexpressing membrane vesicles. Control (EYFP) and BSEP vesicles were exposed to 3H-TCA in the presence or absence of ATP for 5 min. Rifampicin (Rif, 50 μM), a known inhibitor of BSEP, was included as a positive control. The data represent the mean ± SEM of observations from a representative experiment performed 3- to 9-fold (*** indicates P < 0.001).

the vesicles displayed ATP-dependent transport of 3H-TCA. Rifampicin (50 μM), a known BSEP inhibitor, showed a significant reduction in BSEP-mediated 3H-TCA uptake. We did not detect any significant differences between control uptake with AMP in BSEP vesicles and ATP-dependent uptake in EYFP vesicles. Therefore, in the following experiments only EYFP vesicles incubated in the presence of ATP were included as controls. Initial Compounds. An inhibition experiment was performed to identify inhibitors of BSEP-mediated 3H-TCA uptake in the initial compounds set. Compounds 02, 12, 14, 15, 26, and 32 were significant inhibitors of 3H-TCA uptake (P < 0.05, Figure 2). Compounds 02, 12, 15, 26, and 32 were considered to be strong inhibitors of BSEP transport, showing more than 50% inhibition. These structures were used for building the pharmacophore model. Validation Compounds. From the validation data set (Figure 3), baicalein, benzbromarone, glibenclamide, indocyanine green, losartan (potassium salt), mifepristone, pioglitazone hydrochloride, saquinavir mesylate, simvastatin (sodium salt), tacrolimus (FK506), taurocholate, and troglitazone could be identified as strong inhibitors, inhibiting BSEP-mediated 3H-

TP × TN − FP × FN (TP + FP)(TP + FN )(TN + FP)(TN + FN ) (1)

Data Analysis and Statistical Evaluation. Data analysis was performed using the GraphPad Prism (version 5.00, GraphPad Software, San Diego, CA, USA) software package and the R statistical computing software package (version 2.15.2, R Foundation for Statistical Computing, Vienna, Austria). Data are expressed as the mean ± SEM. One-way ANOVA with appropriate posthoc testing was 875

dx.doi.org/10.1021/tx5000393 | Chem. Res. Toxicol. 2014, 27, 873−881

Chemical Research in Toxicology

Article

Information). The plots show a weak correlation for SlogP (r2 = 0.1628, P < 0.0001) and molecular weight (r2 = 0.1554, P < 0.0001), and no correlation for the other properties. For further investigation of this observation, a combined correlation plot was made for SlogP and molecular weight (Figure 4). This combined plot shows that large and hydrophobic compounds are indeed more likely to be BSEP inhibitors. Using our experimental data, we were able to set threshold values for BSEP inhibition. Compounds need to have at least an SlogP of −1.35 (baicalein) and a molecular weight of 355 Da (pioglitazone hydrochloride) to be a BSEP inhibitor. Below these values, no BSEP inhibitors were found. These threshold values define an area (Figure 4) where BSEP inhibitors and also noninhibitors are detected, and on the basis of molecular properties, it is not possible to distinguish between BSEP inhibitors and noninhibitors. Pharmacophore Modeling. To retrieve further insights into other descriptors important for BSEP inhibition, a pharmacophore model containing the molecular features necessary for BSEP inhibition was generated (Figure 5a). A total of two hydrogen bond acceptor (HBA)/anionic features (F2 and F3) in combination with two HBA projection vector features (F4 and F8) and four hydrophobic/aromatic features (F1, F5, F6, and F7) are part of the resulting pharmacophore model (Figure 5b and Table 1). The HBA projection vector features indicate the hydrogen bond direction toward the receptor. For the pharmacophore search, only 7 out of the 8 features have to be matched by each compound. In addition, the large hydrophobic core (F1, F5, F6, and F7) was marked essential. To optimize the model, exclusion volumes were placed to take the shape of the transporter binding site into account (Figure S2, Supporting Information). The final pharmacophore model recognized 5 out of 5 inhibitors (100%) and 11 out of 27 noninhibitors (41%) from the initial compounds (Table 2). In the area of high SlogP and molecular

Figure 2. Effect of the initial compounds on BSEP-mediated 3H-TCA uptake (0.1 μM, 5 min). Rifampicin was included as a control inhibitor. Compounds were tested at a concentration of 50 μM. The data are expressed as a percentage of control and represent the mean ± SEM of observations from 4 experiments each performed in duplicate (* indicates P < 0.05, ** indicates P < 0.01, and *** indicates P < 0.001). The filled bars represent the compounds that inhibit BSEP activity more than 50%.

TCA transport by more than 50%. Pioglitazone hydrochloride was the most potent inhibitor at a concentration of 50 μM, reducing 3H-TCA uptake to 1.6 ± 0.4% of the control. Also, 50 μM vinblastine, tolbutamide, rosuvastatin, quinine, mitoxantrone, glimepiride, fluvastatin, dipyridamole, and cholate demonstrated statistically significant inhibition (P < 0.05). The remaining BSEP activity for these compounds ranged between 57 ± 4% (dipyridamole) and 83 ± 2% (glipizide) of the control. Analysis of Compound Properties. For all tested compounds, the BSEP-mediated 3H-TCA uptake percentages were plotted against SlogP, molecular weight, charge, number of rotatable bonds, number of hydrogen bond acceptors, and number of hydrogen bond donors (Figure S1, Supporting

Figure 3. Effect of the validation compounds on BSEP-mediated 3H-TCA uptake (0.1 μM, 5 min). The data are expressed as a percentage of control and represent the mean ± SEM of observations from a representative experiment performed in triplicate (* indicates P < 0.05, ** indicates P < 0.01, and *** indicates P < 0.001). The filled bars represent the compounds that inhibit BSEP activity for at least 50%. 876

dx.doi.org/10.1021/tx5000393 | Chem. Res. Toxicol. 2014, 27, 873−881

Chemical Research in Toxicology

Article

Figure 4. Plot showing SlogP and molecular weight for initial compounds and validation compounds. The open circles represent the noninhibitors (remaining BSEP activity >50% of the control value), and the filled circles represent the BSEP inhibitors (BSEP-mediated 3H-TCA uptake ≤50% of the control). The crossed circles represent the compounds recognized by the pharmacophore model. The lines demarcate a gray area (upper right quadrant) in which compounds have a molecular weight and lipophilicity that appears to be associated with an increased likelihood for BSEP inhibition.

Figure 5. (a) Alignment of the active initial compounds 02, 12, 15, 26, and 32 (gray) and the consensus pharmacophore. (b) Compound 26 fitted to the BSEP pharmacophore. The pharmacophore consists of 8 features, including 4 essential hydrophobic/aromatic features (F1, F5, F6, and F7 shown in green) and 2 HBA/anionic features (F2 and F3 shown in pink) that are associated with HBA projection vector features (F4 and F8 shown in cyan). The model also includes exclusion volumes (see Figure S2, Supporting Information).

weight (Figure 4), the pharmacophore model predicts 15 out of 26 noninhibitors correctly, which otherwise would have been classified as inhibitors based on their molecular properties. Validating the Pharmacophore Model. The pharmacophore model was further evaluated by performing a pharmacophore search on the set of validation compounds, recognizing 9 out of 12 inhibitors (75%) and 8 out of 47 noninhibitors (17%) (Table 2). Of the 47 noninhibitors, 33 could be classified as noninhibitors using the cutoff values as defined in Analysis of Compound Properties. For the remaining

14 noninhibitors, the pharmacophore model is able to predict 9 noninhibitors correctly (Table 2). Virtual Screening. The pharmacophore model was then used to screen the CoCoCo database containing 6,957,134 unique commercially available compounds, and 608,497 compounds (8.7%) were found to fit the model. An adjusted Lipinski filtration was used to identify all biocompatible compounds.19 The adjustment of the standard filter consisted of the maximal weight of the Lipinski filter being increased to 1000 Da because most known BSEP inhibitors have a high molecular weight. This resulted in 491,465 (7.1%) hit 877

dx.doi.org/10.1021/tx5000393 | Chem. Res. Toxicol. 2014, 27, 873−881

Chemical Research in Toxicology

Article

Table 1. Pharmacophore Featuresa name

typeb

radius [Å]

essential

F1 F2 F3 F4 F5 F6 F7 F8

Aro|Hyd O2|Acc|Ani O2|Acc|Ani Acc2 Aro|Hyd Aro|Hyd Aro|Hyd Acc2

1.26 1.00 1.60 1.10 1.28 1.60 1.60 1.80

yes no no no yes yes yes no

a

For each feature, the type, according to the nomenclature in MOE, the radius, and whether the feature is essential or not are given. bAro = aromatic, Hyd = hydrophobic, O2 = carboxylate centroid, Acc = hydrogen bond acceptor, Ani = anionic atom, and Acc2 = hydrogen bond acceptor projection vector.

compounds, ranked according to their RMSD values. An RMSD value of 0.0 Å indicates a perfect fit to the pharmacophore model; the further the deviation from the centers of the pharmacophore features, the higher is the RMSD value. To select structurally diverse compounds, a visual inspection of the top hits was performed. Restricted by vendor availability, we selected six structurally diverse compounds (selected compounds) differing in charge (formal charge of −2, −1, 0, and +1), molecular weight, and differing in their scaffolds compared to those of initial and validation compounds for experimental validation in the laboratory (Figure 6a and Table S3, Supporting Information). Experimental Testing of the Selected Compounds. Inhibition experiments with the selected compounds showed that all six compounds gave a statistically significant inhibition of BSEP-mediated 3H-TCA uptake (P < 0.001, Figure 6b). Compounds A, B, C, and F reduced uptake levels by more than 50% and were considered strong BSEP inhibitors.

Figure 6. (a) Molecular structures of the compounds selected by virtual screening. (b) Effect of these compounds on BSEP-mediated 3 H-TCA uptake (0.1 μM, 5 min). The data are expressed as a percentage of the control and represent the mean ± SEM of observations from 3 experiments, each performed in triplicate (*** indicates P < 0.001). The filled bars represent the compounds that inhibit BSEP activity more than 50%.

■

DISCUSSION AND CONCLUSIONS A ligand-based pharmacophore model is developed and presented in this article, describing the interactions necessary for human BSEP inhibition. To our knowledge, such a model has not been published before. The building and evaluation of the model was based on extensive uptake experiments in BSEP overexpressing membrane vesicles. The majority of vesicular BSEP interaction studies reported in literature is based on membrane vesicles derived from insect cells.7,20,21 We used a human cell-line (HEK293) to overexpress BSEP and isolate membrane vesicles. It has been shown that insect cells have a lower cholesterol content, which affects the transport activity of the BSEP transporter.22,23 Our test system demonstrated appropriate BSEP-mediated ATP-dependent 3H-TCA transport, and rifampicin, the known BSEP inhibitor, was also capable of inhibiting BSEP transport.24

For the generation of the pharmacophore model, we used the five BSEP inhibitors identified from the initial compounds. The low molecular diversity of the initial compounds led to 8 highly specific pharmacophore features (Figure 5b). To be able to detect a wider range of molecules, while still selecting for a large hydrophobic core, we decided to make the 4 hydrophobic/ aromatic features essential. This means that all of the essential features have to be matched, while only 3 of the 4 remaining features have to be matched. We also discovered that many of the inactive compounds from the initial compounds set were

Table 2. Pharmacophore Screening Results

a

compound set

TPa %

FPb %

precision

recall

specificity

MCCc

initial initial (MW>355 Da, SlogP > −1.35) validation validation (MW>355 Da, SlogP> −1.35)

100 100 75 75

41 42 17 36

0.31 0.31 0.53 0.64

1.00 1.00 0.75 0.75

0.59 0.58 0.83 0.64

0.43 0.42 0.52 0.39

True positives. bFalse positives. cMatthews correlation coefficient. 878

dx.doi.org/10.1021/tx5000393 | Chem. Res. Toxicol. 2014, 27, 873−881

Chemical Research in Toxicology

Article

have reported the BSEP activity to be highly dependent on the molecular weight and the logP of the compounds. Although there is a similar trend visible in our data set, it is not possible to distinguish between inhibitors and noninhibiters based solely on molecular weight and SlogP. On the basis of the correlation shown in Figure 4, we can conclude that the inhibitors should have a molecular weight of at least 355 Da and a SlogP of at least −1.35. The 57 compounds in the marked area of Figure 4 would thus all be classified as BSEP inhibitors, but this prediction is only 30% correct. To correctly classify compounds as BSEP inhibitors, a more specific model is needed, which includes more features than SlogP and molecular weight alone. The BSEP pharmacophore model correctly predicted 38 out of the 57 compounds (67%) in the marked area in Figure 4. This model improves the classification of BSEP inhibitors by 37%. Hirano et al. published a 2D-QSAR model of BSEP in 2006.8 They suggest that a (thio)ester group attached to the carbon of a heterocyclic ring and that carbocyclic systems with at least one aromatic ring are important for BSEP interaction. We found that this is consistent with the results of our model. Both troglitazone and pioglitazone have such an ester incorporated in their thiazolidinedione groups and are potent BSEP inhibitors. We also identified several carbocyclic systems with an aromatic ring among the active BSEP inhibitors. Since the location of the aromatic ring does not seem to be restricted, we incorporated 4 aromatic/hydrophobic pharmacophore features into our model, instead of just 4 hydrophobic features. Ouabain and mycophenolic acid also have the described ester but are not BSEP inhibitors. Our model managed to exclude ouabain from the hits because it was too large to fit inside the area defined by the exclusion volumes. Mycophenolic acid was recognized by the pharmacophore model, as it had a borderline fit to each of the features of the model. It can, however, be excluded by the low molecular weight of 319.33 Da.8,31,32 Traditionally, the type of screening assay presented in this article uses micromolar concentrations of the test compounds. Whether in vivo similar intracellular concentrations will be reached is unknown, particularly in the early stages of drug development, when human dose predictions may not have been performed. Horikawa et al.33 and Dawson et al.20 published sets of maximum inhibitor concentrations at the inlet of the liver. Most of the inhibitor concentrations are