The high-resolution K-Ras G12C mutant protein structure (PDB ID: 6P8Y) that was inhibitor co-crystallized was retrieved from Protein Data Bank [23, 24]. The retrieved protein structure was in the inactive GDP form with 2.31 Å and an Rfree value of 0.269. The mutations found in the PDB structure of G12C were reverted to the actual amino acids (AA) (C51S; C80L; C118S) with the help of the residue mutate module within the Schrödinger. The G12C residue of the K-Ras protein was mutated to G12D using the module mentioned above, and the protein model was engaged for the K-Ras G12D mutation-based module. The Schrödinger suite’s protein preparation wizard was employed to prepare the protein structure [25]. The proteins were subjected to tautomerization and ionization state adjustments, as well as the inclusion of missing residues, hydrogens, formal charge correction, and bond order correction at pH 7. The protein structure was further processed where the water molecules were removed; the structural integrity was mended using the protein preparation module and OPLS3e FF, respectively. While retaining the heteroatoms, the zero-order states were assigned to the metal, and the conformations of the AA side chains were evaluated. In K-Ras, the amino acids of histidine have a different shape. The K-Ras structure was optimized, the protonation states and tautomers were also predicted. The structure was also refined by improving the H-bond and eliminating the steric impediments by utilizing constrained heavy atom minimizations with a root mean square deviation (RMSD) value of 0.3.

For the E-pharmacophore modeling investigation, seven inhibitors were identified from the literature survey and used as templates. The 3D structures of sotorasib (AMG510), ARS1620, Indole lead 1${}^{\text{Amgen}}$, 4_am, bortezomib, and MRTX849 [26, 27, 28, 29, 30] were retrieved from the PubChem database. The molecules were prefiltered based on the Rule of five (ROF) to exclude compounds with no drug-likeness properties [31]. The LigPrep module [32] was used to develop the ligand dataset files, and the OPLS3e force field was then used to minimize energy further. The ligand molecules were reconstituted to their 3D geometries using Epik, the hydrogens were added at an ionization pH of 7.0 $\pm$ 2.0, and the optimum stereochemical conformation for each ligand was calculated. Further, the ligand sets were obtained from the Drug Bank database [33, 34]. The phase database was created using the 3D atomic structures of 2672 FDA-authorized compounds. The LigPrep and Epik modules were initially used to prepare the ligands. The Confgen (OPLS3e FF) was then used to generate the ligands’ chiralities and conformers (50 conformers). The molecules underwent an initiatory screening based on ROF to weed out any false-positive compound. Further, ADME attributes of the molecules were estimated, and maestro files were generated.

The initial docking of ligands to the G12C and G12D proteins is crucial for E-pharmacophore-based modeling, which combines ligand- and structure-based approaches. The Grid Ligand Docking and Energetics (GLIDE) module from the Schrödinger interface performed docking studies on the optimized protein structure. A 3D grid of 20 Å was created around the active sites of the protein structures to accommodate the binding residues before starting the docking process [35, 36, 37, 38]. The seven ligands were docked into the active site of the protein structure using GLIDE XP mode. The PHASE module of the Schrödinger suite was used with the docked files to create an E-pharmacophore model. This module uses six built-in pharmacophore features for model generation: hydrogen bond acceptor (A), hydrogen bond donor (D), hydrophobic (H), negative ionizable (N), positive ionizable (P), and aromatic ring (R). Pharmacophoric sites are created around atoms contributing to GLIDE XP energies during hypothesis development. These sites are ranked by calculating the energy terms of atoms. The hypothesis is used as a query to search for hits in the phase database that are more potent than the reference ligand.

A three-step virtual screening method was performed on the phase database using GLIDE software v12.5, New York, NY, USA, a conventional molecular docking tool that has successfully screened large databases during drug development [39]. The virtual screening process includes High Throughput Virtual Screening (HTVS), Standard Precision (SP), and XP docking. This approach is the first step in the docking process to reduce the number of compounds used for precision docking investigations, as HTVS helps in rapid screening. The G12C and G12D hypotheses were used to screen the phase database for suitable compounds. These compounds were then subjected to HTVS, followed by SP and rigorous XP docking to obtain better binding affinities and eliminate false positives [36, 40]. Based on their binding energies, the docking results were compared to those of the reference compound (sotorasib) [41]. Additionally, the pharmacokinetic properties of the compounds that passed the filtering procedures were evaluated.

