Combined Pharmacophore Modeling, Docking, and 3D-QSAR Studies of PLK1 Inhibitors

Polo-like kinase 1, an important enzyme with diverse biological actions in cell mitosis, is a promising target for developing novel anticancer drugs. A combined molecular docking, structure-based pharmacophore modeling and three-dimensional quantitative structure-activity relationship (3D-QSAR) study was performed on a set of 4,5-dihydro-1H-pyrazolo[4,3-h]quinazoline derivatives as PLK1 inhibitors. The common substructure, molecular docking and pharmacophore-based alignment were used to develop different 3D-QSAR models. The comparative molecular field analysis (CoMFA) and comparative molecule similarity indices analysis (CoMSIA) models gave statistically significant results. These models showed good q2 and r2 pred values and revealed a good response to test set validation. All of the structural insights obtained from the 3D-QSAR contour maps are consistent with the available crystal structure of PLK1. The contour maps obtained from the 3D-QSAR models in combination with the structure based pharmacophore model help to better interpret the structure-activity relationship. These satisfactory results may aid the design of novel PLK1 inhibitors. This is the first report on 3D-QSAR study of PLK1 inhibitors.

structure-based alignments were both found to be statistically valid in terms of the interpretation of interaction mode and the predictability to internal and external compounds. The contour plots obtained from the 3D-QSAR models correlated well, not only with the detailed interactions between the ligands and active-site residues in the crystal structures of PLK1, but also the pharmacophore features directly derived from the receptor-ligand interactions in crystal structures. The developed computational models are expected to help with better understanding of the QSAR of this class of compounds, as well as ensuring the researcher an in-depth analysis about the lead compounds for PLK1 inhibitor in further studies. To the best of our knowledge, this work will be the first 3D-QSAR study of PLK1 inhibitors reported.

Multiple Structure Alignment Analyses
The accuracy and reliability of the CoMFA and CoMSIA model is directly dependent on the structural alignment rule. Therefore, before PLS analyses to construct the 3D-QSAR models, we performed the structure and ligand based alignment to find the effective alignment to this dataset (Table 1). Because the alignments involved were actually based on the co-crystal structure of compound 73 with PLK1, a preliminary analysis on its binding mode was necessary. Figure 1 shows the co-crystal interaction mode of compound 73 with PLK1 (2YAC, resolution: 2.2 Å). The core of compound 73 is sandwiched between Phe183 at the bottom of the ATP binding pocket and Cys67 in the back of the G-loop. An aromatic ring stacking interaction is found between the Phe183 and compound 73, which has an important influence on the conformational equilibrium of the whole compound. The 4-methylpiperazinyl moiety penetrates to the solvent accessible region, which may be involved in hydrophilic interactions. The 2-hydroxyethyl group positions at the same place related to the ribose moiety of ATP, which is a site tolerant to the long chain substituent. In addition to two conservative hydrogen bonds formed with the hinge region residue Cys133, the amide moiety and the trifluoromethoxyl group are engaged in three and one hydrogen bonding interactions, respectively. The knowledge on binding mode will assist in the evaluation of compound alignments as well as QSAR analyses.  As depicted in Figure 2, all 73 compounds were aligned well, using the common substructure based method. GLIDE performed quite well, as most conformations bind in a way analogous to the bound ligand of 2YAC, i.e., compound 73 ( Figure 3). Thus, the alignment derived by GLIDE docking is considered reasonable. Figure 4a,b illustrate two pharmacophore models deduced from the PLK1 crystal structures 2YAC and 3KB7 by LigandScout. It is obvious that they share nearly identical features. To cover the most common features that may be required by PLK1 inhibitory potency, we clustered and subsequently merged them to a new pharmacophore model. This merged model ( Figure 4c) consists of one hydrogen-bond acceptor, one hydrogen-bond donor, three hydrophobic and one ionizable positive, which is simplified by discarding three redundant hydrophobic features. Figure 4d shows the result of pharmacophore mapping of those compounds, which also suggests an excellent alignment.    Features are color-coded with magenta for hydrogen-bond donor, green for  hydrogen-bond acceptor, light-blue for hydrophobic, red for ionizable positive. Taken together, the results from common substructure and the merged pharmacophore based alignments, as well as GLIDE docking, proved to be reasonable and effective. It was difficult to judge which alignment would be more practicable, therefore, they were all subjected to the next step for model generation to further investigate their applicability and gain a more extensive insight to QSAR of pyrazoloquinazoline PLK1 inhibitors.

