Identification of thiazolidin-4-one scaffold-based molecules as PI3Kα inhibitors: Molecular docking, MM/GBSA, IFD, ADMET, water-map analysis, and MD simulations

2.1. Ligand selection and screening of a database

A library of thiazolidine-4-one containing molecules with 17,729 compounds was exported to structure-data format (SDF) format and downloaded from Chemdiv. The downloaded SDF file was imported into the Maestro 2023 interface of the Schrodinger suite and was further processed by the Lig-prep tool. Energy optimization of all the compounds was carried out by giving an ionization target pH of 7.4 ± 0 using the Epik application. Desalting was defined along with determining 3D structure chirality and generating several ligands at 32 under the force field OPLS4.

2.2. Protein preparation and grid generation

An extensive literature review was conducted, by which it was revealed that the α isoform of PI3K will serve as a potential target for the treatment of various types of cancers. An X-ray crystal structure of PI3Kα with a cocrystal ligand NVP-BYL719 (PDB id: 4JPS). The Phosphatidylinositol-4,5-bisphosphate-3-kinase catalytic subunit alpha isoform consisted of chain A with a sequence length of 1,074 and two mutations in the gene PI3KCA. The PI3K regulatory subunit alpha consisted of chain B with a sequence length of 293 and no mutations in the genes. The structure had a resolution of 2.2 Å and an observed R-value of 0.206 and was downloaded from the protein data bank in PDB file format and imported to the Maestro interface of the Schrodinger suite.

Following the import of protein structure to Maestro, it was reviewed, refined, and modified using the protein preparation wizard. The structure was checked for missing loops and sidechains. The prime module helped to fill the missing loops and side chains. The regulatory subunit of the protein, which consisted of chain B, was deleted. The bond orders were assigned, hydrogen atoms were added, waters beyond 5 Å were deleted, and the energy minimization of the protein was carried out in PROPKA at a pH of 7.4 using the OPLS4 force field.

The prepared protein was further subjected to grid generation using the Glide mode in the Schrodinger software. A grid was generated, taking the co-crystallized ligand alpelisib as the reference point. Furthermore, this site was confirmed using the site map analysis tool present in the Schrodinger software. As both these methods resulted in the same sites, it was further used for the docking studies.

2.3. Screening by molecular docking

To overcome the limitations of traditional methods of searching for therapeutically active molecules, the use of computational tools emerged. By using computational techniques, the cost and time required for the process were reduced, and at the same time, large libraries of molecules could be easily screened virtually against a particular biological target or multiple targets [12]. Virtual screening involves the use of simple techniques like molecular docking, as well as advanced statistical and artificial intelligence-driven machine learning techniques [13]. Here, we use virtual screening of a large compound library through various molecular docking protocols and select suitable hits through it.

The prepared ligands were screened against PI3Kα using Glide Ligand Docking of the Maestro, Schrodinger suite. The process of screening and filtering the hits was done by systematically carrying out high-throughput virtual screening (HTVS) initially, followed by standard docking (SP) protocol and extra precision (XP) docking protocol. Initially selected 17,729 ligands were screened against the prepared protein by HTVS. Ligands with considerable docking scores were selected and further screened by SP and XP docking.

2.4. Prime molecular mechanics generalized born surface area analysis

The top 10 molecules after XP docking were further taken for prime-molecular mechanics energies combined with generalized Born and surface area (MM-GBSA) analysis. It is a useful method for determining the binding free energy involved during the binding of small molecules with macromolecules such as proteins [14]. This involves a static approach to predict the stability of protein-ligand interactions [15,16]. During these calculations, the contribution of individual energies such as that of Hydrogen bond, Vander Waals interactions, Electrostatic interactions, and solvation energies is taken into account and is calculated based on the following equations [14]:

ΔG_bind = ΔE_MM + Δ G_GB + ΔG_SA

ΔE_MM = E_ele + E_vdw + E_Hbond

• _·E_MM – molecular mechanical energy

• G_GB – polar contribution towards solvation energy

• G_SA – contribution towards solvation energy

• E_ele – electrostatic energy

• E_vdw – vanderwaals energy

• E_Hbond – Hbond energy

It was carried out using the prime module of the Maestro interface. Ligands and receptors were selected from the project table and workspace entries, respectively. The solvation model was chosen as variable-dielectric surface-generalized born by default settings, and it was run on the OPLS4 force field by selecting minimize as the sampling method.

