Molecular dynamics simulation analysis of Focal Adhesive Kinase (FAK) docked with solanesol as an anti-cancer agent

Focal adhesion kinase (FAK) plays a primary role in regulating the activity of many signaling molecules. Increased FAK expression has been associated in a series of cellular processes like cell migration and survival. FAK inhibition by an anti cancer agent is critical. Therefore, it is of interest to identify, modify, design, improve and develop molecules to inhibit FAK. Solanesol is known to have inhibitory activity towards FAK. However, the molecular principles of its binding with FAK is unknown. Solanesol is a highly flexible ligand (25 rotatable bonds). Hence, ligand-protein docking was completed using AutoDock with a modified contact based scoring function. The FAK-solanesol complex model was further energy minimized and simulated in GROMOS96 (53a6) force field followed by post simulation analysis such as Root mean square deviation (RMSD), root mean square fluctuations (RMSF) and solvent accessible surface area (SASA) calculations to explain solanesol-FAK binding.


Background:
Cancer has affected more than 15 million people across the globe and it is expected to affect 25 million people in the next 20 years [1]. Cancer is regarded as the fastest progressing noncommunicable disease with frequent appearance of numerous new forms every year, which are often resistant [2]. Even though massive hard work have been put by researchers, healthcare professionals, etc. in discovering, developing, and utilizing chemotherapeutic approaches for managing neoplasm, but still it continues to be a serious issue across the world at present. The identification of several unexplored classes of cell-cycle regulators and apoptotic stimuli were the emerging strategy in anti-cancer drug discovery. At present, a number of inhibitors have been discovered which modulates cellular resistance, hormonal milieu, angiogenesis, cellular migration, cell proliferation, and DNA targets [3]. The inhibition of angiogenesis (preventing new blood vessels formation) has been representing as an impressive strategic approach in managing cancer and which in turn prevention of metastasis.
Angiogenesis is the process in which new blood vessels originate from the existing blood vessel system [4]. Innumerable proangiogenic factors such as focal adhesion kinase (FAK), TGF-b, vascular endothelial growth factor (VEGF), angiogenin, fibroblast growth factor (FGF), and many other factors, also known as "angiogenic switch" promote constitution of new vessels [5]. In the case of tumors, these factors get over-expressed, that eventually enhances tumor growth and metastasis [6]. From last few years, various researches have revealed successful restriction of tumor proliferation by preventing accessibility of nutrient supply to the highly metabolizing cell like cancer [7]. The inhibition of these pro-angiogenic factors like FAK will lead to a halt in tumor progression.
Focal adhesion kinase (FAK) is a cytoplasmic tyrosine kinase (PTK2) that plays a primary role in growth factor mediated signaling and mediates significant role in cell migration, proliferation, and angiogenesis [8]. Several vascular growth promoting factors like insulin-like growth factor (IGF)-I, vascular endothelial growth factor (VEGF), and basic fibroblast growth factor (bFGF) activates FAK [9]. FAK have also been identified recently to be the prime factor in retinal angiogenesis [10]. Any mutation or alterations in the expression of FAK results in the generation of tumors promote metastasis and endorse vascular growth [11]. FAK has been associated with several forms of cancers like breast, colon, ovarian, prostate, head and neck, thyroid, oral, stomach, cervical, liver, sarcoma, melanoma, and glioblastoma [12]. It has been reported that blockade of FAK often results in reduction of metastasis and mobility in case of breast cancer [13].
Natural products have a huge reputation as promising anticancer agents and in the modern era, naturally derived products are frequently accepted in therapy among the masses [14]. A lot of dietary isoprenoid-based compounds have come into limelight owing to their chemopreventive and antiproliferative properties [15]. Some less known isoprenoid derivatives are now finding applications in mainstay chemotherapy [16]. Solanesol is believed to be among one of the most eligible candidates to demonstrate anti-cancer activity. Commercially, solanesol gets extracted from the tobacco plant, the richest source. As evidenced and rejuvenating the fact from several famous English literature that "Mother Nature inculcates all miraculous remedies in it", our search for an unexplored class ended in tobacco. In one-half, the American Cancer Society (ACS) sticks to the fact that tobacco is the leading cause of cancer. On the other hand, solanesol, the chemopreventive, anti-angiogenic and anti-tumor principle is obtained from the same plant, tobacco, which is quite mysterious and highlights the immense beauty of nature.
Solanesol is a non-cyclic polyisoprenoid alcohol (composed of nine isoprenoid unit) present mainly in solanaceous crops like tobacco, tomato, potato, eggplant, and pepper plants [17]. Solanesol appears as a waxy white solid at room temperature, it is optically inactive, non-polar in nature and exhibit solubility in like-wise non-polar solvents soluble The present research involves establishing of solanesol as a focal adhesion kinase (FAK) inhibitor by applications of computational in silico methods. Though, solanesol used as a bioactive agent in industries for decades, due to its highly flexible nature, there is no successful in silico protein binding and simulation data available online till date. We also proposed a method of binding of the highly flexible compound to protein targets using an enhanced contact based scoring method. This method scores the residues rather than the conformations. The higher scored residues were then used for more "focused" docking on those residue regions.

