Sotorasib Ras inhibitor

Keywords: anticancer, cancer, G13D, inhibitors, K-Ras, metastasis, pharmacophore model

ABSTRACT

The mutations at the hotspot region of K-Ras result in the progression of cancer types. Our study aimed to explore the small molecule inhibitors against the G13D mutant K-Ras model with anti-cancerous activity from Food and Drug Administration (FDA)-approved drug compounds. We implemented several computational strategies such as pharmacophore-based virtual screening, molecular docking, absorption, distribution, metabolism and excretion features, and molecular simulation to ensure the identified hit compounds have potential efficacy against G13D K-Ras. We found that the FDA-approved compounds, namely azelastine, dihydrocodeine, paroxetine, and tramadol, are potential candidates to inhibit the action of G13D mutant K-Ras. All four compounds exhibited similar binding patterns of Sotorasib (AMG510), and a structural binding mechanism with significant hydrophobic contacts. The descriptor features from the QikProp of all four compounds are within allowable limits compared to sotorasib drug. Consequently, a molecular simulation result emphasized that the dihydrocodeine and tramadol exhibited less fluctuation, minimal basin, significant hydrogen bonds, and potent inhibition against G13D K-Ras. As a result, the current research identifies prospective K-Ras inhibitors that could be further improved with biochemical analysis for precision medicine against K-Ras-driven cancers.

1. INTRODUCTION

The rat sarcoma virus (RAS) GTPases H-Ras, N-Ras, and K-Ras are the critical regulators of EGFR-mediated proliferation signals. Mutations in K-Ras result in an EGFR-independent activation of the proliferation signals. The most common mutation sites of K-Ras-driven cancers are codons 12, 13, and 61. These mutations frequently inhibit GTP hydrolysis due to the loss of vital structural conformations critical for GTP hydrolysis. The functional implications of K-Ras mutations vary, with codon 12 variations having a higher oncogenic and transforming potential than codon 13 variants. The most common codon 13 mutation in K-Ras is G13D and has been reported as the third most common mutation in colorectal cancer (approximately 25%).

Because direct targeting of K-Ras was technically challenging, efforts were concentrated on targeting upstream and downstream mediators of the K-Ras signaling cascade. Current fundamental and applied research breakthroughs have used mutation-specific approaches to explore and tackle individual K-Ras variants. Novel therapies targeting G13D and other K-Ras mutants have been developed and are currently being reviewed clinically. Food and Drug Administration (FDA)-approved drugs Adagrasib (MRTX849) and Sotorasib (AMG 510) are available for treating K-Ras G12C mutants specifically. JDQ4433, a drug targeting G12C, and siG12D-LODER targeting G12D are now in phase 1/2 clinical trials. The mRNA-derived K-Ras targeted vaccine V941, also known as mRNA-5671, is in Phase 1 clinical development targets G12C, G12D, G12V, and G13D driver mutations.

The development of mutation-specific drugs for K-Ras has proven challenging. However, it could be achieved by adding chemical groups that can interact efficiently with the specific mutant features. For example, patients with K-Ras mutations responded poorly to the anti-EGFR drugs panitumumab and cetuximab. The development of cmp4-based drugs (pan-Ras inhibitors) is promising for combination therapies and cetuximab in reducing tumor cell proliferation. Our previous study explored the importance of K-Ras binding site water molecules and dynamical changes of K-Ras mutant structures compared to the native model. Another study focused on developing an engineered chimeric toxin, a pan-Ras biologic inhibitor for use in cancer therapy. A recent study identified several FDA-approved inhibitors targeting G12C and G12D K-Ras models using pharmacophore modeling and a drug-repurpose approach.

Several research studies have utilized contemporary computational techniques to identify the small molecule inhibitors and structural mechanisms that help experimental researchers to test the drugs in vivo and in vitro, for targeted therapy. We integrated pharmacophore modeling with many such approaches for mutation-specific K-Ras drug development to determine the compounds targeting K-Ras G13D. This was done by generating pharmacophore hypotheses utilizing the recognized inhibitory compounds to generate a model with all the characteristics needed for the prospective drug. The pharmacophore hypothesis was then utilized in the PHASE database screening to find the probable compounds exhibiting K-Ras G13D inhibitory activity.

2. MATERIALS AND METHODS

2.1 Ligand Preparation

