Discovery of vascular endothelial growth factor receptor tyrosine kinase inhibitors by quantitative structure–activity relationships, molecular dynamics simulation and free energy calculation

Abstract

Vascular endothelial growth factor (VEGF), along with its receptor tyrosine kinases VEGFR-2 or kinase insert domain receptor (KDR), are targets for development of novel anticancer agents. Accurately predicting the structural characteristics of the target and chemical features of ligands can greatly reduce the cost and shorten the cycle of designing selective KDR inhibitors with desired activity. In this study, a docking strategy and three dimensional holographic vector of atomic interaction field (3D-HoVAIF) were applied in QSAR analysis of KDR inhibitors. The optimal model was constructed by using stepwise regression combined with partial least squares regression (SMR-PLS). Integrating the results of QSAR analysis, ADMET, pharmacophore modeling and a reverse screening strategy, eight derivatives were identified as potential KDR inhibitors. Then molecular dynamics (MD) simulations and free energy calculations were employed to explore the detailed binding process, so as to compare the potential binding modes of inhibitors with different activities. By analyzing the key residues in the binding site, it was found that different KDR–ligand complexes had similar binding modes. The predicted binding affinities were highly correlated with the experimental biological activity. Free energy analysis indicated that van der Waals interactions provided the major driving force for the binding process. Furthermore, key residues, such as Leu840, Val848, Ala866, Lys868, Leu889, Val899, Thr916, Phe918, Cys919, Leu1035, Cys1045, Asp1046, and Phe1047 played a vital role in forming hydrogen bonds, salt bridges, and hydrophobic interactions with the conformation of KDR. The above results will help design more efficient KDR inhibitors.

---

Introduction

Angiogenesis, the formation of new blood vessels from existing vasculature, is an important pathway for normal embryogenesis, growth, and tissue repair.^1–4 It is not only essential to many normal physiological processes, but also contributes to pathological disorders and, in particular, has implications in the growth and metastasis of solid tumors.^5–9 Neovascularization resulted from a complex set of events involved a number of diverse stimuli, however, signaling induced by vascular endothelial growth factor (VEGF) was considered rate-limiting. VEGF drove the angiogenic cascade by promoting endothelial cell activation and proliferation,^2,8–10 and enhanced the migration and invasion of endothelial cells.¹¹ Tumors can up-regulate VEGF expression or secretion via a quantity of pathways, including oncogenes, cytokines, or cell growth factors activation, loss of tumor suppressor function, and change in oxygen levels.¹² VEGF bound to vascular cell surface expressed VEGFR-2 (KDR), a receptor tyrosine kinase (RTK), which upon dimerization/activation by VEGF underwent autophosphorylation and initiated downstream signaling, ultimately leaded to tumor vascularization, survival and migration.^13,14 Up to now, several successful strategies for the inhibition of angiogenesis, in particular, used small molecule inhibitors to interdict the KDR (VEGFR-2) signaling had been proved to effectively inhibit tumor progression and dissemination in preclinical and clinical studies.^15,16 Therefore, the development of KDR inhibitors had become an attractive strategy in the treatment of cancers.¹⁷

In recent years, many studies had indicated that computational approaches can timely provide very useful information and insights for drug development.^18–20 Such as predicting functions of proteins,²¹ structural bioinformatics,²² molecular docking, and predicting drug–target interaction^23–25 were widely welcome by science community. Quantitative structure–activity relationship (QSAR) method, as a powerful quantitative methodology, were proved to be one of the most promising and successful methods for rapidly predicting the biological activity and/or toxicity of chemicals.²² Since the first topological approach was proposed by Wiener in 1948,^26,27 there were many QSAR methods (e.g. ST-scale,²⁸ HESH,²⁹ G-scale,³⁰ and T³¹ etc.) achieved predictable performance on drug development as well as the relevant research areas. In the process of drug discovery and design, researchers often interested in quicksort the ligands of potential pharmaceutical target according to the binding-free energies (or affinity) toward a given protein. Accurately predict the interaction and the binding-free energy (or binding affinity) of protein and ligand was crucial in many fields such as structure-based drug design and lead optimization.^25,32–36

To seek and design new compounds that maintained or enhanced the excellent therapeutic properties, while decreased side effects associated with treatments by using the computer aided drug design (CADD) approaches. Molecular modeling studies of KDR inhibitors were performed by combining molecular docking, 3D molecular structural characterization (MSC), pharmacophore modeling, reverse docking, molecular dynamics (MD) and free energy calculations. These methods were employed to investigate the binding modes between KDR with a series of inhibitors and also to find the key structural features affecting the inhibiting activities. All these information can provide some insights into the further structural modifications and development of new potent inhibitors with high inhibitory activities, which can selectively inhibit KDR.

Methodology

Atomic interaction potential energies

Electrostatic interaction. As an important non-bonded interaction obeying Coulomb's law, electrostatic interactions among all atoms included in a molecule could be given out by eqn (1), and then accumulated them together into each of the 55 interaction items according to their atom-pair attributes. Where r_ij denoted interatomic Euclid distance, with the unit of m; e was the elementary charge (1.6021892 × 10⁻¹⁹ C); ε₀ represented dielectric constant 8.85418782 × 10⁻¹² C² J⁻¹ m⁻¹ in vacuum; Z was the amounts of Mülliken partial charges for atoms; m and n were atomic types.³⁷

Steric interaction. It described by Lennard–Jones formula eqn (2). Amongst, ε_ij = (ε_ii × ε_jj)^1/2 was potential well of atomic pairs, with its value taken from references.

R_ij³ = (C_hR_ii³ + C_hR_jj³)/2, was van der Waals' radius for modified atom-pair, with corrected factor C_h of 1.00 in case of sp³ hybridization, 0.95 sp² hybridization and 0.90 sp hybridization.³⁸

Hydrophobic interaction. In 3D-HoVAIF, hydrophobic interaction which was very important for drug molecules to bind to organisms, indicated information on systematic entropic changes, was indicated by eqn (3). It was simply defined in Hint proposed by Kellogg et al. In that, S was the solvent accessible surface area (SASA) for atoms,³⁹ indicated information on surface area when water-molecule probe roiling its sphere at the atomic surface; a was atomic hydrophobic constant, with the value taken from ref. 40 T was discriminant function, denoted entropy changing orientation in the case of different interatomic interactions.

Data source

The information of 82 compounds (ESI, Table S1†) were obtained from previous publication.^41–43 A naphthamide inhibitor (compound 10) in the X-ray crystal structure was used as the basic skeleton to sketch chemical structures. The molecular modeling and calculations were performed by SYBYL2.0 molecular modeling package. Energy minimization process was performed using the MMFF94 force field with a distance-dependent dielectric and Powell gradient algorithm with a convergence criterion of 0.01 kcal mol⁻¹. Partial atomic charges were loaded using MMFF94 method. Inhibitory activity (IC₅₀, nM) were converted to reciprocal logarithmic values (pIC₅₀ = −log IC₅₀) which had been used for the QSAR analysis. The whole dataset was divided randomly into a training dataset (63 compounds) which was used to construct a predictive model validated by the leave-one-out (LOO) validation, and a test dataset (19 compounds) which was used to validate the reliability of the model (Fig. 1).

