Identification of defactinib derivatives targeting focal adhesion kinase using ensemble docking, molecular dynamics simulations and binding free energy calculations

Abstract Focal adhesion kinase (FAK) belongs to the nonreceptor tyrosine kinases, which selectively phosphorylate tyrosine residues on substrate proteins. FAK is associated with bladder, esophageal, gastric, neck, breast, ovarian and lung cancers. Thus, FAK has been considered as a potential target for tumor treatment. Currently, there are six adenosine triphosphate (ATP)-competitive FAK inhibitors tested in clinical trials but no approved inhibitors targeting FAK. Defactinib (VS-6063) is a second-generation FAK inhibitor with an IC50 of 0.6 nM. The binding model of VS-6063 with FAK may provide a reference model for developing new antitumor FAK-targeting drugs. In this study, the VS-6063/FAK binding model was constructed using ensemble docking and molecular dynamics simulations. Furthermore, the molecular mechanics/generalized Born (GB) surface area (MM/GBSA) method was employed to estimate the binding free energy between VS-6063 and FAK. The key residues involved in VS-6063/FAK binding were also determined using per-residue energy decomposition analysis. Based on the binding model, VS-6063 could be separated into seven regions to enhance its binding affinity with FAK. Meanwhile, 60 novel defactinib-based compounds were designed and verified using ensemble docking. Overall, the present study improves our understanding of the binding mechanism of human FAK with VS-6063 and provides new insights into future drug designs targeting FAK. Communicated by Ramaswamy H. Sarma

Human FAK (UniProt ID (Bateman et al., 2019;Bateman et al., 2017;Bateman et al., 2015;Wu et al., 2006): Q05397) has a molecular weight of 125 kDa and consists of three main structural domains (Figure 1): the amino-terminal 4. 1,ezrin,, which regulates protein-protein interactions such as FAK interactions with integrin (Cooper & Giancotti, 2019), p53 (Lim et al., 2008) and epithelial growth factor receptor (Shen & Guo, 2020); the protein kinase domain (KD) (residues: 422-680), which phosphorylates tyrosine residues on substrate proteins or activation loop of FAK to regulate downstream genes, such as phosphatidylinositol 3-kinase (Xia et al., 2004); and the carboxyl-terminal focal adhesion target (FAT) domain (residues: 707-1052), which is required for binding to focal adhesion proteins (Nowakowski et al., 2002), such as paxillin (Rashid et al., 2017).The kinase domain is formed by an N-terminal lobe (N-lobe) region, C-terminal lobe (C-lobe) region and a linker hinge loop between N-lobe and C-lobe.A conserved glycine-rich (GxGxUG) loop (also known as P-loop) in the N-lobe links the b1-and b2-strands and plays a key role in substrate and inhibitor binding.In addition, a salt bridge is formed between the lysine (K454 of human FAK) within the b3-strand (AxK sequence) and the glutamate (E471 of human FAK) near the middle of the aChelix, which is necessary for the expression of the full protein kinase with catalytic activity (labeled as aC-in conformation) (Roskoski, 2016).In contrast, the conformation was labeled as aC-out without the K454 À E471 salt bridge.The activation loop (also known as the T-loop or A-loop) extends from the 564 DFG 566  (Asp-Phe-Gly) motif to the 590 APE 592 (Ala-Pro-Glu) motif and serves as an important regulatory element for protein kinase activity.The activation conformation of protein kinases has considerable with DFG-in conformation, where the side chain of D564 was orientated into the adenosine 5onformation, w (ATP)-binding pocket and the side chain of F565 was orientated out (Modi & Dunbrack, 2019).Conversely, the side chain of D564 orientated out the ATP-binding pocket and side chain of F565 into was defined as DFG-out conformation.
VS-6063, a second-generation ATP-competitive inhibitor of FAK, has shown considerable antitumor activity (inhibitory concentration of 50%, IC 50 ¼ 0.6 nM) in preclinical studies (Jones et al., 2015;VS-6063 PIs, 2014;Jones et al., 2011).Moreover, VS-6063 has also shown potential at targeting proline-rich tyrosine kinase-2 (Pyk2). 47It displayed over 100-fold greater selectivity for FAK and Pyk2 than for other nonreceptor kinases ( Jones et al., 2015 ).Additionally, VS-6063 inhibits the FAK phosphorylation with IC 50 ¼ 3 nM, as per the A431 epidermoid carcinoma cell study (Jones et al., 2015).VS-6063 also blocks FAK phosphorylation in vivo (Jones et al., 2015 ).The phase I studies of VS-6063 were conducted in patients with advanced solid malignancies (NCT02913716, NCT01943292 and NCT00787033) (Jones et al., 2015;Shimizu et al., 2016).Phase II clinical trials of VS-6063 have also been completed for non-small cell lung cancer patients with KRAS mutation (NCT01951690) (Gerber et al., 2020).In addition, some other clinical trials tested VS-6063, such as for cancers with NF2 genetic changes (NCT04439331), low-grade serous ovarian cancer with or without a KRAS mutation (NCT04625270), metastatic uveal melanoma (NCT04720417), recurrent KRAS mutation and BRAF non-small cell lung cancer (NCT04620330), advanced pancreas adenocarcinoma (NCT04331041) and advanced refractory solid tumors, lymphomas, or multiple myelomas (NCT02465060).Notably, both VS-6063 and its metabolites have been found to be UGT1A1 inhibitors (Shimizu et al., 2016).However, the binding mechanism of VS-6063 with human FAK has not been explained using crystal-structure or theoretical studies.Therefore, in this context, the binding mechanism will be explained from the theoretical methods to provide a reference for the development of new antitumor FAK-targeting drugs.
A reliable ligand/protein complex structure is vital for explaining the binding mechanism between ligand and protein (Waseem et al., 2022;Sen Gupta et al., 2022;Matsukura & Miyashita, 2020).Meanwhile, the complex structure plays an important role in structure-based drug design.In addition, elucidation of inhibitor/protein binding models requires good understanding of the incorporation of the small molecule inhibitor and the macromolecule protein.Commonly, molecular docking can be applied to create a possible binding model of the inhibitor/protein complex and is also widely accepted in designing novel inhibitors targeting the specific protein (Bhagat et al., 2021;Morris & Della Corte, 2021;Jakhar et al., 2020;Chen et al., 2020).Ensemble docking (Amaro et al., 2018;Vilar & Costanzi, 2013), which selects different initial conformations of a protein as the docking receptor temples, has been employed to construct the initial binding models in inhibition and targeting studies (Shi et al., 2021;Aguayo-Ortiz et al., 2022;Yuce et al., 2022).Molecular dynamics (MD), an effective method to sample the possible binding model between ligand and protein, can help elucidate the interactions and find the best binding model to guide the drug design (Salo-Ahen et al., 2020;Sivakumar et al., 2020;Salmaso & Jacobson, 2020;Bera & Payghan, 2019;Du et al., 2020;Yang et al., 2021;Zhan et al., 2016).
Currently, there is no crystalline structure published for VS-6063 binding to FAK.Thus, the absence of a binding model hinders the development of novel potential FAK inhibitors and optimization of VS-6063 to find inhibitors with higher binding affinity than that of VS-6063.In this study, the goal primarily was to understand the VS-6063/ FAK binding mechanism.VS-6063/FAK binding models were constructed using ensemble docking and MD simulations.The binding free energies for VS-6063/FAK were estimated with the molecular mechanics/generalized Born (GB) surface area (MM/GBSA) method.The per-residue energies have been calculated using decomposition analysis to find the key residues involved in VS-6063/FAK binding.Second, 60 novel compounds based on VS-6063 were designed to target FAK.The binding affinity of the designed compounds was estimated using ensemble docking.We believe that these results could provide useful information to develop FAK inhibitors.

