In Silico Quantitative Structure-Activity Relationship Studies on P-gp Modulators of Tetrahydroisoquinoline-Ethyl-Phenylamine Series

Background Multidrug resistance (MDR) is a major obstacle in cancer chemotherapy. The drug efflux by a transport protein is the main reason for MDR. In humans, MDR mainly occurs when the ATP-binding cassette (ABC) family of proteins is overexpressed simultaneously. P-glycoprotein (P-gp) is most commonly associated with human MDR; it utilizes energy from adenosine triphosphate (ATP) to transport a number of substrates out of cells against concentration gradients. By the active transport of substrates against concentration gradients, intracellular concentrations of substrates are decreased. This leads to the cause of failure in cancer chemotherapy. Results Herein, we report Topomer CoMFA (Comparative Molecular Field Analysis) and HQSAR (Hologram Quantitative Structure Activity Relationship) models for third generation MDR modulators. The Topomer CoMFA model showed good correlation between the actual and predicted values for training set molecules. The developed model showed cross validated correlation coefficient (q2) = 0.536 and non-cross validated correlation coefficient (r2) = 0.975 with eight components. The best HQSAR model (q2 = 0.777, r2 = 0.956) with 5-8 atom counts was used to predict the activity of test set compounds. Both models were validated using test set compounds, and gave a good predictive values of 0.604 and 0.730. Conclusions The contour map near R1 indicates that substitution of a bulkier and polar group to the ortho position of the benzene ring enhances the inhibitory effect. This explains why compounds with a nitro group have good inhibitory potency. Molecular fragment analyses shed light on some essential structural and topological features of third generation MDR modulators. Fragments analysis showed that the presence of tertiary nitrogen, a central phenyl ring and an aromatic dimethoxy group contributed to the inhibitory effect. Based on contour map information and fragment information, five new molecules with variable R1 substituents were designed. The activity of these designed molecules was predicted by the Topomer CoMFA and HQSAR models. The novel compounds showed higher potency than existing compounds.

Background MDR, the principal mechanism by which many cancers develop resistance to chemotherapy drugs, is a major factor in the failure of many forms of chemotherapy [1]. In MDR tumour cells, various member of the ABC family of transport proteins can simultaneously be overexpressed: these include P-gp (ABCB1), breast cancer resistance protein (BCRP, ABCG2) and MDR associated protein 1 (MRP1, ABCC family) [2]. These transporters utilize energy from ATP hydrolysis to transport a wide variety of substances out of cells against concentration gradients. The active efflux of substances from cells decreases their intracellular concentration and results in failure of chemotherapy.
Among the 49 identified human ABC transporters, P-gp is most intensively studied [3], and is a member of MDR/TAP (transporter associated proteins) subfamily. P-gp has ability to transport a wide variety of structurally unrelated substances out of cells [4][5][6]. Pgp is extensively distributed and expressed in the intestinal epithelium, hepatocytes, renal proximal tubular cells, adrenal gland and capillary endothelial cells comprising the blood-brain and blood-testis barrier. Pgp transports structurally diverse substrates and most are anticancer drugs such as doxorubicin, daunorubicin, paclitaxel, etoposid, teniposid, vinblastine and vincristine [7]. P-gp does not interact with anionic compounds but does interact with amphipathic compounds with molecular masses between 400-1900 daltons [8,9]. The calcium channel blocker verapamil can overcome MDR in cancer cells [10]. Another drug, cyclosporine-A, was designed as an immunosuppressant, but shows a promising P-gp inhibitory effect. Both these drugs are used as first generation P-gp inhibitors, also called MDR modulators. The use of these modulators has been limited because of low efficacy and higher dose-related toxicity. The second generation modulators, dexverapamil and PSC833, had higher efficacy and lower toxicity, but produce serious drugdrug interactions clinically.
Nowadays, a third-generation of MDR modulators are under investigation. These drugs, which include tariquidar (XR9576), zosuquidar (LY335979) and laniquidar (R1010933), possess selectivity, low toxicity and high efficacy [11]. These modulators are structurally different from the first-and second-generation modulators. Early in their evaluation, these modulators displayed promising activity. However, toxicity was subsequently observed [12]. The toxicities were found not to be mechanism based. These modulators showed some potential as new drugs, but were dropped due to toxicity related to high dose to be effective physiologically. Therefore, there is a clear need to enhance the activity of modulators that would also reduce the required dose. The search for new nontoxic, efficacious, potent modulators without drug-drug interactions has been intensive. The studies include a 3D-QSAR and free Wilson analysis on a series of tariquidar analogues as MDR modulators [13], a QSAR study on anthranilamide derivatives containing the nucleus of XR9576 and a 3D-QSAR study using tariquidar derivatives (WK-X and WK-Y compounds) to develop QSAR CoMFA/ CoMSIA models [14]. With the aim to find out important groups and atoms for P-gp antagonism, the present study selected a third-generation MDR modulator to develop Topomer CoMFA and HQSAR models [15]. This series consisted of same tetrahydroisoquinolineethyl-phenyl-amine nucleus present in two of the aforementioned series.

