https://orcid.org/0000-0003-3165-929X
Figures
Abstract
Aurora kinases (AURKs) have been identified as promising biological targets for the treatment of cancer. In this study, molecular dynamics simulations were employed to investigate the binding selectivity of three inhibitors (HPM, MPY, and VX6) towards AURKA and AURKB by predicting their binding free energies. The results show that the inhibitors HPM, MPY, and VX6 have more favorable interactions with AURKB as compared to AURKA. The binding energy decomposition analysis revealed that four common residue pairs (L139, L83), (V147, V91), (L210, L154), and (L263, L207) showed significant binding energies with HPM, MPY, and VX6, hence responsible for the binding selectivity of AURKA and AURKB to the inhibitors. The MD trajectory analysis also revealed that the inhibitors affect the dynamic flexibility of protein structure, which is also responsible for the partial selectivity of HPM, MPY, and VX6 towards AURKA and AURKB. As expected, this study provides useful insights for the design of potential inhibitors with high selectivity for AURKA and AURKB.
Citation: Alharthy RD, Fatima G, Yousaf N, Iqbal MS, Sattar S, Alanzi AR, et al. (2023) Binding selectivity analysis of AURKs inhibitors through molecular dynamics simulation studies. PLoS ONE 18(12): e0295741. https://doi.org/10.1371/journal.pone.0295741
Editor: Chandrabose Selvaraj, Saveetha University - Poonamallee Campus: SIMATS Deemed University, INDIA
Received: October 5, 2023; Accepted: November 29, 2023; Published: December 19, 2023
Copyright: © 2023 Alharthy et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The raw data can be found at the following link: https://doi.org/10.6084/m9.figshare.24615636.
Funding: The author(s) received no specific funding for this work.
Competing interests: The authors declare there are no competing interests.
1. Introduction
The Aurora kinase group is composed of serine/threonine kinases, known as Aurora kinase A (AURKA), Aurora kinase B (AURKB), and Aurora kinase C (AURKC) [1]. AURKs serve an important role in controlling the cell cycle, with AURKA and AURKB being particularly important during mitosis [1], while AURKC, is crucial for gametogenesis [2]. The kinase domain of AURKs, containing three distinct domains, is highly homologous across all of its members [3]. But the N-terminal region’s sequence vary [3]. The location and spatiotemporal expression of AURKs clearly characterize their functional roles [4]. According to study, the overexpression of AURKs in malignancies results in genomic instability and aneuploidy [5] which causes the development, invasion, and spread of a tumor. Several studies emphasize the significance of AURKA in cancer treatment after extensive research into its functions [6].
AURKA and AURKB are two essential members of the serine/threonine kinases group [1, 7], with AURKA being associated with mitotic commitment, spindle construction, spindle maintenance, and centrosome function [8, 9]. TPX2 plays a crucial role in localizing AURKA to the mitotic spindle by binding and activating it. [10] The association between AURKA and TPX2 relies on the presence of glycine 198 (G198) in the catalytic domain of AURKA [10]. Additionally, the interaction between protein phosphatase 1 and AURKA, regulated by phosphorylation during mitosis, is vital for proper chromosomal segregation [11]. On the other hand, AURKB, encoded by the AURKB gene on chromosome 17, also plays a critical role in regulating the cell cycle. Both AURKA and AURKB phosphorylate histone H3, which is essential for chromosomal segregation during cell division [12]. AURKB is linked to the Chromosomal Passenger Complex, consisting of Survivin, Borealin, and INCENP, which plays a key role in various aspects of mitosis, including chromosome alignment, cytokinesis, and segregation.
Two significant regulators of cell division, aurora kinases A and B, have very similar amino acid sequences. The N-terminal domains, protein kinase domains, and C-terminal domains of Aurora A and B have an impressive conservation rate, with their catalytic domains being 71% identical [3]. Furthermore, the 26 residues around the ATP-binding active regions of both kinases are similar: In Aurora A, L215, T217, and R220, and in Aurora B, R159, E161, and K164, are the only residues differentiating their ATP-binding sites [13]. Despite these similarities, Aurora A and Aurora B have different chromosomal affiliations: Aurora A is linked to chromosome 20q13.2, whereas Aurora B is linked to chromosome 17p13.1 [14]. The average fraction of similar amino acids within the vertebrate Aurora A and B families is significantly larger (0.84 ± 0.5) than within either family alone, suggesting recent vertebrate evolution [15]. The high conservation rate is essential in relation to the distinctive pairing of substrates and inhibitors. These kinases interact with diverse substrates and subcellular localizations with minimal sequence change, despite their dissimilar structures and motifs. This highlights the essential functions that both kinds of kinases carry out in regulating cell division.Many studies have been reported that overexpression of AURKs are responsible in variety of human cancers, and the mutations in Aurora kinases have been identified in a variety of somatic cancer samples, that includes lung cancer, and melanoma [16], this suggests that the function of Aurora kinases in cell transformation and oncogenesis is crucial. In recent decades, there has been increased research on the role of these potentially oncogenic proteins in tumor growth.
The carboxyl terminus catalytic domain of AURKA and AURKB share approximately 70% similarity. Both kinases are essential for mitotic progression, but they have distinct localizations and roles. To investigate the reason for the difference between AURKA and AURKB, studies have used paired shRNA suppression with overexpression of Aurora mutants. Results showed that when the catalytic domain residue, glycine 198 is replaced with asparagine to mimetic the aligned asparagine 142 of Aurora B the AURKA bind to the AURKB binding protein INCENP, instead of TPX2 which is AURKA binding protein [10].
The Aurora B mitotic function is restored by the mutant Aurora A indicating that the binding to INCENP is important for AURKB’s unique functionality. Although AURKA needs G198 for TPX2 binding, and AURKB requires asparagine 142 for INCENP binding and function of AURKB [10].
Previous research has reported that certain compounds, including ZM447439 [17] and VX-680/MK-0457 [18] have demonstrated effects as Aurora kinase A inhibitors and Hesperidin [19] as AURKB. Despite significant experimental study on the interaction of inhibitors with AURKA and AURKB in different studies, decoding the atomic-level conformational changes of these two proteins due to inhibitor interactions is still crucial [20, 21].
With the rapid advancement of simulation and calculation methods [22], several molecular dynamics (MD) techniques, such as traditional MD [23, 24], multiple short molecular dynamics simulations [25, 26], accelerated MD (aMD) simulations [22, 27, 28], have been widely used to carry out target conformational evolution, various techniques for predicting binding free energy, including the Poisson Boltzmann surface area (MM-PBSA) method [29, 30], thermodynamics integration (TI) [31, 32], free energy perturbation (FEP) [33, 34], and solvated interaction energy (SIE) methods [35, 36], are frequently used to assess ligands’ capacity to bind to targets. Furthermore, techniques for machine learning and deep learning are presented to effectively examine the ligand-target binding process and reveal the underlying molecular causes of ligand-target interactions [37, 38]. These modeling techniques have also contributed to successful understandings of the inhibitor-receptor binding process.
The development of small drug like compounds that inhibit the activity of AURKs is still a major area of research. Two inhibitors MPY (2BMC) [39] and HPM (2C6E) [40] were designed to inhibit the activity of AURKA, while a small molecule VX6 (4AF3) [18] was developed to suppress the activity of AURKB. A crucial chemical process for the creation of small molecules that target AURKs can be found by further examining the differences in the binding patterns of VX6, HPM, and MPY to AURKA and AURKB. The plausible binding modes of these inhibitors are shown in Fig 1A and 1B, while the structures of HPM, MPY, and VX6 are shown in Fig 1C–1E. In this study, molecular dynamics simulations were used to enhance the conformational sampling of inhibitor-AURKs complexes, the cross-correlation matrix [41, 42] was used to understand the internal dynamics of inhibitor-bound AURKs, and calculations of residue-based free energy decomposition were used to identify the binding ability of VX6, HPM, and MPY to AURKs by employing MM/GBSA method. MM/GBSA offers a balanced approach, considering both molecular mechanics and solvation effects.
(A) Supxerposition of AURKA (green) and AURKB (blue) to analyze the conformation of bound inhibitors; (B) The plausible binding modes of inhibitors (represented with sticks) in AURKA and AURKB binding pocket (Shown in surface). (C), (D), and (E) represent the structures of HPM, MPY, and VX6, respectively.
2. Material and methods
2.1. Simulated system setup
The crystal structures of AURKA complexed with HPM (PDB ID: 2BMC) and MPY (PDB ID: 2C6E) and AURKB complexed with VX6 (PDB ID: 4AF3) were retrieved from Protein Data Bank (PDB). The co-crystal poses of the inhibitors were extracted from the structures and six complexes i.e., HPM-AURKA, MPY-AURKA, VX6-AURKA, HPM-AURKB, MPY-AURKB, VX6-AURKB were prepared using the PyMOL software. During structure preparation, all the non-inhibitor molecules and crystal water were removed, and the missing hydrogen atoms were added using the Amber’s Leap tool. The simulation parameters of the proteins were prepared by using the ff14SB forcefield [43, 44]. Similarly, the geometries of the inhibitors were optimized at the semiempirical AM1 and Amber’s Antechamber module [45, 46], and Gasteiger changes were allocated to each inhibitor. The general amber forcefield (GAFF) [47] was used to generate the parameters for HPM, MPY, and VX6 inhibitors. After parametrization, the solvation of the complexes was done in a periodic box of 10 Å containing TIP3P water molecules [48]. Then systems were neutralized by the addition of 7 Cl- ions in AURKA complexes and 8 Cl- ions in AURKB complexes. The initial confirmation of the complexes was subjected to 250ns MD simulation with randomly assigned velocities.
2.2. Molecular dynamics simulations
Before subjecting the system to the production run, the prepared systems were optimized through steepest decent minimization of 10000 steps to remove the unfavorable atomic interactions. After that, the solvation system was equilibrated for an additional 10000 steps. Then the temperature of the system was gradually raised from 0 to 300K and then the systems were further optimized at 300K for equilibration. The systems that undergone the process of equilibration were then subjected to the production run for 250 ns long simulation at 310K temperature and 1 atm pressure using NPT ensemble. The SHAKE algorithm was used to constrain the hydrogen bond forming atoms, and the particle mesh Ewald (PME) approach was used to identify the long-range electrostatic interactions [49, 50] at the cutoff range of 10Å. The molecular dynamic simulations were run by using NAMD [51]. The molecular dynamic trajectories were analyzed using the VMD [52] and BIO3D package of R [53].
2.3. Principle component analysis
PCA has proven to be a valuable method for identifying coordinated motions in a collection of conformational structures obtained from either molecular simulations or experimental data. This technique has been widely used to study how changes in conformation affect the function of receptors [54, 55]. To perform PCA, we used the atomic coordinates obtained from the molecular dynamics simulations to construct the covariance matrix. This matrix is then diagonalized to produce a set of eigenvalues and eigenvectors. After diagonalization, the eigenvector of the matrix demonstrates the directions of movement of protein domains, and the correlated eigenvalues indicate the square mean fluctuations along the respective eigenvectors. The first few eigenvectors with high eigenvalues are very helpful in demonstrating the overall movements of proteins BIO3D package [53] of R was used to compute the dynamic movement of the protein complexes.
2.4. Calculations of MM/GBSA
MM/GBSA method provides more reliable binding free energy values than many molecular docking scoring functions [56, 57]. Similarly some studies have also shown that MM/GBSA approach is accurate and reliable enough for predicting the small drug like compounds and their protein targets binding free energies [58, 59]. Considering these studies in view, binding free energies of AURKs complexes were calculated using MM/GBSA method by employing the below mentioned equation.
The ΔGcomp, ΔGpro, and ΔGlig indicate the binding energies of AURKs complexes. The ΔEele and ΔEvdW show the electrostatic and van der Waals interactions of the inhibitors to AURKs. The term ΔGgb presents the polar solvation energy which is solved by using the Generalized Born (GB) model [60] while ΔGnonpol represent the nonpolar free energy terms. Lastly, TΔS indicates the entropy caused by the ligands.
3. Results and discussion
3.1. AURKA and AURKB’s structural fluctuation and flexibilities
Root mean square deviations (RMSDs) of the backbone atoms from the initial optimized configuration were calculated for apo structures of AURKA and AURKB and their complexes to evaluate the extent of structural fluctuations across molecular dynamics simulations (Fig 2A and 2B). The information from RMSD plots revealed that the AURKA systems attained equilibrium at 25 ns. After equilibration, the RMSD of the apo AURKA increased to ~3.5 Å at 100 ns and then attained a stable range of 3 Å at 150 ns. The RMSD of AURKA-HPM complex showed deviations to 2.5–3.5 Å till 100 ns and then attained stability at ~3.5 Å till the end of simulation. The RMSD of AURKA-MPY complex showed a similar trend to HPM complex whereas the RMSD of the AURKA-VX6 complex showed major deviation during 100 to 150 ns but it attained stability after 150 ns in the range of ~2.5-3Å. Among the three complexes, the AURKA-HPM showed higher deviations than the other two complexes. On the other hand, AURKB complexes attained stability at the start of the simulation and their RMSD values remained in the range of ~2–3 Å throughout the simulation time while the RMSD of apo structure was higher than the complexes as indicated in Fig 2B. The RMSD analysis showed that the systems attained equilibrium at specific value and remained stable during the simulation. Furthermore, the conformational changes in the AURKS structures were analyzed by aligning the different snapshots obtained at 0, 50, 100, 150, 200, and 250 ns. The alignment of the snapshots revealed that the ligands remained bound to the proteins pockets and did not dissociate during the simulation (Fig 3).
(A) RMSD plots of AURKA complexes. (B) RMSD plots of AURKB complexes.
To assess the fluctuations of AURKA and AURKB structures upon binding of HPM, MPY, and VX6 inhibitors, the Root Mean Square Fluctuations (RMSF) of the residues were calculated from the trajectories (Fig 4). The RMSF plots of AURKA and AURKB showed similar trends indicating that both structures have same number of rigid residues and flexible regions. Structural fluctuations were observed in five regions consisting of L1, L2, L3, L4 and L5 regions. These results showed that some residues from these regions were situated near the binding sites of the AURKA and AURKB. The RMSF values of AURKA and AURKB bound to VX6 were higher than the AURKs bound to the HPM and MPY, especially at the L2, L3, and L5 regions indicating that the binding of HPM and MPY restricted the motions of these regions. The structural analysis revealed that the regions L2 and L4 were near the binding sites of AURKA and AURKB which showed that some residues in these regions play a significant role in the binding selectivity of HMP, MPY, and VX6 towards the AURKs. While the residues in L1, L3 and L5 were not near the binding pocket and the fluctuations in these regions can play a vital role in the binding of these inhibitors to the AURKA and AURKB. The RMSF values of apo structures revealed that the flexibility in the loop regions was more in apo structure as compared to the complexes indicating the stability of protein structures upon binding of the ligands [61].
RMSF of residues in AURKA and AURKB during MD simulations: (A) RMSF for apo AURKA and its complexes with HPM, MPY, and VX6, (B) the structure of AURKA, (C) RMSF of apo AURKB and its complexes with three inhibitors and (D) the structure of AURKB.
3.2. Dynamics behavior of AURKA and AURKB
The cross-correlation matrices of the Cα atomic coordinates of AURKA and AURKB complexes were computed to find the differences in structural dynamics (Fig 5). The red and pink colors show the positive correlated motions while blue and dark blue indicate the anti-correlated motions. The diagonal part of the matrix shows the correlated motion of domain relative to itself while the correlation in among different domains are depicted by the off-diagonal regions. As shown in Fig 5, the structural dynamics behavior of the AURKA and AURKB was influenced upon binding of the inhibitors, HPM, MPY, and VX6.
Cross correlation matrices derived from Cα atomic coordinates (A), (C), and (E) denoting AURKA with HPM, MPY, and VX6 respectively; (B), (D), and (F) depicting AURKB with HPM, MPY, and VX6.
For AURKB complexes (Fig 5B, 5D and 5F) the anticorrelated motions were observed at three regions R2, R3, and R4, while the positively correlated motions were observed at the diagonal and R1 region. By comparing the HPM bound to AURKB, the binding of HPM to AURKA weakens the positive correlating motions at R1 regions and anticorrelated motions at the R2, R3, and R4 regions (Fig 5A). The binding of MPY to AURKA slightly affected the anti-correlated motions at R2, R3, and R4 regions while the binding of MPY to AURKB affected the anti-correlated motions but did not affect the positive correlated motions at R1 region (Fig 5D). By comparing the AURKA-VX6 and AURKB-VX6, the associations did not alter the motions at R1 region in AURKA, but it strengthened the anticorrelated motions in R2, R3, and R4 regions in AURKA (Fig 5E). According to the above discussion, the binding of identical inhibitors to AURKA and AURKB lead to motion mode difference, indicating that the residues in R1–R4 regions may be involved in binding with HPM, MPY and VX6 [62].
3.3. PCA analysis
The use of principal component analysis (PCA) is widespread in the investigation of concerted movements of protein structural domains. This approach can effectively filter significant collective motions from structural ensembles obtained from experimental or simulation studies. In this study, PCA was utilized to decode the molecular mechanism underlying the binding selectivity of HPM, MPY, and VX6 to AURKA and AURKB. To perform PCA, a covariance matrix was constructed using the Cα atomic coordinates extracted from the starting conformational trajectory (SCT). The proportion of variance of first eigenvalue for the HPM bound to AURKA was 47.4% that was higher than the first eigenvalue of AURKB bound to HPM, indicating the higher structural variation in AURKA-HPM complex (Fig 6A and 6D). The binding of MPY to AURKA and AURKB did not show major difference in the variance of the structures as the eigenvalues for AURKA and AURKB were 22.9% and 22.5%, respectively (Fig 6B and 6E). Lastly, the binding of VX6 showed more variation in AURKA than AURKB. The results indicated that the binding of these inhibitors influence the structural variations of the AURKA and AURKB.
3.4. Binding free energy calculations
The MM/GBSA approach was utilized to determine the binding free energies of HPM, MPY, and VX6 to AURKA and AURKB. This involved calculating the energetic data for three hundred structural frames obtained from a 250 ns trajectory at 2 ns intervals. The results of the MM/GBSA calculation for all three ligands bound to AUKRA and AURKB are presented in Table 1.
The ΔEvdW value for AURKB-HPM was higher than AUKRA-HPM complex. The similar trend was observed in the ΔEvdW value of MPY and VX6 complexes. The ΔEele energy component for AURKA-HPM was -0.01±0.42, while it was -7.59±0.43 for HPM-AURKB complex. The electrostatic energy contributions indicated that the AURKA complexes showed better electrostatic energies than AURKB complexes except for HPM complexes. After analyzing the surface energies, the AURKB-HPM complex showed a better surface energy of -8.98±0.02 than all other complexes. The ΔGtotal values of HPM-AUKRA (-53.09±0.09), HPM-AURKB (-69.34±0.33), MPY-AURKA (-49.34±0.32), MPY-AURKB (-55.28±0.49), VX6-AURKA (-52.79±0.49), VX6-AURKB (-59.32±0.42) indicated that HPM had better selectivity towards AURKB rather than AURKA. Similarly, MPY and VX6 had better binding free energies with AURKB than AURKA. The relative contribution of energy terms in the AURKs complexes is shown in Fig 7.
3.5. Binding selectivity uncovered by inhibitor-residue interaction analyses
To clarify the binding selectivity of the three inhibitors (HPM, MPY, and VX6) to AURKA and AURKB, we used the MM/GBSA method to analyze the interactions between the inhibitors and specific residues. Table 2 Demonstrate the decomposition of ΔGligand–residue values into contributions from the sidechain and backbone of key residues in AURKA and AURKB bound by each inhibitor. We discovered that the sidechains of residues play an important role in contributing energy to inhibitor-residue interactions. Fig 8 displays the key residues of AURKA and AURKB that form important inhibitor–residue interactions with energies stronger than -1 kcal/mol [63]. Additionally, to determine hydrogen bond interactions (HBIs) between the inhibitors and AUKRA, and AURKB, we used the CPPTRAJ tool in Amber 21. The results are summarized in Table 3.
(A) the AURKA-HPM complex, (B) AURKB-HPM complex, (C) AURKA-MPY complex, (D) AURKB-MPY complex, (E) AURKA-VX6 complex and (F) AURKB-VX6 complex.
For the HPM bound to AURKA and AURKB, HPM showed better interactions with L139, V147, L194, L210, Y212, G216, L263, and F275 in AURKA and all the interactions were stronger than -1kcal/mol (Fig 8A). Additionally, a Hydrogen bond interaction (Ala87-H——HPM-N17) was detected with an occupancy of 78.23% (Table 3). The HBI indicated that the hydrogen atom of Ala87 residues engaged in hydrogen bonding with the nitrogen seventeen of HPM. On comparison with the binding of HPM to AURKB, it was observed that the binding modes of HPM with AURKB were like AURKA (Fig 8B). It was observed that the interactions difference of HPM with residues (L139, L83), (V147, V91), (L194, L138), (L210, L154), (G216, G160), (L263, L207) corresponding to AURKA and AURKB was less than 0.30 kcal/mol. The interaction energy difference of HPM with F275 in AURKA and F219 in AURKB was strengthened by 1.20 kcal/mol, indicating that these residues play a significant role in binding selectivity of HPM to AURKA over AURKB.
In case of MPY binding to the AURKA, the interactions stronger than -1.00 kcal/mol were observed in the residues L139, V147, Y212, L215, and L263 (Fig 8C). The hydrogen bond interactions showed that the HA atom of Asp148 engaged in hydrogen bonding with the O36 of MPY with an occupancy of 44.09%. The binding mode of MPY in AURKB was like AURKA as the interacting residues (L139, L83), (V147, V91), (Y212, Y156), and (L263, L207) depicted the better interactions with energy stronger than -1.00 kcal/mol. The binding interactions of L207 and Y212 were strengthened by 0.95kcal/mol and 0.87kcal/mol, indicating these residues were involved in the binding selectivity of MPY towards AURKA and AURKB [64].
Lastly, VX6 bound to AURKA and AURKB interacted with four residues, (L139, L83), (V147, V91), (Y212, Y156), and (L263, L207) as show in (Fig 8E and 8F). VX6 formed HBI with AURKA (Ala87-H-----VX6-N14), and AURKB (Ala88-H-----VX6-N14) with an occupancy of 64.48% and 74.69%, respectively. The interactions of VX6 with L139 and L263 in AURKA were strengthened by -2.29 and -2.32kcal/mol as compared to the binding of VX6 to L83 and L207 with AURKB, which showed that these residues are important for the binding of VX6 to AURKA and AURKB. The binding modes of the inhibitors with the interacting residues are shown in Fig 9.
Further, the distances between the hydrogen bond forming residues with the inhibitors were calculated to analyze the stability of the complexes. The analysis of AURKA complexes showed that the distance between HPM and Ala87 at the start of simulation was 3.25 Å which gradually increased to 4.5 Å at 60 ns but it dropped to 2.5 Å at 75 ns and then remained in this range till 225 ns. The distance increased to 5 Å towards the end of simulation. Similarly, the average distance between MPY and Asp148 3.5 Å during 250 ns simulation. The distance between VX6 and Ala87 was 3.5 Å in most of the frames while in some frames it reached 5 Å (Fig 10). On the other hand, the distance between HPM and Lys38 of AURKB was in the range of 3.5 Å most of the time. The average distance between MPY and Gly91 was 4 Å and the distance between VX6 and Ala88 was 3 Å (Fig 11).
The role of residual contribution in the total binding free energy was further estimated by alanine scanning. Table 2 indicates that F275 play an important role in binding free energy of HPM-AURKA complex. In MPY-AUKRA complex, Y212 showed better contribution in binding free energy while in VX6-AURKA complex, the major energy contribution was exhibited by L263. In the AURKB complexes, F88 showed major contributions in all complexes. These residues were mutated to alanine and then their effect on total binding free energy was calculated by alanine scanning. The reduction in the binding free energies upon mutation is shown in Table 4. The mutation of F275 to alanine in HPM-AURKA complex reduced the binding free energy to -49.07 with an energy loss of -4.01 kcal/mol. The mutation of Y212 in MPY-AURKA caused an energy loss of -2.23 kcal/mol while the mutation of L263 in VX6-AURKA complex resulted in an energy loss of -4.73 kcal/mol. In AURKB complexes, the mutation of F88 to alanine showed the more energy loss in HPM-AURKB complex with a value of -4.07 kcal/mol. The loss of energy due to mutation in these residues indicated the significance of these residues in total binding free energy of complexes.
4. Conclusion
Molecular dynamics simulations were conducted on three AURKA and AURKB systems bound by three inhibitors (HPM, MPY, and VX6) to gain insights into their binding selectivity for anti-cancer drug development. The study revealed that the structural flexibility of two regions L4 and L5 in AURKA was higher than that of AURKB, and these domains exhibited distinct internal dynamics behavior. MM/GBSA calculations showed that the binding free energies of HPM, MPY, and VX6 with AURKB were higher than the AURKA complexes, indicating the superior selectivity and binding ability of inhibitors towards AURKB. Furthermore, residue-based free energy decomposition analysis identified four common residue pairs (L139, L83), (V147, V91), (L210, L154), and (L263, L207) that played a significant role in inhibitor binding affinity, suggesting their crucial role in determining the selectivity of inhibitors towards AURKA and AURKB. These findings could provide a better understanding of the molecular mechanisms and structure-affinity relationships for designing highly selective AURKs inhibitors.
Acknowledgments
The present study was supported by the Researchers Supporting Project (grant no. RSPD2023R885), King Saud University, Riyadh, Saudi Arabia.
References
- 1. Nigg E.A., Mitotic kinases as regulators of cell division and its checkpoints. Nature Reviews Molecular Cell Biology, 2001. 2(1): p. 21–32. pmid:11413462
- 2. Dieterich K., et al., Homozygous mutation of AURKC yields large-headed polyploid spermatozoa and causes male infertility. Nature Genetics, 2007. 39(5): p. 661–665. pmid:17435757
- 3. Carmena M. and Earnshaw W.C., The cellular geography of aurora kinases. Nat Rev Mol Cell Biol, 2003. 4(11): p. 842–54. pmid:14625535
- 4. Li S., et al., Spatial Compartmentalization Specializes the Function of Aurora A and Aurora B. Journal of Biological Chemistry, 2015. 290(28): p. 17546–17558. pmid:25987563
- 5. Marumoto T., Zhang D., and Saya H., Aurora-A—A guardian of poles. Nature Reviews Cancer, 2005. 5(1): p. 42–50. pmid:15630414
- 6. Tang A., et al., Aurora kinases: novel therapy targets in cancers. Oncotarget, 2017. 8(14): p. 23937–23954. pmid:28147341
- 7. Chan C.S. and Botstein D., Isolation and characterization of chromosome-gain and increase-in-ploidy mutants in yeast. Genetics, 1993. 135(3): p. 677–91. pmid:8293973
- 8. Bischoff J.R. and Plowman G.D., The Aurora/Ipl1p kinase family: regulators of chromosome segregation and cytokinesis. Trends in Cell Biology, 1999. 9(11): p. 454–459. pmid:10511710
- 9. Giet R. and Prigent C., Aurora/Ipl1p-related kinases, a new oncogenic family of mitotic serine-threonine kinases. J Cell Sci, 1999. 112 (Pt 21): p. 3591–601. pmid:10523496
- 10. Hans F., et al., Molecular distinctions between Aurora A and B: a single residue change transforms Aurora A into correctly localized and functional Aurora B. Mol Biol Cell, 2009. 20(15): p. 3491–502. pmid:19494039
- 11. Katayama H., et al., Interaction and Feedback Regulation between STK15/BTAK/Aurora-A Kinase and Protein Phosphatase 1 through Mitotic Cell Division Cycle *. Journal of Biological Chemistry, 2001. 276(49): p. 46219–46224. pmid:11551964
- 12. Crosio C., et al., Mitotic phosphorylation of histone H3: spatio-temporal regulation by mammalian Aurora kinases. Mol Cell Biol, 2002. 22(3): p. 874–85. pmid:11784863
- 13. Long L., et al., Structure-based drug design: Synthesis and biological evaluation of quinazolin-4-amine derivatives as selective Aurora A kinase inhibitors. Eur J Med Chem, 2018. 157: p. 1361–1375. pmid:30196060
- 14. Bischoff J.R., et al., A homologue of Drosophila aurora kinase is oncogenic and amplified in human colorectal cancers. Embo j, 1998. 17(11): p. 3052–65. pmid:9606188
- 15. Hall M., Lindholm A.K., and Brooks R., Direct selection on male attractiveness and female preference fails to produce a response. BMC Evol Biol, 2004. 4: p. 1. pmid:14720309
- 16. Dutertre S., et al., The absence of p53 aggravates polyploidy and centrosome number abnormality induced by Aurora-C overexpression. Cell Cycle, 2005. 4(12): p. 1783–7. pmid:16258285
- 17. Ditchfield C., et al., Aurora B couples chromosome alignment with anaphase by targeting BubR1, Mad2, and Cenp-E to kinetochores. J Cell Biol, 2003. 161(2): p. 267–80. pmid:12719470
- 18. Harrington E.A., et al., VX-680, a potent and selective small-molecule inhibitor of the Aurora kinases, suppresses tumor growth in vivo. Nat Med, 2004. 10(3): p. 262–7. pmid:14981513
- 19. Hauf S., et al., The small molecule Hesperadin reveals a role for Aurora B in correcting kinetochore-microtubule attachment and in maintaining the spindle assembly checkpoint. J Cell Biol, 2003. 161(2): p. 281–94. pmid:12707311
- 20. Clegg M.A., et al., Application of Atypical Acetyl-lysine Methyl Mimetics in the Development of Selective Inhibitors of the Bromodomain-Containing Protein 7 (BRD7)/Bromodomain-Containing Protein 9 (BRD9) Bromodomains. J Med Chem, 2020. 63(11): p. 5816–5840. pmid:32410449
- 21. Park S.W. and Lee J.M., Emerging Roles of BRD7 in Pathophysiology. Int J Mol Sci, 2020. 21(19). pmid:32992509
- 22. Wang J. and Miao Y., Mechanistic Insights into Specific G Protein Interactions with Adenosine Receptors. J Phys Chem B, 2019. 123(30): p. 6462–6473. pmid:31283874
- 23. Wu S.L., et al., Insights into interaction mechanism of inhibitors E3T, E3H and E3B with CREB binding protein by using molecular dynamics simulations and MM-GBSA calculations. SAR QSAR Environ Res, 2021. 32(3): p. 221–246. pmid:33661069
- 24. Xue W., et al., What Contributes to Serotonin-Norepinephrine Reuptake Inhibitors’ Dual-Targeting Mechanism? The Key Role of Transmembrane Domain 6 in Human Serotonin and Norepinephrine Transporters Revealed by Molecular Dynamics Simulation. ACS Chem Neurosci, 2018. 9(5): p. 1128–1140. pmid:29300091
- 25. Wang Y., et al., Binding selectivity of inhibitors toward the first over the second bromodomain of BRD4: theoretical insights from free energy calculations and multiple short molecular dynamics simulations. RSC Advances, 2021. 11(2): p. 745–759.
- 26. Chen J., et al., Molecular mechanism with regard to the binding selectivity of inhibitors toward FABP5 and FABP7 explored by multiple short molecular dynamics simulations and free energy analyses. Physical Chemistry Chemical Physics, 2020. 22(4): p. 2262–2275. pmid:31917380
- 27. Duan L., et al., Accelerated Molecular Dynamics Simulation for Helical Proteins Folding in Explicit Water. Front Chem, 2019. 7: p. 540. pmid:31448259
- 28. Chen J., et al., Q61 mutant-mediated dynamics changes of the GTP-KRAS complex probed by Gaussian accelerated molecular dynamics and free energy landscapes. RSC Advances, 2022. 12(3): p. 1742–1757. pmid:35425180
- 29. Chen J., et al., Binding mechanism of inhibitors to p38α MAP kinase deciphered by using multiple replica Gaussian accelerated molecular dynamics and calculations of binding free energies. Comput Biol Med, 2021. 134: p. 104485.
- 30. Wang W. and Kollman P.A., Computational study of protein specificity: The molecular basis of HIV-1 protease drug resistance. Proceedings of the National Academy of Sciences, 2001. 98(26): p. 14937–14942. pmid:11752442
- 31. Chen J., et al., A Comparative Insight into Amprenavir Resistance of Mutations V32I, G48V, I50V, I54V, and I84V in HIV-1 Protease Based on Thermodynamic Integration and MM-PBSA Methods. J Chem Inf Model, 2015. 55(9): p. 1903–13. pmid:26317593
- 32. Tzoupis H., et al., A Comparative Molecular Dynamics, MM-PBSA and Thermodynamic Integration Study of Saquinavir Complexes with Wild-Type HIV-1 PR and L10I, G48V, L63P, A71V, G73S, V82A and I84V Single Mutants. J Chem Theory Comput, 2013. 9(3): p. 1754–64. pmid:26587633
- 33. Aldeghi M., et al., Predictions of Ligand Selectivity from Absolute Binding Free Energy Calculations. J Am Chem Soc, 2017. 139(2): p. 946–957. pmid:28009512
- 34. Chen J., et al., Effect of mutations on binding of ligands to guanine riboswitch probed by free energy perturbation and molecular dynamics simulations. Nucleic Acids Res, 2019. 47(13): p. 6618–6631. pmid:31173143
- 35. Naïm M., et al., Solvated interaction energy (SIE) for scoring protein-ligand binding affinities. 1. Exploring the parameter space. J Chem Inf Model, 2007. 47(1): p. 122–33. pmid:17238257
- 36. Yan F., et al., Understanding conformational diversity of heat shock protein 90 (HSP90) and binding features of inhibitors to HSP90 via molecular dynamics simulations. Chem Biol Drug Des, 2020. 95(1): p. 87–103. pmid:31560152
- 37. Bennett W.F.D., et al., Predicting Small Molecule Transfer Free Energies by Combining Molecular Dynamics Simulations and Deep Learning. J Chem Inf Model, 2020. 60(11): p. 5375–5381. pmid:32794768
- 38. Gainza P., et al., Deciphering interaction fingerprints from protein molecular surfaces using geometric deep learning. Nat Methods, 2020. 17(2): p. 184–192. pmid:31819266
- 39. Yan A., et al., Aurora-A kinase inhibitor scaffolds and binding modes. Drug Discov Today, 2011. 16(5–6): p. 260–9. pmid:21147253
- 40. Backes A., et al., Small-molecule inhibitors binding to protein kinase. Part II: the novel pharmacophore approach of type II and type III inhibition. Expert Opin Drug Discov, 2008. 3(12): p. 1427–49. pmid:23506107
- 41. Chen J., et al., Binding of Inhibitors to BACE1 Affected by pH-Dependent Protonation: An Exploration from Multiple Replica Gaussian Accelerated Molecular Dynamics and MM-GBSA Calculations. ACS Chem Neurosci, 2021. 12(14): p. 2591–2607. pmid:34185514
- 42. Ichiye T. and Karplus M., Collective motions in proteins: A covariance analysis of atomic fluctuations in molecular dynamics and normal mode simulations. Proteins: Structure, Function, and Bioinformatics, 1991. 11(3): p. 205–217.
- 43. Song D., Luo R., and Chen H.F., The IDP-Specific Force Field ff14IDPSFF Improves the Conformer Sampling of Intrinsically Disordered Proteins. J Chem Inf Model, 2017. 57(5): p. 1166–1178. pmid:28448138
- 44. Maier J.A., et al., ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J Chem Theory Comput, 2015. 11(8): p. 3696–713. pmid:26574453
- 45. Jakalian A., Jack D.B., and Bayly C.I., Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation. J Comput Chem, 2002. 23(16): p. 1623–41. pmid:12395429
- 46. Jakalian A., et al., Fast, efficient generation of high‐quality atomic charges. AM1‐BCC model: I. Method. Journal of computational chemistry, 2000. 21(2): p. 132–146.
- 47. Wang J., et al., Development and testing of a general amber force field. J Comput Chem, 2004. 25(9): p. 1157–74. pmid:15116359
- 48. Jorgensen W.L., et al., Comparison of simple potential functions for simulating liquid water. The Journal of chemical physics, 1983. 79(2): p. 926–935.
- 49. Darden T., York D., and Pedersen L., Particle mesh Ewald: An N⋅ log (N) method for Ewald sums in large systems. The Journal of chemical physics, 1993. 98(12): p. 10089–10092.
- 50. Essmann U., et al., A smooth particle mesh Ewald method. The Journal of chemical physics, 1995. 103(19): p. 8577–8593.
- 51. Phillips J.C., et al., Scalable molecular dynamics on CPU and GPU architectures with NAMD. 2020. 153(4): p. 044130.
- 52. Humphrey W., Dalke A., and Schulten K., VMD: visual molecular dynamics. Journal of molecular graphics, 1996. 14(1): p. 33–38. pmid:8744570
- 53. Grant B.J., Skjærven L., and Yao X.Q.J.P.S., The Bio3D packages for structural bioinformatics. 2021. 30(1): p. 20–30.
- 54. Chen J., et al., Free energy profiles relating with conformational transition of the switch domains induced by G12 mutations in GTP-bound KRAS. 2022. 9.
- 55. Yang L.-Q., et al., Substrate-induced changes in dynamics and molecular motions of cuticle-degrading serine protease PL646: a molecular dynamics study. 2017. 7(67): p. 42094–42104.
- 56. Crean R.M., et al., Reliable in silico ranking of engineered therapeutic TCR binding affinities with MMPB/GBSA. 2022. 62(3): p. 577–590.
- 57. Wang C., et al., Recent developments and applications of the MMPBSA method. 2018. 4: p. 87.
- 58. Wang E., et al., End-point binding free energy calculation with MM/PBSA and MM/GBSA: strategies and applications in drug design. 2019. 119(16): p. 9478–9508.
- 59. Sun H., et al., Assessing the performance of MM/PBSA and MM/GBSA methods. 5. Improved docking performance using high solute dielectric constant MM/GBSA and MM/PBSA rescoring. 2014. 16(40): p. 22035–22045.
- 60. Onufriev A., et al., Exploring protein native states and large‐scale conformational changes with a modified generalized born model. 2004. 55(2): p. 383–394.
- 61. de Souza A.S., et al., 3-Acyltetramic acids as a novel class of inhibitors for human kallikreins 5 and 7. 2019. 29(9): p. 1094–1098.
- 62. Wang L., et al., Binding selectivity-dependent molecular mechanism of inhibitors towards CDK2 and CDK6 investigated by multiple short molecular dynamics and free energy landscapes. 2023. 38(1): p. 84–99.
- 63. Wang L., et al., Binding Selectivity of Inhibitors toward Bromodomains BAZ2A and BAZ2B Uncovered by Multiple Short Molecular Dynamics Simulations and MM-GBSA Calculations. 2021. 6(18): p. 12036–12049.
- 64. Wang L., et al., Deciphering Selectivity Mechanism of BRD9 and TAF1 (2) toward Inhibitors Based on Multiple Short Molecular Dynamics Simulations and MM-GBSA Calculations. 2023. 28(6): p. 2583.