Identification of novel dihydroorotate dehydrogenase (DHODH) inhibitors for cancer: computational drug repurposing strategy

TL;DR

This study uses computational methods to find new drugs that can inhibit a cancer-related enzyme called DHODH, with the goal of improving cancer treatment.

Contribution

The study introduces a computational drug repurposing strategy to identify novel DHODH inhibitors with improved potency and reduced toxicity.

Findings

Two FDA-approved molecules, DB09026 and DB00503, were identified as potent DHODH inhibitors.

Molecular docking and simulation studies confirmed strong binding to key DHODH residues and anti-cancer activity.

Ritonavir and Aliskiren are proposed as promising candidates with minimal side effects for cancer treatment.

Abstract

Dihydroorotate dehydrogenase (DHODH) is a crucial enzyme in de novo pyrimidine production, initially sought since its disruption is frequently observed in malignancies. DHODH inhibitors have been demonstrated in multiple trials to effectively destroy tumour cells. For instance, leflunomide, teriflunomide and brequinar are currently in practice for DHODH based therapeutics. However, their usage is hampered due to their less efficiency and toxicity issues. Adding together, no studies have reported drug repurposing efforts targeting DHODH.

To address these challenges, the present study aimed to identify novel and potent DHODH inhibitors through virtual screening, with a distinct focus on repurposing. Initially, 2619 FDA approved molecules were subjected to molecular docking using AutoDock Vina and Molsoft ICM-Pro. Consequently, binding free energy were performed using Uni-GBSA and PRODIGY. Toxicity and cancer cell line activity were assessed using high precision machine learning techniques. In the end, gold standard simulation studies executed to validate the hit compound inhibitory activity against DHODH protein.

The results of our analysis identified two molecules, DB09026 and DB00503, as potent DHODH inhibitors. It is worth noting that the identified compound able to bind with key residues in the DHODH target protein. Moreover, scaffold analysis supports the existence of anti-cancer activity of the identified compounds. In essence, long 100ns molecular dynamic simulation results were also correlates well with the previous results.

Collectively, we hypothesize that both ritonavir and Aliskiren exhibits minimal side effect, it could be of interesting choice for the management of cancer due to its improved potency.

Not applicable.

Figures

Superimposition of docked brequinar and cocrystallized complex in the active site of 1D3G protein a) Superimposition of docked brequinar (red) and cocrystallized complex (green) in the active site of 1D3G protein using AutoDock Vina b) Superimposition of docked brequinar (white) and cocrystallized complex (yellow) in the active site of 1D3G protein using Molsoft ICM pro

Interaction analysis of a) Reference - Brequinar; b) DB09026 Aliskiren and c)DB00503 Ritonavir towards the target protein, blue-coloured lines indicate hydrogen bond interactions and grey colour line signifies hydrophobic interactions

Scaffold analysis of a) Reference-brequinar with Quinoline carboxylic acid moiety; b) DB09026 Aliskiren with 4-methoxy-3-(3-methoxypropoxy) phenyl moiety and c) DB00503 Ritonavir with Thiazole moiety

RMSD profiles of the protein–ligand complexes are shown, with brequinar as the reference (black), Aliskiren (DB09026) in orange, and Ritonavir (DB00503) in violet

Root mean square fluctuation (RMSF) profiles of the protein–ligand complexes. Brequinar (reference) is shown in black, Aliskiren (DB09026) in orange, and Ritonavir (DB00503) in violet

Hydrogen bond (H-bond) profiles of the protein–ligand complexes, with Brequinar (reference) represented in black, Aliskiren (DB09026) in orange, and Ritonavir (DB00503) in violet

Radius of Gyration (Rg) profiles of the protein–ligand complexes. Brequinar (reference) is shown in black, Aliskiren (DB09026) in orange, and Ritonavir (DB00503) in violet

Solvent accessible surface area (SASA) profiles of the protein–ligand complexes. Brequinar (reference) is depicted in black, Aliskiren (DB09026) in orange, and Ritonavir (DB00503) in violet

Full text

Introduction

Pyrimidine nucleotide play a crucial role in cell proliferation and synthesis of DNA and RNA precursors [1]. The production of pyrimidine nucleotides involve two significant pathways namely salvage synthesis and de novo pyrimidine biosynthesis [2]. De novo pyrimidine biosynthesis produces more number of pyrimidines for the cell growth and it is essential in production of RNA and DNA synthesis [3]. One of the key players in de novo pyrimidine biosynthesis is dihydroorotate dehydrogenase (DHODH), a redox enzyme that catalyzes the conversion of dihydoorotate to orotate (rate limiting step) [4, 5]. Recently, literature evidence has shown that DHODH is involved in numerous processes, such as cellular metabolism, growth signalling, ferroptosis, transcription, carcinogenesis, and tumour metastasis, in addition to its role in pyrimidine nucleotide production [2]. Notably, elevated DHODH expression has been implicated in the advancement of numerous diseases, including cancer. Consequently, it represents a crucial therapeutic target [6, 7].

