Molecular Docking, Molecular Dynamics Simulations, and Free Energy Calculation Insights into the Binding Mechanism between VS-4718 and Focal Adhesion Kinase.

Focal adhesion kinase (FAK) is a 125 kDa nonreceptor tyrosine kinase that plays an important role in many carcinomas. Thus, the targeting of FAK by small molecules is considered to be promising for cancer therapy. Some FAK inhibitors have been reported as potential anticancer drugs and have entered into clinical development; for example, VS-4718 is currently undergoing clinical trials. However, the lack of crystal structural data for the binding of VS-4718 with FAK has hindered the optimization of this anticancer agent. In this work, the VS-4718/FAK interaction model was obtained by molecular docking and molecular dynamics simulations. The binding free energies of VS-4718/FAK were also calculated using the molecular mechanics generalized Born surface area method. It was found that the aminopyrimidine group formed hydrogen bonds with the C502 residue of the hinge loop, while the D564 residue of the T-loop interacted with the amide group. In addition, I428, A452, V484, M499, G505, and L553 residues formed hydrophobic interactions with VS-4718. The obtained results therefore provide an improved understanding of the interaction between human FAK and VS-4718. Based on the obtained binding mechanism, 47 novel compounds were designed to target the adenosine 5'-triphosphate-binding pocket of human FAK, and ensemble docking was performed to assess the effects of these modifications on the inhibitor binding affinity. This work is also expected to provide additional insights into potential future target design strategies based on VS-4718.

Introduction

Focal adhesion kinase (FAK, also known as protein tyrosine kinase 2, PTK2) is a nonreceptor tyrosine kinase and is one of the FAK family members. In a biological context, FAK plays a key role in the adhesion, motility, invasion, metastasis, and survival of cancer cells. FAK has been described as a protein that possesses increased tyrosine phosphorylation, which is of particular importance in many carcinomas, including papillary thyroid carcinoma,[1] neck cancer,[2−5] malignant melanomas,[6,7] bladder cancer,[8] intrahepatic cholangiocarcinoma,[9] ovarian cancer,[10] esophageal cancer,[11] breast cancer,[12] and pancreatic ductal adenocarcinoma.[13,14] Additionally, high levels of FAK in cancer patients are generally associated with poor prognosis. Thus, FAK has been considered a promising potential target for cancer therapy.

The FAK (1052 residues) consists of three domains (Figure ), namely, the N-terminal 4.1 protein, ezrin, radixin, moesin (FERM) homology domain (residues 35–355), the middle protein kinase domain (residues 422–680), and the C-terminal focal adhesion target (FAT) domain (residues 707–1052).[15] The kinase domain of FAK (also known as the catalytic domain) possesses a highly conserved amino acid sequence and structure that is formed from an N-terminal region (N-lobe), a C-terminal region (C-lobe), and a linker between the N- and C-lobes, which is referred to as a hinge loop.

Figure 1

Structures of FAK and VS-4718. (A) FAK (1052 residues) consists of three domains: (i) N-terminal 4.1 protein, ezrin, radixin, moesin (FERM) domain (residues 35–355), (ii) middle protein kinase domain (residues 422–680), and (iii) C-terminal FAT domain (residues 707–1052). (B) Kinase domain of FAK is formed from an N-terminal region (N-lobe) and a C-terminal region (C-lobe). (C) Structure of VS-4718 with a 50% inhibitory concentration (IC50) of 1.5 nM in the targeting of FAK.

As previously reported, conformational rearrangements of the kinase domain can take place upon inhibitor binding.[16−19] Generally, such rearrangements are based on one of the following processes: (i) movement of the P-loop (i.e., the glycine-rich loop), (ii) movement of the T-loop (namely, the activation loop or A-loop), or (iii) rotation of the α-C-helix in the N-lobe. It should be noted here that the DFG (Asp-Phe-Gly) motif is highly conserved and follows the T-loop, which serves as an important regulator of the various kinase activities.[20−23] However, partial rearrangement of the T-loop can be induced by different kinase inhibitors for open and/or closed conformations. For example, pyrazolobenzothiazine can bind with the open conformation of the T-loop of FAK (PDB ID: 4I4E(24)), while BI-4464 can bind with the closed conformation (PDB ID: 6I8Z(25)) (Figure S1). In addition, pyrrolo[2,3-d]thiazole can bind with the DFG-in (the side chain of D564 pointing into the active pocket of human FAK) motif of FAK (PDB ID: 3PXK(26)), while pyrazolobenzothiazine can bind with the DFG-out (the side chain of D564 pointing out the active pocket of human FAK) motif (PDB ID: 4I4F(24)) (Figure S2). Thus, the conformational rearrangement of the active pocket of the kinase domain plays an important role in the design of novel kinase inhibitors. However, the conformation that occurs in the binding of VS-4718 with FAK remains unclear. Therefore, for this study, four models of FAK were selected as the initial receptor structures for constructing VS-4718/FAK complexes, namely, FAK-I (DFG-in and T-loop open), FAK-II (DFG-in and T-loop closed), FAK-III (DFG-out and T-loop open), and FAK-IV (DFG-out and T-loop closed).

Till date, several small-molecule inhibitors of FAK have been reported, and some are currently under clinical development, including defactinib (VS-6063, PF-04554878), GSK2256098, VS-6062 (PF-00562271, PF-5662271), CEP-37440, BI-853520 (IN-10018), and VS-4718 (PND-1186) (Figure S3).[2,4,27−32] In addition, a number of effective FAK inhibitors have been proven to inhibit tumor growth and metastasis.[33] For example, VS-6063 is a highly effective FAK inhibitor that has completed phase II clinical trials in patients suffering from KRAS mutant nonsmall cell lung cancer.[34] Furthermore, GSK2256098, as a reversible adenosine 5′-triphosphate (ATP)-competitive inhibitor, is currently undergoing clinical trials for patients with advanced solid tumors.[35,36] Furthermore, VS-6062 has completed phase I clinical trials in the treatment of advanced solid tumors.[37]

