Identification of Novel Mutant (R132H) Isocitrate Dehydrogenase 1 Inhibitors for Glioma Therapy

doi:10.21203/rs.3.rs-1344831/v1

Download PDF

Research Article

Identification of Novel Mutant (R132H) Isocitrate Dehydrogenase 1 Inhibitors for Glioma Therapy

https://doi.org/10.21203/rs.3.rs-1344831/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Neomorphic transformation in isocitrate dehydrogenase 1 (IDH1) are the key mutations prevalently found in various cancers including glioma. Recently identified mIDH1 specific inhibitors such as ivosidenib and Vorasidenib were restricted for use due to its modest brain penetrating potential and dose limiting toxicity respectively. Herein, we elucidate integrated virtual screening strategies to discover persuasive mIDH1 inhibitors from the approved subset of the DrugBank database consisting of 2715 molecules. Initially, structural similarity search identified a total of 1432 lead molecules. The resultant compounds were inspected by molecular docking along with MM-GBSA and ADMET analyses. Altogether, the analyses identified DB12001 (Abemaciclib) as the hit against mIDH1. Notably, Abemaciclib was able to form hydrogen bond interaction with active site residues of mIDH1 protein. In the end, the dynamic behavior of the hit complex was also examined using molecular dynamic (MD) simulation studies. The outcome of the study culminates that the hit complex was stable throughout the simulation period of 100ns. It is worth noting that benzimidazole moiety of Abemaciclib was reported to show inhibitory activity against glioma cells. Overall, these findings highlight that DB12001 has the potential as lead molecule against glioma. Indeed the screened hit compound could be further explored for the development of mIDH1 inhibitor with great brain penetrating ability and low toxicity.

Mutated Isocitrate Dehydrogenase 1

Tanimoto co-efficient

Molecular docking

MM-GBSA

Molecular dynamics simulation.

Gliomas are identified as the primary tumor of the central nervous system that derives from glial cells or neural stem cells [1]. It holds the poorest prognosis when compared to other cancer types with increased mortality [2]. Globally, every year 6 per 100,000 individuals were reported with new glioma incidence and the survival rate was found to be at a median of 15 months after conventional cancer treatments [3]. At present the treatment regimes for glioma is restricted to surgical resection along with adjuvant strategies such as chemotherapy and radiotherapy. The inefficiency of the available treatment methods as imposed an emerging necessity to identify new treatment targets and therapeutics against glioma. Recently, mutations in Isocitrate dehydrogenase (IDH) have transformed our understanding of glioma biology and being viewed as an attractive target for the treatment. It catalyzes the conversion of isocitrate to alpha-ketoglutarate (α-KG) through oxidative decarboxylation reaction to generate reduced NADPH from NADP. Further, the generated NADPH is procured for lipid metabolism. In addition, IDH 1 and 2 are key drivers that function at the crossroad of DNA repair, cellular metabolism, redox states and epigenetic regulations. Prevalence of IDH mutations is commonly found in several types of malignancies such as low-grade gliomas, glioblastomas [5], acute myeloid leukemia [6], chondrosarcomas [7] and intrahepatic cholangiocarcinomas [8]. Notably, these somatic mutations appear in the hotspot region of the active site of the enzyme that drives the production of a neomorphic metabolite called D-2-hydroxyglutarate (D-2HG). Shreds of literature evidence highlight that IDH mutation marks the earliest genetic changes that develop during glioma progression [9, 10].

The commonly identified oncogenic mutations of IDH1 protein are mapped to some key catalytic residues within the active site of the enzymes that are crucial for enzyme binding. For instance, literature evidence have identified various mutations of IDH1 which includes R132H, R132C, R132S, R132L and R132G [11].Of note the most prevalent somatic mutation of IDH1 (> 90%) is observed when arginine at 132 positions is substituted with histidine residue (R132H). The significance of IDH1 mutation as galvanized attention as a potential molecular target for glioma treatment. Although mutated IDH1 specific inhibitors were developed, only a few inhibitors were evaluated in clinical trials [12–14]. Of note, ivosidenib is the only drug approved by FDA (in 2018) as a potent inhibitor of mIDH1 for the first-line treatment of Acute Myeloid Leukemia (AML) [15]. However, the use of the ivosidenib was restricted because of its modest blood-brain penetrating capability. Recently, Vorasidenib was identified as dual inhibitor of both the isoforms of mutated IDH proteins with appreciable brain penetrating property [16]. However the pipelines of molecules are needed to overcome the drug resistance situation and to provide better treatment option. Thus, we employed virtual screening and molecular dynamic simulation studies to identify active hits against theR132H mutant IDH1protein.

Structure based virtual screening methods have remarkably featured towards the identification of anti-glioma compounds that target mIDH1. For instance, the molecular docking and MD based study performed by Zou et al in 2016 identified a highly selective mIDH1 inhibitor, FX-03 [17]. In this study about 200,000 molecules were screened by docking against the allosteric site of mIDH1. The cross docking based virtual screening was carried out by Zou et al (2019) to identify novel inhibitor against mIDH1. It is worth mentioning that the study identified ZX06 as the potent inhibitor with modest toxicity [18]. Another study by Zheng et al in 2017 recognized clomifene as an effective natural product based inhibitor against mIDH1 [19]. The screening was performed against ZINC Drug database with 2924 molecules. Duan et al employed structure-activity relationship based study along with MD simulation to discover DC_H31 as novel mIDH1 inhibitor [20]. Recently, Wang et al combined in-silico and in-vitro methods and reported 10 molecules with significant experimental activity against mIDH1 [21]. Despite the number of available evidence, repurposing the approved drug for identifying mIDH1 (R132H) inhibitor is not reported in the literature. Therefore in the present study approved entities of DrugBank database was explored to identify selective and potent mIDH1 inhibitors.

2.1. Structural Refinement

The crude structures of the mIDH1 protein were procured from Protein Data Bank (PDB). The repository contains 19 PDB codes that correspond to the mutated IDH1 protein structure specifically with R132H mutation. The list of PDB IDs along with their resolution was enlisted in Table S1. Here, the protein structure with the ID 6VEI was used to carry out virtual screeningas it possessed appreciable resolution alongside bound Vorasidenib molecule. The raw structures of the mIDH1 protein and small molecules were subjected for preprocessing and minimization using the 'Protein preparation wizard' and 'LigPrep’ module of Schrödinger software respectively [22]. The approved subset of the DrugBank database with approximately 2715 compounds was utilized for virtual screening process. The schematic workflow of our study is highlighted in Fig. 1.

