2D- and 3D-QSAR studies of a series of benzopyranes and benzopyrano[3,4b][1,4]-oxazines as inhibitors of the multidrug transporter P-glycoprotein

The ATP-binding cassette efflux transporter P-glycoprotein (P-gp) is notorious for contributing to multidrug resistance in antitumor therapy. Due to its expression in many blood-organ barriers, it also influences the pharmacokinetics of drugs and drug candidates and is involved in drug/drug- and drug/nutrient interactions. However, due to lack of structural information the molecular basis of ligand/transporter interaction still needs to be elucidated. Towards this goal, a series of Benzopyranes and Benzopyrano[3,4b][1,4]oxazines have been synthesized and pharmacologically tested for their ability to inhibit P-gp mediated daunomycin efflux. Both quantitative structure–activity relationship (QSAR) models using simple physicochemical and novel GRID-independent molecular descriptors (GRIND) were established to shed light on the structural requirements for high P-gp inhibitory activity. The results from 2D-QSAR showed a linear correlation of vdW surface area (Å2) of hydrophobic atoms with the pharmacological activity. GRIND (3D-QSAR) studies allowed to identify important mutual distances between pharmacophoric features, which include one H-bond donor, two H-bond acceptors and two hydrophobic groups as well as their distances from different steric hot spots of the molecules. Activity of the compounds particularly increases with increase of the distance of an H-bond donor or a hydrophobic feature from a particular steric hot spot of the benzopyrane analogs. Electronic supplementary material The online version of this article (doi:10.1007/s10822-013-9635-9) contains supplementary material, which is available to authorized users.


Introduction
Development of multidrug resistance (MDR) is one of the major challenges in cancer chemotherapy, as it limits the effectiveness of many clinically important agents [1]. One of the basic underlying mechanisms is overexpression of the mdr1 gene product, P-glycoprotein (P-gp) [2], which belongs to the ATP-binding cassette (ABC) family of transporters [3]. It is highly promiscuous in its ligand recognition profile and thus transports a large variety of structurally and functionally diverse compounds out of Electronic supplementary material The online version of this article (doi:10.1007/s10822-013-9635-9) contains supplementary material, which is available to authorized users. tumor cells [4]. Apart from its role in tumor cells it is expressed at epithelial cells of liver, kidney, intestine and colon, as well as at the blood brain barrier. Thus, P-gp not only plays an important role in maintaining a concentration gradient of toxic compounds at these physiological barriers, but also modulates the pharmacokinetics of drugs that are recognized as P-gp substrates.
Within the past decade numerous inhibitors of P-gp mediated drug efflux have been identified [3]. Several compounds entered even phase III clinical studies, such as MS-209 (dofequidar fumarate), tariquidar, valspodar and elacridar [5,6]. However, none made it to the market so far, mainly because of lack of efficacy or severe side effects. In light of our extensive SAR and QSAR studies of propafenones [7,8], benzophenones [9] and dihydrobenzopyrans [10], a new class of conformationally restricted benzopyrano [3,4b] [1,4]oxazines have been synthesized and biologically tested with respect to their ability to block P-gp mediated daunomycin efflux. These new P-gp inhibitors offer the advantage of a remarkably reduced conformational flexibility, which renders them versatile molecular tools for probing stereoselective differences of drug/P-gp interaction [11], as well as for 3D-QSAR studies. These might be performed by utilizing alignment-dependent approaches, such as CoMFA and CoMSIA, or by alignment independent methods using descriptors derived from Molecular Interaction Fields (MIFs), like the GRIND [12]. In particular, the latter allow the analysis of structurally diverse data series. GRID MIFs [13] have been applied to many areas of computational drug discovery, including 3D-QSAR [14], docking [15], high-throughput virtual screening [16], ADME profiling, kinetic [17,18] and metabolism prediction [19] of early drug candidates. In this manuscript we explore the capability of the GRIND approach to derive predictive 3D-QSAR models for a set of diastereomeric benzopyrano [3,4b] [1,4]oxazines. The GRIND based 3D-QSAR models added value in recognition of important pharmacophoric features and their mutual distances. In addition, molecular shape of the P-gp inhibitors has been recognized as an important structural prerequisite for high pharmacological activity.

