Exploration of the selective binding mechanism of protein kinase Aurora A selectivity via a comprehensive molecular modeling study

Preparation of the studied systems

The crystal structures of human Aurora A/AT9283 (PDB code: 6CPG) and Aurora A/Danusertib (PDB code: 2J50) were obtained from The Protein Data Bank (PDB) web server (Smith et al., 2015). The missing side-chains and loop structures in the crystal structures of Aurora A were refined using the Loops/Refine Structure module in UCSF Chimera program (Pettersen et al., 2004). The protonation states of the Aurora A residues were estimated using PDB2PQR web Server (Dolinsky et al., 2004). The charge state of the Gleevec, AT9283 and Danusertib were set as +1, 0, 0, respectively. The initial structure of Aurora A/Gleevec was implemented using molecular docking in the AutoDock 4.2 program (Morris et al., 2009). A cubic box of 25 Å × 25 Å × 25 Å was defined, with ATP-binding sites as the center and a grid spacing of 0.3 Å. The AutoDockTools program was used to distribute Gasteiger partial charges to the atoms of Gleevec (Morris et al., 2009). AutoGrid software was applied to estimate the affinity maps of Aurora A. The docking parameters were set as follows: trials of 300 dockings and clustered according to the RMSD tolerance of 1.0 Å, maximum number of evaluations was set as 25,000,000, and other parameters were set at default values. The highest ranking structure for Aurora A/Gleevec was selected for further molecular dynamics (MD) simulation protocols.

Conventional MD simulation

Conventional MD (cMD) simulations were performed using Assisted Model Building with Energy Refinement (Amber) 16 software. The restrained electrostatic potential (RESP) method was utilized to calculate the partial atomic charges for AT9283, Danusertib and Gleevec based on HF/6-13G* basis set. Then, the Amber ff14SB force field and the General Amber Force Field (GAFF) in the LEaP modules were implemented to describe the Aurora A and ligands in Amber 16 software, respectively (Maier et al., 2015; Wang et al., 2004). Subsequently, these systems were solvated in a cubic box of TIP3P water molecules with the distance between the complex surface and the box boundary set at 20 Å. To ensure the overall charge neutrality, an appropriate number of counter ions were added.

A three-step minimization protocol was implemented to optimize the systems prior to the cMD simulations. Firstly, the 6,000-step steepest descent (SD) method was utilized to minimize the water molecules. This was followed by the implementation of the 6,000-step conjugate gradient (CG) method in order to keep the proteins in a position other than the hydrogens. Then, the same minimization scheme was applied to optimize the protein side chains. At last, the whole system was relaxed for the 10,000-step SD method followed by the 6,000-steps CG method without any restraints. Then, the optimized systems were progressively heated from 0 K to 100 K and then to 310 K with the protein restrained over 300 ps in the canonical (NVT) ensemble. To accommodate solvent density, the whole system was equilibrated over 500 ps at a constant pressure and for 5 ns in the isothermal isobaric (NPT) ensemble. Finally, the unrestrained 500 ns cMD simulations were carried out in the NPT ensemble to produce trajectories with the temperature set at 310 K. During the simulation process, temperature was regulated by Langevin dynamics and the pressure was controlled using Berendsen barostat (Izaguirre et al., 2001; Larini, Mannella & Leporini, 2007). Particle-mesh Ewald (PME) was used to consider the long-range electrostatics with a cutoff value of 10 Å for the non-bonding interactions, and the SHAKE method was selected to maintain the rigidity of all bonds (Essmann et al., 1995; Krautler, Van Gunsteren & Hunenberger, 2001). The time step was set to 2 fs, and the trajectory frames were saved at an interval of 10 ps for further analysis.

Dynamic cross-correlation (DCC) analysis

