Searching for Natural Aurora a Kinase Inhibitors from Peppers Using Molecular Docking and Molecular Dynamics

1. Introduction

The Aurora kinase family, including aurora kinases A, B, and C, is a group of highly conserved serine/threonine kinases that are important for the correct transition through mitosis [1]. Errors in mitosis can lead to genomic changes that are typically associated with tumorigenesis [2]. Although the faithful process of mitosis is connected with the activity of Aurora kinases, its overexpression may be connected with the development and course of some cancers, such as laryngeal, ovarian, breast, colorectal and gastric cancer [3]. This is associated with the partial degradation of the p53 protein, the physiological function of which leads to damage in the apoptotic pathway. Aurora kinase, due to the phosphorylation of p53 (Ser315) protein, inhibits this activity. Moreover, as described earlier, the overexpression of Aurora kinases induces chemoresistance in breast cancer cells [4]. The overexpression of Aurora A may result in resistance to the treatment of breast cancer with paclitaxel via the transactivation of the FOXM1 promoter [5]. Since the therapeutic process requires the inhibition of overexpressed Aurora kinase, the use of several small molecule aurora kinase inhibitors as potential anticancer therapeutic agents has recently been investigated. Although there are numerous studies of synthetic molecules [6,7,8], natural products or directly derived compounds play a crucial role in the discovery of new drugs.

One of many well-known and studied active molecules derived from natural products is capsaicin from chili peppers (Capsicum). The traditional use of peppers and their use in modern societies has raised interest in their biological applications, including cytotoxic and antiproliferative effects. Hui-Chung et al. reported that Aurora A protein increased in response to cisplatin and was degraded upon combined treatment with capsaicin with cisplatin, suggesting that the Aurora A-mediated signaling pathway is responsible for the resistance to cisplatin in cisplatin-resistant gastric cancer cell lines [9].

Therefore, the aim of this study was to find out if the other compounds contained in the Piper, Capsicum and Pimenta (the main genera of peppers consumed worldwide), as well as their analogs, can act as Aurora A (Aurka) inhibitors. Since capsaicin is one of the best-known activators of the ion-channel type TRPV1 receptor, we looked for the connection between TRPV1 and Aurora inhibiting activity.