As the compound of interest in the current study, VS-4718 is a reversible and selective inhibitor that exhibits an IC50 (50% inhibitory concentration) value of 1.5 nM in an in vitro kinase assay.[38] The VS-4718 has also demonstrated a cellular IC50 of ∼100 nM in malignant pleural mesothelioma cell lines and breast carcinoma. Furthermore, it has a median relative cellular IC50 of 1.22 μM against the pediatric preclinical testing program cell line.[38−42] These results indicate that VS-4718 possesses on-target and off-target concentrations of <100 nM and >1.0 μM, respectively, in these cell lines.[39] In terms of its clinical development, VS-4718 has been demonstrated to act as a potential inhibitor for triple-negative breast cancer stem cells.[43] It is currently being evaluated in advanced cancer (NCT02651727), metastatic nonhematologic malignancies (NCT01849744), and acute myeloid or B-cell acute lymphoblastic leukaemia (NCT02215629) (Table S1). Additionally, five off-targets for VS-4718 have been reported with >65% inhibition at 1 μM,[38] thereby indicating that optimization of VS-4718 is necessary to obtain selective FAK inhibitors based on VS-4718. In this context, the binding mechanism of VS-4718 could be used to provide reference for the development of new antitumor drugs to target FAK.

It is known that receptor–ligand interactions play an important role in the elucidation of drug-target binding mechanisms; hence, a reliable receptor–ligand structure is essential to permit structure- or knowledge-based drug development. Because molecular docking can create a static image of the drug–target complex, it has been employed in the area of drug design.[44−46] Furthermore, molecular dynamics (MD) simulations have also been employed to help elucidate the interactions present during drug binding.[47,48]

Till date, no crystal structural data are currently available for the binding of VS-4718 with FAK, which hinders the development of superior FAK inhibitors. We herein report a VS-4718/FAK binding model obtained through the use of molecular docking and all-atom MD simulations. More specifically, the binding free energies are also calculated using the molecular mechanics generalized Born surface area (MM/GBSA) method, which is a valuable and powerful method for carrying out binding free energy calculations.[49−56] In addition, the hot residues present in the VS-4718/FAK binding models, which are likely to alter the binding affinity of VS-4718 with FAK, are identified by analyzing the energy decomposition for each residue. Ultimately, our aim is to carry out simulations that will provide the binding mechanism for VS-4718 with FAK, in addition to useful information that will facilitate the development of innovative FAK inhibitors.

Materials and Methods

Molecular Docking

To date, more than 35 crystal structures have been published for human FAK (UniProt ID: Q05397) in the PDB,[15,24−26,57−72] with over 24 being known to form complex structures between the FAK kinase domain and either an inhibitor or ADP (Table S2). When bound to an inhibitor, the T-loop (i.e., the active, or A-loop) of FAK can form an open or closed conformation; for example, 4-((4-(((1R,2R)-2-(dimethylamino)cyclopentyl)amino)-5-(trifluoromethyl)pyrimidin-2-yl)amino)-N-methylbenzenesulfonamide binds with the open conformation of the T-loop (PDB ID: 6YVY(58)), while N-methyl-N-(3-(((2-((2-oxo-1,2,3,4-tetrahydroquinolin-6-yl)amino)-5-(trifluoromethyl)pyrimidin-4-yl)amino)methyl)pyridin-2-yl)methanesulfonamide binds with the closed conformation (PDB ID: 6YQ1(58)) (Figure S4). In addition, the DFG domain, which is also an important loop for determining the kinase activity, can form the in or out conformation, such as in the case of DFG-in upon binding with N-(3-(((5-cyano-2-phenyl-1H-pyrrolo[2,3-b]pyridin-4-yl)amino)methyl)pyridin-2-yl)-N-methylmethanesulfonamide (PDB ID: 4GU6(65)) or DFG-out upon binding with 1-(4-(6-amino-9H-purin-9-yl)phenyl)-3-(3-(tert-butyl)-1-(p-tolyl)-1H-pyrazol-5-yl)urea (PDB ID: 4K9Y(66)) (Figure S5). However, little information is currently available regarding the conformation of the FAK active pocket with VS-4718. Therefore, we employed four representative conformations of the active groove to construct structures for the VS-4718/FAK complex in our docking study. More specifically, these conformational models are as follows: (i) with an open T-loop and DFG-in, FAK-I (PDB ID: 6YVY(58)), (ii) with a closed T-loop and DFG-in, FAK-II (PDB ID: 6I8Z(25)), (iii) with an open T-loop and DFG-out, FAK-III (PDB ID: 4K9Y(66)), and (iv) with a closed T-loop and DFG-out, FAK-IV (PDB ID: 4EBV(64)) (Figure S6).

The three-dimensional crystal structure of human FAK (UniProt[73] ID: Q05397) cocrystallized with BI-4464 (PDB ID: 6I8Z(25)) was obtained from the PDB,[74,75] as were those of 6YVY,[58]4K9Y,[66] and 4EBV.[64] These four crystal structures were selected as models for molecular docking. PYMOL 2.1[76] was used to prepare the cocrystallized structures of FAK by removing the cocrystallized ligand molecules. The crystallographic water molecules were removed to obtain the final protein structures and maintain chain A within the molecular docking.

