Studies on [5,6]-Fused Bicyclic Scaffolds Derivatives as Potent Dual B-RafV600E/KDR Inhibitors Using Docking and 3D-QSAR Approaches

Research and development of multi-target inhibitors has attracted increasing attention as anticancer therapeutics. B-RafV600E synergistically works with vascular endothelial growth factor receptor 2 (KDR) to promote the occurrence and progression of cancers, and the development of dual-target drugs simultaneously against these two kinds of kinase may offer a better treatment advantage. In this paper, docking and three-dimensional quantitative structure activity relationship (3D-QSAR) studies were performed on a series of dual B-Raf/KDR inhibitors with a novel hinge-binding group, [5,6]-fused bicyclic scaffold. Docking studies revealed optimal binding conformations of these compounds interacting with both B-Raf and KDR. Based on these conformations, comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) 3D-QSAR models were constructed, and the best CoMFA (q2 = 0.542, r2 = 0.989 for B-Raf; q2 = 0.768, r2 = 0.991 for KDR) and CoMSIA models (q2 = 0.519, r2 = 0.992 for B-Raf; q2 = 0.849, r2 = 0.993 for KDR) were generated. Further external validations confirmed their predictability, yielding satisfactory correlation coefficients (r2pred = 0.764 (CoMFA), r2pred = 0.841 (CoMSIA) for B-Raf, r2pred = 0.912 (CoMFA), r2pred = 0.846 (CoMSIA) for KDR, respectively). Through graphical analysis and comparison on docking results and 3D-QSAR contour maps, key amino acids that affect the ligand-receptor interactions were identified and structural features influencing the activities were discussed. New potent derivatives were designed, and subjected to preliminary pharmacological evaluation. The study may offer useful references for the modification and development of novel dual B-Raf/KDR inhibitors.


