Pyrrolo[3,2-d]pyrimidine Derivatives as Type II Kinase Insert Domain Receptor (KDR) Inhibitors: CoMFA and CoMSIA Studies

Kinase insert domain receptor (KDR) inhibitors have been proved to be very effective anticancer agents. Molecular docking, 3D-QSAR methods, CoMFA and CoMSIA were performed on pyrrolo[3,2-d]pyrimidine derivatives as non-ATP competitive KDR inhibitors (type II). The bioactive conformation was explored by docking one potent compound 20 into the active site of KDR in its DFG-out inactive conformation. The constructed CoMFA and CoMSIA models produced statistically significant results with the cross-validated correlation coefficients q2 of 0.542 and 0.552, non-cross-validated correlation coefficients r2 of 0.912 and 0.955, and predicted correction coefficients r2pred of 0.913 and 0.897, respectively. These results ensure the CoMFA and CoMSIA models as a tool to guide the design of a series of new potent KDR inhibitors.


Introduction
Angiogenesis, the form of new blood vessels by capillary sprouting from pre-existing vasculatures, is a normal process for organ development during embryogenesis, wound healing and female reproductive cycling [1,2]. Abnormal regulation of angiogenesis has been shown to get involved in OPEN ACCESS many diseases, such as diabetic retinopathy, psoriasis, rheumatoid arthritis, and cancer. In particular, it is widely recognized that the growth and metastasis of solid tumors is dependent on angiogenesis [3,4]. Out of the many factors that are involved in angiogenesis, vascular endothelial growth factors (VEGFs) are of particular interests [5,6]. The VEGFs are required for vasculogenesis and angiogenic sprouting, and act through receptor tyrosine kinase (VEGFR-1, -2, and -3). Among them, VEGF receptor-2 (VEGFR-2) or kinase insert domain receptor (KDR) play crucial roles in vessel sprouting and new vessel initiation in early stage of angiogenesis.
It is well known that inhibiting of KDR leads to suppression of angiogenesis and tumor growth. A number of preclinical and clinical studies have shown that many small-molecule KDR inhibitors are capable of inhibiting angiogenesis, tumor progression, and dissemination [7][8][9][10][11]. The vast majority of the KDR inhibitors known to date, such as Gefitinib, are ATP-competitive and classified as type I inhibitors. Such inhibitors target the ATP binding pocket in its active conformation of the activation loop. This conformation is normally referred to as DFG "in" based on the position of the conserved triad aspartate-phenylalanine-glycine (DFG) at the entrance of the activation loop. Type I inhibitors typically function in the DFG "in" conformation of KDR through hydrogen bonding with the backbone residues of the hinge region as well as hydrophobic interactions in and around the adenine region. These features, however, make it quite difficult to design highly selective type I inhibitors [12,13]. In addition, ATP-competitive inhibitors have to compete with high levels of intracellular ATP, leading to a significant difference between the in vitro and in vivo activities. In response to these issues, non-ATP competitive kinase inhibitors, such as Imatinib have been identified. These inhibitors, i.e., type II inhibitors bind to and stabilize an inactive kinase form that features the DFG motif in an "out" conformation. The different position of the DFG residues in the "out" form creates a new hydrophobic binding pocket that is adjacent to the ATP-binding site. This pocket, also known as the allosteric site, is characteristic of kinase in an inactive conformation. Type II inhibitors predominantly occupy the ATP binding site, but they also exploit unique hydrogen bonding and hydrophobic interactions with the allosteric site. Compared with type I kinase inhibitors, type II inhibitors have several advantages, including great cellular potency and improved kinase selectivity. In addition, type II inhibitors, because of their interactions with both the ATP pocket and the allosteric site, may provide an avenue to overcome the mutations that induce resistance to the other types of inhibitors [14][15][16][17][18].
Recently, Oguro et al. reported a set of pyrrolo[3,2-d]pyrimidine derivatives 1-52 (Table 1) as potent and selective KDR inhibitors in biochemical and cellular assays [19,20]. These compounds are structurally distinct from the ATP-competitive inhibitors, and the X-ray crystallographic analysis has shown that compound 27 induces an inactive conformation, in which the diphenylurea moiety occupies the hydrophobic pocket created by the conformation change of DFG motif (DFG-out) [19,20]. In the present study, to gain further insight into the interactions of these compounds with the KDR, we employed the molecular docking-guided three dimensional quantitative structure-activity relationship (3D-QSAR) study, comparative molecular field analysis (CoMFA) [21] and comparative molecular similarity indices analysis (CoMSIA) [22] to address how steric, electrostatic, hydrophobic and hydrogen-bonding interactions modulate the inhibitory activities.