A set of ligands was used composed of 16 compounds occurring in Piper, Capsicum and Pimenta genera or their derivatives (see . The compounds were selected according to their reported antitumor activity against breast cancer. Molecular docking calculations were carried out to Aurka and TRPV1.

This is the first step in the study of natural compounds such as Aurka inhibitors. Currently, in silico studies often precede the in vivo and in vitro tests and show the direction in which to conduct further biological research [10]. Molecular docking has become one of the increasingly important tools for drug discovery [11]. This method is less time- and cost-consuming and allows one to choose the best promising compounds for in vitro tests. Molecular docking is widely used in studies of anticancer properties of small molecules. It has been used to describe the disruption of binding between Bcl-2 and Bax proteins involved in the induction of apoptosis [12]. The additional use of molecular dynamics calculations enables the description of the changes in the structural integrity of proteins [13] or to predict the structure of the protein [14]. Our results will be used to plan further molecular dynamics and in vitro studies.

2. Results and Discussion

Firstly, the validation of the molecular docking protocol has been performed. A set of 12 ligands with known affinity to Aurka has been created. Kd values showing affinity to Aurora A protein were taken from PubChem (AID 624919). The structures of ligands, Kd values, and results of docking are included in . Docking scoring functions were compared with the literature Kd values. Then, the linear regression was analyzed, and coefficients of determination (R²) were calculated. The results are presented in . The highest R² value has been observed using Surlex (0.67). AutoDock and Glide gave less significant correlations (0.41 and 0.21, respectively). Therefore, we have chosen Surflex as the software that gives the best results.

In the second step of validation, docking of the Aurka inhibitor present in the 3H10 crystal structure was performed, and its conformation was compared with that of the crystal structure . Both conformations are similar to each other, which is an additional confirmation of the correctness of the simulations.

Finally, we checked the chosen program (Surflex) on ligands with lower affinity to Aurka (Kd > 10 according to PubChem AID 624919). The results show that Total Score values for these ligands are lower than 5.00. It indicates that this model may be used also for compounds with lower affinity to Aurka.

After validation, molecular modeling simulations of a set of 16 ligands to the Aurka protein were performed. In the first stage, docking was carried out using three different programs (Surflex, AutoDock, and Glide). However, the optimization of the program used for docking has been performed earlier, so we wanted to compare results from these three programs. It would be interesting if the comparison of docking scores with biological activity measured as IC50 towards MCF-7 cells and affinity measured as Kd would be comparable. shows a list of ligands, the IC50 values determined on MCF-7 cell lines (breast cancer cells), log p values, and the molar mass of each ligand.

The results obtained from docking are presented in and .

To better understand the contribution of various energy terms to the binding of ligands to Aurka, the MM/GBSA calculations have been performed . It has been found that for the ligands, the binding results mostly from van der Waals interactions and non-polar (lipophilic) contributions to the solvation free energies. Positive values of the polar contributions to the solvation free energies found for all ligands (except 14 and 15) indicates that the desolvation of the ligand is thermodynamically non-favorable as all those compounds are well soluble in water due to the multiple interactions they form with this solvent, which is a direct consequence of multiple hydrogen bond acceptors present in those molecules.

To compare the results obtained from docking using Sybyl, we performed docking using AutoDock. The results of docking are summarized in . The docking binding energy, intermolecular energy, total internal energy, torsional energy and unbound energy allowed to distinguish the real binding conformations of ligands on Aurora kinase from energetically unfavorable ones like in the study of Reghav et al. [23].

The results from AutoDock suggest that all ligands have a similar affinity to Aurka. The overall comparison of binding energy values showed the highest binding energy (−8.97 kcal/mol) for ligand 1. Ligands 1 and 10 have calculated IC in the range of nM, which corresponds to the lowest binding energy. Additionally, ligands 15 and 16 have the lowest affinity to Aurka using Surflex and AutoDock.

The correlation was statistically significant (Pearson correlation, p < 0.05). The results obtained from docking are expressed by the Total Score value related to the stability of the ligand–protein complex (dissociation constant). The obtained values make it possible to rank the compounds according to the supposedly increasing affinity for Aurka.

Then, we compared the Total Score values obtained for the tested ligands and the Total Score values calculated in the same way for Aurka inhibitors with known biological activity. Three Aurka inhibitors were selected: AZD1152, hesperidin and PHA-73935. The results obtained from docking, along with the literature IC50 value (Aurka), are summarized in .

Total Score values obtained for known Aurka inhibitors are in the range of 5.817–8.043 This means that by comparing only the values of the scaling functions, it can be assumed that ligands 1–8 with Total Score values within this range can have the properties inhibiting the activity of Aurora A kinase. As the analysis of the Total Score function alone is insufficient, the interactions exhibited by the ligands in the Aurka active site should also be analyzed.

To verify the correctness of the molecular docking performed in Surflex, we analyzed the arrangement and interactions of hesperadin in the Aurka active site and then compared them with the results published earlier. Talele and McLaughlin published a study on the molecular docking and molecular dynamics of Aurka inhibitors [28].

According to their study, hesperadin forms hydrogen bonds with Ala 213 (O···HN) and Glu211 (-NH···O) backbone atoms. Moreover, hydrophobic contacts with Thr217 and Leu263 are necessary to stabilize the ligand–protein complex. The sulfonamide group is located next to Lys162 within the hydrogen bonding distance. Talele and McLaughlin confirmed the formation of hydrogen bonds using molecular dynamics simulations [28].

Our results are in agreement with the above results. We obtained the conformation of hesperadin stabilized using Ala213 and Lys162 hydrogen bond interaction (2.089 Å and 1.976 Å, respectively) and localized in the active pocket made by the hydrophobic group of Thr217 and Leu263 (Figure 1).

Figure 1. Conformation of hesperadin in the active site of Aurka from molecular docking studies.

Since in earlier stages, we found that the results obtained with Surflex correlate best with IC50 and Kd values, the tested compounds have Total Score values similar to known Aurora A kinase inhibitors, and the configuration and interactions of hesperadin obtained in our simulations correspond to those published previously, we also analyzed the interactions in the active site of ligands 1–3, according to their best Totals Score values . All three ligands were docked in the pocket formed by Arg137, Lys162, Pro214, Ala213, Tyr212, Asp274, Thr 217, Phe275 and Leu263. It was the same pocket in which hesperadin had been docked. All ligands 1–3 form hydrogen bonds with Lys162 by Me-O···NH interactions (1—1.847 Å, 2—2.625 Å, 3—1.980 Å). Ligand 2 is also stabilized in the active site by forming a hydrogen bond with Tyr212 (1.906 Å) and ligand 3 with Asp274 (2.346 Å) and Phe 275 2.677 Å). The ligand–protein interactions are presented in Figure 2.

Figure 2. Conformations of (1) (red), (2) (blue) and (3) (magenta) in the active site of Aurka from molecular docking studies.

The results of the in silico study (molecular docking) suggest that a large group of compounds tested by us can potentially act as Aurka inhibitors. The analysis of the docking site shows the importance of the ligand hydrogen bond formation with the amino group of Lys162. This bond is usually formed by oxygen of the -OMe group attached to the aromatic ring.

Surprisingly, this structure is important for TRPV1 ligands. The structures responsible for the formation of hydrogen bonds are present in all TRPV1 ligands. The general scheme of the TRPV1 ligand structure is shown in Figure 3.

Figure 3. A scheme of TRPV1 receptor ligand.

In the central part, urea (X = N), amide (X = C), or thiourea moiety may be present. The aromatic group can be monocyclic or bicyclic. The π–π interactions between the ring and the corresponding aromatic fragment in the receptor affect the additional stabilization of the molecule in the active site. Such interactions are observed, for instance, for capsaicin, which can form a specific conformation with the bent side chain in the active site of the receptor [29]. Other capsaicinoids can exist in a similar conformation in the active site of TRPV1 [30]. The lipophilic side chain interacts with amino acids of the hydrophobic side chain present in the TRPV1 binding site. The connectors act as a “scaffold” that links both elements together [31].

Since the same moieties are involved in hydrogen bond formation in the active site of Aurka and TRPV1, we thought that there could be a relationship between Aurora– and TRPV1–ligand interactions.

The results of the molecular docking to TRPV1 receptor expressed as Total Score function are presented in in . As we expected, the best affinity was observed for capsaicin. The conformation of capsaicin was similar to those described by Darre and Domene [32]. The interactions include hydrogen bonding to Thr550 and Tyr511 and stabilization in the hydrophobic cluster formed by Leu515, Ile573, Phe587 (aromatic ring) and Met547, Ala665, Phe591 (aliphatic chain; Figure 4).

Figure 4. Conformation of capsaicin (red) in the active site of TRPV1 (5IS0).