Introduction
Over the past decades, the targeted therapy, which aims to develop highly selective inhibitors against a specific drug target, has become the primary paradigm in drug discovery [1][2][3]. However, many complex diseases, such as cancers, are not modulated by a single target, but rather an intricate network of signal pathways. Given the complexity of cancers, there is a general agreement that a multi-target approach might be more effective [4][5][6]. Therefore, drugs simultaneously modulating multiple targets have been extensively explored, which may help solve problems of single-target agents, such as limited efficacy, poor safety, and resistance. One of the multi-target strategies involves developing a single drug that can simultaneously act on two or multiple targets, which is more predictable on the pharmacokinetics and pharmacodynamics [7][8][9][10][11][12]. To date, more and more multi-target drugs have been successfully developed [7].
Members of the Raf family play a key role in the mitogen-activated protein kinase (MAPK) pathway. This pathway is closely related to cellular proliferation, differentiation, and survival [13]. In many cancers, oncogenic mutations of B-Raf are frequently observed and cause abnormal proliferation receptor tyrosine kinases (RTKs) without any control [14]. In addition, many solid tumors require sufficient amount of nourishment and new blood vessels are constructed for their progression in the aid of vascular endothelial growth factor (VEGF), which stimulates adjacent vascular endothelial growth factor receptor 2 (VEGFR-2 or KDR). The VEGFR-2 plays a significant role in the progression of solid tumors [15]. B-Raf with VEGFR-2 works synergistically to promote certain cancers, and to develop multi-target drugs simultaneously against these two kinases may provide a better therapeutic advantage. Thus, dual B-Raf and VEGFR-2 inhibitors could be effective therapeutic agents for cancer treatment [16]. Fortunately, the high degree of sequence and structural homology among the ATP-binding pockets of B-Raf and VEGFR-2 makes it plausible to design dual B-Raf and VEGFR-2 inhibitors [17,18]. Sorafenib is the first kinase drug that targets both Raf/Mek/Erk and VEGFR-2/PDGFR signaling cascades to block the tumor cell proliferation and inhibit the tumor angiogenesis [19]. Sorafenib is a type II inhibitor that interacts with the inactive DFG-out conformations of wild type B-Raf and oncogenic mutant B-Raf V600E as well as VEGFR-2 [20][21][22]. Studies have indicated that Type II inhibitors have certain advantages compared to Type I inhibitors, which includes the improvement of biochemical efficiency and selectivity [23]. Thus, dual inhibition of B-Raf and VEGFR-2 may provide strong antitumor efficacy, especially Type II inhibitors.
To design multi-target agents, theoretical studies on the structural features of the compounds may provide useful information. Molecular docking methodology is robust in the study of protein-ligand interaction, and quantitative structure-activity relationship (QSAR) methodology is considered to be the most effective in understanding the structure and activity relationship and designing new chemical entities [24][25][26]. The three-dimensional quantitative structure activity relationship (3D-QSAR) studies, including comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) [27,28], help to identify the interaction fields surrounding the compounds according to their corresponding activity, which may shed light on designing of novel B-Raf and KDR inhibitors and help with the issue of commonality and selectivity among B-Raf and KDR members. Therefore, combined docking and 3D-QSAR studies [29][30][31][32][33] may offer useful information for designing and obtaining more potent new dual kinase inhibitors.
Recently, molecules with a novel hinge-binding group have been reported to show dual inhibition against both B-Raf and KDR, which form two hydrogen bond interactions with hinge region Cys amino acid residues through various [5,6]-fused bicyclic scaffolds. These derivatives were proven to be the highly potent DFG-out type RAF/VEGFR-2 inhibitors. The representative compound 0JA (Ligand from the cocrystal structure of VEGFR-2 (PDB code: 3VNT) and B-Raf (PDB code: 4DBN)) showed potent inhibitory activity against B-Raf V600E (IC50 = 7.0 nM) and wild-type B-Raf (12 nM), C-Raf (1.5 nM) as well as VEGFR-2 (2.2 nM) ( Figure 1) [34][35][36][37]. In this paper, novel series of [5,6]-fused bicyclic derivatives as dual B-Raf and KDR inhibitors were carried out in a theoretical study using the docking and 3D-QSAR methods. The binding modes and conformations of these dual-target inhibitors with both B-Raf and KDR were determined by docking, and 3D-QSAR models were established based on docked conformations. The model is expected to facilitate deeper understanding of the structure activity relationship of these compounds, and provide useful information for rational design of novel and potent dual B-Raf and KDR inhibitors.

Structure Alignment Analyses
The accuracy and reliability of the 3D-QSAR model is dependent on two important factors: the biological conformation selection and the structural alignment rule. In this study, the 3D-QSAR model is constructed by the docking-based conformation alignment. To ensure molecular docking could recapture feasible conformation, Glide-SP docking algorithm was first validated by re-docking of the crystal ligand. The best docking pose of compound 0JA (good docking score was observed) was only minimally differed from the co-crystal ligand conformation with a root mean-square deviation (RMSD) of 0.4967 and 0.2836 Å (Table 1, Figure 2a,b), suggesting high accuracy and reliability of Glide in reproducing the experimental binding mode. Ligand 0JA showed significant binding interactions with the key residue of the target protein in the docking study. 0JA acted on the DFG-out conformation of B-Raf (PDB code: 4DBN). The [1,3]-thiazolo-[5,4-b]-pyridine-2-amine moiety formed the hydrogen bonds with the carbonyl and NH of Cys532 in the hinge region. The core was positioned in a hydrophobic cleft consisting of Ile463, Val471, Ala481, Trp531, Phe583 and Phe595. The central phenyl occupied the hydrophobic pocket adjacent to the ATP binding site. Carbonyl and NH of the biphenyl amide moiety formed hydrogen bonds with Asp594 and Glu501, respectively. Terminal phenyl extended into a back hydrophobic pocket generated by the Phe595 movement and lined with Val504, Leu505, Ile513, Ile572 and His574.
Docking of ligand 0JA in the X-ray co-crystal structure of KDR (PDB code: 3VNT) was depicted in Figure 2b. Similar to B-Raf, 0JA acted on the DFG-out conformation of KDR. The [1,3]-thiazolo-[5,4-b]-pyridine core was positioned in a hydrophobic cycle lined with Leu840, Ala866, Phe918 and Leu1035. Similar hydrogen bond interactions were observed with the carbonyl and NH of Cys919 in the hinge region. The central phenyl also occupied the hydrophobic pocket adjacent to the ATP binding site. Carbonyl and NH of amide formed hydrogen bonds with Asp1046 and Glu885, respectively. The terminal phenyl extended into a back hydrophobic pocket created through the Phe1047 flip and lined with Ile868, Leu889, Ile892, Leu1019 and His1026. To obtain accurate binding conformations, all compounds ( Table 2) were docked into the active site of B-Raf and KDR kinases by using the Glide-SP docking algorithm. Figure 3a,b show the 3D models of molecular alignments of 40 compounds at the binding site of B-Raf and KDR. Docking conformations of all studied compounds were similar to co-crystal ligand 0JA conformation, providing an excellent alignment. Taken together, results from the Glide docking based alignments were suggested to be reasonable.

