Study of Imidazolium Salt Derivatives as PIK3CA Inhibitors Using a Comprehensive in Silico Method

A series of imidazolium salt derivatives have demonstrated potent antitumor activity in prior research. A comprehensive in silicon method was carried out to identify the putative protein target and detailed structure-activity relationship of the compounds. The Topomer CoMFA and CoMSIA techniques were implemented during the investigation to obtain the relationship between the properties of the substituent group and the contour map of around 77 compounds; the Topomer CoMFA and CoMSIA models were reliable with the statistical data. The protein–protein interaction network was constructed by combining the Pharmmapper platform and STRING database. After generating the sub-network, the phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit α (PIK3CA with protein data bank ID: 3ZIM) was selected as the putative target of imidazolium salt derivatives. A docking study was carried out to correlate interactions of amino acids in protein active pockets surrounded by the ligand with contour maps generated by the structure-activity relationship method. Then the molecular dynamics simulations demonstrated that the imidazolium salt derivatives have potent binding capacity and stability to receptor 3ZIM, and the two ligand-receptor complex was stable in the last 2 ns. Finally, the ligand-based structure-activity relationship and receptor-based docking were combined together to identify the structural requirement of the imidazolium salt derivatives, which will be used to design and synthesize the novel PIK3CA inhibitors.


