Molecular features of interaction between VEGFA and anti-angiogenic drugs used in retinal diseases: a computational approach

Anti-angiogenic agents are biological drugs used for treatment of retinal neovascular degenerative diseases. In this study, we aimed at in silico analysis of interaction of vascular endothelial growth factor A (VEGFA), the main mediator of angiogenesis, with binding domains of anti-angiogenic agents used for treatment of retinal diseases, such as ranibizumab, bevacizumab and aflibercept. The analysis of anti-VEGF/VEGFA complexes was carried out by means of protein-protein docking and molecular dynamics (MD) coupled to molecular mechanics-Poisson Boltzmann Surface Area (MM-PBSA) calculation. Molecular dynamics simulation was further analyzed by protein contact networks. Rough energetic evaluation with protein-protein docking scores revealed that aflibercept/VEGFA complex was characterized by electrostatic stabilization, whereas ranibizumab and bevacizumab complexes were stabilized by Van der Waals (VdW) energy term; these results were confirmed by MM-PBSA. Comparison of MM-PBSA predicted energy terms with experimental binding parameters reported in literature indicated that the high association rate (K_on) of aflibercept to VEGFA was consistent with high stabilizing electrostatic energy. On the other hand, the relatively low experimental dissociation rate (K_off) of ranibizumab may be attributed to lower conformational fluctuations of the ranibizumab/VEGFA complex, higher number of contacts and hydrogen bonds in comparison to bevacizumab and aflibercept. Thus, the anti-angiogenic agents have been found to be considerably different both in terms of molecular interactions and stabilizing energy. Characterization of such features can improve the design of novel biological drugs potentially useful in clinical practice.

Introduction

Vascular endothelial growth factor (VEGF) plays a pivotal role in angiogenesis through activation of specific receptors (VEGFR). Among different VEGF isoforms, VEGFA (VEGF₁₆₅) is the isoform that has principally been involved in angiogenesis. Angiogenesis is the physiological process of formation of new vessels from pre-existing ones and is crucial during development or tissue repair, but detrimental in some disease states. Therefore, in some pathological conditions angiogenesis inhibition is a relevant therapeutic goal (Carmeliet and Jain, 2011). Anti-VEGF agents are used to inhibit primary tumor and metastasis growth (Chung and Ferrara, 2011), both in adjuvant, and more recently, in neoadjuvant settings (Welti et al., 2013). Recently, anti-angiogenic agents have been used to treat ocular pathological conditions such as age-related macular degeneration (AMD) and diabetic macular edema (DME) (Holz et al., 2014; Stewart, 2014a). Diabetic retinopathy is a leading cause of vision loss of working-age adults, and DME is the most frequent cause of vision loss related to diabetes. AMD is a progressive multifactorial neurodegenerative disease that impairs the visual field. Currently, three VEGF inhibitors are used to treat retinal disorders characterized by neovessel formation: ranibizumab (Lucentis^®), aflibercept (Eylea^®) and bevacizumab (Avastin^®), this latter is used off-label (Figure 1). Furthermore, several new anti-angiogenic agents are in clinical development (Teicher, 2011; Hanout et al., 2013; Li et al., 2014). Bevacizumab is a humanized monoclonal antibody, ranibizumab is the mutated Fab (Fragment antigen-binding) of the monoclonal antibody (Ab) originating bevacizumab and aflibercept is a fusion protein that works as decoy VEGF receptor (Figure 1). Ranibizumab has been developed for intravitreal injection and shows improved ocular pharmacokinetics (Xu et al., 2013) compared to bevacizumab. X-ray structures of VEGFA bound to ranibizumab (PDB: 1CZ8) (Chen et al., 1999) or Fab-bevacizumab (PDB: 1BJ1) (Muller et al., 1998) have been solved. The three dimensional structure of aflibercept is not available, albeit widely investigated. Aflibercept, also known as VEGF-trap, is a decoy receptor where two binding domains, the domain 2 (d2) of VEGFR1 and the domain 3 (d3) of VEGFR2 (from N-terminus to C-terminus of primary sequence) are connected to the fragment crystallizable region (Fc) of human immunoglobulin (Ig) (Holash et al., 2002). Throughout the text the binding domain of aflibercept is named “VEGFR1d2_R2d3.”

FIGURE 1

Figure 1. Schematic structures of ranibizumab (A), bevacizumab (B), and aflibercept (C). Ab stands for antibody, Fab stands for fragment antigen binding, Fc stands for fragment crystallizable region. R1d2 stands for domain 2 of vascular endothelial growth factor receptor VEGFR1 and R2d3 stands for domain 3 of VEGFR2. Black bars correspond to inter-chain disulfide bridges.

Iyer et al. (2010) reported the x-ray structure of VEGFR1 domain 2, as binary complex with VEGFB, while Leppänen et al. (2010) have solved the x-ray structure of VEGFR2 domain 2 and 3 (VEGFR2d2_d3) in complex with VEGFC, that we used for homology modeling of aflibercept's binding domain (VEGFR1d2_R2d3). However, in silico comparison of anti-VEGF/VEGFA complexes has not yet been carried out. Therefore, the aim of our study was to test the hypothesis that the three anti-VEGF/VEGFA complexes, involving ranibizumab, bevacizumab or aflibercept have different energy terms corresponding to different molecular and/or atomic interactions. To this end we carried out protein-protein docking with the software PyDock (Cheng et al., 2007; Jiménez-García et al., 2013), simulation of complexes with molecular dynamics (MD, GROMACS, Pronk et al., 2013) and molecular mechanics Poisson-Boltzmann surface area calculation (MM-PBSA) by using the g_mmpbsa tool (Kumari et al., 2014). Furthermore, we have added the analysis of the protein contact networks (Vishveshwara et al., 2009; Di Paola et al., 2013) on MD trajectories, to analyze the correlation of key topological properties over the time. Our results show that the three anti-angiogenic agents considerably differ, both in terms of molecular interactions and stabilizing energy; furthermore our in silico data are consistent with published experimental binding parameters (Papadopoulos et al., 2012).

Methods

Molecular Modeling and Protein-protein Docking