Our study employed 13 K-Ras inhibitors’ structural data from the literature survey as a template for ligand-based pharmacophore modeling. The chemical structures of 7 compounds were retrieved, namely ARS853, sotorasib, ARS1620, ARS1630, bortezomib, ARS1323, and MRTX849 from PubChem. The other inhibitors (2E07, Indol lead 1Amgen, Compound 12, 4_am, 6H05, and Compound 9) required for the present study were sketched (2D) with the help of a sketcher platform. Using default modes with the LigPrep component and OPLS3E force field minimization, the structural data of the inhibitors were used for ligand preparation. The inhibitor sotorasib was used as the reference compound. The ligand was prepared by transitioning the ligands into their 3D geometries and adding hydrogens (pH of 7.0 ± 2.0) using Epik. Following this, the best stereochemical structure computed per ligand was retained. Supporting Information: Figure S1 depicts the workflow of the current study.

2.2 Phase Database Construction

A phase database was created by retrieving coordinates (3D) of 2752 FDA-approved compounds in a suitable format (SDF). The DrugBank database was used to obtain the ligand sets for library construction. The ligand structures were preprocessed during the database construction to ensure the precise bond order, protonation, and tautomeric states. The LigPrep and Epik modules were used to prepare each ligand with chiralities, and 50 conformers were produced for every compound. Furthermore, using the QikProp module, all of the compounds’ absorption, distribution, metabolism and excretion (ADME) properties were created, and also considered Lipinski’s Rule of 5 to eliminate false-positive drug-like molecules. Finally, the compounds were taken in maestro format with respective features of ADME.

2.3 K-Ras Protein Preparation

The PDB structure of G13D mutant (K-Ras) crystallized with an inhibitor was downloaded from the Protein Data Bank (PDB ID: 6P8Y). The mutations such as C80L, C51S, and C118S in the K-Ras structure were restored to native residues. The G13 position was mutated to G13D using the same method, and the protein structure was subsequently employed for K-Ras G13D mutant-based research. Compared to other 3D structures, we chose the 6P8Y, which was used in several computational studies and the missing residues are not present in the 3D structure. In addition, the Protein Preparation Wizard was used to prepare the protein and restore its integrity by employing the OPLS3E force field. The protein was prepared by adding hydrogen atoms and removing water molecules above 5 Å. This was followed by correcting inappropriate bond orders, creating disulfide bonds, and generating tautomers at pH 7. The metal ion was configured to zero-order states, retaining the heteroatoms. The configuration of the amino acid side chains was assessed. The structure was optimized by predicting the protonation of histidine as K-Ras histidine exhibits alternate conformation. The hydrogen bond network was further optimized, and steric hindrances were removed. The Ramachandran plot was used to evaluate the integrity of the prepared protein structure before subjecting them to molecular docking analysis.

2.4 Ligand-Based Pharmacophore Model Generation

Our study opted to generate a ligand-based pharmacophore model using the 13 inhibitors LigPrep optimized. The PHASE module of Maestro was used for developing the common pharmacophore hypothesis (CPH). As the pharmacophore models are developed from a set of active molecules exhibiting a common framework, we classified the inhibitors into three categories. The inhibitor compounds were classed as active, inactive, or intermediate based on their experimental and clinical outcomes. The optimized geometries were matched, but the feature descriptions were left unchanged to determine the inhibitors’ optimal alignment and common pharmacophore features (CPFs). PHASE used six built-in pharmacophore properties to generate the pharmacophore models: a hydrogen bond acceptor (A), hydrophobe (H), and aromatic ring (R). By providing PHASE with user-defined values, at least half of the active set of compounds should be matched by pharmacophores and extract all association of hypotheses, retaining 4-5 CPFs. The Phase HypoScore, survivability score, and BEDROC score were used to assess and rank all combinations of common pharmacophore. Eventually, a pharmacophore hypothesis was created for the following analysis to describe characteristics with equivalent spatial configurations across all ligands.

2.5 Molecular Docking Simulations

The GLIDE module of Schrödinger was used to conduct grid-based flexible-ligand docking to analyze the interaction between K-Ras G13D mutant receptor and sotorasib. The molecular recognition interactions of the K-Ras G13D-sotorasib complexes were assessed using interaction diagrams (2D). The active site’s size and position were determined as per the ligand present in crystallized protein, and a grid was positioned 10 Å. While retaining the docking at their ranges, grid files were used for the XP docking approach to check the binding modalities for the G13D mutant receptor. The resulting values were computed as XP GScores and served as the standard in this study.

2.6 Screening of Potential Candidates and MM/GBSA Calculations

