Structural studies of B-type Aurora kinase inhibitors using computational methods.

AIM
To characterize the structural features of quinazoline-based Aurora B inhibitors that influence its inhibitor activity.


METHODS
Two geometrical methods, Method 1 and Method 2, were used to develop the 3D-QSAR models. The most active ligand was used as the template for the alignment of all the ligands in Method 1, and a conformer of the cocrystal ligand was used as the template for the alignment of all the ligands in Method 2.


RESULTS
The models suggest that highly active ligands can be designed by varying the R1 substituent at position 7 of the quinazoline ring with positively charged, bulky, hydrophobic groups, while bulky and hydrophobic groups around the thiazole ring are desirable for higher activity.


CONCLUSION
This study emphasizes that the bioactive conformer is rather different from the minima. The steric, electrostatic, and hydrophobic field effects contribute to its inhibitory activity.

outcomes in patients with endometrial carcinoma [19] . The forced expression of Aurora B in Chinese hamster embryo cells results in chromosomal instability and increased tumor invasiveness in association with the constitutive expression of phosphorylated (p)-histone H3 on Ser10 [20] , suggesting that Aurora B can act as an oncogene. The development of Aurora kinase inhibitors evidence linking Aurora kinases to malignancies has raised the possibility of targeting these kinases for cancer therapy.
Many Aurora inhibitors in various stages of development have been reported [21] . Among them, MK-0457/VX-680, MLN8054, ZM447439, PHA-739358, AZD1152, and AT9283 are in clinical development. These compounds are of different scaffolds with various specificities to the subtypes of Aurora kinases. VX-680 is a potent inhibitor of all three Aurora kinases with K i values of 0.6, 1.8, and 4.6 nmol/L for the A, B, and C isoforms, respectively [22] , and is undergoing clinical trials for solid tumors and hematological malignancies. PHA-739358, a pyrazole derivative, is a non-selective agent with IC 50 values of 13, 79 and 61 nmol/L for Aurora A, B, and C, respectively [23] . AZD1152 and ZM447439 are both quinazoline derivatives in clinical stages. AZD1152 selectively inhibits Aurora B rather than A or C and resulted from the optimization of a series of 5-acetanilide-3-aminopyrazole (3-pyrazole)substituted quinazolines [24] . Recently, a new class of 1-acetanilide-4-aminopyrazole-substituted quinazoline has been reported as selective Aurora B inhibitors [25] . These findings suggest the possibility of designing more potent quinazolinebased Aurora B inhibitors. The three-dimensional (3D) quantitative structure-activity relationship (3D-QSAR) techniques, such as the comparative molecular field analysis (CoMFA) and comparative molecular similarity analysis (CoMSIA) [26][27][28] , are routinely used in modern drug design to understand the drug-receptor interaction. It has been shown in the literature that these computational techniques can strongly support and aid in the design of novel, more potent inhibitors by revealing the mechanics of the drug-receptor interactions [29][30][31] . To explore the further possibilities of novel drugs, we have developed predictive 3D-QSAR models using quinazoline-based Aurora B inhibitors [24,32] .

Materials and methods
The basic structure of the quinazolines is shown in Figure  1. A series of 48 quinazoline derivatives with their inhibitory activities was retrieved from the literature [24,32] . We have considered only the IC 50 value for the cell-based assay (not an enzyme-based assay), which is the inhibition of histone H3 in SW620 cells. Histone H3 is a cellular substrate of Aurora B and a marker of Aurora B kinase inhibition in vivo. The IC 50 values in the micro-molar (µmol/L) range were converted into the molar (M) range and then into pIC 50 values using the relationship pIC 50 =-logIC 50 values. The pIC 50 values were used as the dependent variables in the statistical modeling. The dataset of Aurora B kinase inhibitors was divided into training (37 compounds) and test (11 compounds) sets considering the fact that the test set molecules must present a range of biological activity and the typical chemical structures similar to those of the training set as reported in Table 1. The test set compounds  are marked by asterisks in Table 1. Molecular modeling and all calculations were performed using SYBYL 7.3 [33] running on Linux cluster. The molecules were aligned using two different geometric methods.