To observe the relationship between the measured variables, PCA (Principal Component Analysis) was applied. Data include not only Total Score functions from molecular docking to Aurka and TRPV1 (only Surflex was included because it showed the best correlation with IC50 values for MCF-7) but also logP, molar mass and IC50 (MCF-7 cells) values.

Figure 5 shows a standard PCA biplot where the first principal component (PC1) describes 26% and the second (PC2) describes 41% of the total subset variance.

Figure 5. A standard biplot PC1 vs. PC2. Aurka, TRPV1—Total Score values from molecular docking; M—molar mass.

These results show an inverse correlation between the affinity to the Aurka and IC50 values. Moreover, the IC50 value can also be connected with the logP value and binding energy to Aurka but not with the affinity to TRPV1 nor molar mass.

The points corresponding to the molecules with the best activity against both MCF-7 and Aurka are located on the right or central side of the plot. The molecules with the lowest activity are on the left side, which corresponds to the highest IC50 and lowest Total Score, but also higher logP values. This suggests that both the inhibition of Aurora A kinase and the anti-tumor activity against MCF-7 cells require lipophilic properties (preferably a logP value below 3). This corresponds to one of Lipinski’s rule of five, stating that an orally active drug should have an octanol–water partition coefficient [33] (logP) that does not exceed 5. In our study, particles with logP values above 4 have the lowest Total Score values for Aurka and the highest IC50 tested on MCF-7.

The analysis of docked ligands shows that they meet the requirements of the Lipinski rule of 5 [34]. These molecules have no more than five hydrogen bond donors, no more than 10 hydrogen bond acceptors, a molecular mass of less than 500 daltons and an octanol–water partition coefficient (logP) that does not exceed 5. However, in compounds with the lowest docking score, the logP value almost exceeds 5, and ligand 7 has a molecular mass of more than 500 daltons.

It is possible that the inhibition of the overexpression of Aurka is one of the mechanisms of anticancer activity observed in MCF-7 cell line studies. Chen et al. studied the anticancer activity of capsaicin against breast cancer and concluded that the mechanism could be related to the down-regulation of the FBI-1-mediated NF-κB pathway [35]. The anticancer activity of capsaicin can also be related to the caspase-independent pathway [36]. Thus, the inhibition of Aurka can be one of the new possible anticancer mechanisms of the compounds obtained from peppers. To confirm these results, the in vitro studies are in progress.

While the ligand–receptor complex obtained from a molecular docking method captures just one moment of the mutual orientation of a ligand and target protein, a more computationally demanding but also much more informative method is molecular dynamics (MD) simulations. This approach enables simulation of the behavior of the modeled system (e.g., ligand–protein complex) in time, allowing to obtain more ligand–protein poses. Consequently, the more poses are being considered, the higher the amount of information provided. Therefore, in this study, we applied both molecular docking and MD simulations to explain the activity profiles of chosen Aurka ligands.

Figure 6 and Figure 7 show the root mean square deviation of atomic positions (RMSD) plots obtained for Cα-atoms of Aurka (Figure 6) and its chosen ligands (Figure 7) in the complexes formed with protein during 100 ns MD simulation. RMSD analysis can indicate if the simulation has equilibrated—its fluctuations towards the end of the simulation are around some thermal average structure. Usually, changes of the order of 1–3 Å are perfectly acceptable for small proteins. From Figure 6, it can be seen that the RMSD has stabilized itself at the level of c.a. 2.5 Å after c.a. 5 ns of simulations, meaning that the protein molecule in all of the studied complexes is stable from the structural point of view.

Figure 6. RMSD plot obtained for Cα-atoms of Aurka in the complexes formed between Aurka and chosen ligands during 100 ns MD simulation.

Figure 7. RMSD plots obtained for the studied ligands in the complexes formed between them and Aurka during 100 ns MD simulation.

The analysis presented in Figure 7, however, shows totally different behavior of 1 in comparison with other ligands. For 2, 3 and AZD1152, the RMSD reached the plateau after a few ns of simulation, at the level of 0.5 Å for 2 and 2.0 Å for 3 and AZD1152. However, the complex between the 1 and Aurka, while initially stable during the first 10 ns of simulation, started to dissociate between 15 and 35 ns and dissociated completely during the second half of the simulation period. This is indicated by the large RMSD value, exceeding 15 Å at some points. Additionally, the initial and final snapshots of the complex confirm the observation about its instability and dissociation (Figure 8).

Figure 8. Initial (A), and final (B), snapshot of the 100 ns MD simulation of the complex formed between 1 and Aurka. During the simulation the dissociation of the complex occurred.

This observation suggests that to better understand the behavior of the complexes formed between Aurka and its ligands, at least 50 ns molecular dynamics simulations should be performed. It is especially important as in the previously published works describing MD simulations of complexes formed between various other ligands and Aurka, significantly shorter total simulation times have been used, i.e., 1 ns [28] and 20 ns [37].

Root mean square fluctuation (RMSF) is a numerical measurement similar to RMSD, but instead of indicating positional differences between entire structures over time, RMSF is a calculation of individual residue flexibility or how much a particular residue moves (fluctuates) during a simulation. shows the RMSF plots obtained for Cα-atoms of Aurka in the complexes formed with the studied ligands during 100 ns MD simulation. It can be observed that the tails (N- and C-terminal) fluctuate more than any other part of the protein. The secondary structure elements like alpha helices and beta strands are more rigid than the unstructured part of the protein and thus fluctuate less than the loop regions. In all of the studied complexes, the maximal RMSF per residue was larger than 0.5 Å and was usually not exceeding 3.0 Å, with the exceptions for the termini, both N- and C-. This indicates that no abnormally high dynamic residues were detected, indicating the correctness of the studied model of macromolecules. The major differences in the behavior of Aurka in the studied complexes include the residues 160–170. This fragment is also characterized by the highest RMSD values. Also, solely for the complex with 3, the large RMSD values were recorded for residues 125–130.

