10.8
CiteScore
 
5.3
Impact Factor
Generic selectors
Exact matches only
Search in title
Search in content
Post Type Selectors
Search in posts
Search in pages
Filter by Categories
Corrigendum
Current Issue
Editorial
Erratum
Full Length Article
Full lenth article
Letter to Editor
Original Article
Research article
Retraction notice
Review
Review Article
SPECIAL ISSUE: ENVIRONMENTAL CHEMISTRY
10.8
CiteScore
5.3
Impact Factor
Generic selectors
Exact matches only
Search in title
Search in content
Post Type Selectors
Search in posts
Search in pages
Filter by Categories
Corrigendum
Current Issue
Editorial
Erratum
Full Length Article
Full lenth article
Letter to Editor
Original Article
Research article
Retraction notice
Review
Review Article
SPECIAL ISSUE: ENVIRONMENTAL CHEMISTRY
View/Download PDF

Translate this page into:

Original article
2021
:14;
202106
doi:
10.1016/j.arabjc.2021.103144

Molecular description of pyrimidine-based inhibitors with activity against FAK combining 3D-QSAR analysis, molecular docking and molecular dynamics

School of Life Science, Linyi University, Linyi 276000, China
Warshel Institute for Computational Biology, School of Science and Engineering, The Chinese University of Hong Kong, Shenzhen 518172, China
School of Biotechnology, University of Science and Technology of China, Hefei 230026, China
Biomedicine Discovery Institute, Monash University, Melbourne, 3800 VIC, Australia
State Key Laboratory of Functions and Applications of Medicinal Plants, College of Basic Medical, Guizhou Medical University, Guizhou 550004, China

⁎Corresponding author. yu100288@163.com (Fangfang Wang)

Disclaimer:
This article was originally published by Elsevier and was migrated to Scientific Scholar after the change of Publisher.

Abstract

Focal adhesion kinase (FAK) is a promising target for developing more effective anticancer drugs. To better understand the structure-activity relationships and mechanism of actions of FAK inhibitors, a molecular modeling study using 3D-QSAR, molecular docking, molecular dynamics simulations, and binding free energy analysis were conducted. Two types of satisfactory 3D-QSAR models were generated, comprising the CoMFA model (R2cv = 0.528, R2pred = 0.7557) and CoMSIA model (R2cv = 0.757, R2pred = 0.8362), for predicting the inhibitory activities of novel inhibitors. The derived contour maps indicate structural characteristics for substituents on the template. Molecular docking, molecular dynamic simulations and binding free energy calculations further reveal that the binding of inhibitors to FAK is mainly contributed from hydrophobic, electrostatic and hydrogen bonding interactions. In addition, some key residues (Arg14, Glu88, Cys90, Arg138, Asn139, Leu141, and Leu155) responsible for ligand-receptor binding are highlighted. All structural information obtained from 3D-QSAR models and molecular dynamics is consist with the available experimental activities. All the results will facilitate the optimization of this series of FAK inhibitors with higher inhibitory activities.

Keywords

FAK
CoMFA
CoMSIA
Molecular docking
Molecular dynamics
1

1 Introduction

Focal adhesion kinase (FAK) is a non-receptor protein-tyrosine kinase, which can control cellular responses (migration, morphology, cell proliferation, and cell survival) to a diverse group of stimuli, such as integrins, cytokines, chemokines, and growth factors (Frisch et al., 1996; Parsons, 2003; Rodriguez-Fernandez and Rozengurt, 1998). Many studies indicate that FAK is overexpressed in some tumors, which appears in head, neck, colon, breast, prostate, liver and thyroid, and so on (Kornberg, 2015; Sood et al., 2004; Judson et al., 1999; Cance et al., 2000). Transgenic mouse experiments have confirmed that specific FAK deletion would result in defect in tumor growth, invasion, and metastasis (Luo et al., 2009; Mclean et al., 2005; Lahlou et al., 2007). In addition, overexpression of FAK is linked with higher tumor stage and metastasis, which is also linked with poor prognosis in some tumors (breast, liver, lung, and esophagus) (Hochwald et al., 2009).

Unlike other non-receptor protein tyrosine kinases, FAK is composed of three domains: a central catalytic domain, N-terminal domain and C-terminal non-catalytic domain. The N-terminal FERM domain (Ceccarelli et al., 2006) mainly interacts with the receptor tyrosine kinase (epidermal growth factor receptor-EGFR, platelet-derived growth factor receptor-PDGFR), shown in Fig. 1 (Sieg et al., 2000). Protein-protein interaction sites are located in the C-terminal region, which includes two proline-rich areas (PR2-PR3) and a focal adhesion targeting (FAT) area. The C-terminal region mediates the binding to molecular regulators and effectors (Schaller, 2010). The central domain, known as an auto-inhibited structure, binds the N-terminal to FERM domain and it’s C-terminal to the proline-rich and FAT protein. Therefore, this domain is crucial for the FAK activity (Brami-Cherrier et al., 2014). However, the activity of FAK is inhibited when binding of the FERM domain to the central kinase domain (Ceccarelli et al., 2006), which can be activated by phosphorylation of some tyrosine residues: Tyr397 and Tyr407 (N-terminal domain), Tyr576 and Tyr577 (central domain), Tyr861 and Tyr925 (C-terminal domain). The auto-phosphorylation of the critical residue Tyr397 is triggered by extracellular stimuli (Mitra et al., 2005), serving as a canonical SH2-domain binding site for SRC-family kinases, for instance, PI3K and GRB7 (Chen et al., 1996; Han and Guan, 1999). Then, FAK would be full activated by phosphorylation of Tyr861 and Tyr925 through the above kinases, further leading to the activation of downstream AKT and ERK oncogenic pathways (Schlaepfer and Hunter, 1996; Xia et al., 2004). Therefore; FAK has been proposed as a potential target for cancer therapy.

The domain organization of FAK.
Fig. 1
The domain organization of FAK.

Numerous researchers have developed several FAK inhibitors, which would bind to the ATP-binding site of the FAK-kinase domain and block catalysis (Cellular Characterization, 2007; Roberts et al., 2008; PND, 2010), for example, PF-562271 (Roberts et al., 2008), PF-573228 (Mabeta, 2016), PF-431396 (Seungil et al., 2009), TAE226 (Liu et al., 2007), Y15 (Cance and Golubovskaya, 2008), Y11 (Golubovskaya et al., 2012), PND-1186 (Tancioni et al., 2014), Defactinib (Shimizu et al., 2016), GSK2256098 (Mak et al., 2019), BI-853520 (Lee and Gan, 2019), CT-707 (Wang et al., 2016), CEP-37440 (Ott et al., 2016), which are in different phases of clinical trials, and the structures of these compounds are listed in Table 1. All these compounds inhibit the proliferation of cancer mainly by means of inhibiting FAK phosphorylation or blocking signal transduction through the cell membrane.

Table 1 FAK inhibitors introduced in this work.
Compound name Structure
PF-562271
PF-573228
PF-431396
TAE226
Y15
PND-1186
Defactinib
GSK2256098
BI-853520 Structure not disclosed
CT-707
CEP-37440

The activities and selectivity of the novel discovered compounds are almost tested in human or calf origin, which procedure is time-consuming and costly. Considering this problem, molecular modeling approaches can be applied, which have proven effective in drug design and discovery (Pourbasheer et al., 2014; Wang et al., 2020; Sigalapalli et al., 2020; Righetti et al., 2020; Cichero et al., 2017; Guariento et al., 2017; Franchini et al., 2017). Therefore, in order to obtain more potent FAK inhibitors and gain an insight into the main molecular interactions involved FAK inhibition at the molecular level, a combination of three dimensional quantitative structure-activity relationship-3D-QSAR (include comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA)) (Hong et al., 2003; Klebe et al., 1994), molecular docking, molecular dynamics (MD), and binding free energies were performed on a series of pyrimidine analogs targeted on FAK. Based on the results obtained from this work, valuable information can be retrieved for further inhibitor design.

2

2 Materials and methods

2.1

2.1 Data sets

In the present work, a series of FAK inhibitors were taken from the works of Zhao et al. (Wang et al., 2020) and the related structures were shown in Table S1. The in vitro inhibitory activity (IC50) was initially converted to negative logarithmic units listed as pIC50 for CoMFA and CoMSIA analysis. Additionally, the whole data set was split into two data set: the training set and the test set, which was used to generate the final 3D-QSAR models and evaluate the predictive power of the derived models, respectively. Furthermore, the test set was selected considering the criteria that the structures of the test set inhibitors should be diverse and their inhibitory activities were similar to those of the training set (Zhang and Zhong, 2010).