In recent years, numerous small-molecule DHODH inhibitors have emerged as promising therapeutic agents, particularly for cancer treatment. Among them, FDA-approved inhibitors such as leflunomide (LEF) and teriflunomide (TFM), originally used to treat autoimmune diseases like rheumatoid arthritis (RA) and multiple sclerosis (MS), respectively, received a black-box warning in 2010 due to their limited anticancer efficacy and risk of acute liver failure [8]. However, brequinar a potent DHODH inhibitor primarily used for cancer treatment [9], demonstrated higher potency than leflunomide (LEF) and teriflunomide (TFM). Although it showed greater efficacy, its clinical development faced setbacks due to concerns over toxicity and limited efficacy in solid tumors, ultimately leading to the failure of clinical trials [10]. Additionally, several small-molecule DHODH inhibitors, including BAY2402234, ASLAN003, PTC299, JNJ74856665, AG-636, and RP7214, have entered clinical trials for cancer treatment [11].

In one of the study, the integration of in vitro and in vivo studies to induce reactive oxygen species production and inhibit DHODH activity. The results showed that naphthol [2,3-d] [1,2,3]-triazole-4-9-dione compounds exhibited antitumor effects [12]. Similarly, the synthesis and evaluation of indoluidin and its derivatives against various cancer cell lines. The study showed that the compounds had a significant antitumor efficacy in A549 xenograft model [13]. Another investigation assessed the teriflunomide derivatives against colorectal cancer using a biphenyl scaffold. The findings revealed that A37 would be an efficient DHODH inhibitor for the management of colorectal cancer [14]. Furthermore, a study reported the link between DHODH enzyme and lysine degradation using a bioinformatics approach and highlighted its importance as a therapeutic target for cancer therapy. Despite the number of recent studies, none of these small molecules have succeeded in clinical trials. These challenges underscore the urgent need for the development of more potent and selective DHODH inhibitors with improved efficacy and safety profiles.

Previous attempts to develop DHODH inhibitors were unsuccessful due to toxicity and low efficacy [8]. Adding together, no studies have reported drug repurposing efforts targeting DHODH. To address these challenges, the present study aimed to identify novel and potent DHODH inhibitors through virtual screening, with a distinct focus on repurposing. In-silico approaches like virtual screening and molecular docking are very significant in the initial stages of discovering new drugs as they facilitate rapid and cost-effective means for the identification of lead compounds [4]. It helps to screen large repositories with high precision. Such an initial screening helps to retrieve potent compounds with less computational cost and ultimately drive the cancer drug therapeutics pipeline [5]. Indeed, this present strategy of combining machine learning with high-end simulation approach offers a promising avenue for overcoming past limitations and advancing DHODH targeted cancer therapy.

Methods

Data set curation and docking validation

The 3D structural information of the DHODH protein (PDB ID: 1D3G) in complex with brequinar at a resolution of 1.60 Å was retrieved from the Protein Data Bank and used as the reference compound. A total of 2619 FDA approved compounds were procured from DrugBank repository and utilized for the virtual screening. Initially, the accuracy and reliability of the docking was assessed by employing superimposition technique. For this purpose, complexes were generated both using AutoDock Vina and Molsoft ICM-Pro and the results are analysed based on the superimposed (RMSD) Root Mean Square Deviation.

Molecular docking

Molecular docking using AutoDock Vina

The binding affinity of the protein-ligand complexes was analyzed through molecular docking using AutoDock Vina v1.1.2 and Molsoft ICM-Pro. Prior to docking, hydrogen atoms were incorporated into the protein using AutoDock 4.2.6. All the water molecules and bounded ligands were systematically eliminated from the PDB structure. Subsequently, protein was assigned with Kollman charges and polar hydrogen bonds [15]. The gasteiger charges were assigned, the torsional root was augmented, and subsequently, brequinar was introduced to the target protein [16]. The grid box was generated with center coordinates (x = 49.872, y = 40.422, z = -4.83) around the active site region. Furtherly the compounds were subjected docking analysis using AutoDock Vina v1.1.2. It utilizes Monte Carlo algorithm combined with BFGS (Broyden-Fletcher-Goldfarb-Shanno) gradient-based optimizer [17] for calculating the docking score.

Molecular docking using Molsoft ICM pro

The docking analysis was also executed using Molsoft ICM-Pro software to avoid the false positive prediction. The protein structure (1D3G) was turned into ICM objects by removing water molecules and optimising hydrogen bonds. Missing side chains were treated before the docking. The Vander Waals contact of the receptor is utilised to locate pockets for a grid construction [18]. Biased probability Monte Carlo was used in docking to optimize the small molecules internal coordinates [19].

Binding free energy calculations

The Molecular mechanics/Generalized-Born (Poisson–Boltzmann) surface area (MM/GB(PB)SA), is one of the most popular method, achieved reasonable correlation with experimental affinities was employed to validate the accuracy of the docking results [20, 21]. In the present study, Uni-GBSA and PRODIGY tools were employed to ascertain the ΔG_bind_ of the docked protein - ligand complexes.

Uni-GBSA is an automatic workflow platform used to calculate MM/GBSA from force field building and structure optimize free energy calculation [22]. The binding free energy is determined using the following equation:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta\text{G}=\Delta\text{G}_{\text{protein-ligand}}-\Delta\text{G}_{\text{protein}}-\Delta\text{G}_{\text{ligand}}$$\end{document}