Data Set
The activity dataset, which was selected from reported literature [15], consisted of 39 molecules (Table 1). For analysis, the given inhibitory concentration values were changed to minus logarithmic scale value (pIC 50 ), as a dependent variable for Topomer CoMFA and HQSAR analysis by using the formula provided below. It is common to convert the biological activity data into a logarithmic scale, because the resulting model behaves The dataset was randomly partitioned into training and test set molecules by considering range of molecules (pIC 50 = 4.52-6.66), so that both the training and test sets consist of high, medium and low activity molecules. The training and test set consist of 31 and 8 molecules, respectively. All the molecules were built using the SYBYL 8.1 molecular modeling package [16]. All the dataset molecules were sketched by the SYBYL sketching program and were minimized by using the Tripos force field. They were then subjected to simulated annealing to get a stable conformation. Simulated annealing was performed for each ligand up to 200 cycles with default parameters, and then conformations were sorted according to the least potential energy value. These conformations were minimized with quantum mechanical semi-empirical AM1 method with precise convergence and full optimization commands with MMOK (Molecular Mechanics Correction to CONH Bonds) keywords. The dataset was then used for Topomer CoMFA and HQSAR analysis.

Topomer CoMFA
A Topomer CoMFA technique merges CoMFA [17] and topomer technology, to overcome the alignment problem of CoMFA [18]. Topomer CoMFA includes alignment of structural fragments. Structural fragments by definition contain a common feature, the "open valence" or "attachment bond". The Topomer methodology overlaps this common feature to provide an absolute orientation for any fragment. A Topomer is an invariant three-dimensional (3D) representation of molecular subunit generated from its two-dimensional (2D) topology by topomer alignment in topomer CoMFA [19]. In Topomer CoMFA analysis, all molecules of dataset were divided into two fragments, shown as R1 (blue) and R2 (red) groups in Figure 1. Each Topomer fragment was applied with topomer alignment to make a 3D invariant representation [20]. In Topomer CoMFA, atomic charges were calculated by the Gasteiger-Marsilli method for the topomer structure. Topomer CoMFA acts in two different ways for the calculation of molecular fields. An 'attenuation factor' reduces the field contributions of fragment atoms more distant from the attachment bond. Finally, the r 2 is calculated by using the same optimum number of component obtained from leave-one-out (LOO) cross validation analysis. Topomer CoMFA steric and electrostatic fields were calculated at a regular space grid of 2 Å, and were fixed automatically into a 1000 point cube to contain a Topomer. A sp 3 hybridized carbon atom was used as a probe atom for the steric field calculation and a negative oxygen atom was used as a probe for electrostatics field.

Partial Least Square (PLS)
The relationship between structural parameters and biological activities of compounds under study has been quantified using a PLS algorithm [21][22][23]. Topomer CoMFA descriptors were used as independent variable and biological activity (pIC 50 ) as a dependent variable. The cross-validation analysis was performed by using the LOO method, in which one molecule is removed from the dataset and its activity is predicted by using the model derived from rest of the molecules in the dataset. The q 2 resulted in an optimum number of components and the lowest standard error of prediction. The q 2 is calculated using the following equation: where γ pred , γ actual and γ mean are predicted, actual and mean values of the target property (pIC 50 ), respectively.