To assess the repeatability of results, the molecular dynamics simulations of the complex formed between the Aurka and AZD1152 using the same computational settings have been repeated three times, and the comparisons of the RMSD obtained from those three trajectories are presented in Figure 9 and Figure 10. No major changes have been observed, both for the ligand and the protein.

Figure 9. RMSD plot obtained for Cα-atoms of Aurka in the complexes formed between Aurka and AZD1152 during each of the three molecular dynamics trajectories (R1, R2 and R3) (A) and average value (red) with the error bars (black) (B).

show the secondary structure (SS) content of the Aurka in its complexes with the studied ligands during the MD simulation. It can be seen that for the vast majority of the residues, their alignment to the particular type of SS remains constant during the simulation, with the exception of the regions characterized by large values of RMSF, which are mostly assigned to unstructured regions of the receptor. In both cases, the SS content remains constant during the simulation, indicating that no major structural transitions can be observed. From the SS analysis point of view, the dynamic behavior of Aurka in all of the studied complexes was similar, with a constant value of SSE at 45% during the whole simulation time.

Specific protein–ligand contacts arise due to multiple weak, low-energy (1–5 kcal/mol), noncovalent interactions such as H-bonds, ionic, and hydrophobic forces at short distance ranges enough for bonding, typically 2.5–3.5 Å.

Protein–ligand interaction diagrams are presented in Figure 11, . It should be noted that interactions observed during the MD simulations differ from those obtained from molecular docking. This is due to the presence of water molecules in the system subjected to MD simulations, which were absent in molecular docking. In addition, the simulation enabled the observation of the conformational changes in the ligand and the formation of intermolecular forces between the different functional groups of ligands, mostly in the case of the complex with AZD1152.

Among the residues that form interactions with all of the studies ligands are Arg137, Leu139, Val147, Ala160, Lys162, Leu194, Leu210, Ty212, Ala213, Pro214, Thr217, Arg220, Leu263 and Val279. Only some of those residues, i.e., Val279, form similar interactions in terms of both type and strength, regardless of the ligand. For most of the residues, significant differences can be observed in terms of their interactions with the studied inhibitors. For example, Lys162, which is very important in terms of binding with Aurka, form either direct hydrogen bond type interactions (with 2) or mixed direct/water bridge mediated H-bonds (with 1 and 3). For AZD1152, the interaction with Leu162 is partially hydrogen bonding via a water bridge and partially hydrophobic. The importance of this introduction is also shown in the case of the complex of Aurka with 1. As long as this interaction was present, the complex remained stable. Once this bond was broken, it resulted in complex dissociation. As can be seen, 1 does not form any long (>30% simulation time) lasting interactions with Aurka, which is a consequence of this complex dissociation. Another very important residue, in terms of complex stabilization, is Ala213. When it forms the direct H-bond with the ligand, such as with 2, 3 or AZD1152, such interaction is strong and stable. On the other hand, if this H-bond is water-mediated, it is eventually broken, resulting in complex dissociation such as in the complex of Aurka with 1.

The average number of interactions occurring simultaneously between Aurka and the ligand was 3–4 for the complexes formed with 2, 3 and 1-until the dissociation occurred . A significantly larger number of interactions were observed between AZD1152 and the studied protein, ranging from 9 to 12. In the case of the complex formed between Aurka and its most potent inhibitor, AZD1152, the phosphate group present in this ligand forms multiple strong interactions with Arg220, Lys224, and even one intramolecular H-bond with charged tertiary amine nitrogen. The other nitrogen atoms present in this ligand molecule form multiple H-bonds, either direct (with Glu211 and Ala213) or water mediated (Lys162, Thr217, Asp274). As can be seen in Figure 11, most of those interactions occur during the whole simulation time.

To analyze the behavior of the ligand in the form of the complex, Figure 12 was created. The RMSD for ligand was found to be the lowest for 2, indicating that this ligand conformation is the most stable. In the cases of 1, 3 and AZD1152, a few conformations of ligand are observed during simulation, which is represented by the different values of RMSD observed during calculations.

The radius of gyration (rGyr) of the ligand was found to be the largest for AZD1152 as it is also the largest of the studied ligands in terms of molecular size. A decrease in rGyr value during the simulation time, observed for 3, implies an increase in ligand structure compactness, thereby suggesting decreased flexibility and higher stability.

In the studied ligands, no intramolecular hydrogen bonds were observed, which is not surprising, taking into consideration their rigid structure and lack of H-bonds donor in the cases of 1 and 2. An exception to this observation is the behavior of AZD1152, in which a single intramolecular H-bond was formed after 10 ns of simulation, forming between one of the oxygen atoms of the phosphate group and a charged tertiary nitrogen amine atom. Interestingly, this interaction was not present at the beginning of the simulation. This may indicate that either this conformation is more stable under dynamic conditions or the MD was necessary to finding this deeper minimum.