General
The crystallographic coordinates of KDR in complex with small-molecule inhibitors were obtained from the Brookheaven Protein Databank as entries 2OH4 [23]. All the molecular modeling and calculations were performed using Sybyl 7.3 molecular modeling package [24].

Data Sets
In this study, 52 pyrrolo[3,2-d]pyrimidine derivatives were taken from the work of Oguro et al. [19,20]. Their structures and inhibitory activities are listed in Tables 1 and 2. The IC 50 values (M) were converted to the corresponding pIC 50 (= −logIC 50 ). On the basis of the diversity in the structures and activities, these 52 compounds were divided into two groups: 44 compounds were used as the training set to build the 3D-QSAR models, and 8 compounds that are marked with an asterisk in Table 1 were used as the test set to evaluate the predictive power of the developed CoMFA and CoMSIA models.  Because the crystal structure of KDR in complex with pyrrolo[3,2-d]pyrimidine was not available in the Brookheaven Protein Databank (PDB), the bioactive conformation was simulated by docking using Surflex-dock program. The crystallographic coordinates of KDR in complex with its inhibitor, which was reported to be in the inactive DFG-out conformation of KDR, were obtained from the PDB as entries 2OH4 [23]. Surflex-Dock program [25,26] has been widely used to calculate the protein-ligand interactions, and to efficiently predict the active conformations [27][28][29][30][31][32][33][34][35]. Surflex-Dock uses a Protomol-based method and an empirical scoring function to dock a ligand into the binding site of a receptor. The Protomol is an idealized representation of a ligand that forms every potential interaction with the binding site. Surflex-Dock's scoring function contains the factors that play crucial roles in the ligand-receptor interaction, including hydrophobic, polar, repulsive, entropic and solvation terms. In this study, the Protomol was generated using a ligand-based approach. During the Protomol generation process, two particular parameters, Protomol_bloat and Protomol_threshold, were specified to form the appropriate binding pocket. The former determines how far the site should extend from a potential ligand, whereas the latter determines how deep the atomic probes that are used to define the Protomol can penetrate into the protein. In the present work, Protomol_bloat and Protomol_threshold default values (0 and 0.5, respectively) were used when a reasonable binding pocket was obtained. During the docking process, the default values of all the other parameters were assigned. The highest-scored conformation of a potent compound 20 based on the Surflex-Dock scoring functions, was selected as the final bioactive conformation.

Molecular Modeling
In the 3D-QSAR study, the selection of active conformations is a key step for CoMFA and CoMSIA studies. The bioactive conformation of compound 20 was simulated using Surflex-Dock. The docked conformation with the highest total score was used as the template to construct the 3D structures of the rest compounds in the data set. Structural energy minimization process was performed using the Tripos force field with a distance-dependent dielectric and Powell gradient algorithm with a convergence criterion of 0.001 kcal/mol. Partial atomic charges were calculated using Gasteiger-Hückel method.

Molecular Alignment
In the 3D-QSAR study, the alignment rule is also a key step. The predictive accuracy of the CoMFA and CoMSIA models and the reliability of the contour maps are directly dependent on the structural alignment rule. The compounds were aligned by the atomfit to the template 20. The aligned compounds are shown in Figure 1.