where the energy estimates of the optimized complex (protein-ligand), optimized free ligand, and optimized free protein [23]. The ΔG is estimated as ΔG = ΔH - TΔS, ΔH is expressed as ΔH = ΔE_MM_ + ΔG_solv_. Δ_EMM_ is the differences between the protein ligand complex as the total minimized energies and ΔG_solv_ is the sum of polar and non-polar concentrations of the protein ligand complex. -TΔS represents the protein-ligand complex’s conformational entropy, where T is the absolute temperature are represented respectively [22].

Similarly, PRODIGY-LIG, an online web server, uses a combination of structural properties like inter-molecular electrostatic energy and number of intermolecular atomic contacts from experimental datasets to predict the ΔG bind of protein-ligand complexes [24]. The ∆G_prediction_ is calculated as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}\Delta{\text{G}}_{\text{n}\text{o}\text{e}\text{l}\text{e}\text{c}\text{t}}\:=&\:0.0354707\:\text{*}\:\text{A}\text{C}\text{N}\text{N}\hspace{0.17em}-\hspace{0.17em}0.1277895\:\text{*}\\&\text{A}\text{C}\text{C}\text{C}\hspace{0.17em}-\hspace{0.17em}0.0072166\:\text{*}\:\text{A}\text{C}\text{C}\text{N}\hspace{0.17em}-\hspace{0.17em}5.1923181\end{aligned}$$\end{document}

where the intermolecular atomic connections (ACs) (within a cutoff of 10.5 Å) is categorised based on the atoms that interact O = Oxygen, C = Carbon, X = other atoms, and N = nitrogen [25, 26].

In-silico toxicity prediction and interaction analysis

The Protox-II algorithm (https://tox.charite.de/protox3/) was used to assess the potential acute toxic endpoints of the screened lead compounds. The server predicts LD_50_ and evaluates a variety of toxicological endpoints, including organ toxicity, oral toxicity, and toxicity targets [27]. ProTox-II uses molecular similarity, pharmacophore-based, fragment propensities, and machine learning models to predict multiple toxicity endpoints. This in-silico prediction platform is expected to improve the hit selection and optimisation process while also providing fresh insights into the toxicity mechanism [28]. Finally, the binding of ligand in DHODH was evaluated using the protein-ligand interaction profiler server [20].

In-silico biological activity prediction

The biological anticancer activity of the reference and screened lead compounds were assessed using the pdCSM-cancer tool (https://biosig.lab.uq.edu.au/pdcsm_cancer/prediction) Notably, pdCSM-cancer employs a graph-based signature approach, which has been demonstrated to efficiently represent chemical and biomolecular datasets, as well as to predict pharmacokinetics, toxicity, and bioactivity. Specifically, pdCSM-cancer forecasts the activity of compounds against cancer cell lines based on their graph-derived structural features. The model was developed and validated using clinical data comprising over 18,000 compounds tested across 74 cancer cell lines and 9 tumor types [29]. The literature evidences highlights that Pearson correlation coefficients (r) ranging from 0.70 to 0.85 and root mean square error (RMSE) values between 0.5 and 0.9 for IC_50_ predictions across various cancer cell lines, based on GDSC and CCLE datasets [30]. These metrics underscore the high fidelity of pdCSM-cancer’s predictions when benchmarked against experimental drug sensitivity data. The GI50% was predicted using the SMILES notation of the compounds as input. It predicts the cell-line growth inhibitory action of substances in logarithmic units (log µg/ml) [31].

Molecular dynamic simulations

The dynamic and adaptable properties of the complexes were investigated at the molecular level using molecular dynamic simulations. In this investigation, MD simulations for the DHODH complexes were carried out using GROMACS version 2020.2. To create a ligand topology, docked ligands were obtained and used with the CHARMM General Force Field (CGenFF). Furthermore, the protein topology was constructed using the CHARMM36 force field. The systems were housed inside a dodecahedron box and solvated using the Simple Point Charge water model. Water molecules were replaced with counter-ions in the solution to provide complete neutrality [20, 27]. Initially, the energy was minimised over 50,000 cycles using the steepest descent approach to resolve steric disputes. Subsequently, the system was further minimised and equilibrated for 1000 ps under the NVT and NPT stages. Using a modified Berendsen thermostat for temperature coupling, a modified Parrinello-Rahman pressure coupling approach, and the V-rescale method, the system was kept at 300 K and 1 bar. After the systems had stabilised, they were subjected to a 100-ns production run at 300 K and 1 bar of pressure. The resultant trajectories were examined using GROMACS tools for several metrics, including radius of gyration (Rg), solvent accessible surface area (SASA), root mean square fluctuation (RMSF), and root mean square deviation (RMSD). The Grace tool was used to graphically visualise these parameters [32].

Results and discussion

Docking validation

The docking approach was validated by re-docking the brequinar into the active region of the DHODH protein. The superimposed fit was assessed using the RMSD value between the crystallographic structure and the redocked complexes. The RMSD value of less than 2.0 Å is considered high correlation with experimentally derived structure. Here, the RMSD value was observed to be 0.60 Å and 0.25 Å for AutoDock Vina and Molsoft ICM-Pro respectively. The super imposition of the redocked complex (red) and the co-crystalised complex (green) using AutoDock Vina is shown in Fig. 1a. Redocked complex (white) and co-crystalised complex (yellow) using Molsoft ICM pro is represented in Fig. 1b.

Fig. 1. Superimposition of docked brequinar and cocrystallized complex in the active site of 1D3G protein a) Superimposition of docked brequinar (red) and cocrystallized complex (green) in the active site of 1D3G protein using AutoDock Vina b) Superimposition of docked brequinar (white) and cocrystallized complex (yellow) in the active site of 1D3G protein using Molsoft ICM pro