Molecular docking
The two-dimension structure of VS-6063 was drawn using InDraw (http://www.integle.com/static/indraw),and the three-dimension structure was generated using OpenBabel 3. 1 (O'Boyle et al., 2011).The VS-6063 structure was optimized at the semi-empirical PM3 method (Stewart, 1989) using MOPAC2016 software (Stewart Computational Chemistry, Colorado Springs, CO; http://OpenMOPAC.net).In addition, the structures of human FAK (UniProt (Bateman et al., 2019) ID: Q05397) were obtained from the Protein Data Bank (PDB) (Berman et al., 2000;Burley et al., 2021).There are more than 24 crystal structures that are known to bind with small molecules (Arold et al., 2002;Nowakowski et al., 2002;Hoellerer et al., 2003;Xia et al., 2004;Lee, 2006;2007;Garron et al., 2008;Lietha and Eck, 2008;Koolman, 2011;Gradler et al., 2013;Heinrich et al., 2013;Iwatani et al., 2013;Tomita et al., 2013;Brami-Cherrier et al., 2014;George et al., 2015;Kadare et al., 2015;Zhou et al., 2015;Yen-Pon et al., 2018;Momin et al., 2019;Popow et al., 2019;Thifault, 2020;Berger et al., 2021).The classes of small molecule protein kinase inhibitors can be defined as types I, II, III, IV, V and VI (Roskoski, 2016).The criterion for type I inhibitor was based on the conformation of protein kinase with DFG-in and aC-in.The classes of FAK inhibitors were crated with the conformation of the DFG motif and aC-helix (Figure S2).Most of the inhibitors belonged to the type I inhibitor class of FAK, such as 1MP8 (Chuang et al., 2022) (PDB ID), 2ETM (Lee, 2006), 4I4E (Tomita et al., 2013), 6I8Z (Mousson et al., 2021), and 6YVY (Berger et al., 2021).In addition, less type IV inhibitors were discovered, including 4EBV (Iwatani et al., 2013), 4EBW (Iwatani et al., 2013) , and 4I4F (Tomita et al., 2013).Type II inhibitors included 4K9Y (Gradler et al., 2013), 4KAO (Gradler et al., 2013) and 4Q9S (George et al., 2015).From the crystal structures of FAK bound with a type I inhibitor (Figure S3), two different states of the kinase domain of human FAK were identified (FAK-I and FAK-II).FAK-I was defined with an open T-loop, aC-in, and DFG-in conformation, such as in PDB ID: 6YVY (Berger et al., 2021).Meanwhile, FAK-II was defined with a closed T-loop, aC-in, and the DFG-in conformation, such as in PDB ID: 6I8Z (Popow et al., 2019).Based on previous studies, VS-6063 was found to be an ATP-competitive inhibitor of FAK and serves as a type I inhibitor binding with maternal embryonic leucine zipper kinase (PDB ID: 5MAH) (Jones et al., 2011;Klaeger et al., 2017).VS-6063 inhibits the activity of MELK with IC50 ¼ 1964 964ibits the activity of MELK with IC50 ¼ or of FAK and serves as a type I inhibitor binding with materFigure S4).Meanwhile, the similarity inhibitors of FAK, such as (E)-N-methyl-N- (3-(((2-(phenylamino)-5-(trifluoromethyl)pyrimidin-4-yl)imino)methyl)pyridin-2-yl)methanesulfonamide (Figure S5), were also type I inhibitors with DFG-in and aC-in conformation of FAK.Thus, we can speculate that VS-6063 is a type I inhibitor for FAK and binds with human FAK with the DFG-in and aC-in conformation.Furthermore, the two conformations of T-loop in type I inhibitors bind with human FAK.Based on the crystal structures for human FAK binding with type I inhibitors, the conformation of T-loop can be found in an open or closed conformation.Currently, there is no crystalline structure published for VS-6063 binding to FAK.Additionally, the T-loop conformation for VS-6063 binding with FAK is also unclear; therefore, both the open and T-closed conformation need to be considered in this study.As such, VS-6063 binding with an open or closed conformation of T-loop was employed to construct the complex structures of VS-6063/FAK in this study.To prepare the receptor structure, PYMOL 2.1 (Schr€ odinger, 2010) was used to remove the co-crystallized ligand molecules and crystallographic water molecules.Subsequently, FAK and VS-6063 were pretreated by adding hydrogen atoms and Gasteiger charging (Gasteiger & Marsili, 1980) and adjusting unreasonable atomic overlap using AutoDockTools 1.5.6 (Sanner, 1999).A grid box with a 40 � 40 � 40 grid size and 0.375 Å grid spacing was estimated with AutoGrid v.4.2 (Morris et al., 2009).Meanwhile, the center of the grid box was defined from the center of the native ligand molecule in the crystal complex.Twenty-hundred conformations were obtained using the Lamarckian genetic algorithm (Fuhrmann et al., 2010) per system in AutoDock v.4.2 (Morris et al., 2009).To identify the best initial binding models, the conformation with the best rational orientation in the active pocket, based on the published crystal complex structures for the type I ATP-competitive inhibitor with human FAK (such as the key hydrogen bonds with the hinge loop of FAK), was selected as the best conformation from each of the docking models from the docking experiment.

MD simulations
Detailed information on the binding mechanism of VS-6063/ FAK was explored using MD simulations.The docking models of the VS-6063/FAK complex structures were applied as the possible initial binding models for the corresponding MD simulations.The standard Amber ff19SB force field (Tian et al., 2020) has been used on human FAK.However, there is no standard force field for the small molecule VS-6063.Therefore, two steps of the general amber force field (version 2; GAFF2) (Wang et al., 2004) procedure were used to generate the force field of VS-6063.In the first step, the structure of VS-6063 was optimized using Gaussian09 program (Frisch et al., 2009) with B3LYP/6-31G* theory level.In the second step, the restrained electrostatic potential (RESP) protocol (Bayly et al., 1993) was performed to fit the partial atomic charges of VS-6063.Subsequently, the VS-6063/FAK complex systems were solvated in a cuboid box (approximately 87 Å � was � was performed to fit th � was � was performed to fit the (Jorgensen et al., 1983) type water, and the distance between the complex systems and solvate edge was more than 15 Å.To maintain the electrically neutral state of the solvated system, sodium chloride was added in the solvated box.Finally, the complex system included human FAK, the inhibitor VS-6063, solvated water, and sodium chloride.
The steepest descent method and conjugate gradient method were applied to optimize the VS-6063/FAK system and avoid the formation of unfavorable interactions due to solvents and ions addition.First, the water molecules and counterions of the VS-6063/FAK system were optimized by restraining FAK and VS-6063, and the weight for the positional restraints was 2.0 kcal/mol/Å (Chuang et al., 2022).Subsequently, the overall complex system was minimized without constraints.For both instances, the 8000 steps of the steepest descent algorithm and 2000 steps of the conjugate gradient algorithm were employed.Next, the temperature was set to 300 K with Langevin dynamics (Feller et al., 1995;Martyna et al., 1994) within 200 ps, and the collision frequency was set to 2.0 ps À .The system pressure was maintained at 1 bar using the isotropic position scaling method for 200 ps.Subsequently, the 200-ps simulation was applied to equalize the system at 300 K and 1 bar in the NPT ensemble (isothermal-isobaric ensemble).To collect the data for analysis, an additional 500-ns simulation was performed for every system in the same situation with an equilibrium step.The coordinate trajectories were saved every 10 ps (5000 steps) for every simulation.
In these simulations, periodic boundary conditions were applied to treat unphysical edge effects.The SHAKE algorithm (Ryckaert et al., 1977) was used to manage the hydrogen atoms of the protein and ligand molecules related to covalent bonds, allowing an integration time step of 2 fs.The particle mesh Ewald algorithm (PME) (Darden et al., 1993) was employed to decrease the long-range electrostatic force.A cutoff value of 12 Å was used to address the shortrange electrostatic force and the van der Waals interactions.The seed for a pseudorandom number was generated on the clock.The system was prepared using the tleap model in AmberTools21 (Case et al., 2020).The analytical process from the MD trajectories was conducted using the CPPTRAJ (Roe & Cheatham, 2018;Roe & Cheatham, 2013) model.However, MD simulations were performed from the PMEMD in AMBER20 (Case et al., 2020).The CUDA version of PMEMD was employed to decrease the time based on the Nvidia GPUs.For every system, the simulation was replicated thrice with the same initial conformation.

Binding free energy calculation
In addition to qualitative analysis of ligand/protein binding, quantitative analysis is also important in determining the binding affinity between ligands and proteins.Several methods, such as molecular mechanics/Poisson Boltzmann (or generalized Born) surface area (MM/PBSA or MM/GBSA) (Srinivasan et al., 1998;Lee et al., 2004), solvated interaction energy (Naim et al., 2007), linear interaction energy (Perdih et al., 2009;Bren et al., 2006;Aqvist et al., 1994), free energy pathway method (Gilson & Zhou, 2007), and linear response approximation (Lee et al., 1992), have been developed to estimate the absolute binding free energy between inhibitors and their target proteins.The MM/GBSA approach is a timesaving and efficient method to evaluate the binding free energy between inhibitors and proteins (Shi et al., 2021;King et al., 2021;Cheng et al., 2018;Shirvani & Fassihi, 2022;Shi et al., 2022).In this study, only a short description of the calculation of the binding free energy (DG binding ) from the MM/ GBSA method is provided in the following formulae: (1) G complex , G protein , and G ligand denote the free energies of VS-  6063/FAK, FAK, and VS-6063, respectively.G can be decomposed into enthalpy (H ¼ E gas þ E sol ) and entropy (TS).
The molecular mechanical energies (E gas ) can also be summarized from the intramolecular energy (E int ), van der Waals forces (E vdW ), and electrostatic forces (E ele ).Meanwhile, the contributions of E int , E vdW , and E ele can be obtained through the statistical average based on molecular mechanics.In addition, solvation free energy (E sol ) can be divided into polar solvation (E polar ) and nonpolar solvation energies (E nonpolar ).E nonpolar is obtained from the favorable van der Waals interactions between the solute and solvent and the unfavorable cost of surface formation.E nonpolar can be calculated using equation ( 5), where c ¼ 0.0072 072.re 2 and b ¼ 0.0 kcal/mol.The linear combination of pairwise overlaps (LCPO) method (Weiser et al., 1999) was employed to estimate the solvent accessible surface area (SA).However, the E polar contribution was calculated from the GB equation (Still et al., 1990;Srinivasan et al., 1999).The dielectric constants for the ligand/protein and water were set to 1 and 80, respectively.In addition, normal model analysis was used to calculate the entropy contribution.Thousand snapshots were extracted from the last 200 ns MD trajectory to calculate the statistical average of the MM/GBSA method.The 100 snapshots from the last 200 ns MD trajectory were employed to estimate the entropy contribution.
The binding free energy between VS-6063 and each residue of human FAK was also decomposed for van der Waals (DG vdW ), electrostatic (DG ele ), polar solvation (DG polar ), and nonpolar solvation energies (DG nonpolar )-but not for entropy-using the MM/GBSA method, and the same parameters were applied in the binding free energy calculation.In addition, the free energy decomposition was calculated for the backbone (BDG subtotal ) and sidechain energies (SDG subtotal ) for each residue.All methods employed during this study are outlined in Figure S6.

Initial models
Generally, the type I protein kinase inhibitor binds in the ATP-binding pocket of an active conformation of protein kinase with DFG-in and aC-in (Roskoski, 2016).The type II  inhibitor binds in the ATP-binding site of an inactive conformation of protein kinase with DFG-out and aC-helix as either aC-in or aC-out.The inhibitors of the human FAK have been classified according to the conformation of DFG motif, aC-helix, and binding site (Figure S2).Most of the inhibitors were type I with DFG-in and aC-in conformation.In previous studies, VS-6063 was identified as a type I inhibitor of MELK (PDB ID: 5MAH) (Klaeger et al., 2017).Meanwhile, the similarity inhibitors of FAK were also classified as type I inhibitors with DFG-in and aC-in conformation of FAK.Thus, it can be speculated that VS-6063 also binds to human FAK as a type I inhibitor with DFG-in and aC-in.The docking power was referred to as the ability of a molecular docking strategy to identify the native binding conformation from the crystal structure (Sanner et al., 2021;Su et al., 2019).In this work, the 'redocking' strategy was employed to evaluate the docking power for the selected parameters, which was applied to construct the VS-6063/FAK complex structures.The crystal complex structures for 4-((4-(((1 R,2R)-2-(dimethylamino)cyclopentyl)amino)-5- (trifluoromethyl)pyrimidin-2-yl)amino)-N-methylbenzenesulfonamide/FAK (FAK-I) and BI-4464/FAK (FAK-II) were then redocked.The root mean square deviation (RMSD) values were 0.87 and 0.97 Å for FAK-I and FAK-II, respectively, obtained between the crystal structure conformation of the ligand and the conformation with the lowest energy (i.e.À 7.88 and À 11.9 kcal/mol for FAK-I and FAK-II, respectively) based on the docking results (Figure S8).Therefore, the same docking strategy was used to construct the VS-6063/ FAK-I and VS-6063/FAK-II complexes.
The possible VS-6063/FAK binding models were obtained from the docking experiments.A total of 2000 conformations were obtained from molecular docking, which were clustered based on the RMSD values with 2.0 Å.For VS-6063/FAK-I, the first cluster with 466 conformations had 23.3% occupancy, and the second cluster had 48.6% occupancy (Figure S9).In addition, the lowest binding energies for those two clusters were similar (-8.02 and À 8.00 kcal/mol for the first and second clusters, respectively).Thus, it can be speculated that these two clusters may act as possible binding models.When the conformation was checked, the main binding pattern was similar between the two clusters, excluding the ethylmethylsulfonamido group.In the following simulation, the lowest binding energy in the first cluster was defined as the possible binding model for VS-6063 with FAK-I.In contrast, the occupancy of the first cluster with VS-6063 binding FAK-II was 49%, approximately half of that of the 2000 conformations (Figure S10).Thus, the conformation with the lowest binding energy in cluster 1 (-10.02kcal/mol) was selected as the possible binding model for VS-6063/FAK-II.Two possible binding models (Figure 2), which were studied in the following simulations, were identified using molecular docking experiments.

System stability
The docking experiments did not consider the additional interactions between VS-6063 and FAK or the rearrangement of the residues in the active site of FAK.Thus, MD simulations were required to obtain additional information related to the binding mechanism of the VS-6063/FAK complex system.During this process, MD simulations for each VS-6063/FAK system over 500 ns were performed.The RMSD values of the heavy atoms of the FAK backbone (C, Ca, O, and N atoms of every residue) and of the nonhydrogen atoms of the VS-6063 were equilibrated after 100 ns and maintained at a plateau until the end of the simulations (Figure 3).Fluctuations of the RMSD values for the protein kinase were observed, with values of 1.91 ± 0.28 and 1.61 ± 0.20 Å for FAK-I-1 and FAK-II-1, respectively (Figure S11).This result indicates that the FAK structures were stable.In contrast, the fluctuations of VS-6063 were 3.18 ± 0.43 and 2.01 ± 0.25 Å for FAK-I-1 and FAK-II-1, respectively.Meanwhile, the simulations for the FAK/VS-6063 complex systems were rerun for three replicates (FAK-I-1, FAK-I-2, and FAK-I-3 for the three replicates of FAK-I; FAK-II-1, FAK-II-2, and FAK-II-3 for the three replicates of FAK-II), which also showed similar results (Table S1).In addition, the radius of gyration of the VS-6063/FAK complex was also calculated to assess the stability when VS-6063 was bound to FAK and indicated that those systems were stable in the last 200 ns (Figure S12).Therefore, these simulations can be used to analyze the interaction between VS-6063 and FAK.
The root-mean-square fluctuation (RMSF) was used to quantify the flexibility of specific residues during the MD simulations.Based on the 500-ns MD trajectory, the RMSF values of all residues were calculated for every VS-6063/FAK complex system (Figure 3C).The amino acids at 564-592 (Tloop) for FAK showed more flexibility than did those of other regions for FAK-I and FAK-II with higher RMSF values.Meanwhile, a similar RMSF fluctuation trend was found for VS-6063 bound with FAK-I and FAK-II.This observation, which is consistent with the RMSD analysis, indicates that the T-loop of FAK is unstable during the simulations (Figure S13).In addition, the relatively flexible T-loop conformations were consistent with those of other kinases, such as saltinducible kinase 2 (SIK2) (Shi et al., 2021;Shi et al., 2021;Shi et al., 2022), microtubule affinity-regulating kinase 2 (MARK2) (Ahrari et al., 2017), and MARK4 (Ahrari et al., 2020).In this study, the initial conformation of FAK-II (closed) and FAK-I (open) were able to form binding models with VS-6063.Moreover, the T-loops exhibited open or closed conformations for FAK that were able to recognize various inhibitors.In addition, other protein kinases also had a similar recognition mechanism, such as the open conformation of the Tloop for Abelson leukemia virus tyrosine kinase 1 (ABL1) to bind to dasatinib (Tokarski et al., 2006) and the closed conformation to bind to imatinib (Cowan-Jacob et al., 2007).The hinge loop was very stable in the simulation for both FAK-I and FAK-II systems.In the contrary, the P-loop and aC were more flexible for FAK-I than for FAK-II.Excluding the T-loop region, the loop regions for b1-b2, b2-b3, b3-b4, and b4-b5 had high values of RMSF, indicating that those regions proved flexible in the simulation (Figure S14).Meanwhile, FAK-II was more stable than FAK-I, especially for the N-lobe regions.However, the flexible regions of FAK were not clearly changed when VS-6063 was bound (Figure S14).Overall, these observations suggest that the protein conformations of FAK should be stable, excluding the T-loop, when VS-6063 binds to the ATP-binding pocket.
The snapshots, which were extracted from 100, 200, 300, 400, and 500 ns, were aligned and depicted for FAK-I (Figures S15-S17) and FAK-II (Figures S18-S20).The conformations of the initial and last 500-ns MD trajectories for the two complex systems were constructed to check the possible changes from initial docking models and the refined binding model from the MD (Figures 3 and S21).FAK in the simulations adapted almost the same conformation as its initial conformation, and the binding modes of VS-6063 were nearly similar to the initial structure from the docking experiments, excluding the methanesulfonamide group.Meanwhile, the activation loop (T-loop) was changed with an induced-fit with inhibitors, which can also be observed from the RMSD value of the T-loop along with the simulation times (Figure S13).However, the T-loop maintained the closed or open conformations, especially the DFG domain (564-566), which can also be observed in the RMSF values (such as 0.65, 0.79, and 1.01 for D564, F565, and G566 of FAK-I-1, respectively).Thus, VS-6063 maintains binding to the ATP-binding site of FAK with the open and closed conformations of the T-loop.

Hydrogen bond analyses
Some potential hydrogen bonds between FAK and VS-6063 had been identified from the initial docking models.However, only a few hydrogen bonds can be found between the inhibitor VS-6063 and FAK protein in the simulations (Figure 4), with an average number of hydrogen bonds of 2.80, 2.63, and 2.64 for FAK-I-1, FAK-I-2, and FAK-I-3, respectively.Meanwhile, there were more than three hydrogen bonds for the FAK-II system, with 3.33, 3.54, and 3.27 hydrogen bonds for FAK-II-1, FAK-II-2, and FAK-II-3, respectively.Therefore, more hydrogen bonds can be formed between VS-6063 and FAK for the closed conformation than for the open conformation of the T-loop.The occupancy of the hydrogen bonds was also assessed in the 500-ns simulation for the VS-6063/FAK complex systems.The N4 atom of VS-6063 acted as a donor atom to form a hydrogen bond with the oxygen atom of the C502 residue of FAK (over 93.00% occupancy excluding FAK-II-3 with 87.07%).The nitrogen atom of the same residue (C502) of FAK also maintained hydrogen bonds with the N6 atom of VS-6063 as an acceptor and 93.00% occupancy of the VS-6063/FAK systems.The two hydrogen bonds formed between FAK and VS-6063 are located within the diamino-pyrimidine group, which is a universal core for FAK inhibitors (Chen et al., 2021;Berger et al., 2021;Qi et al., 2021) and other protein kinase inhibitors to form hydrogen bonds with the hinge loop (Sun et al., 2020;Saha et al., 2021).Meanwhile, the 7H-pyrrolo 2,3-d pyrimidines (Choi et al., 2006;Choi et al., 2006), diamino-pyridine (Mousson et al., 2018;Lu & Sun, 2020;Shah et al., 2014), and thieno 3,2-d pyrimidine (Cho et al., 2021) rings also contribute to the binding model in a similar manner to the diaminopyrimidine ring, thereby indicating that the diaminopyrimidine ring can be replaced with a pyrazole or pyridine ring.
The hydrogen bonds were then examined further, wherein the distances between the donor and acceptor atoms and the angles with the donor, hydrogen, and acceptor atoms were calculated from MD simulations (Table S2).More specifically, the hydrogen bond between the oxygen atom of C502 and N4 of VS-6063 in the FAK-I-1 system was stronger (distance: 2.93 ± 0.16 Å) than that between the nitrogen atom of C502 and the N6 atom of VS-6063 (distance: 3.13 ± 0.18 Å) (Figure S22).Meanwhile, a similar trend was also observed for the hydrogen bonds for the FAK-I-2 and FAK-I-3 systems (Figures S23 and S24).However, the two hydrogen bonds were similarly strong for the FAK-II systems; for instance, the distances for C502_O-LIG_N4 and LIG_N6-C502_N were 3.3.04± 0.22 and 3.09 ± 0.15 Å, respectively, for FAK-II-1 (Figure S25).This was also found for the FAK-II-2 and FAK-II-3 systems (Figures S26 and S27), indicating that these two hydrogen bonds play the same key role in VS-6063 binding to FAK.Therefore, the hydrogen bonds between the pyrimidine group of VS-6063 and the hinge loop of FAK play an important role in orientating VS-6063 into the ATP-binding site.

Binding free energy
As described above, the interactions between VS-6063 and FAK can be depicted using docking experiments and MD simulations.However, it was not possible to identify the binding affinity of VS-6063 to FAK.Thus, MM/GBSA was employed to calculate the absolute VS-6063/FAK binding free energies for the two systems (Tables 1 and S3-S8).All entropy values (DTS total ) and enthalpies (Dð gas þ Da sol ) were negative (less than À 17.22 and À 44.53 kcal/mol, respectively), indicating that the formation of a binding complex is an enthalpy-driven process.For FAK-I, the calculated binding free energies (DG cal bind ) were approximately À 27.30, À 23.15, and À 25.17 kcal/mol for FAK-I-1, FAK-I-2, and FAK-I-3, respectively.Meanwhile, the DG cal bind values for FAK-II-1, FAK-II-2, and FAK-II-3 were À 29.03, À 30.56, and À 27.28 kcal/mol, respectively.The largest binding energy for the FAK-II system (-27.28kcal/mol) was similar to the lowest binding energy for the FAK-I systems.In addition, the average binding free energy of FAK-II was decreased to À 28.96 kcal/mol compared with the À 25.21 kcal/mol of FAK-I.Therefore, the closed conformation of the T-loop can be considered the preferred model when VS-6063 binds to human FAK.The FAK-II-2  The final estimated binding free energy from DE gas þ DE sol -DTS total : � The standard errors for all of terms are included in the parentheses, which were calculated as the root mean square error for all of frames extracted in the MM/GBSA calculation.
system has the lowest binding free energy, based on our calculation.The closed conformation of the T-loop could offer the DFG motif interaction with the side chain of VS-6063.The affinity of VS-6063/FAK was also found to be strong, consistent with the experimental results.
The binding energy can usually be decomposed into its polar (E ele þE polar ) and nonpolar (E vdW þ E nonpolar Þ terms.More specifically, the polar terms were 14.18 kcal/mol for FAK-II-2 and over 12.50 kcal/mol for the other systems.The positive values for the polar contribution indicate that the polar interactions between VS-6063 and FAK are antagonistic to this binding.Meanwhile, the polar terms for FAK-II were larger than those for FAK-I, indicating that they exert a more severe effect on FAK-II systems than on FAK-I systems.In contrast, the nonpolar values were less than À 67.65 kcal/mol for FAK-II systems and less than À 53.89 kcal/mol for FAK-I systems.This negative value of nonpolar terms shows that the nonpolar terms can be seen as the main contributor to VS-6063 binding to human FAK.In other words, the polar term is unfavorable for VS-6063 binding to human FAK, whereas the nonpolar term is favorable.The van der Waals (E vdW ) interactions (-50.40, À 52.29, À 47.70, À 60.16, À 60.72, and À 60.79 kcal/mol for FAK-I-1, FAK-I-2, FAK-I-3, FAK-II-1, FAK-II-2, and FAK-II-3, respectively) are the main nonpolar contributors and are conducive to binding with more than 110.00% of enthalpies.Hence, van der Waals interactions are the predominant factor responsible for VS-6063 binding to human FAK.The difference in the binding free energy of VS-6063 bound to FAK between the FAK-I and FAK-II systems was obtained from the van der Waals interaction with 10.43 kcal/mol.Conversely, the electrostatic interaction offered À 4.23 kcal/mol for binding.Overall, van der Waals interactions play a key role in VS-6063 binding to human FAK.This result also suggests that it is essential to focus on van der Waals interactions to improve the inhibitory potency of FAK inhibitors.

Free energy decomposition
The binding model can be explained from the MD simulation.However, the key residues for VS-6063 bound to human FAK were unclear.Meanwhile, the per-residue free energy decomposition strategy can be applied to analyze the key residues for inhibitor-protein interactions.Thus, in this study, the interaction energies between each residue of human FAK and the small inhibitor VS-6063 were estimated using the MM/GBSA decomposition protocol (Tables S9-S13).
In addition, L553 also provides a contribution of more than À 2.60 kcal/mol for VS-6063 binding to FAK.L553 contributed more energy to FAK-II systems (-2.98,À 2.94, and À 2.97 kcal/mol for FAK-II-1, FAK-II-2, and FAK-II-3, respectively) than to FAK-I systems (-2.67,À 2.77, and À 2.42 kcal/ mol for FAK-I-1, FAK-I-2, and FAK-I-3, respectively).L553 formed stable hydrophobic interactions with the pyridine ring of VS-6063 from the distance between the side chain of L553 and the center of the pyridine ring (Figure S28).The pyrimidine ring has also been widely used for the development of FAK inhibitors, such as CEP-37440 and VS-6062.Therefore, the hydrophobic interaction between VS-6063 and L553 of FAK could be vital for ligand binding.Meanwhile, the pyrimidine ring also forms interactions with A452 of FAK, which contributed more than À 1.00 kcal/mol to VS-6063 binding to FAK.A452 and L553 clipped the pyrimidine ring of VS-6063 to increase the interaction that can be found from the angle among A452, the pyrimidine ring, and L553 (Figure S29).
I428, the most favorable residue on the P-loop, can form hydrophobic interactions with the benzene ring of VS-6063.The side chain of I428 mainly contributes to the interaction with the benzene ring from the closer distance and the vertical angle (Figure S30).In addition, the benzene ring also formed hydrophobic interactions with G505 with À 1.60, À 1.40, À 2.04, À 1.41, À 1.38, and À 1.31 kcal/mol for FAK-I-1, FAK-I-2, FAK-I-3, FAK-II-1, FAK-II-2, and FAK-II-3, respectively.G505 and I428 can form a clip for the benzene ring from the angle between the side chain of G505, the center of the benzene ring, and the side chain of I428 (Figure S31).This indicates that this benzene ring cannot be substituted by other nonaromatic rings.
The DFG motif is highly conserved and follows the T-loop, which can be used to design inhibitors.Meanwhile, the DFG motif acts as an important regulator of kinase activities.D564 was found to form a hydrogen bond with the O1 atom of VS-6063, which contributed À 1.76, À 1.79, and À 1.75 kcal/ mol for FAK-II-1, FAK-II-2, and FAK-II-3, respectively; however, this hydrogen bond formed between D564 and VS-6063 could not be found for FAK-I systems.D564 interacted with trifluoromethyl of VS-6063 through halogen bonds, as widely applied in the medicinal chemistry field (Shinada et al., 2019;Margiotta et al., 2020;Riel et al., 2019;Xu et al., 2011).
Residue L567 on the T-loop also contributed more than À 2.51 kcal/mol to FAK-II systems.This interaction mainly results from the hydrophobic interaction from the pyrazine ring, which was obtained from the distance between L567 and the pyrazine ring (Figure S32).If the conformation of the T-loop is open, L567 will move toward the solvent environment, destroying the key hydrophobic contribution.Thus, the aromatic ring pyrazine needs to remain intact.

Design of novel FAK inhibitors
Based on the MD simulation and binding free energy calculation, FAK-II can be seen as the better binding model than FAK-I for VS-6063 bound with kinase domain of human FAK.Here, the 500 ns frame of FAK-II-1 was selected as the model for analysis.The FAK-II-1 was compared with the published crystal structures for ligand binding with human FAK (Figure S33).This indicates that the obtained binding model is reasonable.The hinge loop of FAK interacts with the aminopyrimidine group of VS-6063, also found in the crystal complex structures and defined as R3, to form two hydrogen bonds (Figure 6).The interactions of the regions with the solventexposed pocket and the non-conserved upper lobe residues were probed as R1 for the amide group and R2 for the benzene group.The amide group R1 mimics the similarity site with amide group of 6YVY.Meanwhile, the benzene group R2 is located at a site similar to that of a benzene ring or substituted benzene ring in crystal complex structures such as 6YVS.Meanwhile, the FAK back-pocket interface, universally occupied by trifluoromethyl in the crystal complexes, was modified as R4.The linker between the aminopyrimidine group and pyrazine ring was modified as R5, influencing the flexibility of the pyrazine group of VS-6063.Additionally, the DFG helix induced in FAK after ligand binding was investigated as R6 for the pyrazine group and R7 for the sulfamide group.The binding model of VS-6063 with FAK was similar to the binding models for aminopyrimidine-based FAK inhibitors (Berger et al., 2021).From the interaction model between VS-6063 and human FAK, 60 novel compounds (Figure S34) were designed to target the ATP-binding pocket of human FAK, and ensemble docking (Amaro et al., 2018;Vilar & Costanzi, 2013) was performed to assess the effects of these modifications on the inhibitor binding affinity.
Due to VS-6063 binding to human FAK with a closed and open T-loop conformation, the ensemble docking method for the design of the inhibitors, which was based on VS-6063, also needed to consider the closed and open conformation.Thus, the representative conformations of the FAK-I and FAK-II systems were selected from the last frame for every system (six frames in total) and applied for ensemble docking.The RMSD among the three frames with 3 Å indicates that the conformational difference was minor between FAK-I and FAK-II (Figure S35).In contrast, major difference between FAK-I and FAK-II from the T-loop conformation with open or closed was observed.In spite of the small difference among FAK-I-1, FAK-I-2, and FAK-I-3 (or among FAK-II-1, FAK-II-2, and FAK-II-3), the six frames were selected as the ensemble docking conformations to find better binding model during novel inhibitors designing.Extended details on the ensemble docking procedure are shown in the Supporting Materials.The power of the docking process was based on the RMSD between the native conformation and the docking conformation of VS-6063.This approach confirmed that the docking process was robust for running ensemble docking for the design of inhibitors (Figure S36).The docking scores for the novel compounds were selected from the lowest binding free energy of each docking (Table S14).The lowest energy for every docking run was compared within the six ensemble dockings for every designed compound, and the lowest binding energy was selected as the score for the compound.The score for VS-6063 (-9.74 kcal/mol), which also re-docked with the six conformations of FAK, was used as a standard to check whether the potential of the novel compounds in targeting human FAK was higher than that of VS-6063.Most scores were similar to that of VS-6063, and few compounds had better scores than that of VS-6063.As per the modification of R2, three design-compounds showed less binding affinity to inhibit FAK activity than VS-6063.The best scores for the other six compound series are shown in Figure 6.As indicated, loss of the R1 region decreased the binding affinity.Thus, it was necessary to retain this region.In addition, certain modifications of the R3 region to keep the key hydrogen bonds to the hinge loop did not affect the binding affinity, which indicates that scaffold hopping can be applied for the R3 region.The R4 region was tolerant of similar groups, such as the four-membered ring cyclopropane, and methyl group.The activity of linker R5 can be decreased using a methylene group but not an amide group.The R6 aromatic ring must remain intact to increase the binding affinity.When a 2H-pyrrol-2-one ring was fused to the pyrazine ring (R7), the binding affinity was enhanced.Conversely, the 2H-pyrrol or 2H-pyrazole rings decrease the binding affinity.Therefore, we inferred that the electronegativity groups need to remain in the R7 region.
The representative inhibitors, including R1d, R2c, R3c, R4f, R5f, R6a, and R7g, binding with human FAK (FAK-II-1 was selected as the initial model) have a similar binding model to that of VS-6063/FAK (Figure S37).Those seven inhibitor/ FAK systems were run using one hundred-nanosecond simulations to check the binding free energies.Generally, the inhibitor/FAK systems were stable after 50 ns from the RMSD vs.Time (Figure S38).Close inspection of the trajectories extracted from 10,20,30,40,50,60,70,80,90, and 100 ns for all seven complex systems revealed that designed VS-6063 derivatives were tightly bound (Figures S39-S45), excluding the R7g/FAK-II-1 system where 2-methylisoindolin-1-one group translated to form more stable interactions with the DFG-motif of FAK.In addition, the binding free energies from the last 20 ns simulation for the seven complex systems were also calculated (Figure 6 and Tables S15-S21).R1d, R2c, R4f, and R6a showed more affinity to bind with FAK-II-1 compared with VS-6063.Nonetheless, the activity of those designed inhibitors has not been validated.However, those designed small molecules were based on VS-6063 and are very similar to VS-6063.Ultimately, the activities of these designed molecules require further testing, and research on more diverse molecules with similar potential must be conducted.

Conclusions
FAK has shown potential as a therapeutic target for various cancers, such as ovarian cancer, breast cancer, and neck cancer.Previously, VS-6063 was reported as a second-generation reversible inhibitor of FAK and inhibited the activity of FAK with 0.6 nM (IC 50 ).However, the binding mechanism of VS-6063 and human FAK was unclear.Thus, we carried out a molecular modeling study of the binding mechanism between VS-6063 and FAK to probe the key interactions and designed novel VS-6063 derivatives.More specifically, ensemble docking, MD simulations, binding free energy calculations, and energy decomposition methods provided critical information regarding the molecular interactions and binding affinities within the VS-6063/FAK complexes, and a reasonable interaction model between the inhibitor and the protein was established.Overall, the obtained results indicated that VS-6063 can be modified to enhance its binding affinity to FAK, and 60 novel compounds were designed.Some potential compounds with better affinity than VS-6063 were identified using ensemble docking.Overall, the present study not only provides a better understanding of the binding mechanism of human FAK to VS-6063 but also provides additional insights into potential future design strategies for this inhibitor.The synthesis and biological evaluation of novel FAK-targeting inhibitors will be further considered in our future work.

Figure 2 .
Figure 2. Initial VS-6063/focal adhesion kinase (FAK) binding models identified using molecular docking experiments.(A) Binding model for FAK-I, (B) binding model for FAK-II, (C) An overlay of the binding models for FAK-I and FAK-II, and (D) T-loop conformation for FAK-I and FAK-II.Human FAK protein is represented in its cartoon format.VS-6063 ligand is shown in stick format.

Figure 3 .
Figure 3. Fluctuations in the focal adhesion kinase (FAK) conformation upon binding to VS-6063.(A) Root-mean-square deviation (RMSD) values of heavy atoms of the backbone for the protein receptor along 500 ns molecular dynamics (MD) simulation for the VS-6063/FAK complex systems.(B) RMSD values of nonhydrogen atoms of the ligand along the 500 ns MD simulation for VS-6063/FAK complex systems.(C) Root-mean-square fluctuation (RMSF) for the FAK backbone residues during the 500 ns simulation.Frames of the VS-6063/FAK complexes at the initial conformation (green) and the last (500 ns) frame (cyan) for the (D) FAK-I-1 and (E) FAK-II-1 systems.The human FAK protein is represented in cartoon form, while the VS-6063 ligand is shown in stick format.T-loop, P-loop, hinge loop, and aC helix domains of human FAK are highlighted.

Figure 4 .
Figure 4. Hydrogen bond analysis for the VS-6063/FAK systems.(A) Distribution of the number of hydrogen bonds for 50000 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 � .Additionally, schematic diagrams of hydrogen bonding are shown in the (C) FAK-II-1 and (D) FAK-II-4 systems.

Figure 5 .
Figure 5. Key residues for VS-6063 binding to focal adhesion kinase (FAK).Key residues for VS-6063 binding with FAK shown with the representative frame from the last one of the (A) FAK-I-1 and (B) FAK-II-1 systems.(C) Binding energy of the key residues from the binding energy decomposition for six VS-6063/FAK systems.Key residues are defined as those that contribute more than -0.50 kcal/mol to VS-6063 binding with human FAK in this study.

Figure 6 .
Figure 6.Potential design of focal adhesion kinase (FAK) inhibitors based on the binding mechanism between VS-6063 and human FAK.
DE vdW : Contribution to the free energy of binding from the van der Waals energy.DE ele : Contribution to the free energy of binding from the electrostatic energy.DE polar : Contribution to the free energy of binding from the polar solvation energies.DE nonpolar : Contribution to the free energy of binding from the nonpolar solvation energies.DE gas : Contribution to the free energy of binding from DE vdW þ DE ele : DE solv : Contribution to the free energy of binding from solvation energies DE polar þ DE nonpolar : DTS total : Contribution to the free energy of binding from the entropy energy.DG cal bind :