Methodology: Protein selection and preparation:
The crystallographic co-ordinates for Focal adhesion kinase (PDB ID: 4Q9S) [25] were retrieved from the Protein Data Bank (PDB). Prior to docking, protein structures were prepared by removing water molecules using UCSF Chimera software [26]. Following which, bond orders were assigned, and hydrogen atoms were added to the crystal structures.

Ligand preparation:
Solanesol exist in both cis and trans states, for this experiment we considered only trans is found in natural sources [27]. The structure of Solanesol was obtained from PubChem compound (CID 5477212). Gaussian 09 program was used to obtain the optimum geometry of the structures using the density function theory at the B3LYP/6-31G (d,p) level [28].

Molecular docking:
All the molecular docking studies of Solanesol to FAK were performed using Autodock 4.2 [29]. Autodock uses a semiempirical free energy force field to evaluate binding conformations of ligand while docking. The AutoDockTools was used for preparing protein and ligand parameters files. For ligands; the hydrogens, compute Gasteiger charges, and nonpolar H atom were added which ensuring the total charge corresponds to the tautomeric state followed; and, at last torsion tree root and rotatable bonds were chosen. For macromolecules; hydrogens, compute Gasteiger charges, and merge non-polar H were added which was followed by Stouten atomic salvation parameters assignment; and at last flexible residues PDBQT in addition to the rigid PDBQT file were created. A RMSD value inferior or close to 2A ° was considered as a successful docking [30].

Binding site analysis:
Solanesol is a 45-carbon chain with 26 rotatable bonds. As it is extremely flexible it is hard to determine the bind mode of it with the protein. The commonly used protocol for determination of binding pocket is "Blind docking", which was initially developed for to determining peptide docking with protein [31]. In this method the constrained ligand (or peptide) is docked with the whole protein surface. The place where it forms a cluster with higher energy determines the binding site. Then these sites were used for "refined docking" where the lowest binding modes for each of these places (in case if there are more than one) where determined by molecular mechanics and molecular dynamics studies. For calculating the possible area of interaction or binding site of a highly flexible ligand we enhanced blind docking using Ligand Contact Based Scoring function for Residues (LCBSR). This is an atomic contact/clash based scoring method, in which the residues are scored based on the higher favorable interactions and probability of formation of a hydrogen bond. As it doesn't depend on the clustering or solely on binding energy, it statistically enhances the probability of finding the possible binding site.
It can be represented as, The "Contact" here is defined as the instance when the difference between the distance of two atoms and the sum of their van der Waal radii is 0.4 A or more or, in other words, the distance is greater than the sum of van der Waal radii (Eq. 1) of two atoms. Whereas "Clash" is defined as the condition where the van der Waal radii of two atoms unfavourably overlaps each other and the distance is lower than the sum of radii (Eq. 2) of the two atoms.
This can be represented as: Where, ΣrVDW(i,j) is the sum of van der Waal radii of interacting atoms of ligand (i) and residue (j) and Dij is the distance between interacting atoms of ligand (i) and residue (j). A Higher value of LCBSR of any residue implies the residue may be a part of the binding pocket for the ligand and, if it is capable, it may also form a hydrogen bond with the ligand. Lesser score implies a lower interaction or high chances of unfavorable clashes or low chances of forming a hydrogen bond. Considering all the conformers may lead to false positives, thus conformers were separated in three criteria: (1) binding energy less than -2.0 kcal/mol; (2) binding energy less than -3.0 kcal/mol; and binding energy less than -4.0 kcal/mol. Provided -4.0 is roughly half the average of the binding energy (Average B.E = -7.10 kcal/mol) of all the three experiments, ensuring only conformations with low B.E were considered. The resultant values from each of these three were added to corresponding residues. This way the residues interacting more with the conformations of better binding energy will have a better score. For getting a statistically significant result, all the scores of the residues from all the three experiments were added to get the final score of each residue. Only those residues, which appeared in more than 2 experiments, were considered.