Molecular Alignment
In the 3D-QSAR model development, molecular alignment is a key step, and usually the molecules are aligned with a suitable conformational template that is supposed to be bioactive so that the inhibitors have comparable conformations and similar orientations in space [26] . In this study, two different geometrical methods, Method 1 and Method 2, were used.

Method 1
The most active molecule (compound 40) was drawn and optimized using Tripos Force Field (TFF) [34] with a convergence criterion of 0.005 kcal/mol/Å and 10 000 iterations. The random search-based minimum energy conformer of the most active compound (compound 40) was obtained and used as the template. This template was modified for further ligands, and all of the molecules were optimized using the TFF level with a convergence criterion of 0.005 kcal/mol and 10 000 iterations. These minimized structures were aligned over the template using the common substructure method ( Figure 2) and subsequently used for the CoMFA/CoMSIA.

Method 2
A recently reported crystal structure for a similar ligand (PDB=2c6e) was obtained from the public domain. The ligand position within the active site is displayed in Figure 3. The docked ligand was extracted and used as template to align all of the ligands. The ligands were all minimized within the receptor site using TFF, but the whole protein active site was fixed during the minimization. All of the minimized structures were aligned over the template using a common substructure method and directly used for CoMFA and CoMSIA. The aligned structures of Method 2 are displayed in Figure 4.

CoMFA
The steric and electrostatic potential fields for the CoMFA were calculated at each lattice intersection of a regularly spaced grid of 2.0 Å. The lattice was defined automatically and extended up to four units past the Van der Waals volume of all molecules in the X, Y, and Z directions. The Van der Waals potential and columbic terms, which represent the steric and electrostatic fields, respectively, were calculated using TFF and a distance-dependent dielectric constant. An sp 3 carbon atom with a Van der Waals radius of 1.52 Å and a charge of +1.0 served as the probe atom to calculate the steric and electrostatic fields. The steric and electrostatic contributions were truncated at ±30 kcal/mol; the electrostatic contributions were ignored at the lattice intersections with the maximum steric interactions. The CoMFA steric and electrostatic fields generated and scaled by the CoMFA standard option are given in SYBYL [33] . The structures A and B correspond to Figure 1  The reported CoMSIA method is based on molecular similarity indices [27] with the same lattice box used for the CoMFA. Molecular similarity expressed in terms of five different properties (steric, electrostatic, hydrophobic, H-bond donors and H-bond acceptors) was calculated using a C + probe atom with a radius of 1.0 Å placed at a regular grid spacing of 2.0 Å. The CoMSIA similarity indices (A F ) for molecule j with atoms i at grid point q are calculated by A q F,K (j)=-∑ω prob,k ω ik e -αr2 iq (1) where k represents the following physicochemical properties: steric, electrostatic, hydrophobic, H-bond donor, and H-bond acceptor. A Gaussian type distance dependence was used between grid point q and each atom i of the molecule. The default value (0.3) was used as the attenuation factor (α). The steric indices are related to the third power of the atomic radii, the electrostatic descriptors are derived from the atomic partial charges, the hydrophobic fields are derived from atom-based parameters [35] , and the H-bond donor and acceptor indices are obtained by a rule based method based on experimental results [36] .
Partial Least Square (PLS) analysis and validation of the QSAR models To derive the 3D-QSAR models, the CoMFA and CoMSIA descriptors were used as independent variables and the pIC 50 values as the dependent variable. The PLS method [37,38] was used to linearly correlate the CoMFA and CoMSIA descriptors to the activity. The CoMFA cutoff value was set to 30 kcal/mol for both the steric and the electrostatic fields, and all of the fields were scaled by the default options in SYBYL. The cross-validation analysis was performed using the leave one out (LOO) method [39,40] in which one compound was removed from the dataset and its activity was predicted using the model derived from the rest of the dataset. The crossvalidated correlation coefficient (q 2 ) that resulted in an optimum number of components and the lowest standard error of prediction was considered for further analysis and calculated using where, y pred , y actual , and y mean are the predicted, actual and mean values of the target property (pIC 50 ), respectively, and PRESS is the sum of the predictive sum of squares. The non-crossvalidated PLS analyses were performed with a column filter value of 2.0 to reduce the analysis time with a small effect on the q 2 values. To have robustness and statistical confidence of the derived models, a bootstrapping analysis was used for 100 runs. To assess the predictive power of the derived 3D-QSAR models, the activity of the test set of eleven molecules was predicted. The predictive ability of the models is expressed by the predictive r 2 value, which is analogous to the crossvalidated r 2 (q 2 ) and is calculated using (4) where SD is the sum of the squared deviations between the biological activities of the test set and the mean activities of the training molecules and PRESS is the sum of the squared deviation between the predicted and observed activity of the test set molecules and is calculated using equation 4.