The polar surface area (PSA) measures the solvent surface area that is produced by polar groups, such as nitrogen and oxygen atoms of the ligand during the simulation [38]. Polar surface area (PSA) was found to oscillate between 70 and 80 A³ for 1 and 2, while it reached 100 A³ for 3. Significantly larger values have been observed for AZD1152, exceeding 250 A³. This is not only due to the larger size of this ligand but also as a result of the presence of multiple basic nitrogen atoms and large phosphate group in this ligand.

The SASA refers to the solvent accessibility that provides the quantitative calculation of ligands with protein-implicit water molecules [38]. The studied ligands differ significantly in terms of this parameter. SASA for 2 oscillated between 50 and 75 A³, while SASA for 3 was found to be almost two times larger, despite the similar molecular masses of those ligands (328.26 and 305.42 Da, respectively). This was due to the significantly larger flexibility of 3 and the less compact shape of this molecule. Concluding, 3 is characterized by significantly larger flexibility and a less compact shape than 2; however, during the simulation, a decrease in the flexibility and an increase in compactness of 3 is also observed. In the case of 1, the 300% increase in SASA during the simulation was caused by the dissociation of the complex and the shift of the ligand from the binding site onto the surface of the protein. This process started after 20 ns of simulation and was completed at 70 ns.

The MolSA measures the water solvent area that ligand molecules acquire. In the MolSA study, the surface area is measured by applying a 1.4 Å probe radius, which is approximately equivalent to the van der Waals surface area and radius of one water molecule [38]. The observed MolSA was found to oscillate around 380, 330, 345 and 510 A³, respectively, for 1, 2, 3 and AZD1152. The decrease in MolSA for 3 was the direct result of switching the conformation by this ligand, as indicated by the decrease in RoG.

3. Materials and Methods

The protocol of molecular docking has been validated in 4 steps:

• In this study, 12 ligands with known affinity to Aurka expressed as Kd (PubChem AID 624919) were docked to Aurora A kinase using 3 programs (Surflex—Sybyl X 1.2, AutoDock 4.2, and Glide 6.6.), ;
• Kd values were correlated with scoring functions from 3 programs. The best results were obtained by Surflex, so it was recommended for further analysis ;
• To check the correctness of docking, ligand from crystal structure 3H10 was redocked, and the low-energy conformation was compared to conformation from the crystal structure ;
• Ligands with lower affinity to Aurka (structures and Kd from PubChem AID 624919) were docked to Aurora A using Surflex .

The set of ligands was selected from compounds obtained from peppers (Piper, Capsicum, and Pimenta), or their analogs with reported anticancer activity on the MCF-7 cell lines. Compounds 1 (piplartine analog), 4 (piplartine analog), 7 (chabamide), 10 (from P. regnelii), 12 (chabamide), 13 (piperine), 16 (from P. regnelii) were obtained from various Piper species. Compounds 6, 9, 11 are chalcones from Pimienta methysticum and Pimienta dilatatum) whereas 2, 5, 8 are their analogs. Compounds 3 (capsaicin), 14 and 15 (analogs) were obtained from the Capsicum species. Well-known Aurka inhibitors: AZD1152, hesperidin, PHA-73935 [23,24,25] were used for docking study as reference molecules. The structures of ligands were drawn in Sybyl-X 1.2 Sketch (Tripos International). Hydrogens were added and the geometry was optimized using Tripos forcefield (gradient 0.05 kcal/mol × A). The ligand preparation was performed by detecting the root, set number of torsions, and aromaticity criteria using AutoDock Tool [39].

Protein structures were taken from PDB entries 3H10 (Aurka) and 5IS0 (TRPV1). Molecular targets were prepared by removing water, adding hydrogens and performing the optimization using the Sybyl structure preparation tools. Then, protomol (an idealized active site ligand) was defined based on the ligand position in the crystal structures from PDB. The parameters were set as maximum conformations per fragment 20, a maximum number of poses per ligand 20 and maximum RMSD between final poses 0.05. Molecular docking was carried out using Surflex-Dock (SFXC) docking mode. The docking simulations by Autodock 4.2 program [40] applied a Lamarckian genetic algorithm to identify possible binding modes for all ligands. Ligands were prepared for docking by adding charges (gasteiger) and setting roots of torsions. The AutoDock grid box was centered on the ligand from the PDB crystal structure with the size set at 40 × 40 × 40 (x, y and z) and spacing between grid points 0.375 Å. The other parameters of AutoDock were used as default values. The three best compounds against Aurka and a known Aurka inhibitor (AZD1152) were chosen for further MD simulation [41].

The verification of docking was performed by comparing the conformation of capsaicin obtained in our study with that described earlier [32] and its interaction with TRPV1. Similarly, as for Aurka, the comparison was drawn between the results of molecular docking and the molecular mechanics calculation of hesperadin. Conformations of ligands and interactions in the active site were comparable with those described earlier; therefore, we assume that our results are reliable.

Molecular dynamics (MD) simulations were performed using Schrӧdinger Maestro 12.8. version (Schrӧdinger, LLC, New York, NY, USA, 2023).

3.1. Structures Preparation

Three-dimensional crystal structure of Aurka bound to the inhibitor 9-chloro-7-(2,6-difluorophenyl)-N-{4-[(4-methylpiperazin-1-yl)carbonyl]phenyl}-5H-pyrimido[5,4-d][2]benzazepin-2-amine (97B) at 2.20 Å resolution (PDB ID: 3H10) was retrieved from RCBS PDB database [42]. The structure was processed prior to docking using the Protein Preparation Wizard to remove unwanted water molecules and metal ions. This procedure also simplifies multimeric complexes, creates disulfide bonds, assigns bond orders properly, adjusts ionization states, and fixes the orientation of misoriented groups. Hydrogen atoms were added to the protein structure, and standard protonation states at pH 7.0 were used. The preprocessed structure was then optimized and minimized to generate geometrically stable structures. The prepared protein structure was then used for further modeling.