2.2. Tanimoto Coefficient

The molecular similarity search method assesses a compound's likelihood of being active against a therapeutic target according to its structural similarity to known active compounds. [23]. In the present study, the similarity between the compounds is measured in terms of the Tanimoto coefficient (Tc). The Tanimoto coefficient is measured using the following equation.

Tc = N_C / N_A + N_B – N_C (1)

Where N_A denotes the number of fingerprints in structure A,

N_Bdenotes the number of fingerprints in structure B and

N_C denotes the number of fingerprints common in both structures A and B.

The DataStructs subpackage from RDKit was utilized to calculate the Tanimoto coefficient of all the compounds from the approved subsets of the DrugBank database [24].

2.3. Virtual Screening

The compounds with appreciable Tanimoto coefficient were screened and all the compounds were further subjected for molecular docking using the Glide module of the Schrödinger suite.The Glide docking is a grid-based technique that identifies the best binding position of the ligand, based on their affinity towards each other. Prior to molecular docking, the grid was generated around the catalytic site of the mIDH1 protein identified from the literature [25]. A cubic grid of 1nm was constructed identifying the centroid of its active site. Additionally, the partial charge cut-off and van der Waals radius scaling factor were set as 0.25 and 1 respectively to lower the receptor’s non-polar parts potential [26]. After the grid generation, coordinates of the grid were used as an input for the molecular docking process. The advantage of using the Glide module is the three-level hierarchical filters that classify the docking modes into high-throughput virtual screening (HTVS), standard precision (SP) and extra precision (XP).The XP docking mode uses an anchor and grow algorithm that refines the predicted docking pose using an extensive scoring function which identifies and rewards the significant structural entities important for binding [27]. Thus, it enhances the overall accuracy of the docking process.

2.4. Prime MM-GBSA Analysis

The binding affinity of all the screened molecules from XP docking was then rescored using Prime Molecular Mechanics Generalized Born Surface Area (MM-GBSA) analysis. Prime MM-GBSA analysis estimates the relative binding free energy of ligand receptor complex using the following equation.

ΔG_bind = ΔE_MM + ΔG_solv + ΔG_SA (2)

Herein, ΔE_MM indicates the energy difference between protein-ligand complex and the energy summation of unligand and liganded protein. ΔG_solv indicates the GBSA solvation energy difference between protein-ligand complex and the energy summation of unligand and liganded protein. ΔG_SA indicates the energy difference in surface area between protein-ligand complex and the energy summation of uncomplexed and complex protein [28].

2.5. Identification of Drug-like Candidate using In-silico Screening

The pharmacological liabilities of the lead compounds serve as one of the salient features in lead optimization. The growing literature evidence highlights that the increasing failure rate of screened molecules in clinical trials was ascribed majorly due to deviation from Absorption, Distribution, Metabolism, Excretion (ADME) and toxicity properties [29]. Therefore, the QikProp module of Maestro was employed to analyze the drug-likeness of screened molecules. The descriptors such as CNS, stars, blood-brain barrier (logBB) were taken into consideration for the selection of lead compounds. Since mIDH1 protein is predominantly found in glioma, the drug molecule should penetrate the central nervous system. Therefore, the compounds with active CNS scores of 0 and above alongside positive valued log BB were identified as crucial threshold for the hit molecules. The potential acute toxicity endpoints associated with the chemical structure of screened lead molecules were assessed using Protox-II algorithm. The server evaluates various levels of toxicological endpoints such as oral toxicity, toxicity targets and organ toxicity [30]. This in-silico prediction platform is believed to enhance the process of hit selection and optimization and it also provides additional insights into the mechanism of toxicity. Towards the end, the biological activity of the screened molecules was measured in terms of P_a and P_i values obtained from the PASS prediction server [31]. The software predicts the pharmacological effect and biological activity spectrum of the chemical compound based on the structural formula.

2.6. Molecular Dynamic Simulation

The molecular dynamics simulation study of the mIDH1 protein and ligand complex was performed using GROMACS 5.1.2. Firstly, the topology files for the screened small molecules were generated using the PRODRG server in the framework of the GROMOS96 43a1 force field. The docked protein-ligand complex was solvated using a simple point charge (SPC) water solvation model enclosed within the dodecahedron box of 10 Å. The electrostatic energy calculation was computed using the Particle Mesh Ewald method and linear constraint solver algorithm (LINCS) was utilized for covalent bond constraints [32, 33]. Prior to minimization, the neutralization step was carried out by adding 3 sodium ions to the system. Further, the energy minimization was performed by using the steepest descent approach. Berendsen coupling and Parrinello-Rahman method were employed to regulate temperature (NVT) and pressure (NPT) at 300k and 1 bar respectively inside the box. SHAKE algorithm was used to constrain the length of the hydrogen bond [34]. Ultimately, 100 nanosecond simulations were performed for the reference and the hit complexes. Various inbuilt gromacs utilities such as gmx rms, gmx hbond, gmx gyrate, gmx sasa, gmx covar and gmx anaeig were exploited to compute the outcomes of MD simulation and the resultant trajectories were visualized using xmgrace.

3.1. Tanimoto Coefficient

The similarity search algorithm based on molecular fingerprint is a significant tool usually employed in retrospective benchmarking studies. According to the similarity-property principle, increased chemical similarity correlates with a greater possibility that two molecules shares similar activity [35]. Therefore the degree of resemblance between the known inhibitor and the small molecules from the approved subset of the DrugBank database were measured using the Tanimoto coefficient (Tc). The overall structural similarity of all the compounds in our study ranged from 0.16 to 0.52and thus the mean value of the tanimoto co-efficient was used as selection criteria for screening compounds. So, compounds with the threshold value of 0.2 and above were considered for further molecular docking. A total of 1432 compounds possessed a structural similarity of 0.2 and above were used for further analyses.

3.2. Virtual Screening