The molecular modelling and 3D-QSAR studies were performed using the Sybyl software (Biotech Ltd, Shanghai, China). A systematic conformational search was initially performed to obtain stable structure of the inhibitors. Then, partial atom charges were added using the Gasteiger-Hückel method (Gasteiger and Marsili, 1980). Then, all inhibitors were minimized using a conjugate gradient minimizer, applying with the standard Tripos force field (Clark et al., 1989) with a distance-dependent dielectric function. The conformation can be achieved when the energy gradient convergence criterion was reached to 0.05 kcal/mol or 1000 step minimization cycle limit.

2.2

2.2 Alignment of the inhibitors

The alignment approach is one of the most sensitive part for CoMFA and CoMSIA analysis (Xia et al., 2002; Cho and Tropsha, 1995). In this work, three different alignment rules were adopted: the template ligand-based alignment (Alignment 1), the docking-base alignment (Alignment 2) and the common scaffold-based alignment (Alignment 3).

Alignment 1: due to the bioactive conformation of the inhibitors was unknown, the most potent inhibitor 24 was selected as the template. The common substructure in compound 24 for alignment is marked in blue (Fig. 2A). The remaining inhibitors were aligned to the template using the “Align Database” function. The result is shown in Fig. 2B.

(A) Compound 24 is used as FAK template for alignment. (B) The alignment result from template ligand-based alignment for FAK.
Fig. 2
(A) Compound 24 is used as FAK template for alignment. (B) The alignment result from template ligand-based alignment for FAK.

Alignment 2: it is assumed that homologous series of inhibitors would bind to the receptor in a similar way. Therefore, the conformations for each inhibitor derived from molecular docking were overlaid to do the docking-based alignment, and the result is shown in Fig. S1A.

Alignment 3: the process of scaffold-based alignment is the same as Alignment 1, while the conformations were originated from Alignment 2. This alignment strategy is listed in Fig. S1B.

2.3

2.3 CoMFA and CoMSIA

The aligned inhibitors were placed in a 3D grid box with a 2 Å spacing, such that the whole dataset compounds were included in it. The CoMFA fields (steric and electrostatic) were generated using a sp3 carbon atom with +1 charge, which were scaled by CoMFA-STD method. The Lennard-Jones potential and Coulomb potentials, representing steric and electrostatic fields, respectively, were calculated using the standard Tripos force fields (Lu et al., 2010; Buolamwini, 2009). To speed up the analysis and reduce noise, the energy cutoff and the minimum column filtering were set to 30 kcal/mol and 2.0 kcal/mol, respectively.

In the case of CoMSIA analysis, similarity indices were produced with the same lattice box that was used in CoMFA (Wold et al., 1984). In addition to steric and electrostatic fields, hydrophobic, hydrogen bond donor and hydrogen bond acceptor parameters were also calculated. The following standard parameters were applied for CoMSIA model: a radius of 1 Å, a charge of +1, hydrophobicity of +1, hydrogen bond donor of +1, hydrogen bond acceptor of +1, the default value of 0.3 as attenuation factor α, which the optimal value normally ranging from 0.2 to 0.4. Larger values of α result in a steeper Gaussian function, and an increasing attenuation of the distance-dependent effects of molecular similarity. Global molecular features become less important, and there is little averaging of local features. If the attenuation factor α is set to 0.3, each property value of a given atom is felt by 74.1% at 1 Å from the atom, by 30.1% at 2 Å, and by 6.7% at 3 Å.

2.4

2.4 Partial least squares (PLS) analysis

PLS method (Clark and Cramer, 2010; Bush and Nachbar, 1993; Cruciani et al., 1992) was employed to develop 3D-QSAR models in which the CoMFA and CoMSIA descriptors and inhibitory activities were used as independent and dependent variables, respectively. Before PLS analysis, the CoMFA and CoMSIA columns were filtered using column filtering. Then, the cross-validation analysis was performed using the leave-one-out (LOO) approach where one inhibitor was removed from the dataset and its activity was predicted using the derived model, and the cross-validated coefficient R2cv and the optimum number of components (ONC) were produced. Subsequently, non-cross-validation analysis (Baroni et al., 1992; Golbraikh and Tropsha, 2002) was employed to generate the final PLS regression model with non-cross-validated correlation coefficient (R2ncv) and standard error of estimate (SEE) parameters.

To assess the predictive power of the developed 3D-QSAR models, the predictive correlation coefficient (R2pred) (Golbraikh and Tropsha, 2000; Lietha and Eck, 2008) is calculated using Eq. (1):

(1)
R pred 2 = SD-PRESS /SD Where SD is the sum of squared deviation between the inhibitory activities of the test set and mean activities of the training set, and PRESS is the sum of squared deviation between the experimental and predicted activities for compounds in the test set.

2.5

2.5 Applicability domain (AD) analysis

The AD is a theoretical region in chemical space, which can be used to predict whether the developed 3D-QSAR models can be used to truly predict a novel compound, further validating the accuracy, reliability and external applicability of the model. Therefore, the AD for the developed 3D-QSAR models was calculated. In the present work, a simple approach for defining outliers (in the case of the training set) and the compounds residing outside the AD (for the test set) was applied, which can be achieved through the following link http://dtclab.webs.com/softwaretools or http://teqip.jdvu.ac.in/QSAR_Tools/.

2.6

2.6 Molecular docking

2.6.1

2.6.1 Receptor preparation

The crystal structure of FAK (PDB code: 2JKK) (Weiner et al., 1984) was retrieved from the protein data bank (http://www.rscb.org/pdb). Initially, the crystalized waters and ligand were removed from the original PDB file. Polar hydrogens, Kollman charges (Morris et al., 1998), atomic solvation parameters and fragmental volumes were assigned using Autodock Tools (ADT).

2.6.2

2.6.2 Ligand preparation

For all the FAK inhibitors, the 3D structures were constructed using the Sketch Molecule function in Sybyl software. Gasteiger partial charges and non-polar hydrogen atoms (Gasteiger and Marsili, 1980) were assigned. All torsions were set to rotatable during molecular docking.

2.6.3

2.6.3 Docking

In order to obtain the probable binding conformations, all inhibitors were docked into the binding site of FAK using AutoDock (The Scripps Research Institute, La Jolla, USA) (Mehler and Solmajer, 1991). Firstly, a grid box was constructed with dimensions of 60 × 60 × 60 Å separated by 0.375 Å, which was situated around the binding site based on the position of the co-crystallized ligand. The hydrogen bonds, van der Waals, electrostatic interactions were calculated using Lennard-Jones parameters 12–10, 12–6, and distance-dependent dielectric permittivity of Mehler and Solmajer, respectively (Tai-Sung et al., 2018). In addition, the number of molecular docking runs, the population in the genetic algorithm, the energy evaluations, and the maximum number of iterations were set to 50, 50, 250,000 and 27,000, respectively. Finally, the conformations with lowest binding energy were aligned together for constructing receptor-based 3D-QSAR models.

2.7

2.7 MD simulation

In order to validate the reliability of molecular docking and provide binding complexes with higher precision, thermo dynamic simulations were performed on FAK-Cpd 1 (the lowest inhibitor) and FAK-Cpd 24 (the potent inhibitor).

2.7.1

2.7.1 MD inputs

The docked complex Cpd 1-2JKK and Cpd 24-2JKK were employed as the starting conformations for MD simulations. For the ligands (Cpd 1 and Cpd 24), the parameters were generated using the antechamber of Ambertools 18 (Wang et al., 2007; Ashvar et al., 1996), the RESP charges calculated by the DFT function of B3LYP/6-31G** from Gaussian 09 were added (Kubicki et al., 2014; Wang et al., 2004), the bond constants were extracted from Amber GAFF (Simmerling et al., 2015). The FAK of the two systems were then defined by AMBER 14SB force field (Schow et al., 2011). The 16 Å TIP3P water layer (Case et al., 2018) was added (the dimensions of each water box were 82.09 × 81.69 × 83.80 Å3 and 82.09 × 81.69 × 83.80 Å3, respectively, for FAK-Cpd 1 and FAK-Cpd 24), and the systems were neutralized by adding counter ions (Cl− and Na+). Finally, each system was composed of approximately 45,824 atoms and 45,841 atoms for FAK-Cpd 1 and FAK-Cpd 24, respectively.

2.7.2

2.7.2 MD process

Classical MD simulations were preformed using the AMBER v18 (Particle Mesh Ewald Molecular Dynamics) (University of California, San Francisco, USA) (Ryckaert et al., 1977). The systems were initially minimized by 10,000 steps to remove bad contacts between water, ions and the complexes. The systems were then equilibrated: 50 ps of heating and 50 ps of density equilibration with weak restraints (5 kcal/mol), 500 ps of constant pressure (1 atm) equilibration at 300 K at a time step of 2 fs. The SHAKE algorithm was applied to constrain all bonds involving at least one hydrogen atom with non-bonded cut-off of 20.0 Å (Abraham et al., 2015). After this point, the simulation was lasted using the NVT ensemble and conformations were sampled every 10 ps. A simulation of 200 ns with 3 replicates was followed.

2.7.3

2.7.3 Local principal component analysis (PCA)

To search stable binding pose with minimum local free energy, the PCA was performed on the replicates of each production phase MD system by Gromacs v 5.11 (Kumar et al., 1992). In each system, the heavy atoms of the receptor and the inhibitor were combined and used for generating the displacement matrix of PCA. In the essential dynamic analysis (EDA), PC1&2 vectors calculated from PCA were used for distinguishing the top two significant movements of the system. The total combined production phase trajectories (600 ns for each system) were then projected onto a subspace which was defined by PC1 and PC2. Finally, the possibility of the microstates found in the subspace was calculated for the following local free energy F a landscape analysis, which was estimated by a weighted-histogram analysis (WHAM) (Marinelli et al., 2009; Wolfram Research, 2018). The free energy of a microstate a was described as follows:

(2)
F a = - T l o g i n a i j e 1 T ( f j - V a j ) Where n a i is the frequency of the microstate a that observed in the trajectory i . The normalization constants f j are self-consistently determined as in the WHAM method (Wolfram Research, 2018). The method described by Marinelli et al. was used to define the variation of the bias over different conformations that assigned to the same cluster a (Wolfram Research, 2018). The V a i is the bias potential acting on microstate a , which is estimated as the time average of the history-dependent potential acting on the trajectory i multiply by S a , the center of microstate a :
(3)
V a i = V G i ¯ S a = 1 t sim - t eq t eq t sim d t ' V G i ( S a , t ' )
Where t sim is the production simulation time, and t eq is the beginning of production phase, which is the last time frame of equilibrium when the bias potentials become stable. Therefore, according to Eqs. (2) and (3), we converted the possibility 2D map to the 2D local free energy landscape (LFEL) heat map and 3D LFEL surface map by using Gromacs v5.11 (Kumar et al., 1992)and Mathematica v11.3 (Kollman et al., 2000).