2.5. Absorption, distribution, metabolism, excretion, and toxicity (ADMET) predictions

The Qik-prop tool of the Maestro interface was utilised to predict the drug likeliness of the top 10 compounds by calculation of various physicochemical and pharmacokinetic properties. Important rules, such as Lipinski’s rule of 5, were checked for the 10 compounds to understand their drug likeliness.

2.6. Induced fit docking

Since the time of recognition of molecular docking as an important technique in drug design and determination of drug-target interactions, there have been many advancements to increase the efficiency of the technique [17]. Induced fit docking (IFD) overcomes the limitations of rigid docking protocols, which are commonly employed in normal docking procedures [18]. There occurs a huge number of rearrangements occur frequently in the ligand binding site of protein, along with conformational changes in the binding site as well as to the whole protein structure [17]. Considering all the conformational changes and incorporating them into the docking protocol is a computationally expensive and time-consuming process [18]. The induced fit approach was developed to overcome these challenges and was accepted as a promising approach due to its ability to cover various conformational spaces of the protein and to generate various possible orientations. IFD takes into consideration the flexibility of ligands as well as proteins and characterises binding site interactions [19]. Molecules obtained after MM-GBSA analysis and drug likeliness studies are subjected to the IFD protocol of the Maestro interface. The top ligands, along with the standard drug, were chosen from the project table, and the box centre was chosen as the centroid of the workspace ligand. Standard protocol was used, which generates up to 20 poses. In prime refinement, residues were refined within 5 Å of the ligand poses.

2.7. Molecular dynamics simulations

With the use of the molecular dynamic (MD) simulation technique, the drug discovery process and the refining process of the hit molecules, to get the best promising ligand out of it, can be accelerated. Classical MD simulation involves the application of Newton’s laws that come under classical mechanics to the movement of atoms and molecules [19]. In the MD simulation process, the millions of atoms in a system are allowed to move in discrete time steps, and during each step, laws of classical mechanics are applied with the help of a force field to calculate the energy of these atoms and the overall energy of the system. Various force fields are developed, which generally consider various non-bonded interactions of these atoms that contribute to the energy of the system. With the help of this simulation technique, the time-dependent behaviour of microscopic systems can be studied easily [19,20]. It finds use in the study of the time-dependent behaviour of the protein-ligand complex and its stability throughout the simulation time.

The top molecules of IFD results were subjected to MD simulation studies with the minimized PI3Kα protein structure using the Desmond application. Initially, for the simulation, the system was built in an OPLS4 Force field using a predefined simple point charge solvent model, and boundary conditions were set with an orthorhombic box shape. The system was neutralized by adding Na⁺ ions. After system building, energy minimization was carried out, and the simulation was run for 100 ns with NPT conditions at 300 K and 1.01325 bar pressure at an approximate number of frames set to be 1,000.

2.8. Water map analysis and advanced MD simulation study

The water map application of Schrodinger was used for refining and selecting the best compound. During water map calculation, the free energies associated with the water molecules that occupy various hydration sites on the apoprotein are calculated. The default scoring function of the water map was utilised, which calculates binding free energies associated with each ligand. When ligand molecules bind to the protein in its binding site, it displaces several water molecules. The free energies associated with this displacement of water molecules are summated to get the binding free energy of each ligand, upon which ligands can be scored [21].

In this study binding site was defined by selecting the atoms of the ligand that are bound to the protein. The water molecules around 10 Å of the binding site were set for analysis. The water molecules associated with the protein were deleted, and a short simulation of 2 ns was run in OPLS4 Force field in the Water map. From the results obtained regarding each hydration site and its associated overlap, occupancy, and free energies of each hydration site was analysed. Also, weighted Monte Carlo/molecular mechanics (WM/MM) scoring was run to calculate the total binding free energies of the ligand, based on which ligands could be scored.

After the water map calculation, the top molecule, which showed the most promising results, was selected and was subjected to an MD simulation of 500 ns run with default settings and the same solvent system as done previously in the OPLS4 force field [22].

3.1. Virtual screening through molecular docking technique