Refined docking:
Binding site for Solanesol was considered using residues with a higher LCBSR score. This binding site was then used for refined docking using Autodock. The experiment was done twice with (1) relaxed parameters, GA maximum energy evaluations 2.5 × 106, for 200 GA runs (2) exhaustive parameters, GA maximum energy evaluations 3.5 × 107, for 200 GA runs.      method was used, with a coulomb cutoff of 1.2 nm, Fourier spacing of 0.16 nm and an interpolation of order 4. Energy minimization of the system was carried out using steepest descent algorithm with a tolerance value of 1000 kJ mol -1 nm -1 . After energy minimization, NVT and NPT equilibrations were done on the system until it reached the room temperature and water density. Production MD was performed for 20 ns time duration for both the simulations.

Molecular dynamics trajectories analysis:
The root mean square deviation (RMSD) and root mean square fluctuations (RMSF) of FAK backbone were calculated using "g_rms" and "g_rmsf" utility commands, respectively. A spherical probe of radius 1.4 Å across the protein surface was used for calculating solvent-accessible surface area (SASA) by "g_sasa" tool of Gromacs. Hydrogen bonds between Solanesol and FAK were calculated using "g_hbond" tool with proton donor and acceptor distance ≤ 3.5 Å and the angle between acceptor-donor-hydrogen ≤ 30.0 degrees.

Binding free energy calculations:
The molecular mechanics Poisson-Boltzmann surface area (MMPBSA) [38] approach was to estimate the binding free energy of protein-ligand interaction. For this purpose, "g_mmpbsa" [39] tool was used. The tool calculates the molecular mechanics potential energy and the free energy of solvation and excludes the entropy calculations. MM-PBSA calculations were performed using 1000 snapshots taken from last 5 ns of trajectories of the complex system. The MM-PBSA based binding affinity (ΔΔG) was calculated using the g_mmpbsa provided script.
ΔΔGBE=ΔGComplex-(ΔGReceptor + ΔGLigand) ΔEMM, ΔGSol and TΔS represent the molecular mechanics component in the gas phase, stabilization energy due to salvation, and a vibrational entropy term, respectively. ΔEMM is the summation of ΔEint, ΔEcol, and ΔEvdw, which are the internal, coulomb, and van der Waals interaction terms, respectively. ΔGsol is the salvation energy and it is divided into an electrostatic salvation free energy (GPB) and a non-polar salvation free energy (GSA).

Result and Discussion: Docking analysis Blinding docking and LBCSR Score:
Solanesol shows very small cluster with insignificant binding affinity towards FAK when docked "blindly" (Table 1). Though, the structure shows quite high binding affinity towards the kinase but the conformations with higher binding energy fails to form any significant cluster.
From all the generated conformations of Solanesol, which were having binding energy lower than -4.0 kcal/mol, all the favorable and unfavorable overlaps of the atoms were calculated using "FindClash" tool of UCSF Chimera. Chimera tool, "FindHbond" was also used to find the hydrogen bonds between the ligand and residue atoms. The default values were kept for all the calculation in Chimera.
An in-house python script was used to calculate the number of "Clashes", "Contacts" and "Hbonds" between all the conformations and residues as well as individual scores. The scores from all the three experiments were added to give a final score for each residue (Table 2).
Based on the LBCSR score calculated for all three experiments, the scores for all the residues were plotted as graph (Figure 3) as well as plotted as false color on the 3D structure of the protein (Figure 2). From the LBCSR score it can be inferred that the ligand shows a high affinity towards the region with Asp564, Asn551, Arg550, Leu553 and Ile428, respectively.

