Abstract
Using in vitro data, we previously built Catalyst 3-dimensional quantitative structure activity relationship (3D-QSAR) models that qualitatively rank and predict IC50 values for P-glycoprotein (P-gp) inhibitors. These models were derived and tested with data for inhibition of digoxin transport, calcein accumulation, vinblastine accumulation, and vinblastine binding. In the present study, 16 inhibitors of verapamil binding to P-gp were predicted using these models. These inhibition results were then used to generate a new pharmacophore that consisted of one hydrogen bond acceptor, one ring aromatic feature, and two hydrophobes. This model predicted the rank order of the four data sets described previously and correctly ranked the inhibitory potency of a further four verapamil metabolites identified in the literature. The degree of similarity in rank ordering prediction by these inhibitor pharmacophore models generated to date confirms a likely overlap in the sites to which the three P-gp substrates used in these studies (verapamil, vinblastine, and digoxin) bind. Alignment of the three substrate probes indicated that they are likely to bind the same or overlapping sites within P-gp. Important features on these substrates include multiple hydrophobic and hydrogen bond acceptor features, which are widely dispersed and in agreement among most of the five inhibitor pharmacophores we have described so far. These 3D-QSAR models will be useful for future prediction of likely substrates and inhibitors of P-gp.
P-glycoprotein (P-gp) is an efflux transporter highly expressed at the interface of many important organs with their environment, where it acts as a barrier limiting exposure to xenobiotics (Wandell et al., 1999a. In addition, P-gp expression in malignant cells has been associated with the multidrug resistance phenomenon resulting from the P-gp-mediated active transport of anticancer drugs from the intracellular to the extracellular compartment (Wandell et al., 1999a. Interestingly, CYP3A4, a drug-metabolizing enzyme with broad substrate specificity, seems to coexist with P-gp in organs such as the intestine and liver. These observations lead to the hypothesis that there may be an interrelationship between these two proteins in the drug disposition process. Wacher et al. (1995) described the overlapping substrate specificity and tissue distribution of CYP3A4 and P-gp. Schuetz et al. (1996) found that modulators and substrates coordinately up-regulate both proteins in human cell lines. Similarly, P-gp mediated transport was found to be important in influencing the extent of CYP3A induction in the same cell lines as well as in mice (Schuetz et al., 1996). More recent data have suggested that there may be a dissociation of inhibitory potencies for molecules against these proteins. Although some molecules can interact with CYP3A4 and P-gp to a similar extent, for the most part the potency of inhibition for CYP3A4 did not predict the potency of inhibition for P-gp and vice versa (Wandell et al., 1999a. Moreover, not all CYP3A substrates, such as midazolam and nifedipine, are P-gp substrates (Kim et al., 1999). Thus, although there seems to be a relationship between the active sites of CYP3A4 and P-gp, it is not absolute.
To account for the observed broad substrate specificities for both CYP3A4 and P-gp, the presence of multiple drug binding sites has been proposed (Greenberger et al., 1990; Ekins et al., 1998; Korzekwa et al., 1998; Houston and Kenworthy, 2000). The first elegant signs of a complex behavior determined experimentally for P-gp appeared in 1996, when cooperative, competitive, and noncompetitive interactions between modulators were found to interact with at least two binding sites in P-gp (Ayesh et al., 1996). The multiple site hypothesis was confirmed by other groups (Dey et al., 1997; Scala et al., 1997; Shapiro and Ling, 1997). Subsequent results have indicated there may be three or more binding sites (Shapiro et al., 1999). Steady-state kinetic analyses of P-gp mediated ATPase activity using different substrates indicates these sites can show mixed-type or noncompetitive inhibition indicative of overlapping substrate specificities (Wang et al., 2000). Other researchers have determined that immobilized P-gp demonstrates competitive behavior between vinblastine and doxorubicin, cooperative allosteric interactions between cyclosporin and vinblastine or ATP, and anticooperative allosteric interactions between ATP, vinblastine, and verapamil (Lu et al., 2001). Clearly allosteric behavior by multiple substrates, inhibitors, or modulators of CYP3A4 or P-gp complicates predicting the behavior and drug-drug interactions of new molecules in vivo and has important implications for drug discovery.
Recently we have described how computational approaches can be used to predict inhibition of P-gp–mediated transport of digoxin, calcein accumulation, vinblastine accumulation, and vinblastine binding (Ekins et al., 2002). From this work, we hypothesized that vinblastine and digoxin are likely to bind a single site because strong correlations between the two models were observed. In the present study, the utility of the newly-derived computational models was tested further using available literature data on verapamil–P-gp binding in vinblastine-induced Caco-2 cells (Neuhoff et al., 2000). In addition, a new P-gp inhibition pharmacophore was constructed using the verapamil binding data and tested with the previously generated data. Finally, a P-gp substrate pharmacophore was produced using an alignment of verapamil and digoxin onto which vinblastine was fitted. Analyses of data using these computational models confirms that verapamil may interact with the vinblastine/digoxin binding site(s) in P-gp.
Materials and Methods
In Vitro Studies.
In vitro procedures for inhibition of P-gp using digoxin, calcein, and vinblastine as substrate probes with different in vitro systems have been described previously in detail (Ekins et al., 2002). Other experimental data was obtained from the literature for inhibition of racemic [3H]verapamil binding to P-gp in vinblastine induced Caco-2 cells (Neuhoff et al., 2000).
Molecular Modeling.
The computational molecular modeling studies were carried out using Silicon Graphics Octane and O2 workstations. Molecular structures were used as SMILES (Simplified Molecular Input Line Entry System) string format (Weininger, 1988) or imported from the ISIS MDDR-3D database (version 2000.2; MDL information Inc., San Leandro, CA).
Modeling with Catalyst.
Briefly, models were constructed as described previously (Ekins et al., 2002) using Catalyst version 4.5 (Molecular Simulations, San Diego, CA) after importing the molecular structures for molecules (Fig. 1) used in the in vitro studies (Ekins et al., 2002) and from the literature (Wandell et al., 1999a Neuhoff et al., 2000). The three-dimensional molecular structures were generated as described previously for CYP3A4 (Ekins et al., 1999a). In this study, the data for 16 inhibitors of verapamil binding (Table 1) were each treated as follows. For each molecule in either the training or test sets, the number of conformers generated using the ‘best’ functionality for each inhibitor was limited to a maximum number of 255 (with an energy range of 20 kcal/mol). Ten hypotheses were generated using these conformers for the 16 inhibitor training sets and the IC50 values. This was possible after selection of the following features for the inhibitors; hydrogen bond donor, hydrogen bond acceptor (defined by Catalyst as nonbasic amines that have a lone pair, sp or sp2 nitrogens, sp3 oxygens or sulfurs, and sp2 oxygens), hydrophobic (aromatic rings and aliphatic chains), and ring aromatic. After assessing all 10 hypotheses generated, the lowest energy-cost hypothesis was considered the best, because it possessed features representative of all the hypotheses. The total energy cost of the generated pharmacophores can be calculated from the deviation between the estimated activity and the observed activity, combined with the complexity of the hypothesis (i.e., the number of pharmacophore features). A null hypothesis can also be calculated that presumes that there is no relationship in the data and the experimental activities are normally distributed about their mean. Hence, the greater the difference between the energy cost of the generated hypothesis and the energy cost of the null hypothesis, the more likely it is that the hypothesis does not reflect a chance correlation.
Molecules used for modeling with an LY number.
Inhibition of verapamil binding (Neuhoff et al., 2000) predicted using previously generated pharmacophore models
The quality of the structure activity correlation between the estimated and observed activity values was estimated by means of an r2 value. Statistical significance of the retrieved hypothesis was verified by permuting (randomizing) the response variable 10 times (i.e., the activities and structures of the training set compounds were mixed 10 times), so that each value was no longer assigned to the original molecule, and the Catalyst hypothesis generation procedure was repeated.
Structural relationships between the P-gp substrates verapamil, vinblastine, and digoxin were initially assessed using the common features function (HIPHOP(Ekins et al., 1999b)] in Catalyst. Ultimately, however, vinblastine required fitting to the result of the digoxin and verapamil alignment. This alignment was then used with fast-fit and best-fit procedures to assess how the substrates might orient within the P-gp binding site. Fast-fit refers to the method of finding the optimum fit of the substrate to the hypothesis among all the conformers of the molecule without performing an energy minimization on the conformers of the molecule. The best-fit procedure starts with fast-fit and allows individual conformers to flex over a given energy threshold of 20 kcal/mol. This allows examination of more conformational space and minimizes the distance between the center of the hypothesis features and their mapping to atoms on the molecule. The result of the best-fit procedure is a conformer whose energy is generally higher than the starting one (Catalyst tutorials release 3.0; MSI, San Diego, CA.).
Validation of the Catalyst Models.
The various test sets contained molecules with IC50 orK i values not included in the initial training sets as described previously. These test set molecules were fit by the fast-fit algorithm to the respective Catalyst models to predict an IC50 value as described previously for CYP3A4 (Ekins et al., 1999a).
Statistical Evaluation of Test Set Predictions.
Observed and predicted inhibition data were graphed (data not shown) and fit using Excel 97 (Microsoft, Redmond, WA) to generate an r2 value. In all cases, the few molecules that were present in the training set were excluded from the subsequent correlation analysis. This data was also analyzed using the nonparametric Spearman's Rho test available in JMP 4.0.2 (SAS Institute Inc., Cary, NC). This test represents a correlation coefficient computed on the relative rank of the data values and not the values themselves.
Results
To elucidate likely binding sites of P-gp, the present study used a computational approach to model in vitro data obtained from the literature and our own laboratories. The Catalyst pharmacophore method generates models of the position of important ligand features in three-dimensional space, which may ultimately relate to features within P-gp itself. These three-dimensional models were achieved by assessing multiple conformations of each ligand alongside the experimental inhibition data. The end result is a computational model that can be used to predict how other molecules might inhibit P-gp solely by providing the structure to the software, which then attempts to fit the numerous conformations of the molecule as close to the centroids of the pharmacophore features as possible. Computational models derived from inhibition of digoxin transport in Caco-2 cells, vinblastine, and calcein accumulation in P-gp–expressing LLC-PK1 cells as well as vinblastine binding in vesicles derived from CEM/VLB100 cells (Ekins et al., 2002) were used to rank the experimental data for inhibition of verapamil binding in Caco-2 cells (Neuhoff et al., 2000). Additionally, these same data sets were used to test the pharmacophore built using literature data for the inhibition of verapamil binding in Caco-2 cells (Neuhoff et al., 2000).
The previously generated digoxin pharmacophore (Ekins et al., 2002) was assessed; the data from the inhibition of verapamil binding was derived from a recent publication (Neuhoff et al., 2000) (Table 1). The digoxin pharmacophore model provides a good rank ordering of P-gp inhibition; all predictions are within a log cutoff of the observed value. When these predictions are analyzed, the r2 value was 0.55 and the Spearman's Rho coefficient of 0.87 (p < 0.0001) suggesting the previously generated Catalyst model for inhibition of digoxin transport produces a significant rank order for these 15 molecules. The Catalyst pharmacophore for inhibition of vinblastine binding was less predictive of the verapamil binding inhibition data; eight molecules fell outside the log unit cutoff and r2 value was 0.15 (Table 1). However, the Spearman's Rho coefficient of 0.70 (p < 0.0034) indicated a significant rank ordering. The vinblastine accumulation pharmacophore (Table 1) predicted the data well; only two molecules were outside the log unit cutoff, yielding an r2value of 0.76 and a Spearman's Rho coefficient of 0.96 (p < 0.0001). The calcein accumulation pharmacophore (Table 1) also predicted the data well, with only two molecules predicted outside the log cutoff, yielding an r2value of 0.72 and a Spearman's Rho coefficient of 0.89 (p < 0.0001).
Next, a P-gp model was built using Catalyst with the inhibition of verapamil binding data (Neuhoff et al., 2000). These data generated a pharmacophore with 4 features: two hydrophobic, one hydrogen bond acceptor, and one ring aromatic feature (Fig.2). The correlation of observed and predicted data possessed an r2 value of 0.96 and fitted the inhibitors well (Fig. 3). The total energy cost of this pharmacophore was 75.2 units, lower than the null cost (101.1 units). Permuting the activity and structures for the published data for inhibition of verapamil binding 10 times resulted in a lower r2 value (0.72) and higher energy cost (87.4).
P-gp verapamil binding inhibition pharmacophore;n = 16 data points from Neuhoff et al. (2000)(orange, ring aromatic feature; green, hydrogen bond acceptor; and cyan, hydrophobic features).
Molecule Lu-47474 fitted to the P-gp verapamil binding inhibition pharmacophore; n = 16 data points from Neuhoff et al., (2000) (orange, ring aromatic feature; green, hydrogen bond acceptor; and cyan, hydrophobic features).
The inhibition of verapamil binding pharmacophore seemed valuable in predicting the inhibition of vinblastine binding (Table2). Although 14 of 21 molecules were predicted outside the log unit cutoff, an r2value of 0.89 and a Spearman's Rho coefficient of 0.55 (p = 0.011) were obtained. Predictions by the inhibition of verapamil binding pharmacophore for the inhibition of digoxin transport data set exceeded the log unit cutoff in most cases (Table 3), yielding an r2 value of 0.28 and a Spearman's Rho coefficient of 0.64 (p = 0.0005). Predictions by the inhibition of verapamil binding pharmacophore for the inhibition of vinblastine accumulation data (Table 4) were comparable with the digoxin transport data [r2 = 0.28 and Spearman Rho coefficient of 0.68 (p = 0.0024)], even though nine molecules were poorly predicted. The ability of the inhibition of verapamil binding pharmacophore to rank the calcein accumulation data (Table 4) was considerably poorer, because the r2 = 0.019 and the Spearman Rho coefficient was 0.4 and not statistically significant. However, this latter model had only four molecules that exceeded the log unit cutoff for the prediction error.
Test set molecules for inhibition of vinblastine binding predicted using inhibition of the verapamil binding model
Digoxin transport IC50 predictions with the inhibition of verapamil binding model
Vinblastine accumulation and calcein accumulation IC50predictions with the inhibition of verapamil binding model
When all five Catalyst pharmacophores generated by our laboratories were used to rank the inhibition of P-gp–mediated verapamil transport by verapamil metabolites (Table 5), all models predicted D-703 and norverapamil as more potent than D-617 and D-620, in keeping with experimental data (Pauli-Magnus et al., 2000). All five P-gp inhibitor pharmacophores generated in this and a previous study were also merged in Catalyst in an attempt to uncover features that could occupy similar regions in space (Fig.4). This analysis suggested the presence of at least four distinct groupings of features, consisting of two hydrophobic domains at the extremes of the figure along with a hydrogen bond acceptor region and ring aromatic region, both near one of the hydrophobic domains.
Prediction of the extent of inhibition of P-gp by the major metabolites of verapamil using pharmacophore models described in this and a previous study
Five merged P-gp inhibition pharmacophores defining common areas of hydrophobicity, hydrogen bonding, and ring aromatic features (orange, ring aromatic feature; green, hydrogen bond acceptor; purple, hydrogen bond donor; and cyan, hydrophobic features)
The common features of pharmacophore alignment of the P-gp substrates verapamil and digoxin and subsequent fitting of vinblastine to this alignment revealed numerous conserved pharmacophore features (Fig. 5, A and B). The fast alignment (Fig. 5A) suggested a less extended conformation for digoxin than using the best-fit approach (Fig. 5B). All three substrates were predicted to possess multiple hydrophobic and hydrogen bond acceptor features at their extremities, pointing to these as important characteristics of substrates binding at this site within P-gp.
Fast fit (A) and best fit (B) to the Catalyst common features model for P-gp substrates. Molecules colored as follows: blue, digoxin; green, verapamil; red, vinblastine. Pharmacophore features colored green, hydrogen bond acceptor; and cyan, hydrophobic features.
Discussion
We have shown previously that structurally diverse P-gp inhibitors can be used to generate three-dimensional quantitative structure-activity relationship models that significantly rank the ability of agents not in the training set to inhibit P-gp (Ekins et al., 2002). Subsequently, these models provide us with a means to computationally predict whether a molecule is likely to interact with P-gp in vitro from its three-dimensional structure alone. Our four distinct pharmacophores derived from different in vitro systems contained a combination of hydrogen bond acceptors, hydrogen bond donors, hydrophobes, and ring aromatic features, and most displayed cross predictivity (Ekins et al., 2002). This work also indicated some degree of overlap for the binding site(s) probed by digoxin and vinblastine. In the present study, we have extended this work using literature data for inhibition of verapamil binding, to produce a fifth P-gp pharmacophore that has been used to predict the other four data sets (Neuhoff et al., 2000). In addition, we used this data and other published data (Pauli-Magnus et al., 2000) to test our models to date. In the process of this work, we have suggested that it is likely that some of the P-gp substrates used may bind at a similar site and this was evaluated using an alignment of substrates.
All four previously generated pharmacophores for P-gp inhibition were useful in producing statistically significant rank ordering of the verapamil binding inhibition data used in this study (Table 1). The Catalyst model built with these 16 P-gp inhibitors of verapamil binding was also predictive using the Spearman's Rho rank-ordering of data generated in our previous work based on inhibition of digoxin transport, vinblastine binding, and vinblastine accumulation (Tables2-4). Our findings suggest that vinblastine, verapamil, and digoxin have an overlapping affinity for similar or identical binding site(s) within P-gp. Merging all five pharmacophores derived to date suggest the presence of clearly defined regions in space occupied by clusters of identical features such as hydrophobes, hydrogen bond acceptors, and ring aromatic features (Fig. 4). This in itself is interesting, because this figure represents the sum of different substrate probes and experimental systems for evaluating P-gp. To some extent, the inhibition pharmacophore for calcein accumulation overlaps with the other four inhibition pharmacophores.
The hypothesis that some of the substrate probes for P-gp bind similar sites was tested further by aligning verapamil and digoxin to assess likely common features followed by fitting vinblastine to this substrate model (Fig. 5, A and B). All three molecules generated extended conformations; hydrophobic and hydrogen bond acceptor features were common to all. These in turn were distributed toward the extremities of the aligned molecules. Therefore, this P-gp substrate pharmacophore (Fig. 5, A and B) seems to be consistent with the features present in some of the inhibitor pharmacophores generated previously (Ekins et al., 2002) and for verapamil in this study (Figs.2 and 3). The alignment of other known P-gp modulators such as trimethoxybenzoylyohimbine (TMBY) and reserpine was evaluated (Fig.6) with the P-gp substrate model derived from the verapamil and digoxin alignment. Both TMBY and reserpine molecules fit to multiple features in the extended conformations; these were found to be similar to verapamil in the case of reserpine. This finding is in agreement with previous work that had aligned reserpine analogs using different computational software (Pearce et al., 1989). In addition, this provides additional support to our suggestion that vinblastine, verapamil, and digoxin are likely to bind similar sites in P-gp. Our findings are also consistent with available data, which show TMBY and verapamil bind to a single or overlapping site in a human leukemic cell line (Shepard et al., 1998) and that TMBY is a competitive inhibitor of vinblastine binding to P-gp (Dantzig et al., 1996). In addition, LY335979 has been shown to competitively block vinblastine binding (Dantzig et al., 1996) and vinblastine competitively inhibits verapamil stimulation of P-gp–ATPase (Shepard et al., 1998). Our recent in vitro and in silico approaches seem to confirm the chemical features important for binding this site as either substrates or inhibitors.
Best fit alignment with the P-gp common features substrate model, molecules colored as follows: green, verapamil; yellow, TMBY; purple, reserpine; pharmacophore features colored green, hydrogen bond acceptor; and cyan, hydrophobic features.
When the P-gp models in this and previous studies are compared with pharmacophores for inhibitors and substrates of CYP3A (Ekins et al., 1999a,b), there is some similarity between these pharmacophores (Ekins et al., 2002). Both proteins share pharmacophores with multiple hydrophobic features and at least one hydrogen bond acceptor feature but in slightly different arrangements. This might explain why potent CYP3A4 inhibitors are not necessarily potent P-gp inhibitors (Wandell et al., 1999a or substrates (Kim et al., 1999). This could also help clarify the independent nature of the relationship between ligands of these two proteins.
The utility of the P-gp inhibitor pharmacophores produced to date can be assessed by using them to predict the rank order of potential P-gp inhibitors not present in the models. A recent publication described data for the varying potency of P-gp inhibition of verapamil metabolites (Pauli-Magnus et al., 2000). We used each pharmacophore model to predict the likely IC50 of the four verapamil metabolites D-617, D-620, D-703, and norverapamil (Table 5). All five pharmacophores were able to correctly identify the least potent P-gp inhibitors as D-617 and D-620. These molecules are more likely to be hydrophilic as suggested previously for low-affinity inhibitors of binding at the verapamil site (Neuhoff et al., 2000). Two of the pharmacophores based on inhibitors of digoxin or verapamil binding provided useful quantitative predictions for the inhibition of verapamil transport with verapamil metabolites. The digoxin pharmacophore performed particularly well with the potent inhibitors D-703 and norverapamil that fit the majority of the features in the models. The third model based on inhibitors of vinblastine binding failed with D-703 in that the prediction exceeded 1 log order of magnitude. These predictions for verapamil metabolites are perhaps not totally unexpected in that all three of these models used verapamil as a training set member. However, the digoxin model used a much higher IC50 value for verapamil making it one of the less potent inhibitors, which might explain its lesser predictive ability in this case. Our models for inhibition of vinblastine and calcein accumulation differed in that the latter, although generating the correct rank order, was less quantitative in the nature of the predictions.
In conclusion, we have confirmed our previous suggestion that digoxin and vinblastine are likely to bind a similar or overlapping P-gp binding sites by using data from inhibition of verapamil binding as an additional test case. It seems increasingly likely that some P-gp substrates possess common structural features dominated by hydrophobic features as well as hydrogen bond acceptor features. These features are naturally present on many of the P-gp inhibitors evaluated in this and other studies (Ekins et al., 2002). The correct alignment of these molecular features with the P-gp binding site will result in a potent inhibitory interaction. We now have a computational means to predict likely inhibition of P-gp substrate probes, and undoubtedly, at some point, quantitatively predictive substrate models will be generated that will expand upon the simple substrate alignment model presented in this study. With the understanding that there may be multiple binding sites within P-gp, we have shown that pharmacophores represent useful tools for classifying and characterizing substrate-binding site interactions. We have also shown this technique has the potential to enhance our understanding of a complex transporter such as P-gp in the absence of a crystal structure.
Acknowledgments
S.E. acknowledges Dr. Homer L. Pearce for stimulating discussions in the early stages of this work.
Footnotes
- Received November 2, 2001.
- Accepted February 13, 2001.
-
This work was supported in part by United States Public Health Service grant GM31304 (to R.B.K.) and National Institutes of Health grant ES08658 (to E.S.).
Abbreviations
- P-gp
- P-glycoprotein
- TMBY
- trimethoxybenzoylyohimbine
- LY335979
- 4-(1,1-difluoro-1,1a,6,10b-tetrahydrodibenzo[a,e]cyclopropa[c]cyclohepten-6-yl)-[(5-quinolinyloxy)methyl]-1-piperazineethanol
- The American Society for Pharmacology and Experimental Therapeutics