To better understand the dynamic behavior of the simulated systems, DCC analysis was conducted to evaluate the cross-correlation displacement of the protein backbone atoms (C_α). The cross-correlation matrix (C_ij) between residues i and j were calculated using the Bio3D package in R (Skjaerven et al., 2014). The C_ij was determined by the following equation: (1)${\displaystyle C_{ij}=\frac{〈\Delta r_{i\cdot }\Delta r_j〉}{\sqrt{〈\Delta r_i^2\Delta r_j^2〉}}}$ where Δr_i or Δr_j represents the displacement from the mean position of the ith or jth atom, and angle bracket denotes an average over the sampled period. In this study, only C_α was applied for analysis by averaging motions of C_α atoms deviating from the mean structure based on last 200 ns trajectories from cMD and GaMD simulations with a total of 5,000 snapshots for each method.

Principal component analysis (PCA)

PCA was performed to study the overall motion of Aurora A protein in each simulated system. It was performed through a commonly used protocol, which involved the elimination of the translational and rotational motions of the protein C_α atoms by aligning the structures from cMD/GaMD simulation trajectories. Then, the 3N ×3N covariance matrix was created using Cartesian coordinates, followed by construction of eigenvectors by diagonalization of the covariance matrix. Herein, PCA was calculated using the trajectories from the 500 ns cMD simulation and 600 ns GaMD simulation trajectories.

Thermodynamic calculations

In accordance with a standard protocol reported in previous studies, the binding free energies and per-residue contribution to ligands were calculated using the Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) method based on the single cMD simulation trajectory (Piao et al., 2019; Sun et al., 2018; Zheng et al., 2018). The ΔG_binding is given by: (2)${\displaystyle G_{binding}=G_{R+L}-\left(G_R+G_L\right)}$ (3)${\displaystyle 〈G〉=〈E_{\mathrm{MM}}〉+〈G_{\mathrm{sol}}〉-T〈S_{\mathrm{MM}}〉}$ (4)${\displaystyle E_{\mathrm{MM}}=E_{\mathrm{int}}+E_{\mathrm{vdW}}+E_{\text{elec}}}$ (5)${\displaystyle G_{\mathrm{sol}}=G_{\mathrm{GB}}+G_{\mathrm{SA}}}$ where ΔG_R+L, ΔG_R, and ΔG_L represent the free energies of the receptor–ligand complex, receptor, and ligand, respectively (Eq. 2). The free energy term is calculated as an average over the ensemble of time-equidistant snapshots (Eq. 3). The molecular mechanical energy (E_MM) can be divided into three terms (Eq. 4): intermolecular energy (E_int), van der Waals energy (E_vdW), and electrostatic energy (ΔE_elec). The solvation term (G_solv) consists of polar (G_SA) and nonpolar contributions (G_GB). In this study, the G_GB was estimated using the Generalized Born (GB) model reported by Onufriev, Bashford & Case (2004) (igb = 2); the dielectric of solute and solvent were set to 1 and 80, respectively. The G_SA is estimated using a fast linear combination of the pairwise overlap (LCPO) model with a probe radius of 1.4 Å. Conformational entropy (−TΔS) was not considered due to the high computational demand and low accuracy (Hou et al., 2011). A total of 2,000 snapshots provided by last 200 ns simulation trajectories were used for binding free energy and per-residue energy contribution calculations.

Gaussian accelerated MD (GaMD) simulations

The equilibrated snapshots extracted from cMD simulations served as the initial structures for the GaMD simulations. In our study, dual-boost GaMD simulations were carried out, which means that both the dihedral energy and total potential energy of the system were added to the simulated systems. The threshold energy was set to the lower bound (E = V_max), and the average and standard deviation (Std) of the simulated system potential energies were calculated every 200 ps. The upper limits of the boost potential Std of both the dihedral and total potential energies were set as 6.0 kcal/mol. Herein, a 10 ns cMD simulation was used to collect the maximum, minimum, average, and Std values of the potential of simulated system for obtaining the acceleration parameters of GaMD. Finally, 600 ns GaMD simulations were carried out with random atomic velocity initializations. During the simulation, temperature was regulated by Langevin dynamics and pressure was controlled using Berendsen barostat (Izaguirre et al., 2001; Larini, Mannella & Leporini, 2007). Particle-mesh Ewald (PME) was used to consider the long-range electrostatics with a cutoff of 10 Å for the non-bonded interactions, and the SHAKE method was adopted to maintain all bonds rigid (Essmann et al., 1995; Krautler, Van Gunsteren & Hunenberger, 2001). The time step was set to 2 fs, and the trajectory frames were saved at an interval of 20 ps for further analysis. After the GaMD simulation, the principal component 1 (PC1) and principal component 2 (PC2) calculated from PCA in conjunction with the potential boost provided by GaMD simulations were performed to obtain the original free energy map (FEM) by the cumulant expansion to the second order method (Miao et al., 2014; Roe & Cheatham 3rd, 2013).

