Abstract

An alignment-free, three dimensional quantitative structure-activity relationship (3D-QSAR) analysis has been performed on a series of -carboline derivatives as potent antitumor agents toward HepG2 human tumor cell lines. A highly descriptive and predictive 3D-QSAR model was obtained through the calculation of alignment-independent descriptors (GRIND descriptors) using ALMOND software. For a training set of 30 compounds, PLS analyses result in a three-component model which displays a squared correlation coefficient () of 0.957 and a standard deviation of the error of calculation (SDEC) of 0.116. Validation of this model was performed using leave-one-out, of 0.85, and leave-multiple-out. This model gives a remarkably high (0.66) for a test set of 10 compounds. Docking studies were performed to investigate the mode of interaction between -carboline derivatives and the active site of the most probable anticancer receptor, polo-like kinase protein.

1. Introduction

The -carboline alkaloids, which are originally isolated from the medicinal plant Peganum harmala [13], comprise a planar tricyclic system with different degrees of aromaticity and various substituents at positions 1, 2, 3, 7, and 9; the presence, location, and nature of these substituents play a crucial role in biological and pharmaceutical properties of these compounds [46]. According to previous investigations, -carbolines exert a broad spectrum of pharmacological effects such as anxiolytic, sedative [79], antimicrobial [10, 11], antiviral [12], antithrombotic [13, 14], and antiparasitic effects [15, 16]; also they are associated with alcohol dependence [17], and neurological diseases such as Parkinson’s disease [18, 19]. Moreover, a large series of -carbolines have been reported for their affinity with several receptors such as imidazoline [20, 21], benzodiazepine (BZ) [2224], 5-hydroxytryptamine (5-HT) [25, 26], and dopamine (DA) [27].

Recently, -carboline alkaloids have drawn attention due to their antitumor activity [5, 2833] and their functions through multiple mechanisms such as their ability to intercalate into DNA helix, inhibiting topoisomerases I and II [34, 35], monoamine oxidase [3639], CDK (cyclin-dependent kinases) [40, 41], and MK-2 [42]. -Carbolines have the potential to be used as anticancer drug leads. One of the major antitumor mechanisms is apoptosis in cultured HepG2 cells induced by -carbolines with downregulate the expression of Bcl-2 gene and upregulate the expression of death receptor Fas,without altering the level of Bax and P53 [43].

Previous data suggest that -carbolines can inhibit the activity of polo-like kinases (PLKs). The investigations of Han et al. show that PLK1 is attractive as a therapeutic target for inhibition of cancer cell growth [44, 45]. Among the four members of PLKs which are a family of serine-threonine kinases, PLK1 plays the most crucial role in cell proliferation and is an important regulatory enzyme in the G2/M transition through phosphorylation of substrates [46]. PLK1 contains an N-terminal catalytic domain and a C-terminal regulatory polo box domain (PBD) which is composed of two polo boxes. Since the PBD is contributing to PLK1 localization and function, the PBDs of PLK1 are promising alternative targets for designing novel antitumor drugs [47, 48].

The docking study was done using CDOCKER algorithm and the picture of the virtual receptor site was validated by investigation interactions between receptor and some inhibitors. CDOCKER is an implementation of a CHARMm-based docking tool using a rigid receptor that generates several prime random ligand orientations within the receptor active site followed by MD-based simulated annealing and final refinement by minimization [49].

This paper intended to get a better pharmacophoric pattern of these compounds and thus we obtained an alignment-independent 3D-QSAR model for the potency of 40 from total 47 -carboline derivatives which have significant pIC50 toward HepG2 human tumor cell lines (Table 1).

2. Materials and Methods

2.1. GRIND Descriptors

A new class of molecular descriptors called GRid-INdependent descriptors has been developed by Pastor et al. [50]. There are many advantages in creating such a predictive and robust 3D-QSAR model with a mathematical description of molecules over the other methods. First, in contrast to other methodologies in which one of the most difficult and time-consuming steps is molecular superimposition, in GRIND there is no need for alignment of compounds. Therefore the descriptors obtained are insensitive to coordinate frame of the space. Moreover, the GRIND descriptors are easily interpretive by referring back to the compounds so the original information can be reconstructed. Finally because of smaller variables/objects ratio in this method in comparison with other 3D-QSAR methods, the models obtained by improper use of variable selection with FFD (fractional factorial design) variable selection procedure are less prone to overfitting.