CoMFA and CoMSIA Statistical Results
Owing to the fact that common substructure, GLIDE docking and the pharmacophore methods all produced acceptable alignments of 73 known PLK1 inhibitors, the corresponding CoMFA and CoMSIA analyses were performed independently for further comparison. The statistical results of PLS analyses for CoMFA and CoMSIA studies are listed in Table 2. The pharmacophore based model yielded q 2 = 0.628 and r 2 = 0.941 for CoMFA, whereas the GLIDE docking and common substructure based model produced a lower q 2 value of 0.283 and 0.578, and r 2 value of 0.420 and 0.867 for CoMFA, respectively. Multiple CoMSIA models were derived based on three types of alignment, with various combinations of steric, electrostatic, hydrophobic, hydrogen bond donor and hydrogen bond acceptor fields. To get a clear view, only parameters of models whose q 2 value are higher than those of other models derived from the same alignment were considered. It is obvious that CoMSIA models from common substructure, GLIDE docking and pharmacophore based alignments showed comparable results. The inhibitory activities (pIC 50 ), the predicted activities using the CoMFA and CoMSIA models, and the corresponding residual values for the training set compounds are listed in Tables 3 and 4. Graphic representations of experimental vs. predicted inhibitory activity of training set for pharmacophore-based CoMFA and typical CoMSIA models are shown in Figure 5. In all, the CoMFA and CoMSIA models we constructed possess high q 2 and r 2 value, indicating that they have good internal predictive ability and that results were not based on any chance correlation. To validate both the predictability and accuracy of the models for external compounds, the predictive correlation coefficient r 2 pred was calculated for the test set. As shown in Tables 3 and 4, the r 2 pred value of pharmacophore-based CoMFA model and CoMSIA models spans from 0.605 to 0.827 and most of the residual values are less than 1.0, revealing that the models are highly reliable and can be used to predict the biological activities of novel compounds; whereas, the r 2 pred value of CoMFA models from GLIDE docking and common substructure based alignment is somehow lower, reflecting poor predictive ability. The plots of experimental vs. predicted inhibitory activity of test set for CoMFA and CoMSIA models are shown in Figure 5, showing that the predicted activities were in good agreement with the original data and the reliable CoMFA and CoMSIA models have a robust external predictive ability.
It can be concluded easily that the best model for CoMFA was obtained from the pharmacophore-based method (model 2) while it was difficult to distinguish the best CoMSIA model because there is no significant difference between PLS statistical results. Hence, we have paid attention to all CoMSIA models (models 4, 5 and 6), considering the representation of different fields, the satisfactory internal and external predictive ability in terms of q 2 and r 2 pred value, respectively.