Molecular docking analysis

The crucial residues in the binding pockets of DHODH were explored from the literature evidences before docking analysis [33]. The active site comprised of hydrophobic residues such as TYR38, MET43, LEU46, LEU50, ALA55, LEU58, ALA59, PHE62, LEU67, LEU68, PRO69, PHE98, MET111, LEU359, and PRO364 alongside ARG136, GLN47, TYR356 and THR360, crucial player in hydrogen bonding interaction [33, 34]. Subsequently, brequinar and 2619 FDA- approved molecules were docked against the binding pocket of DHODH using AutoDock tool and Molsoft ICM-Pro. According to the docking results, the binding affinity of the brequinar was found to be -13.20 kcal/mol in AutoDock Vina and − 16.47 kcal/mol in Molsoft ICM-Pro. The binding affinities of the drug molecules library ranged from − 98.8 to -1.0 kcal/mol (Table 1). The docking results shown that total 920 compounds exhibited better results in the AutoDock Vina and 1027 compounds shown a better results in the Molsoft ICM-Pro. The results of the docking algorithms are integrated to avoid the false positive prediction. This integration yielded a total of 205 compounds that exhibited better binding affinity than brequinar.

Table 1. Docking results of the AutoDock Vina and Molsoft ICM pro of reference brequinar and 17 lead moleculesS.noDrugBank IDNameAutoDockVina (kcal/mol)Molsoft ICM pro (kcal/mol)1DB03480Brequinar-13.2-16.472DB00137Lutein-17.6-20.833DB00481Raloxifene-16.8-17.154DB00503Ritonavir-16-31.165DB00549Zafirlukast-15-41.486DB00654Latanoprost-17-58.477DB00944Demecarium-18.3-31.898DB00947Fulvestrant-17.3-34.649DB01100Pimozide-15.9-22.9410DB06249Arzoxifene-14.3-17.5711DB06401Bazedoxifene-16.7-16.8412DB06695Dabigatran etexilate-14.3-42.8913DB08909Glycerol phenylbutyrate-16.2-56.5014DB08964Gemeprost-15.8-55.3115DB09026Aliskiren-16.1-16.6816DB09030Vorapaxar-16.6-18.1417DB13615Mifamurtide-14.2-48.3218DB14185Aripiprazole lauroxil-15.8-45.67

Binding free energy analysis using Uni-GBSA and PRODIGY

The molecular docking was further rescored by binding free energy analysis using the MM/GB(PB)SA algorithm [35]. The ΔG_bind_ of brequinar was determined to be -67.57 kcal/mol used as a threshold to scrutinise the 205 obtained compounds. It is evident that the 51 compounds exhibits better binding energy profile against DHODH protein.

Further, the binding free energies of 51 compounds were re-evaluated using PRODIGY-LIG tool [26]. The binding free energy of the molecules spanned from − 14.65 to -12.00 kcal/mol. The binding free energy of the brequinar (-10.6 kcal/mol) was used as a threshold to sort the lead molecules. It is important to note that a total 17 molecules showed better binding free energy than the reference. Therefore, these 17 lead compounds were considered for further analysis. The binding free energy calculations are tabulated in Table 2.

Table 2. Binding free energy calculations of reference brequinar and 17 lead moleculesS.noDrugBank IDName(MMGBSA) Molecular mechanics with generalized born and surface area solvationProdigy webserverVan der Waals (Kcal/mol)Electrostatic (Kcal/mol)Solvation (Kcal/mol)ΔG_bind_ (Kcal/mol)ΔG prediction (Kcal/mol)1DB03480Brequinar-63.78-9.655.86-67.57-10.62DB00137Lutein-81.44-1.90-2.77-86.11-15.213DB00481Raloxifene-70.07-3.010.96-72.11-13.254DB00503Ritonavir-91.30-10.035.66-95.67-14.445DB00549Zafirlukast-78.64-0.76-0.75-80.16-12.126DB00654Latanoprost-69.61-6.373.01-72.98-13.037DB00944Demecarium-87.52-2.30-4.58-94.41-12.538DB00947Fulvestrant-78.163.85-4.94-79.24-15.919DB01100Pimozide-66.010.36-2.74-68.40-12.2510DB06249Arzoxifene-69.64-2.520.28-71.89-13.2911DB06401Bazedoxifene-71.53-1.440.90-72.06-13.2812DB06695Dabigatran etexilate-100.12-5.125.00-100.24-12.0013DB08909Glycerol phenylbutyrate-84.14-4.652.07-86.72-14.6514DB08964Gemeprost-63.84-7.973.75-68.05-12.0615DB09026Aliskiren-67.19-11.616.77-72.03-12.5016DB09030Vorapaxar-69.86-3.572.12-71.32-12.7817DB13615Mifamurtide-101.30-18.9820.64-99.64-12.2518DB14185Aripiprazole lauroxil-91.63-2.20-1.73-87.68-13.93