3D-QSAR Analyses
The CoMFA and CoMSIA models for B-Raf and KDR were performed based on docking conformations alignment. The statistical parameters of PLS analyses and external validation results were listed in Table 3.

3D-QSAR Results for B-Raf Kinase
Through PLS statistical analysis, Glide docking based CoMFA model for B-Raf kinase showed a reliable cross-validated correlation coefficient q 2 of 0.542 and optimal number of components of 6, indicating a good internal prediction of model. The CoMFA model also exhibited a high correlation coefficient r 2 of 0.989 with relatively lower SEE of 0.075 and relatively higher F value of 376.623 in the final non-cross-validated model. The contribution of steric fields and electrostatic fields was 45.7% and 54.3%, respectively, which indicated that both steric field and electrostatic field had equally important influences. The above values suggested a good statistical correlation and a good internal predictive ability of the CoMFA model as shown in Figure 4a. In order to further validate the models' predictive ability, activities of test set compounds not included in the construction of the 3D-QSAR models were predicted (shown in Table 4). Both CoMFA and CoMSIA exhibited satisfactory results in term of predictive correlation coefficient r 2 pred of 0.764 and 0.841, respectively. The favorable r 2 pred values indicated good external predictive ability and were able to predict inhibitory activities of new inhibitors. The experimental, predicted and residuals of test set molecules were shown in Table 4 and the correlations between the actual and predicted pIC50 values of test set compounds are presented in Figure 4b,d. These statistical data and graphs indicated that the predicted values were in agreement with the experimental values in the allowable error range.

3D-QSAR Results for KDR Kinase
The 3D-QSAR analysis of KDR was performed on the same training set for B-Raf. The CoMFA model yielded q 2 = 0.768 and r 2 = 0.991, respectively. It was obvious that CoMSIA models also presented well in terms of higher q 2 = 0.849 and r 2 = 0.993 values. Other statistical parameters of PLS analyses were listed in Table 3.
The actual activities (pIC50), the predicted activities and the corresponding residual values of the training set for the CoMFA and optimal CoMSIA models are listed in Table 4. Graphical representations of actual vs. predicted activities of training set are shown in Figure 5a,c. The CoMFA and optimal CoMSIA models possessed high q 2 and r 2 values, indicating that models had good internal predictive ability. To validate the external predictability of the models, the r 2 pred values were calculated for test set. As shown in Tables 3 and 4, the r 2 pred values of CoMFA and optimal CoMSIA models spanned from 0.912 to 0.846, revealing that models had good external predictive ability and could be used to predict the biological activities of novel compounds. The plots of actual vs. predicted activities of test set were shown in Figure 5b,d, showing that the predicted activities were in good agreement with the actual data.