CoMFA Contour Maps
The results of CoMFA analyses from pharmacophore-based alignment (model 2) are displayed as color-coded contours, allowing visual inspection of regions responsible for favorable or unfavorable interactions with PLK1. The green contours indicate regions where bulky substitution enhances binding affinity, and the yellow contours suggest regions where bulky substitution reduces the binding affinity. In the electrostatic interaction map, the blue contours indicate regions where more positively charged substituents are favored and the red contours suggest regions where more negatively charged substituents are favored. The favorable and unfavorable contributions of both fields were plotted as default proportion (80:20). Since the QSAR models were developed based on the information from receptor (docking, structure-based pharmacophore and common substructure from crystal bound conformation), the contour maps produced by CoMFA and CoMSIA could be superimposed onto the PLK1 structure. Thus, to get a straightforward insight into the steric and electrostatic interaction between compounds and PLK1, we introduced the van der waals surface or electrostatic potential surface of protein as background.
The steric and electrostatic fields contribute to 61.9% and 38.1% of the variance, respectively. The steric contour map is shown in Figure 6a with one of the most potent inhibitors, i.e., compound 73 as a reference. A moderate green contour is seen in proximity to the o-trifluoromethoxyl group of phenyl ring, but sandwiched by the protein surface, suggesting that only the medium-sized substituent is favored at this position such as methoxyl, trifluoromethoxyl, and methyl. The large green contours are found around the 4-methylpiperazinyl moiety at 5′ position of the phenyl ring, which penetrate to the solvent accessible region. This indicates that diverse substituents with bulky size at 5′ position of the phenyl ring are favorable to activity and their orientations are tolerable in space, except those extending to the upper side. Those situations for steric favorable substituents are the same for  compounds 35, 42, 44, 47, 48, 51, 53, 58 and 59, all of which show better activity (below 10 nM) and have moderate and bulky groups at 3′ and 5′ positions of the phenyl ring. However, compounds that have only one steric favorable site show only moderate binding affinity, e.g., compounds 5, 12, 14, 15,   22, 25, 28, 29, 38 and 33, demonstrating that both sites are indispensable. The emergence of yellow contours in the front of 2-hydroxyethyl group suggests that a more bulky substituent at position 1 of pyrazole ring would lower the activity. For example, replacement of 2-hydroxyethyl group with longer substituents (compounds 69, 70 and 72) led to a significant reduction of potency. Another yellow contour over the phenyl ring shows that a bulky substituent in this area is not favorable; but when the phenyl ring is flipped vertically by ~90 °, this area can also be occupied by substituents, such as compounds 11, 16, 21, 27, 30 and so on. In all, there is a definite requirement of an appropriate shape to exhibit high potency when designing novel PLK1 inhibitors, and thus it is important to pay attention to the steric characteristics. Electrostatic contour maps are also shown with compound 73 as a reference (Figure 6b). In general, red contour maps are close to heteroatoms such as nitrogen and oxygen, whose partial atomic charges are highly negative. Two main red contours are found close to 4-methylpiperazinyl moiety, of which N atoms bear negative charges, indicating negative potential is preferred in these areas. This trend can be reflected in the activities of compounds 35, 44, 47, 48, 49, 51, 52, 53, 58, 59 and 73, which all have tertiary amine substructures. Confusingly, the simple replacements of 4-methylpiperazinyl moiety with other tertiary amines for compound 53, either open-chain or cyclic, result in decreased potency against PLK1 from 10-fold to more than 300- fold, for compounds 40, 41, 43, 56, 57, 60, 61, 62, 63, 64, 65 and 68. b a This can be explained by taking the steric factor into account, that an unfavorable steric contour exists over the 4-methylpiperazinyl moiety as illustrated in Figure 6a. With respect to the favorable positive potential, the blue contours are distributed around the o-trifluoromethoxyl group, 2-hydroxyethyl group and amide moiety, where the negative charge of protein surface is concentrated, suggesting that positive charged groups such as substituents with electron-withdrawing atoms increase the activity compared to the hydrogen atom. This is consistent with the increase of potency for compound 53 as compared to compound 35, due to the replacement of hydrogen atom with fluorine atom. A similar situation can be observed between compounds 71 and 73, for which the replacement of nitrogen atom with oxygen atom leads to increased activity. Besides, two small blue regions close to the amide moiety also represent the preferred electrostatic interaction around there, indicating that decreasing the electron-withdrawing effect would cause a reduction of binding affinity, such as compounds 1, 2, 3, 4 and 36. Collectively, the distribution of electrostatic contours is well corroborated with the electrostatic potential information of the PLK1 active site, demonstrating its rationality.