With the aim of getting a virtual molecular interaction fields (MIFs) pattern of the regions (fingerprint of receptor) so-called virtual receptor sites (VRS) to reveal the most common structural groups in the active compounds, first we should compute the MIFs of nodes. Therefore for the derivation of MIFs we used the four most recommended probes. To represent steric and hydrophobic interaction, hydrogen bond acceptor groups, and hydrogen bond donor groups, we use DRY (hydrophobic probe), O (carbonyl oxygen), and N1 (amide nitrogen), respectively. These probes represent strong noncovalent interactions between ligands and receptor. In addition, to consider molecular shape effects in the receptor-ligand interaction process, and as complementary to point-based interaction information, we used a supplementary probe, called TIP (shape probe), that extracts each ligand’s isosurface at 1 kcal/mol from the field of a normal GRID calculation. This method has been described elsewhere in more details [50, 51], but it is worth mentioning the whole procedure briefly.

In the absence of binding site knowledge it is possible to explore the process of ligand-receptor interactions with the help of the MIFs [52, 53]. The procedure for obtaining the GRIND starts with calculation of MIFs of the ligands using the program GRID, in order to obtain VRS described above. Due to the vast number of nodes, between 10000 to 100000, from which we obtain the MIFs for a drug molecule, we should filter these MIFs to obtain the most relevant regions and remove redundant information. This filtering uses an optimization algorithm, called Federov-like optimization algorithm [54], that selects from each field a fixed number of nodes (as a default value of 100) with optimizing a scoring function. MIF points that have low energy value and are as far as possible are retained. Also the scoring function can be tuned by giving to each criterion a different relative weight. In the next step the geometrical relationship between the VRS regions should be encoded into GRIND variables; thus there is no dependence upon the position and orientation of the molecular structure in the 3D space. The encoding is based on auto- and cross-correlation transform method called maximum auto- and cross-correlation (MACC) or, more specifically, MACC2. The continuous distances between pairs of filtered MIF points are divided into discrete “distance bins.” For each distance bin the product of the interaction is calculated but only the maximum energy product is conserved while the others are discarded. The output of MACC2 analysis can be represented directly in correlograms, where each point in it represents the product of two particular nodes within the distance bin separating the nodes of a certain compound. Wide peaks represent intense interactions which are produced by contiguous nodes around the compound. Conversely, weaker interactions are represented by narrow peaks. Each block in auto- and cross-correlogram indicates interactions between couples of nodes (node pairs) generated by the same or different probes, respectively. This approach does not require any scaling or pretreatment. The default settings were used for all other parameters.

2.2. Dataset

The dataset was adopted from the work of Cao et al. [55] in which they reported the cytotoxic potential of a number of synthesized -carbolines against a panel of human tumor cell lines (Table 1). The 3D structures of dataset were constructed in SYBYL7.3 molecular modeling package (Tripos Inc., St. Louis, USA), and the energy minimizations were performed using Tripos force field with a distance-dependent dielectric and the Powell conjugate gradient algorithm convergence criterion of 0.01 kcal/mol Å. Gasteiger-Hückel method is used to calculate partial atomic charge of all compounds. All other computational procedures were performed within the software ALMOND (version 3.3.0) running on Red Hat Enterprise Linux 4.7 workstation.

The dataset was partitioned into training and test set using most descriptive compound (MDC) method in which the compounds are weighted according to their population density [56]. The training set is composed of 30 compounds which were used to adjust the parameters of the models. The corresponding pIC50 values which are listed in Table 1 range from 3.63 for the most weak compound 45 to a value of 5.8 for the most potent compound 57 and cover a spectrum of approximately 3-log units. The 3D-QSAR model derived was successfully validated by using a test set of 10 similar compounds, with comparable pIC50.

2.3. ALMOND Model

The -matrix was generated, with compounds as its rows and the GRIND variables grouping into correlogram’s blocks as its columns. The information contained in this matrix can be directly used for chemometric analysis such as principal component analysis (PCA) and partial least square (PLS) analysis. A principal component analysis (PCA) was carried out to make a closer inspection on the distributions of the variables and compounds. The optimum number of PLS components (latent variables, LV) was chosen by monitoring changes in the model’s predicting index ( and leave-one-out) evaluated by applying the cross-validation procedure available in ALMOND. The optimum number was given at three PLS components. In order to select the most informative variables in the -matrix, FFD is applied that yields 145 -variables out of total 437 original -variables. The same data matrix was subjected to PLS analysis. Statistical parameters of the model obtained are summarized in Table 2. PLS analysis resulted in a 3-component model; a correlation coefficient of 0.96 and a standard deviation of the error of calculation (SDEC) of 0.12 were found. Also, the validation of this model using cross-validation method and an external test set shows = 0.85 of the validated model and the for the external test set was 0.66. The plot of experimental versus calculated values of pIC50 in Figure 1 proves the good quality of the PLS model obtained.

2.4. Molecular Docking