One way to decrease the failure rates during drug development is to focus on molecular characteristics such as absorption, distribution, metabolism, excretion, and toxicity (ADMET) [42]. The current study employed the QikProp module to forecast the ADME properties and drug-likeness (DL) [43]. QikProp evaluates molecular properties using a pharmacokinetic descriptor range that covers 95% of existing drugs. To better understand the molecule’s potential as an inhibitor, sotorasib was compared to the predictions generated for other hit compounds.

MDS of nine K-Ras mutant complexes were performed using GROMACS software for 1800 ns [44]. The complexes were prepared with drugs such as sotorasib, afatinib, osimertinib, hydroxychloroquine, sotorasib, doxapram for G12C and fluphenazine, hydroxyzine, and zuclopenthixol for G12D. The force field parameters for the hit molecules were calculated using Charmm27 [45, 46], and the variables and topological files were generated using SwissParam [47]. A dodecahedral box with simple point-charge (SPC) waters was positioned at 1.0 nm from the box’s edge for the mutants complexed with hit molecules, and sotorasib and counterions were added to neutralize the system. Energy minimization was performed using the steepest descent method for 5000 steps to remove steric conflicts. The volume and pressure ensembles were set to default values. The molecular dynamics (MD) trajectories of each complex were analyzed using GROMACS and visualized using ggplot2 in R studio [48]. Principal component analysis (PCA) was used to obtain the spatial structure of the mutant K-Ras inhibitor complexes and sotorasib [49]. The PC1 and PC2 projections and eigenvectors were calculated using the default approach in GROMACS and our previous simulation protocol [50, 51, 52, 53]. The covariance matrices of all complexes were evaluated, and the conformational space was computed using the atom positions from the MD trajectories. The 3D free energy landscape (FEL) was plotted using OriginPro 2022 SR1 software, Northampton, MA, USA, and its projections on eigenvectors 1 and 2 were evaluated and plotted using R studio.

The PHASE module in the Maestro workspace generated an E-pharmacophore hypothesis based on seven pre-existing inhibitor compounds: sotorasib (AMG510), ARS1620, ARS1323, ARS853, bortezomib, MRTX849, and ARS1630. The pharmacophore model was developed using docked poses of these inhibitors and incorporates structural and energy information between G12C, G12D, and the ligands. The model generation aims to identify new active sites for developing novel inhibitors. The seven ligands were docked using XP, and their descriptor files served as input for generating the model. The “E-pharmacophore” option was defined using the docking post-processing module in the Maestro interface. To map the ligands and pharmacophoric sites into the ligand’s atoms, the fragment mode option was used during model building. The pharmacophoric sites retrieved were fixed to two at the start, and a pharmacophore hypothesis was generated in which the features were prioritized according to their energy contributions to binding. For G12C, the top five attributes were selected for the hypothesis, while for G12D, the top three attributes were taken based on an energy score of $>$1.0 kcal/mol (Table 1). The pharmacophoric model generated for G12C featured three aromatic groups (R), one hydrogen bond donor, and one acceptor group, respectively (Fig. 1A), while G12D contained two aromatic rings (R) and one hydrogen bond acceptor group (Fig. 1B). The phase database was then screened using the generated E-pharmacophore models, which identified many active compounds, indicating that the database was effectively screened. This advanced scientific edit provides a detailed overview of how the E-pharmacophore hypothesis was generated and its potential use in identifying new active sites for developing novel inhibitors.

Fig. 1.

The pharmacophore hypotheses were represented using the PHASE module, which included features from Sotorasib (AMG510), MRTX849, ARS1620, Bortezomib, ARS853, ARS1630, and ARS1323. The hypothesis for G12C (with A4, D7, R15, R16, R17) is shown in (A), while the hypothesis for G12D (with A4, R14, R16) is shown in (B). The acceptor groups are represented in blue with arrows indicating h-bond acceptors. The donor groups are represented in pink with arrows indicating h-bond donors. The aromatic rings are represented in orange.