3.2. Molecular Docking

3.3. Molecular Dynamics (MD) Simulations

The obtained ligand–receptor complexes with the best Gscores constituted an input for MD simulations. For that purpose, the Desmond System Builder module was used. The ligand–protein complexes were put in the orthorhombic boxes. The systems were solvated with a water, using TIP3P water model [43], with a buffer distance of 10 Å. In each case, the system was neutralized by the addition of the appropriate number of Cl-ions. The systems were subjected to steepest descent minimization with Desmond’s default protocol prior to performing MD simulations.

The relaxation protocol consists of eight stages that included minimization with restraints on solute heavy atoms, minimization without any restraints, simulation with heating from 0 K to 300 K, H₂O barrier and gradual restraining, simulation under NPT equilibration with H₂O barrier with heavy atoms restrained, NPT equilibration of solvent and lipids, simulation under the NPT ensemble with protein heavy atoms restraint reduced from 10.0 to 2.0 kcal/mol, NPT equilibration with Cα atoms restrained at 2 kcal/mol, and simulation for 1.5 ns under the NPT ensemble with no restraints.

After relaxation, an unrestrained simulation run was performed for 100 ns for each system. The simulations were performed under the NPT ensemble using the Nose–Hoover thermostat to maintain a constant temperature of 300 K and isotropic Martyna–Tobias–Klein barostat to maintain the pressure at 1 atm. The short-range Coulombic interactions were analyzed with a cut-off value of 9.0 Å using the short-range method. A time-reversible reference system propagator algorithm (RESPA) integrator was used with a time step of 2.0 fs. The trajectories were saved at 100 ps intervals for analysis. After simulations were performed, RMSD, RMSF, and protein–ligand contacts were evaluated using a Simulation Interaction Diagram from the Schrödinger Suite. Interactions occurring in each frame of the performed simulations were encoded in the form of interaction fingerprints (IFPs).

The LogP (octanol–water coefficient) and molar mass values were calculated in HyperChem 7.5 for each ligand. PCA was evaluated using Statistica 10 (StatSoft Inc., Tulsa, OK, USA).

4. Conclusions

The conducted research showed that the compounds contained in various varieties of pepper can inhibit Aurora A kinase. Molecular docking allowed us to select the three best compounds showing kinase activity. Analysis of the docking results suggests that there is probably no relationship between the strength of the interaction with the TRPV1 receptor and the Aurora A kinase. The molecular dynamics study confirms that 2 and 3 form a stable complex with Aurora A. Our results show that the duration of the MD simulations is significant for the reliability of the results. In the case of short simulations, it would not be possible to detect that 1 forms an unstable complex. The best compounds selected from our study will be used in subsequent analyses using in vitro activity assays.

References

