Molecular Dynamics Simulation of VEGFR2 with Sorafenib and Other Urea-Substituted Aryloxy Compounds
The binding mode of sorafenib with VEGFR2 was studied using molecular docking and molecular dynamics method. The docking results show that sorafenib forms hydrogen bonds with Asp1046, Cys919, and Glu885 of VEGFR2 receptor. Molecular dynamics simulation suggests that the hydrogen bond involving Asp1046 is the most stable one, and it is almost preserved during all the MD simulation time. The hydrogen bond formed with Cys919 occurs frequently after 6 ns, while the bifurcated hydrogen bonds involving Glu885 occurs occasionally. Meantime, molecular dynamics simulations of VEGFR2 with 11 other urea-substituted aryloxy compounds have also been performed, and the results indicate that a potent VEGFR2 inhibitor should have lower interaction energy with VEGFR2 and create at least 2 hydrogen bonds with VEGFR2.
Vascular endothelial growth factor (VEGF) is an important signaling protein involved in both the growth of blood vessels from preexisting vasculature (angiogenesis) and the formation of the circulatory system (vasculogenesis). VEGF binding to tyrosine kinase receptors (VEGFR) can cause itself dimerization and become activated through transphosphorylation. There are three main subtypes of VEGFR: VEGFR1, VEGFR2, and VEGFR3. Among which VEGFR2 appears to mediate almost all of the known cellular responses to VEGF. Inhibiting the tyrosine kinase VEGFR2 signaling pathway may disrupt the angiogenesis process of solid tumor, thus blocking tumor growth and spread . Therefore, the design of inhibitors targeting VEGFR2 is an attractive approach for the development of new therapeutic agents.
Recently, lots of VEGR2 inhibitors have been reported. Among them sorafenib (codeveloped and comarketed by Bayer and Onyx pharmaceuticals as Nexavar) is one of the potent inhibitors of VEGFR in vitro and is approved by FDA for the treatment of advanced renal cell carcinoma and advanced primary liver cancer [2, 3]. Although sorafenib treatment of melanoma cell lines and tumor xenografts leads to cell death and tumor growth delay [4–6], it has little or no antitumor activity in advanced melanoma patients when used as a single agent . Moreover, its water solubility is not satisfactory . Thus, there are still needs to develop new VEGFR2 inhibitors. With the development of the computational software and hardware, the usage of in silico methods, such as docking , QSAR, pharmacophore modeling , and molecular dynamics method , has become essential to pharmaceutical research. Abreu et al.  have studied the VEGFR2-selective side-chain residue flexibility using AutoDock Vina docking software, and they found that Glu885 flexible conformation could present best scores. An and coworkers  have used molecular dynamics method to investigate the complex of VEGFR2 and its inhibitor sunitinib, and they found that sunitinib could form hydrogen bonds with three residues.
Recently Garofalo and coworkers  have developed a novel series of aryloxy-linked quinazolines with good VEGFR2 inhibitory. They have same quinazoline ring-like gefitinib and vandetanib, among which 11 compounds have the same urea scaffold as sorafenib. Illustrating the binding mode of such inhibitors should be helpful to understand the inhibitory effects caused by different groups. An improved understanding of this kind of VEGFR2 inhibitors may be useful to the design and synthesis of new potential VEGFR2 inhibitors. Thus to this end, molecular docking and molecular dynamics simulations were employed to study the binding mode of sorafenib and other urea-substituted aryloxy compounds with VEGFR2 in this paper.
2. Computational Methods
2.1. Molecular Docking
The structure of VEGFR2 receptor was taken from Protein Data Bank with the ID code 2RL5  (Figure 1). The receptor was prepared using Protein Prepare wizard. The crystal ligand 2RL was used as the docking center with a box size of 12 Å. The structure of sorafenib was sketched in Maestro. LigPrep was used to perform ligand preparation. Glide Module of Schrödinger software  was used to carry out molecular docking computation using default parameters.
2.2. Molecular Dynamics of VEGFR2-Sorafenib
Both the unbound and bound VEGFR2 were used as the initial structures for MD simulation. The topology parameter of sorafenib was obtained through Dundee PRODRG server (http://davapc1.bioch.dundee.ac.uk/programs/prodrg/) . The complex was placed in a dodecahedron box with an edge of 1.2 nm, then the box was solvated with SPC water molecules. Two chloride ions were added to neutralize the whole system. Energy minimization was performed on the complex. Subsequently, 100 ps NVT and 100 ps NPT MD simulation at 300 K were performed with position restrictions on protein and ligand. Finally, MD simulation was carried out on the whole system for 10 ns. Particle Mesh Ewald (PME) method was used to control long range electrostatic interaction. Cutoff of van der Waals was 1.4 nm. The V-rescale and Parrinello-Rahman methods were used for temperature and pressure coupling, respectively. The reference temperature was set as 300 K. Time step was 2 fs and the trajectory was saved every 4 ps. All the simulation was carried out using GROMACS 4.5.4 package [18–21]with GROMOS96 53a6 force field as this force field was validated and successfully used in many systems including protein and drug molecules [22, 23].
2.3. Molecular Dynamics of Other Aryloxy-Quinazolines Derivatives 12–22
The structures of 11 aryloxy-quinazolines derivatives 12–22 are shown in Figure 2. These compounds have been sketched in Maestro with the docking pose of sorafenib as the reference structure. Then the complex structure was minimized using Embrace method of Macromodel module  in Schrödinger software. The obtained complex was subjected to molecular dynamics simulation as for VEGFR2-sorafenib complex.
2.4. Energy Analysis
Compounds 12–22 have the same scaffold as sorafenib. So we divided them into three parts: A, B, and C (Figure 2). Part B was the same for all the studied compounds. The interaction energy between VEGFR2 and the inhibitor was also divided into three parts according to the division of inhibitor. Total interaction energy is the sum of van der Waals and Coulomb energies. The van der Waals and Coulomb energies between certain part of inhibitor and VEGFR2 were computed using the following equations: where is the distance between the and atoms in selected part of drug and in receptor; and are the corresponding coefficients of these two atoms taken from GROMOS force field; and are the corresponding charge numbers of the two atoms in drug and in receptor, respectively. and are the total atom number of certain part and VEGFR2, respectively.
3. Results and Discussion
3.1. Molecular Docking of VEGFR2-Sorafenib
Figure 3 is the obtained docking mode of sorafenib with VEGFR2. The docking score of sorafenib is −7.36 with a glide energy of −61.19 kJ·mol−1. It can be seen from Figure 3 that sorafenib forms hydrogen bonds with three residues: Asp1046, Cys919, and Glu885. Cys919 forms a weak bifurcated hydrogen bond with the nitrogen atom of pyridine ring and carboxyl oxygen atom (the two bond distances are 2.661 and 2.698 Å, resp.), while Asp1046 forms hydrogen bond with the carboxyl oxygen atom of sorafenib (the distance is 2.422 Å). OE2 atom of Glu885 forms bifurcated hydrogen bonds with the two N–H bonds of urea group with a bond distance of 1.863 Å. From Figure 3 we also can see that there exists a strong salt bridge between Lys868 and Glu885. The interaction between Lys868 and Glu885 may weaken the hydrogen bonds involving Glu885. Except hydrogen bond interaction, hydrophobic interaction plays an import role in stabilizing VEGFR2-sorafenib complex, which can be seen from (see supplementary Figure 1 of the Supplementary Material available online at http://dx.doi.org/10.1155/2013/739574). The hydrophobic pocket is mainly made up of Leu840, Val848, Ala866, Ile888, Leu889, Val899, Phe918, Thr916, and Leu1035.
3.2. Molecular Dynamics of VEGFR2-Sorafenib
3.2.1. Reliability of the MD Simulation
Figure 4(a) shows the root mean square deviation (RMSD) profiles of backbone atoms for VEGFR2 and VEGFR2-sorafenib. It can be seen from Figure 4(a) that the overall RMSD of VEGFR2-sorafenib complex after simulation for 10 ns was about 0.4 nm. The RMSD increases rapidly during the first 1 ns, which means that the structure of VEGFR2 deviated from the initial structure rapidly during the 1 ns simulation time. This increase in the RMSD was due to the optimization of interactions with the protein structure, as well as with the water solvent. The overall RMSD was finally fluctuated in a range of 0.1~0.2 nm, indicating that the system was stabilized. Hence the last 2 ns of MD production were satisfactory for analysis. Radius of gyration () of unbound and bound VEGFR2 complexes was analyzed to determine the effect of sorafenib on the folding of VEGFR2. The variation of value is found almost the same after 2 ns, which can be seen from Figure 4(b). All these results indicate that the trajectory of MD simulation after equilibrium was reliable for postanalysis.
3.2.2. Hydrogen Bonds
Figure 5 shows the distance fluctuations of hydrogen bonds (distances between two heavy atoms) between sorafenib and VEGFR2. The distance of OAT and N2 (Asp1046) is almost less than 0.35 nm and the corresponding hydrogen bond angle is less than 30° (see supplementary Figure 2), suggesting that this hydrogen bond is a persistent hydrogen bond. At the beginning of MD simulation the distance between NAJ atom and N1 atom of Cys919 becomes longer and longer, and it reaches 1.0 nm at 3 ns, then it fluctuates about 0.6~0.9 nm, indicating the breakage of this hydrogen bond. But after 6 ns this hydrogen bond reformed (the hydrogen bond distance is about 0.35 nm). The variation of the hydrogen bonds between NAU and NAR atoms with OE2 atom of Glu885 are similar (Figures 5(c) and 5(d)). Both the distance and bond angle fluctuate in a relative large range (0.3~1.0 nm for bond distance and 10°~90° for bond angel), suggesting that these two hydrogen bonds are much more unstable as compared with the other two hydrogen bonds.
Through the bond distance and bond angle variation results, we can see that the hydrogen bond formed with Asp1046 is stable during the MD simulation, which with Cys919 is less stable, while the bifurcated hydrogen bond with Glu885 is the most unstable one.
Figure 6 is the hydrogen bond existence map for VEGFR2-sorafenib complex, from which we can see that one persistent hydrogen bond (hydrogen bond index 0) forms, indicating that the carboxyl group of the sorafenib makes a stable hydrogen bond with Asp1046. As for the hydrogen bond involving Cys919 (hydrogen bond index 1, 2, and 4), frequent breaking and reforming is observed due to the competition of the three acceptor atoms (NAJ, NAB, and OAD) in sorafenib for the same donor of Cys919. The two hydrogen bonds with Gly841 and Thr916 rarely occurred. Glu885 forms unstable hydrogen bond with the two NH of urea group. This coincides with the previous hydrogen bond distance and bond angle results.
3.3. Molecular Dynamics of Aryloxy-Quinazolines Derivatives
From what was previously mentioned we can see that the urea linkage of sorafenib could form hydrogen bonds with ASP1046 and GLU885. Besides, the nitrogen atom of pyridine ring also could form hydrogen bond with CYS919. Holding this in mind we further studied the interaction mode between VEGFR2 and a serial of urea derivatives with such functional groups .
Supplementary Figure 3 lists the hydrogen bond existence of the 11 compounds with VEGFR2 protein, from which we can see that all the compounds forms hydrogen bond with Asp1046 and the occurrence is almost persistent, suggesting that this hydrogen bond is very stable. Among all the 11 complexes, only three form hydrogen bond with Lys868: 12, 13a, and 14a. This hydrogen bond occurs rarely in VEGFR2-13a complex, while it occurs frequently in the last 5 ns MD of VEGFR2-12 and VEGFR2-14a complexes. As for the hydrogen bonds with Thr916, this hydrogen bond is stable in the complexes formed between VEGFR2 with 15 and 18, relatively stable in the complexes formed between VEGFR2 with 13a, 14a, and 17, while rarely occurs in VEGFR2-16 and VEGFR2-20 complexes.
The hydrogen bonding modes of 13a and 14a are somewhat alike. They also have a stable hydrogen bond with Asp1046 and a less stable hydrogen bond with Cys919. The occurrence of hydrogen bond involving Glu885 for 13a, 19, and 20 is like that for sorafenib. Glu885 forms more stable hydrogen bond in VEGFR2 complexes with 15, 16, and 18 as compared with that in VEGFR2-sorafenib complex.
All of the compounds form more than 2 hydrogen bonds with VEGFR2 except 12 and 21. VEGFR2-21 remains one persistent hydrogen bond (Asp1046) while VEGFRR2-12 remains two persistent hydrogen bonds (Asp1046 and Lys868) during most time of MD. From Table 1 we can see that these two compounds have bad inhibitory activity (the inhibitory activities of 21 and 12 are 5.5 μm and 0.40 μm, resp.).
3.4. Energy Decomposition
Table 1 lists the interaction energies of all the VEGFR2 complexes calculated using the last 2 ns MD outcomes. From Table 1 we could see that interaction energies of part A for compounds 12–22 are all lower than sorafenib, suggesting that replacement of N-methylcarbamoyl-4-pyridine to 6,7-dimethoxyquinazoline strengthens the interaction with VEGFR2. From Supplementary Table 1 we can see that Coulomb interaction contributes less to the whole interaction energy between VEGFR2 and the title compounds. The interaction is governed by van der Waals interactions, which is in accordance with the results of Muñoz et al. . Parts A and B have little difference due to same structure for compounds 12–22. The interaction energies of part C for compounds 12 and 22 are −47.82 and −47.69 kJ·mol−1, respectively, which are obviously higher than other compounds. We think that their inhibitory activity against VEGFR2 might be poorer than others. In other words, this indicates that when urea was substituted by an alkyl group (butyl and cyclohexyl), its inhibitory activity decreased, which is in accordance with the experimental results . From Figure 7 we can see that there exists a good correlation () between the calculated interaction energy of part C with experimental pIC50 values. An exception is compound 21. Compound 21 has as lower interaction energy as compound 20. However, 2-naphthyl substituent (20) has an excellent inhibitory activity ( μm) while 1-naphthyl substituent (21) has a poor inhibitory activity ( μm). This may be due to poor hydrogen bond type with VEGFR2 of 21; it only creates one hydrogen bond with Asp1046 of VEGFR2 (see supplementary Figure 3). Thus we can see that a good inhibitor should not only have lower interaction energy with VEGFR2 but also form at least 2 hydrogen bonds with VEGFR2.
Docking and molecular dynamics simulations of VEGFR2-sorafenib complex show that the possible mechanism that sorafenib inhibits VEGFR2 is by forming hydrogen bonds with residues Asp1046, Cys919, and Glu885 of the VEGFR2 binding pocket as well as by hydrophobic contacts with a set of hydrophobic residues. The aforementioned interaction affects the activity of such inhibitors to inhibit VEGFR2. The interaction energy between compounds 12–22 and VEGFR2 has been computed as three parts to consider the effect of different substituent group. Because part A and B are the same for compounds 12–22, part C may play a key role in the interaction with VEGFR2 and the calculated results confirm this. The interaction energy of part C with VEGFR2 and the inhibitory activity do have a good correlation. The obtained results may be helpful in the structural modification and design of VEGFR2 inhibitors, aiding the design and synthesis of new VEGFR2 inhibitors in future.
MD simulations were performed on TianHe-1A supercomputer of National Supercomputing Center in Tianjin. The project was supported by the National Natural Science Foundation of China (21103125) and National Key Project for Innovative Drug (2011ZX09401-009).
Supplementary Figure 1: Hydrophobic surface of VEGFR2-sorafenib.
Supplementary Figure 2: Bond angle variations along MD simulation.
Supplementary Figure 3: Hydrogen bond existence map of compounds 12-22.
Supplementary Table 1: Structure of compounds 12-22 and their VEGFR-2 inhibitory, computed interaction energies.
C. Muñoz, F. Adasme, J. H. Alzate-Morales, A. Vergara-Jaque, T. Kniess, and J. Caballero, “Study of differences in the VEGFR2 inhibitory activities between semaxanib and SU5205 using 3D-QSAR, docking, and molecular dynamics simulations,” Journal of Molecular Graphics and Modelling, vol. 32, pp. 39–48, 2012.View at: Publisher Site | Google Scholar
X. Chen, X. X. Liu, H. Huang, H. H. Hu, and F. C. Jiang, “Construction of pharmaeophore model of EGFR TK inhibitor,” Acta Physico-Chimica Sinica, vol. 24, pp. 281–288, 2008.View at: Google Scholar
F. Luo, J. Gao, Y. H. Cheng, W. Cui, and M. J. Ji, “Interaction mechanisms of inhibitors of glucoamylase by molecular dynamics simulations and free energy calculations,” Acta Physico-Chimica Sinica, vol. 28, no. 9, pp. 2191–2201, 2012.View at: Google Scholar
Glide, Version 5.5, Schrödinger, LLC, New York, NY, USA, 2009.
H. J. C. Berendsen, D. van der Spoel, and R. van Drunen, “GROMACS: a message-passing parallel molecular dynamics implementation,” Computer Physics Communications, vol. 91, no. 1–3, pp. 43–56, 1995.View at: Google Scholar
MacroModel, Version 9.7, Schrödinger, LLC, New York, NY, USA, 2009.