CoMSIA Contour Maps
The CoMSIA contour maps of four models (models 4, 5 and 6), based on different alignment methods or combinations of various fields, are shown in Figure 7, which are depicted with compound 73 as the default from nearly the same viewpoint for a convenient comparison and analysis. The steric and electrostatic contours are colored identically with CoMFA contours map. In addition, the hydrophobic interactions are shown by yellow and white contours, whereas hydrogen bond donor interactions are represented by cyan and purple contours, indicating their favorable and unfavorable regions. The contributions of fields to the variance are listed in Table 2. It is evidenced that the electrostatic field contributes about 1.7-3 times more than the steric field, which is the opposite of corresponding relationships in the CoMFA models. As the alignment results are identical for CoMFA and CoMSIA models, this discrepancy may be explained by the different implementations of the fields for CoMFA and CoMSIA. Thus, we assume that both steric and electrostatic fields play important roles in the binding affinity and should be given equal attention.
Since the steric and electrostatic interactions have been discussed above in detail, a critical eye has been given to the comparison of the contour distributions. As shown in Figure 7a-c, the green contours mainly concentrate around the 4-methylpiperazinyl moiety and the trifluoromethoxyl group, denoting that the bulky substituents are indeed favorable to these regions; a large yellow contour is constantly located between the 2-hydroxyethyl group and 4-methylpiperazinyl moiety, suggesting a potential steric clash may exist between those two substituents. However, there is a distinct difference in the steric contour above the 2-hydroxyethyl group as the green contour in model 4 is conversely yellow in model 6. It can be found that a small sub-pocket is positioned over the 2-hydroxyethyl group, which means bulky substituents are not acceptable, such as 2-(tetrahydro-2H-pyran-2-yloxy)-ethyl and 3-aminopropyl groups, leading to lower binding affinities for compounds 70 and 72, respectively. Therefore, we conclude that a small yellow contour is more appropriate in that position. As for the electrostatic field, it can be observed that the contours of CoMSIA models are more concentrated around the 2-hydroxyethyl group and 4-methylpiperazinyl moiety (Figures 7d-f) in comparison with those of CoMFA model. In spite of this, common characteristics still exist. Commensurate with the CoMFA model, in all three CoMSIA models, a blue contour and large red contour are close to the trifluoromethoxyl group and the 4-methylpiperazinyl moiety, respectively, denoting the electrostatic nature of those two positions are authentically reflected. In addition, a middle-sized blue contour is also proximate to 4-methylpiperazinyl moiety in the CoMSIA model derived from common substructure based method (model 6). Since the most frequently used substituents at position 5′ of phenyl ring are various tertiary amines, the blue contour may account for the electron-deficient methylenes. Therefore, the substituents at position 5′ of the phenyl ring should consist of the electron-withdrawing atom and the electron-deficient atom simultaneously rather than solely the electron-withdrawing atom, such as carbamoyl or sulfoamino group. This conclusion is supported by the fact that the position 5′ of phenyl ring is oriented to the solvent accessible region of PLK1 and the substituents with the ionizable groups at that position will be ionized in the solution, stabilizing interaction and enhancing potency.
Areas favored by hydrogen bond donors are shown in cyan and magenta, respectively (Figures 7g-i). For all three CoMSIA models, two cyan contours are equidistant and close to the amino group of the amide moiety, mirroring the potency of two hydrogen atoms in the NH 2 group to form hydrogen bonds with the residues of receptor in the corresponding orientations, such as Asp194. These contours can be associated with the increment in activity when the NH group of the amide moiety changes from ethoxyl and substituted amino groups in compounds 1, 2, 3, 4 and 36, implying the NH group plays a major role in binding to the PLK1 active sites. The purple contours are found in common around the 4-methylpiperazinyl moiety, denoting the disadvantage of the hydrogen bond donor at this position for activity. This is corroborated with the distribution of electrostatic contours. A small purple contour is shown near the trifluoromethoxyl group that formed a hydrogen bond with the backbone NH group of Arg136 in our models and the guanidine NH group of Arg57 in PLK1 crystal structure (2YAC), respectively. Thus, the hydrogen bond acceptors at this position are favorable. A confusing purple contour is observed near the NH group of the amide moiety. As a crystal structure, and our models confirm the favor of the hydrogen bond donor at this area, this contour cannot be associated to the NH group. From a systematic investigation of the conformations superimposed with contours, we found the substituents at position 1 of pyrazole of compounds 69 and 72 might account for that purple contour, whose hydrogen bond donor groups reach this point due to the flexibility of the alkyl chain.