Umbrella sampling (US) simulations

The equilibrated structures extracted from the trajectories from MD simulations served as the initial structures for the umbrella sampling simulations. The direction of the reaction coordinates (RCs) along the largest binding pocket was determined using the CAVER Analyst 1.0 software (Kozlikova et al., 2014). Subsequently, the distance between one atom in Aurora A (C_α of His-187) and another carbon atom in ligand was selected as the RCs (Fig. S1). In this study, the RCs of the studied systems were extended 20 Å from the initial position unbound from ATP binding pocket. For each system, the RCs were separated into 41 windows, each with a length of 0.5 Å. The initial conformation for each window was from the last snapshot of the previous window and each window was simulated with an 8 ns umbrella sampling simulation to ensure system convergence. A constant force of 10 kcal mol⁻¹ Å⁻² was added to each ligand to pull away from the binding pocket. The weighted histogram analysis method (WHAM) was used to calculate the potential of mean force (PMF), and the biased probability distribution was converted into a normal distribution (Kumar et al., 1992).

Structural equilibrium and flexibility of the studied systems

The binding mode from molecular docking prompted us to explore the dynamic behavior of the simulated complexes. To analyze the stability and flexibility of the system, the root-mean square deviations (RMSDs) of the C_α atoms relative to their original position were monitored. As plotted in Fig. 2, the RMSD evolution from the cMD simulations of the backbone atoms (C_α) of the Apo Aurora A, Aurora A/Gleevec, Aurora A/AT9283, Aurora A/Danusertib protein tended to converge after ∼50–300 ns, indicating that the simulated systems achieved sufficient stability through 500 ns of cMD simulations. Accordingly, the RMSD evolution of the heavy atoms of the Gleevec, AT9283, Danusertib in each simulated system maintained relative stability (RMSD fluctuation <2 Å) during the 500 ns simulation (Fig. 3). Interestingly, the average RMSD values of Aurora A bound with/without different ligands from cMD simulations follow the order of Apo >Gleevec >Danusertib >AT9283. Similarly, ligands that follow the order of Gleevec >AT9283 >Danusertib showed similar values. These results were in direct agreement with the experimental data that selective ligands of AT9283 and Danusertib were more stable in Aurora A than the non- selective ligand of Gleevec. Moreover, the apo Aurora A exhibited more fluctuation during the whole cMD simulation, suggesting that apo Aurora A or Aurora A bound with non-selective ligand may undergo dramatic conformational fluctuations that cannot stabilize the protein conformation. Therefore, the structural and energetic properties of simulated systems were analyzed during the last 200 ns cMD simulation trajectories.

To further investigate the effects of selective and non-selective inhibitors on the structural movements of Aurora A, principal component analysis (PCA) and dynamic cross-correlation (DCC) analysis were carried out. As plotted in Fig. 4, the correlation ranges from −1 to 1, corresponding to the extent of correlation or inverse correlation as indicated by the respective color intensities. Blue color represented positive correlation values ranging from 0.25 to 1; magenta color represented inverse correlation values ranging from −0.25 to −1; and light blue or magenta color represented weak or no-correlation for values ranging from −0.25 to +0.25. As shown in Figs. 4A–4B, the blue and magenta regions for Apo Aurora A and Aurora A/Gleevec were quite similar. A similar pattern was observed for Aurora A bound with AT9283 or Danusertib (Figs. 4C–4D). However, those regions in the Apo Aurora A and Aurora A/Gleevec were larger and more intense than those of Aurora A/AT9283 and Aurora A/Danusertib. These results present evidence that the overall correlated motions of Apo Aurora A and Aurora A/Gleevec are more dramatic than those of Aurora A/AT9283 and Aurora A/Danusertib, suggesting that the correlation or inverse correlation motions were enhanced in the Apo Aurora A and Aurora A/Gleevec. Besides, binding of Apo Aurora A or Aurora A with non-selective inhibitor may be more pronounced than binding with selective inhibitor.