Contour Maps
To visualize the results of the CoMFA and CoMSIA models more directly, the 3D coefficient contour maps of CoMFA (steric and electrostatic fields) and CoMSIA (steric, electrostatic, hydrophobic, and hydrogen bond acceptor fields) were generated (Figures 6-8 and 10), respectively. To facilitate the analysis, ligand 0JA was selected as the reference in the 3D coefficient contour maps. The results of the CoMFA and CoMSIA models were graphically interpreted by the field contribution maps.

CoMFA Contour Maps
The contour maps of CoMFA (steric and electrostatic fields) are shown in Figure 6. In the contour map of steric field, green contour showed sterically favored region while yellow region indicated the area where bulky groups may cause decline in the inhibition activity of compounds. In the contour map of electrostatic field, red contour showed the region where electronegative group was favorable to increase the inhibitory activity while opposite was for blue contours. In the contour map of steric field (Figure 6a), a large green contour was observed around the cyanocyclopropyl group of 2-chloro-3-(1-cyanocyclopropyl)benzene ring (ring-C), suggesting the bulky substituent was favored at this region such as methoxyl, trifluoromethoxyl, and cyanocyclopropyl. This bulky hydrophobic interaction may had an important role in enhancing the cellular activity against B-Raf, which was illustrated by the experimental fact that compound 1-10 exhibited higher activity than corresponding compound 1-9. The small green contours were found at the N atom position on the pyridine ring of [1,3]-thiazolo-[5,4-b]-pyridine scaffold (ring-A), around which the moderate-sized substituent was favored such as cyano. The cyclopropyl group near a hinge was surrounded by two yellow contours, which suggested bulky groups at this position would be unfavorable for binding the B-Raf protein. Because the position existed in a narrow space formed by the indole ring of Trp531 and Gly534, and that the optimal size for the narrow space would be smaller than the pyranyl group. This may explain why compound 1-1 with a bulky group showed significantly decreased activities than other compounds with a minor substituent near a hinge. Therefore, replacement of the pyranyl group in compound 1-1 with a methyl group in compound 1-2 showed significantly increased B-Raf inhibitory activity (IC50 = 4.1 nM). This was also well illustrated by the order of activity for these compounds: 1-4 > 1-5 > 1-6.
Electrostatic contour map was shown in Figure 6b. In general, red contours are close to heteroatoms, whose partial atomic charges were negative such as nitrogen and oxygen. The moderate red contours were found close to the N atom position on the pyridine ring of [1,3]-thiazolo-[5,4-b]-pyridine scaffold, indicating negative charged groups was favored. This trend can be reflected by the activities of compounds 5-2, 5-3, 6-5, 6-6 and 6-7 which all had the cyano group. These compounds were more active than corresponding compound 2-5. The 2-chloro-3-(1-cyanocyclopropyl)benzene ring was surrounded by three small red contours, suggesting that electronegative potential was preferred at these positions, such as F, Cl, OC(CH3)2CN, C(CH3)2OH, and CF3. The red contours above the phenyl of the benzene amide group demonstrated that the negative-favorable substituents induced the positive carbon of phenyl ring. With respect to the favorable positive potential, the blue region around the phenyl urea and benzene amide group showed key positive-favorable property, which reflected an important H-bond donor feature of NH group, consistent with the docking study. The blue contours were distributed around the central benzene ring (ring B) group, suggesting that positive charged groups increase the activity. This was consistent with the increase in the potency of compound 3-3 as compared to compound 3-2, due to the substitution of hydrogen atom with fluorine atom in the central phenyl group.