Comparison of Pharmacophore Model and CoMSIA Model
Considering that the pharmacophore model we have constructed also consists of hydrogen bond related features and a hydrophobic feature, it would be significant to compare it with the CoMSIA models. Hence, the merged pharmacophore model was reproduced using Unity module in SYBYL 6.9 [25]. The graphical interpretation of the superimposition of the features and contours reflecting the hydrophobic and hydrogen bond donor fields is shown in Figure 8 with compound 73 as a reference. The contours reflect the corresponding fields of model 6 (derived from pharmacophore-based alignment). Two large green contours are close to two hydrophobic features located at the trifluoromethoxyl group and in the vicinity of the phenyl ring, suggesting bulky and hydrophobic substituents at these positions are favorable (Figure 8a). The ionizable positive feature is covered by the large red contour, indicating a hydrophilic characteristic is preferred here (Figure 8b). The hydrogen bond donor and acceptor features do not intersect with the contours related to the hydrogen bond donor field (Figure 8c). Despite only partial pharmacophore features overlapping well with corresponding contours, it is still considered reasonable because the features not overlapped with contours belong to the maximum common substructure of compounds or are conservative in most compounds used in this study, while the CoMSIA method mainly focuses on the variable parts for a class of compounds. In this sense, our 3D-QSAR model and pharmacophore model complement each other well in elaborating the interaction mode of compound.  Figure 4. The contours are depicted as mesh.

Structural Insights from 3D-QSAR and Pharmacophore Studies
Our analyses found that the electrostatic, steric and hydrogen bond donor characteristics are highly desirable for potent inhibitory activity. The contour maps show that a moderate bulky substituent with hydrogen bond donor at position 3 of pyrazole ring (R 1 ), a moderately bulky and hydrophobic group with electron-withdrawing heteroatom at position 2′ of phenyl ring (R 3 ), and a bulky and amphoteric substituent with both hydrophilic and hydrophobic moieties at position 5′ of phenyl ring (R 4 ), can play important roles in enhancing binding affinity. Moreover, the length of alkyl chain for the substituent at position 1 of pyrazole ring (R 2 ) should be no more than three carbon atoms. In addition, two hydrogen bonds formed between compounds and Cys133 in hinge region are important which can induce the conformation of the whole compound in the ATP pocket. Thus, influencing these two hydrogen bonds should be avoided when changing substituents at other positions. These insights are consistent with the structural features of ATP pocket, further indicating that our 3D-QSAR models are reasonable. As depicted in Figure 9, R 2 and R 4 groups project into the solvent accessible region, thus allowing a larger extent of variability for the steric, electrostatic and other properties of substituents. Nevertheless, the steric clashes between R 2 and R 4 groups should be avoided. Moreover, R 2 and R 4 groups still have an impact on PLK1 selectivity. The reasonable combinations of substituents at these two positions can increase the selectivity of PLK1 against PLK2-3 up to 5000 times. Although R 1 group can form hydrogen bonds with a water molecule, substituents that discard this interaction may also be positive to the enhancement of binding affinity and target selectivity. This is concluded from the study of Fernandez and coworkers, which indicated that sculpting the shifting hydration patterns of the target would stabilize the protein surface and avoid disfavored induced fit [26]. Figure 9. Schematic representation for the SAR of 4,5-dihydro-1H-pyrazolo [4,3-h]quinazoline derivatives as PLK1 inhibitors. R 1 : medium-sized substituent with hydrogen bond donor and acceptor; R 2 : open-chain alkyl group with less than three carbon atoms or unsubstituted hydrogen; R 3 : hydrophobic group with small size and strong electron-withdrawing atom, especially hydrogen bond acceptor; R 4 : bulky substituents simultaneously with hydrophobic and hydrophilic moiety.