Predictive Correlation Coefficient (r 2
pred ) The predictive power of Topomer CoMFA and HQSAR models were derived from set of eight molecules, which were excluded during model development. In the structural preparation of test set molecules, sketching and optimization was same as the training set molecules.
The activity of the test set was predicted by using model derived from training set. The predictive correlation coefficient for developed model was determined by using following formula: where, PRESS is the sum of the squared deviation between the predicted and actual activity of the test set molecules, and SD is defined as the sum of the square deviation between the biological activity of the test set compounds and the mean activity of the training set molecules.

HQSAR
HQSAR is a technique that employs fragment fingerprints as predictive variables of biological activity or other structural related data [24]. HQSAR does not require a 3D structure of bioactive conformation or molecular alignments. HQSAR model generation deals with the 2D structure directed fragment fingerprints [25]. These molecular fingerprints are broken into strings at fixed intervals as specified by a hologram length (HL) parameter. The HL determines the number of bins in the hologram into which the fragments are hashed. The optimal HQSAR model was derived from screening through the 12 default HL values, which were a set of 12 prime numbers ranging from 53-401. The model development was performed using the following parameters: atom (A), bond (B), connection (C), chirality (Ch), hydrogen (H) and donor/acceptor (DA). The validity of the model depends on statistical parameters such as r 2 , q 2 by LOO, predictive r 2 pred and standard error. The robustness of the model depends on the more challenging test set prediction reflected by its predictive r 2 pred value.

Topomer CoMFA Model Analysis
The Topomer CoMFA model with good predictive ability in terms of r 2 (goodness of fit of the model) and q 2 (internal predictivity of the model) was presently developed. The model displayed a q 2 = 0.536 and r 2 = 0.975, with 0.460 standard error of prediction and 0.110 standard error of estimate. The number of components that provided the highest q 2 was eight. The summary of PLS results is provided in Table 2. The predictive ability of the developed topomer CoMFA model was assessed by the test set (eight molecules) predictions, which were excluded during the Topomer CoMFA model generation. The predictive ability of the test set was 0.604. The actual and predicted activities of the training set and test set molecules, along with the R1 and R2 fragment contributions are given in Table 3. The graph of predicted versus actual activity for training set and test set molecule is shown in Figure 2.

