New inhibitors of VEGFR-2 targeting the extracellular domain dimerization process.

We are reporting the discovery of small molecule inhibitors for vascular endothelial growth factor receptor type 2 (VEGFR-2) extracellular domain. The VEGFR-2 extracellular domain is responsible for the homo-dimerization process, which has been recently reported as a main step in VEGFR signal transduction cascade. This cascade is essential for the vascularization and survival of most types of cancers. Two main design strategies were used; Molecular docking-based Virtual Screening and Fragment Based Design (FBD). A virtual library of drug like compounds was screened using a cascade of docking techniques in order to discover an inhibitor that binds to this new binding site. Rapid docking methodology was used first to filter the large number of compounds followed by more accurate and slow ones. Fragment based molecular design was adopted afterwards due to unsatisfactory results of screening process. Screening and design process resulted in a group of inhibitors with superior binding energies exceeding that of the natural substrate. Molecular dynamics simulation was used to test the stability of binding of these inhibitors and finally the drug ability of these compounds was assisted using Lipinski rule of five. By this way the designed compounds have shown to possess high pharmacologic potential as novel anticancer agents.

Background

Neo-angiogenesis has a crucial role in the progression and survival of most types of cancer beside some other proliferative diseases [1]. Vascular endothelial growth factors (VEGFs) are the most important angiogenesis regulators [2]. Recognition of the crucial role of VEGF pathway in the regulation of angiogenesis has led to the development of VEGF-targeted therapy for the development of selective and safe anticancer agents. Several strategies were taken in consideration in targeting the VEGF signaling pathway. Neutralization of the VEGF or VEGFR using antibodies was investigated and has been reported for treatment of metastatic colorectal cancer [3], non-small cell lung cancer [4] and metastatic breast cancer [5]. Targeting VEGFR using small molecule inhibitors was another strategy as well. The small molecule VEGFR-2 kinase inhibitors sorafenib [6] and sunitinib [7] have been approved by the FDA for the treatment of advanced renal cell carcinoma [8]. A number of other small-molecule tyrosine kinase inhibitors are currently under investigation in Phase III clinical trials. Another approach includes the use of Peptoid ligands for the extracellular domain of VEGFR5 [9]. VEGFR-2 is considered to be the major receptor responsible for mediating physiological and pathological effects of VEGF- A on endothelial cells [10]. Although VEGFR-1 has a tenfold higher binding affinity to VEGF, it exerts less activation of intracellular signaling intermediates than VEGFR2 and consequently it can act as a negative regulator of angiogenesis by binding VEGF and preventing its binding to VEGFR-2 [2]. After the Determination of the crystal structure of the most membrane proximal domain [11], it was revealed that homotypic contacts between them are essential for ligand-induced receptor activation and cell signaling. These contacts were shown to be mediated by salt bridges and van der Waals contacts formed between Arg726 of one protomer and Asp731 of the other protomer which was proved by the fact that ligand-induced auto-phosphorylation and cell signaling via VEGFR1 or VEGFR2 harboring mutations in critical residues (Arg726 or Asp731) are impaired. The newly discovered homotypic contact region of VEGFR2 provides a novel target for pharmacological intervention of pathological activation of VEGFR2 which is essential for cancer development. In this study, we have used the newly discovered crystal structure of the extracellular domain of the VEGFR for the purpose of developing a small molecule anti-angiogenic drug candidate as anticancer agent through inhibition of the VEGFR dimerization process. For this purpose we have used molecular docking techniques for virtual high throughput screen (v-HTS) for an inhibitor of this process [12]. Additionally, fragment based design (FBD) was also employed in contrast as an alternate strategy in developing molecular designs of relatively higher binding affinities [13].

Materials and Methodology