In- silico drug toxicity prediction

The fundamental challenge in lead optimization is distinguishing drug-like molecules from non-drug compounds [27]. The significant toxicological endpoints, such as carcinogenicity, cytotoxicity and mutagenicity, were investigated for both the reference and lead molecules [28]. Table 3 shows the toxicity profile of the brequinar and the lead molecules. The brequinar was seen to be inactive with LD_50_ of 495 mg/kg. The results of 17 lead compounds fall under the toxicity class of II – V with LD_50_ range from 6 mg/kg to 4000 mg/kg. The analysis revealed that 15 lead molecules except DB00944 and DB06249 showed active results in the carcinogenicity and cytotoxicity.

Table 3. Toxicity prediction using PROTOX-II webserver for reference brequinar and17 lead moleculesS.noDrugBank IDNameLD50Toxicity classCarcinogenicityMutagenicityCytotoxicity1DB03480Brequinar495 mg/kgIVInactiveInactiveInactive2DB00137Lutein10 mg/kgIIInactiveInactiveInactive3DB00481Raloxifene400 mg/kgIVInactiveInactiveInactive4DB00503Ritonavir1000 mg/kgIVInactiveInactiveInactive5DB00549Zafirlukast300 mg/kgIIIInactiveInactiveInactive6DB00654Latanoprost50 mg/kgIIInactiveInactiveInactive7DB00944Demecarium6 mg/kgIIactiveInactiveInactive8DB00947Fulvestrant2000 mg/kgIVInactiveInactiveInactive9DB01100Pimozide464 mg/kgIVInactiveInactiveInactive10DB06249Arzoxifene2320 mg/kgVInactiveInactiveactive11DB06401Bazedoxifene200 mg/kgIIIInactiveInactiveInactive12DB06695Dabigatran etexilate200 mg/kgIIIInactiveInactiveInactive13DB08909Glycerol phenylbutyrate731 mg/kgIVInactiveInactiveInactive14DB08964Gemeprost57 mg/kgIIIInactiveInactiveInactive15DB09026Aliskiren4000 mg/kgVInactiveInactiveInactive16DB09030Vorapaxar9 mg/kgIIInactiveInactiveInactive17DB13615Mifamurtide1832 mg/kgIVInactiveInactiveInactive18DB14185Aripiprazole lauroxil800 mg/kgIVInactiveInactiveInactive

In-silico anticancer activity analysis

The pdCSM-cancer, a machine learning based online tool was used to predict the biological activity for 17 lead compounds against the most prevalent cancer types including lung, breast, colorectal, and prostate [36]. The pdCSM-cancer algorithm is trained and validated using data from human cancer cell lines (NCI-60). Here, the cancer cell lines such as A549 (lung), MCF7 (breast), HCT116 (colorectal) and PC3 (prostate) considered in the present analysis. The results were analysed in terms of the GI50% (growth inhibitory concentration by 50%) of the compounds. To benchmark our predictions, brequinar was used as a reference compound. The GI₅₀ values of brequinar were found to be 4.43, 4.41, 5.71, and 5.66 in cell lines respectively A549, MCF7, HCT116 and PC3. Notably, brequinar exhibited least activity against HCT116 as observed from its higher GI₅₀ value than studied cell lines. Interestingly, all the screened compounds demonstrated lower GI₅₀ values in HCT116 compared to brequinar, suggesting enhanced anticancer activity against colorectal cells. Zafirlukast and Glycerol Phenylbutyrate showed potency in lung, colorectal and prostate cancer cell lines. None of the screened compounds showed potency against breast cancer cell line, MCF7 in our analysis (Table 4). Thus, we have taken all the compounds for the subsequent interaction analysis.

Table 4. Anticancer prediction from PdCSM webserver for reference brequinar and 17 lead moleculesS NoDrug Bank IDNameLung (A549)Breast (MCF7)Colorectum (HCT116)Prostate (PC3)1DB03480Brequinar4.434.415.715.662DB00137Lutein4.444.55 4.45

4.63 3DB00481Raloxifene4.835.03 5.15

5.22 4DB00503Ritonavir4.535.44 4.49

5.58 5DB00549Zafirlukast 4.42 5.26 4.61

5.17 6DB00654Latanoprost4.824.97 4.59

4.90 7DB00944Demecarium4.635.26 5.15

5.38 8DB00947Fulvestrant4.745.33 5.35

5.07 9DB01100Pimozide4.815.27 5.07

5.07 10DB06249Arzoxifene5.085.02 5.36

5.29 11DB06401Bazedoxifene4.955.12 5.26

5.28 12DB06695Dabigatran etexilate4.655.60 5.22

5.36 13DB08909Glycerol phenylbutyrate 4.35 5.03 4.53

4.91 14DB08964Gemeprost4.525.25 4.79

5.19 15DB09026Aliskiren5.045.55 5.24

5.38 16DB09030Vorapaxar4.765.34 5.04

5.32 17DB13615Mifamurtide 4.16 5.40 5.02 6.4118DB14185Aripiprazole lauroxil4.815.42 5.22