Refined docking:
All the residues with more than LBCSR score was considered for the active site prediction. A total of 77 residues were found to be above and incidentally which also forms the ATP binding site and the catalytic loop (546-551) and formed between the N and C lobe.
The centre of geometry (coordinates = 10.8, 1.0, 15.0) of these 77 residues was considered for the centre of the binding pocket of Solanesol. A grid of size 60 × 62 × 68 was considered to exactly fit all the 77 residues. Autodock4.2 was again used for docking of Solanesol with FAK with this grid setting for two more times, first time with default setting for 200 conformations and later docked with exhaustive setting for the same number of conformations. All 400 conformations were rescored using DSX online server with CSD settings ( Table 3).

Analysis of Final Docked structure:
The hydroxyl end of Solanesol binds to the binding site of ATP and interacts with Ile428, Val436, Ala452, Lys454 and Leu501 (Figure 4). These residues that forms the binding pocket for ATP, forms Alkyl-Pi hydrophobic interactions with the double bonds of the ligand. Where the isobutyl end of the ligand gets attached to the Catalytic loop and αC helix. The ligand interacts with Phe542, Arg545, Arg550 (Catalytic loop) and Phe478 (αC helix terminal) with Alkyl and Sigma-Pi interactions. Pro585, Ile586 near the ATP active site also forms similar hydrophobic interactions with the ligand. Oxygen of Solanesol forms a very conventional H-bond with OE1 of Gln438, which is very near to G-Loop and ATP binding site.
These 12 hydrophobic interactions of 11 residues with Solanesol stabilize the ligand at the middle of N and C lobe of FAK. As the ligand share interacting residues with both the side it may be act as a better inhibitor for FAK. This structure was used for molecular dynamics studies.

Analysis of Molecular Dynamics:
For the measuring the stability, the RMS deviation in the backbone of the Solanesol bound FAK and only FAK been measured and plotted in accordance to time ( Figure 5). After an initial instability, the backbone changes its shapes linearly with a linear slope increase in RMSD between 4.5ns to 7.7ns after again a short destabilization the system vaguely equilibrates from 9.4ns to 14.2ns. It takes the system almost 15ns to equilibrate, after which the system maintains it position and shape.
Comparison of the RMSD of both the trajectories from the minimized structures shows, Solanesol bound FAK backbone show more stability than that of independent FAK backbone. RMSF (Root mean Square Fluctuation) of Solanesol bond FAK exhibits less fluctuation at the places where Solanesol is bond ( Figure 6).
Total solvent accessible surface area (SASA) was checked by g_sas tool of Gromacs for analyzing the change in surface area with respect to time. It show a quite negative correlation with the RMSD suggesting the better binding of the ligand leads to a lower surface area of the protein thus closing the active site for any further contact. It remains stable for the wider part of 6ns -18ns range after which the backbone gets stabilized and may also affect the surface area.

MM-PBSA Calculation:
Binding site residues for Solanesol were selected by taking 3.5A radius from Solanesol. Molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) was calculated for the last 5 ns with 20 ps steps for the binding site residues versus Solanesol using g_mmpbsa tool. Per residue analysis of the result was done and plotted using the script provided. This analysis suggests that, Ile428, Gly431, Val436, Val484, Met499, Cys502 and Leu553 interacts most favorably with the ligand (Figure 7). Interestingly these all residues are part of ATP binding pocket of FAK. Gly431 is a part of the G-loop, which helps in the phosphorylation of the protein, also interacts favorably along with Leu584, which is part of the activation loop. A very low ΔΔG, of -113.85 kJ/Mol (Table  4), proves solanesol have a very high binding affinity towards the FAK structure. The change in binding energy with time was also plotted (Figure 8). It shows that the ligand gets stabilized with a very high binding affinity after 17 ns of simulation.

Conclusion:
We report the binding of Solanesol with FAK using a GROMOS96 (53a6) force-field simulated docked model of Solanesol-FAK complex. The binding of solanesol at the ATP binding site to inhibit the phosphorylation of FAK is explained. Data help to understand binding of flexible compounds with FAK for potential inhibition.