From the results obtained, all molecules that showed considerably good docking scores (>7.0kcal/mol) were 1,428 in number, which were screened again using SP docking. After SP docking, 1,159 molecules were obtained with docking scores of more than 7.0 kcal/mol. These molecules were further screened by XP docking, and from the results top 10 molecules were chosen. Similarly, XP docking for the prepared standard drug Alpelisib was also done. The XP docking results of the top 10 molecules selected according to their docking score are shown in Table 1.

Table 1. List of the top ten compounds after XP docking, as per their docking score, their respective structure and MMGBSA ΔG scores.

Table 2. 2D interaction diagrams of top 10 molecules (selected from XP docking result) with PI3Kα protein (4JPS) along with the interaction diagram for Alpelisib- PI3Kα complex.

The docking scores and ligand protein interactions of the top 10 molecules were analysed. By analysing ligand interaction diagrams of those molecules (Table 2) and of the standard drug Alpelisib, the common interaction was noted to be that of protein amino acid residue Val851 forming hydrogen bonds with one or more ligand atoms. In top molecules—KPK2, KPK3, KPK5, KPK6, KPK7, KPK9, and KPK10, VAL851 formed two hydrogen bonds with carbonyl oxygen and nitrogen atom of the Thiazolidin-4-one nucleus, respectively. In Alpelisib, VAL851 was only shown to be interacting with the Nitrogen heteroatom of the Thiazole ring. This could be one of the reasons for better docking scores for the ligands with Thiazolidin-4-one moiety over the Thiazole moiety of Alpelisib. But it is also to be noted that the molecule KPK1 with the highest docking score does not show significant interaction with VAL851.

Table 3. Top 5 compounds with their respective IFD scores after IFD.

The other common interaction among the ligands was found to be the π-π interaction of residue TYR836 of the protein with the phenyl rings of the ligands. These interactions may further increase the binding efficiency, as they will result in stable ligand-protein complexes. π-π stacking with TYR836 was observed in KPK2, KPK3, KPK5, KPK6, KPK8, and KPK10. The Sulphur atom of the thiazolidine-4-one ring was exposed to solvent in most of the ligands. It was noted that KPK5 showed a maximum number of hydrogen bonds with various residues of the protein structure. It showed seven hydrogen bonds with residues such as ASP933, ASP810, and LYS802. These residues formed a Hydrogen bond with a water molecule in the protein, and the Nitrogen atom of the thiazolidine ring exhibited H-bond interaction with the water molecule. Along with it, VAL851 and TYR836 also showed H-bond interactions with the ligand. These observations suggest better stability of the KPK5 ligand–PI3Kα protein complex.

3.2. Prime molecular mechanics generalized born surface area analysis

By performing Prime molecular mechanics generalized born surface area (MMGBSA) analysis of the top 10 compounds, selected from XP docking results, the top 5 compounds with good values were selected. The molecule KPK9 showed the best MMGBSA ΔG bind value of −64.15, followed by KPK5, KPK6, KPK2, and KPK8 with values −63.14, −61.65, −58.04, and −57.87, respectively. The standard drug alpelisib showed a value of −10.24. These top 5 molecules were further taken for IFD analysis.

3.3. IFD analysis

After the MMGBSA analysis, the top five compounds that were selected were subjected to IFD against the protein, the results of which are given in Table 3.

Table 4. ADMET analysis of top two molecules and standard drug.

From the IFD scores, the two compounds, KPK5 and KPK3, show better scores than the standard drug Alpelisib. Both KPK5 and KPK3 showed 13 and 18 binding poses, respectively, with differing IFD scores, and among them, the top 3 poses of both KPK5 and KPK3, according to the IFD scores, with their 2D interaction diagrams are shown in Figures 1–3.

Figure 1. IFD results of KPK5-PI3Kα complex- A, B, C represents top three binding poses. [Click here to view]

Figure 2. IFD results of KPK3-PI3Kα complex- A, B, C represents top three binding poses. [Click here to view]

Figure 3. IFD results of Alpelisib-PI3Kα complex- A, B, C represents top three poses. [Click here to view]

From the interaction diagrams for the top three poses of the KPK5-PI3Kα complex, the interactions with the following residues: TYR836, ASP933, and LYS802 were maintained in all three poses. Three of these residues showed hydrogen bond interaction with the ligand. Only the third pose showed H-bond interaction with VAL851 and π-π stacking between TYR836 and the aromatic ring.