Full protein sequences of Fab-bevacizumab binding domain and ranibizumab were retrieved from the DrugBank database (http://www.drugbank.ca, accession numbers: DB00112 and DB01270, respectively). The sequence of aflibercept binding domain (VEGFR1d2_R2d3) was built by connecting the sequences of domain 2 of human VEGFR1 and domain 3 of human VEGFR2:

SDTGRPFVEMYSEIPEIIHMTEGRELVIPCRVTSPNITVTLK KFPLDTLIPDGKRIIWDSRKGFIISNATYKEIGLLTCEATVNG HLYKTNYLTHRQTNTIIDVVLSPSHGIELSVGEKLVLNCTART ELNVGIDFNWEYPSSKHQHKKLVNRDLKTQSGSEMKKFLST LTIDGVTRSDQGLYTCAASSGLMTKKNSTFVRVHEDPIEGR.

Structures of Fab-bevacizumab and ranibizumab binding domains were modeled by means of the SwissModel web server (http://swissmodel.expasy.org/; Schwede et al., 2003) using, respectively, the x-ray structures of bevacizumab (PDB: 1BJ1) and ranibizumab (PDB: 1CZ8) in complex with the dimer of VEGFA as templates. This procedure included the modeling of short loops, which are lacking in the X-ray structures. The structure of VEGFR1d2_R2d3 was modeled using the x-ray structure of VEGFR2 in complex with VEGFC (PDB: 2X1W). Protein-Protein docking was carried out through the pyDockWEB (http://life.bsc.es/servlet/pydock/home/; Jiménez-García et al., 2013). The 3D coordinates of the ligand, the VEGFA dimer, were extracted from the crystal structure of VEGFA bound to ranibizumab (1CZ8). Homology models of Fab-bevacizumab and ranibizumab were set as PDB files of receptors. The complex VEGFR1d2_R2d3/VEGFA was built by using two randomly selected frames, belonging to relative conformational minimum of a preliminary MD of VEGFR1d2_R2d3 (Supplementary Material). Because of the low RMSD (0.06 nm), no differences were found between the two predicted complexes of the two analyzed frames. In order to evaluate the effectiveness of MD on VEGFR1d2_R2d3 structural optimization, the protein-protein docking of VEGFR1d2_R2d3 was also carried out using its starting model.

Given the 3D coordinates of two interacting proteins, pyDockWEB uses FTDock for generation of the predicted complexes by rigid-body docking (Gabb et al., 1997). No restrains have been applied in this stage. The predicted complexes were then evaluated by PyDock scoring function, which includes electrostatics, desolvation energy and limited Van der Waals (VdW) contributions (Cheng et al., 2007).

Molecular Dynamics Simulation (MD)

All-atom MD has been carried out with GROMACS 4.6 (Pronk et al., 2013), taking advantage of graphic processing unit (GPU) acceleration. Disulfide bridges in the object molecules were built as follows: VEGFA (Cys26-68Cys, Cys102-52Cys, Cys102-57Cys) in each monomer; Fab-bevacizumab and ranibizumab (Cys23-88Cys and Cys134-194Cys in L chain, Cys22-96Cys and Cys150-206Cys in H chain); VEGFR1d2_R2d3 (Cys30-79Cys in R1d2 and Cys189-125Cys in R2d3). Ionization state of residues was assigned at pH 7.4, OPLS-AA (Jorgensen et al., 1996) force field parameters were assigned to protein; thereafter, the following simulation protocol was applied: (i) building of the simulation system (explicit solvation with water molecules, SPC-E model (Berendsen et al., 1987), and neutralization with NaCl 150 mM); (ii) 2000 steps of steepest descent energy minimization (time-step = 1 fs); (iii) 500 ps of NVT equilibration (time-step = 1 fs); (iv) 500 ps of NPT equilibration (time-step = 1 fs); v. production runs (time-step = 2 fs). In order to use a 2 fs time-step we applied LINCS constraints (Hess et al., 1997) through the P-LINCS algorithm (Hess, 2008). Long-term electrostatic interactions were treated by applying the Particle Mesh Ewald (PME) method (Cheatham et al., 1995) (pme_order = 4; fourier spacing 0.16). Temperature was kept constant at 300 K with Berendsen thermostat applied to two groups (protein and non-protein), pressure was controlled by the Parinello-Rahman barostat (P = 1 bar; water compressibility 4.5 E-5 bar⁻¹). Trajectory and energy outputs were collected every 10 ps. Three independent replicas for each MD were carried out: VEGFR1d2_R2d3 (40 ns), ranibizumab (20 ns), Fab-bevacizumab (20 ns), VEGFR1d2_R2d3/VEGFA complex (20 ns), ranibizumab/VEGFA complex (20 ns), Fab-bevacizumab/VEGFA complex (40 ns). VEGFA was simulated for just 3 ns, because the RMSD, secondary structure and cosine content indicated an equilibrated simulation. MD runs were launched on the EURORA machine (HPC-CINECA) using an average of 64 CPUs, 8 GPUs and 8 message passing interface (MPI) processes per run. Visual analysis of trajectories was carried out with Visual Molecular Dynamics VMD 1.9 (Humphrey et al., 1996); whereas analysis of MD simulations was carried out with the analysis package of GROMACS 4.6 (Pronk et al., 2013). We carried out Principal Component Analysis (PCA) of trajectories for characterization of principal collective motions of the analyzed molecules; such motions can identify the movement of a domain toward another domain of a given protein (Amadei et al., 1993). Each principal component is also named eigenvector to which an eigenvalue (~frequency) is associated. Eigenvectors were projected in the MD trajectories in order to visualize the collective motions of the protein, thereby reducing the dimensionality of trajectories. Because PCs with high cosine content are indicative of random diffusive movement (Hess, 2002; Neugebauer et al., 2002), in order to assess the correct conformational sampling of trajectories, cosine content of the first two PCs has been determined. A cosine content lower than 0.5 is considered as indicative of satisfactory conformational sampling. The number of contacts at protein-protein interface have been determined with g_mindist tool of GROMACS, residues were defined as interacting at 3.5 Å. High frequency H-bonds at the protein-protein interface have been analyzed using the HBonanza open source software (Feenstra et al., 1999). RMSD (Figures 2, 4) calculations have been carried out by selecting Cα carbons of proteins and excluding the analysis of translational and rotational movements. The splitting of RMSD (Figure 4), of VEGFA/anti-VEGF trajectories, was obtained as follows: we launched the command g_rms using the gromacs index.ndx file in order to carry out the calculation on each component of complexes, the VEGFA and the anti-VEGFA, respectively; therefore RMSD profiles of complexes in Figure 2 are not the sum of split RMSD in Figure 4.

FIGURE 2

Figure 2. Root-mean-square deviation (RMSD) profile of simulated molecules. (A) ranibizumab; (B) Fab-Bevacizumab; (C) VEGFR1d2_R2d3; (D) ranibizumab/VEGFA complex; (E) Fab-bevacizumab/VEGFA complex; (F) VEGFR1d2_R2d3/VEGFA complex; (G) VEGFA replica 1; (H) VEGFA replica 2; (I) VEGFA replica 3. Black line corresponds to the first replica, red line corresponds to the second replica, green line corresponds to the third replica.

Protein Contact Networks

The analysis was applied to 100 frames, 10 ps-spaced, taken from the last 10 ns of each MD replica; whereas MD of VEGFA was analyzed in the whole 3ns-replicas. Each frame was converted into an indirect, unweighted graph, whose nodes are the α-carbons and edges represent the mutual spatial distances between residues when their distance is within 4 and 8 Å, accounting to VdW interactions. This method, applied in a previous work to highlight the allosteric character of protein-ligand complexes (De Ruvo et al., 2012), accounts only for non-covalent bonds. After the definition of the protein contact network, the following global (whole structure) topological descriptors and a chemical-physical descriptor were derived:

(1) adeg (average degree): the average number of contacts involving a single residue;

(2) asp (average shortest path): the shortest path between two residues indicates the minimum number of steps (links) from one residue to another; asp is the average value over all residue pairs;

(3) E (the Graph Energy): couples the graph global connectivity to the interactions in the represented molecular structure (Balakrishnan, 2004);

(4) dG_solv (free solvation energy): quantitative descriptor of protein stability in water (Eisenberg and McLachlan, 1986), takes into account the overall energy gain of atoms passing from protein to water.

In order to assess whether MD reached a relative conformational minimum, the correlation pattern of topological descriptors and energy (dG_solv) was analyzed over the time; the conformational minimum is characterized by a non-significant correlation coefficient of the variables. Furthermore, the protein contact network has been partitioned into clusters, according to a spectral clustering algorithm, which was previously applied to split the protein structure into functional modules. This method is able to detect functional domains in protein structures and complexes, along to the topological role of single residues that account for inter- and intra-module interaction (Tasdighian et al., 2014). Clustering results for each average structure of complexes and protein systems is represented as a partition color map of a two-dimensional matrix of cluster distribution along the sequence. Background (blue) corresponds to residues that do not belong to the same cluster. Residues belonging to the same cluster are represented with the same color. An interruption between cluster-sequence continuity, i.e., a residue shifting to a different cluster, which corresponds to a long-range contact, is represented as a projection termed “whisker.”

MM-PBSA and Energy Decomposition

The MM-PBSA method calculates the three energetic terms of the binding free energy (Equation 1) (Kollman et al., 2000):

$\begin{array}{cc}\Delta {\text{G}}_{\text{binding}}={\text{G}}_{\text{complex}}-({\text{G}}_{\text{protein}}+{\text{G}}_{\text{ligand}})& \left(1\right)\end{array}$

Free energy of either products or reagents is calculated taking in account three terms (Equation 2):

$\begin{array}{cc}{\text{G}}_{\text{x}}=<\text{E}{}_{\text{MM}}>+<{\text{G}}_{\text{solvation}}>-\text{TS}& \left(2\right)\end{array}$

Where E_MM is the vacuum potential energy and G_solvation is the free energy of solvation. E_MM includes E_bonded and E_non−bonded energies; E_non−bonded energy is the summation of Van der Waals (Lennard-Jones potential function) and electrostatic (Coulomb potential function) energy terms. G_solvation, is characterized by the summation of two terms, G_polar and G_apolar, which represent the electrostatic and the non-electrostatic term. G_polaris calculated using a continuum implicit solvent model applying the Poisson-Boltzmann equation (Baker et al., 2001). The G_apolarterm results from the summation of G_cavity and G_VdW terms; G_cavity is the work done by the solute to create a cavity in the solvent, G_VdW is the attractive Van der Waals energy between solvent and solute. G_apolar accounts for the hydrophobic effect (Richmond, 1984). Single trajectory MM-PBSA calculations were carried out on each of the three MD replicas of complexes by using the g_mmpbsa tool (Kumari et al., 2014), which integrates functions from GROMACS and APBS (http://rashmikumari.github.io/g_mmpbsa/). The dielectric relative constant ε has been set to 3 for protein and 80 for water (Kukic et al., 2013). The solvent accessible surface area (SASA) method was used for calculation of G_apolar; the surface tension constant γ was set to 0.022 KJ/mol Ų (Nicholls et al., 1991). The current implementation of the MM-PBSA method in g_mmpbsa does not include calculation of the entropic term (S) in the equation 2; indeed, g_mmpbsa is unable to provide prediction of absolute binding free energy, providing mainly relative binding energies. For this reason, we use through the text the notation “ΔE_binding” instead of “ΔG_binding,” this latter would include entropy. The g_mmpbsa tool predicts the contribution of residues to the binding free energy by means of energy decomposition calculations; it allows the visualization of such results by means of common molecular visualization tools. The energy contribution is written in the b-factor field of the.pdb file and can be mapped in a 3D structure. MM-PBSA was applied to 100 frames, 10 ps-spaced, taken from the last 10 ns of each MD replica. The hardware used for these calculations was a Desktop PC (12 core Intel i7, 64 GB RAM, two GeForce GTX 680-SLI) launching 16 MPI processes per job with maximum available performance (8 h per calculation, 2% load memory).

Statistical and Graphical Analysis

GraphPad (version 6; San Diego, CA, USA) was used to carry out statistical analysis and graph creation. Comparisons between two independent groups were made by unpaired Student's t-test; p < 0.05 were considered significant. RMSD graphs have been created with xmgrace (open GNU license). Figures have been created with OPEN PyMOL Molecular Graphics System, Schrödinger, LLC (New York, NY, USA).

Results

Homology Modeling and Protein-protein Docking

Modeled structures of Fab-bevacizumab and ranibizumab showed low root mean square deviation (RMSD), respectively 0.003 nm and 0.002 nm, upon superimposition on PDB: 1BJ1 and PDB: 1CZ8 x-ray structures, respectively. The model of VEGFR1d2_R2d3 showed low RMSD (0.04 nm) upon superimposition on the corresponding template PDB:2X1W. VEGFR1d2_R2d3 was subjected to short all-atom MD simulation prior protein-protein docking with VEGFA (Supplementary Material, Figure S1). Protein-protein docking predictions were carried out with PyDock. The software was first validated by building complexes of Fab-bevacizumab and ranibizumab with VEGFA. RMSDs between the best scored Fab-bevacizumab/VEGFA and ranibizumab/VEGFA complexes and the correspondent x-ray structures were negligible (0.046 and 0.045 nm respectively). Subsequently we modeled the VEGFR1d2_R2d3/VEGFA complex. When compared to the starting model, MD simulation gave a better docking score for VEGFR1d2_R2d3 (Table 1). The complex VEGFR1d2_R2d3/VEGFA was compared to X-ray structures of VEGFR2 bound to VEGFC and VEGFA (PDB: 2X1W and PDB: 3V2A, respectively; Supplementary Material, Figure S2). Rough energetic evaluation of predicted complexes, obtained with PyDock, is shown in Table 1. Notice that VEGFR1d2_R2d3/VEGFA was stabilized by electrostatic interaction energy compared to Fab-bevacizumab/VEGFA and ranibizumab/VEGFA complexes, which were rather characterized by stabilizing desolvation and VdW energy terms.

TABLE 1

Table 1. Energetic contribution to PyDock score.

Molecular Dynamics Simulation (MD)

Three independent MD replicas of VEGFA, ranibizumab, Fab-bevacizumab, VEGFR1d2_R2d3 and of the corresponding 1:1 complexes with VEGFA were carried out. Ranibizumab and Fab-bavacizumab reached a relative minimum within 10 ns (Figures 2A,B), VEGFR1d2_R2d3 reached a relative minimum within 40 ns (Figure 2C). The great conformational fluctuation of unbound VEGFR1d2_R2d3 is mainly related to rotational freedom of the connecting hinge between domains R1d2 and R2d3, and to inter-conversion of turns to coil and vice-versa, of the loops in R1d2 domain (see also Supplementary Material, Figures S3–S5). Fab-bevacizumab/VEGFA complex was characterized in all three replicas by an RMSD higher than the other complexes (Figure 2E), despite each replica reached a relative minimum after about 10 ns simulation, as did the other complexes (Figure 2); VEGFA reached a relative conformational equilibrium in 3 ns (Figures 2G–I).

The secondary structures of unbound VEGR1d2_R2d3 and Fab-bevacizumab/VEGFA did not significantly change over the time (Supplementary Material, Figures S3–S8).

In order to characterize the principal collective motions of proteins and their respective macromolecular complexes we carried out PCA of covariance matrix of trajectories. The first six eigenvectors explained about 99% of variance of each simulation. The first principal component (PC1) was projected into the MD trajectory of each simulated molecule (videos in Supplementary Material). In order to verify the correct conformational sampling of MD cosine content of the first two eigenvectors of each trajectory was analyzed (Table 2). Because cosine content of PC1 and PC2 was lower than 0.4 (cut-off 0.5), we assumed that conformational sampling of all MD runs was satisfactory.

TABLE 2

Table 2. Cosine content of the first two PCs of molecular dynamics simulations.

Protein Contact Networks Analysis

To confirm that the systems were analyzed at relative conformational minimum, topological descriptors of protein contact network were also analyzed. The results of correlation analysis between topological descriptors and time for complexes are reported in Table 3 (for unbound systems see also Supplementary Material, Table S1). The energetic descriptor dG_solv did not change over the time, such that the correlation dG_solv vs. time was not significant.

TABLE 3

Table 3. Correlation analysis of topological parameters of anti-vegf/VEGFA complexes.

In contrast, some topological descriptors showed correlations with time and/or each other. Indeed, the average shortest path asp positively correlated with time (0.61), while the average degree adeg and the graph energy E showed, respectively, a positive and a negative correlation in ranibizumab complex. The topological descriptors adeg and asp showed negative correlations (-0.60 and -0.54, respectively) with time in VEGFR1d2_R2d3. Finally, strong correlation between E and adeg for all complexes indicated that the graph energy E describes the overall connectivity energy.

After partitioning of the protein contact network into clusters the structure of the complexes was represented as functional modules. Figure 3 reports the clustering of VEGFA (partitioned in two clusters) and all complexes (partitioned in four clusters). Clustering of unbound anti-VEGFs is reported in Supplemental Material (Figure S9).

FIGURE 3

Figure 3. Clustering of protein contact networks. Distribution of clusters along sequence in a matricial space. Partition color maps of (A). VEGFA (200 nodes), (B). VEGFR1d2_R2d3/VEGFA, (C). Fab-bevacizumab/VEGFA, (D). ranibizumab/VEGFA. The first 200 nodes in complexes corresponds to bound VEGFA. In (A) green and light blue clusters corresponds to VEGFA. In (B–D) dark red and orange are clusters of anti-VEGF. Residues (nodes) belonging to the same cluster have the same color, long projections “whiskers” in the map represent residues shifting to a different cluster with respect to that of neighbors in sequence. The background is dark blue and characterizes residues that are not in the same cluster.

The VEGFA partition in two clusters revealed that they were intermingled (Figure 3A). VEGFR1d2_R2d3/VEGFA partition in four clusters (Figure 3B) revealed a conserved network for VEGFA and two distinct domains corresponding to R1d2 and R2d3, with a lot of long-range interactions with VEGFA (Figure 3B). The partition in four clusters of Fab-bevacizumab/VEGFA (Figure 3C) and ranibizumab/VEGFA (Figure 3D), revealed for bound VEGFA a partition that differ from that of unbound VEGFA, while some whiskers projected from VEGFA to Fab-bevacizumab and ranibizumab modules; furthermore, a greater number of long-range interactions was formed between ranibizumab and VEGFA in comparison to Fab-bevacizumab.

MM-PBSA Calculation

An in-depth analysis of MD data was carried out with MM-PBSA calculations. GROMACS output files were directly analyzed with the g_mmpbsa tool. The estimated binding energy ΔE_binding and its contributing terms were compared to experimental binding and kinetic parameters (Table 4; Papadopoulos et al., 2012).

TABLE 4

Table 4. MM-PBSA results compared to experimental binding parameters.

The relationship between K_D and binding energy is given by: ΔG = RT ln K_D at 1M concentration, where R is the ideal gas constant and T is the absolute temperature. In Table 4 ΔE_binding is reported instead of ΔG_binding, because the entropic term is not included. The comparison of experimental K_D with predicted ΔE_binding confirmed a most favorable binding energy for VEGFR1d2_R2d3 compared to ranibizumab and Fab-bevacizumab bound to VEGFA. The differences in ΔE_binding for ranibizumab/VEGFA and Fab-bevacizumab/VEGFA vs. VEGFR1d2_R2d3/VEGFA were significant (t-test respectively p = 0.003 and p = 0.004). Because ΔE_binding results as summation of different energy terms, we analyzed each single energy term; interestingly, there was a correlation between experimental K_D and ΔG_Apolar (apolar contribution to desolvation energy), suggesting that the hydrophobic effect substantially accounts for affinity. As reported in Table 1, a relevant electrostatic stabilization was predicted by PyDock for the VEGFR1d2_R2d3/VEGFA complex, a data confirmed by MM-PBSA (Table 4). Furthermore, the higher K_on of aflibercept (410 M⁻¹ s⁻¹), as compared to ranibizumab (1.6 M⁻¹ s⁻¹) and bevacizumab (5.6 M⁻¹ s⁻¹), was consistent with the favorable electrostatic component of the binding energy, because the association rate of proteins is known to be related to electrostatic forces (Wade et al., 1998; Zhou, 2001; Pan et al., 2013). Polar contribution to the solvation energy (and to the whole ΔE_binding) is positive, i.e., unfavorable, because of the polar/charged residue transition to a more hydrophobic environment; however, despite the ΔG_polar is more positive for VEGFR1d2_R2d3/VEGFA than for ranibizumab/VEGFA and Fab-bevacizumab/VEGFA, ΔE_{electrostatic}-ΔG_polar gives a favorable gain only forVEGFR1d2_R2d3/VEGFA. Compensation of favorable electrostatic and unfavorable polar desolvation energy terms is commonly found in complex formation, while most stabilization arises from non-polar interactions and hydrophobic effect (Ozboyaci et al., 2011; Spiliotopoulos et al., 2012; Kar et al., 2013). Another information provided by the PyDock prediction was the substantial favorable VdW energy term for ranibizumab/VEGFA and Fab-bevacizumab/VEGFA (Table 1); this observation was confirmed by MM-PBSA (Table 4). Furthermore, VdW energy did not appear to be correlated to any kinetic binding parameter.

Dissociation Rate of Complexes

A lower dissociation rate is reported for ranibizumab/VEGFA compared to the other two complexes (Papadopoulos et al., 2012). The number of contacts at protein-protein interface of complexes at 3.5 Å were determined by the g_mindist tool of GROMACS, while the number of H-bond was assessed by Hbonanza. The number of contacts resulted as follows: Ranibizumab/VEGFA, 480.7 ± 0.5; Fab-bevacizumab/VEGFA, 436.5 ± 0.4; VEGFR1d2_R2d3/VEGFA, 289.9 ± 1.8. The number of H-bonds, whose location is reported in Table 5, were: Ranibizumab/VEGFA, 10; Fab-bevacizumab/VEGFA, 5; VEGFR1d2_R2d3/VEGFA, 4. This analysis suggested that the complex Ranibizumab/VEGFA might be more stable than the other two complexes, in terms of residency time. To test this hypothesis we further analyzed the profiles of complexes by splitting the RMSD for each interacting protein. As shown in Figure 4, ranibizumab in the complex ranibizumab/VEGFA had the lowest RMSD, Fab-bevacizumab in the complex Fab-bevacizumab/VEGFA had an intermediate RMSD, VEGFR1d2_R2d3 in the complex VEGFR1d2_R2d3/VEGFA had the highest RMSD. Further analysis was carried out by calculating the residue-based root mean square fluctuation (RMSF) over the whole simulations. The residue-based RMSF of complexes was visualized into 3D structures (Figure 5).

TABLE 5

Table 5. High frequency H-bonds at anti-VEGF/VEGFA interfaces.

FIGURE 4

Figure 4. Representative graphs of split RMSD of VEGFA (black line) and anti-VEGF (red line) binding domains. (A) Split RMSD for the ranibizumab/VEGFA complex; (B) split RMSD for Fab-bevacizumab/VEGFA complex; (C) split RMSD for VEGFR1d2_R2d3/VEGFA complex.

FIGURE 5

Figure 5. Residue based-root mean square fluctuations (RMSF) of complexes. RMSF increases from blue to red color. Bound VEGFA molecule is highlighted with a red circle. (A) VEGFR1d2_R2d3/VEGFA. (B,C) Fab-bevacizumab/VEGFA, respectively from top and side view. (D,E) Ranibizumab/VEGFA respectively, from top and side view.

The visualization of RMSF confirmed less structural fluctuation of ranibizumab/VEGFA compared to Fab-bevacizumab/VEGFA, consistent with their difference in experimental K_off, and RMSD profiles. The representation of residue-based RMSF of VEGFR1d2_R2d3/VEGFA showed that VEGFA is mainly stabilized at the contact surface with domain 2 and domain 3 (Figure 5A); a high degree of fluctuation, however, detectable out of the contact region, may account for the conformational flexibility of VEGFR1d2_R2d3. Ranibizumab/VEGFA showed less conformational flexibility compared to both VEGFR1d2_R2d3/VEGFA and Fab-bevacizumab/VEGFA, suggesting a higher conformational stability of the complex (Zheng et al., 2006). The lower K_off of ranibizumab/VEGFA is in accordance with the dependency of the residence time τ (1/K_off) on the conformational stabilization of the complex (Copeland, 2011).

Energy Decomposition

Anti-VEGF/VEGFA complexes were further analyzed by looking at residues that favorably contribute to ΔE_binding, by means of energy decomposition calculation. Results were visualized in terms of b-factor (Figure 6). “Hot-spots,” i.e., residues that contribute significantly to complex stabilization (Clackson and Wells, 1995), were identified. As expected, the area that shows the highest stabilization was identified in the contact surface between anti-VEGF and VEGFA. Table 6 shows residues of anti-VEGFs that contribute to stabilization of complexes. The mutations (Ser105Thr; His101Tyr, Asn31His) carried out on Fab-bevacizumab to obtain ranibizumab (Yu et al., 2011), resulted in about a two-fold higher energy stabilization for the ranibizumab/VEGFA complex (Table 6). Interestingly, the stabilizing residues of VEGFR1d2_R2d3 bound to VEGFA are all basic amino acids, confirming the substantial contribute of the electrostatic contribution to the overall ΔE_binding, as mentioned above.

FIGURE 6

Figure 6. Three-dimensional projection of energy decomposition results. (A) VEGFR1d2_R2d3/VEGFA. (B) Fab-bevacizumab/VEGFA. (C) Ranibizumab/ VEGFA. Stabilizing effect of residues decreases from blue (stabilizing negative energy) to red (destabilizing positive energy). Green or yellow identify neutral (close to zero) contribution to the binding free energy. VEGFA bound to anti-VEGFA agents is highlighted by a black circle. Arrows highlight the areas of proteins with high RMSF.

TABLE 6

Table 6. Energy decomposition of predicted ΔE_binding.

Discussion

The main object of this study was the computational analysis, at molecular level, of binding between VEGFA and the three available anti-VEGF drugs, namely ranibizumab, bevacizumab (Fab-bevacizumab) and aflibercept (VEGFR1d2_R2d3). We have limited our study to interaction between VEGFA and binding domains of the above mentioned anti-VEGF drugs, excluding the Fc fragment of aflibercept and bevacizumab, because the Fc fragment does not seem to influence the pharmacodynamic properties of these drugs (Stewart, 2014b). Wu et al. (2013) already reported the modeling of cobercept (Li et al., 2014), which binds VEGFA with domains VEGFR1d2_R2d3_R2d4, however, no studies have analyzed in detail the energy components contributing to the complex VEGFR1d2_R2d3/VEGFA or compared complexes of different anti-VEGF agents. Therefore, the present study, for the first time, compares the interaction of the three different anti-VEGF agents with VEGFA; this is particularly relevant, considering that aflibercept is structurally unrelated to the other two agents. The entire computational study was carried out with open source tools and software packages.

Protein-protein docking, carried out with PyDock, was the first step. The rough energetic evaluation of complexes predicted by PyDock showed substantial difference between VEGFR1d2_R2d3, ranibizumab and Fab-bevacizumab, VEGFR1d2_R2d3/VEGFA being stabilized mainly by electrostatic energy, whereas the other two complexes were stabilized by VdW and desolvation energy. MD simulation (GROMACS) combined with MM-PBSA calculation (g_mmpbsa tool) was the second step. All MD simulations have been analyzed over a time sufficient to reach relative conformational minimum. Some systems (unbound VEGFR1d2_R2d3 and Fab-bevacizumab/VEGFA) showed high RMSD fluctuations, though secondary structure was conserved and cosine content of eigenvectors was low, indicating correct conformational sampling. MM-PBSA calculations confirmed most of results obtained with PyDock, such as the contribution of electrostatic energy to stability of VEGFR1d2_R2d3/VEGFA and the contribution of Van der Waals interaction energy to ranibizumab/VEGFA and Fab-bevacizumab/VEGFA. Furthermore, MM-PBSA provided energy contributions to ΔE_binding in good agreement with experimental binding data (Papadopoulos et al., 2012). MM-PBSA calculation carried out with g_mmpbsa seems at least as successful as other existing tools for analysis of protein-protein docking and MD (Spiliotopoulos et al., 2012; Corrada and Colombo, 2013), but has the advantage of being implemented for GROMACS output files. We obtained a good correlation between experimental K_D and apolar desolvation energy (ΔG_apolar), in accordance to the leading role of solvent exclusion, which is strictly related to the hydrophobic effect of inter- and intra-molecular assembly (Richmond, 1984; Chandler, 2005). Breaking down of ΔE_binding helped the identification of several features of complexes. The kinetics of aflibercept/VEGFA binding has been found to be characterized by fast K_on, which is consistent with substantial favorable electrostatic forces (Papadopoulos et al., 2012). Ranibizumab/VEGF-A had shown low experimental K_off (Papadopoulos et al., 2012), i.e., long lasting binding, and this could be related to higher number of contacts and H-bonds and less conformational freedom compared to the other two analyzed complexes (Copeland, 2011). The residue-based RMSF confirmed that upon binding ranibizumab stabilizes VEGFA in comparison to Fab-bevacizumab. VEGFR1d2_R2d3 also stabilizes VEGFA; however, the C-terminal and the N-terminal of VEGFR1d2_R2d3 are characterized by high RMSF, which may account for the higher K_off of aflibercept. Furthermore, the g_mmpbsa tool allowed us to identify the residues that favorably contribute to binding energy. A similar approach was carried out by Corrada and Colombo (2013) who have studied correlation of energetic parameters with affinity maturation of 17 variants of bevacizumab bound to VEGFA. A previous study reports two single nucleotide polymorphisms (SNPs) that may influence the binding at VEGFR2 (Wang et al., 2007):

• SNP1192G/A (rs2305948, in exon 7), that corresponds to a mutation Val279Ile in domain 3 of VEGFR2;

• SNP1719A/T (rs1870377, in exon 11), that corresponds to a mutation Q472H in domain 5 of VEGFR2.

SNP1192G/A is more interesting for our study because is located in the domain 3 of VEGFR2, included in aflibercept, and precisely in a beta sheet (one of the two anti-parallel beta sheets) where the residue interacts with other hydrophobic residues. Valine and isoleucine are both hydrophobic, however the isoleucine is more bulky than valine, and Val279Ile mutation might influence the stability of the beta sheets of domain 3. This mutation, along to other SNPs of VEGFRs, is worthy to be studied with the computational approach hereby described.

A recent binding study (Yang et al., 2014) reported an affinity of ranibizumab for VEGFA higher than that of aflibercept. At variance with computational studies, however, experimental binding data appear significantly influenced by the methodology used (Wang and Yang, 2013; Yadav et al., 2014). We have found correspondence between predicted ΔE_binding and contributing energy terms to the kinetic and affinity parameters of ranibizumab, bevacizumab and aflibercept, measured by Surface Plasmon Resonance (SPR) by Papadopoulos et al. (2012), who used captured anti-VEGF agents in the matrix and free ligands (VEGFA, VEGFB and PlGF) in the eluent. In contrast, the other SPR setup has used captured VEGFA and free anti-VEGF agents. In this latter setup, the finding that affinity of ranibizumab for VEGFA results greater than that of aflibercept is likely to depend on the blocked access of aflibercept to both ends of VEGFA (Yang et al., 2014).

Protein contact network approach provides a complementary analysis of evolution time for different topological properties of complexes. Fluctuations around a relative conformational minimum were characterized by properties that varied according to given constraints, whereas the energetic descriptor of proteins (dG_solv) did not change over the time, indicating that MD frames represented relative conformational minimum. The clustering of protein contact network revealed in the VEGFR1d2_R2d3/VEGFA that interaction wiring of VEGFA was conserved in comparison to the other two complexes, where the VEGFA network was altered. Furthermore, clustering of protein contact network of ranibizumab/VEGFA and Fab-bevacizumab/VEGFA was similar but not identical, because VEGFA showed a greater number of long-range interactions toward ranibizumab, in comparison to Fab-bevacizumab. Protein contact networks are built mainly considering VdW interactions; the observation that clustering of VEGFA in ranibizumab/VEGFA and Fab-bevacizumab/VEGFA was altered may be accounted for the higher contribution of VdW in such complexes, as observed with docking and MM-PBSA.

Some controversy exists about the correct length of MD (Dror et al., 2010; Genheden and Ryde, 2010, 2012). Long time scale simulation, in the micro- to millisecond range, is necessary whenever the phenomena observed are in evolution (e.g., binding and unbinding processes, protein folding, conformation transition) (Dror et al., 2010). In the present study we have simulated preformed macromolecular complexes in a water environment in order to carry out MM-PBSA, i.e., analysis of binding energy. To this end, the nanosecond scale seems adequate. Within this time scale, it has been shown that MM-PBSA carried on short (20 ns) independent replicas gives results comparable to a single longer simulation (Genheden and Ryde, 2010, 2012). Furthermore, longer simulations certainly increase sampling size of MD, but would not affect MM-PB(GB)SA results (Genheden and Ryde, 2012).

Our data may provide a theoretical basis for understanding differences in affinity for VEGFA of clinically used anti-VEGFs. Such differences may impact the pharmacodynamics and the therapeutic effectiveness of these biological drugs. However, the clinical outcomes of the VEGF Trap-Eye: Investigation of Efficacy and Safety in Wet AMD (VIEW) studies (Schmidt-Erfurth et al., 2014) indicate that ranibizumab and aflibercept are comparable with respect to the number of treatments and visual acuity gains when dosed in identical regimens and concentrations (0.5 mg). The clinical outcomes of Comparison of Age-related Macular Degeneration Treatments Trials (CATT) studies (Martin et al., 2012) indicate that ranibizumab and bevacizumab are equally effective on visual acuity over a 2-year period when used in same regimens but at different doses (0.5 and 1.25 mg for ranibizumab and bevacizumab, respectively). The same studies indicate that the proportion of patients with 1 or more systemic serious adverse events was higher with bevacizumab than ranibizumab (39.9 vs. 31.7%; adjusted risk ratio, 1.30; 95% CI, 1.07–1.57; P = 0.009). In general, affinity alone, even when it is assessed with high accuracy, does not always straightly correlate with clinical outcome. In fact, between affinity and clinical outcomes there are potency, dose, and regimen parameters to take into account in order to translate in silico, in vitro, and in vivo data from bench to bedside.

Here we presented in silico data of anti-VEGF/VEGFA complexes showing that they considerably differ both in terms of molecular interactions and stabilizing energy. Detailed understanding of such drug-target interactions may help in developing novel biological drugs.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported in part by the National grant PON01-00110. The authors thank Dr. Elisa Muscianisi for the fruitful scientific discussion. Authors acknowledge the CINECA Award N. HP10C8LAAA, 2013 for the availability of high performance computing resources and support. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphar.2015.00248

Abbreviations

VEGF(A,B,C…), Vascular endothelial growth factor and its isoforms; VEGFR (1, 2…), VEGF receptor and its isoforms; anti-VEGF, Drugs that block VEGF; AMD, Age related macular degeneration; DME, Diabetic macular edema; VEGFR1d2_D2d3, Aflibercept's binding domain; Fab-bevacizumab, Bevacizumab's binding domain; Fab, Fragment antigen-binding; Fc, Fragment crystallizable region; Ig, Human immune globulin; MD, Molecular dynamics; g_mmpbsa, Tool for molecular mechanics Poisson Boltzmann Surface Area calculations (MM-PBSA).

References

Amadei, A., Linssen, A. B., and Berendsen, H. J. (1993). Essential dynamics of proteins. Proteins 17, 412–425. doi: 10.1002/prot.340170408

PubMed Abstract | CrossRef Full Text

Baker, N. A., Sept, D., Joseph, S., Holst, M. J., and McCammon, J. A. (2001). Electrostatics of nanosystems: application to microtubules and the ribosome. Proc. Natl. Acad. Sci. U.S.A. 98, 10037–10041. doi: 10.1073/pnas.181342398

PubMed Abstract | CrossRef Full Text | Google Scholar

Balakrishnan, R. (2004). The energy of a graph. Linear Algebra Appl. 387, 287–295. doi: 10.1016/j.laa.2004.02.038

CrossRef Full Text | Google Scholar

Berendsen, H. J. C., Grigera, J. R., and Straatsma, T. P. (1987). The missing term in effective pair potentials. J. Phys. Chem. 91, 6269–6271. doi: 10.1021/j100308a038

CrossRef Full Text | Google Scholar

Carmeliet, P., and Jain, R. K. (2011). Molecular mechanisms and clinical applications of angiogenesis. Nature 473, 298–307. doi: 10.1038/nature10144

PubMed Abstract | CrossRef Full Text | Google Scholar

Chandler, D. (2005). Interfaces and the driving force of hydrophobic assembly. Nature 437, 640–647. doi: 10.1038/nature04162

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheatham, T. E., Miller, J. L., Fox, T., Darden, T. A., and Kollman, P. A. (1995). Molecular-dynamics simulations on solvated biomolecular systems - the particle mesh ewald method leads to stable trajectories of DNA, Rna, and proteins. J. Am. Chem. Soc. 117, 4193–4194. doi: 10.1021/ja00119a045

CrossRef Full Text | Google Scholar

Chen, Y., Wiesmann, C., Fuh, G., Li, B., Christinger, H. W., McKay, P., et al. (1999). Selection and analysis of an optimized anti-VEGF antibody: crystal structure of an affinity-matured Fab in complex with antigen. J. Mol. Biol. 293, 865–881. doi: 10.1006/jmbi.1999.3192

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, T. M., Blundell, T. L., and Fernandez-Recio, J. (2007). pyDock: electrostatics and desolvation for effective scoring of rigid-body protein-protein docking. Proteins 68, 503–515. doi: 10.1002/prot.21419

PubMed Abstract | CrossRef Full Text | Google Scholar

Chung, A. S., and Ferrara, N. (2011). Developmental and pathological angiogenesis. Annu. Rev. Cell Dev. Biol. 27, 563–584. doi: 10.1146/annurev-cellbio-092910-154002

PubMed Abstract | CrossRef Full Text | Google Scholar

Clackson, T., and Wells, J. A. (1995). A hot spot of binding energy in a hormone-receptor interface. Science 267, 383–386. doi: 10.1126/science.7529940

PubMed Abstract | CrossRef Full Text | Google Scholar

Copeland, R. A. (2011). Conformational adaptation in drug-target interactions and residence time. Future Med. Chem. 3, 1491–1501. doi: 10.4155/fmc.11.112

PubMed Abstract | CrossRef Full Text | Google Scholar

Corrada, D., and Colombo, G. (2013). Energetic and dynamic aspects of the affinity maturation process: characterizing improved variants from the bevacizumab antibody with molecular simulations. J. Chem. Inf. Model. 53, 2937–2950. doi: 10.1021/ci400416e

PubMed Abstract | CrossRef Full Text | Google Scholar

De Ruvo, M., Giuliani, A., Paci, P., Santoni, D., and Di Paola, L. (2012). Shedding light on protein-ligand binding by graph theory: the topological nature of allostery. Biophys. Chem. 165–166, 21–29. doi: 10.1016/j.bpc.2012.03.001

CrossRef Full Text | Google Scholar

Di Paola, L., De Ruvo, M., Paci, P., Santoni, D., and Giuliani, A. (2013). Protein contact networks: an emerging paradigm in chemistry. Chem. Rev. 113, 1598–1613. doi: 10.1021/cr3002356

PubMed Abstract | CrossRef Full Text | Google Scholar

Dror, R. O., Jensen, M. Ø., Borhani, D. W., and Shaw, D. E. (2010). Exploring atomic resolution physiology on a femtosecond to millisecond timescale using molecular dynamics simulations. J. Gen. Physiol. 135, 555–562. doi: 10.1085/jgp.200910373

PubMed Abstract | CrossRef Full Text | Google Scholar

Eisenberg, D., and McLachlan, A. D. (1986). Solvation energy in protein folding and binding. Nature 319, 199–203. doi: 10.1038/319199a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Feenstra, K. A., Hess, B., and Berendsen, H. J. C. (1999). Improving efficiency of large time-scale molecular dynamics simulations of hydrogen-rich systems. J. Comput. Chem. 20, 786–798.

Google Scholar

Gabb, H. A., Jackson, R. M., and Sternberg, M. J. (1997). Modelling protein docking using shape complementarity, electrostatics and biochemical information. J. Mol. Biol. 272, 106–120. doi: 10.1006/jmbi.1997.1203

PubMed Abstract | CrossRef Full Text | Google Scholar

Genheden, S., and Ryde, U. (2010). How to obtain statistically converged MM/GBSA results. J. Comput. Chem. 31, 837–846. doi: 10.1002/jcc.21366

PubMed Abstract | CrossRef Full Text | Google Scholar

Genheden, S., and Ryde, U. (2012). Will molecular dynamics simulations of proteins ever reach equilibrium? Phys. Chem. Chem. Phys. 14, 8662–8677. doi: 10.1039/c2cp23961b

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanout, M., Ferraz, D., Ansari, M., Maqsood, N., Kherani, S., Sepah, Y. J., et al. (2013). Therapies for neovascular age-related macular degeneration: current approaches and pharmacologic agents in development. Biomed Res. Int. 2013:830837. doi: 10.1155/2013/830837

PubMed Abstract | CrossRef Full Text | Google Scholar

Hess, B. (2002). Convergence of sampling in protein simulations. Phys. Rev. 65:031910. doi: 10.1103/physreve.65.031910

PubMed Abstract | CrossRef Full Text | Google Scholar

Hess, B. (2008). P-LINCS: a parallel linear constraint solver for molecular simulation. J. Chem. Theory Comput. 4, 116–122. doi: 10.1021/ct700200b

CrossRef Full Text | Google Scholar

Hess, B., Bekker, H., Berendsen, H. J. C., and Fraaije, J. G. E. M. (1997). LINCS: a linear constraint solver for molecular simulations. J. Comput. Chem. 18, 1463–1472.

Google Scholar

Holash, J., Davis, S., Papadopoulos, N., Croll, S. D., Ho, L., Russell, M., et al. (2002). VEGF-Trap: a VEGF blocker with potent antitumor effects. Proc. Natl. Acad. Sci. U.S.A. 99, 11393–11398. doi: 10.1073/pnas.172398299

PubMed Abstract | CrossRef Full Text | Google Scholar

Holz, F. G., Schmitz-Valckenberg, S., and Fleckenstein, M. (2014). Recent developments in the treatment of age-related macular degeneration. J. Clin. Invest. 124, 1430–1438. doi: 10.1172/JCI71029

PubMed Abstract | CrossRef Full Text | Google Scholar

Humphrey, W., Dalke, A., and Schulten, K. (1996). VMD: visual molecular dynamics. J. Mol. Graph. Model. 14, 33–38. doi: 10.1016/0263-7855(96)00018-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Iyer, S., Darley, P. I., and Acharya, K. R. (2010). Structural insights into the binding of vascular endothelial growth factor-B by VEGFR-1(D2): recognition and specificity. J. Biol. Chem. 285, 23779–23789. doi: 10.1074/jbc.M110.130658

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiménez-García, B., Pons, C., and Fernández-Recio, J. (2013). pyDockWEB: a web server for rigid-body protein-protein docking using electrostatics and desolvation scoring. Bioinformatics 29, 1698–1699. doi: 10.1093/bioinformatics/btt262

PubMed Abstract | CrossRef Full Text | Google Scholar

Jorgensen, W. L., Maxwell, D. S., and Tiradorives, J. (1996). Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 118, 11225–11236. doi: 10.1021/ja9621760

CrossRef Full Text | Google Scholar

Kar, P., Lipowsky, R., and Knecht, V. (2013). Importance of polar solvation and configurational entropy for design of antiretroviral drugs targeting HIV-1 protease. J. Phys. Chem. B 117, 5793–5805. doi: 10.1021/jp3085292

PubMed Abstract | CrossRef Full Text | Google Scholar

Kollman, P. A., Massova, I., Reyes, C., Kuhn, B., Huo, S., Chong, L., et al. (2000). Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 33, 889–897. doi: 10.1021/ar000033j

PubMed Abstract | CrossRef Full Text | Google Scholar

Kukic, P., Farrell, D., McIntosh, L. P., García-Moreno, E. B., Jensen, K. S., Toleikis, Z., et al. (2013). Protein dielectric constants determined from NMR chemical shift perturbations. J. Am. Chem. Soc. 135, 16968–16976. doi: 10.1021/ja406995j

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumari, R., Kumar, R., and Lynn, A. (2014). g_mmpbsa-A GROMACS tool for high-throughput MM-PBSA calculations. J. Chem. Inf. Model. 54, 1951–1962. doi: 10.1021/ci500020m

PubMed Abstract | CrossRef Full Text | Google Scholar

Leppänen, V. M., Prota, A. E., Jeltsch, M., Anisimov, A., Kalkkinen, N., Strandin, T., et al. (2010). Structural determinants of growth factor binding and specificity by VEGF receptor 2. Proc. Natl. Acad. Sci. U.S.A. 107, 2425–2430. doi: 10.1073/pnas.0914318107

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Xu, G., Wang, Y., Xu, X., Liu, X., Tang, S., et al. (2014). Safety and efficacy of conbercept in neovascular age-related macular degeneration: results from a 12-month randomized phase 2 study: AURORA study. Ophthalmology 121, 1740–1747. doi: 10.1016/j.ophtha.2014.03.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, D. F., Maguire, M. G., Fine, S. L., Ying, G. S., Jaffe, G. J., Grunwald, J. E., et al. (2012). Ranibizumab and bevacizumab for treatment of neovascular age-related macular degeneration: two-year results. Ophthalmology 119, 1388–1398. doi: 10.1016/j.ophtha.2012.03.053

PubMed Abstract | CrossRef Full Text | Google Scholar

Muller, Y. A., Chen, Y., Christinger, H. W., Li, B., Cunningham, B. C., Lowman, H. B., et al. (1998). VEGF and the Fab fragment of a humanized neutralizing antibody: crystal structure of the complex at 2.4 A resolution and mutational analysis of the interface. Structure 6, 1153–1167. doi: 10.1016/S0969-2126(98)00116-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Neugebauer, J., Reiher, M., Kind, C., and Hess, B. A. (2002). Quantum chemical calculation of vibrational spectra of large molecules - Raman and IR spectra for buckminsterfullerene. J. Comput. Chem. 23, 895–910. doi: 10.1002/jcc.10089

PubMed Abstract | CrossRef Full Text | Google Scholar

Nicholls, A., Sharp, K. A., and Honig, B. (1991). Protein folding and association - insights from the interfacial and thermodynamic properties of hydrocarbons. Proteins 11, 281–296. doi: 10.1002/prot.340110407

PubMed Abstract | CrossRef Full Text | Google Scholar

Ozboyaci, M., Gursoy, A., Erman, B., and Keskin, O. (2011). Molecular recognition of H3/H4 histone tails by the tudor domains of JMJD2A: a comparative molecular dynamics simulations study. PLoS ONE 6:e14765. doi: 10.1371/journal.pone.0014765

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, A. C., Borhani, D. W., Dror, R. O., and Shaw, D. E. (2013). Molecular determinants of drug-receptor binding kinetics. Drug Discov. Today 18, 667–673. doi: 10.1016/j.drudis.2013.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Papadopoulos, N., Martin, J., Ruan, Q., Rafique, A., Rosconi, M. P., Shi, E., et al. (2012). Binding and neutralization of vascular endothelial growth factor (VEGF) and related ligands by VEGF Trap, ranibizumab and bevacizumab. Angiogenesis 15, 171–185. doi: 10.1007/s10456-011-9249-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Pronk, S., Páll, S., Schulz, R., Larsson, P., Bjelkmar, P., Apostolov, R., et al. (2013). GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 29, 845–854. doi: 10.1093/bioinformatics/btt055

PubMed Abstract | CrossRef Full Text | Google Scholar

Richmond, T. J. (1984). Solvent accessible surface area and excluded volume in proteins. Analytical equations for overlapping spheres and implications for the hydrophobic effect. J. Mol. Biol. 178, 63–89. doi: 10.1016/0022-2836(84)90231-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmidt-Erfurth, U., Kaiser, P. K., Korobelnik, J. F., Brown, D. M., Chong, V., Nguyen, Q. D., et al. (2014). Intravitreal aflibercept injection for neovascular age-related macular degeneration: ninety-six-week results of the VIEW studies. Ophthalmology 121, 193–201. doi: 10.1016/j.ophtha.2013.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwede, T., Kopp, J., Guex, N., and Peitsch, M. C. (2003). SWISS-MODEL: an automated protein homology-modeling server. Nucleic Acids Res. 31, 3381–3385. doi: 10.1093/nar/gkg520

PubMed Abstract | CrossRef Full Text | Google Scholar

Spiliotopoulos, D., Spitaleri, A., and Musco, G. (2012). Exploring PHD fingers and H3K4me0 interactions with molecular dynamics simulations and binding free energy calculations: AIRE-PHD1, a comparative study. PLoS ONE 7:e46902. doi: 10.1371/journal.pone.0046902

PubMed Abstract | CrossRef Full Text | Google Scholar

Stewart, M. W. (2014a). Anti-VEGF therapy for diabetic macular edema. Curr. Diab. Rep. 14:510. doi: 10.1007/s11892-014-0510-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Stewart, M. W. (2014b). Pharmacokinetics, pharmacodynamics and pre-clinical characteristics of ophthalmic drugs that bind VEGF. Expert Rev. Clin. Pharmacol. 7, 167–180. doi: 10.1586/17512433.2014.884458

PubMed Abstract | CrossRef Full Text | Google Scholar

Tasdighian, S., Di Paola, L., De Ruvo, M., Paci, P., Santoni, D., Palumbo, P., et al. (2014). Modules identification in protein structures: the topological and geometrical solutions. J. Chem. Inf. Model. 54, 159–168. doi: 10.1021/ci400218v

PubMed Abstract | CrossRef Full Text | Google Scholar

Teicher, B. A. (2011). Antiangiogenic agents and targets: a perspective. Biochem. Pharmacol. 81, 6–12. doi: 10.1016/j.bcp.2010.09.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Vishveshwara, S., Ghosh, A., and Hansia, P. (2009). Intra and inter-molecular communications through protein structure network. Curr. Protein Pept. Sci. 10, 146–160. doi: 10.2174/138920309787847590

PubMed Abstract | CrossRef Full Text | Google Scholar

Wade, R. C., Gabdoulline, R. R., Lüdemann, S. K., and Lounnas, V. (1998). Electrostatic steering and ionic tethering in enzyme-ligand binding: insights from simulations. Proc. Natl. Acad. Sci. U.S.A. 95, 5942–5949. doi: 10.1073/pnas.95.11.5942

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., and Yang, J. (2013). Analysis of the binding affinity of vascular endothelial growth factor A (VEGF) to ranibizumab, aflibercept and bevacizumab. Invest. Ophthalmol. Vis. Sci. 54. Abstract retrieved from ARVO Annual Meeting Abstracts database. (Accession No. 1961).

Google Scholar

Wang, Y., Zheng, Y., Zhang, W., Yu, H., Lou, K., Zhang, Y., et al. (2007). Polymorphisms of KDR gene are associated with coronary heart disease. J. Am. Coll. Cardiol. 50, 760–767. doi: 10.1016/j.jacc.2007.04.074

PubMed Abstract | CrossRef Full Text | Google Scholar

Welti, J., Loges, S., Dimmeler, S., and Carmeliet, P. (2013). Recent molecular discoveries in angiogenesis and antiangiogenic therapies in cancer. J. Clin. Invest. 123, 3190–3200. doi: 10.1172/JCI70212

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Z., Zhou, P., Li, X., Wang, H., Luo, D., Qiao, H., et al. (2013). Structural characterization of a recombinant fusion protein by instrumental analysis and molecular modeling. PLoS ONE 8:e57642. doi: 10.1371/journal.pone.0057642

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, L., Lu, T., Tuomi, L., Jumbe, N., Lu, J., Eppler, S., et al. (2013). Pharmacokinetics of ranibizumab in patients with neovascular age-related macular degeneration: a population approach. Invest. Ophthalmol. Vis. Sci. 54, 1616–1624. doi: 10.1167/iovs.12-10260

PubMed Abstract | CrossRef Full Text | Google Scholar

Yadav, S., Demeule, B., Liu, J., and Shire, S. J. (2014). The relative affinities of ranibizumab and aflibercept for vascular endothelial growth factor A (VEGF) – an analytical ultracentrifugation study. Invest Ophthalmol. Vis. Sci. Abstract retrieved from ARVO Annual Meeting Abstracts database. (Accession No. 1942).

PubMed Abstract | Google Scholar

Yang, J., Wang, X., Fuh, G., Yu, L., Wakshull, E., Khosraviani, M., et al. (2014). Comparison of binding characteristics and in vitro activities of three inhibitors of vascular endothelial growth factor A. Mol. Pharm. 11, 3421–3430. doi: 10.1021/mp500160v

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, L., Liang, X. H., and Ferrara, N. (2011). Comparing protein VEGF inhibitors: in vitro biological studies. Biochem. Biophys. Res. Commun. 408, 276–281. doi: 10.1016/j.bbrc.2011.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, J., Ma, B., Tsai, C. J., and Nussinov, R. (2006). Structural stability and dynamics of an amyloid-forming peptide GNNQQNY from the yeast prion sup-35. Biophys. J. 91, 824–833. doi: 10.1529/biophysj.106.083246

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, H. X. (2001). Disparate ionic-strength dependencies of on and off rates in protein-protein association. Biopolymers 59, 427–433. doi: 10.1002/1097-0282(200111)59:6<427::AID-BIP1047>3.0.CO;2-7

PubMed Abstract | CrossRef Full Text | Google Scholar