Interspecies comparison of putative ligand binding sites of human, rat and mouse P-glycoprotein.

Prior to the clinical phases of testing, safety, efficacy and pharmacokinetic profiles of lead compounds are evaluated in animal studies. These tests are primarily performed in rodents, such as mouse and rats. In order to reduce the number of animal experiments, computational models that predict the outcome of these studies and thus aid in prioritization of preclinical candidates are heavily needed. However, although computational models for human off-target interactions with decent quality are available, they cannot easily be transferred to rodents due to lack of respective data. In this study, we assess the transferability of human P-glycoprotein activity data for development of in silico models to predict in vivo effects in rats and mouse using a structure-based approach. P-glycoprotein (P-gp) is an ATP-dependent efflux transporter that transports xenobiotic compounds such as toxins and drugs out of cells and has a broad substrate and inhibitor specificity. Being mostly expressed at barriers, it influences the bioavailability of drugs and thus contributes also to toxicity. Comparison of the binding site interaction profiles of human, rat and mouse P-gp derived from docking studies with a set of common inhibitors suggests that the inhibitors share potentially similar binding modes. These findings encourage the use of in vitro human P-gp data for predicting in vivo effects in rodents and thus contributes to the 3Rs (Replace, Reduce and Refine) of animal experiments.

Introduction

The efflux transporter P-glycoprotein (P-gp) is a protein of high interest among other major anti-targets (Cramer et al., 2007). It is expressed in tissues such as intestine, liver, kidney, placenta, testis, and in the capillary endothelial cells of the brain (Seelig, 1998; Thiebaut et al., 1987), and plays an important role in the absorption, distribution and excretion of many drugs. Overexpression of P-gp has been implicated in resistance to multiple chemotherapeutic drugs and is a widely accepted mechanism underlying multidrug resistance (Aller et al., 2009; Fojo et al., 1987; Widmer et al., 2003). Co-administration of a P-gp inhibitor with a drug can lead to altered disposition of the latter, resulting in elevated plasma levels of the drug which could lead to adverse effects (Bussey, 1982; Tsuji, 2002; Verschraagen et al., 1999). Furthermore, the partial blockade of P-gp expressed in the blood–brain barrier or placenta could lead to an increased distribution of a co-administered drug in the corresponding organs. Thus, concomitant administration of substrates and P-gp inhibitors may lead to adverse drug reactions and organ toxicities (Balayssac et al., 2005). In this respect, the United States Food and Drug Administration (US FDA) guidance requires new drug candidates to be routinely screened against P-gp as part of the clinical drug interaction studies (“Clinical Drug Interaction Studies — Study Design, Data Analysis, and Clinical Implications Guidance for Industry,”, 2017; Klepsch et al., 2011). Therefore, computational methods that characterize P-gp interactions and thus guide the prioritization of compounds in the early phase of the drug discovery process are of considerable interest (Schneider, 2010).

In early stages of drug development, pharmacokinetic and toxicity profiles of a candidate drug are evaluated in animal models (typically rats or mouse) prior to the clinical phases of testing in humans. A substantial amount of experimental data against human P-gp is already available and has been utilized for the development of in silico models (see e.g. livertox.univie.ac.at). However, besides developing in silico models for the prediction of ligands for human P-gp, it would be beneficial to also establish models for rat and mouse P-gp in order to predict the outcomes of preclinical animal studies. Unfortunately, limited availability of experimental data for rat and mouse P-gp restricts the development of such models. In this context the question arises, whether predicted interaction profiles of ligands with human P-gp could be transferred to rodent P-gp. This would require a comprehensive comparison of the putative binding sites of the P-gp structures across species. Literature sheds little light on this, suggesting the need for exploration of species-related differences in P-gp mediated drug transport activity (Martignoni et al., 2006; Schwab et al., 2003; Suzuyama et al., 2007).