But in the case of the KPK3-PI3Kα complex, all three poses showed H-bond interaction of ligand atoms with Val851 and π-π stacking between TYR836 and the aromatic ring. SER854 and ASN853 also showed H-bond interactions with the ligand atoms in all three poses. Therefore, in the case of KPK3, interactions were all maintained throughout the top three poses.

Looking at the interaction diagrams of the top three poses of the Alpelisib-PI3Kα complex, the VAL851 H-bond interaction with the ligand was only present in the second and third poses, not in the top pose. H-bond interaction of GLN859 with the thiazole ring of the ligand was observed in all three poses of the complex. ARG770 also showed interaction with the ligand in the first and third binding poses.

It is seen that the IFD scores for all poses of KP8, KPK9, and KPK6 come under the IFD scores for Alpelisib. Hence, only the top two molecules, i.e., KPK5 and KPK3, were selected for running the MD simulation with the PI3Kα protein. Before the MD simulation process, the two compounds, along with alpelisib, were checked for their drug likeness and ADMET properties.

3.4. ADMET analysis and drug likeness prediction

The results of ADMET analysis of the top 2 compounds, obtained after IFD, were done and the results are shown in Table 4.

Table 5. Water map results of the top 2 molecules and the standard drug.

According to various molecular descriptors shown in Table 4, the ADMET properties of the two molecules can be roughly approximated. The molecular weights of the two compounds were found to be below 500 daltons, which is ideal for drug likeness. The QPlog_o/w (Partition coefficient of the molecule in an octanol–water mixture) for both molecules was found to be in the acceptable range of −2.0 to 6.5. QPlogS value shows the predicted aqueous solubility, and ideally, it ranges from −6.5 to 0.5. The values for the two molecules KPK5 and KPK3 were found to fall within this range. The predicted apparent Caco-2 cell permeability in nm/s was found to be more than 25, which shows that the molecules have a considerable permeability. The predicted brain/blood partition coefficient also falls within the recommended values, –3.0 to 1.2. Human oral absorption of 1, 2, or 3 signifies low, medium, or high absorption, respectively. Both KPK5 and KPK3 were found to have high human oral absorption. The QPlogHERG values show the possibility of the compound inhibiting HERG channels. It represents the potential for cardiotoxicity, and if the value goes below −5.0, it is of concern. The molecule KPK5 shows a QPlogHERG value slightly less than −5.0. However, it was noted that the standard drug Alpelisib showed a QPlogHERG value of −2.124, which shows a much worse condition than the slight possibility of cardiotoxicity exhibited by KPK5. Thus, in terms of cardiotoxicity, KPK5 is comparatively better and less cardiotoxic than the standard drug Alpelisib. None of the compounds violated Lipinski’s rule of five, which shows that both compounds have druglike properties.

3.5. MD simulations

MD simulation was performed for the top 3 molecules obtained after IFD analysis, which are: KPK5, KPK3, and KPK4, and for the standard drug Alpelisib. All simulations were carried out for 100 ns, and the results obtained were analysed.

Initially, the simulation results of the Alpelisib- PI3Kα complex were analysed. In the protein-ligand random mean square deviation (RMSD) (P-L RMSD) (Fig. 6A), it was observed that the ligand RMSD increases from 2. 4 Å and reaches values of 6.4 Å in the initial 0–10 ns period, then stabilizes and fluctuates around 4.8–6.4 Å. On nearing 40 ns, a sudden hike of ligand RMSD to 7.2 Å is observed, which again comes down and fluctuates around 4.8–6.4 Å. The protein RMSD, on the other hand, fluctuates in a range of 1.8 Å (at the beginning of the simulation) to 3.8 Å (at the end of the simulation). High RMSD for ligand and protein around the first 0–10 ns period shows decreased stability of the protein-ligand complex during that time. The protein random mean square fluctuation (RMSF) (Fig. 6B) also shows several fluctuations or peaks in the graph, indicating the alterations in the protein structure due to local changes in the protein structure throughout the simulation time. These fluctuations indicate constant changes in the protein stability pattern throughout the simulation, affecting the biological efficiency. While analysing the ligand-protein contacts (Fig. 6C and D), the crucial interaction of the ligand with Val851 was only sustained for a short period of simulation time and was not retained throughout the simulation. On the other hand, there were significant interactions observed during the whole simulation time between ligands and residues such as ARG770, LYS802, TRP780, and GLY859. Contact with LYS 802 and ARG770 was maintained by the ligand quite stably throughout the simulation, as shown protein-ligand contact graph.