CoMFA and CoMSIA Studies
Standard CoMFA and CoMSIA procedures were performed. A 3D cubic lattice was created automatically by extending at least 4 Å beyond all the aligned molecules in X, Y and Z directions with 2.0 Å grid spacing. The CoMFA steric (Lennard-Jones potential) and electrostatic (Coulomb potential) fields at each lattice were calculated using the standard Tripos force field method. A distance dependent dielectric constant of 1.0 was used, and an sp 3 hybridized carbon atom with one positive charge and a radius of 1.52 Å served as a probe atom to calculate the steric and electrostatic fields. The default cutoff value of 30.0 kcal/mol was adopted.
Compared with CoMFA, CoMSIA methodology has the advantage of exploring the impacts of more fields. In addition to the steric (S) and electrostatic (E) fields used in CoMFA, the CoMSIA method defines hydrophobic (H), hydrogen bond donor (D), and hydrogen bond acceptor (A) descriptors. The CoMSIA fields were derived, according to Klebe et al. [22], from the same lattice box that was used in the CoMFA calculations, with a grid spacing of 2 Å and a probe carbon atom with one positive charge and a radius of 1.0 Å as implemented in Sybyl. Arbitrary definition of cutoff limits was not required in CoMSIA method, wherein the abrupt changes of potential energy near the molecular surface were taken into account in the distance dependent Gaussian type functional form. The default value of 0.3 was used as the attenuation factor.

PLS Regression Analysis and Validation of QSAR Models
Partial least squares (PLS) approach was used to derive the 3D QSAR models. The CoMFA and CoMSIA descriptors were used as independent variables and the pIC 50 values were used as dependent variables. CoMFA and CoMSIA column filtering was set to 2.0 kcal/mol to improve the signal-to-noise ratio. The leave-one-out (LOO) cross-validation was carried out to obtain the optimal number of components (N) and the correlation coefficient q 2 . The obtained N was then used to derive the final QSAR model and to obtain the non-cross-validation correlation coefficient r 2 , standard error of estimate (SEE), and Fischer (F) ratio value.
To assess the predictive power of the derived 3D-models, a set of test compounds that had known biological activities, was used to validate the obtained models. The predictive correlation r 2 pred value was calculated using Wherein SD is the sum of the squared deviations between the biological activities of the test compounds and the mean activities of the training compounds, and PRESS is the sum of the squared deviations between the experimental and the predicted activities of the test compounds.

Binding Modes of Pyrrolo[3,2-d]pyrimidine Derivatives
To determine the probable binding conformations of these compounds, Surflex-Dock was used to dock one potent compound 20 into the active site of KDR (PDB code: 2OH4). First, the docking reliability was validated by a known inhibitor 53 ( Figure 2) that was reported to bind in the DFG-out inactive conformation of KDR [23]. The co-crystallized 53 was re-docked into the binding site, and the docked conformation having the highest total score was selected as the most probable binding conformation ( Figure 3). The low root mean-square deviation (RMSD) of 0.58 Å between the docked and the crystal conformations demonstrated the high reliability of Surflex-dock in reproducing the experimentally observed binding mode for these KDR inhibitors. As shown in Figure 3a, redocked 53 was almost in the same orientation with co-crystallized 53 at the active site of KDR. Therefore, Surflex-Dock docking protocol and the used parameters were extended in search for the binding conformations of KDR inhibitors. The binding mode of compound 20 at the active site of KDR in its inactive conformation is shown in Figure 3b. It is clear that compound 20 adopts an overall conformation that is very similar to that of compound 53, and forms multiple hydrogen bonds with KDR. Specifically, the 1-nitrogen of the pyrrolo[3,2-d]pyrimidine subunit forms one hydrogen bond with the backbone-NH of Cys917 in the hinge region of KDR at the angle and distance of 130.76° and 3.253 Å, respectively. In the allosteric site, the two NH groups of the urea moiety form two hydrogen bonds with the side chain carboxylate of Glu883 at the angles of 99.12° and 97.85°, and the distances of 2.566 Å and 2.873 Å, respectively, whereas the CO moiety interacts with the backbone-NH of Asp1044 in the DFG motif through hydrogen bonding at the angle and distance of 151.02° and of 2.666 Å, respectively. These interactions are typical characteristics of the interactions between the inhibitors and the inactive conformation of the kinase. In addition, the terminal phenyl moiety occupies the hydrophobic pocket that is formed from residues Ile886, Leu887, Ile890, Val896 and Leu1017.