The approved set of compounds with Tc ≥ 0.2 was retrieved and they were preprocessed using the LigPrep module of the Schrödinger suite. Subsequently, molecular docking was performed for all the compounds to identify the potential lead molecules against mIDH1 protein. Initially, all the compounds (1432 compounds) were subjected to the HTVS mode of docking. As the HTVS method indulges in faster prediction of the appropriate binding mode and ranks the ligand-based on its empirical scoring functions [27]. The top-scoring ligands were selected from HTVS level and used for SP docking. Towards the end, XP docking process was initiated with the high scoring compounds respectively. For instance, 299 molecules were subjected for XP docking and the binding score of each compound was compared with the reference compound Vorasidenib that had a binding score of -5.403 kcal/mol. In comparison, about 59 molecules showed increased binding affinity towards mIDH1 protein. The binding score of the screened compound was ranging from − 5.433 kcal/mol to − 9.964 kcal/mol. The Glide XP GScore of the docked complexes was depicted in Table S2.

3.3. Binding Free energy Analysis

The accuracy of docking was further validated by binding free energy analysis using MM-GBSA algorithm. The reference molecule Vorasidenib exhibited binding free energy (ΔG_bind) of -56.44kcal/mol. Therefore, ΔG_bind of Vorasidenib was used as the threshold to screen molecules with an enhanced binding affinity towards the target protein mIDH1. A total of 16 molecules from the approved set of DrugBank database showed lower binding energies than the reference compound. The results of the compounds were represented in Table 1. It is evident from the table that the binding energies of the compounds vary from − 57.42 kcal/mol to -77.33 kcal/mol. Further the major contributions to the enhanced binding of the hits are due to the exceptionally strong lipophilic and van der Waals interaction. For instance, shreds of literature evidence also highlight the significance of Coulomb, lipophilic energy, electrostatic solvation and van der Waals interactions in the effective binding of protein-ligand complex [36, 37]. Therefore, in addition to ΔG_bind, the above mentioned parameters were also included in selection criteria. Notably, only 10 molecules namely, DB01076, DB13751, DB04868, DB00229, DB01698, DB06769, DB12001, DB00183, DB06590 and DB00492 possessed appropriate energy values. Other screened molecules such as DB01326, DB05039, DB03310, DB13874, DB00293 and DB12789 were eliminated because of its disfavored Coulomb, van der Waals and solvation energies with the values less than − 7.97 kcal/mol, -27.84 kcal/mol, and 11.88 kcal/mol respectively. These selected 10 molecules were then examined for its drug-likeliness screening.

Table 1

The energetic of the reference and the hits from MM-GBSA analysis
S.No.	Compound	MMGBSA dG Bind (kcal/mol)	MMGBSA dG Bind Coulomb (kcal/mol)	MMGBSA dG Bind Lipo (kcal/mol)	MMGBSA dG Bind Solv GB (kcal/mol)	MMGBSA dG Bind vdW (kcal/mol)
1.	Vorasidenib	-56.44	-7.97	-27.84	11.88	-29.58
2.	DB00492	-77.33	-13.60	-41.43	25.26	-48.79
3.	DB06590	-75.12	-34.31	-35.1	42.11	-47.68
4.	DB12789	-71.77	-20.69	-47.58	18.68	-29.50
5.	DB00183	-69.44	-8.64	-38.60	20.3	-46.73
6.	DB12001	-69.26	-9.23	-41.92	20.91	-44.45
7.	DB00293	-68.84	-28.75	-21.27	20.04	-31.48
8.	DB06769	-68.83	-18.14	-40.93	22.08	-30.83
9.	DB01698	-68.69	-19.00	-39.01	17.07	-31.45
10.	DB00229	-67.25	-8.27	-31.51	14.82	-36.34
11.	DB13874	-66.98	-16.66	-28.68	11.46	-35.72
12.	DB03310	-65.49	-34.80	-24.7	44.61	-48.42
13.	DB04868	-64.92	-8.88	-34.56	18.95	-43.26
14.	DB05039	-62.92	-21.96	-25.49	13.55	-24.02
15.	DB01326	-61.75	-19.46	-28.16	12.69	-28.13
16.	DB13751	-61.09	-18.42	-28.98	24.47	-34.85
17.	DB01076	-57.42	-19.21	-34.80	20.95	-30.75

3.4. In-silico drug bioavailability screening

In search of an efficient drug candidate for glioma treatment, the resultant compounds with high binding affinity and binding energies were analyzed for their drug-likeliness and pharmacokinetic properties using the QikProp module of the Schrödinger suite. QikProp module is equipped with a varying range of descriptors that broadly categorizes the absorption, distribution, metabolism and excretion properties of the chemical structure. The inclusion criteria considered for the screening of desired small molecules include the blood-brain barrier (log BB), stars, human oral absorption (HOA) and CNS. The FDA-approved small molecule ivosidenib was restricted from use only due to its poor blood-brain penetrating ability with the logBB value of -1.112 [38]. Thus, the screened molecules with increased logBB value were taken into account. The star descriptor identifies a number of physiochemical features each screened molecule failed to abide within the acceptable limits. The compounds with less than 5 star values were selected as hit compounds. The HOA values of reference and the hits were analyzed and compounds with an HOA value of 3 were scrutinized.The CNS value of the molecules lies between the range of -2 to 2. In glioma the neoplastic cells are predominantly present in the neuronal cells; therefore an effective drug molecule should penetrate into the brain cells to inhibit the activity of tumor cells. Hence, an active CNS score of zero and above was considered for screening. The PK/PD properties of all the screened molecules are shown in Table 2. Amongst the 10 hit molecules, only DB12001 was found to have satisfactory QPlogBB, stars, HOA and CNS properties with the value of 0.14, 1, 3 and 1 respectively. Further, the biological activity of the hit compound was predicted using the PASS server [31]. The server predicts the biological activity by calculating the probability of active and inactive. The potent hit molecule should possess a higher probability of active (P_a) score when compared with the probability of inactivity (P_i). Interestingly, the identified hit was found to have anti-cancerous activity with a P_avalue of 0.5 (Table S3). Finally, the toxicological endpoints of the Vorasidenib and DB12001 were examined using the ProTox-II server [30]. From Table 3 we infer that the reference compound was found to exhibit organ toxicity, especially hepatotoxicity. Interestingly, DB12001 was predicted to have a lethal dose level of 2000 mg/kg and categorized under class IV toxicity. Further, the molecular insights of the reference and the hit were inspected through the interaction pattern and structural scaffolds.