Results
The 3D-QSAR analyses were performed by two different geometric methods using ligand-based and receptor-guided techniques.
A comparison of the two geometric methods The structures of conformers A (green) and B (red) of the most    Figure 5. Structure A was obtained from a minimum energy conformer of the most active compound, while structure B was drawn from a modification of the crystallized ligand (gray) and further minimized within the active site. It is clear ( Figure 5) that the receptor guided conformer (red) is closer to residues 211-223, and the key interaction with residue 161 is more significant for compound 40. For this interaction, the distance between the carbonyl carbon and the hydrogen of Lys was reported as 1.78 Å for the co-crystallized ligand, while this distance for the most active ligand (compound 40) was 1.2 Å. The minimum energy conformer (green) has a distance of 3.5 Å. The molecular field analyses and this comparison of different conformers emphasize that the bioactive conformer is not the minimum energy conformer.

Comparative molecular fields analyses (CoMFA)
A Gasteiger Huckel charge was used for the aforementioned different geometrical methods. Six models (CoMFA 1 to CoMFA 6 ) were developed. CoMFA 1 to CoMFA 3 were based on minimum energy geometries with steric, electrostatic, and both fields combined, respectively. The statistical summary of these three models is presented in Table 2. Model 3 (CoMFA 3 ) was better than the other two models because it was based on both fields having high LOO values of q 2 =0.54 and r 2 =0.86. The electrostatic contribution was dominant in the interaction. The model was successfully validated for internal predictivity with r 2 bs =0.90 and the test-set predictivity of 11 compounds with r 2 predictive =0.74. CoMFA 4 to CoMFA 6 were based on the co-crystal geometry with steric, electrostatic, and both fields combined, respectively. The result of both geometries was good; the co-crystal geometry-based model was statistically significant ( Table 2). The best-fit model (CoMFA 6 ) was based on the X-ray geometry and both fields jointly have LOO values of q 2 =0.56 and r 2 =0.96. Even though this model involved both steric and electrostatic interactions, the electrostatic contribution was dominant in the interaction. The model was successfully validated for the internal predictivity with r 2 bs =0.97 and for the test set predictivity of 11 compounds with r 2 predictive =0.76. The q 2 value of model CoMFA 3 was better than that of model CoMFA 6 , but the other values indicated the significance of model CoMFA 6 . This model is more realistic because it is based on the co-crystallized binding mode in which all molecules have nearly the same conformation as in the co-crystal structure, and all have a key contact hydrogen bonding between residue 161 and the carbonyl carbon. All molecules acquired the same binding site and an almost equal interaction with the hinge region. The predicted activities by the CoMFA are summarized in Table 3.
The trend of the observed and predicted activities using model CoMFA 6 is displayed in Figure 6.