Contour Map Analysis
Topologically aligned R1 and R2 fragments are shown in Figure 3. Topomer CoMFA steric and electrostatic  contour maps for the R1 and R2 fragments of the most active molecule (28) are shown in Figure 4. Contour level along with color scheme and estimated volume of contour are summarized in Table 4. In the steric contour map, the green color denotes sterically bulky groups favoured for activity and the yellow color indicates sterically bulky groups unfavoured for activity. In the electrostatics contour map, red indicates electronegative-favoured groups and blue indicates electropositive-favoured group. The steric contour map for the R1 fragment ( Figure 4A) indicated that the bulky 2-nitro group on the phenoxy ring was favourable for activity. The electrostatics contour map ( Figure 4B) indicated that the 2-position was favourable for polar electronegative substituent and, hence, compound 28 displayed the highest activity (pIC 50 = 6.66) among all the compounds. In compound 29 (pIC 50 = 5.47), the 2-nitro group of compound 28 was replaced by the 2-amino group, which produced less of an inhibitory effect. The electropositive nature of the amino group was not favourable for activity when located at the 2-position on phenyl ring, and showed decreased inhibitory potency. For compound 27 (pIC 50 = 5.51), the unsubstituted phenyl ring was not favoured for activity; the compound displayed less inhibitory potency than compound 28. For the R2 fragment, the dimethoxy substitution was important for the inhibitory effect ( Figures 4C and 4D). The green contour near the methyl group indicated that the bulky group was favoured for the inhibitory effect and the red color near the oxygen atom indicated that an electronegative substitution could retain molecular activity. The data supports the idea that the methyl group allows a hydrophobic interaction with a receptor, with an oxygen atom acting as a hydrogen-bond acceptor. This indicates the hydrophilic interaction with receptor, which echoes earlier observations [14]. Steric and electrostatic contour maps for compound 30a (pIC 50 = 6.48) are shown in Figure 5. Fragment R1 consisted of a urea substituted derivative. The steric contour ( Figure 5A) revealed that a phenyl ring with a 2-nitro substitution is favoured for activity. The electrostatic contour map ( Figure 5B) indicated that the electron withdrawing nature of the nitro group was favourable for increasing inhibitory potency. For compound 30b (pIC 50 = 6.11), removal of the dimethoxy group from the R2 fragment produced a decreased activity compared to compound 30a, indicating the necessity of the dimethoxy substitution for the inhibitory effect. For compound 31 (pIC 50 = 6.11), and 32 (pIC 50 = 5.82), the position of nitro group played a major role in the inhibitory effect, as a para-substituted phenyl ring (compound 32) displayed less inhibitory potency than a meta-substituted phenyl ring (compound 31). The yellow contour corresponding to the para position of the phenyl ring indicated that a bulky group at this position was unfavourable for activity, perhaps due to the steric restriction of the receptor pocket. For compound 33a (pIC 50 = 4.68), 34 (pIC 50 = 4.85) and 35 (pIC 50 = 4.51), amino group substitutions at the ortho, meta and para positions led to the least-active compounds, indicating that, for P-gp inhibitory effect, a polar electronegative group (nitro) was more favourable than an electropositive group (amino).
Steric and electrostatic contour maps for compound 13 (pIC 50 = 6.33) are shown in Figure 6. Fragment R1 consisted of an amide-substituted derivative. The steric contour map ( Figure 6A) revealed that the bulky nature of the quinoxalinyl ring was favourable for inhibitory potency. The electrostatics contour map ( Figure 6B) revealed that the electropositive nature of quinoxalinyl ring was favoured for increased inhibitory potency. Replacement of the quinoxalinyl ring by a 2-naphthyl   Steric and electrostatic contour maps for compound 5b (pIC 50 = 4.85) are presented in Figure 7. The steric contour map for fragment R1 (Figure 7A), which displayed a yellow contour near the 2-nitro position, was  indicative of the lower favourability of that particular position for the inhibitory effect. The electrostatic contour map for the R1 fragment ( Figure 7B), which displayed a red contour near the nitro groups, indicated that the position was favorable for an electronegative group. For compound 17 (pIC 50 = 5.48), compound 18 (pIC 50 = 5.78) and compound 19 (pIC 50 = 5.62), substitutions with the less bulky and electronegative bromine at ortho, meta and para positions improved the potency over compound 5b. Fragment R2 (Figures 7C and 7D) displayed green and red contours near positions 6 and 7 of the 3,4-dihydroisoquinoline ring, indicating the importance of the dimethoxy group for the inhibitory effect. In contrast, the presence of a dimethoxy group in compound 5a (pIC 50 = 5.27) improved potency over compound 5b.