We incorporated a three-step virtual screening procedure to identify novel K-Ras G13D inhibitors. An HTVS procedure followed two precision docking techniques: glide SP and XP docking. The HTVS approach ensured that the total molecules submitted for precision protocols were reduced to a minimum. The binding effectiveness of 50% of the top-scoring molecules from HTVS docking was investigated further. Consequently, 50% of the compounds that cleared the prerequisites of SP were taken to XP protocol to exclude false positives from the screened molecules. Eventually, the compounds from XP were evaluated using the molecular mechanics with generalised born and surface area solvation (MM-GBSA) approach to predict their free energies (ΔGbind) during binding with the help of the Prime module. The pharmacokinetic properties of the substances passed through the filtration procedures were studied in greater depth.

2.7 ADME Profiling

The pharmacokinetic properties of the hit compounds were evaluated using the QikProp of Schrödinger. Distinct ADME descriptors were employed to develop more likely compounds to produce adequate ADME performance during clinical trials. QikProp provides comparison ranges for assessing the properties of the compounds to those of 95% of existing drugs. QPlogBB and central nervous system (CNS) were calculated to assess the compound’s ability to cross the BBB. QPlogPo/w quantified the lipophilicity of the lead compounds. QPlogS measured the aqueous solubility of the compounds. QPlogPw and QPlogPoct quantitatively assessed the free energy solvation in water and octanol, respectively. QPlogKp predicted skin permeability. Solvent-accessible surface area and hydrophobic component of the SASA calculated the solvent-accessible surface area and a hydrophobic component of the solvent-accessible surface area of the lead compounds. These ADME descriptors to screen compounds helped eliminate unnecessary compounds that would fail subsequent experimental procedures.

2.8 Molecular Simulations and Trajectory Analysis

The molecular simulations of 200 ns were carried out using the GROMACS package to examine the binding efficiency of the hit compounds (Azelastine, Dihydrocodeine, Paroxetine, and Tramadol) against the K-Ras G13D mutant. The forcefield parameters were computed using the Charmm27. SwissParam built the parameter and topological data for all hit compounds using the Merck molecular force field, allowing Charmm27 to detect molecules in GROMACS during the simulation process. The output files generated by SwissParam were appropriate for GROMACS topology. To solvate the G13D complexed with sotorasib and identified compounds, the SPC waters were situated at least 1.0 nm from the dodecahedral box’s edge. Sodium and chloride ions were added to neutralize the system. The complexes were subjected to a 5000-step steepest descent energy minimization to eliminate any steric clashes. All the complexes underwent 50,000 ps run at 300 K for NPT and NVT ensembles. GROMACS modules were used to explore the trajectories after the simulation, and each complex was explored individually. Principal component analysis (PCA) was carried out using GROMACS to analyze the atomic motion patterns. A standard procedure was used to calculate the first PC1 and PC2 projections. The covariance matrix was examined, and the dynamic motion of the protein was computed. The QtGrace was utilized to plot the graphs.

3. RESULTS AND DISCUSSION

3.1 Pharmacophore Generation

The pharmacophore model was created by utilizing 13 existing training set molecules encompassing active and moderately active drug molecules in the PHASE module of the Maestro workspace. The recent evidence that was suppressing the K-Ras mutation (G13D) and drugs studied in trials for G12C and G12D, we chose the active drugs, and some classified as intermediates (ARS1323, ARS853, and ARS1630), while other compounds were deemed inactive (6H05, compound 9, 2E07, and compound 12). As a result, only the active chemicals involved in the production of CPFs are considered. To eliminate speculations that could not tell the difference between inactive and active molecules, actives and inactives were used. As a result, the active chemicals entangled in the CPFs formation are considered. The hypothesis characteristics and the anticipated values from the intermediates, actives, and inactives are listed in Table 1.

Four pharmacophore hypotheses were generated (AHHR_4, AHHR_5, AHHR_1, and AHRR_2) with survival scores of 4.759, 4.689, 4.829, and 4.737, respectively. Consequently, AHHR_4 was chosen as the hypothesis for further study; it includes four features (one aromatic group, one hydrogen bond donor group, and two hydrogen bond acceptors). The superimposition of active molecules with the selected AHHR_4 hypothesis is shown in Figure 1B–F. The generated pharmacophore model was then utilized to screen the phase database, which revealed many active chemicals, showing that the database had been screened efficiently.

3.2 Database Screening with Pharmacophore Hypothesis