The E-pharmacophore model (DDRRR) was used to screen the phase database and identify a large group of potential inhibitors for G12C and G12D. During phase database construction, Ligand conformers were generated using the ConfGen option, so on-the-fly conformer generation was unnecessary. The ligand conformers were then compared to the E-pharmacophore model to identify those that matched its features. Reorienting the ligand conformers also helped to determine if they matched the pharmacophoric sites. The basic criteria for ligand site matching with the E-pharmacophore model were set to 4 out of 5 sites to ensure no potential molecule was missed. During database screening, the “hit treatment” option was set to “Return at most 1000 hit molecules”. PHASE identified 110 hit molecules for G12C and 189 for G12D that could superimpose the hypothesis when the process was completed. All 299 hit compounds were then subjected to virtual screening processes.

In the first screening, 299 hit compounds were identified out of 2672 FDA-approved compounds. These compounds were expected to have properties from the implied E-pharmacophore and were labeled as possible hit compounds if all the properties were satisfied. After the initial PHASE filtering, docking with the Glide module was used. The molecules with higher glide scores and energies than the reference inhibitor Sotorasib (G12C: –4.648 kcal/mol, G12D: –3.449 kcal/mol) were screened. The FDA-approved inhibitors (sotorasib-AMG510) were used as a reference drug for both mutants. In HTVS, 110 and 189 compounds were docked with the mutant models (G12C and G12D). After HTVS docking, the compounds were subjected to SP docking, based on the Glide G-score ($>$3.0 kcal/mol and 4.0 kcal/mol), 67 and 128 hit compounds, respectively, for G12C and G12D were selected for the XP docking module. As a result, 29 and 113 compounds were obtained as hits from XP docking. The final hits (G12C: 3, G12D: 4) that obtained a higher G-score and MM-GBSA binding energies than the reference compound (sotorasib) (Table 2) were selected for ADME analysis.

The Schrodinger’s (v12.5, New York, NY, USA) QikProp module was used to evaluate the ADME properties of hit compounds for the G12C (Afatinib, Osimertinib, Hydroxychloroquine) (Supplementary Fig. 1A–C) and G12D (Hydroxyzine, Zuclopenthixol, Fluphenazine, Doxapram) (Supplementary Fig. 2A–D) models. The CNS descriptors were a key focus of the study, as they can determine a molecule’s ability to cross the blood-brain barrier. The QikProp results for all hits and their characteristics for the K-Ras mutations are listed in Supplementary Tables 1,2. Compounds that met the description while remaining within acceptable bounds were selected as potential inhibitors. The pharmacokinetic properties of the hit compounds were evaluated for drug-like behavior and found to be within acceptable limits (Table 3). The screened compounds with favorable predicted properties were examined as a group, and the results indicate that the hit compounds satisfied all the ADME descriptors.

The interaction of the compounds with mutant models was visualized using Schrödinger’s LID and Discovery Studio to gain further insights into the binding patterns. Figs. 2,3 show the hit molecules for mutant models in 3D, highlighting their interactions. The grey surface model for G12C (Fig. 2A–D) depicts the binding pocket, while the hit molecules are represented as sotorasib (violet), afatinib (green), osimertinib (blue), and hydroxychloroquine (red). In G12D (Fig. 3A–E), the hit molecules are represented as sotorasib (violet), hydroxyzine (green), zuclopenthixol (yellow), fluphenazine (red), and doxapram (blue) in a ball and stick model. The interaction mechanism of the G12C-sotorasib complex in 2D revealed three H-bonds with carbon atoms C12, E62, and R68, Pi-alkyl with V9, and Pi-Pi stacking at Y96 (Supplementary Fig. 3A,B). The G12D-sotorasib complex had one H-bond with R68, a Pi-cation with K88, and two Pi-Pi stacks at Y96 (Supplementary Fig. 3C,D). The binding patterns of top compounds were compared with the G12C-sotorasib model. The G12C-afatinib complex had two H-bonds with G10 and D69, a Pi-Pi stack with Y96, and a Pi-alkyl with V9, M72, and Y96 (Fig. 4A,B). For G12C-osimertinib, it was observed that five H-bonds were formed with G10, C12, T58, R68, and Y96, with two Pi-Pi stacking with H95 and Y96 (Fig. 4C,D). For G12C-hydroxychloroquine, it was observed that two H-bonds formed with T58 and D92, two Pi-Pi stacking at V9 and H95, and a halogen bond with T58 (Fig. 4E,F). Finally, the binding patterns of top-hit molecules were compared with the reference compound for G12D. For the G12D-hydroxyzine complex, two H-bonds formed with A59 and Q61 and Pi-alkyl with V7, V9, and M72 (Fig. 5A,B). Three H-bonds were formed for G12D-zuclopenthixol with K16, D12, and D69; a Pi-cation with R68; and three Pi-alkyls with R68, H95, and Y96 (Fig. 5C,D). For the G12D-fluphenazine complex, it was observed that a single H-bond formed with G10 and a Pi-PI stacked with Y96 (Fig. 5E,F). A single H-bond formation was observed for the G12D-doxapram complex with K16, a single Pi-Pi stacking at Y96 (Fig. 5G,H).