CoMFA and CoMSIA Results
The CoMFA and CoMSIA 3D-QSAR models were derived from a training set of 44 compounds. The statistical results of the CoMFA and CoMSIA 3D-QSAR models are presented in Table 3. The CoMFA model gave a cross-validated correlation coefficient q 2 of 0.542, an optimal number of principal components (N) of 4 and a non-cross-validated correlation coefficient r 2 of 0.912. The corresponding contributions of steric and electrostatic fields were 52.5% and 47.5%, respectively. The CoMSIA model gave a cross-validated correlation coefficient q 2 of 0.552, an optimal number of principal components of 5 and a non-cross-validated correlation coefficient r 2 of 0.955. The corresponding contributions of steric, electrostatic, hydrophobic, hydrogen bond donor and acceptor fields were 18.4%, 22.8%, 34.3%, 6.3% and 18.2%, respectively. Both the CoMFA and CoMSIA models were satisfactory from the viewpoint of statistical significance. The activities of the 44 training compounds were predicted with the constructed CoMFA and CoMSIA models. The predicted pIC 50 values are shown in Table 2 and Figure 4. It can be seen that the predicted pIC 50 values were in good agreement with the experimental values, indicating that the obtained CoMFA and CoMSIA models had strong predictive ability.

Validation of the 3D-QSAR Models
The predictive powers of the CoMFA and CoMSIA models were validated by the eight test compounds. The predicted pIC 50 values were found to be in good agreement with the experimental data within an acceptable error range (Table 2 and Figure 4). The predictive correction coefficients of the CoMFA and CoMSIA models were 0.913 and 0.897, respectively. This result indicates that the CoMFA and CoMSIA models may be used to predict the inhibitory activities of novel pyrrolo [3,2-d]pyrimidine derivatives as type-II KDR inhibitors.

Contour Analysis
To visualize the results of the CoMFA and CoMSIA models, the 3D coefficient contour maps were generated. The CoMFA and CoMSIA results were graphically interpreted by the field contribution maps using the STDEV*COEFF field type. The contour maps of CoMFA (steric and electrostatic) and CoMSIA (steric, electrostatic, hydrophobic, hydrogen bond donor and acceptor fields) are shown in Figures 5 and 6, respectively. Compound 20 was displayed in the map in aid of visualization. All the contours represented the default 80% and 20% level contributions for favorable and unfavorable regions, respectively, except 90% and 10% level contributions in Figure 6c, respectively. The CoMFA contour maps of the steric and electrostatic fields are shown in Figure 5. In the map of steric field, the green contours represent the regions in which bulky groups confer an increase in the activity, whereas the yellow ones represent the regions where bulky groups may lead to a decrease in the activity. Similarly, in the map of electrostatic field, the blue contours indicate the regions where electropositive substitution increases the inhibitory activity, whereas the red contours indicate the regions where electronegative substitution increases the activity.
In the CoMFA steric contour map (Figure 5a), a large green contour near the 3'-position of the terminal 3'-trifluoromethylphenyl group of compound 20 suggests that introducing of bulky groups at this position would increase the activity. In consistent with this, compounds bearing bulky groups at this position, for example, compounds 18, 20, 21 and 23 showed high activities, whereas the ones bearing small groups at the same position, for example compounds 16 and 22, showed low activities. Large yellow contours near the 4'-and 5'-positions of the terminal 3'-trifluoromethylphenyl group suggest that introducing of bulky groups at these positions would decrease the activity. This is in agreement with the fact that compounds 44-52 with bulky groups at 4'-or 5'-positions of the terminal 3'-trifluoromethylphenyl showed decreased activities. In addition, big yellow contours near the pyrrolo [3,2-d]pyrimidine subunit suggest that steric bulkiness is unfavorable by the model. When pyrrolo [3,2-d]pyrimidine subunit is conjugated to the central phenyl group via a sulfur or nitrogen atom (X), the subunit falls into the region of the yellow contour, which suggests that a decrease in the inhibitory activity will be observed. For example, sulfur-linked derivative 12 (IC 50 = 110 nM) and amine-linked derivative 13 (IC 50 = 1400 nM) were 3-and 40-fold less active than oxygen-linked derivative 16 (IC 50 = 33 nM), respectively.
In the CoMFA electrostatic contour map, a small red contour near the 3'-position of the 3'-trifluoromethylphenyl group of compound 20 indicates that introduction of electronegative groups around this position would increase the inhibitory activity. This, together with the green contour discussed above, suggests that electronegative and bulky groups near the 3'-position of the terminal phenyl group are favored by the CoMFA model. A blue contour near the 4'-position of the 3'-trifluoromethylphenyl group indicates that introducing of electropositive groups around this position would increase the inhibitory activity. For example, compound 24 with electropositive hydrogen at the 4'-position showed higher activity than the corresponding compounds 45-48 with electronegative substituents.