Table 2

The Pharmacokinetics and drug-likeliness properties of the reference and the screened hit molecules
S.No.	Compound
		Lipinski rule of five (Ro5)					Jorgensen’s rule of three (Ro3)				Others
		MWⁱ	HBAⁱⁱ	HBDⁱⁱⁱ	QplogPo/W^iv	NRB^v	QPPCaco^vi	QPlogS^vii	#Star^viii	CNS^ix		HOA^x	QPlogBB^xi
1.	Vorasidenib	414.74	5.00	2.00	4.786	7	2832.41	-6.728	3	1		1	0.328
2.	DB00492	563.67	12.00	1.00	4.13	17	86.59	-5.21	0	-2		2	-1.77
3.	DB06590	420.49	9.70	5.00	-1.427	9	0.906	-3.42	0	-2		2	-2.52
4.	DB00183	767.89	13.25	4.75	3.31	28	0.55	-6.00	12	-2		1	-4.66
5.	DB12001	506.60	9.00	1.00	4.51	7	103.70	-5.81	1	1		3	0.14
6.	DB06769	358.26	4.50	1.00	4.49	10	182.89	-5.36	1	-1		3	-0.75
7.	DB00229	525.61	13.25	3.25	-1.68	12	0.36	-2.99	1	-2		1	-2.91
8.	DB04868	529.52	8.00	2.00	5.90	8	617.10	-9.33	3	-2		1	-1.04
9.	DB13751	822.94	21.30	6.00	1.92	15	0.009	-6.034	12	-2		1	-5.47
10.	DB01076	558.64	6.90	3.00	6.55	16	48.853	-7.521	3	-2		1	-1.99
11.	DB01698	610.52	20.55	9.00	-2.38	16	1.10	-2.70	9	-2		1	-4.89
MWⁱ -Molecular weight
HBAⁱⁱ- Hydrogen Bond acceptor
HBDⁱⁱⁱ - Hydrogen Bond donor
QplogPo/w^iv- Predicted water/octanol partition co-efficient
NRB^v– Number of rotatable bonds
QPPCaco^vi- Predicted apparent Caco-2 cell permeability in nm/sec
QplogS^vii- Predicted aqueous solubility
#stars^viii - Number of property or descriptor values that fall outside the 95% range of similar values for known drugs
CNS^ix - Predicted central nervous system activity on a -2 (inactive) to + 2 (active) scale
HOA^x- Predicted human oral absorption
QPlogBB^xi- Predicted brain/blood partition coefficient

Table 3

Toxicity endpoint analysis of the reference and hit using ProTox-II server
S. No	Compound	Hepatotoxicity	Mutagenicity	Carcinogenicity	Predicted LD₅₀(mg/kg)	Predicted toxicity class
1.	Vorasidenib (reference)	Active	Inactive	Inactive	1700	4
2.	DB12001	Inactive	Inactive	Inactive	2000	4

3.5. Interactions of the Protein-ligand complex

Initially, the binding scheme of all the available PDB codes of mIDH1 protein along with their respective ligand was analyzed using PLIP [39]. Figure 2 highlights the important interacting residues of mIDH1 protein along with their PDB codes. It is evident from the Fig. 2 that most of the mIDH1 proteins were forming hydrophobic interactions with residues such as TRP124, VAL255 which are highlighted in red color. Predominantly, hydrogen bond interactions were formed by the residues GLN277, ASP279, SER280. Similarly, the binding pattern of the Vorasidenib and DB12001 was analyzed and key interactions were identified. The atomic interactions of the hit and the reference molecules were highlighted in different color representations. The interaction of Vorasidenib and DB12001 is highlighted in Fig. 3. The blue color line depicts the existence of hydrogen bond interactions. From the results, we infer that the Vorasidenib was capable of forming two hydrogen bond interactions with the residue GLN277. Furthermore, Vorasidenib was able to form halogen bond interaction between ASP273 residue and chloro-pyridine moiety which is highlighted in cyan color. Comparably, DB12001 was also observed to create two hydrogen bond interactions with the residue GLN277. In addition, DB12001 was forming hydrophobic interactions (dashed line) with some of the conserved residues present in the catalytical site of mIDH1 protein such as TRP124, ASP252, VAL276 and SER280 [16]. Overall, from the above results, we deduce that the additional hydrophobic interaction might contribute to the increased binding affinity of the hit molecule against the mIDH1 protein.

3.6. Scaffold Analysis

Shreds of literature evidence report that amino triazine and chloropyridine moiety of Vorasidenib enhances its binding affinity towards mIDH1 protein [16, 38]. Similarly, the structure of DB12001 was analyzed to identify the important structural moiety that contributes to the increased binding affinity. The benzimidazole group of DB12001 was reported to exhibit anti-cancerous activity especially against glioma [40]. The 2-dimensional structure of both the reference and the hit are depicted in Fig. 4. The generic name of DB12001 is Abemaciclib. It is an anti-cancerous agent that is reported to dually inhibit cyclin-dependent kinase 4 and 6 (CDK 4 and CDK 6). Abemaciclib is an FDA-approved drug used for the treatment of HR-positive and HER-negative metastatic breast cancer in combination with Fulvestrant [41]. Further, the drug has been used in various other trials including for the treatment of melanoma, lymphoma and solid tumors [42, 43]. All these literature evidence highlights that Abemaciclib could also be repurposed for the inhibition of mIDH1 in glioma treatment.Thus, the structural dynamics of Vorasidenib and Abemaciclib were investigated using molecular dynamics simulation studies to enrich the prediction accuracy.

3.7. Molecular Dynamics Simulation

3.7.1. Root Mean Square Deviation (RMSD)

The conformational stability of the protein-ligand complexes were evaluated using MD simulation. The gmx rms utility of gromacs was employed to quantitatively estimate the conformational changes and the stability of the system that occur within the stipulated time boundaries [44]. In our study, the average root mean square deviation (RMSD) of the reference and the hit complexes were calculated for backbone atoms of the protein. Figure 5 highlights the RMSD plots of mIDH1-Vorasidenib (reference) and mIDH1-Abemaciclib (hit) complex. From Fig. 5, it is evident that the reference complex exhibited least deviation from 0 to 27ns. However, the average RMSD of the reference complex increased from ~ 0.41 to ~ 0.65 nm between 30 to 100ns. Although the hit complex showed increased deviation from 0 to 60 ns with the RMSD value of ~ 0.49 nm, the complex attained the state of equilibrium within 70 ns time frame. Interestingly, at the end of the 100ns simulation the complex showed the average RMSD value of ~ 0.48 nm. Thus, from RMSD plot we infer that the mIDH1-Abemaciclib complex proclaimed lesser backbone deviation when compared to mIDH1-Vorasidenib complex.