2.7.4

2.7.4 Binding free energy calculation

The binding free energies were further estimated by MM-GBSA and MM-PBSA method (Wang et al., 2001; Miller et al., 2012) with MMPBSA.py.MPI (Hawkins et al., 1996). The binding free energy Δ G binding is computed as follows:

(4)
Δ G binding = Δ G com - Δ G rec + Δ G lig
(5)
Δ G com / r e c / l i g = Δ H - T Δ S
(6)
Δ H = Δ E gas + Δ G sol
(7)
Δ E gas = Δ E int + Δ E vdw + Δ E ele
(8)
Δ G sol = Δ G PB / G B + Δ G NP
(9)
Δ G NP = γ S A S A + β
where Δ G com , Δ G rec and Δ G lig is the free energies of the complex, the receptor and the ligand, respectively. Δ G com / r e c / l i g is the discrepancy between enthalpy contribution ( Δ H ) and the conformational entropy ( T Δ S ). T is the temperature of the simulated environment. The Δ S is the entropy of the compound, estimated from a normal-mode analysis of harmonic frequencies calculated at the molecular mechanics (MM) level (Wang et al., 2001). The enthalpy ( Δ H ) is the sum of internal energy from gas phase ( Δ E gas ) and the solvation free energy ( Δ G sol ). Δ E gas is the standard gas phase energy, which includes internal energy ( Δ E int ), van der Waals interactions ( Δ E vdw ), and electrostatic energies ( Δ E ele ). Since complex MD simulations were only performed here, the Δ E int to the binding free energy was zero. The solvation energy Δ G sol is the sum of non-ploar energy ( Δ G NP ) and electrostatic energy ( Δ G PB / G B ). The Δ G PB is calculated from Poisson-Boltzmann function with the default cavity radii from AMBER palmtop files. The dielectric constant is set to 1 for the interior solute and 80 for the exterior solvent. While the Δ G GB is an alternative part for Δ G GB which uses Hawkins, Cramer, and Truhlar pairwise generalized Born model (Hou et al., 2006; Onufriev et al., 2000) with the parameters described by Tsui and Case (Weiser et al., 1999). The LCPO approach (Hermann, 1972) was used to calculate the solvent accessible surface area (SASA) for the estimation of the non-polar solvation. The γ (0.00542 kcal/mol*Å2) and β (0.92 kcal/mol) were taken from a linear regression of the solvation free energy of a set of small polar molecules in water (Sitkoff et al., 1994).

3

3 Results and discussion

3.1

3.1 CoMFA and CoMSIA studies

The results of the developed models (Table S2-S4) suggest that the template ligand-based alignment is superior to the other two alignments. For alignment 1, all possible field combinations were computed (Table S2), the optimum models were chosen for further analysis (Table 2).

Table 2 Statistical data of optimal QSAR models.
Parameters CoMFA CoMSIA
R2cv 0.528 0.757
R2ncv 0.987 0.971
SEE 0.047 0.067
F 193.472 116.840
R2pred 0.7557 0.8362
SEP 0.425 0.427
ONC 5 4
Field contribution
S 0.690 0.310
E 0.310
H 0.690
D
A

R2cv = Cross-validated correlation coefficient using the leave-one-out methods;

R2ncv = Non-cross-validated correlation coefficient; SEE = Standard error of the estimate; F = Ratio of R2ncv explained to unexplained = R2ncv/(1 − R2ncv);

R2pred = Predicted correlation coefficient for the test set of compounds; SEP = Standard error of prediction; ONC = Optimal number of principal components.

3.1.1

3.1.1 CoMFA analysis

For CoMFA model, PLS analysis shows an R2cv of 0.528. A non-cross-validated PLS analysis results in a conventional R2ncv of 0.987, F value of 193.472, and standard error of the estimate (SEE) of 0.047 with 5 components. The contributions of steric and electrostatic fields are 69% and 31%, respectively, indicating that the steric field has a greater influence on the inhibitory activities. The CoMFA model is further validated by external test set compounds with R2pred = 0.7557, and the fitness plot (Fig. 3A) using the experimental activity on the X-axis and predicted activity on the Y-axis proves that the CoMFA model is statistically significant and predictive.

Graphs of the predicted versus the experimental pIC50 values of the optimal models for FAK. (A) CoMFA model. (B) CoMSIA model.
Fig. 3
Graphs of the predicted versus the experimental pIC50 values of the optimal models for FAK. (A) CoMFA model. (B) CoMSIA model.

3.1.2

3.1.2 CoMSIA analysis

In the CoMSIA model, the best combination is steric and hydrophobic fields (R2cv = 0.757, R2ncv = 0.971, F = 116.84, SEE = 0.067, ONC = 4), indicating that this model is reliable and predictive. Furthermore, the steric and hydrophobic fields contribution are 31% and 69%, respectively, indicating that hydrophobic filed has important influence on the ligand-receptor interactions. At the same time, the decrease of the contribution of steric field in the CoMSIA model compared with CoMFA model can be explained by the incorporation of other more significant fields (hydrophobic) in the model. To validate the efficacy of the established CoMSIA model, the external validation using the test set compounds gives R2pred of 0.8362, further demonstrating that the CoMSIA model possesses satisfactory predictive ability. Fig. 3B illustrates the correlation plot of the experimental versus the predicted inhibitory activities. Clearly, all points are uniformly distributed around the regression line, suggesting the satisfactory predictive ability of the CoMSIA model. Besides, compared with the CoMFA model, the internal prediction and external prediction capabilities all superior to the former.

3.2

3.2 Contour maps analysis

The best CoMFA and CoMSIA models were selected to construct the 3D contour maps to view the effect of different fields on the target features. The contour maps around the inhibitors would identify the key features contributing to the binding activities. Furthermore, the contribution values of favored and disfavored regions were set to 80% and 20% by default. The most potent compound 24 was also selected as the template compound and displayed along with the contour maps.

3.2.1

3.2.1 CoMFA contour maps

The steric contour maps of the CoMFA model are shown in Fig. 4A, where green and yellow contours refer to steric favored and unfavored areas for the inhibitory activity, respectively. A large and a small yellow contours mapped near the 6-position of ring D (Fig. 2A), suggest that small groups are favored at this position, which is consistent with the fact that the activities of Cpds 1426 (–H) are higher than those of Cpds 113 (–CONHCH3), actually, it all comes down to the fact that the possible rotation around NH-ring D, positions 2 and 6 are equivalent, thus in 14–26, the hydrogen group is present on 6-position not a sulfonamide substituent. Therefore, the activities of Cpds 1426 are higher than Cpds 113. Several yellow polyhedrons located around ring C, indicate that minor substituents at this area might have favorable steric interactions with surrounding residues. For instance, Cpds 24 and 26 containing five-membered ring show higher inhibitory activities than Cpds 16, 17, 18, 19 and 21 which have six-membered ring at the same position. The steric favored green contour covering the 4-position of ring C reveals that bulky groups at this location would improve the activity. Take Cpds 15 and 14 as an example, the activity of Cpd 15 (pIC50 = 7.5436) is relatively higher than that of Cpd 14 (pIC50 = 7.2396), because Cpd 15 (3,4-dimethoxybenzene) has larger substituent compared with Cpd 14 (benzene). Additionally, a yellow contour map is observed at the methylene group of ring D, suggesting that minor substituent is favorable for the inhibitory activity. This is consistent with the activity order of Cpd 11 (pIC50 = 7.2612) and Cpd 24 (pIC50 = 8.1487) (Cpd 24 > Cpd 11). It is mainly due to the lack of methylene group in Cpd 11 at this position, which leads to the –NH directly linking with benzene ring, at the same time, the benzene ring with larger volume extends to the yellow contour map, thus reducing its activity.