5.17 The bold highlights the better anticancer activity of the compounds. The predicted GI₅₀ values are expressed in − log₁₀[M]

Interaction analysis of the protein-ligand complex

The post-molecular docking study was examined to comprehend the binding mechanism of the reference and 17 lead compounds in the active site of DHODH. Figure 2 represents the protein-ligand interaction pattern where hydrogen bond interactions are indicated by the blue, while hydrophobic interactions are indicated by the grey dashed lines [27]. The interaction profiling of DHODH-brequinar showed that it forms hydrogen bonds with GLN47 and THR360 residues, respectively. Besides, LEU46, GLN47, ALA55, HIS56, ALA59, PHE62, LEU68, TYR356, THR360 and PRO364 were able to form hydrophobic interactions. ARG136 forms a salt bridge with the carboxylate group of DHODH. It is likely that TYR356, provides the protons required to from a neutral hydroquinone following electron transfer from FMN [33, 34]. Although there are many binding interactions, hydrogen bonding is one of the key non-covalent interactions in stabilizing the protein-ligand complex [37]. Indeed, these patterns of interactions were explored for all the investigated lead compounds. Among which DB09026, DB00503 and DB13615 showed a high number of hydrogen bonds compared to other compounds. Hydrogen bonds have significant interactions because they can be crucial for drug partition and permeability, structural stability, molecular recognition, and enzyme catalysis. A drug’s solubility makes significant interactions with its biomolecular targets, which can be enhanced by the presence of functional groups that can form hydrogen bonds, promoting strong binding [38]. Further, DHODH-DB09026 was able to form nine hydrogen bonds with ARG136, ALA96, THR283, LYS100, LYS255, TYR356 and ASN145. The presence of TYR356 and ARG136 indicates a strong binding, conforming to the orientation of active site residues closer to the hit molecule. The binding pattern of DHODH-DB00503 involves interactions with ASN145, TYR356, LYS255, GLY335 and THR357, forming seven hydrogen bonds. Ten hydrophobic bonds form with ALA55, ALA59, VAL143, THR360, TYR356, TYR147, and PHE149. Primarily, TYR356 forms a higher number of hydrogen bonds, showing strong hydrogen bonding and facilitates the binding of small molecule with minimal strain of the hits.

Fig. 2. Interaction analysis of a) Reference - Brequinar; b) DB09026 Aliskiren and c)DB00503 Ritonavir towards the target protein, blue-coloured lines indicate hydrogen bond interactions and grey colour line signifies hydrophobic interactions

Although DB13615 forms two hydrogen bonds with ARG136, GLN47, and ASP51, due to its higher molecular weight DB13615 is not considered a hit compound. Higher molecular weight FDA-approved compounds can lead to issues such as reduced solubility, lower permeability, inconsistent pharmacokinetics and low oral bioavailability [39]. Overall, DB09026 and DB00503 showed a greater number of hydrogen and hydrophobic interactions and showed enhanced binding with DHODH protein.

Scaffold analysis

Past studies indicate that brequinar has a quinoline-4-carboxylic acid moiety as the biological activity and is present in different natural and synthetic products. Many secondary metabolites contain quinoline, a heterocyclic molecule made up of pyridine and benzene rings [40]. The screened hit aliskiren (DB09026) contains a phenyl group with 3-methoxypropoxy group showing antibacterial and antifungal activities that enhance the microbial activity [41]. Similarly, Aliskiren is a mono methoxybenzene-based molecule, characterized by the presence of a 3-methoxypropoxy group. Mono methoxybenzene scaffolds are well-recognized for their broad spectrum of bioactive properties, including antiviral and anticancer activities. Recent studies had reported that a mono methoxybenzene derivative isolated from Zingiber officinale exhibited significant anticancer and anti-inflammatory effects. These findings further support the hypothesis that Aliskiren, through its structural features and scaffold similarity, holds promise as a potential anticancer agent, justifying its repurposing for oncology applications [42]. Furthermore, this drug is used as a direct inhibitor to reduce blood pressure by inhibiting renin and is also used for the treatment of cancer cachexia [43, 44]. Similarly, thiazole moieties of Ritonavir have been validated for diverse activities including antiviral, anticancer and anti-inflammatory properties [45–47]. The 2D structures of reference and the hit compounds emphasizing its scaffolds were visualized in Fig. 3.

Fig. 3. Scaffold analysis of a) Reference-brequinar with Quinoline carboxylic acid moiety; b) DB09026 Aliskiren with 4-methoxy-3-(3-methoxypropoxy) phenyl moiety and c) DB00503 Ritonavir with Thiazole moiety

Molecular dynamic simulation

Conformational stability analysis

Molecular Dynamic simulations are performed to identify the conformational properties of protein-ligand complex [48]. RMSD analysis is performed to examine the stability profiles of the two lead compounds in complex with DHODH, as well as the DHODH- brequinar complex. The RMSD plots of DHODH-brequinar (black), DHODH-DB09026 (orange), DHODH-DB00503 (violet) complexes are given in Fig. 4. Brequinar (black) showed a deviation at 20 ns and 40 ns, stabilizing at 0.25 nm until 75 ns. Similarly, DB09026 (Orange) exhibited a 0.2 nm deviation at the beginning, with a maximum deviation of 0.25 nm observed at 75 ns. DB00503 exhibited the minimum deviation at 60 ns. At the end of 100 ns simulation, the average RMSD values for brequinar, DB09026, and DB00503 were found to be 0.25 nm, and 0.15 nm, respectively. Thus, the binding of hits resulted in the least conformational change.