Inhibition of P-gp activity as a result of drug interactions has been reported in both animals and humans (Bussey, 1982; Choo et al., 2000; Pedersen, 1985), but only a few studies discussed species-related differences in the inhibitory effects on the P-gp function (Chu et al., 2013; Suzuyama et al., 2007; Zolnerciks et al., 2011). A few studies proposed moderate species differences, human vs. rat (Molden et al., 2000), human vs. mouse (Adachi et al., 2001; Lin and Yamazaki, 2003) and also among the three species (human vs. rat vs. mouse) (Katoh et al., 2006), while a few other studies reported no significant differences between human, rat and mouse P-gp (Chu et al., 2013; Feng et al., 2008; Hsiao and Unadkat, 2012). Nicklisch et al. (2016) compared the binding site residues of mouse and Thunnus albacares (Yellowfin tuna) P-gp. Despite the low similarity, of the 15 residues of mouse P-gp that interacted with the inhibitor BDE-100, 13 were conserved in the tuna P-gp sequence, and the remaining two did not play a major role in the binding site. However, it must be noted that only a small number of compounds were tested in these studies. It might thus well be that the inhibitory effects on P-gp-mediated drug transport are subjective to both the chemical structure of substrates/inhibitors and to the species. Moreover, it is not yet clear if the possible species differences in the inhibitory effects of P-gp activity are due to differences in binding site residues of P-gp, which is therefore worth investigating.

To the best of our knowledge, no computational study compared the binding site interaction profiles of P-gp across different species (human, rat and mouse) so far. In this study, we used a structure-based approach to compare their binding sites in order to derive information concerning potential species differences in P-gp-mediated drug transport. Since an X-ray crystal structure is available for mouse P-gp alone, homology modeling was performed to construct the models for human P-gp and for rat P-gp. Subsequently, docking of common inhibitors of rat, mouse and human P-gp was performed. Next, known inhibitors of human P-gp were docked into the models of the three species followed by an analysis of the interactions between the inhibitors and binding site residues. The interaction profiles of the P-gp binding sites of the three species were then compared to evaluate the transferability of in vitro human P-gp data for development of models to predict effects in rat and mouse.

Methodology

Dataset

A substantial amount of human P-gp data is made publicly available through previous literature reports (Broccatelli et al., 2011; Chen et al., 2011; Klepsch et al., 2014). However, due to the limited availability of rat P-gp data in public domain bioactivity databases such as ChEMBL (Gaulton et al., 2012; Willighagen et al., 2013) and BindingDB (Liu et al., 2007), an exhaustive literature search was performed. A total of 18 rat P-gp inhibitors could be identified that are known to also inhibit both human P-gp and mouse P-gp. Due to the inconsistencies in the assay conditions, these compounds unfortunately could not be utilized to compare inhibitory profiles across the species. Suzuyama et al. (2007) studied the species differences (human, monkey, canine, rat and mouse) in the inhibitory effects of the prototype P-gp inhibitors quinidine and verapamil. These two drugs served as the starting point for in silico comparison of binding site interaction profiles across the species. Further, we also extracted the human P-gp data from Broccatelli et al. (2011) in order to perform protein-ligand interaction fingerprint (PLIF) analysis and to identify the common functional group residue interactions among the three species. The dataset was standardized according to the procedure described in Pinto et al., 2012. (Pinto et al., 2012) The final dataset contained a total of 1161 compounds (612 inhibitors and 549 non-inhibitors).

Homology modeling

Based on sequence identity and resolution, the corrected mouse P-gp structure (mdr1a; PDB ID: 4M1M; UNIPROT ID: P21447) was selected as the most structurally related template protein for constructing the homology models for human P-gp (MDR1; UNIPROT ID: P08183), rat P-gp (MDR1a; UNIPROT ID: Q9JK64 and MDR1b; UNIPROT ID: P43245) and mouse P-gp (mdr1b; UNIPROT ID: P06795). Rat and mouse P-gp proteins are encoded by two paralogous genes namely MDR1a and MDR1b that show a sequence identity of 83% (Chu et al., 2013; Devault and Gros, 1990). Therefore, we constructed in total four homology models to consider the paralogs too. Homology models were constructed using MODELLER 9.13 (Eswar et al., 2007) and the Prime module in Maestro (Schrödinger, Inc. V-10.1.013) (Jacobson et al., 2004, 2002). The energy minimized models were further evaluated using DOPE score (Shen and ) and GA341 score (John and Sali, 2003; Melo et al., 2002). Quality of the stereochemical parameters and the normality of the structures were checked using the PROCHECK program, included in the PDBsum analysis (Laskowski et al., 1993). Ramachandran plot (Zhou et al., 2011) and G-factor (Engh and Huber, 1991), and finally the Q-score (Benkert et al., 2008, 2009) values were evaluated to identify the best homology models. The electrostatic potential surface (EPS) of each of the three best models for the three species was also calculated and compared using MOE 2013 ().