StDev*Coeff contour plots for FAK inhibitors in combination of Cpd 24. (A) The steric contour map, where the green and yellow contours represent 80% and 20% level contributions, respectively. (B) The electrostatic contour map, where the blue and red contours represent 80% and 20% level contributions, respectively.
Fig. 4
StDev*Coeff contour plots for FAK inhibitors in combination of Cpd 24. (A) The steric contour map, where the green and yellow contours represent 80% and 20% level contributions, respectively. (B) The electrostatic contour map, where the blue and red contours represent 80% and 20% level contributions, respectively.

For electrostatic contour maps shown in Fig. 4B, the red and blue contours represent the negative charge and positive charge favorable regions, respectively. There is a small red region near ring D, indicating that the introduction of negatively charged groups near this position is conducive to increase the inhibitory activity, therefore, modifications can be made at this position to improve the activity. A big and a small red contours near the 2-position of ring D indicate the need of electronegative groups for improving potency. Cpds 1416 with electronegative groups at that position are more active than Cpds 113 with hydrogen substituent. Red and blue contour maps are distributed near ring C, indicating that the substituents here need to be carefully selected. Additionally, another red contour near the 4-position of ring C suggests the need of electronegative atoms at this area can increase the activity. For example, the order of activity of Cpds 1416 (pIC50 = 7.2396, pIC50 = 7.5436, pIC50 = 7.6421, respectively) is Cpd 16 (3,4,5-trimethoxybenzene) > Cpd 15 (3,4-dimethoxybenzene) > Cpd 14 (benzene).

3.2.2

3.2.2 CoMSIA contour maps

The CoMSIA steric contour maps (Fig. 5A) are similar to those obtained from the CoMFA model. However, other small yellow contour map is observed at the sulfonamide group of ring D, indicating that minor groups are beneficial to the activity.

CoMSIA StDev*Coeff contour plots for FAK inhibitors in combination of Cpd 24. (A) The steric contour map, where the green and yellow contours represent 80% and 20% level contributions, respectively. (B) The hydrophobic contour map, where the yellow and white contours represent 80% and 20% level contributions, respectively.
Fig. 5
CoMSIA StDev*Coeff contour plots for FAK inhibitors in combination of Cpd 24. (A) The steric contour map, where the green and yellow contours represent 80% and 20% level contributions, respectively. (B) The hydrophobic contour map, where the yellow and white contours represent 80% and 20% level contributions, respectively.

Fig. 5B shows the CoMSIA hydrophobic contour maps, where yellow contours represent hydrophobically favored regions and white contours indicate hydrophilically favored regions. There is a yellow contour around ring D, indicating that hydrophobic substituents in this region are favored. This may be the reason why Cpds 1426 with relatively hydrophobic groups at this position show higher potencies than Cpds 113. Furthermore, several yellow contours are observed near ring C, illustrating that hydrophobic groups would benefit the inhibitory activity. For instance, Cpd 2 (pIC50 = 7.2233) shows improved potency than Cpd 1 (pIC50 = 6.7597). Moreover, the activity discrepancies for Cpd 22 (pIC50 = 7.8153) and Cpd 24 (pIC50 = 8.1487) can also be explained by these yellow contour maps. Another yellow contour map is located around the methylene group of ring D, which indicate that the hydrophobic groups are important for the inhibitory activity. For example, Cpd 24 exhibits higher activity than Cpd 11. This phenomenon is originated from the fact that the methylene is present in this position of Cpd 24, which increases the flexibility of the bond rotation and makes the benzene ring extend to the yellow contour map, while Cpd 11 in the position has –NH directly connecting to the benzene ring, which reduces the rotation flexibility, thus the –CONH group extends to the yellow contour and reduces the activity. In addition, a white contour map is located at the oxygen atom at 4-position of ring C, indicating that this region favors hydrophilic groups. It can be validated by Cpd 24 (pIC50 = 8.1487) and Cpd 23 (pIC50 = 7.4260). Here, Cpd 24 has oxygen atom while Cpd 23 has fluorine atom at this position, therefore, the activity of Cpd 24 is higher than Cpd 23.

3.3

3.3 AD analysis

In case of the CoMFA model, it is interesting to point out that there is no outlier for the training set and the test set, indicating that the reliability of the derived model from another aspect. For the CoMSIA model, no outliers are detected for the training set and test set. Overall, the results of AD further indicate that the developed models are reliable and applicable for predicting novel compounds with common structure.

3.4

3.4 Molecular docking

The accuracy of molecular docking was validated before docking. The crystallized ligand B19 was extracted and re-docked into the binding site of FAK protein. Fig. 6 displays the superimposed conformation between the co-crystallized and the docked conformation with Root Mean Square Deviation (RMSD) value of 0.365 Å, suggesting accurate pose prediction. Therefore, all FAK inhibitors were docked into the binding site using the obtained parameters, and the docked poses with lowest binding energies were used as starting conformations for following MD simulations.

Re-docking pose with the RMSD value of 0.365 Å (Cyan = Original, Red = Docked).
Fig. 6
Re-docking pose with the RMSD value of 0.365 Å (Cyan = Original, Red = Docked).

3.5

3.5 MD simulations

In order to examine the reliability of molecular docking and further investigate the significant residues involved in inhibitor-receptor interactions, MD simulations were performed on the docked complexes of Cpd 24 and Cpd 1.

3.5.1

3.5.1 Plasticity of the MD systems

The resilience of FAK and ligands were initially investigated by calculating the RMSD of the backbone atoms (Cα, N, O, C) and the heavy atoms of the docked ligands. The residue flexibility was further examined by calculating the Root Mean Square Fluctuation (RMSF) of the backbone atoms, which can be summarized to exhibit the contribution on the residue-level. The trajectories were aligned by the Cα atoms of the starting structures of MD simulations beforehand to eliminate the rotation and transition of the whole systems.

The RMSD analysis shows that FAK fluctuates a little bit more (2.43 ± 0.52 Å) in the Cpd 24 system than in the Cpd 1 system (2.32 ± 0.44 Å), whereas Cpd 24 is more stable (0.82 ± 0.25 Å) than Cpd 1 (2.42 ± 0.35 Å) (Fig. 7). In run 2 of the FAK-Cpd 1 system, the fluctuation of the ligand is even reached 3.4 ± 0.35 Å, indicating that the binding of Cpd 24 to FAK is more stable than Cpd 1.

RMSD of protein backbone and heavy atoms of Cpd 1 and Cpd 24. (A). The RMSD of the three replicates of FAK-Cpd 1. (B). The RMSD of the three replicates of FAK-Cpd 24.
Fig. 7
RMSD of protein backbone and heavy atoms of Cpd 1 and Cpd 24. (A). The RMSD of the three replicates of FAK-Cpd 1. (B). The RMSD of the three replicates of FAK-Cpd 24.

Additionally, RMSF analysis (Fig. 8) suggests that the major fluctuations are mainly located at the flexible loop residues 570–584 (8.87 ± 2.25 Å in FAK-Cpd 1, 6.24 ± 1.25 Å in FAK-Cpd 24) for Cpd 1 and Cpd 24 systems. This is indeed reasonable as the loop misses 17 residues in the original crystal structure, and we have modelled it in the protein. Additionally, it can be seen that those residues (Pro444-Pro447) situated at the rim of the protein fluctuate a lot more (2.54 ± 0.32 Å) in both systems. In addition, the ligand binding domain is more stable in Cpd 24 system (0.89 ± 0.14 Å) than that in the Cpd 1 system (1.23 ± 0.14 Å), suggesting that the fluctuation of Cpd 1 is higher than Cpd 24 system. Therefore, the binding of Cpd 1 is not as stable as that of Cpd 24.

The RMSF plot of the two systems. Replicate trajectories were combined.
Fig. 8
The RMSF plot of the two systems. Replicate trajectories were combined.

3.5.2

3.5.2 Essential dynamics(ED) and free energy landscape (FEL)