Introduction
Cancer is the second leading cause of death, globally, and was responsible for 8.8 million deaths in 2015. Globally, nearly 1 in 6 deaths are due to cancer [1]. Chemotherapy is generally the main treatment for various cancers. Chemotherapeutic agents (anti-cancer drugs) have a range of side-effects such as immunosuppression, myelosuppression, anemia, teratogenicity, infertility, and even secondary neoplasm [2][3][4][5][6]. The major goal of oncology scientists is to design a selective and effective anticancer agent that is only sensitive in normal cancer cells, as well as the ability to predict, alter, or block the hallmark of cancer cells and is likely to improve the therapeutic index [7]. Therefore, the search for a targeted, effective drug with minimum toxicity is urgently necessary [3,7].
Previous studies have shown that the phosphatidylinositol-3-kinase (PI3K) signaling pathway is a crucial one for many aspects of cell growth and survival. Abnormalities in the PI3K pathway are common in cancer and have a role to play in neoplastic transformation [8]. The most frequent genetic aberrations in cancer are linked to somatic missense mutations in the gene encoding PIK3CA (p110α) [9]. Given the important role of the PI3K signaling pathway, some selective inhibitors-PX-866 and PEG Wortmannin-have entered into preclinical status [10,11]. The correlation coefficient (R 2 ) between the Experimental pIC50 and predicted pIC50 of all 77 compounds is shown in Figure 2; the correlation coefficients of the topomer CoMFA and CoMSIA models were found to be 0.9027 and 0.9204, respectively, showing that the topomer CoMFA and CoMSIA models were reliable and precise for the prediction of activity.  The correlation coefficient (R 2 ) between the Experimental pIC 50 and predicted pIC 50 of all 77 compounds is shown in Figure 2; the correlation coefficients of the topomer CoMFA and CoMSIA models were found to be 0.9027 and 0.9204, respectively, showing that the topomer CoMFA and CoMSIA models were reliable and precise for the prediction of activity.  The correlation coefficient (R 2 ) between the Experimental pIC50 and predicted pIC50 of all 77 compounds is shown in Figure 2; the correlation coefficients of the topomer CoMFA and CoMSIA models were found to be 0.9027 and 0.9204, respectively, showing that the topomer CoMFA and CoMSIA models were reliable and precise for the prediction of activity.  The steric and electrostatic contour maps of topomer CoMFA are shown in Figure 3a-d, and the hydrophobic and hydrogen-bond acceptor contour maps of CoMSIA are shown in Figure 3e,f. The structure of the most active compound-compound 72-was selected as the reference structure for the generation and visualization of the topomer CoMFA and CoMSIA contour maps.  The steric and electrostatic contour maps of topomer CoMFA are shown in Figure 3a-d, and the hydrophobic and hydrogen-bond acceptor contour maps of CoMSIA are shown in Figure 3e,f. The structure of the most active compound-compound 72-was selected as the reference structure for the generation and visualization of the topomer CoMFA and CoMSIA contour maps. CoMFA steric contour maps of fragment 1 and 2, the green and yellow regions indicate the sterically favorable and unfavorable regions, respectively; (c,d) Topomer CoMFA electrostatic contour map of fragment 1 and 2, the blue and red regions are favorable to positively and negatively charged substituents, respectively; (e) CoMSIA hydrophobic contour map, the yellow and white regions are favorable and unfavorable to hydrophobic substituent groups, respectively; (f) CoMFA hydrogen bond acceptor field, magenta and red indicate regions favorable and unfavorable to hydrogen bond acceptor (HBA) atoms, respectively. The steric contour maps of the results of topomer CoMFA for fragment 1 and fragment 2 are shown in Figure 3a,b, respectively. In Figure 3a, the green contour map around the R1 substituent of imidazole/triazole ring indicates that this region was favorable for bulky groups. The results can be proved by the fact that compounds 72, 73 and 74 (with IC 50 values of 0.45, 0.68, and 0.58 µM, respectively) with 5,6-Dimethyl-benzimidazole exhibit more potent cytotoxic activities than compound The configuration change of the substituent group in R2 occurred when the hydroxyl in this place (in the upper right corner of Figure 3a) and the yellow region neared the substituent groups directly connected to the hydroxyl, which revealed that this place was sterically unfavorable for functional groups; for example, compounds 40 and 51, with 4-Bromobenzyl (with IC 50 values of 1.73 and 1.09 µM, respectively) in the R2 substituent group, were more positively charged, resulting in more potent activity than compounds 36 and 53 with 4-Bromophenacyl (with IC 50 values of 7.31 and 1.9 µM, respectively). Figure 3b shows the steric contour map of fragment 2, which is composed of two entirely different types of tricyclic substituent groups: flexible hexahydropyrrolo [2,3-b]indole substituents and rigid 9H-fluorene substituents. The green and yellow regions above the tricyclic substituent groups neared the 9H-fluorene and hexahydropyrrolo [2,3-b]indole, respectively; the combined region indicates that the compounds with 9H-fluorene substituents ( The topomer CoMFA electrostatic contour map of fragment 1 and 2 are shown in Figure 3c,d, respectively, in Figure 3c, the blue regions are found around the R1 and R3 regions, demonstrating that these positions were favorable for the electropositive group. Compounds in which the electropositive 5,6-Dimethyl-benzimidazole and 2-Methyl-imidazole moieties were present (compounds 4 and 72, with IC 50 values of 0.47 and 0.45 µM, respectively) will demonstrate higher cytotoxic activity than compounds with imidazole, triazole and benzimidazole (compounds 46, 57 and 30, with IC 50 values 3.49, 2.80 and 8.29 µM, respectively). The blue region near and around the benzene ring revealed that, whether or not there was hydroxyl in R2, the electropositive group was favorable for the activity, hence the compounds with electropositive methoxy and naphthyl benzene ring groups substituent in the benzene ring demonstrated more potent cytotoxic activity than those with electronegative fluorine and bromine substituents. It can be proved that compounds 36 and 40, with 4-Bromophenacyl and 4-Bromobenzyl (IC 50 = 7.31 and 1.73 µM) exhibited weaker cytotoxic activity than the compounds 37, 38 and 42, with 4-methoxyphenacyl, 2-naphthylacyl and 2-Naphthylmethyl (with IC 50 values of 6.23, 1.6 and 1.35 µM, respectively). In addition, the more the electropositive group is replaced in the benzene ring, the more the activity will decrease. For example, compound 75 (with IC 50 values of 1.78 µM) with fluorine substituent in the benzene ring demonstrated weaker activity than compound 74 (with IC 50 values of IC 50 = 0.58 µM) with bromine substituent. There is no distinct red region in the electrostatic contour map of fragment 1. In Figure 3d, the blue region near the N-benzyl group indicates that this place is favorable for electropositive groups, and compounds with benzyl groups demonstrated more potent activity than compounds without substituents in the N-benzyl group, which is consistent with the steric contour map of fragment 2.
The CoMSIA hydrophobic contour is shown in Figure 3e, the yellow region indicates that R1 was favorable for hydrophobic groups, which can be validated by the fact that the compounds with bicyclic Benzimidazole skeletons exhibited more potent activity than the monocyclic imidazole skeletons. The yellow region can also be seen in the benzyl at the R2 and R3 positions, so the more hydrophobic group naphthyl in R2 position and methyl in R3 position demonstrated higher cytotoxic activity (compound 71, with an IC 50 value of 0.57 µM) than other compounds that did not have this moiety. The white region around hexahydropyrrolo [2,3-b]indole indicates that the introduction of this hydrophobic skeleton was unfavorable for the cytotoxic activity, which is in accordance with the steric contour map of fragment 2; substituent groups with weak hydrophobicity were more beneficial for the activity than hydrophobic groups.
The CoMSIA hydrogen bond acceptor contour map is shown in Figure 3f, the red regions near the two sites of imidazole indicate that the introduction of the hydrogen bond acceptor carbonyl in this position was unfavorable for cytotoxic activity, which corresponds with the results of the topomer CoMFA steric contour map. There was no distinct magenta region around the compound, because the structural modifications of the hydrogen bond acceptors of compounds did not exhibit good biological activity.
We received 3000 protein targets for 10 active compounds from the PharmMapper result list. After removing duplicates, 722 targets were used for screening cancer-related proteins; only proteins with clear cancer drug ligands can be used for further PPI analysis. Finally, 27 cancer-related proteins were identified, as shown in Table 2. Then all the proteins were uploaded to the STRING database to find their direct and functional partners and to obtain the primary PPI network of each protein. The software Cytoscape 3.5.0 (U.S. National Institute of General Medical Sciences, Bethesda, MD, USA) was utilized to merge the PPI network and analyze the merged network using its functional plugins. Finally, a network with 104 nodes and 496 edges was obtained, as shown in Figure 4a. skeletons. The yellow region can also be seen in the benzyl at the R2 and R3 positions, so the more hydrophobic group naphthyl in R2 position and methyl in R3 position demonstrated higher cytotoxic activity (compound 71, with an IC50 value of 0.57 μM) than other compounds that did not have this moiety. The white region around hexahydropyrrolo [2,3-b]indole indicates that the introduction of this hydrophobic skeleton was unfavorable for the cytotoxic activity, which is in accordance with the steric contour map of fragment 2; substituent groups with weak hydrophobicity were more beneficial for the activity than hydrophobic groups. The CoMSIA hydrogen bond acceptor contour map is shown in Figure 3f, the red regions near the two sites of imidazole indicate that the introduction of the hydrogen bond acceptor carbonyl in this position was unfavorable for cytotoxic activity, which corresponds with the results of the topomer CoMFA steric contour map. There was no distinct magenta region around the compound, because the structural modifications of the hydrogen bond acceptors of compounds did not exhibit good biological activity.
We received 3000 protein targets for 10 active compounds from the PharmMapper result list. After removing duplicates, 722 targets were used for screening cancer-related proteins; only proteins with clear cancer drug ligands can be used for further PPI analysis. Finally, 27 cancer-related proteins were identified, as shown in Table 2. Then all the proteins were uploaded to the STRING database to find their direct and functional partners and to obtain the primary PPI network of each protein. The software Cytoscape 3.5.0 (U.S. National Institute of General Medical Sciences, Bethesda, MD, USA) was utilized to merge the PPI network and analyze the merged network using its functional plugins. Finally, a network with 104 nodes and 496 edges was obtained, as shown in Figure 4a.
EGFR ERBB ERBB1 HER1 Figure 4. Network of cancer-related targets (a) and sub-network with essential targets (b), as generated by Cytoscape software and the CytoNCA plugin.  The CytoNCA plugin was utilized to calculate the Subgraph Centrality, Betweenness Centrality and Closeness Centrality of all 104 nodes. After calculating the eight centrality measures of all nodes (see Table S1), all nodes were sorted by three centrality measures in descending order, the top 10% of the three centrality measures are colored black, blue and white respectively, and the overlapping nodes are colored with color mixtures, with the colored network being shown in the Figure. The top 10% of the three centrality measures were merged together to generate the sub-network with essential nodes, as shown in Figure 5b; in the merged list (see Table 3), the nodes are sorted by the three comprehensive centrality measures in descending order. The PIK3CA had the highest values for the three centrality measures, which indicates that PIK3CA produces some interactions with other proteins, and serves as a more essential node than the other proteins obtained from PharmMapper, hence the PIK3CA (with PDBID: 3ZIM in the PharmMapper results) was selected for the further docking study.  Table S2). The ligand interaction of compound 04, 72 and 74 is exhibited in Figure 5b-d, analyzing the interaction between the substituents in the compound skeleton and the key amino acids in the active pocket. Figure 5a exhibits the combinations of selective targeted covalent inhibitor CNX-1351 and PIK3CA, this ligand-receptor interaction was selected as the reference for analyzing the docking results of the three compounds, because the ligand-binding affinity to the receptor has been verified by in vitro biological experiment [29].