1. Giet, R.; Prigent, C.J. Aurora/Ipl1p-related kinases, a new oncogenic family of mitotic serine-threonine kinases. J. Cell Sci. 1999, 112, 3591–3601. [Google Scholar] [CrossRef] [PubMed]
2. Kollareddy, M.; Zheleva, D.; Dzubak, P.; Brahmkshatriya, P.S.; Lepsik, M.; Hajduch, M. Aurora kinase inhibitors: Progress towards the clinic. Investig. New Drugs 2012, 30, 2411–2432. [Google Scholar] [CrossRef]
3. Mountzios, G.; Terpos, E.; Dimopoulos, M.A. Aurora kinases as targets for cancer therapy. Cancer Treat. Rev. 2008, 34, 175–182. [Google Scholar] [CrossRef]
4. Lee, H.H.; Zhu, Y.; Govindasamy, K.M.; Gopalan, G. Downregulation of Aurora-A overrides estrogen-mediated growth and chemoresistance in breast cancer cells. Endocr.-Relat. Cancer 2008, 15, 765–775. [Google Scholar] [CrossRef]
5. Yang, N.; Wang, C.; Wang, J.; Wang, Z.; Huang, D.; Yan, M.; Kamran, M.; Liu, Q.; Xu, B. Aurora kinase A stabilizes FOXM1 to enhance paclitaxel resistance in triple-negative breast cancer. J. Cell. Mol. Med. 2019, 23, 6442–6453. [Google Scholar] [CrossRef]
6. Meulenbeld, H.J.; Mathijssen, R.H.; Verweij, J.; de Wit, R.; de Jonge, M. Danusertib, an aurora kinase inhibitor. Expert Opin. Investig. Drugs 2012, 21, 383–393. [Google Scholar] [CrossRef]
7. Manfredi, M.G.; Ecsedy, J.A.; Meetze, K.A.; Balani, S.K.; Burenkova, O.; Chen, W.; Galvin, K.M.; Hoar, K.M.; Huck, J.J.; LeRoy, P. Antitumor activity of MLN8054, an orally active small-molecule inhibitor of Aurora A kinase. Proc. Natl. Acad. Sci. USA 2007, 104, 4106–4111. [Google Scholar] [CrossRef] [PubMed]
8. Qi, W.; Cooke, L.S.; Liu, X.; Rimsza, L.; Roe, D.J.; Persky, A.M.D.O.; Miller, T.P.; Mahadevan, D. Aurora inhibitor MLN8237 in combination with docetaxel enhances apoptosis and anti-tumor activity in mantle cell lymphoma. Biochem. Pharmacol. 2011, 81, 881–890. [Google Scholar] [CrossRef]
9. Huh, H.C.; Lee, S.Y.; Lee, S.K.; Park, N.H.; Han, I.S. Capsaicin induces apoptosis of cisplatin-resistant stomach cancer cells by causing degradation of cisplatin-inducible Aurora-A protein. Nutr. Cancer 2011, 63, 1095–1103. [Google Scholar] [CrossRef] [PubMed]
10. Aziz, N.; Kim, M.Y.; Cho, J.Y. Anti-inflammatory effects of luteolin: A review of in vitro, in vivo, and in silico studies. J. Ethnopharmacol. 2018, 225, 342–358. [Google Scholar] [CrossRef]
11. Meng, X.-Y.; Zhang, H.-X.; Mezei, M.; Cui, M. Molecular docking: A powerful approach for structure-based drug discovery. Curr. Comp. Aided Drug Des. 2011, 7, 146–157. [Google Scholar] [CrossRef] [PubMed]
12. Raghav, P.K.; Kumar, R.; Kumar, V.; Raghava, G.P. Docking-based approach for identification of mutations that disrupt binding between Bcl-2 and Bax proteins: Inducing apoptosis in cancer cells. Mol. Genet. Genom. Med. 2019, 7, e910. [Google Scholar] [CrossRef] [PubMed]
13. Raghav, P.K.; Singh, A.K.; Gangenahalli, G. A change in structural integrity of c-Kit mutant D816V causes constitutive signaling. Mutat. Res. 2018, 808, 28–38. [Google Scholar] [CrossRef] [PubMed]
14. Raghav, P.K.; Verma, Y.K.; Gangenahalli, G.U. Molecular dynamics simulations of the Bcl-2 protein to predict the structure of its unordered flexible loop domain. J. Mol. Model. 2012, 18, 1885–1906. [Google Scholar] [CrossRef] [PubMed]
15. Puunganuru, S.R.; Madala, H.R.; Venugopal, S.N.; Samala, R.; Mikelis, C.; Srivenugopal, K.S. Design and synthesis of a C7-aryl piperlongumine derivative with potent antimicrotubule and mutant p53-reactivating properties. Eur. J. Med. Chem. 2016, 107, 233–244. [Google Scholar] [CrossRef] [PubMed]
16. T Tang, Y.; Li, X.; Liu, Z.; Simoneau, A.R.; Xie, J.; Zi, X. Flavokawain B, a kava chalcone, induces apoptosis via up-regulation of death-receptor 5 and Bim expression in androgen receptor negative, hormonal refractory prostate cancer cell lines and reduces tumor growth. Int. J. Cancer 2010, 127, 1758–1768. [Google Scholar] [CrossRef]
17. Bley, K.; Boorman, G.; Mohammad, B.; McKenzie, D.; Babbar, S. A comprehensive review of the carcinogenic and anticarcinogenic potential of capsaicin. Toxicol. Pathol. 2012, 40, 847–873. [Google Scholar] [CrossRef]
18. Rao, V.R.S.; Suresh, G.; Babu, K.S.; Raju, S.S.; Ramakrishna, S.; Rao, J.M. Novel dimeric amide alkaloids from Piper chaba Hunter: Isolation, cytotoxic activity, and their biomimetic synthesis. Tetrahedron 2011, 67, 1885–1892. [Google Scholar] [CrossRef]
19. Longato, G.B.; Rizzo, L.Y.; de Oliveira Sousa, I.M.; Tinti, S.V.; Possenti, A.; Figueira, G.M.; Ruiz, A.L.T.G.; Foglio, M.A.; de Carvalho, J.E. In vitro and in vivo anticancer activity of extracts, fractions, and eupomatenoid-5 obtained from Piper regnellii leaves. Planta Med. 2011, 77, 1482–1488. [Google Scholar] [CrossRef]
20. Rama Subba Rao, V.; Suresh, G.; Ranga Rao, R.; Suresh Babu, K.; Chashoo, G.; Saxena, A.; Madhusudana Rao, J. Synthesis of piperine–amino acid ester conjugates and study of their cytotoxic activities against human cancer cell lines. Med. Chem. Res. 2012, 21, 38–46. [Google Scholar] [CrossRef]
21. Fernandes, T.B.; de Azevedo, R.A.; Yang, R.; Teixeira, S.F.; Goulart Trossini, G.H.; Marzagao Barbuto, J.A.; Ferreira, A.K.; Parise-Filho, R. Arylsulfonylhydrazone induced apoptosis in MDA-MB-231 breast cancer cells. Lett. Drug Des. Discov. 2018, 15, 1288–1298. [Google Scholar] [CrossRef]
22. Sawasdee, K.; Chaowasku, T.; Lipipun, V.; Dufat, T.-H.; Michel, S.; Likhitwitayawuid, K. Neolignans from leaves of Miliusa mollis. Fitoterapia 2013, 85, 49–56. [Google Scholar] [CrossRef] [PubMed]
23. Raghav, P.K.; Verma, Y.K.; Gangenahalli, G.U. Peptide screening to knockdown Bcl-2’s anti-apoptotic activity: Implications in cancer treatment. Int. J. Mol. Biol. Macromol. 2012, 50, 796–814. [Google Scholar] [CrossRef]
24. Schellens, J.; Boss, D.; Witteveen, P.; Zandvliet, A.; Beijnen, J.; Voogel-Fuchs, M.; Morris, C.; Wilson, D.; Voest, E. Phase I and pharmacological study of the novel aurora kinase inhibitor AZD1152. J. Clin. Oncol. 2006, 24, 3008. [Google Scholar] [CrossRef]
25. Saiprasad, G.; Chitra, P.; Manikandan, R.; Sudhandiran, G. Hesperadin induces apoptosis and triggers autophagic markers through inhibition of Aurora-A mediated phosphoinositide-3-kinase/Akt/mammalian target of rapamycin and glycogen synthase kinase-3 beta signalling cascades in experimental colon carcinogenesis. Eur. J. Cancer 2014, 50, 2489–2507. [Google Scholar] [CrossRef] [PubMed]
26. Katayama, H.; Sen, S. Aurora kinase inhibitors as anticancer molecules. Biochim. Biophys. Acta Gene Regul. Mech. 2010, 1799, 829–839. [Google Scholar] [CrossRef] [PubMed]
27. Selleck Chemicals. Available online: https://www.selleckchem.com/ (accessed on 12 August 2022).
28. Talele, T.T.; McLaughlin, M.L. Modelling. Molecular docking/dynamics studies of Aurora A kinase inhibitors. J. Mol. Graph. Model. 2008, 26, 1213–1222. [Google Scholar] [CrossRef]
29. Siudem, P.; Paradowska, K.; Bukowicki, J. Conformational analysis of capsaicin using 13C, 15N MAS NMR, GIAO DFT and GA calculations. J. Mol. Struct. 2017, 1146, 773–781. [Google Scholar] [CrossRef]
30. Siudem, P.; Bukowicki, J.; Wawer, I.; Paradowska, K. Structural studies of two capsaicinoids: Dihydrocapsaicin and nonivamide. 13 C and 15 N MAS NMR supported by genetic algorithm and GIAO DFT calculations. RSC Adv. 2020, 10, 18082–18092. [Google Scholar] [CrossRef]
31. Siudem, P.; Paradowska, K. Structure, function, and mechanism of action of the vanilloid TRPV1 receptor. Postępy Hig. I Med. Doświadczalnej 2020, 74, 481–488. [Google Scholar] [CrossRef]
32. Darré, L.; Domene, C. Binding of capsaicin to the TRPV1 ion channel. Mol. Pharm. 2015, 12, 4454–4465. [Google Scholar] [CrossRef] [PubMed]
33. Leo, A.; Hansch, C.; Elkins, D. Partition coefficients and their uses. Chem. Rev. 1971, 71, 525–616. [Google Scholar] [CrossRef]
34. Pollastri, M.P. Overview on the Rule of Five. Curr. Protocol. Pharmacol. 2010, 49, 9–12. [Google Scholar] [CrossRef] [PubMed]
35. Chen, M.; Xiao, C.; Jiang, W.; Yang, W.; Qin, Q.; Tan, Q.; Lian, B.; Liang, Z.; Wei, C. Capsaicin inhibits proliferation and induces apoptosis in breast cancer by down-regulating FBI-1-mediated NF-κB pathway. Drug Des. Dev. Ther. 2021, 15, 125. [Google Scholar] [CrossRef]
36. Chou, C.-C.; Wu, Y.-C.; Wang, Y.F.; Chou, M.J.; Kuo, S.J.; Chen, D.R. Capsaicin-induced apoptosis in human breast cancer MCF-7 cells through caspase-independent pathway. Oncol. Rep. 2009, 21, 665–671. [Google Scholar]
37. Ying, Y.; Yulin, S.; Shuyan, L.; Nengzhi, J.; Huanxiang, L.; Xiaojun, Y. Molecular dynamics and free energy studies on Aurora kinase A and its mutant bound with MLN8054: Insight into molecular mechanism of subtype selectivity. Mol. BioSyst. 2012, 8, 3049–3060. [Google Scholar] [CrossRef]
38. Ansari, W.A.; Rab, S.O.; Saquib, M.; Sarfraz, A.; Hussain, M.K.; Akhtar, M.S.; Ahmad, I.; Khan, M.F. Pentafuhalol-B, a Phlorotannin from Brown Algae, Strongly Inhibits the PLK-1 Overexpression in Cancer Cells as Revealed by Computational Analysis. Molecules 2023, 28, 5853. [Google Scholar] [CrossRef]
39. Khan, M.F.; Ansari, W.A.; Ahamad, T.; Khan, M.A.; Khan, Z.A.; Sarfraz, A.; Khan, M.A. Bioactive components of different nasal spray solutions may defeat SARS-Cov2: Repurposing and in silico studies. J. Mol. Model. 2022, 28, 212. [Google Scholar] [CrossRef]
40. Morris, G.M.; Huey, R.; Lindstrom, W.; Sanner, M.F.; Belew, R.K. AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility. J. Comp. Chem. 2009, 30, 2785–2791. [Google Scholar] [CrossRef]
41. Ansari, W.A.; Khan, M.A.; Rizvi, F.; Ali, K.; Hussain, M.K.; Saquib, M.; Khan, M.F. Computational screening of plant-derived natural products against SARS-CoV-2 variants. Future Pharmacol. 2022, 2, 558–578. [Google Scholar] [CrossRef]
42. Berman, H.M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T.N.; Weissig, H.; Shindyalov, I.N.; Bourne, P.E. The Protein Data Bank. Nucleic Acids Res. 2000, 28, 235–242. [Google Scholar] [CrossRef] [PubMed]
43. Jorgensen, W.L.; Chandrasekhar, J.; Madura, J.D.; Impey, R.W.; Klein, M.L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 26. [Google Scholar] [CrossRef]