For studying the overall dynamics that contributes to the ligand binding domain (LBD), we performed Principal Components Analysis (PCA) on the main-chain atoms and the heavy atoms of the ligands of the two MD systems, respectively. PC1 together with PC2 in the two PCA took more than 50% of the overall covariance of the atom displacements. Thus, the major overall movements were obtained by analyzing the essential dynamics (ED) that projected only along PC1 and 2 per MD system. Then, the probability based kinetic free energy landscape (FEL) was calculated to identify the most stable binding states using the PC1 and 2 defined subspace.

In the FAK-Cpd 1 MD system, it can be seen from PC1 that domain 1 experiences a swing movement by 5°, and domain 2 is relatively stable, except the modeled missing loop (possessing huge swing movement) (Fig. 9A), which makes the movements in the LBD quite subtle. In the PC2 of FAK-Cpd 1 system, we can only observe little rotation movement of the modeled missing loop with the other parts behaving rigidly (Fig. 9B). In the PC1 and PC2 defined FEL, with the docked and minimized structure of 4.1 kT, three metastable states are discovered (Fig. 9C): −5.5, 4, 2.1 kT, −5.1, 4.1, 1.7 kT and 10.4, 0.5, 0.2 kT (the most stable state).

The ED, FEL of FAK-Cpd 1. (A) The ED was obtained along PC1. (B) The ED obtained along PC2. Red (from PC1) and green arrows (from PC2) were proportional to the movement scale of the Cα atoms in FAK (C) The FEL on the same PC1&PC2 defined subspace of FAK-Cpd 1.
Fig. 9
The ED, FEL of FAK-Cpd 1. (A) The ED was obtained along PC1. (B) The ED obtained along PC2. Red (from PC1) and green arrows (from PC2) were proportional to the movement scale of the Cα atoms in FAK (C) The FEL on the same PC1&PC2 defined subspace of FAK-Cpd 1.

In the FAK-Cpd 24 MD system, domain 1 only experiences little swing (<1 Å) movement along PC1 ED movements (Fig. 10A). Only a small scale (around 0.3 Å) of rotating movements are observed in the LBD. Large swing movements are found in the modeled missing loop as observed in PC1 of FAK-Cpd 1 system. The movements along PC2 is much smaller than that in the PC1, except for the mild rotation of the modeled missing loop. In addition, we literally observe no motions in the LBD for Cpd 24. The converted FEL from PC1 and PC2 defined subspace indicates that most stable meta-stable states locate around (5.5, 2.9), with a reference free energy of 0.1 kT.

The ED, FEL of FAK-Cpd 24. (A) The ED obtained along PC1. (B) The ED obtained along PC2. Color representations were kept same with Fig. 9. (C) The FEL on the same PC1&PC2 defined subspace of FAK-Cpd 24.
Fig. 10
The ED, FEL of FAK-Cpd 24. (A) The ED obtained along PC1. (B) The ED obtained along PC2. Color representations were kept same with Fig. 9. (C) The FEL on the same PC1&PC2 defined subspace of FAK-Cpd 24.

Overall, the structures derived from the meta-stable states were taken to the binding free energy calculations.

3.5.3

3.5.3 Stable conformation and activity discrepancy analysis

The conformation of the most active Cpd 24 retrieved from the meta-stable states is employed for analyzing the detailed interactions between ligand and receptor, which can also be conducted to facilitate validation of the developed 3D-QSAR models and derive key residues in the binding site of FAK. The interactions between Cpd 24 and FAK are shown in Fig. 11. Cpd 24 is placed into the FAK binding pocket lined by Arg14, Ile16, Gly17, Glu18, Gly19, Val24, Met87, Glu88, Leu89, Cys90, Thr91, Leu92, Gly93, Glu94, Arg138, Asn139, Leu141, Gly151, Asp152, Leu155, and Ser156 (Fig. 11A). As shown in Fig. 11B, Cpd 24 displays three hydrogen bonds with Glu88 and Cys90. One of the hydrogen bond is observed between the –NH of ring A and residue Glu88 (–O⋯HN, 1.82 Å, 160.5°) (H-1). Another two hydrogen bonds are formed between the nitrogen atom of ring B, the –NH group between ring B and ring D and Cys90 (–N⋯HN, 2.39 Å, 166.3°) (H-2), (–O⋯HN, 2.06 Å, 142.1°) (H-3).

(A) The residues in FAK active site around Cpd 24. (B) The enlargement for Cpd 24 in the binding site after molecular docking, which is displayed in stick, hydrogen bonds are shown as dotted black lines, and the nonpolar hydrogens were removed for clarity.
Fig. 11
(A) The residues in FAK active site around Cpd 24. (B) The enlargement for Cpd 24 in the binding site after molecular docking, which is displayed in stick, hydrogen bonds are shown as dotted black lines, and the nonpolar hydrogens were removed for clarity.

In addition, from Fig. 11A, we can clearly see that ring C locates on proximity of residues Cys90, Thr91, Leu92 and Gly93, suggesting there is a steric hindrance for the receptor, which is collaborated with the CoMFA yellow contour map at this area (Fig. 4A and Fig. 5A). The groups at 4-position of ring C are occupied by green contour map (Fig. 4A and Fig. 5A), which is consistent with the interaction environment where a big cavum is appeared.

It can be seen from Fig. 11A that the substituent at 2-position of ring D is closer to electropositive amino acids Arg138 and Asn139, indicating that negatively charged groups are favorable for the ligand-receptor interactions. This is in agreement with the red contours shown in Fig. 4B. MD simulations also display that the only amino acid closer to the substituent at 4-position of ring C is positively charged Arg14, indicating that groups at this position should be electronegative to fit into the active site, proved by the red contour map listed in Fig. 4B.

As shown in Fig. 11A, there is a hydrophobic pocket formed by Leu141 and Leu155, illustrating that hydrophobic groups at ring D are favored for the inhibitory activity, evidenced by the presence of a yellow contour (hydrophobic favorable) around this area by the CoMSIA model (Fig. 5B). Additionally, the oxygen atom at 4-position of ring C is adjacent to hydrophilic amino acid Arg14, therefore, hydrophilic substituent at this position is needed, confirmed by the white contour (Fig. 5B).

The above analysis of the ligand-receptor interactions further confirms that the established 3D-QSAR models are accurate and effective.

Furthermore, the interactions of Cpd 1 and FAK showing various level of inhibitory activity were selected to explore the binding mode and activity discrepancy. The lowest inhibitor Cpd 1 presents interactions with surrounding residues (Fig. 12A): Ile16, Glu18, Gly19, Val24, Glu94, Arg138, Asn139, Val140, Leu141, Gly151, Asp152, Leu155, Arg157, Tyr158, Met159, Glu160, Ser162, Thr163, Tyr164 and Tyr165. As can be seen from Fig. 12B, the –NH group between ring B and ring D (Fig. 12C) forms hydrogen bond with residue Leu155 (–O⋯HN, 1.79 Å, 163.4°) (H-1), and the substituent of ring D also forms two hydrogen bonds with Ser156 (–N⋯HO, 2.94 Å, 111.3°) (H-2), (–O⋯HN, 1.82 Å, 161.0°) (H-3).

(A) The residues in FAK active site around Cpd 1. (B) The enlargement for Cpd 1 in the binding site after molecular docking, which is displayed in stick, hydrogen bonds are shown as dotted black lines, and the nonpolar hydrogens were removed for clarity. (C) Cpd 1 is used for analysis.
Fig. 12
(A) The residues in FAK active site around Cpd 1. (B) The enlargement for Cpd 1 in the binding site after molecular docking, which is displayed in stick, hydrogen bonds are shown as dotted black lines, and the nonpolar hydrogens were removed for clarity. (C) Cpd 1 is used for analysis.

The predicted binding modes of Cpd 1 and Cpd 24 derived from molecular docking and MD simulations are shown in Fig. 13. As shown in Fig. 13A and 13B, we can see that Cpd 1 and Cpd 24 could interact with FAK well both in the molecular docking and MD simulation. However, several differences still exist between the conformations obtained from the two procedures for Cpd 1 and Cpd 24 (Fig. 13C and Fig. 13D). In the case of FAK-Cpd 24, the main difference between docking and MD is the orientation of the substituents at ring C and ring D. The main difference in FAK-Cpd 1 is that the whole Cpd 1 is moved outside of the pocket, indicating that the binding mode of these inhibitors originated from molecular docking is not very stable. Additionally, Cpd 1 and Cpd 24 are aligned together for molecular docking and MD simulation, respectively. Due to the different substituents, Cpd 1 and Cpd 24 exhibits a huge difference in inhibitory activity. This phenomenon occurs because of the altered conformations of these compounds in the binding pocket, and the overall structures are changed, which makes the distribution of the surrounding amino acids altered accordingly. It finally leads to the changes in the surrounding electrostatic, hydrophobic, and steric fields.

