Molecular modelling approaches predicted 1,2,3-triazolyl ester of ketorolac (15K) to be a novel allosteric modulator of the oncogenic kinase PAK1

P21-activated kinases (PAKs) are serine/threonine protein kinase which have six different isoforms (PAK1–6). Of those, PAK1 is overexpressed in many cancers and considered to be a major chemotherapeutic target. Most of the developed PAK1 inhibitor drugs work as pan-PAK inhibitors and show undesirable toxicity due to having untargeted kinase inhibition activities. Selective PAK1 inhibitors are therefore highly desired and oncogenic drug hunters are trying to develop allosteric PAK1 inhibitors. We previously synthesized 1,2,3-triazolyl ester of ketorolac (15K) through click chemistry technique, which exhibits significant anti-cancer effects via inhibiting PAK1. Based on the selective anticancer effects of 15K against PAK1-dependent cancer cells, we hypothesize that it may act as an allosteric PAK1 inhibitor. In this study, computational analysis was done with 15K to explore its quantum chemical and thermodynamic properties, molecular interactions and binding stability with PAK1, physicochemical properties, ADMET, bioactivities, and druglikeness features. Molecular docking analysis demonstrates 15K as a potent allosteric ligand that strongly binds to a novel allosteric site of PAK1 (binding energy ranges – 8.6 to – 9.2 kcal/mol) and does not target other PAK isoforms; even 15K shows better interactions than another synthesized PAK1 inhibitor. Molecular dynamics simulation clearly supports the stable binding properties of 15K with PAK1 crystal. Density functional theory-based calculations reveal that it can be an active drug with high softness and moderate polarity, and ADMET predictions categorize it as a non-toxic drug as evidenced by in vitro studies with brine shrimp and fibroblast cells. Structure–activity relationship clarifies the role of ester bond and triazol moiety of 15K in establishing novel allosteric interactions. Our results summarize that 15K selectively inhibits PAK1 as an allosteric inhibitor and in turn shows anticancer effects without toxicity.

which defines a quantitative classification of global electrophilic nature of a compound 18 .