To highlight the differences of Aurora A’s collective motion in different systems, PCA was performed to further characterize the conformational transformations over time. PCA was commonly utilized to reduce the dimensionality of the MD simulation data and to identify conformational transformations. In theory, the first two principal components (PC1 and PC2) from PCA capture the majority of the variance in the original distribution of molecule conformational ensembles. Herein, the conformational ensembles of the systems under study were analyzed by projecting the trajectories of PC1 and PC2 into a two-dimensional space. When PC1 and PC2 are mapped onto each other, structures with a high degree of similarity cluster together. Therefore, each cluster represents a different protein conformational state. As illustrated in Fig. 5, the simulated conformations in each system were dynamic and fluctuant during 500 ns cMD simulations, and eventually stabilizes into a dominant state. The conformational changes of the Apo Aurora A and Aurora A/Gleevec were greater than those of the Aurora A/AT9283 and Aurora A/Danusertib. This demonstrated that the conformational distributions of Aurora A bound without/with non-selective inhibitor (Apo Aurora A and Aurora A/Gleevec) were remarkably different from those bound with selective inhibitors. Then, the frequencies of PCA scatter plots were quantified and the highest-frequency structures were aligned to their initial structures for the four simulated systems as plotted in the Figs. S2 and S3. Surprisingly, these representative structures did not undergo large conformational changes (except the flexible loops, such as A-loop, Fig. S3). These results indicated that Aurora A bound with selective inhibitors had a different conformational fluctuation compared to that bound with non-selective. Apo Aurora A or Aurora A bound with non-selective ligand may undergo dramatic conformational fluctuation that cannot stabilize the protein conformation. In comparison, Aurora A bound with selective ligand could stabilize its conformation, which is congruent with the DCC analysis results.

Binding free energies calculations to determine the key residues

