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  1), 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 Cterminal 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.
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 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  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 IC 50 (50% inhibitory concentration) value of 1.5 nM in an in vitro kinase assay. 38 The VS-4718 has also demonstrated a cellular IC 50 of ∼100 nM in malignant pleural mesothelioma cell lines and breast carcinoma. Furthermore, it has a median relative cellular IC 50 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) 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- 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: 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 AutoDock-Tools 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 equations where G complex , G protein , and G ligand are the free energies of binding for VS-4718 + FAK, FAK, and VS-4718, respectively. The value of ΔG binding can be summarized based on the enthalpic term (ΔH = ΔE gas + ΔE sol ) and the entropic term (TΔS), while the molecular mechanical energy (E gas ) can be divided into van der Waals forces (E vdW ), electrostatic forces (E ele ), and the intramolecular energy (E int ), which can themselves be obtained from a statistical average based on the ff19SB force field. In addition, the solvation free energy (E solv ) can be decomposed into electrostatic (E GB ) and nonelectrostatic (E surf ) components. More specifically, the E surf 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, E surf 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 E GB 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 2), 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.
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 3A). This observation, which is consistent with the rmsd results, indicates that the Tloop 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 saltinducible 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 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 3B,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 4). 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 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.
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 Nmethylbenzamide 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 Mg 2+ 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 5). 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 4B) 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.
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 1, S3, and S4). It was found that all entropy values (TΔS total ) and enthalpies (E gas + G sol ) 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 (ΔG bind cal ) was approximately −33.13 kcal/mol, whereas for FAK-II-4, the ΔG bind cal 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 IC 50 value of 1.5 nM. 38 The binding energy can usually be decomposed into its polar (E ele + E GB ) and nonpolar (E vdW + E surf ) 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 (E vdW ) 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 1 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 6) 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.
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 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.

ACS Omega
http://pubs.acs.org/journal/acsodf Article 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 7). Thus, based on the interaction model between VS-4718 and human FAK, 47 novel compounds were designed to target the ATPbinding 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.
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 IC 50 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.