Molecular docking analysis.
Three different crystal structures of PAK1 (PDB ID: 3FXZ, 3Q52, and 5DEW) were retrieved from the RCSB protein data bank (PDB) based on the lowest resolution. PAK1 crystal structures 3FXZ, 3Q52, and 5DEW were resolved at 1.64, 1.80, and 1.90 Å, respectively. Proteins were then prepared for docking with removing water, heteroatoms, and co-crystallized ligands using Discovery Studio 4.5 (Accelrys, San Diego, CA, USA). To prepare the ligands for docking analysis, the drugs were saved as protein data bank (pdb) format after minimizing the energies with DFT using B3LYP/6-311G (d,p) level of theory. Molecular docking was employed with AutoDock Vina suite which is widely used for molecular docking as it significantly improves the accuracy of the binding mode predictions compared to its previous version AutoDock 4 19 . The grid box was prepared in the hinge region of PAK1 surrounding the active site residues, allosteric residues, gatekeeper residue, and DFG-motif residue 20,21 . The size of the grid boxes was differed with crystal structures, and the dimension and the center of the grid boxes for each protein structure were presented in Table 1.
The binding affinity of ligands was observed by kcal/mol as a unit for a negative score 19 , and the binding pose with the highest negative value was selected as the best pose for corresponding ligand. Then, the best pose was visualized and analyzed to comprehend the ligand-protein interactions through PyMOL Molecular Graphics System 2.0 (DeLano Scientific LLC, San Carlos, CA, USA) and Discovery Studio 4.5. Predicted protein-ligand interactions were validated by cross-checking with LigPlot + software (obtained from https:// www. ebi. ac. uk/). Three crystal structures of PAK4 (5VEF, 5XVA and 5ZJW) and PAK6 (2C30, 4KS7 and 4KS8) were also retrieved from the RCSB PDB and used them for molecular docking simulation against 15K using the similar method described above. However, in PDB, there are no full crystal structures of PAK2, PAK3, and PAK5. Therefore, homology modelling was carried out using SWISS-MODEL server (https:// swiss model. expasy. org/ inter active) to predict 3D structures for PAK2, PAK3, and PAK5. Their FASTA sequences were copied from the UniProt protein data bank (https:// www. unipr ot. org), under the access code Q13177 (for PAK2), O75914 (for PAK3), and Q9P286 (for PAK5). PAK2, PAK3, and PAK5 consist of 524, 559, and 719 amino acid sequences, respectively. The crystal model for each protein was built based on the template showing the highest sequence identity and query coverage. Ligands were removed from the models and hydrogen atoms were added using Discovery Studio 4.5. Energy minimization was done for these three selected models with Swiss-PdBViewer prior to using for molecular docking with the similar methods described earlier.

Molecular dynamics simulation.
For validating the binding behavior of 15K with PAK1 and comparing with the same shown by inhibitor 2, molecular dynamics (MD) simulations were performed with YASARA Dynamics suite 22 using AMBER14 force field. NPT ensemble simulation was used with normal temperature (298 K) and 1 bar of pressure. Berendsen thermostat was used to control the simulation temperature, and for short-range van der Waals and Coulomb interactions, a cut-off radius of 8.0 Å was considered. The long-range electrostatic interactions were calculated using the particle-mesh Ewald method. The periodic boundary condition was included for performing the simulation, where the cell size was 20 Å larger than the protein in all cases. These simulations were performed in water which involved the simulation of predicted model inside a trajectory box filled with solvent (including NaCl and water) and protein-ligand complex at 298 K. In order to reduce the contact area difference between the complex model and the solvent molecule, the initial model was minimized with steepest-descent algorithm. As time step, 1.25 fs was chosen, and 100 ns MD simulation was performed for each system where trajectories were saved after every 100 ps. The MD trajectories were analyzed by a macro program written in YANACONDA language. In this analysis, the root mean square deviation (RMSD) was computed for the protein backbone and residues with a view to checking the stability of the trajectories. In order to study the flexibility of the trajectories, the root mean square fluctuation (RMSF) was calculated per residue. Further analysis was performed by calculating the radius of gyration (Rgyration) and solvent accessible surface area (SASA) of the enzyme within the 100 ns of production time.  23,24 , which manifests the variability exist in MD trajectories. The bond distances, bond angles, dihedral angles, planarity, van der Waals energies, and electrostatic energies represented structural and energy profile. The PCA analysis assisted to reflect the variations among different groups. The 100 ns trajectory data were pre-processed using centering and scaling for the analysis. The multivariate factors were arranged in the X matrix and expressed them into a product of two new matrices by the following equation: where Tk, matrix of scores symbolizes how sample is relative to each other; Pk, matrix of loadings manifests how variables correlate to each other, k is the number of factors integrated in the model and E is the matrix of residuals. The analysis was performed using R programming language and the plots were visualized by ggplot2 package 25 .
Multiple receptor conformers-based molecular docking. Two methods were used to generate multiple conformers of PAK1 crystal. Firstly, the crystal structure of PAK1 (PDB ID: 3Q52) was used in 100 ns MD by YASARA Dynamics program as stated above. Then, 100 snapshot conformers were generated and saved in PDB format after every 5 ns of 100 ns MD simulation. In another approach, 20 crystallographic structures of PAK1, except for 3Q52, were collected from the RCSB PDB database (PDB ID: 1YHV, 1YHW, 3FY0, 3Q53, 4DAW, 4EQC, 5DFP, 3Q4Z, 4O0R, 4O0T, 4P90, 4ZJI, 4ZJJ, 4ZLO, 4ZY4, 4ZY5, 5DEY, 6B16, 5KBR, 5IME). Flexible docking with 15K against the multiple conformers of the PAK1 was performed and subsequently analyzed using the same protocol by AutoDock Vina. Inhibitor 2 was also docked against multiple conformers in this experiment to compare its binding affinity with that of the 15K.
Prediction of ADMET properties. Pharmacokinetic parameters related to drug absorption, distribution, metabolism, excretion, and toxicity of 15K were predicted with pKCSM online tool 26 and compared to the predicted values of inhibitor 2. pKCSM is a novel approach to predict pharmacokinetic and toxicology parameters using graph-based signatures, which develop and implement 14 quantitative regression models with actual numeric outputs and 16 predictive classification models with categorical outputs for predicting a wide range of ADMET properties for novel diverse molecules. We also used ADMETlab 2.0 (https:// admet mesh. scbdd. com/) to compare the ADMET results.
Prediction of bioactivities, physicochemical properties and druglikeness of 15K. Prediction of activity spectra for substances (PASS) program was utilized to predict biological activities of 15K and inhibitor 2. It is an online-based computer program for the evaluation of overall biological potentials according to the structure-activity relationship (SAR) 27 . Through this program, prediction is carried out by comparing the desired molecules with a training set including more than 205,000 compounds revealing more than 7200 types of biological activities 28 . This tool predicts a compound with its probability to be active (Pa) and a probability to be inactive (Pi). The easy, efficient, and login-free web tool SwissADME 29 was here used to predict the physicochemical properties and druglikeness of 15K.
In vitro toxicity evaluation of 15K. Cytotoxic effects of 15K were evaluated with two different approaches.
Firstly, using the normal fibroblast 3T3L1 cells, cytotoxicity was investigated with MTT assay 30 . The cells were seeded for overnight in a 24-well plate with the density of 2.5 × 10 4 cells/well and then treated with 15K at different concentrations ranging up to 20 µM for 48 h. After that, the cells were washed with ice-cold PBS, and an aliquot of MTT solution (0.5 mg/ml in PBS) was added to each well. After incubating the cells at 37 °C in humidified conditions, formazan crystals formed in each well were dissolved with 400 μL dimethyl sulfoxide (DMSO), and the plate was shaken for 10 min at room conditions. The absorbance was measured at 570 nm wavelength using a microplate reader. Cytotoxicity was calculated based on the absorbance of the treated versus untreated control cells. In the second approach, toxic effects of 15K were evaluated using a small animal model brine shrimp (Artemia salina) 31 . Brine shrimp eggs (200 mg, Tetra Brine Shrimp Eggs, Spectrum Brands Holdings, Inc., Kanagawa, Japan) was hatched in artificial seawater supplemented with dried yeast with suitable aeration in room conditions for 48 h. Active brine shrimp nauplii were then transferred to a glass petri dish-containing artificial seawater and then transferred to each well of a 96-well plate. Each well was contained 10-15 larvae in 100 μl artificial seawater. 15K dissolved in DMSO was added in each well except for the control wells in different concentrations, and then the well plate was incubated at 20 °C. The final concentration of DMSO in the treated wells was 1% (v/v), and DMSO at 1% was used in the control wells. The number of dead larvae in each well was counted using magnifying lenses (5 × or 10 ×) after 24 h and 48 h of incubation. Finally, methanol (50 μl) was added to each well, and the total number of larvae in each well was counted after 1 h. The survival percentage of brine shrimps for each concentration was then calculated.
Protein-protein interaction (PPI) network and enrichment analysis. PPI network was prepared using the default options of STRING database and then it was imported to the Cytoscape for further analysis. Total number of nodes and edges were 11 and 50, respectively. Enrichment analysis was also done with the STRING database, which shows a total of 235 gene ontology (GO) terms (component/process/function) enrich-

Results
Geometry optimization and the features of frontier molecular orbitals. Full geometry optimization of ketorolac, 15K, and inhibitor 2 was carried out by DFT calculation, and their most stable optimized structures are presented in Fig. 1. Table 2 represents the electronic energy, enthalpy, Gibbs free energy, and dipole moment of the compounds. 15K shows different structural and thermodynamic properties compare with its precursor ketorolac. It exhibits higher electronic energy, enthalpy, and Gibbs free energy than ketorolac. Dipole moment is also higher in 15K than ketorolac and inhibitor 2, demonstrating that the polarity of 15K is higher than both other drugs. Frontier molecular orbital analysis reveals that 15K has the lower HOMO-LUMO energy gap (0.144 Hartree) compared with ketorolac and inhibitor 2, and in turn 15K shows low hardness and high softness value (Table 3). Figure 2 represents the molecular orbital distribution in 15K and the energy gap between HOMO and LUMO. The electron distribution in both HOMO and LUMO is observed in the precursor molecule not in its triazole moiety.
Molecular electrostatic potential is highly informative concerning the nuclear and electronic charge distribution of the molecules. It is a tool for interpretation and prediction of chemical reactivity. We calculated molecular electrostatic potential (MEP) to predict the reactive sites for electrophilic and nucleophilic attack of all three optimized structures. Red color indicates maximum negative area (favorable site for electrophilic attack), blue color represents the maximum positive area (favorable site for nucleophilic attack) and zero potential area is represented by green color (Fig. 3). Here, 15K potentially shows negative − 0.2689 a.u. and positive + 0.1745 a.u.  www.nature.com/scientificreports/

Molecular interactions of drugs with different PAK isoforms. Molecular docking analysis shows
that 15K can strongly bind with all the three PAK1 crystals (3FXZ, 3Q52, and 5DEW) compared with ketorolac and inhibitor 2. The binding affinities for ketorolac were ranging − 6.7 to − 7.0 kcal/mol, while the binding affinities for 15K were − 8.6 to − 9.2 kcal/mol (Table 4). 15K exhibited higher binding energies than those shown by inhibitor 2.
15K interacted with all three PAK1 crystals by forming three hydrogen bond interactions, while ketorolac formed only one hydrogen bond with both 3Q52 and 5DEW but no hydrogen bond with 3FXZ ( Fig. 4). However, inhibitor 2 interacted differentially with different crystals forming one hydrogen bond with 3FXZ and 5DEW and four hydrogen bonds with 3Q52. In every case, 15K formed two hydrogen bonds with the αC helix residue Glu315 having the bond distance 2.5 Å or less. It also showed another hydrogen bond interaction with the DFG residue Asp407 of 3FXZ. One electrostatic interaction with Glu315 was observed in 15K-5DEW complex. The hydrogen bond interactions observed by ketorolac were not with any active site or even allosteric residues, while it showed an electrostatic interaction with Glu315 of 5DEW. Inhibitor 2 showed a hydrogen bond interaction with Glu315 of 5DEW crystal and DFG residue Asp407 of 3Q52. Both ketorolac and inhibitor 2 formed hydrophobic interaction with the gatekeeper residue Met344 of two crystals, however 15K only showed this interaction with one crystal. Moreover, hydrogen bond surface shows that the residues such as Ser281, Arg299, Asp407, Gly279, and Ala280 can create strong donor regions on the drug-protein interaction surface (Fig. 5).
When 15K was docked against other PAK isoforms, PAK2-6, it did not show stronger binding affinities than what it showed for PAK1. It only showed better binding energies with the modelled PAK2 crystal (− 8.4 kcal/mol) and one PAK4 crystal (− 8.3 kcal/mol) (Tables 5 and 6). The binding affinities with other PAK isoform crystals were low (− 0.8 to − 7.7 kcal/mol). Most importantly, 15K did not show any hydrogen bond interaction with the active site residues or allosteric residues of PAK2-6, even no hydrophobic and electrostatic interactions were observed with those residues. Ligand-protein molecular interaction results from LigPlot + coincide with those obtained from Discovery Studio (Fig. 4E-G).
Analysis of molecular dynamics simulation. We carried out MD simulation and explored the stability of each MD system to determine the promising PAK1 inhibitors. MD simulation for each complex of PAK1 with two ligands (15K and inhibitor 2) is performed for 100 ns. We also conducted an MD run for PAK1 apo which shows stability for RMSD value during the whole run (average 1.6 Å). Although the RMSDs (0.4-2.7 Å) for α-carbon atoms is slightly higher in 15K than inhibitor 2 (0.3-2.5 Å), our results still suggest that 15K can be stable in physiological conditions. RMSD values of 15K are slightly increased to 2.7 Å after 75 ns and show www.nature.com/scientificreports/ few smaller fluctuations at 35 and 73 ns (Fig. 6A). The trajectories generated otherwise during the whole run are found to be stable. Thus, higher structural stability in the PAK1-15K complex is detected. Inhibitor 2 exhibits high fluctuations from 73 to 74 ns, but is stable in the later phase with the average of 1.7 Å. The average RMSD for 15K is 1.9 Å, indicating that it can be a stable and potent inhibitor. The Rgyration is a parameter that indicates structural compactness of protein. A lower degree of fluctuation with its consistency throughout the simulation run suggests the greater compactness and rigidity of the system 32 . A similar Rgyration is detected for apo (average 20.11 Å), 15K (average 20.05 Å), and inhibitor 2 (average 19.99 Å), which suggests comparatively tight packing of the ligand-protein complex during the whole run. This result indicates that the complexes are relatively compact and not much changed after ligand binding (Fig. 6B). The expansion of protein volume is indicated by higher SASA value and a low fluctuation is expected over the simulation time. SASA shows the lowest value for 15K-PAK1 complex (average 14,013 Å 2 ) compared to inhibitor 2 (average 14,106 Å 2 ) (Fig. 6C). Molecular surface area (MolSA) was evaluated for both ligand-protein complexes. Inhibitor 2 complex shows the highest MolSA (average 14,649 Å 2 ). 15K and apo show slightly lower value of MolSA with the average of 14,571 Å 2 and 14,339 Å 2 , respectively (Fig. 6D).
RMSF property is highly used to capture the conservation of protein dynamics 33 . During the protein-ligand interaction, a significant fluctuation of RMSF values of residues (415-420) is observed (Fig. 6E). In 15K-PAK1 complex, most of the residues show a low RMSF (average 1.29 Å) than inhibitor 2-PAK1 complex (average 1.34 Å) or the apo form (average 1.31 Å). The RMSF of hotspot residues Glu315, Met344, and Asp407 in 15K-PAK1 complex is lower than inhibitor 2-PAK1 complex, reflecting that 15K can stimulate more strong interactions with these key residues 34 . Furthermore, the visual analysis of MD simulation trajectories suggests that both ligands engaged in significant binding interactions with the hotspot residues mentioned earlier.
Stability of the protein molecule largely depends on hydrogen bonds. The intermolecular hydrogen bond formed during the whole 100 ns run is observed and the maximum number of hydrogen bond is predicted for 15K (average 531), whereas the number of hydrogen bond in inhibitor 2-PAK1 and apo is almost similar (523 and 521) (Fig. 6F). The snapshot of conformers displays that both ligands were mostly in the hinge region of PAK1 throughout the entire simulation process (Fig. 7A,B).
Principal component analysis on MD simulation. The PCA model was built to evaluate the structural and energy data retrieved from MD simulation analysis. There are three training sets accounted for the analysis. For the PCA model, total 97.1% of the variance is expressed by PC1 and PC2, wherever PC1 exposes 80.5% and PC2 exposes 16.6% of the variance (Fig. 8). The score plot shows different clusters for apo-PAK1 (Blue), 15K-PAK1 complex (Red) and inhibitor 2-PAK1 complex (Green). The fluctuations during MD simulation are responsible for different cluster formation among three groups. Here, 15K-PAK1 complex exhibits significant www.nature.com/scientificreports/ differences than both apo-PAK1 and inhibitor 2-PAK1 complex. The divergent behavior of 15K-PAK1 complex infers that the complex has exerted substantial variations in the structural and energy profile, which may favor the strong interaction of 15K with PAK1. www.nature.com/scientificreports/ Multiple receptor conformers based molecular docking. Two types of PAK1 conformers generated from MD simulation and PDB were used to detect the difference in binding affinity of 15K and inhibitor 2 ( Fig. 9). Figure 9A represents that 15K exhibits more binding affinity than inhibitor 2 against 20 different MD trajectories (0.1, 5, 10, 15,… 0.100 ns) obtained from 100 ns run. Inhibitor 2 demonstrates an average binding affinity of − 7.9 kcal/mol, whereas 15K shows − 8.2 kcal/mol. Similarly, the binding affinity was higher for 15K than that for inhibitor 2 when they were docked with 20 crystal structures of PAK1 other than 3Q52 (Fig. 9B).

ADMET, physicochemical, and pharmacological properties, and druglikeness of 15K.
To screen the most effective compound during drug development process, different factors like absorption, distribution, metabolism, excretion, and toxicity (ADMET) play important roles. Through an in silico ADMET analysis, 15K shows poor aqueous solubility (log mol/L = − 4.26) (Table 7), indicating poor absorption via oral administration. The ADMET table (Table 7) also expresses that 15K could be distributed in blood stream only and is not able to cross the blood brain barrier (BBB) and in turn it could not pose any psychological side effect. www.nature.com/scientificreports/ In terms of metabolism, 15K has been found to be a non-substrate and non-inhibitor of CYP2D6 enzyme. A noninhibitor of CYP2D6 means that it wouldn't interfere in any kind of biotransformation of any compound metabolized by CYP2D6 enzyme. However, it has been predicted to be a substrate and inhibitor of CYP3A4 enzyme demonstrating that it might interfere the metabolism of drugs that are metabolized by this enzyme. Similarly, it may work as the inhibitor of CYP1A2, CYP2C19, and CYP2C9 enzymes. Predicted total clearance of 15K implies moderate level of clearance through kidney, liver, and lungs. This result also indicates a medium level of protein binding capacity of this drug.
Toxicity of 15K was predicted based on the AMES toxicity test, hERG receptor binding, and skin sensitivity test, and the predicted results demonstrate that it can be a non-mutagenic drug without showing toxicity to the heart and irritation to the skin. Surprisingly, most of the predicted ADMET values from pKCSM server comply with the results obtained from ADMETlab 2.0.
In silico pass prediction server was used to predict the highest probable pharmacological properties of 15K and the results show that it can be highly active against hypoxia-inducible factor 1A ( Table 8). The other pharmacological properties include anti-inflammatory and analgesic effects. Predicted physicochemical properties support 15K as a good orally bioavailable drug (Table 9) which follows all the druglikeness rules introduced by Lipinski et al. 35 , Ghose et al. 36 , Veber et al. 37 , Egan et al. 38 , and Muegge et al. 39 .
In vitro toxic effects of 15K. Toxicity of 15K was evaluated by using brine shrimp lethality test. 15K was tested using 4 different concentrations and as shown in Fig. 10, it did not exhibit toxicity against brine shrimp at the highest 200 nM concentration after 24 and 48 h treatment. It also did not induce notable toxicity up to 20 nM on 3T3L1 fibroblast cells, while it significantly inhibited cells viability at 50 nM. PPI network and enrichment analysis. The prepared PPI network by PAK1 shows that it consists of 11 nodes and 50 edges (Fig. 11). Five nodes-PAK1, Rac1, Cdc42, Ras homolog family member A (RHOA), and cortactin (CTTN) are significantly involved in this network as proved by their highest number of degree (10). The degree value for all other nodes except for filamin-A (FLNA) is 9. PAK1, Rac1, Cdc42, RHOA, and CTTN also show the highest value for betweenness centrality and stress, while for other proteins both the betweenness centrality and stress are zero. The closeness centrality for that five proteins is also higher than those for other proteins and the lowest closeness centrality is shown by FLNA.
GO enrichment analysis based on 11 targets from the predicted PPI network demonstrates that most of the GO components/processes/functions are relevant to the cellular structure organization and signal transduction regulation. The highly enriched GO terms include actin cytoskeleton organization, cell projection organization, www.nature.com/scientificreports/ regulation of signal transduction, plasma membrane bounded cell projection, plasma membrane bounded cell projection organization, lamellipodium, etc. (Fig. 12A). All these enriched GO terms can eventually enrich several pathways in cancer proved by the KEGG enrichment analysis prediction (Fig. 12B). The highly enriched pathways are regulation of actin cytoskeleton, Ras signaling pathway, mitogen-activated protein kinase signaling pathway, focal adhesion, exon guidance, etc.

Discussion
Protein kinases catalyze phosphorylation reaction and activate a wide array of signaling pathways responsible for cell growth, survival, differentiation, and migration. Small molecule kinase inhibitors therefore attracted attention as promising cancer therapeutics, but most of the developed inhibitors are ATP-competitive that may target other kinases due to having conserved binding features. Thus, cancer drug hunters started looking for non-ATP competitive allosteric inhibitors which have significant advantages over ATP-competitive inhibitors with greater specificity and lower off-target toxicity. In this study, we evaluated 15K, which is a PAK1 inhibitor, using computational techniques and found that it can selectively bind with the allosteric site of PAK1 without targeting other PAK isoforms (PAK2-6) and thereby not showing any toxicity to human body so does pan-PAK inhibitors. Thermodynamic properties of 15K show that it has greater negative values for electronic energy, enthalpy, and Gibbs free energy than those shown by ketorolac, indicating that 15K has improved thermodynamic features compared with ketorolac 30 . In addition, the higher dipole moment demonstrates high polarity features of 15K than its precursor ketorolac and the tested compound inhibitor 2. Dipole moment is an indicator of drug-receptor interaction, which may facilitate hydrogen bond interactions with the target protein 40 . Molecular orbital properties reveal that the HOMO-LUMO gap is lower in 15K and thus its hardness is low, and softness  www.nature.com/scientificreports/ is high compared to both ketorolac and inhibitor 2. Drugs with higher softness are highly desirable as they can be degraded promptly after showing their pharmacological effects without exerting any unwanted toxicity 41 . Our group synthesized 15K using "Click Chemistry" technique and proved it to be a more water-soluble and cell permeable azo derivative of ketorolac 9 . Click chemistry technique was introduced by Kolb et al. 42 , based on the principal that coupling of 1,2,3-triazole ring with the carboxylic group in presence of copper catalyzer forms a water-soluble ester. In silico ADMET analysis reveals that 15K can be a poorly water-soluble drug with the log mol/L value of − 4.26. ADMET analysis predicts that it, after ingestion, can be transmitted through intestinal barrier via passive transcellular diffusion. Using multi-drug resistant breast cancer cell line-EMT6, we proved that 15K has the higher cell permeability than ketorolac 9 , which is also justified by the calculated polar surface area (PSA) (88.24 Å 2 ). A small molecule drug with less than 100 Å 2 PSA can show good membrane permeability features 43 . However, it cannot pass the BBB due to having a higher molecular weight than 400 Da 44 ; our ADMET analysis also predicted the same. The BBB is a dynamic hindrance ensuring the brain against undesirable substances, and 15K would therefore not show any physiological side effects.
Molecular interaction analysis explores that 15K shows strong binding affinity with PAK1 crystal (binding energy − 8.6 to − 9.2 kcal/mol), forming the highest number of hydrogen bond interactions (three) regardless of the PAK1 crystal structures, while ketorolac showed a varied interaction with different PAK1 crystals. Most importantly it showed two hydrogen bond interactions with the allosteric residue Glu315, similar to inhibitor 2. Several previous reports have already proved that Glu315 is a novel allosteric residue in αC helix, which is responsible for binding the selective PAK1 inhibitors 14,21 . 15K also can strongly interact with DFG motif residue Asp407 which is another allosteric residue adjacent to the ATP-binding site 14 . During our analysis, 15K did not show notable interaction with other PAK isoforms (PAK2-6). Although the binding affinities with a PAK2 and PAK4 crystal were a bit higher than the normal range, it didn't show any hydrogen bond interaction either with the active site residues or allosteric residues, which clarifies its weak interaction with PAK2 and PAK4 and thereby demonstrating that 15K has selectivity towards PAK1. When 15K was docked with multiple PAK1 conformers, it exhibited stronger binding affinities than inhibitor 2.   www.nature.com/scientificreports/ Drug-protein complexes for both 15K and inhibitor 2 were chosen for MD simulation to investigate their stability over the time. Though the average RMSD values of 15K and inhibitor 2 are higher compared to the apo form, they still indicate structural stability by other parameters, such as by the Rgyration and SASA (Fig. 6B,C). RMSF data demonstrate that 15K-PAK1 complex residues show lower fluctuation than other complexes. Importantly, the higher number of hydrogen bonds formation in 15K-PAK1 complex during MD simulation indicates the higher stability of 15K-PAK1 complex over inhibitor 2-PAK1 complex and apo form. The hotspot allosteric residues-Glu315, Met344, and Asp407-show less fluctuation in 15K-PAK1 complex than what they show in case of two other systems. PCA analysis categorizes 15K-PAK1 complex into a distinct group, clearly indicating that 15K-PAK1 complex bears a substantial structural and energy variation favoring its strong interaction with PAK1. MD simulation results firmly support the stable molecular interactions of 15K with PAK1 and therefore it can be considered as a potent allosteric ligand for PAK1 inhibition. SwissADME prediction reveals that 15K follows all the druglikeness rules according to Lipinski et al. 35 , Ghose et al. 36 , Veber et al. 37 , Egan et al. 38 , and Muegge et al. 39 , and the main features behind its druglikeness are low molecular weight (< 500), less than 5 hydrogen bond donors, and less than 10 hydrogen bond acceptors. Likewise, it would not cause any mutation, toxicity to the heart and skin as the prediction reveals. In vitro toxicity evaluation confirms it to be not toxic to the brine shrimp and fibroblast cells as well. At 50 nM concentration, a 2 × higher concentration than the IC 50 value against A549 lung cancer cells, it doesn't cause the death of C. elegans rather extending the lifespan of worms 12 . All these experiments together reinforce that 15K could be used as a safe drug, while ketorolac poses hepatotoxicity and renal toxicity 45 . Toxicity of 15K is therefore more specific towards cancer cells through inhibiting PAK1. Network and enrichment analysis also direct that it can breakdown the PAK1-driven protein network and in turn induces potent anticancer effects.
Based on the significant anticancer and anti-PAK1 performance of 15K, we designed and synthesized 1,2,3-triazolyl amide derivative of ketorolac and also synthesized an array of 1,2,3-triazolyl ester derivatives using other NSAIDs having COOH moiety including aspirin, naproxen, ibuprofen, mefenamic acid and indomethacin. But   www.nature.com/scientificreports/  Table 8. Pharmacological properties of 15K obtained through in silico PASS prediction, and the properties on which 15K showed more than 70% probability to be active are presented here. HIF1A Hypoxia-inducible factor 1-alpha, Pa Probability to be active, Pi Probability to be inactive.  www.nature.com/scientificreports/ unfortunately, all those derivatives did not show any anticancer effects on A549 lung cancer cells and B16F10 melanoma cells even at high concentration (200 nM) when tested by both trypan blue and MTT assays. Thus, they were not further assessed for anti-PAK1 activities. These results support that the ester bond of 15K has a significant role in its binding with PAK1 and it also has some novel features to be a potent anti-cancer drug. We here propose an SAR of 15K in Fig. 13. 1,2,3-Triazol moiety and the ester bond play a role in identifying and binding with the allosteric pockets, while the terminal ketorolac and methoxy benzene moieties are responsible for establishing hydrophobic bonds with PAK1. Although this study is limited for not showing any experimental evidence in support of strong binding characteristics of 15K, it warrants co-crystallization studies of 15K with PAK1 to confirm its binding to an allosteric  www.nature.com/scientificreports/ site and to reveal novel interactions. Future studies are also urgent to categorize 15K whether it is type III or type IV allosteric inhibitor.
In conclusion, our computational analysis clearly depicts that 15K can be a highly reactive and safe small molecule drug having all druglikeness features and it can strongly bind with the allosteric site of PAK1. It could act as a selective allosteric PAK1 inhibitor without showing any toxic side effects and can therefore be incorporated in clinical trials for anticancer drug development. Our findings could be used as a guide to study the crystallization of 15K-PAK1 complex to explore the novel interactions of 15K with different PAK1 domains. Besides, this study can be an important guideline in organic synthesis pipeline to develop more potent anticancer drugs, particularly the information that the ester bond and triazole moiety should be kept intact.