Chemistry
Synthesis of the benzopyrane common scaffold was achieved in analogy to the procedure reported by Godfrey et al. [20] and following our strategy outlined recently [11]. Besides a set of diastereomeric esters and benzopyrano [3,4b] [1,4]oxazines, also a set of corresponding ethers were synthesized. Respective procedures and experimental details are provided in the supplementary material. In general, compounds showing 4aS,10bR stereochemistry are denoted as (a)-series, whereas the respective 4aR,10bS-analogues are assigned as (b)-series ( Table 1).

Calculation of physicochemical parameters
Hansch analysis Molecular descriptors supplied by the program MOE (atom and bond counts, connectivity indices, partial charge descriptors, pharmacophore feature descriptors, logP (o/w), calculated physical property descriptors) were computed for Hansch analysis. QSAR-Contingency [21], a statistical application in MOE, was used for the selection of relevant descriptors. PLS analysis was performed to determine the relationship between these 2D molecular descriptors and biological activity of the compounds. The predictive ability of the model was determined by classical leave one out (LOO) and leave one pair out cross validation procedures (SM Table 1). In order to remove any bias, the final model was externally validated by using a test set of already published dihydrobenzopyrans [10]. GRIND 3D conformations of the molecules in the data set were obtained from their 2D coordinates by using the program CORINA [22]. Molecular Interaction Fields (MIF) were calculated as GRID based fields in Molecular Discovery software Pentacle [23] using four different probes: DRY probe to represent hydrophobic interactions, O sp 2 carbonyl oxygen probe to represent H-bond donor feature of the molecules, N1 probe to represent -NH which is a neutral flat probe as an H-bond acceptor in the molecules and the TIP probe that represents the shape of the molecule, in terms of steric hot spots. The regions with the most relevant MIF were extracted by applying the AMANDA algorithm [24] that uses the intensity of the field at a node and the mutual nodenode distances between the chosen nodes. At each point, the interaction energy (Exyz) was calculated as a sum of Lennard-Jones energy (Elj), Hydrogen bond (Ehb) and Electrostatic (Eal) interactions.
Default values of probe cutoff (DRY = -0.5, O = -2.6, N1 = -4.2, TIP = -0.74) was used for discretization of MIF. Nodes with an energy value below this cutoff were discarded. The Consistently Large Auto and Cross Correlation (CLACC) algorithm [23] was used for encoding the prefiltered nodes into GRIND thus producing most consistent variables as compared to MACC [25]. The values obtained from this analysis were represented directly in correlogram plots, where the product of node-node energies is reported versus the distance separating the nodes. Highest energy product can be defined for the same probe  cross validation as already described in the 2D-QSAR section (SM Table 2) as well as by an external test set composed of previously published compounds.

Pharmacology
Biological activity of target compounds 5a-22b was assessed using the daunorubicin efflux protocol as described previously [26]. Briefly, multidrug resistant CCRF-CEM vcr 1,000 cells were incubated with daunorubicin and the decrease in mean cellular fluorescence in dependence of time was measured in presence of various concentrations of the modulator. IC 50 values were calculated from the concentration-response curve of efflux V max /K m versus concentration of the modulator. Thus, the effect of different modulators on the transport rate is measured in a direct functional assay. Values are given in Table 1 and are the mean of at least three independently performed experiments. Generally, inter experimental variation was below 20 %.