Fig. 1 Distribution of inhibiting activities for training set and test set compounds in the QSAR study.

Molecular docking

The crystal structure of VEGFR2 (PDB entry code: 3B8Q) was recovered from RCSB Protein Data Bank (PDB http://www.rcsb.org/pdb). Missing loop regions that were not observed in the crystal structure were modeled and refined by Accelrys Discovery studio 3.0. In order to determine the probable binding conformations of inhibitors, the Surflex-Dock program interfaced with SYBYL2.0 was applied to study molecular docking.⁴⁴ Surflex-Dock used an empirically derived scoring function that was based on the X-ray structures and binding affinities of protein–ligand complexes. A patented search engine to dock ligands into the protein's binding site.⁴⁵ Protomol, a computational representation of the intended binding site to which putative ligands are aligned, was used to guide molecular docking.⁴⁴ The protein structure was utilized in subsequent docking experiments without energy minimization. The co-crystal ligand and water molecules from crystal structure were removed. Polar hydrogen atoms were added. Then ligand docking was employed. Other parameters were established by default in software.

Structure characterization

Common atoms in organic molecules, including H, C, N, P, O, S, F, Cl, Br and I, which were mainly dispersed in 5 families in periodic table of elements respectively (i.e. group IA, IVA, VA, VIA, and VIIA). According to the hybridization state, atoms were subsequently subdivided into 10 subtypes, which were regarded as the key to present diverse chemical properties. Thus, there were 55 atomic interactions in a molecule (Table S2†). In this paper, three kinds of potential energies, electrostatic, steric, and hydrophobic interactions, which were deemed to significantly contribute to bioactivities, were employed in 3D-HoVAIF descriptors.^46–49 Here, an organic molecule will generate 165 (3 × 55) variables.

Descriptor selection and PLS modeling

In a QSAR data set, not all the structural variables were relevant to biological activity. Those redundant variables should be deleted from the model in order to promote the robustness and predictive capability. Here, variables selection method was carried out by a classic statistical method “SMR-PLS”. SMR, less time-consuming and easier to implement, which ensure the smallest possible set of predictor variables included in the model was used to optimize variables. Partial least square (PLS) regression was a widely used modeling method, which had advantage of effectively overcoming harmful effects in modeling due to multicollinearity. It was especially suits for regression when the number of observation was less than the number of variables. In addition, PLS had the desirable property that the precision of model parameters was improved with the increasing number of relevant variables and observations.^31,50 It combined basic functions of regression model, principal component analysis (PCA) and canonical correlation analysis.²⁹ An excellent model should have both favorable estimated abilities for any internal samples and outstanding predictive abilities for any external samples.⁴⁹ SMR-PLS was completed by SPSS10.0 and SIMCA-P 10.0, respectively.

ADME and toxicity predictions

Computational approaches were being used nowadays in drug discovery program to assess the ADME and toxicity properties at the early stages of the drug discovery. In order to select potential drug candidates, ADME were carried out by employing the ADME descriptors algorithm of Accelrys Discovery studio 3.0, which stands for absorption, distribution, metabolism and excretion studies. Using this algorithm can predict various pharmacokinetic parameters, such as hepatotoxicity levels, aqueous solubility, blood–brain-barrier penetration (BBB), plasma protein binding (PPB) and human intestinal absorption. These were significant descriptors of drug likeness. Toxicity profiling of designed compounds were performed by toxicity prediction extensible protocol.^51,52

Pharmacophore modelling and reverse screening

In order to design new compounds that maintaining or enhancing the excellent therapeutic properties, while decreasing side effects associated with treatments. A pharmacophore based on an automated computational alignment technique were developed to explore the structural features of inhibitors. The reverse molecular docking was used to explore the underlying toxicity. GALAHAD (Genetic Algorithm with Linear Assignment of Hypermolecular Alignment of Datasets),⁵³ a module in the Sybyl 2.0 modeling environment, aligned a set of compounds that bound to a common target site by generating pharmacophore hypothesis.⁴⁰ In this process, new designed compounds and structurally representative compounds with high activities were separately collected to develop pharmacophore models. GALAHAD was run for 100 generations with a population size of 80 and a tournament pool size of 250. Default parameters were used for other settings. When models for all compounds with contribution to the consensus feature were considered. Based on Pareto ranking, the model with low strain energy (SE), steric overlap (SO) and pharmacophoric similarity (PhS) values was regarded as the best.⁵¹ The reverse screening strategy was carried out with the aid of in silico platform PharmMapper Server⁵⁴ that can find the potential target candidates for the query compounds.⁵⁵ Eight designed compounds as an input were uploaded to this tool to find the best pharmacophore mapping poses. The option of “Energy Minimization” was set to “Yes”, and “human protein targets only” was selected as the targets set. Genetic algorithm was used to optimize the pharmacophore mapping. Other options were the default.

Molecular dynamics simulations

The MD simulation was performed in AMBER 14.0 software package with Amber ff99SB force field.^56,57 Docked structures of KDR with compound 10, the most active compound 12, 14, 30, 64, the least active compound 47 and the new designed compound N3, N12, N13, N15, N16, N17, N18, N22 were used as the initial structures for MD calculations. Force field parameters for ligands were defined by general Amber force field (GAFF) using the Antechamber program.⁵⁸ The topology and parameters for KDR–ligand complexes were constructed by Tleap program.^56,57 Then all hydrogen atoms of the receptor were added by using the leap module, considered ionizable residues set in their default protonation states at a neutral pH value.⁵⁹ Chloridion were added to maintain charge neutrality of system. Each complex was solvated in a 10.0 Å truncated octahedron periodic box with TIP3P water molecules. Particle-Mesh Ewald (PME) algorithm⁵⁶ was used for the long-range coulombic interactions, with a cutoff of 8.0 Å. SHAKE algorithm was employed to astrict all atoms covalently bonded involving hydrogen atoms. Periodic boundary conditions were used to avert edge effects in all calculations.⁶⁰ Before MD production, solvent molecules, counterions, the lateral chain of protein and the entire model system were successively minimized by 2000 steps of steepest-descent, and 8000 steps of conjugated gradient. To obtain the stable system, a 200 ps (picoseconds) of equilibration was applied. Each system was gradually heated from 0 to 300 K over a period of 200 ps with a weak harmonic potential solutes restrained in NVT ensemble. Finally, a production run for 100 ns (nanoseconds) was performed using NPT ensemble at 300 K and 1.0 atm pressure with a time step of 2 fs (femtoseconds). Coordinate trajectories were recorded every 10 ps throughout all equilibration and production runs.

Binding free energy calculations

Free energy calculations played an important role in exploring the dynamical interaction between ligands and the protein. The binding free energies (ΔG_bind) of KDR–ligand complexes were calculated by using the MM-PBSA procedure in AMBER 14.0. Average 100 snapshots were selected from the last 10 ns MD simulation trajectory with a 100 ps interval for further analysis. The methodology details of binding free energy had been described in previous reports.⁶¹ The final ΔG_bind value was the average value from the last 10 ns of MD simulations.

Results and discussion

Docking analysis

Surflex-Dock was used to perform molecular docking studies of analogues described in Table S1† to explain the structure–activity relationship. To validate the docking reliability of this modeling, the small molecular ligand (compound 10) was re-docked to the binding site of protein, and the docked conformation corresponding to the lowest free energies was selected as the most possible binding conformation. Root mean squared deviation (RMSD) of the docked conformation to the experimental conformation was 0.5920 Å, and similarity (a measure of similarity between solution coordinates and reference coordinates) was 0.7190. Fig. 2(a) showed the superposition of Surflex-Dock result with the crystal structure. As can be seen that the binding conformation of compound 10 derived by Surflex-Dock superposes well with the crystal structure. Thus demonstrated that Surflex-Dock method may predict reliable conformations for KDR inhibitors. Then all molecules were docked into the active site by the same way. The docking scores versus activity for the top-ranking pose of each molecule were shown in Fig. 2(b), it could be seen that scores yielded by Surflex-Dock had a good correlation with the inhibitory activity. Fig. 2(c) displayed the 3D model of 82 compounds bound to the surface of KDR groove. Most of these inhibitors had a similar conformation to compound 10 in the X-ray structure co-crystallized with KDR (Fig. 2(d)).

Fig. 2 (a) Binding conformations of the co-crystallized compound 10 and re-docked compound 10 at the active sites of KDR. (b) Docking scores versus pIC₅₀. (c) Binding conformations of docked compounds at the surface of the groove. (d) Schematic representation of interactions between compound 10 and KDR.

QSAR study on KDR inhibitors

Based on the set of binding conformations and their atomic partial charges which were performed by using MMFF94 method, the forms of Cartesian coordinates and partial charges were took, respectively. Spatial position for each atom in a molecule and atomic charges were input into C-edited program Focus_hovf, giving rise to 3D-HoVAIF descriptors of the molecule. All obtained 165 vectors were analyzed by using SMR-PLS. According to the empirical rule, variables introduced into the model should not more than one-third of the number of samples. The 1–8 parameters, including V₇₅, V₅₈, V₅₆, V₆₄, V₁₀₀, V₇₇, V₈₂ and V₁₇ were arranged as a vector in the order of significance.

Fig. 3 was plots of correlation coefficient (R_train²) and leave one out cross validation values (Q_train²) versus the step (M) in SMR. As can be seen from the figure, R_train² and Q_train² gradually increased with the introduction of variables. The optimal variable number was determined by cross-validation Q² of the constructed PLS model. When m = 8, Q_train² reached at the maximum. It was known that this QSAR model with eight variables was the optimal model, which had the best predictive ability and good estimation capacity. It was given below:

Y = 3.566 + 0.116V₇₅ − 0.274V₆₄ + 0.303V₅₆ + 0.368V₇₇ − 0.791V₅₈ + 0.297V₁₇ − 0.184V₈₂ + 0.139V₁₀₀

Fig. 3 Value of fitness R² and cross-validation Q² changed with SMR introduced variables in the PLS model.

The two-component PLS model explained 80.1% of the sum of squares in Y-variance with a predictive ability of 75.9% for training set. R_train², Q_train², and root mean square error for estimation (RMSEE) of the optimal PLS model were 0.801, 0.759, and 0.328, respectively. The results of other modeling methods⁶² were given in Table 1. As shown in this table, the value of R² and Q² based on the process, i.e., SMR-PLS modeling, was obviously superior to others. So, it could be concluded that the model built here had a good performance for bioactivity prediction of drug molecules and it could also be used to explain the mechanism of interaction between the target and ligands. Fig. 4 showed the plots of calculated value against observed value for KDR inhibitors. Most of samples were uniformly scattered around the origin-passed diagonal. Therefore, it was confirmed that PLS models were stable and generalized.

Table 1 Comparison of the different descriptors models

No.	Modeling method	Sum of descriptors	A	Qtrain^2	Rtrain^2	Rtest^2
a Co-crystallized conformer-based alignment (CCBA) model.b Docked conformer-based alignment (DCBA) model.c Optimized molecular structures model.d Docked conformer model. Rtrain^2: cumulative multiple correlation coefficient; Qtrain^2: cross-validated correlation coefficient after leave-one-out procedure; Rtest^2: predicted correlation coefficient for test set of compounds; A: optimal number of principal components.
1	CoMFA^a	2	5	0.484	0.898	0.263
2	CoMSIA^a	5	6	0.711	0.934	0.433
3	CoMFA^b	2	5	0.546	0.936	0.673
4	CoMSIA^b	5	6	0.715	0.961	0.797
5	3D-HoVAIF^c	8	2	0.718	0.774	0.710
6	3D-HoVAIF^d	8	2	0.759	0.801	0.844

---

Fig. 4 Plots of calculated and observed activities for KDR inhibitors.

Fig. 5 showed the scoring distribution scatter of 67 samples at the first and second PLS principal component spaces, of which, circle marked samples with the pIC₅₀ greater than 9, while triangle marked samples with the pIC₅₀ between 8 and 9, and square marked samples with the pIC₅₀ less than 8. It could be seen that the distribution of KDR inhibitors exhibited an increasing trend from the left to the right in terms of their activities. Except for 2 samples (no. 45 and 68), most compounds with similar bioactivities were favorably assembled together. This suggested that the top two principal components of the 8-variable model were sufficient to reflect the essential structural characteristics for activity distribution. In addition, most samples (excluding no. 11 and no. 30) were falling into Hotelling T2 confidence ellipse with 95% confidence, which indicated that high dimensional properties of descriptors variables could construct the stable model.

Fig. 5 Scoring distribution scatters of KDR inhibitors samples at the first and second PLS principal components.

The importance of a given X-variable for Y was proportional to its distance from the origin in the variable importance in the projection (VIP) and corresponds to the PLS regression coefficients.⁵⁰ Fig. 6 delineated the PLS regression coefficients and VIP for 3D-HoVAIF descriptors. The VIP values of V₅₆, V₅₈, and V₇₅ were greater than 1.00, it was evident that steric interaction, such as the steric interaction between the first type of hydrogen atoms (Hs1) and the first type of hydrogen atoms (Hs1), the steric interaction between the first type of hydrogen atoms (Hs1) and the third type of carbon atoms (Csp2), the steric interaction between the third type of carbon atoms (Csp2) and the third type of carbon atoms (Csp2), was more important than electrostatic and hydrophobic interaction. V₅₆, V₅₈, V₇₅, V₇₇, and V₁₀₀ had great impact on the improvement of activity, while V₁₇, V₆₄ and V₈₂ played important roles in the decrease of activity for KDR inhibitors. According to values of VIP and coefficients, it could be concluded that addition of large bulk group will help to improve activity, but this group should not be too large. Because the steric interaction between s-unhybridized H atoms and s-unhybridized H atoms had a positive influence on the activity. Altering these substituents can obtain high biological activities KDR inhibitors. This result was in good agreement with the previous modelling.⁶²

Fig. 6 PLS regression coefficients and VIP for KDR inhibitors.

Designing of potent analogs

The analysis results as discussed in the previous section provided several useful information on the structural requirements for inhibitory activities. From the docking analysis, it could be found that inhibitors force the protein to adopt a “DFG-out” conformation. This conformational changed enabled inhibitors to project into the extended hydrophobic pocket. While the lipophilic pocket could not be accessed when the enzyme was in “DFG-in” conformation. The naphthyl or structurally similar unit could serve as the scaffold to build a novel series of KDR inhibitors.^17,63 To maintain pivotal interactions with the protein, the structure of designed compounds should preserve the pyridine moiety to bind hinge regions of the kinase. The amide functionality should also be retained to serve the engage of Glu 885, Asp 1046 and projected the lipophilic moiety into the extended hydrophobic pocket. Employed these information new analogs were obtained (Table S3†), which were generated by substituting the basic skeleton for compound 14 (the most potent molecule).

ADME and toxicity predictions

Most of drug candidates failed in clinical trials due to poor ADME and high toxicity profile. Thus, an important aspect was to evaluate the critical physio-chemical as well as toxicity properties in initial stages of drug discovery.⁶⁴ It was generally considered that in silico prediction of ADMET properties had a potentially significant contribution to simplification and shortening of a drug design process. So, in present work, in silico ADMET (Absorption, Distribution, Metabolism, Excretion and Toxicity) studies were performed by using ADME descriptors algorithm and toxicity prediction protocol of Accelrys Discovery studio 3.0 to determine true potential of compounds. After strictly applied all ADMET profiles, eight compounds (Fig. 7) were confirmed for further research. The calculated ADMET parameters were presented in Table S4.† All compounds showed significant values for various parameters analyzed and good drug-like characteristics based on Lipinski's rule of five.

Fig. 7 Structures and predicted pIC₅₀ values of designed compounds.

Pharmacophore modeling and reverse screening

Pharmacophore models were predicted to interact with the target binding pocket, which contained steric, electronic and hydrophobic characteristics. To explore the structural features of designed analogs, pharmacophore were developed using the GALAHAD program. In the presented figure (Fig. 8), different colored spheres represented different pharmacophore features: magenta sphere represented hydrogen-bond donor atom (D); green one represented acceptor atom (A); cyan one represented hydrophobic center (H) and red one represented positive nitrogen (P).⁵³ Among 10 best pharmacophore models, pharmacophore features such as hydrophobe (HY), hydrogen bond acceptor (HBA) and hydrogen bond donor (HBD) were included. Hydrophobe and hydrogen bond acceptor played a crucial role in determining activity against KDR. Even more startlingly, this result showed a good fit with features for previous reported modeling that was generated by potent inhibitors (Table S5†) of KDR with high inhibitory activity (Fig. 8(b)).

Fig. 8 (a) Pharmacophore model for designed analogs. (b) Pharmacophore model for potent inhibitors with high inhibitory activity. The colored spheres correspond to pharmacological features in three-dimensional representations with the following colors: magenta, hydrogen-bond donor atom; green, acceptor atom; cyan: hydrophobic center; red, positive nitrogen.

The freely accessible web server Pharmmapper served as a valuable tool for exploring the underlying toxicity of novel designed compounds. Eight analogs were aligned to 2241 human protein targets, from which the top 5 targets ranked by fit values and z-scores were analyzed (Table S6†). As can be seen from these results, most of the designed compounds with good specificity to KDR. Individual compounds were better matched with other targets, which possessed anticancer potential and played crucial roles in cancer progression mainly through four pathways: the cell cycle, angiogenesis, apoptosis, and androgen receptor.

Molecular dynamics simulations

In order to ensure the rationality of MD simulated results, the representative snapshots of KDR-10 complex took from the last 10 ns MD trajectory and the co-crystal structure were compared. As illustrated in Fig. S1,† the binding modes of compound 10 from MD simulated results were nearly the same as that extracted from the X-ray co-crystal structure. This indicated the reasonability of MD results. As shown in Fig. 9, designed compounds retained the important common binding mode known for KDR inhibitors, whose core possessed close resemblance to other inhibitors. Such as hydrophobic nitrogen or oxygen heterocycles. Most compounds targeted the ATP binding site. The hydrogen bond which was defined by distance (<3.5 Å) and orientation (the angle A⋯H–D > 120.0°)^59,65 played an important role in inhibitor binding to receptor. For the fourteen protein-inhibitor systems, hydrogen bond interactions with key residues in the binding pocket were listed in Table 2. Inspection of the docked structures showed that some of compounds fitted the binding site well and the nitrogen atom of dimethoxy quinolone ring participated in a highly conserved hydrogen-bond interaction with the backbone amide-NH of Cys919 for the hinge region. This cysteine residue was unique for KDR, which played a critical role in the recognition of ligands. The naphthyl core that was oriented orthogonal to the quinoline ring was accommodated in the mostly hydrophobic ATP-binding cleft.⁶² In addition, the carbonyl and NH groups of the amide had two interactions with the backbone NH of Asp1046 and the side chain of Glu885, respectively, which were crucial to the interaction between inhibitors and the conformation of receptor. Herein, types of hydrogen bonds included C=O⋯H–N and H–N⋯H–N. The binding pocket were surrounded by residues Leu840, Val848, Ala866, Lys868, Leu889, Thr916, Phe918, Gly922, Leu1035 and Phe1047 through hydrophobic interaction. It was evident that hydrogen bonding interactions, hydrophobic interactions and lipophilic interactions affect activity of KDR inhibitors.

Fig. 9 Predicted binding mode between KDR and inhibitors at the active site. (a) KDR-10, (b) KDR-12, (c) KDR-14, (d) KDR-30, (e) KDR-47, (f) KDR-64, (g) KDR-N3, (h) KDR-N12, (i) KDR-N13, (j) KDR-N15, (k) KDR-N16, (l) KDR-N17, (m) KDR-N18 and (n) KDR-N22.

Table 2 Hydrogen bond analysis from the results of MD simulation

System	Acceptor	Donor	Occupancy (%)	Distance (Å)	Angle (°)
KDR-10	GLU885 OE2	Ligand N12–H	82.3	2.8326	159.0112
	Ligand O20	ASP1046 NH	62.3	2.8611	155.111
	Ligand N25	CYS919 NH	36.3	2.9329	159.8882
KDR-12	Ligand O20	LYS868 NH	83.1	2.8436	157.6212
	Ligand N25	CYS919 NH	34.9	2.9218	155.1886
KDR-14	GLU885 OE2	Ligand N12–H	86.8	2.8315	162.749
	Ligand O19	ASP1046 NH	48	2.8876	159.426
	Ligand N24	CYS919 NH	34.8	2.9211	156.5414
KDR-30	GLU885 OE1	Ligand N12–H	88.1	2.8249	163.11
	Ligand N23	CYS919 NH	37.6	2.9247	156.1
	Ligand O18	ASP1046 NH	37.3	2.8946	161.0181
KDR-47	Ligand O17	ASP1046 NH	58	2.8826	156.5514
	GLU885 OE1	Ligand N12–H	54.2	2.8627	155.8766
	Ligand N22	CYS919 NH	32.4	2.9266	162.6853
KDR-64	Ligand O10	THR916 OH	83.2	2.9241	160.117
KDR-N3	GLU885 OE1	Ligand N12–H	84.4	2.8389	161.3739
	Ligand O20	ASP1046 NH	63.3	2.8796	161.3523
	Ligand N25	CYS919 NH	37.7	2.9274	162.7948
KDR-N12	GLU885 OE1	Ligand N12–H	85.1	2.8219	160.571
	Ligand O19	ASP1046 NH	64.5	2.8728	159.709
	Ligand O24	CYS919 NH	23.9	2.9114	158.642
KDR-N13	THR916 OG1	Ligand N12–H	61.2	2.8978	159.7859
KDR-N15	Ligand F20	CYS919 NH	79.1	2.7426	158.0272
	THR916 OG1	Ligand N12–H	18.2	2.9234	163.7641
KDR-N16	Ligand O20	ASP1046 NH	71.5	2.8703	161.1452
	GLU885 OE1	Ligand N12–H	59	2.828	158.5464
	Ligand O25	CYS919 NH	24.7	2.9189	158.1599
KDR-N17	THR916 OG1	Ligand N12–H	55.8	2.9183	153.1568
	Ligand O31	ASN923 NH	27	2.8886	160.1951
	Ligand N10	CYS919 NH	26.4	2.9278	156.4941
KDR-N18	THR916 OG1	Ligand N11–H	58.5	2.8697	152.6839
KDR-N22	GLU885 OE1	Ligand N11–H	60.4	2.8074	163.9582
	Ligand O25	CYS919 NH	29.8	2.8065	164.0067

---

To further understand the probable binding modes of KDR–ligand complexes, MD simulations which had been widely applied to study protein–ligand interactions were used to discuss the dynamical structural features. At first, MD simulations were carried out on KDR-compound 10 (KDR-10) binding structure obtained from molecular docking. Fig. 10(a) showed the system kept stable on 100 ns molecular dynamics simulation. The RMSD values of KDR-10 complex obtained equilibration at 2.0 Å after 20 ns. The ligand exhibited different dynamics behaviors in the active site of KDR. The RMSD values fluctuate around 0.7 Å during the first 75 ns. Afterward, a sharp rise was observed and the trajectories kept stable with average RMSD values of 1.8 Å for about 10 ns. During about 85 ns, a sudden drop was observed and then the trajectories kept stable for about 3 ns. At 88 ns, a sharp rise was also observed and then the trajectories kept stable with average RMSD values of 1.8 Å in the rest of the simulation. This indicated that the conformation of compound 10 underwent a notable change along the simulation. The behavior was characteristic of the molecule itself and not an effect of the protein conformation. The average structures of compound 10 from 0 to 75 ns, 75 to 85 ns, 85 to 88 ns and 88 to 100 ns MD trajectories were displayed in the ESI (Fig. S2†). It was obvious that the flip took place in the terminal 4-Cl phenyl. Furthermore, analyses of root-mean-square fluctuation (RMSF) versus the residue number for KDR-10 complex were illustrated in Fig. 10(b). It could be found that residues with higher fluctuation values were those in the flexible loops, while alpha helixes and beta strands were more rigid.

Fig. 10 (a) RMSDs of backbone atoms (C, Cα, and N) for the KDR-10 complex system and the heavy atoms in the ligand. (b) RMSF of each residue for the KDR-10 complex system.

It was not surprised that compound 12, 14, 30 and 64 with the same inhibitory activity shared the similar MD trajectories (Fig. 11). For complexes with compound 14, 30 and 64, the RMSD values obtained equilibration at 1.8 Å and 2.0 Å after 20 ns MD simulation. Experiencing the sharp rise and sudden drop, complex with compound 12 had no obvious fluctuation after 50 ns MD simulation with RMSD equilibration at 1.8 Å. Fig. 11(b) showed the RMSD values of four ligands. Except the compound 14 had fluctuations at 22 ns and 58 ns, it could be noted that other ligands were stable throughout the simulations, with average of 0.88 Å and 2.0 Å RMSD fluctuations, respectively. As shown in Fig. 11(c), the RMSF profiles of all complexes had a similar trend.

Fig. 11 (a) RMSDs of backbone atoms (C, Cα, and N) for KDR–ligand complex systems. (b) RMSDs of heavy atoms for ligands. (c) RMSF of each residue for KDR–ligand complex systems.

In order to explore the differences between the most active and the least active inhibitors in the binding process, MD simulations for compounds 14 and 47 bound KDR in explicit aqueous solution were also ran for 100 ns. The plots in Fig. 12 revealed that studied systems reach equilibrium within 30 ns. The protein and ligand in compound 14 bound system were stable after equilibrium, with average 1.8 Å RMSD fluctuation. The fluctuations were very small. While the corresponding RMSD values for the compound 47 bound system were significantly larger with longer time to equilibrium. Compound 47 showed the highest fluctuations during the MD simulations (Fig. 12(b)). This indicated that the conformation of this ligand underwent a notable change, and showed more flexible to a degree. Whereas the compound 14 was more stabilized. The RMSF values of two systems showed that fluctuations were higher in 47-bound system than that in 14-bound system (Fig. 12(c)). Residues with higher fluctuation values were located in turn regions between α helices and β sheets, the activation loop region, or the protein terminus.

Fig. 12 (a) RMSDs of backbone atoms (C, Cα, and N) for KDR-14/47 complex systems. (b) RMSDs of heavy atoms for ligands. (c) RMSF of each residue for KDR-14/47 complex systems.

In this computational study, novel analogs with maintained or improved activity were designed. After molecular docking, ligand binding modes need deeply simulation optimization. As can be seen from the plots (Fig. 13(a)), the RMSD of each system except the N16 complex tended to convergence after 5 ns MD simulation, indicating the systems were equilibrated and stable. For complexes with N3, N13, N17, N18 and N22, the RMSD values obtained equilibration at 2.0 Å. Complexes with N12 and N15 had no obvious fluctuation with RMSD equilibration at 1.8 Å. RMSD values of the N16 bound system fluctuated around 1.8 Å during the first 40 ns. Then this system experienced an obvious drift and kept stable with average RMSD values of 2.0 Å in the rest of simulations. This implied that the conformation of ligand N16 underwent a notable change (Fig. 13(b)). During the whole simulation, designed compound N13, N15, N18 and N22 remained highly stable at 0.7 Å, 1.5 Å, 2.5 Å and 2 Å, respectively. RMSD of ligand N3 and N17 stabilized around 1.5 Å. The ligand N12 had some conformation changes with the RMSD changing from 0.5 Å to 1.5 Å. Interestingly, the RMSD value of KDR-N16 backbone atoms stabilized at 1.8 Å for the first 40 ns, but the ligand occurred large fluctuations at 10 ns. It hardly changed the conformation during the 10 to 75 ns with RMSD value equilibrating at 1.5 Å. Then it occurred a sudden drop and kept stable with average RMSD values of 0.5 Å. Fig. 13(c) illustrated the RMSF values versus the residue number for complexes. It could be seen that all systems shared similar distributions and trends for dynamic features. This suggested that all inhibitors had similar interaction mechanism with KDR.

Fig. 13 (a) RMSDs of backbone atoms (C, Cα, and N) for KDR-analogs complex systems. (b) RMSDs of heavy atoms for ligands. (c) RMSF of each residue for KDR-analogs complex systems.

Binding free energy analysis

Binding free energies of all systems were calculated using the MM-PBSA method in AMBER 14.0.⁶⁶ As indicated in Table 3, the calculated ΔG was greatly compatible with the experimental data (ΔG_exp), which was determined based on IC₅₀ values (ΔG ≈ −RT ln IC₅₀).⁶⁷ MM-PBSA calculations confirmed that ligands 12, 14, 30 and 64 had the same activities against KDR (ΔG_KDR-12 = −15.50 kcal mol⁻¹, ΔG_KDR-14 = −15.79 kcal mol⁻¹, ΔG_KDR-30 = −15.52 kcal mol⁻¹ and ΔG_KDR-64 = −15.04 kcal mol⁻¹). The results were in good agreement with experimental datum. There was a big difference of the ΔG for compounds 14 and 47, which corresponded to the difference in the binding affinity. The IC₅₀ value of compound 47 was about 7589 times higher than that of compound 14. The ΔG_exp was estimated to be 5.29 kcal mol⁻¹ difference for compounds 14 and 47. Furthermore, the ΔG were −15.19 to −18.71 kcal mol⁻¹ for the complexes of KDR-N3, N12, N13, N15, N16, N17, N18 and N22. The compound N3 bound to KDR was 2.92 kcal mol⁻¹ more stable than ligands 14 bound to KDR. It was noteworthy that the ΔG of some complexes were lower than that derived from the experiment. This was caused by the low prediction accuracy of entropic calculation and the approximation of experimental free energy. It was still of great significance to compare the relative value.⁶⁰

Table 3 Binding free energies and its components for potent KDR inhibitors^a

Compound	ΔGvdW	ΔGele	ΔGele,sol	ΔGnonpol,sol	TΔS	ΔGpred	ΔGexp
a All energies were in kcal mol^−1. TΔS: the entropy changes. ΔGpred: the calculated binding free energy by MM-PBSA method. ΔGexp: the experimental binding free energy calculated according to the IC50 by ΔG ≈ −RT ln IC50.
10	−63.63	−33.31	49.70	−44.02	−76.50	−14.76	−12.71
12	−68.44	−19.25	40.25	−44.23	−76.17	−15.50	−13.25
14	−61.61	−28.91	44.53	−42.56	−72.76	−15.79	−13.25
30	−60.68	−33.92	50.23	−40.88	−69.73	−15.52	−13.25
47	−62.11	−29.60	46.41	−43.20	−73.61	−14.89	−7.95
64	−60.67	−28.48	43.70	−43.39	−73.80	−15.04	−13.25
N3	−61.46	−30.65	44.98	−42.79	−71.21	−18.71	−14.24
N12	−53.82	−30.28	44.68	−40.04	−63.36	−16.10	−13.59
N13	−66.63	−26.80	45.82	−45.67	−77.29	−15.99	−13.51
N15	−66.56	−27.81	47.70	−45.83	−77.29	−15.21	−13.27
N16	−65.61	−33.36	50.68	−44.55	−77.45	−15.39	−13.26
N17	−67.88	−34.13	52.98	−46.27	−80.11	−15.19	−13.20
N18	−61.99	−26.76	42.56	−38.43	−69.18	−15.44	−13.30
N22	−63.07	−21.01	39.27	−44.06	−72.95	−15.92	−13.33

---

Enable to understand the protein–ligand complex binding process in detail, compared which energy term had more impact on the binding affinity. The total binding free energy was deconvolute into individual components. In the fourteen complexes, the van der Waals interactions (ΔG_vdw) and the nonpolar solvation energies (ΔG_nonpol,sol) were the basis for favorable binding free energies, which were responsible for imbedded hydrophobic groups for inhibitors. The favorable Coulomb interactions within the protein–ligand complexes were counteracted by detrimental electrostatics of desolvation. Consequently, the contribution of electrostatic interaction were unfavorable for binding in all systems.⁵⁹ The contributions of entropy changes (−TΔS) to free energies impaired the binding of each inhibitor to KDR. From Table 3, it was evident that the value of ΔG_vdw was greater than that of ΔG_ele. The ΔG_vdw was dominant in ΔG. Examined the contributions to each binding energy, it could be found that the summation of electrostatic energy including electrostatic interaction (ΔG_ele) and polar contribution to solvation (ΔG_ele,sol) was similar. The summation of ΔG_vdw and ΔG_nonpol,sol for KDR-47 was less than that of the 12, 14, 30 and 64-bound. This difference could partial interpretation why the conformational change of KDR-47 at binding pocket led to its reduced binding affinity. The ΔG_vdw of KDR-N3 complex (−61.46 kcal mol⁻¹) was approximately equal to KDR-14 complex (−61.61 kcal mol⁻¹). Obviously, the ΔG_vdw part could not account for the activity difference between ligands N3 and ligands 14. ΔG_nonpol,sol for these two complexes (−42.79 kcal mol⁻¹ and −42.56 kcal mol⁻¹) were approximately equivalent. While the difference of electrostatic energy between KDR-N3 and KDR-14 was −1.29 kcal mol⁻¹. Therefore, the electrostatic energy played the most important role in differentiating the binding affinity between KDR-N3 and KDR-14.

Binding free energy deconvolute

To identify the key residues related with the binding process, binding free energies between the protein and inhibitors were deconvoluted into the contribution of each residue. As can be seen from Fig. S3,† the major favorable energy contributions originated from the residues Leu840, Val848, Ala866, Lys868, Leu889, Val899, Thr916, Phe918, Cys919, Leu1035, Cys1045, Asp1046, and Phe1047 in all complexes. Most of these important residues were hydrophobic, which could form strong van der Waals interactions with inhibitors. The contributions of van der Waals and nonpolar solvation energies of residues Leu840, Val899, Cys919, Leu1035 and Phe1047 were mainly responsible for the differences between N3, N12, N22 and the others. And the total electrostatic energies of residues Leu840, GLU885, GLU917, and LYS920 participated in yielding energetic differences for ligand 47, compared with the others. From these results, it could be concluded that hydrogen bond and hydrophobic interactions played an important role in the binding affinity with residues Leu840, GLU885, Val899, Cys919, Asp1046, and Phe1047 (Fig. 14).

Fig. 14 Comparison of per-residue energy decomposition for key residues of fourteen systems: (a) and (c) the sum of vdW and nonpolar solvation energy, ΔG_vdW + ΔG_nonpol,sol, and (b) and (d) the sum of electrostatic and polar solvation energy, ΔG_ele + ΔG_ele,sol.

Conclusions

KDR, alternatively referred to as VEGFR-2, was a receptor for VEGFs. It was directly linked to tumor angiogenesis and blood vessel-dependent metastasis.^13,16 Inhibitors of KDR had been demonstrated to induce tumor regression and reduce metastatic potential in preclinical models.^55,68 In this study, a series of CADD processes, such as molecular docking, 3D MSC method, pharmacophore modeling, reverse docking, molecular dynamics and MM/PBSA calculations were used to build models and explore the probable binding modes between ligands and the receptor protein. Meanwhile, the CADD technology was also employed to design new potent inhibitors with high activities. The model showed a good internal and external predictive ability, with the Q_train² was 0.759, the R_train² of training datasets was 0.801 and the R_test² of test datasets was 0.844, respectively. Quantitative models constructed here indicated that steric interaction played more important roles in their bioactivities than electrostatic interaction. Adding large bulk groups will help to improve activity. According these information eight novel analogs with desired activity were designed. Then the complexes binding modes gained from molecular dynamics were used to analyze the affinity of ligands in the binding pocket. Calculations of binding free energy verified that ligands 12, 14, 30 and 64 had almost the same activities against KDR. This conclusion was greatly compatible with the experimental observation. MD simulations and binding free energy calculations of ligand 47 bound to KDR could give a good explanation of the drastically diminished potency. In addition, free energy calculations suggested that new designed compounds bound to and stabilized the conformation of KDR. The decomposition of binding free energy into each interaction type indicated van der Waals interactions that supplied the substantial driving force for the binding process. The most favorable contributions were came from Leu840, Val848, Ala866, Lys868, Leu889, Val899, Thr916, Phe918, Cys919, Leu1035, Cys1045, Asp1046, and Phe1047. So, the combined strategy presented by this paper could helpful to understand the structural features of target and chemical features of ligands, and provide a basis for further rational design of new potent inhibitors.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (81171508, 31170747, 31400667), Key Project of Natural Science Foundation of Chong Qing (CSTC 2013JJB10004), and Science and Technology project of Chongqing Education Commission (KJ1400941, KJ1400932).

References

1. N. Ferrara, Nat. Rev. Cancer, 2002, 2, 795–803.
2. N. Ferrara, H. P. Gerber and J. LeCouter, Nat. Med., 2003, 9, 669–676.
3. J. M. Cherrington, L. M. Strawn and L. K. Shawver, Adv. Cancer Res., 2000, 79, 1–38.
4. D. W. Siemann, D. J. Chaplin and M. R. Horsman, Cancer, 2004, 100, 2491–2499.
5. S. D. Zhang, C. M. McCrudden and H. F. Kwok, Oncol. Lett., 2015, 10, 1893–1901.
6. S. D. Zhang, K. L. Leung, C. M. McCrudden and H. F. Kwok, J. Cancer, 2015, 6, 812–818.
7. J. Folkman, Ann. Surg., 1972, 175, 409–416.
8. J. Rak, Y. Mitsuhashi, L. Bayko, J. Filmus, S. Shirasawa, T. Sasazuki and R. S. Kerbel, Cancer Res., 1995, 55, 4575–4580.
9. G. Breier and W. Risau, Trends Cell Biol., 1996, 6, 454–456.
10. A. A. Friedman, A. Amzallag, I. Pruteanu-Malinici, S. Baniya, Z. A. Cooper, A. Piris, L. Hargreaves, V. Igras, D. T. Frederick, D. P. Lawrence, D. A. Haber, K. T. Flaherty, J. A. Wargo, S. Ramaswamy, C. H. Benes and D. E. Fisher, PLoS One, 2015, 10, e0140310.
11. L. Scotti, D. Abramovich, N. Pascuali, G. Irusta, G. Meresman, M. Tesone and F. Parborell, J. Steroid Biochem. Mol. Biol., 2014, 144(Pt B), 392–401.
12. D. Mukhopadhyay, B. Knebelmann, H. T. Cohen, S. Ananth and V. P. Sukhatme, J. Steroid Biochem. Mol. Biol., 1997, 17, 5629–5639.
13. M. L. Maitland, C. F. Xu, Y. C. Cheng, E. Kistner-Griffin, K. A. Ryan, T. G. Karrison, S. Das, D. Torgerson, E. R. Gamazon, V. Thomeas, M. R. Levine, P. A. Wilson, N. Bing, Y. Liu, L. R. Cardon, L. N. Pandite, J. R. O'Connell, N. J. Cox, B. D. Mitchell, M. J. Ratain and A. R. Shuldiner, Clin. Cancer Res., 2015, 21, 365–372.
14. D. Shweiki, M. Neeman, A. Itin and E. Keshet, Proc. Natl. Acad. Sci. U. S. A., 1995, 92, 768–772.
15. M. B. Nilsson, U. Giri, J. Gudikote, X. Tang, W. Lu, H. Tran, Y. Fan, A. Koo, L. Diao, P. Tong, J. Wang, R. Herbst, B. E. Johnson, A. Ryan, A. Webster, P. Rowe, I. I. Wistuba and J. V. Heymach, Clin. Cancer Res., 2015 DOI:10.1158/1078-0432.ccr-15-1994.
16. X. Zhang, Y. L. Ge, S. P. Zhang, P. Yan and R. H. Tian, Cell. Mol. Biol. Lett., 2014, 19, 527–541.
17. D. G. Kim, Y. Jin, J. Jin, H. Yang, K. M. Joo, W. S. Lee, S. R. Shim, S. W. Kim, J. Yoo, S. H. Lee, J. S. Yoo and D. H. Nam, mAbs, 2015, 7, 1195–1204.
18. A. K. Halder, S. Mallick, D. Shikha, A. Saha, K. D. Saha and T. Jha, RSC Adv., 2015, 5, 72373–72386.
19. M. Lv, S. Ma, Y. Tian, X. Zhang, H. Zhai and W. Lv, RSC Adv., 2015, 5, 462–476.
20. S. K. Saha and P. Banerjee, RSC Adv., 2015, 5, 71120–71130.
21. L. Hu, T. Huang, X. Shi, W. C. Lu, Y. D. Cai and K. C. Chou, PLoS One, 2011, 6, e14556.
22. A. Cherkasov and B. Jankovic, Molecules, 2004, 9, 1034–1052.
23. Z. He, J. Zhang, X. H. Shi, L. L. Hu, X. Kong, Y. D. Cai and K. C. Chou, PLoS One, 2010, 5, e9603.
24. A. Genoni, M. Pennati, G. Morra, N. Zaffaroni and G. Colombo, RSC Adv., 2012, 2, 4268–4282.
25. M. Zhou, H. Luo, R. Li and Z. Ding, RSC Adv., 2013, 3, 22532–22543.
26. H. Wiener, J. Phys. Colloid Chem., 1948, 52, 1082–1089.
27. H. Wiener, J. Phys. Colloid Chem., 1948, 52, 425–430.
28. L. Yang, M. Shu, K. Ma, H. Mei, Y. Jiang and Z. Li, Amino Acids, 2010, 38, 805–816.
29. M. Shu, Y. Jiang, L. Yang, Y. Wu, H. Mei and Z. Li, Protein Pept. Lett., 2009, 16, 143–149.
30. X. Y. Wang, J. Wang, Y. Lin, Y. Ding, Y. Wang, X. Cheng and Z. Lin, J. Mol. Model., 2011, 17, 1599–1606.
31. J. Wang, X. Y. Wang, M. Shu, Y. Q. Wang, Y. Lin, L. Wang, X. M. Cheng and Z. H. Lin, Protein Pept. Lett., 2011, 18, 956–963.
32. T. Simonson, G. Archontis and M. Karplus, Acc. Chem. Res., 2002, 35, 430–437.
33. H. Gohlke and G. Klebe, Angew. Chem., 2002, 41, 2644–2676.
34. Y. Deng and B. Roux, J. Phys. Chem. B, 2009, 113, 2234–2246.
35. T. Steinbrecher and A. Labahn, Curr. Med. Chem., 2010, 17, 767–785.
36. N. Singh and A. Warshel, Proteins: Struct., Funct., Genet., 2010, 78, 1705–1723.
37. M. Levitt, J. Mol. Biol., 1983, 170, 723–764.
38. M. Levitt and M. F. Perutz, J. Mol. Biol., 1988, 201, 751–754.
39. G. E. Kellogg, S. F. Semus and D. J. Abraham, J. Comput.-Aided Mol. Des., 1991, 5, 545–552.
40. A. R. Katritzky, U. Maran, V. S. Lobanov and M. Karelson, J. Chem. Inf. Comput. Sci., 2000, 40, 1–18.
41. J. C. Harmange, M. M. Weiss, J. Germain, A. J. Polverino, G. Borg, J. Bready, D. Chen, D. Choquette, A. Coxon, T. DeMelfi, L. DiPietro, N. Doerr, J. Estrada, J. Flynn, R. F. Graceffa, S. P. Harriman, S. Kaufman, D. S. La, A. Long, M. W. Martin, S. Neervannan, V. F. Patel, M. Potashman, K. Regal, P. M. Roveto, M. L. Schrag, C. Starnes, A. Tasker, Y. Teffera, L. Wang, R. D. White, D. A. Whittington and R. Zanon, J. Med. Chem., 2008, 51, 1649–1667.
42. M. M. Weiss, J. C. Harmange, A. J. Polverino, D. Bauer, L. Berry, V. Berry, G. Borg, J. Bready, D. Chen, D. Choquette, A. Coxon, T. DeMelfi, N. Doerr, J. Estrada, J. Flynn, R. F. Graceffa, S. P. Harriman, S. Kaufman, D. S. La, A. Long, S. Neervannan, V. F. Patel, M. Potashman, K. Regal, P. M. Roveto, M. L. Schrag, C. Starnes, A. Tasker, Y. Teffera, D. A. Whittington and R. Zanon, J. Med. Chem., 2008, 51, 1668–1680.
43. D. S. La, J. Belzile, J. V. Bready, A. Coxon, T. DeMelfi, N. Doerr, J. Estrada, J. C. Flynn, S. R. Flynn, R. F. Graceffa, S. P. Harriman, J. F. Larrow, A. M. Long, M. W. Martin, M. J. Morrison, V. F. Patel, P. M. Roveto, L. Wang, M. M. Weiss, D. A. Whittington, Y. Teffera, Z. Zhao, A. J. Polverino and J. C. Harmange, J. Med. Chem., 2008, 51, 1695–1705.
44. A. N. Jain, J. Comput.-Aided Mol. Des., 2007, 21, 281–306.
45. A. N. Jain, J. Med. Chem., 2003, 46, 499–511.
46. F. Tian, P. Zhou, F. Lv, R. Song and Z. Li, J. Pept. Sci., 2007, 13, 549–566.
47. L. N. Xu, G. Z. Liang, Z. L. Li, J. Wang and P. Zhou, J. Mol. Graphics Modell., 2008, 26, 1252–1258.
48. M. Shu, Y. R. Zhang, F. F. Tian, L. Yang and Z. H. Lin, Chin. J. Struct. Chem., 2012, 31, 443–451.
49. P. Zhou, F. F. Tian and Z. L. Li, Chemom. Intell. Lab. Syst., 2007, 87, 88–94.
50. M. Henningsson, E. Sundbom, B. A. Armelius and P. Erdberg, Scand. J. Psychol., 2001, 42, 399–409.
51. M. Arooj, S. Thangapandian, S. John, S. Hwang, J. K. Park and K. W. Lee, Int. J. Mol. Sci., 2011, 12, 9236–9264.
52. K. Meraj, M. K. Mahto, N. B. Christina, N. Desai, S. Shahbazi and M. Bhaskar, Bioinformation, 2012, 8, 1139–1146.
53. N. J. Richmond, C. A. Abrams, P. R. Wolohan, E. Abrahamian, P. Willett and R. D. Clark, J. Comput.-Aided Mol. Des., 2006, 20, 567–587.
54. X. Liu, S. Ouyang, B. Yu, Y. Liu, K. Huang, J. Gong, S. Zheng, Z. Li, H. Li and H. Jiang, Nucleic Acids Res., 2010, 38, W609–W614.
55. T. Usha, A. K. Goyal, S. Lubna, H. Prashanth, T. M. Mohan, V. Pande and S. K. Middha, Asian Pacific Journal of Cancer Prevention, 2014, 15, 10345–10350.
56. M. Kolar, K. Berka, P. Jurecka and P. Hobza, ChemPhysChem, 2010, 11, 2399–2408.
57. L. K. Lindorff, S. Piana, K. Palmo, P. Maragakis, J. L. Klepeis, R. O. Dror and D. E. Shaw, Proteins: Struct., Funct., Genet., 2010, 78, 1950–1958.
58. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
59. Y. Yang, J. Qin, H. Liu and X. Yao, J. Chem. Inf. Model., 2011, 51, 680–692.
60. Y. Yang, Y. Shen, H. Liu and X. Yao, J. Chem. Inf. Model., 2011, 51, 3235–3246.
61. A. Jakalian, D. B. Jack and C. I. Bayly, J. Comput. Chem., 2002, 23, 1623–1641.
62. J. Du, B. Lei, J. Qin, H. Liu and X. Yao, J. Mol. Graphics Modell., 2009, 27, 642–654.
63. L. Ding, F. Tang, W. Huang, Q. Jin, H. Shen and P. Wei, Bioorg. Med. Chem. Lett., 2013, 23, 5630–5633.
64. S. Ray, P. B. Madrid, P. Catz, S. E. LeValley, M. J. Furniss, L. L. Rausch, R. K. Guy, J. L. DeRisi, L. V. Iyer, C. E. Green and J. C. Mirsalis, J. Med. Chem., 2010, 53, 3685–3695.
65. X. Wu, S. Wan, G. Wang, H. Jin, Z. Li, Y. Tian, Z. Zhu and J. Zhang, J. Mol. Graphics Modell., 2015, 56, 103–112.
66. K. Kitamura, Y. Tamura, T. Ueki, K. Ogata, S. Noda, R. Himeno and H. Chuman, J. Chem. Inf. Model., 2014, 54, 1653–1660.
67. J. Wang, P. Morin, W. Wang and P. A. Kollman, J. Am. Chem. Soc., 2001, 123, 5221–5230.
68. B. Vijayakumar, S. Parasuraman, R. Raveendran and D. Velmurugan, Pharmacogn. Mag., 2014, 10, S639–S644.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra03743g

---