PDB ID Gene Names
In Figure 5a, the inhibitor CNX-1351 formed several important interactions with amino acids Asp810, Val851, and Ile932. Similarly, in the three docking compounds, the Val851 formed arene-H interactions with the benzyl in N-benzyl group; this docking result was consistent with the steric, electrostatic and hydrophobic contour maps of topomer CoMFA and CoMSIA, which revealed that the benzyl in this position was beneficial for the improvement of the activity. In compounds 72 and 74, the Val 851 and Ile932 formed arene-H interactions with the 9H-fluorene moiety of the compound skeletons, while the hexahydropyrrolo[2,3-b]indole moiety of compound 04 did not form any interaction with the amino acid, indicating that compounds with rigid 9H-fluorene substituents were more similar to inhibitor CNX-1351 than compounds with flexible hexahydropyrrolo[2,3-b]indole rings; this also confirms the topomer CoMFA steric contour map and the CoMSIA hydrophobic contour map. In the docking results for compounds 04 and 72, the Val851 and Lys702 formed arene-H interactions with phenyl in the R2 and R1 positions, respectively, and as the result of CoMFA and CoMSIA, bulky, electropositive and hydrophobic groups were favorable for the cytotoxic activity. The imidazolium moieties of compounds 04, 72 and 74 formed ion contacts with acidic amino acids Asp933 and Glu849; the imidazolium served as the core moiety in a series of imidazolium salt derivatives. While performing structural modification of the imidazolium moieties of compounds 45, 46, 47, 49 and 50 (with IC 50 values of 2.17, 3.49, 1.75, 1.10 and 1.01 µM, respectively), we found that the remodeled compounds 56, 57, 58, 59 and 60, which possess the triazolium moiety (with IC 50 values of 2.05, 8.29, 2.07, 2.55 and 1.70 µM, respectively) did not exhibit any improvement in cytotoxic activity. Additionally, some greasy amino acids, basic amino acids and acidic amino acids around compounds 04, 72 and 74 were also similar to the inhibitor CNX-1351. the remodeled compounds 56, 57, 58, 59 and 60, which possess the triazolium moiety (with IC50 values of 2.05, 8.29, 2.07, 2.55 and 1.70 μM, respectively) did not exhibit any improvement in cytotoxic activity. Additionally, some greasy amino acids, basic amino acids and acidic amino acids around compounds 04, 72 and 74 were also similar to the inhibitor CNX-1351. The molecular dynamics simulation results for the 04-3zim complex and the 72-3zim complex at different temperatures is shown in Figure 6. The values of root mean square deviations (RMSD) indicate that the transformation of the ligand-protein complex backbone forms the initial structure. In results of the compound 04-3zim, the RMSD values increased to 2.2 Å in 2.3 ns, and then maintained at between 2.5-3 Å until 5 ns; similarly, in the 72-3zim complex, the RMSD increased to 2.5 Å in 2.3 ns and then maintained at between 2.6-3.2 Å until 5 ns. The mean RMSD values of the 04-3zim complex and the 72-3zim complex were 2.115 and 2.253, respectively. The MD results revealed that these two active compounds had potent binding capacity and stability to the putative receptor.  The molecular dynamics simulation results for the 04-3zim complex and the 72-3zim complex at different temperatures is shown in Figure 6. The values of root mean square deviations (RMSD) indicate that the transformation of the ligand-protein complex backbone forms the initial structure. In results of the compound 04-3zim, the RMSD values increased to 2.2 Å in 2.3 ns, and then maintained at between 2.5-3 Å until 5 ns; similarly, in the 72-3zim complex, the RMSD increased to 2.5 Å in 2.3 ns and then maintained at between 2.6-3.2 Å until 5 ns. The mean RMSD values of the 04-3zim complex and the 72-3zim complex were 2.115 and 2.253, respectively. The MD results revealed that these two active compounds had potent binding capacity and stability to the putative receptor. Finally, the structural requirement of a series of imidazolium salt derivatives was identified by a comprehensive method that combined ligand-based QSAR study and receptor-based docking-MD simulation (see Figure 7).  Finally, the structural requirement of a series of imidazolium salt derivatives was identified by a comprehensive method that combined ligand-based QSAR study and receptor-based docking-MD simulation (see Figure 7).  The molecular dynamics simulation results for the 04-3zim complex and the 72-3zim complex at different temperatures is shown in Figure 6. The values of root mean square deviations (RMSD) indicate that the transformation of the ligand-protein complex backbone forms the initial structure. In results of the compound 04-3zim, the RMSD values increased to 2.2 Å in 2.3 ns, and then maintained at between 2.5-3 Å until 5 ns; similarly, in the 72-3zim complex, the RMSD increased to 2.5 Å in 2.3 ns and then maintained at between 2.6-3.2 Å until 5 ns. The mean RMSD values of the 04-3zim complex and the 72-3zim complex were 2.115 and 2.253, respectively. The MD results revealed that these two active compounds had potent binding capacity and stability to the putative receptor. Finally, the structural requirement of a series of imidazolium salt derivatives was identified by a comprehensive method that combined ligand-based QSAR study and receptor-based docking-MD simulation (see Figure 7).