The binding free energy of Gleevec, AT9283 and Danusertib associated with Aurora A was calculated using Molecular Mechanics/Generalized Born Surface Area (MM/GBSA), which has been proved to be an important complement to the evaluation of protein-ligand interactions (Chen et al., 2018; Sun et al., 2014a; Sun et al., 2014b; Weng et al., 2019; Xu et al., 2013). In this study, although prediction of binding free energies based on multiple independent trajectories is more rigorous theoretically, the intramolecular energies (angle energy, bond energy etc.) cannot be well canceled by the multiple independent trajectories strategy. This might result in great uncertainty or noise of the calculated binding free energies (ΔG_binding). Therefore, a single trajectory was performed following the procedure in previous studies (Sun et al., 2014b; Xue et al., 2014). As summarized in Table 1, the predicted ΔG_binding for Gleevec, AT9283 and Danusertib bound with Aurora A were -21.39 ± 3.11, −40.47 ± 2.90, and −38.99 ± 3.82 kcal/mol. These results were consistent with both of the reported experimental conclusions. Then, a detailed description of the binding of free energy to each residue in the receptor ligand recognition pattern was analyzed. Per-residue energetic differences between Aurora A bound with Gleevec and bound with AT9283 or Danusertib are plotted in Fig. 6. Positive values suggested that the interaction with Aurora A residue bound with Gleevec was stronger than the Aurora A residue bound with AT9283 or Danusertib. Conversely, negative values suggested that the interaction with Aurora A residue bound with Gleevec was weaker than the Aurora A residue bound with AT9283 or Danusertib. As shown in Fig. 6A, the residues of Glu-211, Tyr-212 and Arg-220 interacted more significantly with AT9283 in Aurora A than in Gleevec. Similarly, the differences between Aurora A/Gleevec and Aurora A/Danusertib were quite similar. When bound to Aurora A, the residues of Glu-211, Tyr-212 and Leu-263 interacted more significantly with AT9283 than with Gleevec (Fig. 6B). Notably, the residues of Glu-211 and Tyr-212 were located in the conserved structural elements of hinge region required for ligand binding (Figs. 6C–6E). Structural analysis of Aurora A/Gleevec revealed that at a distance of 3.4 Å, the nitrogen group of Gleevec formed a hydrogen bond with the backbone carbonyl of Ala-213. In contrast, binding of AT9283 or Danusertib with Aurora A revealed a different picture (Fig. 6C). The nitrogen groups of AT9283 and Danusertib formed important hydrogen bonds with the residues of Glu-211 and Ala-213. In particular, the Glu-211 was not formed for Gleevec (Figs. 6D–6E). To further analyze the stability of these hydrogen bonds, the distance evaluations of hydrogen bonds were monitored during the whole cMD simulation process and plotted (Figs. 6F–6H). As shown in Figs. 6F–6H, the distance of hydrogen bond for Gleevec to Ala-213 were relative stable within 3.5 Å during the 500 ns cMD simulations. Similarly, these for AT9283 and Danusertib were also quite stable during the simulations. However, the distances of hydrogen bond between Gleevec and Glu-213 were more than 4 Å during the whole cMD simulations. That is to say, Gleevec cannot form stable hydrogen bond with Glu-213 during the cMD simulation and cannot stabilize the hinge region. As a comparison, the distance of hydrogen bond for AT9283 and Danusertib to Glu-211 were quite stable within 3.5 Å during the whole cMD simulations. These observations suggested that the hydrogen-bond networks for AT9283 and Danusertib bound to Aurora A were more stable than the networks bound to Gleevec, which were congruent with the DCC analysis results that the hinge region were more dramatic (Fig. 3).

GaMD simulations supported the results from cMD simualtions

To better understand the free energy map (FEM) associated with the conformational changes of Aurora A coupled with different ligands, the conformation ensemble was sampled in more detail and the possible hidden energy barriers in the cMD simulation were detected. Subsequent to the GaMD simulations, the RMSDs of the protein backbone and the heavy atoms of the ligand were calculated. As shown in Fig. 7, the RMSDs of protein backbones in all the simulated systems reached equilibrium after ∼50–100 ns of GaMD simulations. The RMSDs of the heavy atoms of ligands in each simulated system remained dynamic (Fig. 8). Similar to the results from cMD simulations, the average RMSD values of protein backbone of Aurora A from GcMD simulations follow the order of Apo >Gleevec >Danusertib >AT9283, and the heavy atoms of ligands follow the order of Gleevec >AT9283 >Danusertib. In addition, DCC analysis indicated that the correlated and anti-correlated motions were stronger in the Apo Aurora A and Aurora A/Gleevec than those in Aurora A bound with AT9283 and Danusertib, which is consistent with the results from cMD simulations (Fig. 9). The above results indicated that Aurora A, when bound with non-selective inhibitors, induces a larger conformational fluctuation and greater variability between protein subunits. In contrast, when bound with selective inhibitors, Aurora A stabilizes the conformational ensembles.

PCA was then performed to visualize the conformational changes. As shown in Fig. 10, the conformational distributions of Apo Aurora A and Aurora A/Gleevec were significantly different from Aurora A bound with AT9283 or Danusertib. Compared to the Apo Aurora A and Aurora A/Gleevec, the systems of Aurora A/AT9283 and Aurora A/Danusertib showed fewer structural clusters, suggesting that the selective inhibitors could stabilize the protein conformation. Thereafter, the FEM was performed to assess the relationship between the conformational changes and energy changes (Fig. 11). Theoretically, more energy wells (dark blue regions) indicated that the protein underwent larger conformational changes during MD simulations. As shown in Figs. 11A–11B, three and two major deep energy wells were observed for the Apo Aurora A and Aurora A/Gleevec. However, Aurora A bound with AT9283 or Danusertib were confined to a single deep energy well with high similarity, demonstrating the unstable nature of Aurora A bound with non-selective inhibitors during GaMD simulations (Figs. 11C–11D). These results are also supported by previous findings from cMD simulations, which showed that hinge region may directly regulate the conformational ensemble to ligand specifically binding.