The PHASE database searched putative inhibitors for G13D K-Ras mutant using the CPFs identified from the actives (AHHR_4). The properties of the hypotheses were matched to a ligand group from DrugBank. The PHASE database was created by utilizing 2752 FDA-approved molecules, with 50 conformers for each molecule. The absent/incorrect residues were rectified during G13D structure preparation. Consequently, we were able to identify 1611 hit compounds from a total of 2752. The evaluated compounds were projected to contain four properties from the proposed CPFs with no partial entries. If they satisfied all requirements, they were classified as prospective hit molecules. Subsequently, after the initial database screening, the GLIDE module from the Schrödinger software was utilized for further docking. Sotorasib (AMG510), an FDA-approved inhibitor drug, was kept as the reference for the G13D; after docking of sotorasib to the allosteric site of G13D, the XP GScore was estimated to be −4.441 kcal/mole.

HTVS, standard precision (SP) docking, and XP docking are the three systematic stages of the virtual screening process. These three docking precisions performed virtual screening for the G13D mutant model. In the HTVS process, the G13D model was docked to 616 compounds. Further, the top 318 molecules were selected for another docking method (Glide SP) according to their binding score, and GScore should be greater than 3.0 kcal/mol. We found 176 hits for the G13D model according to their docking and glide GScore. At last, the 176 hits were docked using the precision module (XP docking). Consequently, 176 hit molecules were obtained after the XP screening process and had elevated XP GScores than the G13D-sotorasib complex.

FIGURE1 Selected pharmacophore hypothesis and superimposition of actives onto a hypothesis. (A) AHHR_4 hypothesis; (B) Sotorasib onto AHHR_4; (C) 4_am onto AHHR_4; (D) ARS1620 onto AHHR_4; (E) Indole lead IAMGEN onto AHHR_4; (F) Bortezomib onto AHHR_4. Actives are shown in gray stick models.

3.3 Binding Free Energy Calculation

The Prime MM-GBSA was utilized to evaluate the free energy of the G13D mutant model with the top hit molecules obtained from the XP docking process. The MM-GBSA calculations were performed on the docked G13D mutant complexes, and the hits were shortlisted from the XP docking for the G13D model. Here, we considered the ΔGbind and GScores from the XP, the hits that scored lower energies were ruled out. Consequently, the four selected molecules had more significant binding affinities than the reference molecule (sotorasib). The MM-GBSA, ΔGbind score with sotorasib for G13D was estimated at −17.92 kcal/mol.

3.4 ADME Property Analysis

The K-Ras mutations have been linked to the metastases of the brain, intracranial disease, and distant brain failure, CNF descriptors were one of the primary characteristics investigated in our study. The ADME features for the hit molecules of the G13D mutant model (Azelastine, Dihydrocodeine, Paroxetine, and Tramadol) were determined using the Schrödinger software’s QikProp module. The hit molecules that met the description’s requirements within tolerable limits were considered potential inhibitor compounds. The study employed properties like stars, CNS, percent HOA, and HOA as essential standards for drug-like substances. The ADME criteria principally evaluated for the ADME evaluation were QPPCaco (25 poor to > 500 great), QPlogBB (−3 to 1.2), CNS (−2 to +2), and HOA (1–3). The hit compounds were more effective and met the ADME features in the G13D mutant model compared to sotorasib. From the results of XP GScores, ADME features, and MM-GBSA from the compounds, we found four compounds that satisfy all the descriptor values.

3.5 Interaction Between Ligand-G13D Complex

Interactions of the hit molecules were explored to learn more about the unique binding tactics with the K-Ras G13D model by employing Schrödinger’s LID and Discovery Studio. The binding postures of the hit molecules for the G13D mutant model are shown in Figure 3. In the surface model (color) G13D (Figure 3A–E), the reference molecule is represented as sotorasib (magenta), and the remaining hit molecules are represented in the ball and stick (blue) model with interactions. The G13D-sotorasib interaction was investigated; it resulted in the formation of one hydrogen bond at A59, one halogen (fluorine) at A11, and four Pi-alkyl bonds at K16, Y71, M72, and Y96 respectively.

For the G13D-azelastine complex, it was observed that two Pi-cations with R68 and Y96, one Pi-sigma with T58, three Pi-alkyl with V9, Y71, and M72, and one amide-Pi stacking with A59 were formed. For G13D-dihydrocodeine, one hydrogen bond with Y96, one Pi-Alkyl with Y96, two attractive charges with E37, E62, and two Pi-Pi T-shaped with Y96. For G13D-paroxetine, one hydrogen bond was observed to have one unfavorable donor-donor interaction with K16, G10, one halogen interaction (fluorine) at R68, two Pi-alkyl at R68 and Y96. For G13D-tramadol, two hydrogen bonds formed with G10, Y96, and four Pi-alkyl with V7, V9, and M72. Also, at R68, an unfavorable donor-donor interaction was identified.