Fig. 2.

The 3D representation of the G12C K-Ras mutant complexes’ docked poses with hit compounds. The complexes include: (A) G12C-Sotorasib (Purple). (B) G12C-Afatinib (Green). (C) G12C-Osimertinib (Blue). (D) G12C-Hydroxychloroquine (Red). The hit compounds are represented using a ball and stick model, with Sotorasib in purple, Afatinib in green, Osimertinib in blue, and Hydroxychloroquine in red. The oxygen and nitrogen atoms of the hit compounds are shown as red and blue balls, respectively, while the chloride and fluoride ions within the hit compounds are depicted as green and light green balls. The interacting amino acid residues are shown as green sticks, with blue, yellow, pink, and green dashes indicating Pi-Pi stacked hydrogen bonds, attractive charges, and Pi-Cation interactions. The names of the interacting residues are written in light green fonts. The binding pocket is shown as a grey surface model.

Fig. 3.

The 3D representation of the G12D K-Ras mutant complexes’ docked poses with hit compounds. The complexes include: (A) G12D-Sotorasib. (B) G12D-Hydroxyzine. (C) G12D-Zuclopenthixol. (D) G12D-Fluphenazine. (E) G12D-Doxapram. The hit compounds are displayed in a ball and stick model, with Sotorasib in violet, Hydroxyzine in green, Zuclopenthixol in yellow, Fluphenazine in red, and Doxapram in blue. The oxygen and nitrogen atoms are shown in red and blue ball models. The chloride and fluoride ions are depicted in green and light green ball models. The sulfur atom is shown in a yellow ball model. The interacting amino acid residues are displayed in a green stick model. The blue, yellow, pink, and green dashes represent Pi-Pi stacked hydrogen bonds, attractive charges, and Pi-Cation interactions. The names of the interacting residues are written in light green. The binding pocket is shown as a grey surface model.

Fig. 4.

The image depicts multiple interactions between the ligand and the G12C K-Ras protein. The diagrams (A), (C), and (E) show the ligand interactions as visualized by the Maestro suite, while (B), (D), and (F) show the 2D ligand interactions as visualized by Discovery Studio. The complexes shown in (A) and (B) are G12C-Afatinib, those in (C) and (D) are G12C-Osimertinib, and those in (E) and (F) are G12C-Hydroxychloroquine. In the Maestro suite’s ligand interaction diagrams, the red, green, and purple residue interactions represent Pi-Cation, Pi-Pi stacked/T-shaped, and hydrogen bonds, respectively.

Fig. 5.

The image depicts multiple interactions between a ligand and the G12D K-Ras protein. The diagrams (A), (C), (E), and (G) are from the Maestro suite, while (B), (D), (F), and (H) are 2D ligand interactions from Discovery Studio. The G12D-Hydroxyzine, G12D-Zuclopenthixol, G12D-Fluphenazine, and G12D-Doxapram complexes are represented by (A) and (B), (C) and (D), (E) and (F), and (G) and (H) respectively. In the Maestro’s ligand interaction graphs, the red, green, and purple residue interactions represent Pi-Cation, Pi-Pi stacked/T-shaped, and hydrogen bonds, respectively.