HQSAR analysis
The HQSAR model with good predictive ability in terms of r 2 and q 2 was presently developed. The model shows q 2 = 0.777 and r 2 = 0.956 with 0.302 standard error of prediction and 0.105 standard error of estimate. The model was developed with bond (B), connection (C) and donor acceptor (DA) parameters with BHL = 97. The number of components that provided the highest q 2 was six. Table 5 summarizes the PLS results. The predictive ability of the developed HQSAR model was assessed by the test set (eight molecules) predictions, which were excluded during HQSAR model generation. The predictive ability of the test set was 0.730. The actual and predicted activities of the training set and test set molecules are given in Table 6. The graph of predicted versus actual activities for the training set and test set molecules is shown in Figure 8. Numerous models with combinations of the A/B/C/H/ Ch/DA parameters with 4-7 atom counts are given in Table 7. Parameters such as B/C/DA were important for model generation. Bond and connection considered the bond order and hybridization states within fragments, respectively, and DA yielded information about the donor and acceptor atoms. The highest q 2 value was obtained for parameters B/C/DA; for further improvement of q 2 , an optional atom count (1-10) was explored.
A significant difference was noticed in the statistical parameters with different atom counts for model B/C/ DA (Table 8).
A standard color coding system was used to indicate atomic contributions in the HQSAR model. Red, redorange and orange designated unfavourable and negative contribution to the activity, while yellow, green-blue and green denoted favourable or positive contribution to the activity. White indicated an intermediate contribution to activity. For study of atomic contribution, molecules were selected randomly. The positive and negative atomic contributions for the selected molecules are shown in Additional file 1. All the molecules in the dataset had a common substructure and varied only in R 1 and R 2 substructure. The contribution map for        linker. This fragment was present in all molecules of the data set, indicating the importance of tertiary nitrogen for the inhibitory effect. This fragment may be preferred when designing a new scaffold for P-gp antagonism. Fragment F2, which contained an ethyl phenyl ring with a coefficient value of 0.001, was also present in all molecules of the dataset. Fragment F3, which had the same coefficient value, was present in some molecules. Fragment F3 contained a methoxy group on the aromatic ring; the fragment plays an important role in a favourable hydrophobic interaction, consistent with previous results [14]. The oxygen residue of the methoxy group acts as a hydrogen bond acceptor, which supports hydrophilic interaction with receptor, also consistent with previous observations [14]. In this particular series of compounds, the tetrahydroisoquinoline moiety was either unsubstituted or was substituted by a 6, 7-dimethoxy group. The difference in inhibitory effect of the dimethoxy substituted and unsubstituted compounds suggest that this substructure might be important in the inhibitory effect. Fragment F4 displayed a positive contribution of 0.004 for activity and was present in all molecules as a central phenyl ring. It might act through hydrophobic interactions with receptors, as has been previously suggested [13]. Fragment F6 was present in molecule 11 (pIC 50 = 5.33) and contributed  negatively for activity, with a coefficient of -0.002. This result indicates that the 4-quinoline substituent contributes negatively to activity and decreases inhibitory potency. Fragment F7 consisted of a 3-quinoline nucleus and was present in compound 9 (pIC 50 = 6.24) and 9b (pIC 50 = 6.37); the fragment contributed positively to activity (coefficient of 0.008). The aromatic (3-quinolinyl) 'N' acts as a hydrogen-bond acceptor with the appropriate amino acid of the receptor [14]. Fragment F8 was present in molecules 27 (pIC 50 = 5.51), 28 (pIC 50 = 6.66) and 29 (pIC 50 = 5.47), where it contributed positively to the inhibitory effect (coefficient of 0.006); these results indicate that an ether linker is more important for an inhibitory effect. Fragment F9 was present in the urea derivative compound; its' contribution coefficient of 0.004 was indicative of an inhibitory effect, highlighting the importance of urea substituted derivatives for the inhibitory effect. The model generated by both Topomer CoMFA and HQSAR agreed well with each other. The HQSAR analysis showed that tertiary nitrogen with an ethyl phenyl linker was essential for activity and that a dimethoxy group was necessary in inhibition of P-gp. Urea and ether linker were most important for the inhibitory effect and contributed profoundly. The fragment from the isoquinoline ring was also vital for activity and for the inhibitory effect.
Klinkhammer paper is about structure activity relationship (SAR) and this manuscript describes quantitatively structure activity relationship (QSAR), which is consistent with the previous paper. The present approach is in-depth study, i.e., contour map analyses of CoMFA and fragment analyses of HQSAR provided guidelines concerning compound modification.
The presently-developed QSAR models yielded similar findings as XR9576 reports [26]. The third-generation MDR modulators have the same scaffold as XR9576, with different side chains. The sterically bulky side chain of XR9576 could be responsible for the higher inhibitory effect [26]. This view is supported by the present findings of favourable steric (green) and electrostatic (blue) contour maps around the R1 fragment. It also highlights the importance of fused benzo rings. In the current model, the fused ring containing compounds (such as quinolinyl, quinoxalinyl and naphthyl) are predicted to show higher activity. The XR9576 report also highlighted the importance of the positional effect of nitrogen in quinolinyl ring for inhibitory potency [26]. This positional effect was presently confirmed quantitatively.