Fig. 4RMSD profiles of the protein–ligand complexes are shown, with brequinar as the reference (black), Aliskiren (DB09026) in orange, and Ritonavir (DB00503) in violet

Residual flexibility analysis

The flexibility of protein-ligand complexes was evaluated using RMSF analysis. Figure 5 shows that the backbone fluctuation pattern of the protein-ligand complex. Notably, the binding residues (ARG136, GLN47, TYR356, THR360, LEU46, ALA55, HIS56, ALA59, PHE62, ASN145, LYS100, LEU68, LYS255 and PRO364) showed a lower average RMSF in DB09026 (0.12 nm) and DB00503 (0.11 nm) complexes than brequinar (0.15 nm) system at the end of a 100 ns. The research suggests that a higher RMSF number indicates less stability of the complex, whereas a lower RMSF value indicates greater stability [48]. Overall, the results indicate that the hits showed less fluctuation adapts the favourable position than brequinar at the binding site of the DHODH protein.

Fig. 5. Root mean square fluctuation (RMSF) profiles of the protein–ligand complexes. Brequinar (reference) is shown in black, Aliskiren (DB09026) in orange, and Ritonavir (DB00503) in violet

Hydrogen bond interaction

Molecular recognition depends on the specificity and directionality of interactions made possible by the hydrogen bonding connection between the protein-ligand complexes. The gmx hbond tool of GROMACS was used to find the total number of hydrogen bonds between the DHODH ligand complexes [49]. The hydrogen bond plot (Fig. 6) revealed that the DHODH-brequinar, forms 3 hydrogen bonds. Meanwhile, DHODH-DB09026 and DHODH-DB00503 observed to exhibit exhibits 5–6 hydrogen bond interaction. Based on these results, we can conclude that both DB09026 and DB00503 complexes actively forms hydrogen bonds with active site residues of DHODH protein.

Fig. 6. Hydrogen bond (H-bond) profiles of the protein–ligand complexes, with Brequinar (reference) represented in black, Aliskiren (DB09026) in orange, and Ritonavir (DB00503) in violet

Compactness evaluation

The compactness level of hits in the DHODH binding pocket were determined using degree of compactness (Rg). The weighted RMSD of each atom, as established by its centre of mass, is commonly denoted by Rg. Of note, an increase in Rg value indicates a decrease in protein structural compactness, leading to greater flexibility and instability [50]. Figure 7 shows the compactness of DHODH-brequinar, DHODH-DB09026 and DHODH-DB00503, which were found to be 1.93 nm, 1.94 nm and 1.92 nm, respectively, at the end of the 100 ns simulation. Finally, compactness evaluation reveals that hits adapts favourable conformation in the DHODH protein than reference molecule in our analysis.

Fig. 7. Radius of Gyration (Rg) profiles of the protein–ligand complexes. Brequinar (reference) is shown in black, Aliskiren (DB09026) in orange, and Ritonavir (DB00503) in violet

Solvent accessible surface area

The solvent-accessible surface area (SASA) method is used to calculate the surface area of polar and non-polar molecules in order to comprehend how residues interact with the surrounding solvent [15]. The SASA for the brequinar, DB09026, and DB00503 complexes is displayed to be 208 nm^2^, 209 nm^2^, and 208 nm^2^ respectively (Fig. 8). It is observed that both DB09026 and DB00503 complex with protein have stable conformational changes as of brequinar. There is no significant variation observed in terms of solvent accessible surface area analysis.

Fig. 8. Solvent accessible surface area (SASA) profiles of the protein–ligand complexes. Brequinar (reference) is depicted in black, Aliskiren (DB09026) in orange, and Ritonavir (DB00503) in violet

Discussion

The DHODH is the plays a vital role in the de novo pyrimidine biosynthesis pathway and mitochondrial function. The overexpression of DHODH contributes to carcinogenesis and thus is considered as potential target for the cancer therapeutics. Here, libraries of 2619 FDA-approved compounds were exploited for its activity against DHODH through high-throughput strategy. Notably, two tiers of docking and binding free energy analysis were implemented to improve the prediction accuracy. The results yielded 51 compounds with satisfactory docking and binding energy values. It is also note that van der Waals interaction provides adequate involvement in the binding of inhibitors. MM/GB(PB)SA, is one of the most popular method, achieved reasonable correlation with experimental affinities was employed to validate the accuracy of the docking results [21]. Notably, Uni-GBSA, the MM/GBSA implementation used in our study, has demonstrated a Pearson correlation coefficient of 0.64 and RMSE of 1.70 kcal/mol in comparison to experimental affinities, supporting its robustness for predicting ΔG_bind_. Such approaches help narrow down viable candidates for downstream experimental validation and therapeutic exploration [51]. Following this rescoring, 17 lead compounds were selected based on their projected binding free energy. The toxicity and anti-tumor activity assessments for these 17 lead compounds were also found to be in the acceptable range.