The MDS of 200 ns were performed on each complex to examine the stability of the docked hit molecules compared to sotorasib. The Charmm27 force field was used for the simulations. The structural stability of the nine complexes was assessed using the RMSD. The G12C-afatinib and G12C-hydroxychloroquine complexes showed more deviations (more than 0.2 nm) (Supplementary Fig. 4A), while the G12D-zuclopenthixol complex showed the most fluctuations (0.25 nm) (Supplementary Fig. 4B). The G12C-osimertinib and G12D-hydroxyzine complexes exhibited the least deviations compared to the G12C and G12D-sotorasib complexes. The intermolecular hydrogen bonds were calculated for all nine complexes to determine which bonds contribute to stable complex formation. All three inhibitors in G12C had a greater number of hydrogen bonds than sotorasib (1–2), with a total of 5–6 bonds in each case (Supplementary Fig. 4C). In G12D, all four inhibitors showed a moderate number of hydrogen bonds (1–4) in comparison to sotorasib (Supplementary Fig. 4D). The root mean square fluctuations (RMSFs) values were calculated for each of the complexes using ggplot2. In G12C, the highest fluctuations were observed in G12C-hydroxychloroquine (0.38 nm), followed by G12C-afatinib (~0.35 nm) (Fig. 6A). In contrast, in G12D, the highest fluctuations were seen in G12D-doxapram (~0.34 nm), followed by G12D-fluphenazine (~0.29 nm) in comparison to sotorasib and other inhibitors (Fig. 6B). The compact nature of protein was assessed by visualizing the radius of gyration (Rg) from the MD trajectory. The osimertinib complex showed a steady graph in comparison to the sotorasib complex (Fig. 6C), while the G12D-doxapram, zuclopenthixol, and hydroxyzine complexes all displayed a consistent Rg trend in comparison to the sotorasib complex (Fig. 6D). Across the simulation process, the mean Rg for all nine complexes was 1.53–1.59 nm.

Fig. 6.

The plots represent the RMSF and Rg of the G12C and G12D models complexed with reference (Sotorasib) and shortlisted inhibitor compounds. The first plot (A) shows the RMSF of the G12C model complexed with Sotorasib, Osimertinib, Afatinib, and Hydroxychloroquine. The second plot (B) shows the RMSF of the G12D model complexed with Sotorasib, Doxapram, Fluphenazine, Hydroxyzine, and Zuclopenthixol. The third plot (C) shows the Rg of the G12C model complexed with Sotorasib, Afatinib, Osimertinib, and Hydroxychloroquine. The fourth plot (D) shows the Rg of the G12D model complexed with Sotorasib, Doxapram, Fluphenazine, Hydroxyzine, and Zuclopenthixol.

The movements of the C atoms in each of the nine complexes were analyzed using PCA, as the atomic dynamics of the macromolecule are collective. The PCA plot was created using the first and last PCs, representing the models’ basic aggregate motion. When compared to sotorasib, osimertinib showed less separation of motion in the G12C mutant (Fig. 7A), while hydroxyzine and zuclopenthixol showed less separation of motion in G12D (Fig. 7B). The 3D Gibbs FEL plots show a lower dimensionality in the G12C-osimertinib complex and G12D-zuclopenthixol compared to other inhibitors and reference compounds (Supplementary Figs. 5,6). The stability of the complexes is represented by the volume and shape of the lowest energy region, shown in blue.

Fig. 7.

The figure shows the PC1 and PC2 contact maps of the G12C and G12D models complexed with inhibitor molecules, illustrating the C$\alpha$ motions relative to the reference molecule Sotorasib. Essentially, these plots depict the alterations in the protein structure of the G12C and G12D models upon complexation with inhibitor molecules, relative to their structure when complexed with Sotorasib. The initial plot (A) illustrates the contact maps for PC1 and PC2 of the G12C model when complexed with Sotorasib, Afatinib, Osimertinib, and Hydroxychloroquine. The subsequent plot (B) displays the contact maps for PC1 and PC2 of the G12D model when complexed with Sotorasib, Doxapram, Fluphenazine, Hydroxyzine, and Zuclopenthixol.