CoMSIA Contour Maps
The contour maps of CoMSIA steric and electrostatic fields are shown in Figure 7a,b, which were found to be quite similar to the contour maps derived from the above CoMFA model except several slight differences, and thus were not discussed here. The following analyses focus on the CoMSIA hydrophobic and hydrogen bond acceptor contour maps. In hydrophobic contour map of CoMSIA, yellow and gray contours indicated hydrophobic and hydrophilic groups were favored, respectively. Figure 7c showed hydrophobic contour map of CoMSIA, gray contour was located around the linker between the ring-A and ring-B, which suggested hydrophilic contributions had effect on the activity. This may explain why compounds 3-6, 3-7 and 3-8 with an N-Me substituent were less active than compounds only with NH or O substituent here. The gray contour near the cyclopropyl position of the terminal phenyl group of compound 0JA indicated that hydrophilic substitution at this position was favorable. This may explain why compound 4-2 possessing a more hydrophobic substituent (-CH2C(CH3)2CN) showed significantly decreased activities than other compounds (4-1, 4-3) with a less hydrophobic substituent (e.g., -CH2CH2CN, -OC(CH3)2CN).
As indicated in Figure 7c, the white region was close to the polar CONH group and R-position of the ring-B, which interacted with hydrophilic amino acid residues, such as Glu501, Thr508, and Asp594. Hydrophobic substituents at R-position of the ring-B such as -F, -Cl were not beneficial to the improvement of inhibitory activity. For example, compounds 3-2, 6-2 and 6-6 with hydrogen at this position were more active than compounds 3-4, 6-4 and 6-8 with a Cl or F group at the same position. The yellow regions were situated around the ring-C. It was shown that the position points to the hydrophobic pocket of the binding site, interacting with the Ile513, His574, and Ile572 residues. It was shown that N1 of the ring-A is surrounded by the yellow region. This can be explained by the fact that compounds with more hydrophilic substitution (e.g., COOH) had decreased activities than compounds 5-3 and 5-4 with CO2Me and CH2OH as substituents.
The magenta contours indicate hydrogen bond-accepting groups increase the inhibitory activity, whereas the red contours indicate hydrogen bond-accepting groups decrease the activity. A magenta contour located on the carbonyl oxygen of the urea or NHCO moiety suggested that hydrogen bond-accepting groups were favored. This is obvious from the fact that the carbonyl oxygen participated in hydrogen-bonding interaction with the Asp594. The docking study also showed that the carbonyl forms a hydrogen-bond with Asp594. Another magenta contour was found near the N atom position on the pyridine ring of [1,3]-thiazolo- [5,4-b]-pyridine scaffold, which indicated that H-bond acceptor groups were favored there. Thus, it was not surprising that compounds 5-2 and 5-3 with -CN and -CO2Me as substituents had higher activities than corresponding compound 2-5. A red contour located on the S atom position of the thiazole pyridine ring of [1,3]-thiazolo- [5,4-b]-pyridine scaffold, indicating that H-bond acceptor groups were unfavorable there. The remaining red contour was around the H atom of the pyridine ring indicating that the hydrogen bond acceptor substituent was inessential for the increased inhibitory activity. Hence, compounds 4-1, 4-2, 4-3 did not show improved potencies than corresponding unsubstituted derivatives. Moreover, there were red acceptor-unfavorable contours over the 2-chloro-3-(1-cyanocyclopropyl)benzene ring of the template molecule, suggesting that the hydrogen bond acceptor group was not tolerated at this position. It was also in line with the above result that there was a hydrogen bond donor site near this region. For example, compound 1-11 with a -OH as substituent showed higher activity than corresponding compounds 1-10 and 1-12.