Dissociation pathways of Gleevec, AT9283 and Danusertib from Aurora A

A growing body of research has suggested that dissociation process of a ligand from its target can reveal the dynamic mechanisms of interactions, which is vital to the development of a novel selective ligand. Even though recent advances in computing power have allowed for the use of long-time-scale cMD simulations, simulation of the dissociation pathway of a ligand from its target remains difficult. Among the various enhanced sampling methods, the umbrella sampling simulation may be the most classic and commonly used method to explore the dissociation pathway of a specific ligand from the binding pocket from its target. In this study, the dissociation process of Gleevec, AT9283 and Danusertib from Aurora A was characterized by the umbrella sampling method. In order to guarantee the sampling convergence of the umbrella sampling simulations, each simulated window was submitted to 8 ns umbrella sampling simulations, and the convergence of potential of mean force (PMF) was checked every nanosecond. As illustrated in Fig. 12, eight curves were plotted for each simulated system, and the PMFs achieved satisfactory coincidence after 7 ns umbrella sampling simulations for each window (difference of PMFs, 0.5 kcal/mol). Thereafter, the binding affinities (PMF depth, ΔW_PMF) were obtained by averaging the last 5 ns US simulations (21–25 ns). The ΔW_PMF of Aurora A/Gleevec, Aurora A/AT9283 and Aurora A/Danusertib were predicted by umbrella sampling simulations were, −12.73 ± 0.51, −17.83 ± 0.62, and −17.50 ± 0.60 kcal/mol, respectively. The ΔW_PMF can be congruent with the reported experimental data and correctly ranked. Thus, the PMF curves were utilized to clarify the different dissociation processes from the binding pocket of Gleevec, AT9283 and Danusertib toward Aurora A. The processes of the Gleevec dissociation from Aurora A were illustrated in Fig. 13 (yellow line). These results suggested that the PMF curves of AT9283 and Danusertib in Aurora A were obviously different from those of Gleevec, while the PMF values of AT9283 and Danusertib toward Aurora A were not significantly different. Initially, all the inhibitors were located at the binding sites of Aurora A (RC was 0 Å). When the inhibitors were dissociated from the binding site, the PMF values increased rapidly (RC was 0–10 Å). Thereafter, the PMFs along with the RC was on a steady increase toward equivalence (RC was 10–20 Å). During the rapid increase of the PMF values, the PMF curves of AT9283 and Danusertib showed an obvious peak and a valley. This phenomenon highlights the fact that the dissociations of AT9283 and Danusertib from Aurora A have an obvious energy barrier that needs to be overcome, and that no obvious energy barrier was observed for Gleevec. These observations indicated that the dissociation of AT9283 and Danusertib from Aurora A needs to overcome a higher free energy barrier than the Gleevec. Therefore, the AT9283 and Danusertib combined with Aurora A have a longer residence time. Further structural analysis showed that no hydrogen bond between Aurora A and Gleevec was observed when the RC were nearly 4 Å (Fig. 13B). In comparison, some important hydrogen bonds were observed for AT9283 and Danusertib at the peaks of the RCs. In addition, most of these hydrogen bonds were triggered by Glu-211, and Ala-213. These results are in agreement with previous findings based on the per-residue free energy decomposition and structural analysis, where the hydrogen-bond networks for Gleevec bound to Aurora A were not as stable as those for AT9283 and Danusertib, which means that AT9283 and Danusertib, but not Gleevec, could selectively bound to Aurora A.