Figure 4. MD simulation data for KPK5-PI3Kα complex. [A] Protein-Ligand RMSD. [B] Protein RMSF. [C] Protein-Ligand Contacts. [D] Ligand-Protein contacts. [Click here to view]

Figure 5. MD simulation data for KPK3- PI3Kα complex. [A] Protein-Ligand RMSD. [B] Protein RMSF. [C] Protein-Ligand Contacts. [D] Ligand-Protein contacts. [Click here to view]

Figure 6. MD simulation data for Alpelisib- PI3Kα complex. [A] P-L RMSD. [B] Protein RMSF. [C] Protein-Ligand Contacts. [D] Ligand-Protein contacts. [Click here to view]

In the P-L RMSD of the KPK5–PI3Kα complex (Fig. 4A), the protein as well as the ligand RMSD fluctuates to a high value in the initial nanoseconds of simulation. The protein RMSD increases from 1.5 Å to about 3.5 Å during the 0–10 ns period and then stabilizes. The protein RMSD fluctuations thereafter remain quite stable and are observed in the range of 3–3.5 Å. The ligand RMSD also increases rapidly from 2.4 to 7 Å and comes down immediately to a stable range of 2.4 Å. This initial increase and decrease of protein and ligand RMSD during the first few nanoseconds of the simulation shows that the protein-ligand complex attains stability over time. The ligand RMSD fluctuates in the range of 3.2–4.2 Å for a period of 60 ns, and again it comes down to 2.4 Å, increasing and fluctuating in the range of 3.2–4.8 Å. During the 100 ns simulation time, the ligand RMSD never went above 5 Å. This shows that the range at which RMSD fluctuates for KPK5 is less than that of the standard drug Alpelisib, which indicates better stability of the KPK–PI3Kα complex over the Alpelisib-PI3Kα complex. The protein RMSF gives the areas of protein that fluctuate during the simulation time. It was observed that the fluctuations observed were low compared to those of standard drugs (Fig. 4B). This shows that the protein remained quite stable comparatively through the simulation time. The protein-ligand contacts (Fig. 4C and D) were analysed, and it was observed that the crucial interaction of the ligand and protein residue—VAL851 was maintained for more than 90% of the simulation time. The residue showed H-bond interaction for more than 60% of the simulation time and then hydrophobic (π-π stacking) and water bridge interactions, respectively. This water bridge interaction is an H-bond between the residue and the ligand, mediated by a water molecule. It also showed significant interactions of TYR836 and the ligand for more than 70% of the simulation time. A total of five interactions were seen between various residues of protein and ligand that sustained for more than 70% of simulation time, which shows good interaction between the ligand and protein and thereby better stability of the protein-ligand complex.

In the KPK3–PI3Kα complex, the RMSD of simulation (Fig. 5A) shows that ligand RMSD fluctuates between 5 and 7 Å, and initially a rapid increase is observed, like the other ligands, from 1 to 7 Å, and then almost stabilized fluctuations between 5 and 7 Å. These fluctuations are much more than those of KPK5. The protein stability during this simulation is better than the protein stability that was observed during the simulation with alpelisib. The RMSF plot (Fig. 5B) of KPK3 with protein showed comparatively fewer fluctuations than that of alpelisib. In the protein-ligand contacts, three main residues were found to play a significant role, which were GLU259, SER275, and PHE794 (salt bridges in the case of the first two and hydrophobic in the case of the third residue, respectively). The interaction with VAL851 and TYR836 was not shown by this ligand.

From the analysis of the simulation results of the two ligands and that of the standard drug with the protein, it was seen that KPK5 shows much more efficiency in the stability of protein and protein-ligand complex, and shows significant protein-ligand contacts. For further confirmation, the water map analysis results are used to refine the two ligands and the standard drug, and select the best potential inhibitor of the protein structure.

3.6. Water map analysis