(A) Superimposition of FAK-Cpd 24 and FAK-Cpd 1 derived from molecular docking. (B) Superimposition of FAK-Cpd 24 and FAK-Cpd 1 derived from MD simulation. (C) Superimposition of Cpd 24 derived from molecular docking and MD simulation. (D) Superimposition of Cpd 1 derived from molecular docking and MD simulation.
Fig. 13
(A) Superimposition of FAK-Cpd 24 and FAK-Cpd 1 derived from molecular docking. (B) Superimposition of FAK-Cpd 24 and FAK-Cpd 1 derived from MD simulation. (C) Superimposition of Cpd 24 derived from molecular docking and MD simulation. (D) Superimposition of Cpd 1 derived from molecular docking and MD simulation.

As shown in Fig. 14A and Fig. 14C (the contour maps for Cpd 1), some yellow contours are covering ring C, indicating that minor substituents are favorable for higher inhibitory potency. Is is also the reason why the activity of Cpd 1 is lower than that of other compounds containing five-membered rings in the same series. Furthermore, the region of green contour map distributed around the left of ring C indicates a sterically bulky substituent extended to this field is favorable. It also shows that the presence of larger groups on ring C extending to the green contour map is beneficial to the higher activity of the compound. For example, the activities of Cpds 26 are higher than Cpd 1. A yellow region is also observed around the amido group of ring D, suggesting that bulky substituent in this region is unfavorable for the inhibitory activity. This could be part of the reason why Cpd 24 with non-additionally substituent here has higher activity than Cpd 1. Meanwhile, the sulfonamide group of Cpds 1426 is located at the pocket formed by Glu94, Arg136, Arg138, Glu18, and Gly19, which facilitates the formation of hydrophilic interactions, while the amido group of Cpds 113 is close to form hydrogen bonds with amino acids Asp152 and Ser156. This situation suggests that the two groups extend in different directions in the receptor binding pocket, leading to the activity discrepancy.

CoMFA StDev*Coeff contour plots for FAK inhibitors in combination of Cpd 1. (A) The CoMFA steric contour map, where the green and yellow contours represent 80% and 20% level contributions, respectively. (B) The CoMFA electrostatic contour map, where the blue and red contours represent 80% and 20% level contributions, respectively. (C) The CoMSIA steric contour map, where the green and yellow contours represent 80% and 20% level contributions, respectively. (D) The hydrophobic contour map, where the yellow and white contours represent 80% and 20% level contributions, respectively.
Fig. 14
CoMFA StDev*Coeff contour plots for FAK inhibitors in combination of Cpd 1. (A) The CoMFA steric contour map, where the green and yellow contours represent 80% and 20% level contributions, respectively. (B) The CoMFA electrostatic contour map, where the blue and red contours represent 80% and 20% level contributions, respectively. (C) The CoMSIA steric contour map, where the green and yellow contours represent 80% and 20% level contributions, respectively. (D) The hydrophobic contour map, where the yellow and white contours represent 80% and 20% level contributions, respectively.

The interpretation of the electrostatic contour is listed in Fig. 14B. A red and a blue contours are found at ring C, suggesting that modifications at this position should be made for Cpd 1 for higher activity. A red contour map is located at the left of ring C, which indicates that electronegative groups in this region are favorable for potent activity. This explains why Cpd 24 is more active than Cpd 1.

In the hydrophobic field (Fig. 14D), there are several yellow contours near the substituent of ring D, indicating that hydrophobic groups here would increase the activity, however, Cpd 1 possesses hydrophilic group at this point, leading to the reduced inhibitory activity.

3.5.4

3.5.4 MM/PB&GBSA analysis

The above results are in agreement with the experimental data, but more details are still needed for further analysis. The binding free energy between FAK receptor and the two ligands was elaborately calculated using the snapshots extracted from the most stable binding states discovered by FEL.

The results of binding free energy calculation using MM-GBSA and MM-PBSA are shown in Table S5. The binding free energy for Cpd 1 and Cpd 24 is −26.29 ± 2.34 kcal/mol and −13.33 ± 1.90 kcal/mol, respectively (GBSA), and −14.04 ± 2.42 kcal/mol, −4.88 ± 1.89 kcal/mol for PBSA (Cpd 1 and Cpd 24), indicating that Cpd 24 is a quite favorable inhibitor for FAK compared with Cpd 1. The van der Waals energy is −57.56 kcal/mol and −40.67 kcal/mol for Cpd 24 and Cpd 1, respectively, which are much bigger than other energies, suggesting that hydrophobic feature is the main contribution for ligand-receptor interactions. The electrostatic energy is −27.59 kcal/mol and −16.58 kcal/mol for Cpd 1 and Cpd 24, respectively, further illustrating that electrostatic interaction is also important during the binding process. It is noted that the polar contributions are unfavorable to the binding free energy, which is due to the large volume of the ligand binding pocket, leading to the ligands are exposed in the solvent.

Interestingly, the ligand-receptor binding energies of the studied compounds (Cpd 1 and Cpd 24) correlates well with the predicted values calculated by the developed 3D-QSAR models, which in turn, further validates the reliability of the constructed models.

3.5.5

3.5.5 Binding free energy decomposition

We also decomposed the binding free energy (PBSA) after MD simulations to locate the key residues contribute to ligand binding (Fig. 15). Comparing Cpd 1 with Cpd 24, the key residues contributing to the free energy terms are Leu141, Leu155, Cys90, Leu89, Glu88, Ile16. Interestingly, in both systems, residue Asp152 has the negative contribution to the ligand binding, which possesses 3.82 kcal/mol unfavorable energy contribution mostly to Cpd 24, suggesting that Cpd 24 can still be improved at ring D for further inhibitor design (Fig. 2A).

MM-PBSA binding free energy decomposition with key residues contribution.
Fig. 15
MM-PBSA binding free energy decomposition with key residues contribution.

3.6

3.6 Design for novel FAK inhibitors

According to the above analysis, the main characteristics affecting the inhibitory activity were derived in the present work and the substituents were modified at related parts on the basis of the contour maps. The novel designed inhibitors were aligned in the same way as the dataset compounds and were fitted to the 3D-QSAR models. The structures and predicted activities are shown in Table 3.

Table 3 The compounds of newly designed and the predicted values.
Cpd 4-position of ring C ring C 2-position of ring D ring D Predicted activity (pIC50)
CoMFA CoMSIA
1 8.1496 8.1501
2 8.2435 8.3210
3 8.2563 8.2745
4 8.3011 8.3102
5 8.3021 8.2987

4

4 Conclusions

In the present work, a series of computational studies (3D-QSAR, molecular docking and MD simulations, binding free energy analysis) was carried out on FAK inhibitors. CoMFA and CoMSIA approaches were employed to derive 3D-QSAR models, which possess strong reliability and predictive capability. The produced contour maps delineate the structure-activity relationship of the studied inhibitors, which are capable of predicting the activity of the external compounds in a reliable manner. Then, molecular docking was applied to provide initial conformation for following MD simulations. Afterwards, MD simulations reveal that the most favorable contributions to ligand-receptor interactions are mainly come from Arg14, Glu88, Cys90, Arg138, Asn139, Leu141, and Leu155. Comparison between Cpd 1 and Cpd 24 confirms that it is significant to have suitable substituents in ring C and ring D for developing more active FAK inhibitors. Meanwhile, the binding free energy calculations indicate that the hydrophobic, electrostatic and hydrogen bonding interactions are significant for maintaining the activities of these inhibitors. The 3D contour maps correlate well with the detailed binding structures derived from MD simulations, and the detailed structure-activity relationship is concluded as follows: (1) minor and hydrophobic groups at ring C are favorable for the inhibitors of FAK; (2) bulky, electronegative, and hydrophilic groups at the 4-position of ring C are favorable; (3) electronegative and hydrophobic substituents at ring D are beneficial for the improved activity; (4) electronegative groups located at 2-position of ring D will enhance the activity; (5) small groups at 6-position of ring D are favored. Overall, the developed models and elucidation of ligand-receptor interactions will be helpful for rational design of novel and potent FAK inhibitors.

Acknowledgements

The study was supported by the National Natural Science Foundation of China (No. 32001699) and the Shandong Provincial Natural Science Foundation of China (No. ZR2019BC102).

Declaration of Competing Interest

The authors declare that there are no conflicts of interest.