QSAR Study
77 imidazolium salt derivatives were collected from two literature sources and summarized in Table 4, including 63 imidazolium salt derivatives and triazolium salt derivatives. The IC 50 values of HL-60 tumor cell lines in the two literature sources were measured over a similar period of time and in the same laboratory [18,19], so we were able to merge these comparable IC 50 values for the two series of compounds with common nuclear structures for further 3D-QSAR analysis. The IC 50 values of all 77 compounds were converted to pIC 50 ; these pIC 50 values served as the dependent variable, while the descriptor values of CoMSIA and topomer CoMFA served as the dependent variables. The sketch function in software Sybyl X2.0 (Certara, Princeton, NJ, USA) was utilized to draw the structure of the compound, and these 77 molecules were divided in a ratio of 75:25 by Sybyl X2.0 in such a way that both datasets consisted of a balance of active and less-active molecules [30,31]. After minimizing the energy of the compound in the tripos force field with 0.05 kcal/(mol·A) termination and 1000 iterations, the compound was put into the database for topomer CoMFA and CoMSIA analysis in Sybyl X2.0 [32,33].

QSAR Study
77 imidazolium salt derivatives were collected from two literature sources and summarized in Table 4, including 63 imidazolium salt derivatives and triazolium salt derivatives. The IC50 values of HL-60 tumor cell lines in the two literature sources were measured over a similar period of time and in the same laboratory [18,19], so we were able to merge these comparable IC50 values for the two series of compounds with common nuclear structures for further 3D-QSAR analysis. The IC50 values of all 77 compounds were converted to pIC50; these pIC50 values served as the dependent variable, while the descriptor values of CoMSIA and topomer CoMFA served as the dependent variables. The sketch function in software Sybyl X2.0 (Certara, Princeton, NJ, USA) was utilized to draw the structure of the compound, and these 77 molecules were divided in a ratio of 75:25 by Sybyl X2.0 in such a way that both datasets consisted of a balance of active and less-active molecules [30,31]. After minimizing the energy of the compound in the tripos force field with 0.05 kcal/(mol·A) termination and 1000 iterations, the compound was put into the database for topomer CoMFA and CoMSIA analysis in Sybyl X2.0 [32,33].  Topomer CoMFA included steric and electrostatic fields, while CoMSIA included steric, electrostatic, hydrophobic, hydrogen bond donor (HBD) atom, and hydrogen bond acceptor (HBA) atom fields. PLS (partial least squares) techniques associated these field descriptors with the activity value [30]. Many statistics, such as the values Leave One Out (LOO), optimal number of components (ONC), Standard Error of Estimation (SEE), mean absolute error (MAE), cross-validated coefficients (q 2 ), and conventional coefficient (r 2 ), were important in the evaluation of the 3D-QSAR model, and could be worked out using the PLS method. The QSAR model is said to be good when the q 2 value is greater than 0.5 and the r 2 value is greater than 0.6; because q 2 and r 2 values reflect the model soundness, the best QSAR model will have the highest q 2 and r 2 value, the lowest SEE, and optimal number of components [31]. In the case of topomer CoMFA analysis, the PLS leave-one-out (LOO) method with CoMFA standard options for variable scaling was implemented to investigate the topomer CoMFA model [32].
To build the topomer CoMFA QSAR model, the topomer technique was applied to split the molecules into two fragments (the cutting method is shown in Figure 4a) [32]. Database alignment was used to build the CoMSIA QSAR model [31,33]; the imidazole ring was identified as the common core moiety, and the most active compound-72-was selected as the template.
In the CoMSIA analysis, the values of ONC, SEE, and q 2 were worked out by LOO validation, turning on USE SAMPLES, and components set to ONC. In the process of calculating r 2 , USE SAMPLES was turned off, and column filtering was set to 2.0 kcal·mol −1 to speed up the calculation without sacrificing information content [30][31][32][33], components were set to ONC, which was the optimal number of components calculated by performing a SAMPLES run. SEE and r 2 were utilized to assess the non-cross-validated models.