CoMFA Contour Maps
The contour maps of the CoMFA steric and electrostatic fields for KDR are shown in the Figure 8a,b. They were similar to the contour maps of the CoMFA steric and electrostatic fields of B-Raf (Figure 6a,b). However, substitution of the N-acyl group (R1) had negligible impact on the KDR inhibition, with compounds 1-1, 1-2, 1-3, 1-4, 1-5, and 1-6 showing nano-molar order inhibitory activity against KDR. For example, in the steric contour map there was a small yellow contour beside N position at the C-7 position on the pyridine ring of [1,3]-thiazolo- [5,4-b]-pyridine scaffold, indicating that introduction of bulky groups here would decrease the KDR inhibitory activity. This was the opposite with steric contours of B-Raf. This can explain the favorable activity of compounds 2-5, and the unfavorable activity of compounds 5-1, 5-2, and 5-3. The reasons accounting for the aforementioned phenomenon may be as follows: there was a prominent difference in the Phe residue conformations of the DFG motif for B-Raf and KDR (Figure 9). The benzene ring of Phe595 in B-Raf (green) was located under the thiazolo-[5,4-b]-pyridine scaffold, whereas that of KDR (yellow) was located in front of the N atom position of pyridine. Consequently, the introduction of a bulky substituent could reduce the KDR inhibitory activity by steric repulsion. Studies found the C-7 position should be very important in modulating the ratio of KDR IC50 to B-Raf IC50. On the basis of the notable difference, we could design inhibitors to increase B-Raf selectivity over KDR.  In the electrostatic field of CoMFA model of KDR, there is a blue contour located on the S atom position of [1,3]-thiazolo-[5,4-b]-pyridine scaffold, indicating that positive charged groups increased the activity there. With respect to electrostatic contours of B-Raf, no blue contour located on the S position of thiazole, whereas there was blue contour around the amide linker group for B-Raf, which showed a key positive/favorable property and reflected an important H-bond donor feature of NH. According to the compound data, the NHCO substructure connected to ring-B and ring-C could apparently form two significant hydrogen bonds with Glu501 and Asp594. In order to study the structure-activity relationship, the linker between ring-B and ring-C was modified to explore appropriate linkers. An important modification for the linker was insertion of X (NH or CH2). The insertion of amine (X = NH) of ureides 6-1, 6-2, 6-3, 6-4 could form an additional hydrogen bond with Glu501, and the insertion of methylene (X = CH2) of the acetamides 6-5, 6-6, 6-7, and 6-8 was considered flexible enough so that ring-C could adjust the conformation to be accommodated into the hydrophobic back pocket of B-Raf. However, the impact of these chemical modifications of the linker for KDR was negligible.

CoMSIA Contour Maps
As shown in Figure 10, the CoMSIA steric and electrostatic contour maps (Figure 10a,b) were almost identical with CoMFA contour maps. Thus, the following discussion focuses on the hydrophobic and hydrogen bond acceptor fields. The contour maps of the hydrophobic and hydrogen bond acceptor fields of the CoMSIA model can be seen clearly in Figure 8c,d. They were also similar to the CoMSIA hydrophobic and hydrogen bond acceptor contour maps of B-Raf (Figure 7c,d). The hydrophobic fields are presented in Figure 10c, yellow and white contours indicate hydrophobic and hydrophilic groups were favored, respectively. Compared with the hydrophobic contour of B-Raf, a yellow polyhedron that covered cyclopropane group near the hinge region suggested that a moderate-sized hydrophobic group could help increase the biological activity. This could be validated by the fact that the activity of compounds 1-2, 1-3, and 1-4 had better biological activity than compound 1-1. While no white contours are close to the central phenyl ring, indicating that hydrophilic substitution at this position was not necessary. Figure 10d showed contour maps of CoMSIA hydrogen bond acceptor field. Compared to the hydrogen bond acceptor field of B-Raf, the magenta contour near o-Cl of the terminal phenyl ring-C revealed that hydrogen bond acceptor groups might be beneficial for the potency. The results were supported by the higher potency of compounds 1-8, 6-1, and 6-5 with CF3 at this position than corresponding compound 1-4 as well as compounds 1-7, 6-2, 6-6 and 6-7. Without magenta contour regions around urea or amide groups, the CO of urea or amide groups could provide one additional H-acceptor, which leads to a lower inhibitory activity.