Two-dimensional structures of VS-4718 and the other compounds were sketched using InDraw software and were converted into three-dimensional structures using OpenBabel 3.1.[77] These structures were then minimized with the semiempirical PM3 method[78] using MOPAC2016 (Stewart Computational Chemistry, Colorado Springs, CO; http://OpenMOPAC.net). The structures of the protein, VS-4718, and all other compounds were pretreated using AutoDockTools 1.5.6[79] with hydrogen atoms included, followed by Gasteiger charging[80] and unreasonable atomic overlap adjustment. A 40 × 40 × 40 grid map with a 0.375 Å grid spacing was generated using AutoGrid v.4.2;[81] this grid map was based on the center of the ATP-binding groove for FAK. Two hundred conformations per system were generated using the Lamarckian genetic algorithm[82] in AutoDock v.4.2.[81] Finally, the optimal conformation for each docking model was selected based on the docking experiment that gave the best rational orientation in the active pocket, as referenced to previously reported crystal complex structures for inhibitor/FAK.[25,58]

MD Simulations

The MD simulations have become increasingly important in the context of understanding the interactions between receptors (e.g., proteins, enzymes, or cyclodextrins) and ligands (e.g., inhibitors, stabilizers, and supermolecules).[54,55,83−86] Therefore, MD simulations were carried out to explore the binding between VS-4718 and the FAK protein. The initial VS-4718/FAK complex structures obtained from the docking model were used as the initial complex conformations for subsequent simulations. Because the standard force field of a small molecule cannot be obtained, the general Amber force field generation procedure (version 2; GAFF2)[87] was used to generate the force fields of the various ligand molecules examined in this work, including that of VS-4718. For this purpose, the ligand structure was drawn using InDraw and translated into three dimensions using OpenBabel v3.1 prior to structural optimization at the B3LYP/6-31G* level of theory using Gaussian09 software.[88] Then, electrostatic potential was also calculated using the B3LYP/6-31G* method. Subsequently, the restrained electrostatic potential protocol[89] was used to fit the partial atomic charges of the small molecule. Meanwhile, the standard Amber ff19SB force field[90] was employed to create the topology parameters for the FAK protein. The standard residue protonation approach was employed based on the residue obtained in the ff19SB force field (pH = 7.0). Subsequently, the ligand/FAK complex systems were solvated with a cuboid box of TIP3P water[91] at a relative distance of 15 Å from all protein atoms, and a single chloride ion was added to ensure that the entire system was in an electrically neutral state. The final system included the FAK protein, the ligand molecule, and the solvent water molecules.

To avoid the unfavorable interactions produced by additional solvents and ions, the system was minimized by initially restraining the atoms of the protein and ligand molecules to optimize the coordinates of the water molecules and counterions, wherein the weight for the positional restraints was set at 2.0 kcal/mol/Å2. Subsequently, minimization of the overall system was carried out without any constraints. Langevin dynamics with a 2.0 ps–1 collision frequency was used to increase the overall system temperature from 0 to 300 K in an NVT ensemble, while the isotropic position scaling method was applied to maintain the system pressure at 1 bar. Subsequently, the NPT ensemble was applied to equalize the system at 300 K and 1 bar. Simulations were then performed for heating, the application of a constant pressure, and the process of equilibrium maintenance for 200 ps in each case. Finally, the two systems with different conformations of FAK were subjected to a 500 ns MD simulation under the above conditions to collect the necessary data required for analysis.

In these simulations, periodic boundary conditions were used to avoid unphysical edge effects, while the SHAKE algorithm[92] was employed to constrain the covalent bonds of the protein and ligand hydrogen atoms. To mitigate long-range electrostatic forces, the particle mesh Ewald algorithm[93] was applied during the simulation, wherein a cutoff distance of 12 Å was employed to simplify the short-range electrostatic forces. This cutoff was also used to simplify the van der Waals interactions. The seed for the pseudo-random number generator was based on the current date and time. System preparation and analysis of the simulation results were conducted using AmberTools21.[94] However, it should be noted that the simulations were performed using the CUDA version of PMEMD in AMBER20[94] to decrease the required simulation time through the use of Nvidia GPUs. The coordinates of the simulation were saved every 10 ps to obtain the trajectories for the complex systems. The CPPTRAJ[95,96] model was employed to analyze the data obtained from the MD trajectories.

Binding Free Energy Calculations

Qualitative and quantitative analyses are both important when determining the binding free energies of inhibitor–protein binding interactions. Currently, several methods are available for estimating the absolute binding free energy between a protein and an inhibitor, such as the linear interaction energy,[97−99] the linear response approximation,[100] the solvated interaction energy,[101] the free energy pathway,[102] and the molecular mechanics Poisson Boltzmann (or generalized Born) surface area (MM-PB/GBSA)[103,104] approach. However, the MM/GBSA approach is considered the most efficient means for evaluating ligand and enzyme systems,[50,105] and its framework can be represented by the following equationswhere Gcomplex, Gprotein, and Gligand are the free energies of binding for VS-4718 + FAK, FAK, and VS-4718, respectively. The value of ΔGbinding can be summarized based on the enthalpic term (ΔH = ΔEgas+ ΔEsol) and the entropic term (TΔS), while the molecular mechanical energy (Egas) can be divided into van der Waals forces (EvdW), electrostatic forces (Eele), and the intramolecular energy (Eint), which can themselves be obtained from a statistical average based on the ff19SB force field. In addition, the solvation free energy (Esolv) can be decomposed into electrostatic (EGB) and nonelectrostatic (Esurf) components. More specifically, the Esurf component is the combined effect of the unfavorable cost of surface formation and the favorable van der Waals interactions between the solute and the solvent. In this term, Esurf can be evaluated using γ·SA + b, where γ = 0.0072 kcal/Å2 and b = 0.0 kcal/mol and where the LCPO method[106] is used to estimate the solvent accessible surface area (SA). Furthermore, the GB equation[107,108] was applied to calculate the EGB contribution. For the purpose of this work, the dielectric constants for the solute and the exterior were set to 1 and 80, respectively. One thousand snapshots were extracted from the final 200 ns of the MD trajectory to give a statistical average for the MM/GBSA method, and an entropic term was added to improve the accuracy for the binding of VS-4718 with FAK. Normal model analysis coupled with the quasi-harmonic model was applied to estimate the entropic contribution, which was based on 100 snapshots from the final 200 ns of the MD trajectory.[109]

To determine the essential residues involved in the binding interactions between VS-4718 and FAK, the contribution of each residue was evaluated at the atomic level using the energy decomposition method to obtain the binding free energy. The electrostatic contribution to the solvation energy was determined from the GB model, as in the case of the binding free energy calculation, while the single trajectory simulation was used to calculate the binding free energy; the internal energy calculation was excluded. Moreover, based on the corresponding SA, the nonpolar solvation energy per atom was obtained, and the entropic contribution of each residue was included for the energy decomposition calculations. The contribution of each given residue was estimated at the atomistic level by summing the contribution of the overall atoms present in the residue. Similarly, the contributions of the backbone and the side chain were determined by summing the relevant atoms in each case. One thousand snapshots extracted from the 200 ns simulation were also used to estimate the energy decomposition. Finally, all energies involved in the binding of VS-4718 with FAK were calculated using the AMBER MMPBA.py program.[110] All methods employed during this study are outlined in Figure S40.

Results and Discussion

Initial Models

Prior to carrying out the docking experiment, a redocking strategy was employed to evaluate the docking power of the docking procedure. The docking power is defined as the root-mean-square deviation (rmsd) between the conformation from molecular docking and the crystal conformation of the ligand molecule. More specifically, the selected crystal complex structures for 6YVY (PDB ID, FAK-I), 6I8Z (FAK-II), 4K9Y (FAK-III), and 4EBV (FAK-VI) were redocked, and the rmsd between the crystal structure conformation of the ligand and the conformation with the lowest energy (i.e., −8.73, −11.90, −13.25, and −10.03 kcal/mol for FAK-I, FAK-II, FAK-III, and FAK-IV, respectively) for docking was found to be <1.0 Å in each case (Figure S7). It should be noted here that the ATP binding site of FAK was also considered in the FAK-VI model because the original ligand was bound to the allosteric site (Figure S8). Thus, as an example of one redocking experiment, the lowest binding free energy conformation of redocking for the FAK-II model was determined to be −13.25 kcal/mol, and the conformation number of this cluster was 1960, which gives an occupancy of 98% based on the 2000 docking conformations (Figure S9). These results indicate that the docking strategy is suitable for docking the ligand molecule at the ATP binding site. In addition, it was found that VS-4718 was bound to the ATP binding site of FAK, and therefore the same docking strategy was used to construct the VS-4718/FAK-I, VS-4718/FAK-II, VS-4718/FAK-III, and VS-4718/FAK-IV complexes.

The initial conformations for VS-4718 binding with FAK were obtained from docking experiments, wherein the lowest binding free energies were −8.74, −10.11, and −8.61 kcal/mol for the FAK-I, FAK-II, and FAK-III systems, respectively (Figure S10). Although the binding models based on FAK-I and FAK-II (i.e., with DFG-in) were similar (Figure S11), the binding models with DFG-out (i.e., FAK-III and FAK-IV) were not considered in the context of the defined standards for competitive ATP inhibitors.[65,66,111,112] More specifically, there was a large spatial similarity between the ligands and a competitive ATP inhibitor, and hydrogen bonds were not present between the ligands and the hinge loop of FAK. These results suggest that the DFG-in conformation is the most plausible in the binding of VS-4718 with FAK. Furthermore, the greater binding affinity between VS-4718 and FAK-II (cf., FAK-I, based on the docking score) indicates that the closed conformation of the T-loop plays an important role in enhancing the binding of VS-4718 with FAK. Thus, the model of VS-4718 binding with FAK-II was selected for further experiments.

To obtain a more accurate conformation of VS-48718 binding with FAK-II, the subsequent docking experiment incorporated 2000 different conformations that could be divided into >5 clusters (Figure S12), wherein the lowest free energy conformation within each cluster was selected as its representative conformation. The binding models of cluster 1 (defined as FAK-II-1) and 4 (defined as FAK-II-4) were considered to meet the above-described standards for competitive ATP inhibitors. Thus, two binding conformations (i.e., the conformation with the lowest binding score in cluster 1 or cluster 4) were selected for further examining the binding of VS-4718 with human FAK (Figure ), wherein the two nitrogen atoms (N3 and N5) of the pyrimidine–amine moiety form two hydrogen bonds with the nitrogen and oxygen atoms in the C502 residue (Figure S13). In addition, an unfavorable acceptor–acceptor interaction was found between the C502 residue and the methoxy oxygen atom of VS-4718. Furthermore, in the case of FAK-II-1, the O1 atom of the ligand forms a hydrogen bond with the D564 residue in the DFG domain of FAK. Together, these data suggest that the pyrimidine-amine moiety of the inhibitor plays an important role in the binding characteristics. Moreover, the pyrimidine side chains, such as the 2-methoxy-4-morpholinophenyl group, are responsible for orienting the conformation of the inhibitor when bound with FAK. Similarly, the trifluoromethyl group occupies the hydrophobic cavity, while the N-methylbenzamide group points toward the T-loop region.

Figure 2

Docking model for the binding of VS-4718 with human FAK. (A) Binding model for FAK-II-1, (B) binding model for FAK-II-4, and (C) overlay of the binding models for FAK-II-1 and FAK-II-4. (D) Rotated from (C). The human FAK protein is represented in its cartoon form wherein the helical, sheet, and loop structures are colored cyan, magenta, and orange, respectively. The VS-4718 ligand is shown in the stick format in cyan for FAK-II-1 and pink for FAK-II-4.

System Stability

Because the docking experiments did not consider any additional interactions between VS-4718 and FAK or rearrangement of the residues present in the active site of FAK, MD simulations were carried out for the two VS-4718/FAK complex systems to obtain additional information relating to the binding mechanism. More specifically, MD simulations were performed over 500 ns for the various protein–ligand complexes, at which point the rmsd values of the heavy atoms in the protein backbone and in the ligands had reached a plateau (Figure S14). Only small fluctuations of the rmsd values were observed for the protein kinase domain, with values of 1.77 ± 0.25 and 1.77 ± 0.19 Å being determined for FAK-II-1 and FAK-II-4, respectively, thereby indicating that the FAK structures were stable. In contrast, the fluctuations for the ligands were 1.51 ± 0.26 and 2.00 ± 0.24 Å for FAK-II-1 and FAK-II-4, respectively, which indicates that VS-4718 was reorientated in the active pocket for the FAK-II-4 system. This result agrees with the radius of gyration (Figure S15) and the surface area (Figure S16) determined for the VS-4718/FAK complex system, and therefore our observations indicate that these MD simulations are suitable for analyzing the interactions between VS-4718 and FAK.

The root-mean-square fluctuation (RMSF), which is based on the fluctuation of residues, was then used to quantify the stabilities of specific residues during the MD simulations. For this purpose, the RMSF values were analyzed based on the 500 ns MD trajectory for each VS-4718/FAK complex system (Figure S17), and it was found that for both FAK-II-1 and FAK-II-4, the amino acid residues at FAK positions 564–592 (i.e., the T-loop) exhibited a greater degree of fluctuation than those present in other regions (Figure A). This observation, which is consistent with the rmsd results, indicates that the T-loop of FAK is unstable during these simulations. In addition, the A-loop conformations (564–592 for human FAK) were in good agreement with those of other kinases, such as salt-inducible kinase 2,[83,105] microtubule affinity-regulating kinase 2,[113] and microtubule affinity-regulating kinase 4.[114] In our simulations, the initial conformation of FAK-II was closed, which promoted a greater degree of VS-4718 binding. However, as previously mentioned, the T-loops can form open or closed conformations, which are able to recognize various inhibitors; other protein kinases also have a similar recognition mechanism.[115,116]

Figure 3

Fluctuation in the FAK conformation during binding with VS-4718. (A) RMSF for the FAK backbone residues during the 500 ns simulation. Frames of the VS-4718/FAK complexes at 0 nm (green) and 500 ns (cyan) for the (B) FAK-II-1 and (C) FAK-II-4 systems. The FAK protein is represented in its cartoon form, while the VS-4718 ligand is shown in the stick format. The T-loop of FAK is highlighted.

Based on the obtained docking scores, the closed conformation of the FAK T-loop was identified as the optimal conformation for VS-4718 binding with FAK, as indicated by the snapshots obtained at 100, 200, 300, 400, and 500 ns from the MD trajectory. These snapshots were aligned and are shown for the FAK-II-1 and FAK-II-4 system in Figures S18 and S19, respectively. The initial and final structures of the two complexes were also extracted to check the possible structural changes during the MD simulation (Figure B,C), and it was found that the structure of the FAK A-loop (i.e., the T-loop) adopted an induced-fit conformation in the presence of an inhibitor. However, in the DFG domain, the T-loop maintained a closed conformation, thereby indicating that overall, VS-4718 is bound to the ATP binding site of FAK via a closed conformation of the T-loop.

Analysis of the Hydrogen-Bonding Network

The initial docking models identified a hydrogen bonding network between VS-4718 and FAK. However, only a few hydrogen bonds were formed between VS-4718 and the FAK protein (Figure ). Some potential residues for hydrogen bonds were identified in the active site, but less than two hydrogen bonds were found in the majority of simulation times, with an average number of 2.43 and 1.88 hydrogen bonds being determined for FAK-II-1 and FAK-II-4, respectively (Figure S20). The occupancy of the hydrogen bonds was also assessed in the 500 ns simulation for the VS-4718/FAK complexes, wherein the N3 atom of VS-4718 acted as a donor atom to form a hydrogen bond with the C502 residue of FAK (82.98 and 61.75% occupancies for FAK-II-1 and FAK-II-4, respectively). The oxygen atom of the C502 residue also took part in hydrogen bonding with the N5 atom of VS-4718, acting as an acceptor, with occupancies of 96.53 and 66.64% being determined for the FAK-II-1 and FAK-II-4 systems, respectively. Furthermore, it was determined that two of the hydrogen bonds formed between FAK and VS-4718 were located within the diamino-pyridine ring; these interactions represent the universal hydrogen bonds that are commonly found for FAK inhibitors[30,31,117] and other protein kinase inhibitors.[118,119] Moreover, the 7H-pyrrolo[2,3-d]pyrimidine,[120,121] diamino-pyrimidine,[58,122−124] and thieno[3,2-d]pyrimidine[125] rings also contribute to the binding model in a manner similar to the indazole ring, thereby indicating that the indazole ring can be replaced with a pyrazole or pyridine ring.

Figure 4

Hydrogen bond analysis for the VS-4718/FAK systems. (A) Distribution of the number of hydrogen bonds for total 50 000 frames during the 500 ns simulation. (B) Occupancy of each hydrogen bond as a percentage of the investigated period (500 ns), during which specific hydrogen bonds were formed. A hydrogen bond was defined when the distance between the acceptor and donor atoms was <3.5 Å and the internal acceptor···H-donor angle was >120°. Also shown are schematic diagrams of hydrogen bonding in the (C) FAK-II-1 and (D) FAK-II-4 systems.

In addition, the D564 residue present in the DFG region formed hydrogen bonds with VS-4718, although these were weaker than those formed at the C502 residue. Interestingly, in the FAK-II-1 system, the hydrogen bond between this residue and the O1 atom of VS-4718 (46.11% occupancy) was not observed in the FAK-II-4 system. However, greater hydrogen bonding occupancy was found between the D564 residue and the trifluoromethyl group of FAK-II-4 compared to that of the FAK-II-1 system, which was attributed to translation of the N-methylbenzamide group conformation to orientate the conformation of VS-4718 into the ATP binding pocket. Additionally, the D564 residue plays an important role in the coordination of a Mg2+ ion to stabilize the Mg–ATP interactions at the catalytic subunit of the cyclic adenosine monophosphate-dependent protein kinase.[126] As a result, repositioning of the D564 side chain would prevent ATP binding at the binding site.

Subsequently, MD simulation trajectories were employed to calculate the distances between the donor and acceptor atoms, in addition to the angles between the donor, hydrogen, and acceptor atoms. More specifically, the hydrogen bond between the C502 residue of FAK and the N3 atom of VS-4718 in the FAK-II-1 system (distance = 3.15 ± 0.24 Å) was stronger than the FAK-II-4 system (distance = 3.42 ± 0.40 Å) (Figure S21), and a similar trend was also observed for the hydrogen bond between the C502 residue and the N5 atom (distances = 3.15 ± 0.17 and 3.44 ± 0.33 Å for the FAK-II-1 and FAK-II-4 systems, respectively) (Figure S22). In addition, similar hydrogen bonds were found between the N3 atom and the C502 residue, indicating that these two hydrogen bonds play the same key role in the binding of VS-4718 with FAK, as confirmed experimentally for the interactions between FAK and BI-4464 (distances = 3.05 and 3.06 Å, PDB ID: 6I8Z(25)) (Figure S23). Furthermore, hydrogen bonding between the D564 residue and both the trifluoromethyl group and the O1 atom of VS-4718 were examined (Figures S24 and S25). It was found that replacing the O1 atom in the FAK-II-1 system with a trifluoromethyl group in the FAK-II-4 system reduced the binding affinity, as confirmed by measurement of the corresponding angles (Figures S26 and S27). Overall, these results suggest that the FAK-II-1 system is a superior model for investigating the binding of VS-4718 with FAK, wherein hydrogen bonding plays an important role in orienting VS-4718 into the ATP binding site. Our results also suggest that the methyl group of the N-methylbenzamide moiety can also be substituted for ethyl, cyclopropyl, and propyl chains to retain a similar binding affinity.

Analyses of the Interaction Fingerprints

To quantitatively characterize the interactions between VS-4718 and the human FAK ATP binding pocket, the interaction fingerprints of the FAK-II-1 and FAK-II-4 complexes were calculated using IChem[127−130] during the final 200 ns of the MD simulations (Figure ). It was found that the I428, V436, A452, V484, M499, L501, E506, S509, N551, L553, L567, and S568 residues of FAK form hydrophobic interactions with the inhibitor in both systems, wherein the A452, V484, M499, and L567 residues provide the greatest stabilization with a proportion of ∼1.00 [defined as the proportion of frames exhibiting interactions divided by the total (10 000) frames obtained in the final 200 ns of the simulation]. In addition, the L553 residue can form hydrophobic interactions with the pyridine moiety of VS-4718, while the I428 residue forms hydrophobic interactions with the methoxybenzene ring. Additionally, the L567 residue forms hydrophobic interactions with the benzene ring of the N-methylbenzamide group. It was also found that the S568 and N551 residues contributed toward hydrophobic interactions in the FAK-II-4 system but not in the FAK-II-1 system, and therefore the hydrogen bond attributed to the D564 residue in the FAK-II-1 system (see Figure B) can be considered to compensate for these missing interactions. As mentioned above, fingerprint analysis confirmed that the D564 residue forms a hydrogen bond with VS-4718 in the FAK-II-1 system (proportion = 0.64) but not in the FAK-II-4 system (proportion = 0.02). However, the hydrophobic segment of D564 contributes similarly to binding in both systems (proportions = 0.97 and 1.00 for the FAK-II-1 and FAK-II-1 systems, respectively). Moreover, C502 was found to form hydrogen bonds in both systems through interaction with the diaminopyridine group of the inhibitor (proportion ∼1.00), and this result is consistent with the hydrogen bond analyses section in this work. The above results therefore indicate that hydrogen bonding between the inhibitor and the C502 and D564 residues is key to stabilizing the orientation of the VS-4718 molecule within the binding pocket.

Figure 5

(A) Interaction fingerprints between FAK and VS-4718 in the final 200 ns of the MD simulations of the two systems over 10 000 frames. (B,C) Key residues involved in the hydrogen bonding (green) and hydrophobic (gray) interactions with VS-4718 (yellow) for the FAK-II-1 (B) and FAK-II-4 (C) systems.

Binding Free Energies

As described above, the interactions between VS-4718 and the FAK can be depicted by docking experiments and MD simulations. However, these procedures did not allow the determination of binding affinities of the inhibitor to the two model systems. Thus, the MM/GBSA method was employed to calculate the absolute VS-4718/FAK binding free energies for both systems (Tables , S3, and S4). It was found that all entropy values (TΔStotal) and enthalpies (Egas + Gsol) were negative (i.e., less than −23.10 and −55.80 kcal/mol, respectively), which indicates that the formation of a binding complex is an enthalpy-driven process. For FAK-II-1, the calculated binding free energy (ΔGbindcal) was approximately −33.13 kcal/mol, whereas for FAK-II-4, the ΔGbindcal value was −32.70 kcal/mol, thereby confirming that both the FAK-II-1 and FAK-II-4 systems constituted the preferred models. This conformation of FAK-II-1 was also previously observed in the crystal structure of a methanesulfonamide diaminopyrimidine inhibitor bound with FAK (PDB ID: 3BZ3(67)) (Figure S28). Overall, the perfect binding model for VS-4718/FAK involves the side chain (N-methylbenzamide group) of the ligand pointing toward the DFG region, which permits stable binding at the ATP binding site of FAK. The obtained binding free energy therefore indicates that VS-4718 can bind strongly with FAK, as confirmed by the experimentally obtained IC50 value of 1.5 nM.[38]

Table 1

Binding Free Energy for the VS-4718/FAK Complex and Decomposition into Electrostatic Interactions, van der Waals Interactions, Solvation Free Energies, and Entropy Valuesa

energy (kcal/mol)	FAK-II-1	FAK-II-4
ΔEvdW	–61.39 (2.79)*	–60.89 (2.58)
ΔEele	–14.44 (3.21)	–11.86 (3.42)
ΔEGB	25.34 (2.57)	24.28 (2.83)
ΔEsurf	–7.52 (0.27)	–7.33 (0.27)
ΔEgas	–75.83 (4.05)	–72.75 (4.44)
ΔEsolv	17.81 (2.57)	16.95 (2.81)
ΔEgas + ΔEsolv	–58.02 (3.19)	–55.80 (3.04)
ΔTStotal	–24.89 (5.41)	–23.10 (5.48)
ΔGbindcal	–33.13 (6.28)	–32.70 (6.27)

ΔEvdW: contribution of the van der Waals energy to the free energy of binding; ΔEele: contribution of the electrostatic energy to the free energy of binding; ΔEGB: contribution of the polar solvation energies to the free energy of binding; ΔEsurf: contribution of from the nonpolar solvation energies to the free energy of binding; ΔEgas: contribution of ΔEvdW + ΔEele to the free energy of binding; ΔEsolv: contribution of ΔEGB + ΔEsurf to the free energy of binding; ΔTStotal: contribution of the entropy energy to the free energy of binding; ΔGbindcal: the final estimated binding free energy from ΔEgas + ΔEsolv – ΔTStotal. *The uncertainties (shown in parentheses) were calculated as the root-mean-square error for each frame extracted during the MM/GBSA process.

The binding energy can usually be decomposed into its polar (Eele + EGB) and nonpolar (EvdW + Esurf) terms. More specifically, for the FAK-II-1 and FAK-II-4 systems examined herein, the polar terms were determined to be 10.90 and 12.42 kcal/mol, respectively, using the MM/GBSA method. It should be noted here that a positive value for the polar contribution indicates that the polar interactions between VS-4718 and FAK are antagonistic to this binding. In contrast, van der Waals (EvdW) interactions (−61.39 and −60.89 kcal/mol for FAK-II-1 and FAK-II-4, respectively) acted as the main nonpolar contribution and were conducive to binding. Furthermore, the nonpolar terms were determined to be −68.91 and −68.22 kcal/mol for FAK-II-1 and FAK-II-4, respectively. These similar values for the two systems support our previous observation that they exhibited a similar binding affinity toward the inhibitor. Overall, these results indicate that hydrophobic (i.e., nonpolar) interactions are predominantly responsible for the binding of VS-4718 with human FAK.

Free Energy Decomposition

The calculated binding free energies presented in Table show that the nonpolar (i.e., hydrophobic) term plays the most important role in complex formation. Because the per-residue free energy decomposition strategy is known to enable analysis of the inhibitor–protein interactions,[55,131−134] the interaction energies between the various residues of FAK and VS-4718 were computed using the MM/GBSA decomposition protocol (Tables S5 and S6). As indicated by the obtained results, several hydrophobic residues possessed substantial subtotal binding free energies. In the FAK-II-1 system, the L501 and C502 residues present in the hinge loop made a favorable contribution to binding (i.e., −1.85 kcal/mol, Figure ) because the C502 residue can form two hydrogen bonds with VS-4718. It should be noted that the residues of the hinge loop that form hydrogen bonds with the ATP-competitive inhibitor are known to be conserved for FAK[28,31,32] and for other protein kinases.[135−140] In addition, the hydrophobic interaction between the side chain of L501 and the methoxy group of VS-4718 was the main contributor for this residue, giving a value of −1.75 kcal/mol. However, the methoxy group of the inhibitor generated an electrostatic repulsion with C502 due to the proximity of this group to the carbonyl oxygen atom of C502 (Figure S29). Variation in the position of this methoxy group has therefore been used to increase the inhibitor selectivity toward different protein kinases, such as in the case of dasatinib.[141] Therefore, the presence of a substituent at this position is of particular importance. Ideally, the methoxy group could be substituted by a halogen atom or a small alkyl chain to increase the binding affinity and prevent electrostatic repulsion with C502.

Figure 6

Key residues involved in VS-4718 binding with human FAK. The energy was decomposed into its backbone and side-chain components for each residue in the FAK-II-1 (A) and FAK-II-4 (B) systems. The energy also decomposed into its nonpolar solvation, polar solvation, electrostatic, and van der Waals components for each residue in the FAK-II-1 (C) and FAK-II-4 (D) systems. The binding energy was decomposed using the MM/GBSA method.

In addition, the L553 residue also provided a contribution of more than −2.90 kcal/mol due to a stable hydrophobic interaction with the pyridine ring of VS-4718 (Figure S30). Based on previous studies, it would be expected that substitution of the pyridine moiety with pyrimidine would reduce the distance between the ligand and the residue to increase the binding affinity[58] (Figure S31). Indeed, the pyrimidine ring has been incorporated into other FAK inhibitors, such as VS-6063 and VS-6062,[67] suggesting that the hydrophobic interaction between the VS-4718 inhibitor and the L553 residue of FAK could be vital for ligand binding. Furthermore, the I428 residue on the P-loop is of particular importance due to its ability to form hydrophobic interactions with both the pyridine ring and the methoxybenzene ring of VS-4718. In such cases, the side chain of I428 is responsible for the interaction with the methoxybenzene ring, as observed from the distance and angle of the interacting structure (Figure S32). However, the CD1 and CG1 atoms of the I418 side chain exhibit angles >90° with the methoxybenzene ring, and as previously reported, this ring could be modified with aliphatic groups, as observed for the thieno[3,2-d]pyrimidine derivatives.[125] The methoxy group points toward the hinge loop of FAK, while in contrast, the methyl group of dasatinib has been reported to result in the opposite orientation.[105] It would also be expected that if the pyridine and methoxybenzene rings were to be linked by a five- or six-membered ring system, favorable hydrophobic interaction would form with the CD1 or CG1 atoms of the I428 residue to increase the binding affinity. The report of an FAK inhibitor containing a tricyclic pyrimidothiazolodiazepinone core (i.e., BJG-03-025) also indicates that a six-membered ring and a five-membered ring can be linked by a seven-membered ring (Figure S33).[142] However, due to the nonaromatic nature of the resulting seven-membered ring, the affinity of BJG-03-025 was reduced to 20 nM,[142] thereby confirming that the presence of an aromatic ring is essential at this site. Moreover, the tricyclic benzopyrimidodiazepinone core has been shown to act as a privileged scaffold for the generation of potent and selective kinase inhibitors.[142] Therefore, it can be inferred that other tricyclic cores may be used to increase the selectivity of FAK inhibitors. Additionally, the methoxybenzene ring also formed hydrophobic interactions with the G505 residue (−1.60 kcal/mol), and it has been shown that the G505 and I428 residues can form a “clip” to bind the methoxybenzene ring (Figure S34), thereby suggesting that this moiety cannot be substituted by other nonaromatic rings.

It should also be noted here that the DFG motif serves as an important regulator of kinase activities. In addition, it is highly conserved and follows the T-loop and therefore can be used in the design of novel inhibitors. In this motif, the D564 residue forms a hydrogen bond with the O1 atom of VS-4718, which contributed an energy of −1.95 kcal/mol. Furthermore, the L567 residue on the T-loop was found to contribute an energy of −2.39 kcal/mol to the binding of VS-4718 with FAK, and this interaction can be attributed to hydrophobic binding with the benzene ring of the N-methylbenzamide moiety (Figure S35). Importantly, if the conformation of the T-loop is open, L567 will point toward the solvent environment and destroy this key hydrophobic contribution, and therefore the aromatic N-methylbenzamide ring must be retained in the structure. Furthermore, we found that the DFG-motif forms an α-helical structure in VS-4718/FAK (Figure S36), which enables multiple interactions, such as those between the D564/L567 residues and VS-4718. Indeed, this α-helical conformation has been demonstrated to provide selectivity between FAK and PYK2 (61% sequence identity with FAK in the kinase domain).[58] More specifically, the benzene ring of VS-4718 provides a favorable hydrophobic interaction with L567 to decrease the off-rate, which can be used to increase the selectivity between FAK and PYK2.[58] These results therefore suggest that a DFG-in conformation with a closed T-loop is necessary in the FAK structure to obtain selective FAK inhibitors.

With the exception of the D564 residue, we found that the key residues involved in VS-4718 binding with human FAK in the FAK-II-4 system were similar to those in the FAK-II-1 system. More specifically, D564 is unable to form a hydrogen bond with the O1 atom of VS-4718 in the FAK-II-4 system, which ultimately changes the backbone contribution from −1.17 kcal/mol in FAK-II-1 to −0.31 kcal/mol in FAK-II-4. Overall, the I428, A452, L501, C502, G505, L553, D564, and L567 residues appear to be the key residues for VS-4718 binding with human FAK.

Design Strategies

Based on the central aminopyrimidine hinge, VS-4718 was modified at four different regions. More specifically, the interactions with the solvent-exposed pocket and the interactions with the nonconserved upper lobe residues were probed (R1), and the FAK back–pocket interface containing an induced helix was modified (R2). In addition, the DFG helix induced in FAK after ligand binding was investigated (R3), while the influence of the central aminopyrimidine group was evaluated (R4) (Figure ). Thus, based on the interaction model between VS-4718 and human FAK, 47 novel compounds were designed to target the ATP-binding pocket of human FAK (Figure S37), and ensemble docking[143,144] was performed to assess the effects of these modifications on the inhibitor binding affinity.

Figure 7

Representation of the key interactions between VS-4718 and FAK, and the four modification regions. (A) Interactions between VS-4718 and the key residues of human FAK determined herein. (B) Modification of VS-4718 at four different regions (R1, R2, R3, and R4).

As shown in Figures S38 and S39, the representative conformations of the various clusters of the FAK-II-1 and FAK-II-4 systems were analyzed using the rmsd values obtained for the Cα of the protein structure, and similar conformations were obtained for two complex systems in each case. In addition, based on the cluster analysis, the first cluster of the MD simulation was determined to be the representative conformation for this MD simulation because its occupancy was >38.2% for FAK-II-1 and >70.7% for FAK-II-4. The centroid frames in each cluster were therefore selected as representative conformations for ensemble docking.

Overall, 10 conformations of the complex systems were employed in the ensemble docking calculations, and the docking scores for the novel compounds were selected from the lowest binding free energy of each docking (Table S7). As indicated, modification of the R1 region had little influence on the binding affinity, although it was necessary to retain this region to prevent loss of the inhibitor activity. In addition, certain modifications of the R2 region increased the binding affinity (e.g., in the cases of R2a, R2d, and R2g–R2k), while in the R3 region, an improved binding affinity was obtained when the N-methylbenzamide group was present (R3a–R3j). In contrast, aromatic heterocycle group replaced with the benzene ring of benzamide group (R3k–R3p) decreased the binding affinity. Finally, modification of the R4 segment appears key for optimizing the main hydrogen bonding network, and it was found that the R4d and R4e structures constituted novel frameworks to improve binding with the C502 residue.

Conclusions

The inhibition of FAK has shown potential as a therapeutic treatment of various carcinomas, such as breast, ovarian, and neck cancers. Previously, VS-4718 was reported as a selective, reversible inhibitor of FAK, with an IC50 value of 1.5 nM. Such inhibition can be attributed to a conformational rearrangement of the kinase domain of human FAK upon inhibitor binding. Therefore, we carried out a molecular modeling study of the binding mechanism between VS-4718 and FAK to probe the key interactions responsible for this action. More specifically, molecular docking studies, MD simulations, binding free energy calculations, and energy decomposition studies provided critical information regarding the molecular interactions and binding affinities within the VS-4718/FAK complexes, and a reasonable interaction model between the inhibitor and the protein was established. Overall, the obtained results indicated that VS-4718 can be modified to enhance its binding affinity with FAK based on the following strategies: (i) maintenance of the N-methylbenzamide moiety for binding with the Asp-Phe-Gly (DFG) motif of FAK, (ii) enhancement of the interactions with the back pocket of FAK, and (iii) optimization of the hydrogen-bonding interactions from the diaminopyrimidine group. Overall, the present study not only facilitates a better understanding of the binding mechanism of human FAK with VS-4718 but also provides additional insights into potential future design strategies for this inhibitor. The synthesis and biological evaluation of novel inhibitors targeting FAK will be considered in the near future, and the results will be reported in due course.