Sequence alignment

Sequence alignment was performed using ClustalX (Larkin et al., 2007) and verified by including secondary structure predictions. Subsequently, the alignment was analyzed using Jalview (Supplementary Figs. S1–S4) (Clamp et al., 2004; Waterhouse et al., 2009).

Binding site identification and molecular docking

In order to avoid any bias, the binding site for all five structures (human MDR1, rat MDR1a, rat MDR1b, mouse mdr1a, and mouse mdr1b) was defined as the complete transmembrane region, taking 20 Å around the coordinate of the center point to allow subsequent flexible docking of a series of P-gp inhibitors. The protein was prepared using the Protein Preparation Wizard of the Schrödinger Suite (2015) (Sastry et al., 2013; ). Hydrogen atoms were added, and optimal protonation states and ASN/GLN/HIS flips were determined. To assess their correct protonation states, ligands were prepared using the LigPrep module of the Schrödinger Suite, (, ) which produces low-energy 3D structures that can be used for docking. The OPLS_2005 force field was used for minimization of the structures. Different ionization states were generated by adding or removing protons from the ligand at a target pH of 7.0 ± 2.0 using Epik version 3.1., (Greenwood et al., 2010; Shelley et al., 2007) and tautomers were generated for each ligand. To generate stereoisomers, the information on chirality from the input file for each ligand was retained as is for the entire calculation. All docking runs were performed in high-throughput mode with GlideXP (Friesner et al., 2006; Halgren et al., 2004) docking in Maestro. We also used the genetic algorithm-based GOLD suit (version 5.2.0) (Jones et al., 1997; Verdonk et al., 2003) for docking.

Protein ligand interaction fingerprint (PLIF)