Comparison of the Protein Structures
To compare these two proteins, the sequences were aligned and the crystal structures were superposed for B-Raf (PDB code: 4DBN) and KDR (PDB code: 3VNT). The comparison of the sequences and proteins are shown in Figures 9 and 11.
The sequence identity and similarity were 25.0% and 46.2%, respectively. Identity and similarity of B-Raf and KDR chains were not very high, but B-Raf and KDR had a lot of common key points. As displayed in the alignment of the sequences, the same sequence ID for B-Raf and KDR was marked in blue shadows. The key amino acids in the active sites of B-Raf (PDB code: 4DBN) and KDR (PDB code: 3VNT) were marked in red rectangles. It can be observed that the docking active sites for B-Raf and KDR may have almost the same active sequence domains. Figure 9 depicted the structural superposition of B-Raf (PDB code: 4DBN) and KDR (PDB code: 3VNT), where the binding conformation and interaction were illustrated, with the root mean square deviation (RMSD) of 1.964 Å. The binding conformation and interaction of ligand 0JA with the two kinase receptors B-Raf and KDR were strikingly similar. Almost all key amino acids in the active sites (such as Glu501, Cys532, and Asp594 for B-Raf, and Glu885, Cys919, and Asp1046 for KDR) interacting with the [5,6]-fused bicyclic derivatives were well aligned for B-Raf and KDR, which may explain why these compounds are dual B-Raf /KDR inhibitors. Figure 11. Sequence alignment of 4DBN (B-Raf) and 3VNT (KDR). Deep Blue color regions represent that the amino acid residues in the individual column were identical in the sequence alignment. The key amino acids in the binding sites were marked in red rectangles.

Comparison of the 3D-QSAR Results
From the statistical results listed in Table 3, the CoMFA and CoMSIA models for B-Raf and KDR kinases exhibited excellent predictive powers. The statistical results of the CoMFA and CoMSIA models for B-Raf and KDR were very similar. For CoMFA models, both the steric and electrostatic fields both had equally important influences on the ligand-receptor interactions, with approximately 50% contributions to the steric and electrostatic fields, respectively. Furthermore, optimal CoMSIA models of B-Raf and KDR were constructed by using the same combination of steric, electrostatic, hydrophobic and hydrogen acceptor fields, indicating that the ligand binding to both B-Raf and KDR may be the similar. In addition, the electrostatic and hydrophobic fields were observed to play crucial roles for the interaction, with contributions of 61.2% (electrostatic: 35.6%, hydrophobic: 25.6%) and 60.5% (electrostatic: 33.4%, hydrophobic: 27.1%) in the CoMSIA models for B-Raf and KDR, respectively.

Comparison of the 3D-QSAR Contour Maps
Many similarities in sequences and protein structures of B-Raf and KDR may lead to very similar structural features of the B-Raf and KDR inhibitors. Through analysis and comparison of the B-Raf and KDR 3D-QSAR contour maps (Figures 6-8 and 10), very similar structural requirements and small differences could be found for potent ligands of B-Raf and KDR kinases ( Figure 12). From the discussion above, we could conclude common structural requirements for more potent dual inhibitors against B-Raf and KDR as follows: (i) Ring-A and linker between ring-B and ring-C form key H-bonds with kinase.
(ii) Ring-A, -B, and -C are hydrophobic to interact with the Adenine site, DFG motif toward the allosteric site; and the cyanocyclopropyl-position should be hydrophilic-favorable.
However, some small differences between the contour maps of B-Raf and KDR could also be easily found. These differences may contribute to better understanding of their specific structural requirements for kinase selectivity.

Designing Potent Derivatives for Further Verification
Based on the results of established receptor-ligand interactions and 3D-QSAR above, and for further verification and design consideration, some novel compounds based on the 1H-benzo[d]imidazole moiety [38] with potential inhibitory activity were designed as B-Raf V600E /KDR dual inhibitors. The chemical structures and predicted activity values by the CoMFA and optimal CoMSIA models are shown in Table 5. The designed molecules were synthesized and the preliminary pharmacological evaluation would be used to verify the biological activity and our design strategy. The experimental inhibitory activities of designed new compounds were also listed in Table 5. The preliminary biological tests showed new compounds have dual B-Raf V600E /KDR inhibitory activity. Among them, compound D2 had the best inhibitory effect with the IC50 of 420 and 210 nM for B-Raf V600E and KDR kinases, respectively. The results showed these theoretical results would provide the useful references and guidelines for optimizing current molecules and designing novel inhibitors with new scaffolds.