Molecular docking studies were carried out by using CDOCKER (CHARMm-based DOCKER) [57]. The crystal structure of PLK1 was retrieved from the RCSB Protein Data Bank (PDB entry code: 2OWB). For the preparation step of ligands all structures sketched and modified in SYBYL 7.3 molecular modeling package were transferred into Discovery Studio 2.5 (Accelrys Inc, San Diego, CA, USA) and typed with CHARMm force field and Momany-Rone partial charges calculation method [58]. Energy minimization was performed by SMART minimization algorithm, followed by conjugate gradient minimization [49]. The active site of the target protein is created as a spherical region with a radius of 9.18 Å. The next step is protein preparation in which all ligands and water molecules in PLK1(2OWB) were deleted and hydrogen atoms were added to the original crystal structure.

The CDOCKER score (-CDOCKER ENERGY) as a negative value includes receptor-ligand interaction energy and internal ligand strain energy [58]. Default values were assigned to the other parameters. High value of 25.45 for compound 57 as its dock score indicates favorable binding between this active inhibitor and kinase domain of PLK1.

3. Results and Discussion

3.1. The Results of ALMOND Model

PLS coefficients for the model are calculated using the DRY, O, N1, and TIP probes [59]. PLS coefficient plot representing the contribution of each single descriptor to the model with respect to the value of is shown in Figure 2. Since the regions of the MIF corresponding to each important variable can be visualized with the 3D plots included in ALMOND, the interpretation of these variables is straightforward. From PLS coefficient plot it can be seen that the most important variables which have positive impact on inhibitory activity are DRY-DRY: 24 and 25, TIP-TIP: 52 and 56, N1-TIP: 23 and 24, and N1-N1: 28. On the contrary, the analysis of all the distances at higher PLS coefficient showed that the variables TIP-TIP: 42, 43, and 44, N1-TIP: 10 and 7, and N1-N1: 8 correlate negatively with the activity. According to this plot the probe O has a marginal contribution.

3.2. ALMOND and Docking Interpretation

Variables 11–24 and 11–25 which have the strongest positive impact on antitumor potency are within the block of DRY-DRY node pairs. The variable 24 represents a distance of 9.6 Å between hydrophobic nodes due to ring A and phenyl ring in 2-benzyl substituent in the two compounds 57 and 58 (Figure 3(b)) (the bold letter faces show the number of compounds). The variable 25 is the same for compound 57 but in a distance of 10 Å in compound 58 this variable is due to hydrophobic properties of ring B and 2-benzyl in this series of -carboline derivatives.

According to Figure 2 all DRY-DRY interactions are positively related to antitumor potency and differ in their relative intensity. These variables referring to favorable regions for hydrophobic groups mainly arise from interactions of triple ring system of all -carbolines and the aromatic moiety far from template. However the long aliphatic saturated chains are hydrophobic regions in many compounds (e.g., the compounds 33, 47, 48, 50, and 51). These variables can be seen clearly in the active compounds such as the compounds 56, 57, 58, and 60. Those variables are also present in compound 45, due to its pentafluorobenzyl moiety, although it is one of the most inactive compounds.

The results of molecular docking -carboline derivatives of PLK1 are overlaid in Figure 4; it is clearly observed that compound 57, which is one of the most active compounds, interacts with the active site of PLK1 through an extensive hydrophobic region, including mainly the following residues: PHE58, LEU59, GLY60, ALA80, VAL114, LEU130, LEU132, PHE183, and GLY193; thus hydrophobic groups would benefit the potency.

Compound 57 also exhibits a strong binding affinity through stacking with the PHE183 residue of PLK1. The stacking interaction between PHE183 of PLK1 and aromatic rings either in the backbone or one of substitutions of compounds like pentafluorobenzoxyl moiety has crucial role in enhancing the binding affinity since it is absent from most of the compounds with lower pIC50.

Within the block of TIP-TIP node pairs, the largest positive impact on antitumor potency is attributed to variables 52 and 56 (node-pair distances of 20.8 and 22.4 Å, resp.). The variables 44–52 and 44–56 are due to the terminal benzyl ring and the aromatic or ether group on the other end of the molecule in compounds 56, 57, 58, and 60. Compound 50 is active but since it has no benzyl group this variable is zero according to the correlogram. Actually these variables are the most discriminating variables that can differentiate active from inactive compounds (Figure 5). On the other hand, GRIND variables 42, 43, and 44 that correspond to TIP-TIP node-pair distances of 16.8, 17.2, and 17.6 Å, respectively, exhibit strong inverse correlation to the biological activity (). So these variables are present in compounds 3 and 45 (Figure 3(a)).