Dataset
All compounds used in the present study were taken from the literature [18,20,21]. Of the 73 compounds, 52 ones (unasterisked in Table 1) were selected randomly as training set for model construction and the remaining 21 ones (asterisked in Table 1) were used as test set for model validation, according to biological activity range and structural diversity. The IC 50 values of all compounds for PLK1 inhibition were normalized and converted to the logarithm unit of molar grade (pIC 50 = −log IC 50 ), which spanned 4 orders of magnitude (5.00-8.70). The distribution of activity data and the number of compounds were shown in Figure 10 to confirm with the test set as a true representative of the training set. The X-ray crystal structures of this class of compounds bound with PDK1 are available from the protein data bank (PDB). The bound conformation of compound 73 (PDB code: 2YAC) [18] was used as a template to build the 3D structures for both training and test set compounds. The partial charge was calculated with Gasteiger-Hückel method. The common structure was constraint for each compound and only the varying parts were energy minimized by using conjugate gradient method and Tripos force field until an energy gradient of 0.05 kcal/mol was reached. These works were all done in SYBYL 6.9.

Conformational Alignment
Structure alignment is considered as an important and critical step in CoMFA and CoMSIA analyses because this affects the reliability of the models. In order to avoid bias towards a particular alignment method, the structure-based and ligand-based alignments were both used in this study. It should be noted that a study that specifically seeks to understand the influence of alignment methods on the predictive performance of 3D-QSAR model is an important direction but extended in the work presented here. Herein, the common substructure, molecular docking and pharmacophore-based alignment were performed to build the 3D-QSAR models. Meanwhile, the docking and pharmacophore studies would also provide beneficial insight into ligand-receptor interactions to help better understand the structure-activity relationship.

Common Substructure Based Alignment
The key assumptions of CoMFA and CoMSIA are that compounds with common substructure always adopt a similar conformation when binding with the target and the common substructure in each compound contributes equally. Therefore, we selected the co-crystal structure of compound 73 from 2YAC as the template to align the remaining compounds using the "align database" command in SYBYL 6.9. The common substructure used for the alignment is shown in Figure 11.

Molecular Docking Based Alignment
Molecular docking was carried out to obtain reasonable molecular alignments for developing receptor-based 3D-QSAR models. At the beginning, we tested the applicability of three well-known docking software, viz. CDOCKER [27,28] in Discovery Studio 2.5, GOLD 5.0 [29,30] and GLIDE 4.5 [31,32] in Maestro 8.0, by checking if the conformation of the bound ligand in PLK1 crystal structure can be reproduced, and whether the common substructure of all compounds in both training and test sets can overlap well with each other in a way analogous to the bound ligand in PLK1 crystal structure. Docking conformations output by both CDOCKER and GOLD overlapped in a chaotic state, suggesting a failure of alignment. In contrast, GLIDE performed quite well. Thus, GLIDE was eventually selected as the docking tool.
The 3D structure of PLK1 (2YAC) in docking study was downloaded from Protein Data Bank. For GLIDE, the PDB structure was prepared using the "protein prepare wizard" automatically and subsequently its grid file was generated in Maestro 8.0. The initial conformation of compound used was obtained by conformational search in water with force filed of OPLS_2005 based on mixed torsional/low-mode sampling method in Maestro 8.0. The binding site was defined by the co-crystal ligand (compound 73) itself for all three docking software. The XP mode (extra precision) was selected and post-docking minimization was conducted to penalize highly strained ligand geometries and eliminate poses with eclipsing interactions. Finally, other options not mentioned above were kept as default.

Pharmacophore Based Alignment
The structure based pharmacophore model can be derived directly from ligand-protein co-crystal structure and thus can reflect more reliable combination of the essential features required for the relating biological potency [33]. As the compounds used in 3D-QSAR analyses belong to the same class and the co-crystal structures of PLK1 are available, the structure-based pharmacophore was generated utilizing LigandScout 2.02 [34], which is based on a sophisticated ligand-protein complex interpretation algorithm. Two PLK1-ligand co-crystal structures (2YAC and 3KB7) [18,20] available were chosen. When creating pharmacophore model, the "Phase" mode was selected with waters and other heteroatom ignored due to their non-conservation in crystal circumstance. This produced two pharmacophore models. Considering that pharmacophore should contain the most common features and these two models indeed share some identical features, we compared and clustered them in Discovery Studio 2.5 to draw a new pharmacophore model. This model was eventually used to align compounds in Discovery Studio 2.5, during which the conformations of compounds were generated with "best" option and the fitting method was "flexible" with the maximum omitted features of 3.