The 17 lead compounds were categorized based on their functions concerning the DHODH and their hydrogen bond interactions with the active residues. Compounds such as DB00654 (Latanoprost), DB00944 (Demecarium), DB01100 (Pimozide), DB14185 (Aripiprazole Lauroxil), DB00137 (Lutein), DB06695 (Dabigatran Etexilate), DB09030 (Vorapaxar), DB00481 (Raloxifene) and DB06401 (Bazedoxifene) fail to form hydrogen bond interactions with important residues including ARG136, GLN47, THR360, and TYR356 in the protein’s binding pocket. Hence these nine compounds are not taken to the further scaffold analysis. However, DB00947 (Fulvestrant), DB06249 (Arzoxifene), DB08909 (Glycerol Phenylbutyrate), DB08964 (Gemeprost) and DB00549 (Zafirlukast) are able to maintain interaction only with TYR356. Hence, the lack of involvement of other key residues limits their binding strength. Further, DB13615 (Mifamurtide) higher molecular weight significantly impacts the consideration of the compound for further analysis due to poor permeability and low solubility. Adding together, the pharmacokinetic and dynamic properties are also not satisfactory.

On the other hand, key residues such as ARG136 and TYR356 interact with aliskiren, demonstrating its potential to inhibit DHODH. Precisely, TYR356 forms two hydrophobic interactions, indicating stronger binding with the protein and suggesting its potential as a promising target for DHODH in cancer treatment. Furthermore, the binding pattern of DHODH-DB00503 complex reveals extensive molecular interactions. Specifically, it forms seven hydrogen bonds with ASN145, TYR356, LYS255, GLY335 and THR357. In addition, ten hydrophobic interactions are established with ALA55, ALA59, VAL143, THR360, TYR356, TYR147, and PHE149. Among these, TYR356 emerges as key residue due to its higher number of hydrogen bonds, contributing significantly to the stability of the complex.

Aliskiren (DB09026) modulates blood pressure through the RAAS system (renin-angiotensin-aldosterone system), a hormonal system involved in regulating blood pressure and fluid balance system. Interestingly, RAAS, is known to promote cell proliferation, hypertrophy, and inflammation by activating signaling pathways linked through various mechanisms, particularly by influencing cell proliferation, mitochondrial function and nucleotide metabolism [52]. Since DHODH is dependent on mitochondrial electron transport, RAAS activation may modulate pyrimidine synthesis through mitochondrial dysfunction or oxidative stress.

A HIV protease inhibitor called Ritonavir (DB00503) inhibits viral proteinase enzyme that normally cleaves the structural and replicative proteins that are produced by the main HIV genes, including gag and pol [47]. Usually, viruses replication depends on pyrimidine nucleotides for their genome replication. Ritonavir also has off-target effects on cellular metabolism, particularly mitochondrial pathways which could interact with DHODH mechanism through inhibition [53].

In recent studies shows that the ritonavir act as a pharmacokinetic booster due to its potent and irreversible inhibitory effects on the CYP3A4 and also it inhibits drug transporters such as breast cancer resistance protein [54]. Aliskiren also shows anti-atherosclerotic and cardioprotective effects. It also has ability to reduce the reactive oxygen species (ROS) production. Aliskiren reduced tumour necrosis factor-α and interleukin-6 levels and inhibited myeloperoxidase activity, providing anti-inflammatory actions. Aliskiren also demonstrated antiproliferative efficacy by lowering epidermal hyperplasia and proliferating nuclear antigen levels. Similarly, ritonavir combined with docetaxel have increased antitumour effects in prostate and in breast cancers. It also acts as the antiproliferative agent in the ovarian, lung, bladder and pancreatic cancers [55–58]. Ours is the first report on ritonavir inhibitory effects in colorectal cancer. All these information highlights the hit compounds potential to inhibit DHODH in cancer cells. We, further hypothesize that combining the hit molecules with DHODH inhibitors could amplify anticancer effects by restricting pyrimidine biosynthesis.

The potential side effects of ritonavir have seen in the epithelium kidney cells induces endoplasmic reticulum stress and mitochondria stress. Similarly, aliskiren is known to cause kidney injury due to its RAAS blockade by angiotensin receptor blockers [58]. Although reports suggest that both compounds may exhibit minimal side effect, it could be of interesting choice for the management of cancer due to its improved potency. Indeed, the key limitation is lies in the lack of experimental validation of these compounds, which remains essential and could serve as an interesting future direction towards the development of novel DHODH inhibitor.

Conclusion

The present study aimed to identify novel and potent DHODH inhibitors through virtual screening, with a distinct focus on repurposing. Initially, the binding poses of the target protein upon the ligand binding were analyzed using molecular docking and binding free energy calculations. Furthermore, the 17 lead compounds showed a higher binding free energy and was investigated for its toxicity, anti-cancer activity and scaffold analysis. Collective evidences from our analysis highlight satisfactory profile of DB09026 and DB00503 in all the investigation. The scaffold analysis reiterated the fact of anti-cancer activity of the identified compounds. Towards the end, the molecular dynamic simulation revealed comfortable positioning of the compounds in the binding pocket DHODH protein. However, experimental evidences on these compounds are certainly needed and likely to be an interesting future direction to confirm their activity against cancer cell lines.