By performing water map analysis, the hydration sites were identified in the region around the ligand binding site. These hydration sites can be classified into three based on the free energies associated with the site: low-energy (ΔG ≤ 1.5 kcal/mol), medium-energy (1.5 < ΔG < 3.5 kcal/mol), or high-energy (ΔG ≥ 3.5 kcal/mol). The simulation was run without retaining the ligand, so the hydration sites shown will be those sites that exist in the protein in the absence of the ligand [22]. Those sites that are overlapped by ligand atoms will have an overlap factor >0.1. An overlap factor of 1 is associated with a complete overlap of the hydration sites with ligand atoms. Hydration sites with highly negative free energy (i.e., highly negative dg values) will have water molecules that are highly stable at that site, which are not replaceable or displaceable as they are tightly bound to the protein. The hydration sites at medium energy and low energy will have comparatively more easily replaceable water molecules, with the lower energy ones being highly displaceable. As the number of sites with low energies and overlapping factors close to one increase, the easiness of replacing the water at that site and binding with the protein increases for the ligand.

In the water map analysis of the KPK5–PI3Kα complex (Fig. 7A), it was observed that out of all observed hydration sites, nine sites on the protein completely overlapped with that of ligand, with an overlap factor of one. Out of those nine sites, site 26, site 45, and site 71 showed very high positive dg values of 5.33, 3.88, and 4.77, respectively, indicating that they are very easily replaceable. There were also five medium-energy sites and one low-energy hydration site of the protein that completely overlapped with that of the ligand. By observing the 2D interaction diagram associated with the analysis, an unstable hydration site is identified in the regions where significant interaction between ligand and Val851 residue of protein happens; there is a complete overlap of the ligand with that site. The presence of highly replaceable water molecules at this site signifies the easiness with which the ligand replaces this water and forms hydrogen bonds with the protein, which leads to better stability and increased biological efficiency of the protein-ligand complex. There is another replaceable water molecule observed in site 104, where overlap with the ligand is complete, which is involved in mediating the other important hydrogen bond interaction between the ligand and TYR836.

Figure 7. Water map results of hit compounds. [Click here to view]

The water map analysis of the KPK3–PI3Kα complex (Fig. 7B) showed nine sites on the protein structure that completely overlapped with the ligand, and three out of the nine: site20, site33, and site 72 were with high positive free energies of 5.25, 3.97, and 4.93, respectively. The benzene ring that is fused with the piperido–pyrrolidine ring is involved in π-π stacking interaction with the residue TYR836. This benzene ring can easily displace a highly unstable water molecule at hydration site 63, as observed in the interaction images, and this ease helps it bind and establish interaction with the protein efficiently. Also, there is hydration site 20 with easily displaceable water molecules, which overlap with the carbonyl oxygen atom of the ligand. The replacement of water molecules in this site would result in the formation of the important Val851-Ligand interaction.

In the case of the standard drug Alpelisib–PI3Kα complex (Fig. 7C), the hydration sites that showed complete overlap with the ligand were 11 in number, and out of these 11 sites, 3 sites were having highly positive free energies—site 13, site 26, and site 47 with 5.17, 5.17, and 3.99 kcal/mol, respectively. It can be seen that site 26 overlaps the nitrogen atom of the thiazolidine ring in the ligand, and the displacement of water molecules at this hydration site by the ligand will be easy due to the high positive free energy of water molecules at that site. The ligand displaces this water and results in the formation of an H-bond interaction between Val851 and the nitrogen atom of the thiazolidine ring. Therefore, this unstable water molecule at this site helps the ligand to establish the significant Val851-ligand interaction easily.

By seeing the hydration site, its free energies, and overlap factor, filtering of the ligands is difficult as all three ligands showed similarity in those properties. For this, WM/MM scoring for water map results associated with each ligand protein complex was done. Based on the results obtained from the scoring, ligands were selected based on their MM-GBSA ΔG, water map ΔG, and WM/MM scores. The more negative the scores, the better their binding efficiency to the protein. The results are shown in Table 5. KPK5 and KPK3 both have shown good WM/MM ΔG bind scores compared with the standard drug Alpelisib. The water map entropy value for the compound KPK5 is −23.706 Kcal/mol, which is higher compared to both the standard drug as well as the KPK3 compound. From the results, it is also evident that the MMGBSA ΔG Values are higher than the standard for the compounds KPK3 and KPK5. Though the values such as water map ΔG and water map enthalpy are noticed to be less than the standard, the values are not so less to reject them. Therefore, by comparing the results, it can be clearly noted that compound KPK5 shows a good water map analysis result compared to KPK3 and can be used for further analysis.

3.7. Advanced simulation study