CoMSIA Contour Maps
The CoMSIA steric and electrostatic contour maps of compound 20 are shown in Figure 6a,b, respectively. These contours are quite similar to those of CoMFA. Therefore, our following discussion will focus on the hydrophobic, hydrogen bond donor and acceptor fields. Figure 6c shows the hydrophobic contour maps in which yellow and gray contours indicate the regions where hydrophobic and hydrophilic groups are favored by the model, respectively. A yellow contour near the 3'-position of the terminal 3'-trifluoromethylphenyl group indicates that hydrophobic substituent at this position would increase the activity. This has been noticed in compounds 18, 20, 21 and 23 bearing hydrophobic Cl, CF 3 , Br and CH 3 , respectively. These compounds exhibited increased activities. This hydrophobic interaction may play a crucial role in improving of the binding affinity, since it is also observed in the CoMFA and CoMSIA steric contour maps. The large gray contour near the urea indicates that hydrophilic urea at this position is favorable.
The CoMSIA hydrogen bond donor and acceptor contour plots are shown in Figure 6d,e respectively. The cyan contours represent the regions where hydrogen bond-donating groups increase the activity, whereas the purple contours represent the regions where hydrogen bond-donating groups decrease the activity. Similarly, the magenta contours indicate the regions where hydrogen bond-accepting groups increase the inhibitory activity, whereas the red contours indicate the regions where hydrogen bond-accepting groups decrease the activity.
The cyan contours near the NH of urea indicate that hydrogen bond-donating groups are favored. This is well consistent with the observations that NH group in this region forms stable hydrogen bonds with the residue of Glu883 as hydrogen bond donor, and that the amide derivatives 6 and 7 having a CH 2 to replace the NH of the urea 16 led to more than 300-and 100-fold decreases in the activities, respectively. A large magenta contour located on the carbonyl oxygen of the urea moiety suggests that hydrogen bond-accepting groups are favored in this region. This is evident from the fact that thiourea derivative 8 was about 100-fold less active than the corresponding urea derivative 16, and that the carbonyl oxygen of the urea moiety was engaged in hydrogen bonding interaction with the Asp1044 in the DFG motif.

Design of New Inhibitors
As shown above, the CoMFA and CoMSIA have provided detailed insight into the key structural requirements for potent activities of the inhibitors of this class. Specifically, the urea plays a crucial role in the inhibitory activity-its replacement with thiourea or with an acetamide leads to a complete loss of the activity. Substituting of the oxygen linker with sulfur or nitrogen atoms affords less active compounds. Introducing of appropriately bulky and strongly hydrophobic groups at the 3'-position of the terminal phenyl group may significantly increase the activity. To demonstrate the practical values of these structure-activity relationships, a series of new inhibitors were designed and their pIC 50 values were predicted with the established CoMFA and CoMSIA models (Table 4).

Concluding Remarks
In this study, 3D-QSAR analyses, CoMFA and CoMSIA, have been applied to a set of recently synthesized pyrrolo [3,2-d]pyrimidine derivatives as type II KDR inhibitors. The binding mode of the template molecule 20 was clarified by Surflex-dock. The CoMFA and CoMSIA models showed statistically significant results in terms of cross-validated coefficients and conventional coefficients. Their predictive capabilities were verified by the test compounds. The derived CoMFA and CoMSIA models showed predictive cross-validated coefficients of 0.913 and 0.897, respectively, and the activities of the training and test compounds were predicted with good accuracy. Based on the obtained structure-activity relationships, a series of new inhibitors were designed to have excellent activities that were predicted with the developed CoMFA and CoMSIA models. Thus, these models may be expected to serve as a tool to guide the future rational design of pyrrolo [3,2-d]pyrimidine-based novel type II KDR inhibitors with potent activities.