3.7.2. Root Mean Square Fluctuation (RMSF)

The root mean square fluctuations of individual amino acids were analyzed using gmx rmsf tool of gromacs. In the mIDH1-Vorasidenib complex, PRO147, PRO149, GLY150, ILE154, THR155, ASP160, GLY161, PRO384 and SER415 were few majorly flexible residues. While residues such as GLU17, LYS72, ALA193, SER195, SER196, PHE197, GLN198 and MET199 were relatively stable and expressed subtle fluctuations. In case of mIDH1-Abemaciclib complex LYS3, LYS81, VAL146, PRO147, GLY148, PRO149, THR157, PRO158, SER159, ASP160, GLY161, GLU174, GLY175, GLY176 and GLY177 were found to exhibit high flexibility. On the other hand the residues MET13, GLN14, ASP16, ILE26, LYS27, ILE31, GLU110, ILE113, ILE128, ASP137, SER196, MET199, ILE266, TRP267, THR292, GLU306 and TYR316 were found to exhibit less flexibility as evident from its RMSF value less than ~ 0.03 nm. Of note, active site residues such as ILE113, MET199, TRP267 and THR292 exhibited restricted fluctuation which highlights its involvement in the binding of Abemaciclib (Fig. 6).

3.7.3. Hydrogen bond

The gmx hbond was employed to ascertain the specific inter-molecular interactions between protein-ligand complexes. The stability of the H-bond created between reference and the hit complex was deduced by extracting the time dependent hydrogen bond pattern observed throughout 100ns simulation. The results from the trajectory (Fig. 7) revealed that the reference complex formed an average of ~ 6 hydrogen bonds during the simulation. Whilst the hit complex was capable of forming 3 hydrogen bond interaction with the target protein. From interaction studies, we understand that the hit complex was actively forming hydrophobic interactions with the crucial binding site residues of mIDH1 protein (Fig. 3).

3.7.4. Radius of Gyration

The structural compactness of the reference and hit complex were analyzed using gmx gyrate. The inbuilt gyrate (Rg) tool of gromacs calculates the weighted root mean square distance of collective Cα atoms from the center of mass. Thus, Rg imparts insights on the overall dimensions and the folding state of the target protein [45]. Increased fluctuation in Rg value highlights the unfolding of the target protein. Figure 8 illustrates the radius of gyration plot for the Cα atoms of the mIDH1-Vorasidenib and mIDH1- Abemaciclib complex. Initially, both the complexes were found to have increased RG value of ~ 2.31 nm. However, at the end of 100ns simulation period the Rg value of Vorasidenib and Abemaciclib complex was found to be ~ 2.08 nm and ~ 2.11 nm respectively. It is evident from the results that the reference and hit complexes were stable throughout the simulation time [46].

3.7.5. Solvent Accessible Surface Area (SASA)

SASA estimates the interacting surface area of target protein along with its solvent molecules [47]. The gmx sasa tool was employed to measure the average SASA value of mIDH1-Vorasidenib and mIDH1-Abemaciclib complex throughout the time period of 100ns. The results from SASA plot (Fig. 9) illustrates that the mean SASA value of Vorasidenib and Abemaciclib was 197.39 nm² and 200.71 nm² respectively. The increased SASA value of the hit complex signifies that the internal residues of the mIDH1-Abemaciclib complex are disclosed to the solvent molecules for interactions.The free energy of solvation for hit complex is similar to that of the reference. Thus, the results from SASA analysis highlight the stable binding of the hit complex.

3.7.6. Principal Component Analysis (PCA)

Essential dynamics / PCA aids in identifying the most dominant and probable conformational changes that occurs in the target protein at the time ligand binding. This study allows us to quantify the effect of functionally critical movement upon ligand binding [48]. The important conformational subspace with crucial amount of collective motions of the reference and hit complex are confined within first few eigenvectors. Notably, the first two principal components (PCs) were used for detailed study. Each point on the subspace illustrates a specific conformation of the protein ligand complex. The flexibility of mIDH1-Vorasidenib and mIDH1-Abemaciclib complex was evaluated using the trace value. The trace value for the reference and the hit complex was found to be 26.08 nm² and 14.12 nm² respectively. The higher trace value of the reference compound suggests that there is an increased flexibility and expansion in the collective motion than the hit complex. The 2d projection of the MD trajectories was depicted in Fig. 10 (a) illustrates the concreted motion of system in the phase space spanned by the first two principal components of the reference and hit complex. Interestingly, the results revealed that the Abemaciclib complex occupied smaller phase space along the PC when compared to the reference complex implicating higher flexibility of the reference compound [49]. Notably, the results of 2d projection were found to be in agreement with trace values of the covariance matrix.

Furthermore, covariance calculation also aids to correlate the collective motion of the protein. For instance, when two atoms moves unidirectional they are termed as positively correlated and if the atoms moves in opposite directions they are termed as anti-correlated motion. The positively correlated motion of atoms is highlighted in red regions and the negatively correlated motion in blue region. From Fig. 10, we conclude that both the complexes showed net anti-correlated motion. Overall, less trace value and smaller conformational phase space of the hit molecule suggest its stable binding than the reference compound.

3.7.7. Free Energy Landscape (FEL)

Finally, the conformational states of the PCs were examined in terms of free energy landscape to gain insights on the protein folding. The FEL is merely post-processing of PCs where the difference between free energy is evaluated from probability of energy state occupancy [50]. The conformational states of each complex are represented in different color codes. For instance, the global energy minima are represented in blue and the meta-stable states are depicted in green and cyan. For mIDH1-Vorasidenib, a single narrow energy minima basin was identified, while mIDH1- Abemaciclib complex showed a broader free energy minima basin. Figure 11 highlights the FEL plot of reference and hit complex. From the FEL analysis, we infer that the Abemaciclib complex is thermodynamically stable than mIDH1-reference complex.