CoMFA maps
The 3D-CoMFA contour maps using the best-fit model CoMFA 6 are displayed in Figures 7 and 8 with both kinds of field effects on compound 40. The green contours indicate the area in which steric bulk substations might beneficially affect the activity and the yellow region is favorable for small groups ( Figure 7). The blue contour indicates the region where a positive group is required for high activity, while the red zone indicates the region favorable for negative groups ( Figure 8).
The green contour was evident around the pyrrolidine ring, which is directed toward residues 211 to 223, and around the thiazole ring, which is directed toward 207-209 and 161-162, indicating that a bulkier group around these regions will be favorable for higher activity. Particularly, the contact with residue 161 was described as a key interaction for activity. A small yellow contour was evident just near the methoxy group, indicating that a small group would be favorable for higher activity. Blue contours appear around the methoxy and pyrrolidine rings, as well as around the thiazole ring (Figure 8), which indicated that except for the methoxy region, a positive group with steric bulk would be desirable at this position, which is in support of previous findings [41] . A small positive group would be helpful around the methoxy substituent. There is a red contour in the far region of the thiazole ring indicating that a negative group would have a good effect on the activity, which is also clear from the data of the most active compound (compound 40), which has halogens at the same place that are probably key to its higher activity.  Figure 6. The trend of observed and predicted activities for the training and test set using CoMFA based model.      Table 4. The model (CoMSIA 9 ) based on the steric, electrostatic, and hydrophobic effects provided good results, with LOO values of q 2 =0.75 and r 2 =0.98. For this model, the correlation coefficient for the internal predictivity was r 2 bs =0.98 and the correlation coefficient for the test set predictivity of 11 compounds was r 2 predictive =0.73. Similarly, the models CoMSIA 17 to CoMSIA 32 are based on the X-ray geometries and combinations of different field effects. The best-fit model (CoMSIA 25 ) was based on the X-ray geometry and the steric, electrostatic and hydrophobic fields effects, which provided LOO square of cross-validated correlation coefficients of q 2 =0.77 and r 2 =0.98. The model was tested for the internal predictivity (r 2 bs =0.98) and test set predictivity of 11 compounds (r 2 predictive =0.77). The statistical value of model CoMSIA 25 was better than that of model CoMSIA 9 , which indicated the significance of model CoMFA 25 . The predicted activities by model CoMSIA 25 are summarized in Table  3. The trend of the observed and predicted activities using model CoMSIA 25 is displayed in Figure 9.

CoMSIA maps
Like the CoMFA, the CoMSIA contour maps were also developed using the best-fit model CoMSIA 25 , which was based on steric, electrostatic, and hydrophobic interactions. Each  (Figures 10-12) with compound 40. The green contours indicate the area in which steric bulky substations might beneficially affect the activity and the yellow region is favorable for small groups ( Figure 10). The map is quite similar to the CoMFA steric map because the green contour appears around the pyrrolidine ring and other green contours appear around the thiazole ring, which indicates that a bulkier group will be favorable for higher activity. The blue contour indicates the region where a positive group is required for high activity, while the red zone indicates the region favorable for negative groups (Figure 11). Blue and red contours appeared around the thiazole ring, which is similar to the CoMFA electrostatic map. These contours indicate a suitable site for positive and negative substituents. The CoM-SIA map with the hydrophobic field effect is displayed in Figure 12. The white contour favorable for hydrophobic groups was evident around the pyrrolidine and near the thiazole ring. The small yellow contour favorable for the hydrophilic group appeared at a certain distance from the thiazole ring. A comparison indicated that there is a dominance of the hydrophobic and bulky groups around the same region (Figures 10 and  12).