Results and discussion
Structure activity relationships (SAR) Biological activity values of the data series cover a range of more than three orders of magnitude (Table 1) with the two phenylalanine esters 7a and 7b being the most active compounds (7a: 0.55 lM; 7b: 0.77 lM), followed by N-methylated L-valine analogues 9a (0.96 lM) and 9b (1.35 lM), which are by a factor of 2 more active than the corresponding unsubstituted analogs 6a (2.40 lM) and 6b (2.70 lM). The same trend could be observed for the respective D-valine derivatives. This observation is even more pronounced for the alanine derivatives (compare methylated analogs 8a (3.96 lM) and 8b (3.72 lM) versus respective secondary amines 5a (29.85 lM) and 5b (14.55 lM). This most probably is due to a logP effect with more lipophilic compounds showing higher biological activity, which has been shown for numerous classes of P-gp inhibitors [27]. It has to be noted that for all seven diastereoisomeric pairs showing a bicyclic scaffold almost no differences in biological activity exist. However, this pattern changes remarkably upon ring closure to the tricyclic benzopyrano [3,4b] [1,4]oxazines. While all stereoisomers containing a valine moiety (13a,b; 16a,b, 18a,b, 19a,b-22a,b) are still within one order of magnitude, both the alanine and phenylalanine derivatives exhibit remarkable differences in their P-gp inhibitory potency (IC 50 ) values. Interestingly, in case of alanine, the 4aS,10bR-isomer 12a is by a factor of 15 less active than the diastereomeric 4aR,10bS analogue 12b, whereas in case of the phenylalanine derivatives this behavior reverses with the 4aS,10bR-isomer 14a being by two orders of magnitude more active than 14b. This difference in their biological activities might be due to difference in mode of interaction of diastereoisomeric pairs as has been indicated in a preceding publication [11].
Hansch analysis 3D structures of all diastereoisomers were built with the builder function of MOE 2011. 10 and energy minimised using the MMFF94 force field which uses a bond charge increment method to set the electrostatic partial charges [28]. In order to determine the influence of physicochemical properties of the compounds on their biological activity, QSAR analyses were performed by using the software package MOE version 2011. 10 and MOE's contingency analysis tool for identification of the most important descriptors. The multiple linear regression analysis produced an equation solely based on the hydrophobic van der Waals surface area (vsa_hyd) (Eq. 1). Interestingly, descriptors related to electrostatic properties, such as topological polar surface area and molar refractivity, did not show significant contributions to the model. series additional factors other than lipophilicity might play a role. Vsa_hyd describes the sum of van der Waals surface areas of hydrophobic atoms (Å 2 ). This is perfectly in line with previous studies which showed that distribution of hydrophobicity within the molecules influences their mode of interaction with P-gp [29] and lipophilicity needs to be considered as a space directed property [30,31]. This space-directedness might be indicative for different orientations of molecules within the binding area of P-gp, which is mainly hydrophobic [32]. The QSAR model was further validated by using an external test set of already published dihydrobenzopyrans and tetrahydroquinolines [10]. All compounds are predicted well, with the residuals being less than one log unit from their experimental inhibitory potencies (log IC 50 ). This further strengthens the reliability of the final QSAR model (Table 2). Additionally, 18 different models were developed by taking one pair of diastereoisomer out at each step. All models showed q 2 values in the range of 0.57-0.70, which further demonstrates the consistency of the QSAR model (SM Table 1).

GRID Independent molecular descriptor (GRIND) analysis
The previously computed molecular structures along with their activity values (expressed as log1/IC 50 ) were loaded into the software package Pentacle (v 1.06) [23] to derive 3D-QSAR model using GRIND descriptors. According to previous findings for propafenone analogs, all compounds were modeled in their neutral form [33]. Structural variance of the data was analyzed with principal component analysis (PCA) performed on the complete set of GRIND descriptors. The first two principal components explain about 32 % of the descriptor variance in the data set. Principal component analysis (PCA) on the data matrix showed that the series is organized in three different clusters (Fig. 2). Molecules in the cluster on the right hand side (cluster 1) do not contain any H-bond donor, while the second cluster (cluster 2) contains one H-bond donor group. The 3rd cluster (cluster 3), located on the upper left corner of the plot, contains compounds with two H-bond donor groups in their structures. Furthermore, rigid and smaller compounds (cluster 1 and cluster 2) are separated from the flexible ones (cluster 3). Overall, compounds in cluster 3 are more potent than compounds in cluster 1 and 2, suggesting that an elongated structure is an important perquisite for high P-gp inhibitory potency. In order to identify the more important pharmacophoric features of ligand-protein interaction, a PLS model was built, using the complete set of active variables (450) generated by Pentacle (v 1.06). This resulted in a one-latent variable (LV1) model with an r 2 of 0.51 and a cross-validated (LOO) q 2 value of 0.27, which was quite unsatisfactory. Thus, variable selection was applied to reduce the variable number using the FFD variables selection algorithm [34] implemented in Pentacle. This resulted in a decrease from 450 to 196 variables and a large increase of the model quality (r 2 of 0.72, q 2 of 0.58, standard error of prediction 0.52).  With the exception of compounds 14b, 5a, and 12a, all compounds are within one order of magnitude from their predicted values (14b: obs 259.78, pred 23.21; 5a: obs 29.85, pred 2.80; 12a: obs 1241.65, pred 58.30 lM) (Fig. 3). The outlier behavior of these three compounds might be due to potential different interaction behavior of the two diastereomeric series as reported by Jabeen et al. [11]. However, building two separate QSAR models composed of compounds of series (a) and series (b) in two separate training sets showed an analogous picture and did not improve the results (data not shown). Thus, although GRIND descriptors are able to capture different configurations, they were not able to extract the differences of the two diastereomeric series. This might be due to the fact that the molecules are quite compact (Fig. 4) and GRIND is considering MIFs within a grid step of 0.5 Å .
All compounds of the external test set are predicted within one log unit from the experimental inhibitory potencies (log IC 50 ), except (1c), where the residual is slightly more than one log unit ( Table 2). The low activity of 1c mainly might be due to its low logP value, which is not properly reflected in the GRIND based pharmacophoric features. Thus, GRIND is over predicting the compound. The overall good predictive ability and model statistics of all 18 leave one pair out GRIND models further demonstrates the consistency and validity of the GRIND based 3D-QSAR model (SM Table 2).
Analysis of the PLS coefficients profile of the GRIND model allows to identify those descriptors which exhibit the largest contribution to the model. According to the bar plot shown in Fig. 5, certain distances of the N1-N1, O-N1, and O-TIP probes are participating most in explaining the variance in the biological activity values ( Table 3).
The sum of the van der Waals surface areas of hydrophobic atoms (vsa_hyd) has emerged as an important determinant for high biological activity of benzopyranetype P-gp inhibitors (Eq. 1). The 3D-QSAR model using GRIND descriptors further refines this general property and identified two hydrophobic regions (DRY-DRY) separated by a certain distance range in all active compounds. These represent the aromatic ring of the benzopyrane ring system and R1. In the most active phenylalanine derivatives (7a,b and 14a,b) the two regions are separated by a distance of 13.2-13.6 Å , which is considered optimal according to the GRIND model. Thus, adding a large hydrophobic group (large vsa_hyd) at the position of R1 might lead to a further increase of the biological activity.
Previous QSAR studies on propafenone derivatives have demonstrated the importance of H-bond acceptors and their distance from the central aromatic ring [35,36].  acceptors at a distance of 8.0 Å apart from each other in P-gp inhibitors [39], whereas a distance of 11.5 Å between two H-bond acceptors, along with the importance of shape descriptors, have been reported by Cianchetta et al. [40] for substrates of P-gp. However, despite of similarities in the number of H-bond acceptors necessary for a high biological activity, a direct comparison of distance matrices thereof across different chemical scaffolds reveals some differences. This most probably is due to the fact that the large binding site of P-gp has multiple spots able to contribute to H-bond interactions and that different chemical series utilize different H-bond interaction patterns.
Apart from a certain number of H-bond acceptors, one H-bond donor along with hydrophobic regions distribution have been identified as important pharmacophoric features of P-gp inhibitors/substrates [8,36]. It is worth noting that a very similar MIF based pharmacophore of P-gp inhibitors was recently published by Broccatelli et al. [41]. They identified one H-bond acceptor and two large hydrophobic regions, together with an optimal molecular shape, as being important for high activity, and successfully used their model for virtual screening to identify new P-gp inhibitors. The results are further in line with Boccard et al. [42] outlining an optimal shape and hydrophobicity as major physicochemical parameters responsible for the affinity of flavonoid derivatives for P-gp [43,44].
Also in our GRIND model, shape based probes (TIP) defining steric hot spots exhibit a significant contribution. Especially the 9-carbonitrile group in the benzopyrane scaffold encodes an important molecular boundary (steric hot spot) and serves as anchor for defining optimal distance ranges to an H-bond donor (O-TIP correlogram) as well as to a hydrophobic feature (DRY-TIP correlogram). The O-TIP combination of probes encodes the shape of the molecules (steric hot spots) together with an H-bond donor group. Interestingly, O-TIP coefficients are negative for a distance between 5.6 and 6.0 Å , but become positive for larger distances (12.8-13.2 Å ). These distances (12.8-13.2 Å ) are present in benzopyranes bearing tert-butyl esters (5a-11b) and amino alcohol derivatives (19a-20b and 22a,b) as shown in Fig. 6c. In tricyclic diastereoisomers (12a-14b and 17b) these descriptors are linked to shorter distances and mark (-NH) as an H-bond donor at a distance of 5.6-6.0 Å apart from the cyano group, which is the main group contributing to the TIP MIF, and seems to be related with a negative influence for the biological activity of this group (Fig. 6d). This indicates that, in general, the most potent P-gp inhibitors show extended conformations and have an H-bond donor group far from regions with a strong TIP probe related field.
Analyzing the DRY-TIP correlogram it becomes evident that a hydrophobic group at a distance of 15.2-15.6 Å from one of the ''edges'' of the molecule (steric hot spot, cyano group) positively contributes to biological activity. In tert-butyl esters (5a-11b) and 14a,b these two probes map the distance between a hydrophobic group (R 1 ) (14a,b) or tert-butyl group in 5a-11b from the cyano group (Fig. 6b). In analogy to the O-TIP correlogram, DRY-TIP shows a negative contribution towards biological activity for shorter distances (7.6-8.0 Å ) of these probes.
Finally, the O-N1 correlogram (H-bond donor-H-bond acceptor) points towards two positive contributions at a distance of 9.6-10.0 and 2.4-2.8Å , respectively (Fig. 6e, f). The first distance is linked to the hydroxyl and carbonyl group in 5a-11b and is complementary to the N1-N1 correlogram as already discussed. The second distance refers to the -NH and carbonyl group. O-N1 probes at both distance ranges are well pronounced in tert-butyl esters (5a-11b) as well as in amino alcohol substituted derivatives (19a-20 b and 22a,b). However, in all tricyclic compounds (12a-18b) the two probes do not fit either of the distance ranges.
To summarize, the presence of two H-bond-acceptor groups and one H-bond donor at a particular distance from each other and from a particular ''edge'' or steric hot spot of the molecule is associated to an increase of the biological activity in benzopyrane-type P-gp inhibitors (Fig. 6).

Conclusions
Benzopyrano- [3,4b] [1,4]oxazines are versatile molecular tools for probing the stereoselectivity of P-glycoprotein. For a distinct substitution pattern, different pairs of diastereoisomers exhibit a large difference in their potency to inhibit P-gp mediated drug efflux. Unfortunately, GRINDbased 3D-QSAR models were unable to link these differences to concrete differences of distances between pharmacophoric hot spots, even if the GRIND analysis provided a reasonably well performing 3D-QSAR model outlining a set of important pharmacophoric features. Two H-bond-acceptor groups, one H-bond donor at a particular distance from each other as well as distinct distances of these probes to steric hot spots seem to play a major role in the interaction of benzopyrane-type P-gp inhibitors. The activity particularly increases when increasing the distance between an H-bond donor or a hydrophobic feature and a particular steric hot spot of the benzopyrane analogs. This not only further highlights the importance of H-bonding, but also indicates that a certain shape/configuration of the molecules is important for high activity. Further analyses will focus on a generalisation of this finding in other series of P-gp inhibitors. thanks the Higher Education Commission of Pakistan for financial support.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.