Within the N1-N1 block in PLS coefficient plot, it is worth mentioning that the negatively correlated bars are located on the left side (representing smaller node-node distances) and the positively correlated variables are positioned on larger node-node distances, that is, on the right side. It is due to the fact that structural elements which exert negative impact on potency are located closer to each other in molecules than those having positive impact. According to the N1-N1 autocorrelogram, regions favorable for hydrogen bond acceptor groups separated by 3.2 Å (i.e., variable 33–8) seem to be characteristic of inactive compounds. They are associated to compounds 4, 23, 3, 17, 11, 55, and 45 in respect to their energy products value. These variables correspond to the interactions of the polar carbonyl and carboxyl groups as hydrogen bond acceptors separated by 3.2 Å, being of similar energy product for the two compounds 3 and 23.

Figure 6 shows that the carboxyl group at position 3 of compound 3 binds via a hydrogen bond to amino group of ARG136 which indicates that the presence of H-bond donor groups in this region would emphasize the lower affinity of this compound to the active site of enzyme.

Also there is a similarity in node-node energy products between compounds 11 and 17 due to primary amine nitrogen, and also between compounds 45 and 55 due to pentafluoro benzoxyl ring. The highest value of energy interaction product for the variable 33–8 in compound 4 which is clearly evident in correlogram means that the N1 probe interacts more strongly with nitrogen atom in –C(O)NHCH2CH2OH than the groups described above.

According to the docking results there are two hydrogen bonds in the –C(O)NHCH2CH2OH moiety. One is between oxygen of CONH and alcoholic hydrogen within the molecule itself, and the other is between alcoholic oxygen at the end of the moiety and LYS61 amino acid.

Overall, the hydrogen bond acceptor group in the side chain of ring A at position 7 forms a hydrogen bond with LYS82 or ASP194, such as oxygen of ether group in compounds 56, 32, and 43 or fluorine of pentafluorobenzoxyl in compound 57, which plays an important role in the PLK1 inhibitory activity of -carboline. By contrast, presence of hydrogen donor in this position forms unfavorable hydrogen bond with LYS61 like in compounds 40, 42, and 43. The next hydrogen bond formation represents negative influence on inhibitory potency observed for compounds 49 and 53 due to pyridinium nitrogen (ring C) so this nitrogen atom in its protonated form would be more desirable for potent inhibitory activity.

The interpretation of all relevant peaks is summarized in Table 3.

At a deeper insight, in order to determine the most relevant descriptors which can properly distinguish between active and inactive compounds, 6 compounds of training set (1, 3, 45, 57, 58, and 60) and 2 compounds of test set (21 and 59) were selected; then existence of all variables which have high positive or negative PLS coefficients described above is carefully inspected for these 8 compounds, including 4 active and 4 inactive -carboline derivatives (Figure 5). Variables with positive PLS coefficient should be present in active compounds and should be absent from inactive compounds and vice versa. Green boxes in Figure 5 represent presence or absence of variables as expected, and yellow boxes represent unexpected results. As it is shown, variables 44–52, 33–8, and 34–10 are the most discriminative variables. Since more than 65% of boxes are highlighted in green it can be concluded that the obtained model is very predictive and efficient. Also the yellow cells are a consequence of attempting a univariate interpretation of a multivariate model.

3.3. Summary of Structural Insights from 3D-QSAR Models and Docking Studies

The final results of analysis of 3D-QSAR models and docking studies are summarized in Figure 7. As discussed above, the presence, spatial shape and chemical characteristics of R2 and R7, and their mutual distance have high impact on -carbolines inhibitory activity. Pyridinium nitrogen of ring C can form hydrogen bond with ASP194 and LYS82; however there is negative association between presence of this moiety and the activity. Also there should be an electron donating group like benzyl to form N1-N1:28. This substituent should be bulky enough to form the TIP-TIP variable in a distance of 20.8–22.4 Å from R7. Otherwise it may represent the TIP-TIP variable in a distance of 16.8–17.6 Å, with negative impact. At position 3, carbonyl or carboxyl groups as hydrogen bond acceptor are disfavored, because they form the N1-N1: 8 variable with negative impact.

4. Conclusion

Within this paper we aimed at development and validation of a ligand-based 3D-QSAR model in order to get a deeper insight into molecular structure-antitumor potency relationship in a series of -carboline derivatives. The method used in this paper is based on alignment independent descriptors whose particular value is the ability of backtracking of all interactions to the original filtered GRID field.

The hydrophobicity (associated with the DRY probe), shape effects (associated with TIP probe), and hydrogen bond acceptor-donor interactions (associated with N1 probe) are the main factors that determine antitumor potency toward HepG2 cell lines, within studied set. In addition, molecular docking is carried out to map the binding pocket of the PLK1 and its characteristics. And the interactions described with docking studies are similar to those described through analysis of 3D-QSAR model. Results obtained in the present study can be used as guidance for design of novel drugs.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.