References

  1. , , , , , , , . GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX. 2015;1:19-25.
    [Google Scholar]
  2. , , , , , . Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra: 6,8-Dioxabicyclo[3.2.1]octane. J. Phys. Chem.. 1996;100:9262-9270.
    [Google Scholar]
  3. , , , , , , . Predictive ability of regression models. Part II: Selection of the best predictive PLS model. J. Chemometrics. 1992;6:347-356.
    [Google Scholar]
  4. , , , , , , , , , , . FAK dimerization controls its kinase-dependent functions at focal adhesions. EMBO J.. 2014;33
    [Google Scholar]
  5. , . CoMFA and CoMSIA 3D-QSAR studies on S6-(4-nitrobenzyl)mercaptopurine riboside (NBMPR) analogs as inhibitors of human equilibrative nucleoside transporter 1 (hENT1) Bioorg. Med. Chem. Lett. 2009
    [Google Scholar]
  6. , , . Sample-distance partial least squares: PLS optimized for many variables, with application to CoMFA. J. Comput. Aided Mol. Des.. 1993;7:587-619.
    [Google Scholar]
  7. , , . Focal Adhesion Kinase Versus p53: Apoptosis or Survival? Sci. Signaling. 2008;1:pe22.
    [Google Scholar]
  8. , , , , , , , , . Immunohistochemical Analyses of Focal Adhesion Kinase Expression in Benign and Malignant Human Breast and Colon Tissues: Correlation with Preinvasive and Invasive Phenotypes. Clin. Cancer Res. Off. J. Am. Assoc. Cancer Res.. 2000;6:2417-2423.
    [Google Scholar]
  9. I.Y.B.-S. D.A. Case, S.R. Brozell, D.S. Cerutti, T.E. Cheatham, III, V.W.D. Cruzeiro, T.A. Darden, R.E. Duke, D. Ghoreishi, M.K. Gilson, H. Gohlke, A.W. Goetz, D. Greene, R Harris, N. Homeyer, S. Izadi, A. Kovalenko, T. Kurtzman, T.S. Lee, S. LeGrand, P. Li, C. Lin, J. Liu, T. Luchko, R. Luo, D.J. Mermelstein, K.M. Merz, Y. Miao, G. Monard, C. Nguyen, H. Nguyen, I. Omelyan, A. Onufriev, F. Pan, R. Qi, D.R. Roe, A. Roitberg, C. Sagui, S. Schott-Verdugo, J. Shen, C.L. Simmerling, J. Smith, R. Salomon-Ferrer, J. Swails, R.C. Walker, J. Wang, H. Wei, R.M. Wolf, X. Wu, L. Xiao, D.M. York and P.A. Kollman, AMBER 18, University of California, San Francisco., 2018.
  10. , , , , , . Crystal Structure of the FERM Domain of Focal Adhesion Kinase. J. Biol. Chem.. 2006;281:252-259.
    [Google Scholar]
  11. Cellular Characterization of a Novel Focal Adhesion Kinase lnhibitor, J. Biol. Chem. (2007).
  12. , , , , , , , , . Phosphorylation of tyrosine 397 in focal adhesion kinase is required for binding phosphatidylinositol 3-kinase. J. Biol. Chem.. 1996;271(42):26329-26334.
    [Google Scholar]
  13. , , . Cross-validated R2-guided region selection for comparative molecular field analysis: a simple method to achieve consistent results. J. Med. Chem.. 1995;38:1060-1066.
    [Google Scholar]
  14. , , , , , , , , . Benzimidazole-based derivatives as privileged scaffold developed for the treatment of the RSV infection: a computational study exploring the potency and cytotoxicity profiles. J. Enzyme Inhib. Med. Chem.. 2017;32:375-402.
    [Google Scholar]
  15. , , , . Validation of the general purpose Tripos 5.2 force field. J. Comput. Chem.. 1989;10:982-1012.
    [Google Scholar]
  16. , , . The Probability of Chance Correlation Using Partial Least Squares (PLS) QSAR Comb. Sci.. 2010;12:137-145.
    [Google Scholar]
  17. , , , , , , . Predictive ability of regression models. Part I: Standard deviation of prediction errors (SDEP) J. Chemom.. 1992;6:335-346.
    [Google Scholar]
  18. , , , , , , , , , , . Structure-Activity Relationships within a Series of σ1 and σ2 Receptor Ligands: Identification of a σ2 Receptor Agonist (BS148) with Selective Toxicity against Metastatic Melanoma. ChemMedChem. 2017;12:1893-1905.
    [Google Scholar]
  19. , , , , . Control of adhesion-dependent cell survival by focal adhesion kinase. J. Cell Biol.. 1996;134:793-799.
    [Google Scholar]
  20. , , . Iterative partial equalization of orbital electronegativity—a rapid access to atomic charges. Tetrahedron. 1980;36:3219-3228.
    [Google Scholar]
  21. , , . Predictive QSAR modeling based on diversity sampling of experimental datasets for the training and test set selection. Mol. Diversity. 2000;5:231-243.
    [Google Scholar]
  22. , , . Beware of q2! J. Mol. Graph. Model.. 2002;20:269-276.
    [Google Scholar]
  23. , , , , , , , , , , . A small molecule focal adhesion kinase (FAK) inhibitor, targeting Y397 site: 1-(2-hydroxyethyl)-3, 5, 7-triaza-1-azoniatricyclo [3.3.1.1(3,7)]decane; bromide effectively inhibits FAK autophosphorylation activity and decreases cancer cell viability, clonoge. Carcinogenesis. 2012;33:1004-1013.
    [Google Scholar]
  24. , , , , , , , . Exhaustive CoMFA and CoMSIA analyses around different chemical entities: a ligand-based study exploring the affinity and selectivity profiles of 5-HT1A ligands. J. Enzyme Inhib. Med. Chem.. 2017;32:214-230.
    [Google Scholar]
  25. , , . Association of Focal Adhesion Kinase with Grb7 and Its Role in Cell Migration. J. Biol. Chem.. 1999;274:24425-24430.
    [Google Scholar]
  26. , , , . Parametrized models of aqueous free energies of solvation based on pairwise descreening of solute atomic charges from a dielectric medium. J. Phys. Chem. – US. 1996;100:19824-19839.
    [Google Scholar]
  27. , . Theory of hydrophobic bonding. II. Correlation of hydrocarbon solubility in water with solvent cavity surface area. J. Phys. Chem.. 1972;76:2754-2759.
    [Google Scholar]
  28. , , , , , , , , , , . A novel small molecule inhibitor of FAK decreases growth of human pancreatic cancer. Cell Cycle. 2009;8:2435-2443.
    [Google Scholar]
  29. , , , , , , . Comparative molecular field analysis (CoMFA) model using a large diverse set of natural, synthetic and environmental chemicals for binding to the androgen receptor. SAR QSAR Environ. Res.. 2003;14:373-388.
    [Google Scholar]
  30. , , , , , . Computational analysis and prediction of the binding motif and protein interacting partners of the Abl SH3 domain. PLoS Comput. Biol.. 2006;2:e1.
    [Google Scholar]
  31. , , , , . Overexpression of focal adhesion kinase, a protein tyrosine kinase, in ovarian carcinoma. Cancer. 1999;86
    [Google Scholar]
  32. , , , . Molecular similarity indices in a comparative analysis (CoMSIA) of drug molecules to correlate and predict their biological activity. J. Med. Chem.. 1994;37:4130-4146.
    [Google Scholar]
  33. , , , , , , , , , , , , , , , . Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res.. 2000;33:889-897.
    [Google Scholar]
  34. , . Focal adhesion kinase expression in oral cancers. Head Neck. 2015;20:634-639.
    [Google Scholar]
  35. , , , , . Quantum mechanical calculations on cellulose–water interactions: structures, energetics, vibrational frequencies and NMR chemical shifts for surfaces of Iα and Iβ cellulose. Cellulose. 2014;21:909-926.
    [Google Scholar]
  36. , , , , , . The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method. J. Comput. Chem.. 1992;13:1011-1021.
    [Google Scholar]
  37. , , , , , , , . Mammary epithelial-specific disruption of the focal adhesion kinase blocks mammary tumor progression. Proc. Natl. Acad. Sci. USA. 2007;104:20302-20307.
    [Google Scholar]
  38. , , . BI 853520, a FAK-Simile of Prior FAK Inhibitors? Targeted Oncol.. 2019;14:39-41.
    [Google Scholar]
  39. , , . Crystal structures of the FAK kinase in complex with TAE226 and related bis-anilino pyrimidine inhibitors reveal a helical DFG conformation. PLoS ONE. 2008;3:e3800.
    [Google Scholar]
  40. , , , , , , , , , . Inhibition of both focal adhesion kinase and insulin-like growth factor-I receptor kinase suppresses glioma proliferation in vitro and in vivo. Mol. Cancer Ther.. 2007;6:1357-1367.
    [Google Scholar]
  41. , , , . CoMFA and CoMSIA studies on HIV-1 attachment inhibitors. Eur. J. Med. Chem.. 2010;45:1792-1798.
    [Google Scholar]
  42. , , , , , , , , . Regulation of heterochromatin remodelling and myogenin expression during muscle differentiation by FAK interaction with MBD2. EMBO J.. 2009;28:2568-2582.
    [Google Scholar]
  43. Mabeta, Peace, PF573,228 inhibits vascular tumor cell growth, migration as well as angiogenesis, induces apoptosis and abrogates PRAS40 and S6RP phosphorylation. Acta Pharm. 66, 2016.
  44. , , , , , , , , , , . A phase Ib dose-finding, pharmacokinetic study of the focal adhesion kinase inhibitor GSK2256098 and trametinib in patients with advanced solid tumours. Br. J. Cancer 2019
    [Google Scholar]
  45. , , , , . A kinetic model of trp-cage folding from multiple biased molecular dynamics simulations. PLoS Comput. Biol.. 2009;5:e1000452.
    [Google Scholar]
  46. , , , , , , , , , , , , , , , , , , , , . Frame MCSpecific deletion of focal adhesion kinase suppresses tumor formation and blocks malignant progression. Genes Dev.. 2005;18:2998-3003.
    [Google Scholar]
  47. , , . Electrostatic effects in proteins: comparison of dielectric and charge models. Protein Eng. 1991:903-910.
    [Google Scholar]
  48. , , , , , , . MMPBSA. py: an efficient program for end-state free energy calculations. J. Chem. Theory Comput.. 2012;8:3314-3321.
    [Google Scholar]
  49. , , , , , . Schlaepfer DDFocal adhesion kinase: in command and control of cell motility. Nat Rev Mol Cell Biol 6:56–68. Nat. Rev. Mol. Cell Biol.. 2005;6:56-68.
    [Google Scholar]
  50. , , , , , , , . Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function. J. Comput. Chem.. 1998;19:1639-1662.
    [Google Scholar]
  51. , , , . Modification of the generalized Born model suitable for macromolecules. J. Phys. Chem. B. 2000;104:3712-3720.
    [Google Scholar]
  52. , , , , , , , , , , . Discovery of Clinical Candidate CEP-37440, a Selective Inhibitor of Focal Adhesion Kinase (FAK) and Anaplastic Lymphoma Kinase (ALK) J. Med. Chem.. 2016;59:7478-7496.
    [Google Scholar]
  53. , . Focal adhesion kinase: the first ten years. J. Cell Sci.. 2003;116:1409-1416.
    [Google Scholar]
  54. PND-1186 FAK inhibitor selectively promotes tumor cell apoptosis in three-dimensional environments, Cancer Biol. Therapy 9, 2010, 764–777.
  55. , , , , , , . 2D and 3D Quantitative Structure-Activity Relationship Study of Hepatitis C Virus NS5B Polymerase Inhibitors by Comparative Molecular Field Analysis and Comparative Molecular Similarity Indices Analysis Methods. J. Chem. Inf. Model.. 2014;54:2902-2914.
    [Google Scholar]
  56. , , , , , , , , . New Insights into the Binding Features of F508del CFTR Potentiators: A Molecular Docking. Pharmacophore Mapping QSAR Analy. Approach, Pharm.. 2020;13:445.
    [Google Scholar]
  57. , , , , , , , , , , . Antitumor Activity and Pharmacology of a Selective Focal Adhesion Kinase Inhibitor, PF-562,271. Cancer Res.. 2008;68:1935-1944.
    [Google Scholar]
  58. , , . Bombesin, Vasopressin, Lysophosphatidic Acid, and Sphingosylphosphorylcholine Induce Focal Adhesion Kinase Activation in Intact Swiss 3T3 Cells. J. Biol. Chem.. 1998;273:19321-19328.
    [Google Scholar]
  59. , , , . Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys.. 1977;23:327-341.
    [Google Scholar]
  60. , . Cellular functions of FAK kinases: insight into molecular mechanisms and novel functions. J. Cell Sci.. 2010;123:1007-1013.
    [Google Scholar]
  61. , , . Evidence for in vivo phosphorylation of the Grb2 SH2-domain binding site on focal adhesion kinase by Src-family protein-tyrosine kinases. Mol. Cell. Biol.. 1996;16:5623-5633.
    [Google Scholar]
  62. , , , , , , , . Arginine in Membranes: The Connection Between Molecular Dynamics Simulations and Translocon-Mediated Insertion Experiments. J. Membr. Biol.. 2011;239:35-48.
    [Google Scholar]
  63. , , , , , . Structural Characterization of Proline-rich Tyrosine Kinase 2 (PYK2) Reveals a Unique (DFG-out) Conformation and Enables Inhibitor Design. J. Biol. Chem. 2009
    [Google Scholar]
  64. , , , , , , , , , , . A first-in-Asian phase 1 study to evaluate safety, pharmacokinetics and clinical activity of VS-6063, a focal adhesion kinase (FAK) inhibitor in Japanese patients with advanced solid tumors. Cancer Chemother. Pharmacol.. 2016;77:997-1003.
    [Google Scholar]
  65. , , , , , . FAK integrates growth-factor and integrin signals to promote cell migration. Nature Cell Biol. 2000;2:249-256.
    [Google Scholar]
  66. , , , , , , , . Insights into the pharmacophore-based 3D-QSAR modeling, molecular dynamics simulation studies of certain dihydroxy pyrrolidine/piperidine and aza-flavanone derivatives as α-glucosidase inhibitors. J. Mol. Struct.. 2020;1223
    [Google Scholar]
  67. , , , , . ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015
    [Google Scholar]
  68. , , , . Accurate calculation of hydration free energies using macroscopic solvent models. J. Phys. Chem.. 1994;98:1978-1988.
    [Google Scholar]
  69. , , , , , , , , , . Biological significance of focal adhesion kinase in ovarian cancer: role in migration and invasion. Am. J. Pathol.. 2004;165:1087-1095.
    [Google Scholar]
  70. , , , , , , , , , , . GPU-accelerated molecular dynamics and free energy methods in Amber18: performance enhancements and new features. J. Chem. Inform. Model. 2018 acs.jcim.8b00462-
    [Google Scholar]
  71. , , , , , , , , , , . FAK Inhibition disrupts a β5 integrin signaling axis controlling anchorage-independent ovarian carcinoma growth. Mol. Cancer Ther.. 2014;13:2050-2061.
    [Google Scholar]
  72. , , , , , , , , , , . CT-707, a Novel FAK Inhibitor, Synergizes with Cabozantinib to Suppress Hepatocellular Carcinoma by Blocking Cabozantinib-Induced FAK Activation. Mol. Cancer Ther. 2016:2916-2925.
    [Google Scholar]
  73. , , , , . Use of MM-PBSA in reproducing the binding free energies to HIV-1 RT of TIBO derivatives and predicting the binding mode to HIV-1 RT of efavirenz by docking and MM-PBSA. J. Am. Chem. Soc.. 2001;123:5221-5230.
    [Google Scholar]
  74. , , , , , . Development and testing of a general amber force field. J. Comput. Chem.. 2004;25:1157.
    [Google Scholar]
  75. , , , , . Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graph. Model.. 2007;25:247-260.
    [Google Scholar]
  76. , , , , , . Identify of promising isoquinolone JNK1 inhibitors by combined application of 3D-QSAR, molecular docking and molecular dynamics simulation approaches. J. Mol. Struct.. 2020;129127
    [Google Scholar]
  77. , , , , , . Discovery of 7H-pyrrolo[2,3-d]pyridine derivatives as potent FAK inhibitors: Design, synthesis, biological evaluation and molecular docking study. Bioorg. Chem.. 2020;102:104092.
    [Google Scholar]
  78. , , , , , , , , . A new force field for molecular mechanical simulation of nucleic acids and proteins. J. Am. Chem. Soc.. 1984;106:765-784.
    [Google Scholar]
  79. , , , . Approximate atomic surfaces from linear combinations of pairwise overlaps (LCPO) J. Comput. Chem.. 1999;20:217-230.
    [Google Scholar]
  80. , , , , . The Collinearity Problem in Linear Regression. The Partial Least Squares (PLS) Approach to Generalized Inverses. Siam J. Sci Statist. Comput.. 1984;5:735-743.
    [Google Scholar]
  81. I. Wolfram Research, Mathematica, Wolfram Research, Inc., Champaign, Illinois, 2018.
  82. , , , . Synthesis, Fungicidal Activity, and QSAR of Pyridazinonethiadiazoles. J. Agric. Food. Chem.. 2002;50:1451-1454.
    [Google Scholar]
  83. , , , , , . Focal Adhesion Kinase Is Upstream of Phosphatidylinositol 3-Kinase/Akt in Regulating Fibroblast Survival in Response to Contraction of Type I Collagen Matrices via a β1 Integrin Viability Signaling Pathway. J. Biol. Chem.. 2004;279:33024-33034.
    [Google Scholar]
  84. , , . Docking and 3D-QSAR studies of 7-hydroxycoumarin derivatives as CK2 inhibitors. Eur. J. Med. Chem.. 2010;45:292-297.
    [Google Scholar]

Appendix A

Supplementary material

Supplementary data to this article can be found online at https://doi.org/10.1016/j.arabjc.2021.103144.

Appendix A

Supplementary material

The following are the Supplementary data to this article:

Supplementary data 1

Show Sections