The 3D X-ray crystal structure of monomeric and homodimeric forms of membrane-proximal Ig-like domain of the ectodomain (D7) of VEGFR-2 [PDB ID: 3KVQ] was exploited as target for vitrual High-throughput screening (v-HTS) and fragment-based design (FBD). Both the proteins were prepared, undergone rigid and flexible docking, energetically probed, and subjected to molecular dynamics (MD) simulations. All computational analysis was carried out using Discovery Studio (DS) 2.5 (Accelrys Software Inc., San Diego; http://www.accelrys.com). The hardware used was the workstation cluster of pharmaceutical chemistry department, at Ain Shams University, Faculty of Pharmacy.

Receptor/Ligands preparation

The target domain monomeric form as well as the homodimerized complex (biological assembly coordinates file) were obtained, structure cleaned, the hydrogen atoms were added, typing was carried out by all-atom CHARMm forcefield (version 35b1) (Momany-Rone partial charges method) [14] and followed by Smart Minimizer algorithm – an algorithm that performs 1000 steps of Steepest Descent with a root mean square (RMS) gradient tolerance of 3, followed by Conjugate Gradient minimization – until the RMS gradient for potential energy was less than 0.05 kcal /mol/Å. Using the ‘Binding Site’ tools available in DS 2.5 the binding site was defined to be the homotypic contact region inscribed inside a sphere of 13.5 Å radius. Ligands for virtual screening were obtained from Zinc database ( http://zinc.docking.org) from commercially available lead-like compounds, filtered and then prepared using the ‘Prepare Ligands’ protocol in order to standardize charges, enumerate ionization states and generate tautomers at physiological pH range (the eventual count of the library reached around 110,000 molecules). The ligands were typed similarly by CHARMm for partial charges set up. The last two steps were additionally carried out to the ‘De Novo’ designed ligands prior to their flexible docking.

High-throughput docking

Rigid molecular docking is the computational method used to predict the binding of the ligand to the receptor binding site by varying position and conformation of the ligand keeping the receptor fixed. Initially LibDock, [15] a relatively fast algorithm that conducts ‘HotSpots’ matching of ligand conformations with rigid binding site's HotSpots map that is well-suited for largesized libraries, was used to filter the obtained library according to their binding capacity into the proposed active site in order to decrease the library size. With specified 200 hotspots and 3 saved for each ligand with conformational sampling using the ‘BEST’ algorithm. By the utilization of CDocker [16], the protocol that employs a CHARMm-based MD scheme, the ligands were docked into receptor binding site, to further narrow down the number of top hits on a more precise scale. Thus, refining each hit through simulated annealing between 300 and 700 K. Finally the top ligands were selected for flexible docking by using the corresponding protocol in DS 2.5 that allows for receptor flexibility [17]. The flexible residues were determined to be those of ASP710, SER711, ILE713, ILE724, ARG725, ARG726, VAL727, ARG728 and ASP731, which are thought to be the critically interactive residues during the dimerisation process, and post-docking ChiRotor and ligandannealing refinement were allowed.

De Novo design

A set of novel ligands were designed by the help of the ‘De Novo Evolution’ protocol, that constructs larger molecules in the binding site of the receptor starting from an initial scaffold molecule employing a genetic algorithm. The core scaffold was manually chosen and further manual modifications were carried out on the resulting structures. The protocol was run in ‘Full Evolution’ mode starting with a population of 40 at each of the 6 generations with the selection of 15 survivors per generation using the top LUDI 3 [18] – an empirical function that provides a fast and accurate score describing the binding affinity – for scoring structures. Designed ligands were first heated in the receptor site up to 300 K and trajectories analyzed to assure the mode of binding at physiological temperature for 5 fs with Generalized Born Simple Switching (GBSW) implicit solvation conditions – the methodology is well suited with molecular dynamics trajectory [19], in order to initially filter the candidate structures.

Scoring and Energetics

The ligand-receptor complexes of selected clustered poses from each flexibly docked ligand, as well as the designed ligands were further subject to a minimization process using the ‘Minimization’ protocol in DS 2.5 through the ‘Smart Minimizer’ – steepest-descent followed by conjugate gradient method – until potential energy RMS gradient was less than 0.05 kcal/mol/Å, with generalized born simple switching implicit solvation at salt concentration of 0.145 M to better mimic physiological conditions [20]. The same was done to the homodimerised complex for acquiring an energy estimate of the low-affinity D7 dimerization. The last followed binding energy calculation by the ‘Calculate Binding Energy’ protocol yet through the more robust Poisson-Boltzmann solver. The binding energies were used for scoring due to the lack of prior arts correlating scoring functions reproducibility on the current binding site.

Molecular Dynamics

The ‘Standard Dynamics Cascade’ Protocol [21] of DS 2.5 was implemented with a time step of 1 fs. The initial minimization of the implicitly solvated protein–ligand complex was carried out in two steps on residues of a selected sphere of a 13 Å radius inscribing the ligand, binding site and its surrounding residues and loops, keeping about two thirds of the protein domain fixed throughout the whole dynamics simulation cascade. Bad contacts were corrected using the steepest-descent algorithm without major distortion in the structure down to 0.01 kcal/mol/Å RMS gradient. Further minimization was carried out using the conjugate gradient algorithm with a lowenergy starting point down to 0.0001 kcal/mol/Å final RMS gradient. Initially the temperature of the system was raised from 50 to 330 K for 5000 heating steps at each temperature. The system was then equilibrated for 200 ps, and finally the production run was carried out for another 1500 ps.

Results and Discussion

The detailed investigation conducted by Yang et al., [11] concluded the critical role of the dimerization process of D7 through the homotypic contacts in ligand-induced (VEGF-A) autophosphorylation of VEGFR-2 and that latter's eventual activation. Such interaction was regarded to be of low affinity, the fact proven by analytical centrifugation that determined the dissociation constant for (Kd) dimerisation of the isolated D7 region was found to exceed 10-4 M. Relating the total binding energy to the dissociation constant, the following formula was used:

ΔG(lgd-ptn)•ln Kd(ptn-ptn)/ΔG(ptn-ptn) = ln Kd(lgd-ptn) Where ΔG is the binding free energy, and Kd is the dissociation constant. The reference value was calculated from the molar dimerisation free energy of D7 that approximately corresponds to a Kd of 100 µM. The above formula was then used to calculate dissociation constants of the ligand-protein complexes. Eventually, showing potentially efficacious leads resulting from the design procedure as opposed to the virtual screening [11].

Virtual High-throughput Screening

The aforementioned selected library of commercially available compounds (approximately 110,000 compounds) was initially docked with LibDock into the binding pocket (Figure 1) to reduce its size down to approximately 35,000 poses − with a maximum of three poses per molecule, which was appropriate for the more robust CDOCKER protocol, that was allowed to perform final annealing, has shrunk down the size into 751 poses. Based on the ‘CDOCKER Energy’ score which includes internal ligand strain energy and receptor-ligand interaction energy, the 19 non-redundant structures of the top 50 poses were chosen for an additional flexible docking simulation in order to better investigate their mode of binding. The resulting complexes output by the latter process were retrospectively investigated, and five poses for each ligand-flexible receptor poses were selected as candidates for the energetics simulation. The purpose of the re-ranking the compounds based on thorough energetic probing was to better account solvation and entropy-related factors that could strongly affect the binding efficacy. So, the complexes were minimized with implicit solvation and then the total binding energies (regarding receptor and ligand entropy) were calculated to select the best ligand. Fifteen out of the nineteen ligands have shown negative total binding energies as shown in (Table 1, see Table 1). The results show a range from very poor to modest binding affinities. Hence, FBD is capable of generating more robust solutions provided chemical compatibility and synthetic feasibility are addressed.

Figure 1

Interpolated charge surface of the binding pocket residues.

De Novo Design of Small Molecule Inhibitors

Firstly, we have selected two pharmacophoric groups as scaffolds from the designed prototype D3 that have been manually placed into the active site and then minimized for the Ludi algorithm to start with in ‘Full Evolution’ mode, namely, guanidinium group at the acidic pocket of the active site − by natural binding mode mimicry − and tetrazolyl sandwiched between the two guanidinium groups of the ARG725 and ARG726. The standard Ludi fragment libraries − ‘Link’ and ‘Receptor’ − developed by Böhm [22] were used. The 6th generation has proved to supply optimally sized candidate ligands with sufficiently high LUDI scores. The selected ligands were those of the acceptable synthetic feasibility. Later, manual optimization for binding enrichment was carried out, mainly through replacement of weakly interacting or misplaced moieties and addition of electron-rich annular systems in the vicinity of the ionisable guanidinium groups of the active site in aim to the formation of multiple cation-pi interations, the kind of non-covalent bonding that has proved capable of competing with full aqueous solvation as well as baring binding energies beyond ­20 kcal/mol between various protonated amines and activated ring systems [23]. Other condition-specific interactions were also sought during the design process like cyclodione-arginine reversible covalent interaction [24], distributed formal and densely charged groups. Preliminary heating to 300 K was conducted to discard poorly bound ligands, and then the stable complexes were re-cooled to decrease the gradient after relieving steric clashes. Then, the complexes of the five remaining compounds 1D-5D (Figure 2) were minimized and binding energies were calculated (Table 2, see Table 2).

Figure 2

The binding modes of the designed compounds 1D- 5D (from left to right). Hydrogen bonds visualized as dashed black lines; Cation-pi interactions visualized as orange lines.

Molecular dynamics simulation

To account for the effect of solvent on binding stability; GBSW implicit solver was used with the same parameters of prior minimization. The total energy and simulation temperature were found to remain steady with little fluctuation during the production stage time interval (1500 ps), which was preceded by heating followed by equilibration. The conformational sampling of the trajectory at 0, 250, 500, 750, 1000, 1250, and 1500 ps of the production run are shown in Figure 3. Conformations of both the ligand and the protein showed that they stayed conserved with slight perturbation, at the solvent accessible region by the residues of Arg726 and Arg725, especially in the case of the ligands' 2D rotatable tail. Due to the mandatory role played by hydrogen bonds and cation-pi interactions, they were closely monitored and their existence sampled through the trajectory (Table 3, see Table 3). From the table it is clear that the bond showed a high temporal stability throughout the production dynamics phase. Calculating the molecular properties (Table 2, see Table 2) of the designed ligands could predict good pharmacokinetic properties. This should lead to the expectation of a valuable in vivo performance of such agents. Also in most of the cases, Lipiniski's ‘rule of 5’ was satisfied, which suggests a good oral bioavailability of the compounds [25]. Due to the high potentials of these compounds, we are planning to synthesize and report their activities in due course.

Figure 3

Molecular dynamics trajectory for the complexes of the ligands 1D-5D (from left to right). Snapshots of the designed ligands and the side chains of the binding site residue conformers extracted from the production dynamics trajectory at times 0, 250, 500, 750, 1000, 1250, 1500 ps.

Conclusion

A novel homotypic dimerization region that is believed to play a critical role in VEGFR signal transduction was exploited as target in structure-based drug design. A large sized library of commercially available compounds were virtually screened and showed Kd values higher than that of D7 dimerisation. Nevertheless, the de novo design process has resulted in many promising ligands in the sub-micromolar range (binding affinity exceeding 340-fold that of dimerisation). In addition designed compounds could be utilized in construction of a pharmacophore model to screen much larger databases of compounds for possible hits. This study could be a promising in developing lead compounds in the discovery of anticancer drugs.