A PLIF summarizes the interactions between a ligand and a protein using a molecular fingerprint scheme. We generated two types of PLIFs that differ in the information encoded. The first PLIF encodes residues involved in an interaction with the ligand at each bit position. The second type encodes the functional group of the ligand that interacts with the residue. For this, the substructure patterns of 100 functional groups (in SMARTS notation) were extracted from the Daylight website (http://www.daylight.com/dayhtml_tutorials/languages/smarts/smarts_examples.html#GROUP). All PLIF bits were calculated with the MOE 2013 () built-in function CalculateRawInteractions using a 1% threshold for molecular interactions and a 20% threshold for surface contacts. The function was embedded in an in-house SVL script and was post-processed to enable calculation of functional group PLIFs.

Results and discussion

Predicting interactions of small molecules with membrane protein structures has always been challenging. Nevertheless, visualization of the 3D models contributes to the comprehension of the physical and chemical properties of these biomolecules, and of their intermolecular interactions with endogenous and exogenous compounds. Due to the lack of crystal structures for human and rat P-gp, homology modeling and computational ligand docking were used to generate structure-based hypotheses for protein-ligand-interactions.

ABC transporters are transmembrane proteins that are in general difficult to be resolved via crystallization (Klepsch et al., 2010). In such cases, homology modeling is the method of choice for structure-based studies. The homology models generated in this study resemble the open-inward (or apo/ground) state of P-gp. This state was considered because it resembles the first step of the basic catalytic cycle for drug-binding in P-gp (Wilkens, 2015).

Since January 2014, a refined X-ray structure of a eukaryotic ABC efflux pump, ABCB1 (mouse) is available (Li et al., 2014) (PDB code: 4M1M, resolution: 3.8 Å). High sequence identities with human MDR1 (86%), rat MDR1a (94%), rat MDR1b (82%), mouse mdr1b (83%) and a moderate resolution of 3.8 Å renders 4M1M a reasonable template for homology modeling (Pajeva et al., 2009). Moreover, the secondary structure elements are also conserved among the species. When only the transmembrane domain (TMD) was analyzed, the sequence identity is greater than 85% for all structures (Supplementary Fig. S5). The best models had a normalized Dope score of less than −0.6, G-factors less than −0.12, and Qmean scores of greater than 0.60 (see Table 1). For all modelled structures, the Ramachandran plot (Supplementary Figs. S6–S9) showed excellent results with less than 1.9% of residues in generously allowed or disallowed regions. All of these residues are located in the nucleotide binding domains (NBD) or extracellular loops (ECL) and are therefore not involved in drug binding (Supplementary Figs. S10–S13). Table 1 summarizes the model assessment details for the best structure. The X-ray crystal structure and site directed mutagenesis studies on ABCB1 serve as validity tests for both helix orientation in the template (Ward et al., 2007), and the alignment used for ABC transporter modeling (Supplementary Figs. S1–S4). The homology models as well as the crystal structure displayed a V-shaped structure with analogos domain orientations.

Table 1

Results from the stereochemical validation of the homology models.

Model	Dope score	G-factor	Qmean score	Residues in generously allowed or disallowed regions (%)
Human MDR1	−0.633	−0.13	0.68	1.7
Rat MDR1a	−0.795	−0.03	0.70	1.7
Rat MDR1b	−0.703	−0.16	0.65	1.8
Mouse mdr1b	−0.808	−0.06	0.69	2.0

Sequence alignment and binding site analysis

The amino acid sequence is highly conserved among the three species (Supplementary Fig. S5), suggesting a high structural similarity (see Table 2).

Table 2

Sequence identity/similarity [%] between human, rat and mouse P-gp.

	Human-MDR1	Rat MDR1a	Rat MDR1b	Mouse mdr1a	Mouse mdr1b
Human-MDR1	100	87/93	80/90	86/92	80/90
Rat MDR1a	87/93	100	84/91	94/97	84/92
Rat MDR1b	80/90	84/91	100	82/91	93/97
Mouse mdr1a	86/92	94/97	82/91	100	83/91
Mouse mdr1b	80/90	84/92	93/97	83/91	100

Experimental techniques such as cysteine and arginine scanning and photoaffinity labeling were previously employed to determine the drug binding sites of P-gp (Loo and Clarke, 2008; Pleban et al., 2005; Seeger and van Veen, 2009; Shilling et al., 2006). Multiple binding sites were identified and binding at different sites could lead to different inhibitory effects. Well characterized binding sites are the ones of Hoechst 33,342 and Rhodamine, the so called H-sites and R-site (Loo and Clarke, 2002; Qu and Sharom, 2002). Studies also suggest the presence of an allosteric regulatory site as well as a progesterone and prazosin binding region (Martin et al., 2000; Shapiro et al., 1999). The H-site and R-site residues, (characterized by Ferreira et al. (2013)) of the three species were compared and showed a high sequence identity. This would indicate the similar arrangement of the binding sites residues and thus further pointing to the presence of a similar binding/interaction profile of the inhibitors. Mostly identical or similar residues were present in the five structures. The H-site and R-site had 77% and 65% residues identical within the three species. Those residues which show a difference, have mostly similar properties. For example, Glu180 in mouse mdr1a is replaced with Aspartic Acid in mouse mdr1b and in ratMDR1b. Both residues are charged and have acidic properties. In a few instances, charged (basic) amino acids are replaced by polar (neutral) or hydrophobic (aliphatic) amino acids in the other species but most of these residues did not participate in interactions with docked ligands. In general, the H-site has a higher percentage of charged residues (lysine, histidine, and glutamic acid residues), while the R-site has a high number of glycine, glutamine, and proline residues. Interestingly, threonine and tyrosine were not found in the H-site and R-sites, respectively. A detailed comparison of the H-site and R-site residues of the five structures is shown in Supplementary Table S1. These observations signify the harmony of electrostatic properties and molecular features of the drug recognition site (central binding cavity) in the three species. Supplementary Figs. S14–S18 represents the electrostatic potential surface (EPS) of the substrate recognition area of each of the ABC-transporter models. The EPS of the substrate recognition area in the TMDs of the human model is neutral with negative and weakly positive areas, similar to the EPS of rat and mouse models.

Molecular docking

In order to analyze the putative binding pocket of the transport protein in the three species, we proceeded with docking of a set of inhibitors. Ligand docking is a commonly used approach to identify ligand-protein interactions. However, in case of P-gp, this appears to be challenging due to various reasons. Firstly, P-gp possesses a high degree of flexibility with a large binding cavity consisting of multiple binding sites. Secondly, it can harbor more than one ligand simultaneously (Loo et al., 2003a; Lugo and Sharom, 2005). And finally, lack of a high resolution crystal structure of human P-gp necessitates the use of homology models, which add additional layers of uncertainty. A large binding pocket could also be seen in a recent structure (PDB id: 4M1M) wherein large cyclopeptides bind at different sites with partially overlapping residues (Li et al., 2014). Some of these residues are identical to those involved in rhodamine or verapamil binding (Loo et al., 2006; Loo and Clarke, 1997). Other studies reported different prazosin binding sites in hamster (Isenberg et al., 2001) and human P-gp (Ambudkar et al., 2003). Overall, it is understood that P-gp possesses a huge binding pocket with at least four distinct binding sites, with TM 6 as the helix primarily involved in binding (Klepsch et al., 2010). Therefore, we considered the complete TMD as drug binding site (DBS) and generated a large number of docking poses to prevent any bias introduced by scoring functions.

We started with docking of verapamil and quinidine into the binding pocket (complete TMD) of all models. These two compounds were chosen since IC50 values measured under the same assay conditions were available for all three species. Our study revealed that the top ranked docking poses of verapamil were found in the R-site of P-gp in all three species, which is in agreement with previous reports (Ferreira et al., 2013). The top scored docking pose for each of the five models was found in the same region of the binding pocket (R-site) are shown in Fig. 1.

Fig. 1

Best scored docking pose of verapamil: green (human MDR1), yellow (rat MDR1a), pink (rat MDR1b), red (mouse mdr1a), blue (mouse mdr1b), secondary structure of human P-gp. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We used the GlideXP scoring function from Maestro (Friesner et al., 2006; Halgren et al., 2004) to evaluate the binding poses. GlideXP docking also provides the per residue interaction energies for a particular docking pose. For each model, the residue interaction energy (RIE) for the top scored docking poses was calculated. Phe303, Tyr307, Tyr310, Phe336, Phe343, Phe728, Phe983, Met986 and Gln990 (numbering according to human-MDR1) are residues that showed more negative interaction energy values in all three species, indicating their higher contribution to binding. Residue interaction energies for all residues which are involved in interactions with verapamil can be seen in Supplementary Table S2. For example, the residues corresponding to Tyr307 in human MDR1 are Tyr306 (RIE:-3.229 kcal/mol), Tyr306 (RIE:-3.714 kcal/mol), Try303 (RIE:-6.96 kcal/mol) and Tyr299 (RIE:-5.1 kcal/mol) in rat MDR1a, rat MDR1b, mouse mdr1a and mouse mdr1b, respectively. Each of these residues contributes to the binding with more negative interaction energy. We observed less negative RIE values with residues which are different within the species (e.g. human MDR1: Met68, rat MDR1a: Leu67, rat MDR1b: Leu66, mouse mdr1a and mdr1b: Met67), suggesting their small influence on - and involvement in - interactions with verapamil. Thus, the comparison of the per residue interaction energies of the best docking pose of the three species revealed that similar binding site residues (as per the alignment) are involved in strong interactions with the ligand (see Fig. 2).

Fig. 2

Residue interaction energy for common interaction residues in human MDR1, rat MDR1a, rat MDR1b, mouse mdr1a and mouse mdr1b. x-axis denotes residue number in the order human MDR1, rat MDR1a, rat MDR1b, mouse mdr1a and mouse mdr1b, y-axis denotes the corresponding residue interaction energy (kcal/mol).

In case of quinidine, the human MDR1 residues Phe336 (RIE:-2.41 kcal/mol), Gln725 (RIE:-11.577 kcal/mol), Phe728 (RIE:-2.479 kcal/mol), Ser979 (RIE:-1.535 kcal/mol), Phe983 (RIE:-8.114 kcal/mol) and Met986 (RIE:-2.162 kcal/mol) interacted with greater negative interaction energies. The corresponding residues in rat MDR1a, rat MDR1b, mouse mdr1a and mouse mdr1b that demonstrated more negative interaction energies can be found in Supplementary Table S3. Supplementary Fig. S19 shows the RIE for common residues involved in interaction in human MDR1, rat MDR1a, rat MDR1b, mouse mdr1a and mouse mdr1b with the top scored docking pose of quinidine. Replacing a phenyl alanine in human (Phe303) and mouse (Phe299 in mdr1a and mdr1b) with another hydrophobic residue, for instance Tyr302 in MDR1a and MDR1b, still showed negative RIE values. Supplementary Fig. S20 shows the top scored docking poses of quinidine in the five models.

Thus, similar amino acids, as observed with verapamil and quinidine, also confirm the homogeneous nature of the binding site residues in different species. This would further support the hypothesis of similarity in their binding sites. Site directed mutagenesis studies on human ABCB1 also indicated that Ile306 (TMH5) (Loo et al., 2006; Loo and Clarke, 2005), Ile340 (TMH6) (Loo and Clarke, 2002), Phe343 (TMH6) (Loo et al., 2003b, 2006), Phe728 (TMH7) (Loo et al., 2006), and Val982 (TMH12) (Loo and Clarke, 2002, 2005) may participate in ligand binding. As shown in Fig. 3, these residues may form a substrate recognition site in the human ABCB1 model. The involvement of these residues in ligand binding was also confirmed by Li et al. (Aller et al., 2009; Li et al., 2014).

Fig. 3

Key residues of the substrate recognition site in the human ABCB1 model from literature (Loo et al., 2003b, 2006; Loo and Clarke, 2002, 2005).

In our previous work, we demonstrated that the Chemscore scoring function from the GOLD docking suit facilitated docking-based classification of inhibitors and non-inhibitors for P-gp (Klepsch et al., 2014) and the bile salt export pump (BSEP) (Jain et al., 2017) with reasonable accuracies. Therefore, we used the Chemscore scoring function to perform docking of all human P-gp inhibitors into human MDR1, rat MDR1a, rat MDR1b, mouse mdr1a, and mouse mdr1b structures in order to compare the interaction profiles of the binding site residues in the three species via PLIF analysis.

Maestro allows computation of different molecular interactions between binding site residues and a ligand in a specific pose. A PLIF summarizes the interactions between a ligand and a protein using a fingerprint scheme. It provides a detailed picture of the binding modes of different inhibitors. We retrospectively analyzed the PLIFs of complexes of verapamil and quinidine with structures of all three species (human MDR1, rat MDR1a, rat MDR1b, mouse mdr1a and mouse mdr1b) derived from docking, in order to compare their interaction profiles. In case of verapamil, in human MDR1, around 70% of the poses showed hydrophobic interactions with Phe336, Ile340, Phe343, Phe728 and Met986. Also, over 85% of the poses displayed interaction with Met69, Tyr310, Tyr 953, Phe983 (Fig. 4). In the rat structure, more than 73% of the residues showed interaction with Phe328, Phe335, Phe720, Met978 and over 85% residues showed interaction with Phe295, Ile298, Tyr299, Tyr302, Phe975 (Fig. 4). The percentage of binding poses in which specific residues are involved in interactions with verapamil in rat MDR1b, mouse mdr1a, and mouse mdr1b models can also be seen in Fig. 4. Supplementary Fig. S21 provides the same information for quinidine.

Fig. 4

Hydrophobic interactions common in human MDR1, rat MDR1a, rat MDR1b, mouse mdr1a and mouse mdr1bfor verapamil. X-axis denotes residue number in the order human MDR1, rat MDR1a, rat MDR1b, mouse mdr1a and mouse mdr1b, Y-axis denotes frequency of interacting residues (%).

For both verapamil and quinidine, the interacting residues that contributed to a significant number of binding poses were at similar positions in the 3D structures of the five transporters (Fig. 4, Supplementary Fig. S21). Specifically, for verapamil, Tyr310 (human MDR1), Tyr302 (rat MDR1a), Tyr310 (rat MDR1b), Tyr309 (mouse mdr1a and mdr1b) showed an interaction in more than 95% of the poses in all five docked structures. Interestingly, when an amino acid in one species was replaced with another amino acid in another species, a similar percentage of docking poses interacted with this residue. For instance, the exchange of Ile340 in human MDR1 with Leucine in rat MDR1b and mouse mdr1b showed hydrophobic interactions in almost 80% of the poses for verapamil, indicating that the interaction pattern did not change when two hydrophobic residues were interchanged. Similar PLIF-based observations could be inferenced after evaluation of the docking poses of quinidine for the three species. However, due to the lower degree of freedom (flexibility) of quinidine, relatively fewer docking poses could be obtained.

We also identified the interacting residues for a set of 612 human P-gp inhibitors that were docked into these five structures. For human P-gp, more than 70% of the inhibitors interacted with Ile306, Tyr310, Phe336, Phe343, Tyr953, Phe983 and Met986. Fig. 5 shows interacting residues common to human, rat and mouse structures. Supplementary Table S4 lists the occurrence of commonly interacting residues in the three species. PLIFs obtained from docking of a diverse set of human P-gp inhibitors into the five models revealed that similar residues were involved in the interactions, thereby further strengthening the existence of analogous binding site residues and interaction profiles in the three species.

Fig. 5

Hydrophobic interactions common in human MDR1, rat MDR1a, rat MDR1b, mouse mdr1a and mouse mdr1b for human P-gp inhibitors. X-axis denotes residue number in the order human MDR1, rat MDR1a, rat MDR1b, mouse mdr1a and mouse mdr1b, y-axis denotes frequency of interacting residues (%).

Analysis of the interactions of functional groups with protein residues

Investigation of the functional group-residue interactions for the set of 612 P-gp inhibitors docked revealed similar interaction patterns with all three species. Functional groups ether, carbonyl, alkyl carbon, nitrogen and arene showed more prominent interactions with Tyr310, Phe343, Phe983, and Met986 (numbering as per human MDR1). Corresponding residues in other species that participated in interactions are shown in Supplementary Table S5. To illustrate the outcomes, a heat map (Fig. 6, Supplementary Figs. S22–S25) plotting the interacting residues against the functional groups of ligands was generated. The color scale of the heat map represents the number of inhibitors involved in a particular interaction between a specific residue and a specific functional group. Thus, the most significant interactions could be visually identified.

Fig. 6

Heat map illustrating the PLIF analysis of the human P-gp inhibitors for human MDR1. X-axis denotes contact residues. Y-axis denotes functional groups of the ligand which are showing an interaction with the residue. Color scale signifies the number of interacting ligands.

In the three species, the H-site and R-site showed high sequence identity, suggesting the presence of similar residues at specific positions in the 3D structure. Additionally, similar functional group-residue interaction patterns were observed for human (MDR1), rat (rat MDR1a, rat MDR1b) and mouse (mouse mdr1a, and mouse mdr1b) (Supplementary Figs. S22–S25). This further strengthens the idea of utilizing human P-gp activity data, collated from in vitro studies, for structure-based modeling of rodent ligand-target interactions.

Differences in the expression, inhibition and substrate specificity of P-gp across species have already been reported (Chu et al., 2013). As stated earlier, unlike the human P-gp, rat and mouse proteins are encoded by two paralogous genes. These genes are expressed in different organs in the respective species. In humans, MDR1 mRNA is predominantly expressed in the liver, adrenal gland and kidney. In mouse, mRNA analysis suggested that Mdr1a is predominantly expressed in intestine, brain, and testis while Mdr1b is prominently expressed in adrenal gland, ovaries, placenta, and kidney. However both these genes (Mdr1a and Mdr1b) are expressed in liver and heart (Chen et al., 2003; Schinkel et al., 1994). Expression of MDR1/Mdr1 genes is less pronounced in the liver and kidney as compared to other organs in the mouse. On the other hand, enriched expression in the gastro-intestinal tract was reported in rodents, when compared to humans (Bleasby et al., 2006). Ito et al. (2011) and Uchida et al. (2011) reported that the P-gp protein expression was approximately 3-fold higher in mouse as compared to the human brain microvessels. Further, Yamazaki et al. (2001) and Baltes et al. (2007) demonstrated that the antiepileptic drugs phenytoin and levetiracetam were directionally transported by mouse, but not by human P-gp. This would indicate a difference in the P-gp substrate susceptibility between the two species (Baltes et al., 2007). For instance, the lack of P-gp-mediated transport of levetiracetam in rats but not in mouse could be explained by the species differences (Baltes et al., 2007).

A study by Schwab et al. (2003) reported comparable IC50 values in a calcein-AM assay for human MDR1, mouse mdr1a and mouse mdr1b for 28 reference compounds. Zolnerciks et al. (2011) also observed comparable IC50 values for a set of compounds against human and rat P-gp transporter and also suggested that multiple P-gp substrates would be needed to accurately predict clinically significant P-gp drug interactions, in both in vitro and in vivo (including human) drug-drug interaction studies. As mentioned earlier, Suzuyama et al. (2007) evaluated the inhibitory effects of quinidine and verapamil on P-gp-mediated drug transport using MDR1 transfected cell lines of different species. As a common observation, although the IC50 values differ between the species, it was less than 10-fold. This along with our molecular docking and PLIF analysis results signify the possibility of similar interaction profiles in the three species (human, rat and mouse), suggesting the usability and transferability of in vitro human data for development of prediction models for rat and mouse.

Conclusion

P-glycoprotein is a transmembrane efflux transporter that plays an important role in drug absorption, disposition, metabolism, and toxicity. It is essential to investigate the interactions of P-gp with candidate drugs not only to understand the contribution of P-gp to the pharmacological properties of candidate drugs, but also to evaluate their drug-drug interaction (DDI) profiles and thereby their clinical implications. In this regard, it is important to understand the binding site interaction profiles of P-gp in rodents which is poorly addressed so far due to the limited availability of experimental data. In this communication, we compared the P-gp binding sites across human, rat and mouse using molecular docking and protein-ligand interaction fingerprint analysis. The suitability of heterogeneous in vitro and in vivo data to model different aspects of human toxicity (or activity), that remained quite challenging, was recently demonstrated (Taškova et al., 2018). However, transferability of human activity data for modeling in vitro and in vivo effects in rodents presents a contrasting perspective and has been so far unexplored. To the best of our knowledge, this is the first in silico study of its kind that compares the binding sites across three different species with emphasis on their inhibitory interaction profile. It is widely acknowledged that the wealth of preclinical toxicity data aggregated by pharmaceutical industry is not completely exploited. With significant increase in the drug development costs and the pressure to reduce the use of animals for experimental testing, development of in silico systems for predicting organ and in vivo toxicity was regarded as the primary goal of the eTOX consortium (www.etoxproject.eu) (Briggs et al., 2012; Hartmann and Pognan, 2017) that ultimately aimed at early identification of potential toxic effects and risk assessment during drug development. In the light of this, exploiting human in vitro data on transporters to model in in vivo effects in animals could significantly reduce the attrition rates as well as the number of animal experiments in preclinical studies.

Our results show a significant overlap between the binding site interacting residues across the three species. This strengthens the likelihood of similar binding mode of human, rat and mouse P-gp inhibitors, thus supporting the transferability of in vitro human P-gp data for development of computational models to predict effects in rat and mouse. As shown recently, the incorporation of predicted ligand transporter interaction profiles increases the performance of selected in vivo toxicity prediction models. The transferability of human P-gp data to rodent in silico models might thus increase the predictivity of rodent in vivo toxicological outcomes. This will subsequently improve the quality of drug candidates while lowering the attrition rate during subsequent phases of drug development, and, most remarkably, reduce the number of animal experiments in preclinical studies.

Supplementary Material

Supplementary data to this article can be found online at https://doi.org/10.1016/j.ejps.2018.06.022.