The advanced simulation time helps to further determine the stability and duration of that stability for the ligand-protein complex. Longer duration of stability of the complex shows the binding efficiency and inhibitory potential of the ligand [23]. The data regarding P-L RMSD, protein RMSF, and protein-ligand contacts during the 500 ns simulation of the KPK5–PI3Kα complex and Alpelisib- PI3Kα complex are depicted in Figures 8 and 9. It can be observed that the P-L RMSD for KPK5–PI3Kα complex (Fig. 8A) remains in an acceptable range for most of the simulation time but shows fluctuations beyond 5.6 at three points between 200 and 300 ns period. No fluctuation is observed beyond 6.4. But in the case of the Alpelisib-PI3K complex, the P-L RMSD (Fig. 9A) shows fluctuations beyond 6.4 multiple times, and a maximum deviation of 7.2 is observed. The average RMSD value lies in the range of 4.8–6.4 for the standard drug, while for KPK5, this value lies in the range of 2.4–4.8 during most of the simulation time. This difference in the P-L RMSD values of the standard drug and KPK5 shows the superiority of KPK5 over Alpelisib in making a stable and long-lasting complex, which is necessary for eliciting the required biological action.

Figure 8. 500 ns MD simulation data for KPK5-PI3Kα complex. [A] Protein-Ligand RMSD. [B] Protein RMSF. [C] Protein-Ligand Contacts. [D] Ligand-Protein contacts. [Click here to view]

By observing the P-RMSF of both KPK5-PI3Kα (Fig. 8B) and Alpelisib-PI3Kα complex (Fig. 9B), in the latter, the RMSF deviates beyond 6.4 for three residues, while in the case of KPK5, it is observed only for one residue. Also, it can be observed that the fluctuations are comparatively less in the case of KPK5 than in Alpelisib, when they bind with the protein. The fewer the fluctuations, the better the stability. So, these two observations suggest that our ligand KPK5 causes a less conformational change in the protein PI3Kα than the standard drug Alpelisib when they are complexed with it.

Figure 9. 500 ns MD simulation data for Alpelisib-PI3Kα complex. [A] Protein-Ligand RMSD. [B] Protein RMSF. [C] Protein-Ligand Contacts. [D] Ligand-Protein contacts. [Click here to view]

The protein-ligand contacts throughout the 500 ns simulation for the KPK5-PI3Kα complex (Fig. 8C and D) show an interesting fact. Throughout the study, VAL851-ligand interaction was assumed to be a significant one, as it was present in almost all ligands that showed good binding efficiency against the protein. Also, in the 100 ns simulation study, it was observed that the ligand-Val851 interaction was sustained for more than 70% of the simulation time, suggesting its importance in maintaining the stability of the ligand-protein complex. However, as observed in Figure 8C, the interactions that are mostly sustained throughout the 500 ns simulation are TYR836-Ligand interaction and ASP 810-Ligand interaction. Earlier also during the docking studies and the 100 ns MD simulation, the importance of TYR836 was noted and was an important interaction that may play a crucial role in facilitating required biological action along with Val851. Both were observed in the standard drug-protein complex. But by this longer simulation time, it is revealed that TYR836 is the residue that sustains the stability of the complex for a longer duration and, hence plays a more crucial role than Val851 in sustaining the inhibitory action of the drug against target. TYR836 shows H-bond interaction with the ligand almost towards the end of simulation time and shows hydrophobic interaction towards the end of the simulation. ASP810 also shows H-bond interaction with the ligand during most of the simulation time, and salt bridge formation with the ligand towards the end of the simulation time.

In the case of the Alpelisib-PI3Kα complex, only one residue-ARG 770 showed a significant interaction (i.e., hydrogen bond and water bridge interaction) with the ligand till the end of simulation time. The other interactions include the hydrophobic interaction of TRP780 and the interaction of GLN859, which was a hydrogen bond initially and was later converted to a salt bridge interaction. They too remained for less than 70% of the simulation time. From Figures 8D and 9D, it is evident that when two residue ligand interactions in the case of KPK5–PI3Kα complex, TYR836 and ASP810 remained for 99% and 96% of the 500 ns simulation time, respectively, for Alpelisib-PI3Kα complex no interaction between protein residue and ligand was sustained for more than 60% of simulation time. This again shows the superiority of the ligand KP5 over Alpelisib in creating a much more stable and long-lasting protein-ligand complex, which signifies its superiority in the inhibitory action of the protein PI3Kα over Alpelisib.