The present study focused on exploring plausible mIDH1 protein inhibitors using ligand similarity-based repurposing strategy. The compounds with the similarity coefficient of 0.2 and above were identified from the approved subset of the DrugBank database and used for the molecular docking process. After the hierarchical docking process, about 229 small molecules were screened and their respective binding affinities were compared with the reference compound. The binding energies of the resultant compounds were examined using MM-GBSA analysis. About 16 small molecules with optimum binding affinity were scrutinized and their physicochemical properties such as CNS, QPlogBB, Ro5 and Ro3 were also investigated. Collectively, our analysis resulted DB12001 (Abemaciclib) as a potent lead molecule with binding energy of -69.26 kcal/mol, QlogBB of 0.14 together with active CNS score of 1. The higher Pa value of the Abemaciclib compound highlights its anti-cancerous activity. Interestingly, Abemaciclib exhibits similar binding pattern to that of Vorasidenib with the active site of mIDH1. From the scaffold analysis, we hypothesize that the increased activity of Abemaciclib is mainly due to the existence of benzimidazole moiety in its structure. Our study correlates well with the literature evidence that the benzimidazole scaffold reported to have antineoplastic activity against glioma.In essence, the MD simulation studies also highlights that the Abemaciclib complex maintained a stable conformation throughout the simulation period of 100 ns than the reference molecule.Indeed, we believe that DB12001 (Abemaciclib) portrays desirable attributes to be an effective anti-cancer drug, that could be repurposed for glioma treatment in the impending years.

Ethics approval and consent to participate - Not applicable

Consent for publication - Not applicable

Availability of data and materials- Not applicable

Competing interests- Not applicable

Funding- Not applicable

Authors' contributions - P.M. carried out the experiments and wrote the first draft of the manuscript. R.K.conceived and planned the present study, contributed to the interpretation of the results. and supervised the entire study.

Acknowledgements- The authors thank the management of Vellore Institute of Technology for providing the facilities to carry out this research work.

Disclosure of potential conflicts of interest- We wish to confirm that there are no known conflicts of interest associated with this publication.

Research involving Human Participants and/or Animals -Not applicable

Informed consent - Not applicable