CoMFA and CoMSIA Methodology
The CoMFA and CoMSIA analyses were carried out with the RHEL 4.0 operating system using SYBYL 6.9. In CoMFA study, the aligned compounds were placed in the 3D cubic lattice with grid spacing of 2.0 Å. The standard CoMFA steric and electrostatic fields were calculated using a sp3 carbon atom as steric probe and a +1 charge as electrostatic probe, with Lennard-Jones potential and the Coulombic potential, respectively. The cut off value for both fields was set to 30 kcal/mol and the minimum-sigma (column filtering) was set to be 2.0 kcal/mol to reduce the noise by omitting those lattice points. The five fields of CoMSIA (steric, electrostatic, hydrophobic, hydrogen bond donor and acceptor) were calculated for each lattice with a grid of 2 Å by employing a common probe atom with 1 Å radius, +1 charge, and hydrophobic and hydrogen-bond property values of +1 [24]. The attenuation factor was set to the default value of 0.3 for the Gaussian function.

Partial Least Squared (PLS) Analyses and Validation
The relationship between the structures and the biological activities derived by the PLS algorithm. CoMFA and CoMSIA descriptors were used as independent variables and pIC 50 values were used as dependent variables in PLS to generate corresponded 3D-QSAR model. The predictive ability of the models was evaluated by leave-one-out (LOO) algorithm, which gave the optimal number of component (ONC), the lowest standard error of prediction and cross-validation coefficient (q 2 ), calculated with Equation 1, where Y pred , Y exp and Y mean are the values (pIC 50 ) of the predicted activity, experimental activity and mean activity for compounds in training set, respectively. The analysis of non-cross validation was performed to calculate the conventional r 2 using the ONC obtained from the LOO analysis. Validation of the utility of the model as a predictive tool was carried out by predicting the activity of an external test set of 21 compounds. The predictive r 2 (r 2 pred ), reflected the predictive power of the CoMFA and CoMSIA models, was calculated using Equation 2, where SD is the sum of the squared deviations between the experimental activities of the compounds in the test set and the mean activity of the compounds in the training set, PRESS is the sum of the squared deviations between predicted and experimental activities for every compound in the test set.

Conclusions
The 4,5-dihydro-1H-pyrazolo [4,3-h]quinazoline derivatives are a class of novel, potent, selective and orally bioavailable PLK1 inhibitors with reasonable SAR and strong quantitative correlations. The CoMFA and CoMSIA studies were performed on these compounds based on three different alignment methods to build 3D-QSAR models. Most of the models showed good q 2 and r 2 pred values and revealed a good response to the test set validation. The CoMFA model generated from the pharmacophore-based method was found to be superior (model 2, q 2 = 0.628, r 2 pred = 0.785) to those obtained from GLIDE docking-based and common substructure based methods. All the CoMSIA models derived from three different alignment methods gave good results, whose q 2 and r 2 pred values were greater than 0.5 and 0.6, respectively. The q 2 value of the best CoMSIA model was only a little larger than that of other models. In view of that, three CoMSIA models were selected for further comparisons and analyses so that more valuable information for the structural requirements can be obtained. From our studies, it was found that the pharmacophore-based alignment produced the best model for CoMFA and the common structure-based alignment for CoMSIA. This indicated that the discovery of the optimal alignment method should depend on the statistical performances of 3D-QSAR models generated from the alignments based on those methods. In addition, suitable alignment methods for CoMFA and CoMSIA studies might be different. The comparative studies among the best CoMFA and the CoMSIA models were also demonstrated in the crystallographic environment of PLK1 and high consistency was found in steric, electrostatic and hydrogen bond donor fields. Furthermore, the contours of CoMSIA model 6 were compared with the structure-based pharmacophore model and the key factors related to binding affinity were reconfirmed. These satisfactory insights identified in the present study can be utilized to design and predict new potent compounds as PLK1 inhibitor candidates, and to discover compounds with novel scaffolds that can act as PLK1 inhibitors via similar mechanisms.