Discussion
In earlier publications, a number of quinazoline-based series have been described as Aurora kinase inhibitors [2,22,[42][43][44] . Particularly, the p-benzamidoanilinoquinazoline series was shown to inhibit Aurora A and Aurora B equipotently. The investigation of the SAR around the anilino ring linking the key hinge (quinazoline) and the lipophilic pocket (benzamide) binding groups showed that this ring could be replaced with a number of heterocyclic systems to give improvements in both the potency and the physicochemical properties [42] . The linker group also has an important role for the cellular activity. For example, the ether and thiomethyl groups in compound 35 and 36 gave poor cellular potency compared to the acetanilide series. Other linkers like the reverse amide and sulfonamide in compounds 37 and 38 also are responsible for their lower cellular activity. Comparing all the linkers, it is clear that acetanilide shows the potent cellular activity that is present in the most active compound 40. Furthermore, it was found that Aurora A is able to accommodate highly polar functional groups within the side chain at the quinazoline C-7 position. Compounds with elaborated C-7 side chains show good cellular activity and benefit from much better solubility and free drug levels. The quinazoline C-7 side chain might result in a further enhancement of the cellular potency [42] . In the thiazole series, the length of the aliphatic spacer, between the oxygen atom on the C-7 position of the quinazoline ring and the basic nitrogen atom within the side chain of the R1 substituent, governs the ability of the polar group to solvate [42] . Both cyclic and noncyclic amide chains are tolerated by the Aurora kinases.    Among the five and six membered ring compounds, the prolinol side chains have greater cellular activity. The preference for five membered rings can be helpful to support the findings of the 3D-CoMFA model, which suggest a bulkier substituent in this position. The structure-activity relationships of the pyridino-and pyrimidinoquinazoline compounds have been described recently and include one of the first crystal structures with a small inhibitor bound to the Aurora kinase [45] . It was helpful in rationalizing the observed high potency and high specificity for the Aurora kinases. Both the CoMFA and the CoMSIA models suggested that small and lipophilic groups were particularly favorable in the terminal lipophilic binding pocket. Fluoro groups at the ortho and meta positions were found to be particularly favorable for better activity. The 2,3-difluoro and 3-fluoro acetanilide groups generally give rise to better cellular potency compared to larger and more polar functional groups. These findings can be observed from the comparison of compounds 5 and 3 with compound 8 in the pyrazole series. The comparison of compound 47 with the mentioned compounds provides evidence for the binding conformation for both (pyrazole and thiazole) series. These studies suggest that a further hydrogen bond donor in the pyrazole ring (compared to the thiazole ring) may allow additional binding to the protein backbone. The improved physical properties, together with reasonable enzymatic and cellular activities, made the pyrazole series an attractive proposition for further medicinal chemistry aimed at optimizing the cellular potency, while at the same time maintaining the drug-like properties.
A number of crystal structures of Aurora A in different conformations with ATP analogues have been published [45][46][47][48] . A crystal structure of Aurora A complexed with a small molecule inhibitor from the pyrimidinoquinazoline series adopts an inactive conformation with the displacement of the DFG motif (DFG-out) [45] . In the crystal structure ( Figure 5), the distance between the carbonyl oxygen (black skeleton) and the hydrogen of the amino group was 1.8 Å. While using molecular docking, the best active compound acquired the same site as the X-ray structure and apparently the noted distance was 1.2 Å, which facilitates the interaction and results in a more potent inhibitor. Similar binding interactions may exist with the acetanilide chain; however, the interposed methylene and the inversion of the amide may lead to a different binding mode within the hydrophobic pocket that could translate into higher cellular potency.
Compound 40 binds in the ATP-binding site. The quinazoline ring of compound 40 fit into the hinge region ( Figure  13A), and the N-atom formed a hydrogen bond with the main chain of Ala212. In addition, the quinazoline ring also formed hydrophobic interactions with the surrounding residues, including Leu138, Val146, Ala159, and Leu262. In particular, Leu262 made extensive interactions with the quinazoline ring. The prolinol moiety extended into the solvent-exposed region of the binding site. The amide carbonyl near the thiazole shows a hydrogen bond to the conserved Lys161. The terminal aromatic ring fits into a hydrophobic pocket of the kinase and is mainly surrounded by hydrophobic residues, Trp276, Leu163, Leu177, and Leu207. The pyrazole ring and amino nitrogen near the quinazoline of compound 8 form hydrogen bonds ( Figure 13B) with the hinge residues (Glu210 and Ala212, respectively). The O-atom of 3-methoxy acetamide shows a hydrogen-bonding interaction. This may cause inhibitor 8 not to occupy the terminal lipophilic binding pocket fully, which is important for cellular potency.
The Thr217/Glu161 difference has been proposed recently to explain the selectivity of indirubin-based inhibitors for Aurora B [49] . A docking study was performed to observe the ligand-receptor interaction regarding the selectivity for Aurora B. The docking data showed that compound 40 could fit into the catalytic cleft of Aurora B (PDB ID 2VGP) and bind strongly to the residues surrounding the cleft ( Figure 14A). The elaborate docking indicates that compound 40 occupies the ATP binding pocket of the active center of Aurora B with multiple interactions with the residues around the pocket: the quinazolino nitrogen and amino hydrogen are hydrogen bonded to Ala173 and the amide hydrogen to Pro174. At the same time, the hydroxyl of prolinol binds to Glu177 with a hydrogen bond, occupying the entry site to the catalytic cleft in Aurora B ( Figure 14B). This binding mode is similar to the recently reported [50] Aurora B inhibitor. This interaction can be explained by considering the property of the glutamate residue. Glutamate is an acidic residue with a longer and more flexible side chain and is located near the catalytic edge of the catalytic pocket, in a subsite referred to as the ribose binding region. It is a hydrophilic part of the binding site, where interactions between the receptor and the sugar moiety of the natural substrate are formed, as demonstrated in crystal structures such as adenosine complexed with cAMP-dependent kinase (1FMO) or ATP complexed with CDK2 (1HCK). Moreover, crystallographic studies revealed that the direct or water mediated interaction of the inhibitor with residues located in the ribose-binding site enhances the ligand affinity. The superimposition of Aurora A (PDB ID 2C6E) and B (PDB ID 2VGP) shows a clear difference between the active site residues (Figure 15). One of the Aurora B residue equivalents (Lys 122) to Aurora A (Lys 161) lacks a hydrogen bond interaction because of its different orientation. Also, the position of the nonconserved residue Glu177 explains the distinct interaction responsible for the selectivity. The Glu177 residue is closer to the hydroxyl of prolinol to make the hydrogen bond interaction. This type of important interaction revealed from the docking study can be a useful guideline for future selective inhibitor design.
Modification of the R1 substituent at the C7 position of the quinazoline inhibitors is known to profoundly influence their activity and selectivity [51][52][53][54] . The R1 substituent accesses the solvent channel, explaining why a wide range of strongly polar groups is tolerated, provided that a suitable linear spacer is included between the solubilizing group and the quinazoline. Variation of the R1 substituent at the quinazoline C7 side chain allowed the fine-tuning of the overall properties of the compounds. From the above discussion, it can be said that www.chinaphar.com Neaz MM et al Acta Pharmacologica Sinica npg variation of the R1 substituent at the quinazoline C7 position along with an optimized linker and terminal aromatic substituent may help to design more potent and selective Aurora B inhibitors. Beside these findings, a progressive scrambling technique was employed to determine the possibility of a chance correlation in the QSAR model. The cSDEP values are almost at a minimum and the Q 2 values are at a maximum (Tables 5 and 6). The values for the slope in the 2, 3, and 4 component models are acceptable, while a value greater than 1.2 is indicative of a bad model. Thus, we used the CoMFA and CoMSIA models having the least number of components for a stable QSAR model.

Conclusion
The overall study implies that steric and hydrophobic effects, along with a positive electrostatic interaction, play crucial roles in the inhibitory activity of quinazoline derivatives. The activity is quite sensitive to the conformation, and the contact between the carbonyl oxygen and the Lys hydrogen is important to the activity. The importance of steric bulk with electronic interaction is evidenced by the statistics. The ligand-   npg based 3D QSAR model has proven significant, but a more definitive conclusion requires consideration of the receptor site. The receptor-guided CoMFA model has a high value of LOO, q 2 =0.56, r 2 =0.96 and r 2 predictive =0.76; similarly, the CoMSIA model has a high value of LOO, q 2 =0.77, r 2 =0.98 and r 2 predictive =0.77, suggesting that a highly active ligand can be designed by variations of the R1 substituent at position 7 of the quinazoline ring with positively charged, bulky and hydrophobic groups, while around the thiazole ring, bulky and hydrophobic groups are desirable for higher activity. Also, the formation of a hydrogen bond with Glu177 can be useful for selectivity.