Louis DN, Perry A, Wesseling P, Brat DJ, Cree IA, Figarella-Branger D, Hawkins C, Ng HK, Pfister SM, Reifenberger G, Soffietti, R.(2021) The 2021 WHO classification of tumors of the central nervous system: a summary. Neuro Oncol23:1231-1251.https://doi.org/10.1093/neuonc/noab106
Rouse C, Gittleman H, Ostrom QT, Kruchko C, Barnholtz-Sloan JS (2015)Years of potential life lost for brain and CNS tumors relative to other cancers in adults in the United States, 2010. Neuro-Oncol18:70-77.https://doi.org/10.1093/neuonc/nov249.
Weller M, van den Bent M, Preusser M, Le Rhun E, Tonn JC, Minniti G, Bendszus M, Balana C, Chinot O, Dirven L, French P (2021) EANO guidelines on the diagnosis and treatment of diffuse gliomas of adulthood. Nat Rev ClinOncol 18:170-186.https://doi.org/10.1038/s41571-020-00447-z.
Louis DN, Perry A, Reifenberger G, Von Deimling A, Figarella-Branger D, Cavenee WK, Ohgaki H, Wiestler OD,Kleihues P, Ellison DW (2016) The 2016 World Health Organization classification of tumors of the central nervous system: a summary. ActaNeuropathol131:803-820. https://doi.org/10.1007/s00401-016-1545-1
Molenaar RJ, Verbaan D, Lamba S, Zanon C, Jeuken JW, Boots-Sprenger SH, Wesseling P, Hulsebos TJ, Troost D, VanTilborg AA, Leenstra S (2014) The combination of IDH1 mutations and MGMT methylation status predicts survival in glioblastoma better than either IDH1 or MGMT alone.Neuro-Oncol 16:1263-1273.https://doi.org/10.1093/neuonc/nou005.
Figueroa ME, Abdel-Wahab O, Lu C, Ward PS, Patel J, Shih A, Li Y, Bhagwat N, Vasanthakumar A, Fernandez HF, Tallman MS (2010) Leukemic IDH1 and IDH2 mutations result in a hypermethylation phenotype, disrupt TET2 function, and impair hematopoietic differentiation. Cancer cell18:553-567.https://doi.org/10.1016/j.ccr.2010.11.015.
Pansuriya TC, Van Eijk R, d'Adamo P, Van Ruler MA, Kuijjer ML, Oosting J, Cleton-Jansen AM, Van Oosterwijk JG, Verbeke SL, Meijer D, Van Wezel, T (2011) Somatic mosaic IDH1 and IDH2 mutations are associated with enchondroma and spindle cell hemangioma in Ollier disease and Maffucci syndrome. Nat Genet43:1256-1261.
Borger DR, Tanabe KK, Fan KC, Lopez HU, Fantin VR, Straley KS, Schenkein DP, Hezel AF, Ancukiewicz M, Liebman HM, Kwak EL (2012) Frequent mutation of isocitrate dehydrogenase (IDH) 1 and IDH2 in cholangiocarcinoma identified through broad‐based tumor genotyping. The oncologist17:https://doi.org/10.1634/theoncologist.2011-0386
Dang L, White DW, Gross S, Bennett BD, Bittinger MA, Driggers EM, Fantin VR, Jang HG, Jin S, Keenan MC, Marks KM (2009) Cancer-associated IDH1 mutations produce 2-hydroxyglutarate. Nature 462: 739-744. https://doi.org/10.1038/nature08617
Suzuki H, Aoki K, Chiba K, Sato Y, Shiozawa Y, Shiraishi Y, Shimamura T, Niida A, Motomura K, Ohka F, Yamamoto T (2015) Mutational landscape and clonal architecture in grade II and III gliomas. Nat Genet47: 458-468. https://doi.org/10.1038/ng.3273
Lu VM, McDonald KL (2018) Isocitrate dehydrogenase 1 mutation subtypes at site 132 and their translational potential in glioma. CNS Oncol7: 41-50.https://doi.org/10.2217/cns-2017-0019
Davis MI, Gross S, Shen M, Straley KS, Pragani R, Lea WA, Popovici-Muller J, DeLaBarre B, Artin E, Thorne N, Auld DS (2014) Biochemical, cellular, and biophysical characterization of a potent inhibitor of mutant isocitrate dehydrogenase IDH1. J Biol Chem 289: 13717-13725.https://doi.org/10.1074/jbc.M113.511030
Deng G, Shen J, Yin M, McManus J, Mathieu M, Gee P, He T, Shi C, BedelO, McLean LR,Le-Strat F (2015) Selective inhibition of mutant isocitrate dehydrogenase 1 (IDH1) via disruption of a metal binding network by an allosteric small molecule. JBiol Chem 290: 762-774.https://doi.org/10.1074/jbc.M114.608497
Foran JM, DiNardo CD, Watts JM, Stein EM, De Botton S, Fathi AT, Prince GT, Stein A, Stone RM, Patel P, Tallman MS(2019)Ivosidenib (AG-120) in patients with IDH1-mutant relapsed/refractory myelodysplastic syndrome: Updated enrollment of a phase 1 dose escalation and expansion study. Blood 134:4254.https://doi.org/10.1182/blood-2019-123946
Dhillon S (2018)Ivosidenib: first global approval. Drugs78:1509-1516. https://doi.org/10.1007/s40265-018-0978-3
Konteatis Z, Artin E, Nicolay B, Straley K, Padyana AK, Jin L, Chen Y, Narayaraswamy R, TongS, Wang F, Zhou D (2020) Vorasidenib (AG-881): a first-in-class, brain-penetrant dual inhibitor of mutant IDH1 and 2 for treatment of glioma. ACS MedChem Lett11: 101-107.https://doi.org/10.1021/acsmedchemlett.9b00509
Zou F, Pusch S, Eisel J, Ma T, Zhu Q, Deng D, Gu Y, Xu,Y, von Deimling A, Zha X (2016) Identification of a novel selective inhibitor of mutant isocitrate dehydrogenase 1 at allosteric site by docking-based virtual screening. RSC advances. RSC Adv 6: 96735–96742. https://doi.org/10.1039/C6RA21617J
Zou F, Pusch S, Hua J, Ma T, Yang L, Zhu Q, Xu Y, Gu Y, von Deimling A, Zha X(2018) Identification of novel allosteric inhibitors of mutant isocitrate dehydrogenase 1 by cross docking-based virtual screening. Bioorganic Med Chem Lett 28: 388–393. https://doi.org/10.1016/j.bmcl.2017.12.030
Zheng M, Sun W, Gao S, Luan S, Li D, Chen R, Zhang Q, Chen L, Huang J, Li H (2017) Structure based discovery of clomifene as a potent inhibitor of cancer-associated mutant IDH1. Oncotarget 8: 44255–44265. doi: 18632/oncotarget.17464
Duan Z, Liu J, Niu L, Wang J, Feng M, Chen H, Luo C(2019) Discovery of DC_H31 as potential mutant IDH1 inhibitor through NADPH-based high throughput screening. Bioorganic Med Chem 27: 3229–3236.https://doi.org/10.1016/j.bmc.2019.05.040
Wang Y, Tang S, Lai H, Jin R, Long X, Li N, Tang Y, Guo H, Yao X, Leung ELH(2020) Discovery of novel IDH1 inhibitor through comparative structure-based virtual screening. Front Pharmacol11, 1–11. https://doi.org/10.3389/fphar.2020.579768
Sastry GM, Adzhigirey M, Day T, Annabhimoju R, Sherman W (2013) Protein and ligand preparation: parameters, protocols, and influence on virtual screening enrichments. J Comput Aided Mol Des 27: 221-234. https://doi.org/10.1007/s10822-013-9644-8.
Liu R, AbdulHameed MDM, Wallqvist A (2019) Teaching an Old Dog New Tricks: Strategies That Improve Early Recognition in Similarity-Based Virtual Screening. Front Chem7:701.https://doi.org/10.3389/fchem.2019.00701
Vogt M, Bajorath J (2020) ccbmlib–a Python package for modeling Tanimoto similarity value distributions. F1000Research9.doi: 12688/f1000research.22292.2
Poonan P, Agoni C, Soliman M (2021) Dual‐Knockout of Mutant Isocitrate Dehydrogenase 1 and 2 Subtypes Towards Glioma Therapy: Structural Mechanistic Insights on the Role of Vorasidenib. ChemBiodivers 18: e2100110.https://doi.org/10.1002/cbdv.202100110
Anusuya S, Velmurugan D, Gromiha MM (2016) Identification of dengue viral RNA-dependent RNA polymerase inhibitor using computational fragment-based approaches and molecular dynamics study. J BiomolStructDyn 34:1512-1532.https://doi.org/10.1080/07391102.2015.1081620
Repasky MP, Shelley M, Friesner RA(2007)Flexible ligand docking with Glide. CurrProtocBioinform 18:8-12.https://doi.org/10.1002/0471250953.bi0812s18
Lyne PD, Lamb ML,Saeh JC (2006) Accurate prediction of the relative potencies of members of a series of kinase inhibitors using molecular docking and MM-GBSA scoring. J Med Chem 49:4805-4808.https://doi.org/10.1021/jm060522a
Egan WJ, Merz KM, BaldwinJJ (2000) Prediction of drug absorption using multivariate statistics. Journal of medicinal chemistry. J MedChem43:3867-3877.https://doi.org/10.1021/jm000292e
Banerjee P, Eckert AO,Schrey AK, Preissner R (2018) ProTox-II: a webserver for the prediction of toxicity of chemicals. Nucleic acids research. Nucleic Acids Res46: W257-W263.https://doi.org/10.1093/nar/gky318
Lagunin A, Stepanchikova A, Filimonov D Poroikov V (2000) PASS: prediction of activity spectra for biologically active substances. Bioinformatics16: 747-748.https://doi.org/10.1093/bioinformatics/16.8.747
Wang H, Dommert F, Holm C (2010) Optimizing working parameters of the smooth particle mesh Ewald algorithm in terms of accuracy and efficiency. JChem Phys133:034117.https://doi.org/10.1063/1.3446812
Amiri S, Sansom MS, Biggin PC (2007)Molecular dynamics studies of AChBP with nicotine and carbamylcholine: the role of water in the binding pocket. Protein Eng Des Sel 20: 353-359.https://doi.org/10.1093/protein/gzm029
Kumar A, Rajendran V, Sethumadhavan R, Purohit R (2013)Molecular dynamic simulation reveals damaging impact of RAC1 F28L mutation in the switch I region. PloS one8:e77453.https://doi.org/10.1371/journal.pone.0077453
Wang T, Wu MB, Lin JP, Yang LR (2015)Quantitative structure–activity relationship: promising advances in drug discovery platforms. Expert Opin Drug Discov10: 1283-1300.https://doi.org/10.1517/17460441.2015.1083006
Li J, Zhou N, Luo K, Zhang W, Li X, Wu C, Bao J (2014) In silico discovery of potential VEGFR-2 inhibitors from natural derivatives for anti-angiogenesis therapy. Int J Mol Sci. 15:15994-16011. https://doi.org/10.3390/ijms150915994
Kellici TF, Ntountaniotis D, Liapakis G, Tzakos AG, Mavromoustakos T (2019)The dynamic properties of angiotensin II type 1 receptor inverse agonists in solution and in the receptor site. Arab J Chem12:5062-5078.https://doi.org/10.1016/j.arabjc.2016.11.014
Golub D, Iyengar N, Dogra S, Wong T, Bready D, Tang K, Modrek AS, Placantonakis DG (2019) Mutant isocitrate dehydrogenase inhibitors as targeted cancer therapeutics. Front Oncol9: 417.https://doi.org/10.3389/fonc.2019.00417
Adasme MF, Linnemann KL, Bolz SN, Kaiser F, Salentin S, Haupt VJ, Schroeder M (2021) PLIP 2021: Expanding the scope of the protein–ligand interaction profiler to DNA and RNA. Nucleic Acids Res 1:1-5.https://doi.org/10.1093/nar/gkab294
YurttasL, Demirayak S, AkalınCiftci G (2015) Cytotoxic, antiproliferative and apoptotic effects of new benzimidazole derivatives on A549 lung carcinoma and C6 glioma cell lines. Anti-Cancer Agents Med Chem15: 1174-1184.
Wender IO, Haines K, Jahanzeb M (2020) Response to abemaciclib after 10 lines of therapy including palbociclib in metastatic breast cancer: a case report with literature review. OncolTher8: 351-358. https://doi.org/10.1007/s40487-020-00126-0
Patnaik A, RosenLS, Tolaney SM, Tolcher AW, Goldman JW, Gandhi L, Papadopoulos KP, Beeram M, Rasco DW, Hilton JF, Nasir A (2016) Efficacy and safety of abemaciclib, an inhibitor of CDK4 and CDK6, for patients with breast cancer, non–small cell lung cancer, and other solid tumors. Cancer Discov6: 740-753.DOI:1158/2159-8290.CD-16-0095
Tate SC, Sykes AK, Kulanthaivel P, Chan EM, Turner PK, Cronier DM (2018) A population pharmacokinetic and pharmacodynamic analysis of abemaciclib in a phase I clinical trial in cancer patients. ClinPharmacokinet57:335-344.https://doi.org/10.1007/s40262-017-0559-8
Schreiner W, Karch R, Knapp B, Ilieva N (2012) Relaxation estimation of RMSD in molecular dynamics immunosimulations. Comput Math Methods Medhttps://doi.org/10.1155/2012/173521
Baig MH, Sudhakar DR, Kalaiarasan P, Subbarao N, Wadhawa G, LohaniM, Khan MKA, Khan AU (2014)Insight into the effect of inhibitor resistant S130G mutant on physico-chemical properties of SHV type beta-lactamase: A molecular dynamics study. PLoS One 9: e112456.https://doi.org/10.1371/journal.pone.0112456
Baig MH, Danishuddin M, Khan S Khan AU (2012)Screening of inhibitors for S130G inhibitor resistant mutants of TEM type beta-lactamase. Bioinformation8: 1225.doi: 6026/97320630081225
Ausaf Ali S, Hassan I, Islam A, Ahmad F (2014) A review of methods available to estimate solvent-accessible surface areas of soluble proteins in the folded and unfolded states. Curr Protein PeptSci15:456-476.
Wolf A, Kirschner KN (2013) Principal component and clustering analysis on molecular dynamics data of the ribosomal L11· 23S subdomain. JMol Model19:539-549.https://doi.org/10.1007/s00894-012-1563-4.
Pandey B, Grover S, Tyagi C, Goyal S, Jamal S, Singh A, Kaur J,Grover A (2016)Molecular principles behind pyrazinamide resistance due to mutations in panD gene in Mycobacterium tuberculosis. Gene581: 31-42.https://doi.org/10.1016/j.gene.2016.01.024
Tavernelli I,Cotesta S, Di Iorio EE (2003)Protein dynamics, thermal stability, and free-energy landscapes: a molecular dynamics investigation. Biophys J85: 2641-2649.https://doi.org/10.1016/S0006-3495(03)74687-6.

INDSupplementaryfile.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Identification of Novel Mutant (R132H) Isocitrate Dehydrogenase 1 Inhibitors for Glioma Therapy

Status:

Version 1

Abstract

Figures

1. Introduction

2. Methodology

2.1. Structural Refinement

2.2. Tanimoto Coefficient

2.3. Virtual Screening

2.4. Prime MM-GBSA Analysis

2.5. Identification of Drug-like Candidate using In-silico Screening

2.6. Molecular Dynamic Simulation

3. Results And Discussion

3.1. Tanimoto Coefficient

3.2. Virtual Screening

3.3. Binding Free energy Analysis

3.4. In-silico drug bioavailability screening

3.5. Interactions of the Protein-ligand complex

3.6. Scaffold Analysis

3.7. Molecular Dynamics Simulation

3.7.1. Root Mean Square Deviation (RMSD)

3.7.2. Root Mean Square Fluctuation (RMSF)

3.7.3. Hydrogen bond

3.7.4. Radius of Gyration

3.7.5. Solvent Accessible Surface Area (SASA)

3.7.6. Principal Component Analysis (PCA)

3.7.7. Free Energy Landscape (FEL)

4. Conclusion

Declarations

References

Supplementary Files

Status:

Version 1