PPI Network Construction and Analysis
PharmMapper serves as a valuable tool for identifying potential targets for novel synthetic compounds, newly isolated natural products, compounds with known biological activity, or existing drugs whose mechanism of action is unknown [23,34].
The top 10 compounds with the highest pIC 50 values of all imidazolium salt derivatives were selected for subsequent study using PharmMapper, and the mol2 formats of these compounds were uploaded to the PharmMapper server, the parameters Generate Confomers and Maximum Generated Conformations were set to ON and 300, respectively, the target set Druggable Pharmacophore Models (v2017, 16159) served as pharmacophore mapping, and all other options conformed to the default settings [23,34].
Protein-protein interactions (PPIs) participate in many metabolic processes that occur in living organisms, such as cellular communication, immunological response, and gene expression control [24]. Systematic description of these interactions aids elucidation of the interrelationships among targets. Thus, targeting PPIs with small-molecule compounds is becoming an essential step to comprehending even major-target mechanisms [24]. Cancer-related targets of imidazolium salt derivatives were gathered by screening the results of PharmMapper, these proteins were used as the input protein in the STRING (https://string-db.org/) database [35] to search for related proteins or pathways that had previously been reported. Finally, we merged all the PPIs of cancer-related targets of imidazolium salt derivatives into an integrated network to identify the critical protein in Cytoscape 3.5.0 software.
Cytoscape is a free open-source Java application for visualizing molecular networks and integrating gene-expression profiles [36]. Plugins are available for network and molecular profiling analyses, new layouts, additional file format support and connection with databases and searching in large networks [36,37].
CytoNCA is a plugin in Cytoscape that integrates calculation, evaluation and visualization analysis for multiple centrality measures. There are 8 centrality measures provided by CytoNCA: Betweenness, Closeness Centrality, Degree Centrality, Eigenvector Centrality, Local Average Connectivity-based Centrality, Network Centrality, Subgraph Centrality, Information Centrality [38]. The centrality measurement analysis was conducted in order to identify the essential proteins in an already-built PPI network.

Molecular Docking and Dynamics
MOE (Molecular Operating Environment) is a comprehensive Computer-Aided Drug Design (CADD) software program that incorporates the functions of QSAR, molecular docking, molecular dynamics, ADME (absorption, distribution, metabolism, and excretion) and homologous modeling; all these functions can be regarded as conducive instruments in the field of drug discovery and biochemistry. The molecular docking and dynamics technology were performed in MOE2016 software to detect the stability and affinity between the ligands and predictive targets [39,40].
In the case of the docking study, the pdb formats of essential proteins were downloaded from the RCSB Protein Data Bank (http://www.rcsb.org/pdb/home/home.do), which served as the receptors, while the mol2 formats for the top 10 compounds with highest pIC 50 were obtained from Sybyl software and served as the ligands. All the ligands will be docked at the active site of receptors in the manner of Triangle Matcher placement and Induced Fit refinement; the five docking conformations are listed in descending order by London dG functional score in the results panel. Based on the London dG functional score, the interaction between the structure of the compounds and the key amino acids in the active pocket can be visualized using the Ligand interaction function [40].
Integrated MD simulations for the ligand-receptor complex were performed using the program dynamics simulations, and all complexes were generated from the docking results. MD simulations were implemented in three phases: preparation, restrictive MD simulations, and non-restrictive MD simulations. During the preparation phase, the minimization and equilibration of constrained complexes were performed in a 15 Å3-sized water box [41], and an Amber12 EHT force field file was applied for energy minimization and equilibration with Gasteiger-Huckel charges using Boltzmann initial velocity. Restrictive MD simulations were carried out after the preparation phase, while gradually raising the system temperature from 300 to 310 K, and maintaining pressure at 101.325 kPa within 100 ps [42,43]. Finally, a 5 ns MD non-restrictive MD simulation with a 2 fs timestep, was carried out to investigate the stability and affinity of the ligand-receptor complex by contrasting the root mean square deviations (RMSD) [41,43].

Conclusions
A novel comprehensive method was carried out to identify the target and structural requirements of imidazolium salt derivatives. Cancer-related protein PIK3CA was identified as the putative target for a series of imidazolium salt derivatives by the PharmMapper platform and PPI network. In the following QSAR and docking study, the important substituent groups related to the cytotoxic activity formed molecular interaction directly with the key amino acids around the active pocket; further MD simulations showed the docking conformation to have potent binding affinity and stability. After correlating these ligand-based and receptor-based results, we found that bulky, electropositive and hydrophobic groups were beneficial for cytotoxic activity, and that compounds with imidazolium and 9H-fluorene moieties exhibited more potent activity than those with hexahydropyrrolo[2,3-b]indole and imidazolium moieties; additionally, the HBA and electronegative atoms on the two sides of the imidazole ring were adverse to the activity. This comprehensive method not only serves as a guideline for designing and synthesizing a novel imidazolium salt derivative PIK3CA inhibitor, but also provides new ideas for the structural modification of known compounds.