Data Set
For 3D-QSAR analyses, 40 compounds from the literature [34][35][36][37] were employed, which covered a range of more than three magnitudes (1.3-3700 nM for B-Raf, 0.31-730 nM for KDR) for their activity. The IC50 values of all compounds for B-Raf and KDR inhibition were converted to the −log IC50 (pIC50) as dependent variables in the 3D-QSAR analysis. The 3D structures for all 40 compounds were built and minimized, and then Gasteiger-Hückel charge was calculated. All works were all done in SYBYL 6.9 (Tripos, Inc. Normally, the whole dataset for 3D-QSAR study would be divided into the training set and test set randomly. However, here, we performed an alternative way by a data mining methods of clustering [40,41], which may bring us more reliable results. In the clustering process, the dataset was clustered into 7 classes (shown in Table 2) according to the molecular descriptors and structural difference. The compounds were chosen from the above 7 categories with a ratio of about 4/1 of the training set and test set, according to biological activity and structural diversity. So 33 compounds as training set and the other 7 ones (asterisked in Table 2) as test set were selected to construct and validate the model, respectively.

Conformational Alignment
To construct reliable 3D-QSAR models, structure alignment is considered the most important and critical step [42]. The most common method is to align compounds onto a common scaffold. For our 3D-QSAR study, molecular docking based alignments were shown as the most effective alignment by fully considering the binding mode of compounds. Molecular docking can be used for alignment because it can dock compounds into the target binding site with similar modes. Meanwhile, the docking studies would also be conducive to the view of ligand-receptor interaction, assisting in further understanding the structure-activity relationship.
Glide 5.5 was selected to generate biological conformation as it has shown good reproducibility and accuracy in molecular docking and scoring [43,44]. The crystal structures of B-Raf and KDR (PDB code: 4DBN and 3VNT) were prepared with the Protein Preparation Wizard workflow in Schrödinger suite 2009 software (Schrödinger, LLC, New York, NY, USA, 2009). The grid files were generated, which were defined by the co-crystal ligand itself. All remaining parameters were kept as default. The compounds were flexibly docked into the binding site with standard precision (SP) docking mode. The optimal conformation of each compound was retained based on the Glide score and interaction modes.

Generation of 3D-QSAR Model
The CoMFA and CoMSIA analyses were carried out using the QSAR module in SYBYL6. 9. In CoMFA study [27], the standard steric and electrostatic fields were calculated in the cubic lattice by a sp 3 hybridized carbon with a +1.0 charge and grid spacing of 2.0 Å. The cut off value was set to 30 kcal/mol and the column filtering was set to be 2.0 kcal/mol for both fields. In CoMSIA study [28], steric, electrostatic, hydrophobic, and hydrogen bond donor and acceptor fields were calculated at each lattice by a sp 3 carbon probe atom with a charge of +1.0 and a grid of 2 Å. The attenuation factor was set to the default value of 0.3.

Partial Least Squares (PLS) Analysis and Validation
The Partial Least-Squares (PLS) [45,46] algorithm was employed for the structure-activity relationship. Biological activities (pIC50 values) as dependent variables and the field descriptors as independent variables were used to construct 3D-QSAR models. In cross validation, leave-one-out (LOO) algorithm was adopted to yield the optimal number of components (ONC), and cross-validation coefficient (q 2 ), calculated with Equation (1): where Ypred, Yexp and Ymean are the predicted activity values, experimental activity values and mean activity values of training set compounds, respectively. In the non-cross validation, the correlation coefficient (r 2 ) was calculated using the obtained ONC. Additionally, the statistical significance of the models was also described by the SEE and F probability