Inhibitors and druggable sites for K-Ras have recently been identified in several investigations, notably the switch I/II pocket between switches I and II. Consequently, switch I/II interacting residues were identified in our interaction study for G13D complexed with azelastine (T58) and dihydrocodeine (E37 and E62).

3.6 Molecular Dynamics Simulation

The stability of the docked hit molecules with G13D in contrast to sotorasib was assessed using MD trajectories obtained from 200 ns MD simulations for each complex; a total of 1000 ns of MD simulations were performed using the Charmm27 force field. The backbone root mean square deviation was employed to determine the five complexes’ structural stability. G13D-paroxetine complex showed the highest number of fluctuations (approximately 0.33 nm), followed by azelastine (approximately 0.3 nm). In comparison to sotorasib, dihydrocodeine and tramadol showed a lesser number of deviations. It can be observed that the complexes showed an equilibrated state after 160 ns.

By examining the radius of gyration the protein compactness was quantified (Rg). We obtained a stable graph for the G13D-dihydrocodeine, G13D-paroxetine, and G13D-tramadol compared to the sotorasib complex. The study revealed that G13D-dihydrocodeine and G13D-tramadol showed a lower Rg value (approximately 1.51–1.56 nm). The five complexes showed average Rg variations of 1.51–1.61 nm. According to the Rg findings, earlier complexes were more compact and had an improved binding capacity.

The RMSF’s of Cα atoms for all five complexes were plotted using the QtGrace. The highest fluctuations were observed in G13D-paroxetine and G13D-sotorasib (approximately 0.4 nm), followed by G13D-azelastine and G13-dihydrocodeine (approximately 0.35 nm), particularly in the regions such as switch I/II. The G13D-tramadol complex exhibited a comparatively minimal deviation in comparison to other complexes. The minimal deviations can be correlated to the inactivation of structure, which can be caused due to a compound binding to the allosteric site of K-Ras, ultimately favoring the inhibition of G13D when bound to hit molecules.

To gain a deeper insight into the enzyme-ligand interaction, the hydrogen bonds that help form stable complexes were calculated. The hydrogen bond plot of the G13D mutant model showed that a higher number of hydrogen bonds were seen in sotorasib and dihydrocodeine (approximately 1–4). The remaining hit molecules showed an average (approximately 1–3) hydrogen bonds. As the protein functions through cumulative atomic motions, the PCA study was employed to understand the movements of the Cα atoms structurally for the five complexes. Based on the cumulative motion of the complexes, the first and last principal components were employed to construct the PCA graph. It was observed that sotorasib, dihydrocodeine, and paroxetine showed comparatively fewer motion separations in the G13D inhibitor complexes.

The matrix (covariance) values of G13D-sotorasib, azelastine, dihydrocodeine, paroxetine, and tramadol complexes are provided in Table 4. The stability of the complexes is indicated by the blue color representing the size and the form of the minimum energy region. The Gibbs free energy landscape (FEL) was used to calculate the conformational states of complex structures, with the first two major PCA components serving as reaction coordinates, where the ΔG ranged from 0 to 16.5 kJ/mol respectively. The FEL plot illustrated that G13D-azelastine and G13D-dihydrocodeine have a minimal basin compared to other inhibitors and the reference compound (sotorasib), indicating an increase in conformation stability.

4. CONCLUSION

Our study presented the potential inhibitor candidates against the K-Ras G13D mutant. The 13th position in K-Ras is among the three-hotspot mutations (12 and 61) frequently identified in cancer patients. Our developed hypothesis model (AHHR_4) was screened against the FDA-approved drugs created with the PHASE module’s help. A series of analyses such as XP docking, QikProp, and MM-GBSA analysis was performed for the hit compounds. We found four compounds (azelastine, dihydrocodeine, paroxetine, and tramadol) that act as potent inhibitors. The MD trajectory results found that dihydrocodeine and tramadol exhibited less fluctuation, minimal basin, hydrogen bonds, and potent inhibition. In contrast, the identified compounds could efficiently bound to the H95 groove of K-Ras and hence showing its inhibitory action. Furthermore, the in vivo and in vitro studies on these compounds might be a viable approach to construct a G13D mutant inhibitor with higher drug potency.