The complexes of Serine/Threonine Kinase 16 (STK16) with hit molecules NPC132329, NPC160898, and the approved drug Neratinib were subjected to 500 ns MD simulations. In addition to these complexes, the complex with the highly selective STK16 inhibitor (STK16-IN-1) as a positive control was also investigated in the MD simulation studies. The stability of STK16 complexes with control, Neratinib, NPC132329, and NPC160898 was analyzed in terms of RMSD in protein backbone atoms, where lower RMSD indicates reasonable stability in the protein structure [35, 36]. In the analysis of the RMSD in ligand atoms, the backbone atoms of STK16 were used for least square fitting, and RMSD was calculated from ligand atoms, which effectively resulted in analyzing how ligand conformations remained in the binding site and consequent conformational stability of ligands [37].
Further, the RMSD in STK16 ligand atoms was also analyzed by choosing the protein backbone atoms for least square fitting, which effectively eliminated the translational and rotational motions in the protein and was sufficient to analyze the position of ligands relative to the protein. The results showed that the average RMSD in protein backbone atoms was lowest in the STK16-NPC132329 complex with an average of 0.1950 nm, while the STK16-Neratinib complex had the highest average RMSD of 0.2388 nm (Table 2, Fig. 3A). The RMSD in the protein backbone atom remained almost stable until around 300 ns in the STK16-NPC160898 complex with an average of around 0.22 nm and thereafter deviated to slightly higher RMSD, reaching around 0.3 nm. The RMSD almost remained stable in the case of the STK16-control complex with an average of 0.2277 nm. Overall, the results showed that the ligand NPC132329 resulted in the most stable conformation of STK16 compared to the control and Neratinib.
The RMSD in the atoms of control was stable and lowest with an average of 0.1150 nm, while the NPC160898 had the highest with an average of 0.4911 nm (Fig. 3B). The RMSD in NPC132329 was lower during the first 100 ns and thereafter rose to higher RMSD with an overall average of 0.3668 nm. The RMSD in Neratinib atoms remained stable throughout the simulation period with an average of 0.3252 nm. The RMSD in ligand atoms revealed that the ligands, control, and Neratinib remained intact in the binding site. Amongst the ligands NPC132329 and NPC160898, the ligand NPC160898 remained reasonably stable in the binding site.
The overall RMSF in the side chain atoms of residues was lowest for the STK16-NPC132329 complex with an average of 0.1201 nm, while highest for the STK16-Neratinib complex with an average of 0.1446 nm (Fig. 3C). A larger magnitude of fluctuations was observed in the residues in the ranges 25–30 and 225–250 in the case of the STK16-Neratinib complex. The magnitude of fluctuations was minimal in the case of STK16-control and STK16-NPC132329 complexes. The RMSF in STK16-NPC160898 was intermediate and almost similar, comparable to STK16-control and STK16-NPC132329 complexes. The RMSF results suggest that the RMSF in the STK16 bound to Neratinib was higher than all other complexes, suggesting major conformational changes at the binding site. At the same time, the RMSF in the remaining complexes with control, NPC132329 and NPC160898, was reasonably lower and only had major fluctuations in the binding site residues.
The radius of gyration, a measure of the rotation of protein structure around its centre of mass, provides insights into the compactness of the structure, where lower Rg indicates a compact and stable structure [38]. Almost all complexes showed a very stable and compact structure of STK16. The radius of gyration in STK16 was lower in the case of STK16-control and STK16-NPC132329 complexes with averages of 2.0415 and 2.0475 nm, respectively (Fig. 3D). Meanwhile, the Rg in STK16 was slightly higher for STK16-Neratinib and STK16-NPC160898 complexes with averages of 2.0619 and 2.0349 nm, respectively. However, the Rg in STK16-NPC160898 was found reasonably stable during the first 350 ns simulation, and a slightly larger magnitude of deviations was observed after 350 ns until the end of the simulation period. Overall, the control and NPC132329 stabilized the STK16 better than Neratinib and NPC160898.
The solvent-accessible surface area analysis showed an almost stable solvent-accessible surface on STK16 in all the studied complexes with averages of 161.22, 159.27, 159.59, and 160.94 nm2 for STK16-control, STK16-Neratinib, STK16-NPC132329, and STK16-NPC160898 complexes, respectively (Fig. 3E). Overall, the solvent-accessible surface area analysis further confirmed that the solvent-exposed surfaces on STK16 remained relatively uniform in all the complexes representing the solvent-accessible surface area between 150 and 170 nm2.
The stability of protein-ligand complexes can be further gauged by analyzing the extent of non-bonded interactions. The non-bonded interactions, such as hydrogen bonds between protein and ligand, help rationalize the consequent stability of protein-ligand complexes [39]. The hydrogen bond analysis showed an average of around 3 hydrogen bonds between the control and STK16 (Fig. 4A). In order to investigate the key residues involved in hydrogen bond formation, the contact frequency analysis was performed, which showed that residue Phe101 had the highest contact frequency, while the residues Pro99 and Phe100 had more than 50% contact frequency, and the residues Leu 26, Gly 27, and Met 165 had more than 25% contact frequency. The trajectories extracted at different time intervals also confirmed the key hydrogen bond with Phe101.
In the case of the STK16-Neratinib complex, Neratinib formed an average of around 3 hydrogen bonds (Fig. 4B). The residue Phe101 was found to have the highest contact frequency, while the residues Lys49, Gly104, and Thr105 had more than 75% contact frequency, and the residues Asn108, Phe100, Pro99, Leu98 had more than 25% contact frequency. The trajectory at the equilibrium state showed a hydrogen bond between Neratinib and Lys49. Other trajectories at 100, 300, 400, and 500 ns time intervals showed this hydrogen bond. The hydrogen bond with Thr105 was also found in trajectories at 100, 200, and 300 ns time intervals.
NPC132329 formed an average of around 2 hydrogen bonds in the STK16-NPC132329 complex (Fig. 5A). The contact frequency for the Phe101 was around 75%, and for Pro99, it was around 55%. The residues Met165, Ala47, Phe100, and Leu26 showed a contact frequency of more than 25%. The equilibrium state trajectory showed hydrogen bonds with Phe101 and Pro99. The hydrogen bond with Phe101 remained intact, as seen in the trajectories at 200 and 400 ns time intervals. Meanwhile, the hydrogen bond with Pro99 remained intact, as seen in the trajectories at 200- and 300-ns time intervals. In the case of the STK16-NPC160898 complex,
NPC160898 formed an average of around 3 hydrogen bonds (Fig. 5B). The contact frequency for the residue Phe101was more than 75%, while for Phe100 was more than 50%, and for the residues Lys49, Met165, Pro99, Ala47, Leu78, and Leu26 was more than 25%. The trajectory at the equilibrium state and 100 ns showed hydrogen bonds between NPC160898 and Phe101, Asn108, and Pro99. The trajectories at 200 and 300 ns showed hydrogen bonds with Lys49. The trajectories at 400 and 500 ns showed hydrogen bonds with Phe101.
Among all ligands, the common hydrogen bond was seen with the residue Phe101, which suggests the most influential role of this hydrogen bond in the final binding affinity of ligands with STK16. When the contact frequency was investigated to get further clues of key residues that establish the highest contact frequency in the subset of MD trajectories, the result revealed that the ligands control and Neratinib establish more than 75% contact frequency with the residues Phe101, Pro99, Phe100 for control and Phe101, Lys49, Gly104, and Thr105 for Neratinib. Except for the common Phe101, these ligands' key contact frequencies slightly differ. Meanwhile, in the case of ligands NPC132329 and NPC160898, the only prominent contact frequency was observed with Phe101. The flexible loop spanning from the residue numbers 98–106 and 25–34 seems to be influenced by these contact frequencies. However, the hydrogen bond interaction is one of the non-bonded interactions which are responsible for the favourable binding affinities of the ligands and a few other non-bonded interactions, such as van der Waal’s interactions and other hydrophobic interactions, might be crucial in the binding affinities of these ligands.
The principal component analysis derived Gibb’s free energy analysis showed that most of the metastable STK16-control conformations occupy the lowest energy basin between − 5 to 9 on PC1 and − 7.5 to 5 on PC2 (Fig. 6A). Meanwhile, a few metastable conformations were found in another small energy basin occupying the region between − 12 to -10 on PC1 and 1 to 4 on PC2. The representative metastable conformations from these lowest energy basins showed the hydrogen bond between the control and STK16. However, major conformational flexibility was found in two domains of STK16, which were made up of residues in the ranges 24–34 and 98–106 (shown in dark blue). The conformations from the boundary region with ΔG between 2–4 kJ/mol considerably different occupancy of the ligand as well as the orientation of the flexible loop region. In the case of the STK16-Neratinib complex, two lowest energy basins were observed (Fig. 6B). The largest energy basins occupied the region between − 4 and 1 on PC1 and − 4 and 1 on PC2. The representative metastable conformations from this region showed a hydrogen bond between Neratinib and Thr105, Pro99, and Lys49.
The smaller low energy basin occupying the region between 3 and 4.5 on PC1 and − 1 to -2 on PC2 showed only hydrogen bonds with Thr105 and quite a different orientation of the loop region. The trajectories from the region with slightly higher ΔG showed a hydrogen bond with Lys49; however, the loop region was slightly away from the ligand position.
The complex STK16-NPC132329 showed the lowest energy metastable conformations occupying the energy basin between 1 and 3 on PC1 and − 1 and 3 on PC2 (Fig. 7A). The representative metastable conformation from this energy basin showed a hydrogen bond with Phe101 and Pro99 from the loop region of residues 98–106. At the same time, the conformations from other slightly higher energy regions differ in the position of NPC132329 and the orientation of the loop region. The complex STK16-NPC160898 showed a major and more significant energy basin between 0 to 3 on PC1 and − 4 to 2 on PC2 and one smaller energy basin between − 3 to -1 on PC1 and 0 to 1 on PC2 (Fig. 7B).
The metastable conformations from the larger energy basin showed the hydrogen bonds between Phe101, Lys49, and Asp166, while the conformations from the smaller energy basin showed the hydrogen bonds between Phe101 only. It is further confirmed from the results of Gibb’s free energy analysis that the metastable conformations of STK16-control show the key hydrogen bond with Phe101. However, the energy barrier between two states of metastable conformations in this complex revealed the orientation of the control ligand in the binding pocket and the flexibility of the loop regions. The metastable conformations of Neratinib mainly existed in the conformations such that the hydrogen bonds were formed with Thr105, Pro99, and Lys49. Other states with small energy basins with metastable conformations or slightly higher energy states showed hydrogen bonds with either Thr105 or Lys49, respectively.
The metastable conformations of NPC132329 revealed that the propensity to form hydrogen bonds with Phe101 and Pro99 due to the flexibility of the loop region encompassing these residues is important. Similarly, the case of NPC160898 suggested that the propensity of hydrogen bond formation with the residues Phe101, Lys49, and Asp166 is the deciding factor. The slightly higher-energy conformations lack hydrogen bond interactions with Lys49 and Asp168.
The cluster analysis showed four unique clusters for the STK16-control complex (Fig. 8A). The conformations of cluster 1, the major cluster, originated throughout the simulation period, showing the unique position of the control ligand at the binding site. The cluster analysis for the STK16-Neratinib complex showed two unique clusters, with cluster 2 being the major cluster (Fig. 8B).The conformations of this major cluster originated after a simulation period of around 150 ns. In the case of the STK16-NPC132329 complex, 3 clusters were found, out of which cluster 1 had the maximum number of conformations (Fig. 9A).
The conformations from this major cluster originated from the simulation period 225 ns onwards. In the case of the STK16-NPC160898 complex, 6 clusters were observed (Fig. 9B).Cluster 5 had the maximum number of conformations, originating at different time intervals during the simulation. The cluster analysis revealed that the STK16-control exists in 4 unique clusters differing in the ligand position at the binding site due to the flexible loop regions, where most of the conformations exited in cluster 1. STK16-Neratinib complex showed only two unique conformations where the conformations in cluster 2 are predominant. The major conformations of STK16-NPC132329 existed majorly after around 200 ns showed unique occupancy of NPC132329 in the binding pocket. At the same time, the STK16-NPC160898 existed in six clusters, where the major cluster 5 showed the unique binding of NPC160898 with the flexibility of the loop regions.
The MM-GBSA calculations wherein the enthalpic contributions are considered while calculating the binding free energies (ΔG). The relative binding free energies were found most favorable for Neratinib with ΔTotal of -49.17 kcal/mol. Neratinib's lowest relative binding free energy was due to the lowest and favorable van der Waals energy of -47.66. However, Neratinib had a higher entropic energy (-TΔS) of 12.56 kcal/mol, which resulted in the binding free energy (ΔG) of -36.61 kcal/mol. Comparably, the control had the ΔTotal of -33.92 kcal/mol, and the ΔTotal of NPC160898 was almost close to the control with ΔTotal of -32.85 kcal/mol.
Both control and NPC160898 had almost similar van der Waals energies, and NPC160898 had a slightly higher entropic energy of 7.41 kcal/mol. After accounting for the entropic energies, the binding free energies for control and NPC160898 were − 29.81 and − 25.44 kcal/mol, respectively. The NPC132329 had a ΔTotal of -28.68 due to a slightly higher van der Waals energy of -28.99 kcal/mol. After accounting for the entropic energy, the binding free energy for NPC132329 was − 18.61 kcal/mol, which was slightly higher than that of other compounds.
The MM-GBSA analysis revealed that Neratinib had the most favourable binding free energy, considering the enthalpic energy. The ligand NPC160898 has binding free energy comparable to control, while the binding free energy of NPC32329 is slightly higher. In summary, the ligands NPC132329 and NPC16098 showed promising interactions compared to either the control ligand or Neratinib.
Table 2
MM-PBSA calculations for STK16 complexes
Energy component (kcal/mol)
|
Averages for the STK complexes
|
Control
|
Neratinib
|
NPC132329
|
NPC160898
|
ΔVDWAALS
|
-32.68
|
-47.66
|
-28.99
|
-32.02
|
ΔEEL
|
-3.88
|
-35.67
|
-1.98
|
-4.16
|
ΔEGB
|
6.75
|
40.83
|
5.60
|
7.53
|
ΔESURF
|
-4.12
|
-6.66
|
-3.32
|
-4.21
|
ΔGGAS
|
-36.56
|
-83.34
|
-30.97
|
-36.18
|
ΔGSOLV
|
2.64
|
34.16
|
2.29
|
3.32
|
ΔTOTAL (ΔH)
|
-33.92
|
-49.17
|
-28.68
|
-32.85
|
-TΔS
|
4.11
|
12.56
|
10.07
|
7.41
|
ΔG
|
-29.81
|
-36.61
|
-18.61
|
-25.44
|
ΔVDWAALS: van der Waals energy; ΔEEL: Electrostatic energies; ΔEGB: Polar solvation energy; ΔESURF: Nonpolar solvation energy; ΔGGAS = ΔVDWAALS + ΔEEL; ΔGSOLV = ΔEPB + ΔENPOLAR; ΔTOTAL = ΔGSOLV + ΔGGAS; ΔG = ΔH-TΔS |