Designing New Compounds
Ligand-based methods such as Topomer CoMFA are not computationally intensive and can lead to the rapid generation of QSARs, from which the biological activity of newly designed compounds can be predicted. In contrast, an accurate prediction of activity of untested compounds based on the computation of binding free energies is both complicated and lengthy. The Topomer CoMFA contour maps provide clear indicators for designing novel molecules with improved P-gp inhibitory potency. The careful analyses of contour maps and HQSAR results led to the identification of the structural requirements responsible for compounds having improved potency. The information obtained from the contour maps of the most potent molecule (28; Figure 4) was utilized to design new R1 fragment containing compounds.
It was assumed that the 2-nitro group on the phenyl ring is favourable and responsible for retaining higher potency of inhibitors. The electrostatically favourable blue color at the para position of the 2-nitrophenoxy ring indicates that substituents with an electron-rich group at this position might increase the activity. Presently, a chlorine group was substituted at the para position and the phenoxy-acetamide group was replaced by N'-keto-benzohydrazide group. The resulting compound (28B) displayed improved inhibitory potency over compound 28. Here, using information obtained from the Topomer CoMFA model, five new R1 fragment containing compounds were designed, which displayed significant increases in activity. These newly designed compounds were included in the test set and their activity was predicted by the Topomer CoMFA and HQSAR models ( Table 9). The Topomer CoMFA contour map was predicted for the designed compound 28B shown in Figure 10. The (E)-N-3-(2-nitrophenyl)acrylamide (R1 fragment of compound 24) was replaced by a (2Z,4E)-5-(2-nitrophenyl)penta-2,4-dienenitrile (28A fragment) substituent, which resulted in equipotent compound as 28. The 2-nitrophenyl group remained unchanged because of its delicacy and importance at particular position for inhibitory effect, which corresponded to the sterically and electrostatically favourable contour maps. Another highly potent compound 28 derivative (28C) was designed by replacing the 2-nitrophenoxy-acetamide group of compound 28 by 2-(5-(4-chlorocyclohexyl)-1methyl-1H-imidazol-2ylthio)-N-acetamide. Here, the essential 2-nitrophenyl group was replaced by methyl substituted imidazole group by considering its positional and steric effects. The sterically favourable bulky nature of the nitro group was replaced by a methyl group, which corresponded to the favourable green contour map. The fifth position of the imidazole ring was substituted with a 4-chlorocyclohexyl group. The electronrich chloro group at this position was favourable for improved inhibitory potency, which corresponded to the blue contour map in the vicinity.
To check the importance of ortho substituents on the inhibitory effect, the nitro group was removed from the ortho position of compound 28. Decreased inhibitory potency was observed in compound 28D, indicating that retention of inhibitory potency of an inhibitor requires substitution of the phenyl ring with an ortho nitro group. Steric and electrostatic contour maps for the R2 fragment of 28B are displayed in Figure 10, and demonstrate the importance of a sterically bulky methoxy group for the inhibitory effect. We analyzed individual atomic contribution map of newly designed molecules for inhibitory effect. It indicates that the R1 fragment of designed molecules contribute positively towards the inhibitory effect, as denoted by the blue, green-blue and yellow color in HQSAR.
In summary, utilizing information obtained from the Topomer CoMFA and HQSAR analyses, we designed novel fragment containing compounds, which had higher inhibitory potency than the reported compound. From the overall analyses, we conclude that the 2-nitrophenyl group, which has a steric as well as polar nature, is responsible for the higher affinity of the molecules. Additionally, designed fragments underscore the importance of electron-rich substituents at the para position of the phenyl and cyclohexyl ring system.

Conclusion
We derived Topomer CoMFA and HQSAR models with good statistical values. The robustness of these models was confirmed using a test set. Topomer CoMFA analysis provided great insight into the structural requirements for improved potency over existing compounds. The information obtained from HQSAR model shows the importance of bond, connection and donor/acceptor parameters. The overall study indicates that, in HQSAR analysis, fragments containing information about the dimethoxy group are important for an inhibitory effect; this was supported by the findings of the Topomer CoMFA contour map. The contour map for ether and urea linking fragments (R1) indicate that substitution of bulkier and polar group to the ortho position of benzene ring enhances the inhibitory effect, and explains why the compounds with nitro group have good inhibitory potency. Contour map analysis also revealed that bulky and more electropositive substituents on the amide linker are responsible for higher potency; this was supported by the HQSAR atomic contribution map. A central phenyl ring could hydrophobically interact with a receptor and so is important in the inhibitory effect. In summary, both HQSAR and Topomer CoMFA underscore the importance of the aromatic dimethoxy and nitro groups for the inhibitory effect. Molecular modeling techniques like Topomer CoMFA and HQSAR aid the identification of the functional groups and atoms important for the inhibitory potency. Together, these data can be utilized to design more potent compounds than the present series of compounds.

Additional material
Additional file 1: Contribution map. Positive and negative contribution map for few molecules obtained by HQSAR analysis.