Identification of potential inhibitors of dihydrofolate reductase (DHFR) through blocking the folate biosynthetic pathway of Mycobacterium tuberculosis utilizing structure-based virtual screening

Abstract

Tuberculosis (TB) has emerged as a leading cause of death due to a single infectious agent—Mycobacterium tuberculosis (Mt). This situation is exacerbated by delayed diagnosis, inadequate administration of effective TB medications, prolonged duration of treatment, shortage of toxin-free TB drugs, and frequent increases in resistance to most TB drugs. In an urge to find potential drug candidates for the treatment of fatal infectious TB disease, we targeted the folate biosynthetic pathway that involves the ubiquitous enzyme dihydrofolate reductase (DHFR), which catalyzes the NADPH-dependent reduction of dihydrofolate with the generation of tetrahydrofolate (THF). Blocking the enzymatic activity of DHFR exhausts the cellular pool of THF, which results in cessation of DNA synthesis in rapidly proliferating cells and ultimately cell death. Herein, a total of 1026 drug-like molecules with antibacterial activities were tested using several in silico tools for determining drug-likeness features, ADMET (absorption, distribution, metabolism, excretion, and toxicity) profiling, binding affinity, and conformation analysis using Autodock Vina and Schrodinger Suite. This exhaustive investigation identified CHEMBL577, CHEMBL161702, and CHEMBL1770248 as potential drug candidates for the inhibition of M. tuberculosis DHFR protein. Root mean square deviation, root mean square fluctuation, hydrogen bond, and MMGBSA evaluation by 100 ns molecular dynamics simulation (MDS) confirmed their molecular stability with the target protein. All of these drug-like compounds outperformed the control drugs trimethoprim and methotrexate in molecular docking and molecular dynamics simulation tests. Therefore, our study suggests these M. tuberculosis DHFR inhibitors as promising drug candidates. However, additional wet-lab experiments are required to verify their potential therapeutic potency as novel drugs against M. tuberculosis.

---

1 Introduction

According to the 2020 World Health Organization (WHO) report, tuberculosis (TB) has become a major menace to humanity, affecting 10 million people worldwide and leading to 1.5 million fatalities each year.¹ It is caused by Mycobacterium tuberculosis, which is a distinct Gram-positive bacterium with a unique composition and arrangement of the mycobacterial cell envelope that makes treating TB even more challenging.^2–4 Along with poor diagnosis, occasionally inadvertent delivery of TB drugs with toxic properties and prolonged chemotherapy with the existing TB drugs trigger the upsurge of multi-drug-resistant tuberculosis (MDR-TB) (resistant to at least isoniazid and rifampicin) and extensively drug-resistant tuberculosis (XDR-TB) (resistant to at least isoniazid and rifampicin, in addition to capreomycin, fluoroquinolone, amikacin or kanamycin among second-line anti-TB drugs).⁵ Consequently, new anti-TB drugs with no resistance or toxicity are desperately required to fight against MTB. Therefore, we evaluated the potential of and employed the enzyme dihydrofolate reductase as a target for the development of a suitable drug in this study.

Dihydrofolate reductase (DHFR, EC 1.5.1.3) is the most investigated, validated, and targeted member of all the enzymes that comprise the folate cycle. It is a key enzyme in the folate metabolic cycle that is found in all dividing prokaryotic and eukaryotic cells.^6–8 DHFR yields tetrahydrofolate (THF) by catalyzing the NADPH-dependent reduction of dihydrofolate. THF serves its function as a critical enzymatic cofactor that shuttles methyl and formyl groups from one molecule to another one in several one-carbon unit transfer reactions of various biosynthetic and degradative processes.^9–13 THF is necessary for the production of purine and thymidylate nucleotides, pantothenate, and amino acids (e.g., methionine and glycine-to-serine conversion), all of which are required for DNA, RNA, and protein synthesis.^9,11 Moreover, THF is required for the synthesis of formyl methionyl tRNAfMet, which is critical for initiating bacterial protein synthesis.¹¹ Since rapidly dividing cells are in high demand for DNA, RNA, and protein synthesis, blocking the enzymatic activity of DHFR is the basis for the design and development of new agents against DHFR.¹¹ As a result, the loss of DHFR activity depletes the cellular pool of the reduced folate cofactor, THF, and its derivatives, resulting in the cessation of DNA synthesis, cell growth and proliferation, and, finally, cell death. Because of its tiny size (18–22 kDa), availability of the pure form, and the presence of a small-molecule binding pocket with established biochemical characteristics, DHFR attracts researchers' interest as an attractive candidate for utilizing the scope as a therapeutic target.^14–16 As of now, a number of selective DHFR inhibitors are significant in the chemotherapy for several human diseases as anti-protozoal,^17–19 anti-bacterial,^20–22 anti-fungal,^23–25 and anti-cancer agents,^11,26,27 and against psoriasis, autoimmune²⁸ and neoplastic diseases as well.¹¹ For example, methotrexate (MTX) is one of the most widely known and potent inhibitors that bind to both human DHFR (h-DHFR) and Mycobacterium tuberculosis DHFR (mt-DHFR) with no discernible selectivity.^29,30 Both trimethoprim (TMP) and pyrimethamine (PYR) are effective inhibitors of bacterial and protozoal DHFR, respectively, however not potent either for mammalian or mt-DHFRs.^30–32 Although TMP has a low affinity for mt-DHFR,³⁰ the synergistic effect of TMP/sulfamethoxazole DFHR inhibitor medication makes Mtb susceptible.³³ WR99210, another potent DHFR inhibitor, functions both as an anti-malarial agent and as an anti-tubercular agent that has had its clinical development halted due to unacceptable side effects.^34,35

X-ray crystallographic studies unveiled the overall similarities between M. tuberculosis DHFR (mtDHFR) and human DHFR (hDHFR) with a sequence homology of 26%.^36,37 We used the wild-type DHFR (PDB ID: 1DG5), as there are no experimentally resolved mutant DHFR structures of M. tuberculosis currently deposited in the Protein Data Bank. Their high-resolution crystal structures propose some of the crucial structural differences between host and pathogen DHFR enzymes (Fig. 1). Both types of DHFR enzymes exhibit the same general fold, even though hDHFR is significantly larger with 187 amino acid residues compared with 159 for the mtDHFR.^37,38 DHFR is a small protein with an α/β structure folded into a central eight-stranded β-sheet of seven parallel strands with a C-terminal antiparallel strand and four flanking α-helices.³⁸ Most of the structural characteristics are identical in pathogen and host DHFR, which occurs at the C-terminal of the β sheet, whereas structural alterations occur at the N-terminal of the β-sheet. A striking difference located in the proximity of the active region is the extension of a loop domain that puts an end to the interactions between this loop and the loop domain found in hDHFR. This alignment of the hDHFR loop domain is equivalent to the “Met 20 loop” in the E. coli DHFR, which influences the accessibility of the active region.^39,40 Another noticeable difference found from the analysis of the crystal structures of MTX binding to mtDHFR and hDHFR implies that mtDHFR combines with MTX, glycerol (GOL), and NADPH (PDB ID: 1DF7).³⁷ The sequence homology is approximately 55% in the active and ligand binding site than the sequence homology of 17% for the rest of the chain in mtDHFR.^37,39,40

Fig. 1 X-ray crystallographic structure of dihydrofolate reductase of Mycobacterium tuberculosis (mt-DHFR) (PDB: 1DG5). (A) Helix, sheet, loop, NADPH, and native ligand trimethoprim shown in pink, yellow, green, blue, and magenta color, respectively, in the three-dimensional structure. (B) Major amino acid involved in the interaction pattern of the native ligand.

Our in silico technique aims to alleviate the load of the prevailing antibiotic resistance situation by investigating novel therapeutic usage. A total of 1026 drug-like substances were obtained as part of this study to evaluate their efficacy in the treatment of tuberculosis. This investigation followed by measuring the binding energy and binding interaction of 1026 compounds against the dihydrofolate reductase of M. tuberculosis (PDB ID: 1DG5) was performed using two distinct docking programs, Autodock Vina and Maestro Schrodinger Suite. The characterization of pharmacokinetics profiling of retrieved compounds helped to validate their drug-like features and toxic effects. Lastly, structural dynamics simulation and MM-GBSA analysis were performed to assess the structural stability and binding energy of the lead drug-like molecules with their particular binding proteins (Fig. 2).

Fig. 2 Complete methodology of this study in a concise flowchart.

2 Materials and methods

2.1. Protein and ligand preparation

The three-dimensional structure of the receptor dihydrofolate reductase of M. tuberculosis in complex with NADPH and trimethoprim (PDB ID: 1DG5)³⁷ was obtained from the RCSB Protein Data Bank⁴¹ (https://www.rcsb.org/search). All 1026 ligands possessing antibacterial activity have been collected from the ChEMBL Database⁴² (https://www.ebi.ac.uk/chembl/). Each of these selected ligands was converted into the PDBQT format employing Open Babel of the PyRx software.⁴³ Trimethoprim was used as a control drug against dihydrofolate reductase enzyme. Prior to docking analysis, the PDB protein was prepared by subtracting the water molecules from the structure employing the PyMOL program,⁴⁴ and then refined and optimized using Swiss-PdbViewer.⁴⁵

2.2. Molecular docking

Molecular docking is a crucial part of rational drug development,⁴⁶ and it was used in this study to evaluate the binding interaction and affinity of all the ligands with the dihydrofolate reductase. A total of 1026 ligands were docked with the enzyme employing Autodock Vina 4.2.⁴⁷ This tool estimates the binding energy between ligands and the protein employing the Lamarckian Genetic Algorithm. Following protein and ligand preparation, AutoDockTools-1.5.6rc3 was used to convert them into the vina-compliant PDBQT format.⁴⁸ Grid boxes with required dimensions were generated applying the same tool. The center of the grid box for the dihydrofolate reductase (PDB ID: 1DG5) was settled at Center_x = 8.69, center_y = 19.838, Center_z = 16.301 and size_x = 86, size_y = 82, size_z = 70, exhaustion = 10.

2.3. Evaluation of drug-likeness properties

Drug-likeness is the structural or physicochemical properties of a drug-like molecule. A drug-like molecule can be considered as a drug candidate assessing the drug likeness properties. The freely accessible Swiss-ADME server (http://www.swissadme.ch/index.php) was used in this study to find out the key physicochemical, pharmacokinetic, drug-like, and related parameters of the filtered 452 drug-like compounds.⁴⁹ This tool followed Lipinski's rule of five criteria,⁵⁰ Ghose's rule,⁵¹ Veber's rule,⁵² Muegge's rule,⁵³ TPSA, and the number of rotatable bonds to evaluate the properties.

2.4. ADMET prediction

The evaluation of chemical absorption, distribution, metabolism, excretion, and toxicity (ADMET) of the ligand is inevitable before designing a final drug. ADMET profile of 248 ligands was performed in admetSAR (http://lmmd.ecust.edu.cn/admetsar2/).^54,55 The ChEMBL database (https://www.ebi.ac.uk/chembl/) was used to obtain the canonical smiles of the ligands, which were then used as input data for the study of drug-likeness properties.⁴²

2.5. Glide docking and MM-GBSA analysis

A total of 30 ligands were selected from previous analyses. Then they were prepared and refined using the LigPrep option of Maestro Schrodinger Suite.⁵⁶ Precise and energy-minimized 3D molecular structures of those ligands have been produced utilizing Epik2.2 and sustaining pH 7.0 ± 2. Energy minimization of those ligands was carried out exploiting the OPLS3e force field.⁵⁷ The protein DHFR was preprocessed once again using the Protein Preparation module of Maestro Schrodinger.⁵⁸ Hydrogen bonds and bond ordering were distributed, and missing side chains and loops were added to the protein employing the protein Preparation module. Before removing water beyond 3.0 Å, protein structures were optimized using PROPKA (pH 7). Subsequently, the OPLS3e force field was applied for energy minimization.

One of the essential parts of the docking study is defining grid box as it guides the drug-like molecule to the protein's ligand-binding region. Receptor Grid Generation module was used to acquire the grid map utilizing baseline parameters such as keeping a van der Waals radius scaling factor of 1.0 and a charge cutoff value of 0.25. XP Glide methodology was performed utilizing the Ligand docking module of Maestro Schrodinger Suite. XP (extra precision) mode of docking using Glide results in reproducible experimental binding affinity and accurate binding poses. A cutoff score of 0.15 and van der Waals radius scaling factor of 0.80 were assigned for ligand atoms. The glide docking score was determined based on ligand-protein binding complexes within the binding sites of receptors.

The free binding energy of each ligand with the DHFR protein was computed employing the Prime MM-GBSA module of Maestro Schrodinger Suite. The module integrates the OPLS3e force field, VSGB solvation model to calculate the MM-GBSA scores using protein-ligand complexes.⁵⁸

2.6. Evaluation of interaction patterns and visualization

BIOVIA Discovery Studio 4.1 Visualizer was implied to identify and visualize multiple non-bonded interactions within each ligand-protein docked complex.⁵⁹ The hydrogen-bond interaction, hydrophobic contacts between the ligands, and electrostatic interaction patterns among the amino acid of the protein and the ligand are depicted in the 2D and 3D schematics of the docked complex.

2.7. Molecular dynamics simulation

The MD simulation was performed in Schrodinger's Desmond module.⁶⁰ Using the System Builder tool, the best 3 ligands and control Trimethoprim-Mtb dihydrofolate reductase complexes were positioned in the orthorhombic box with a buffer distance of 10 Å in order to add water, and single point charge (SPC) was used to generate a water model. The systems were neutralized by the addition of Na+ and Cl-ions until a 0.15 M concentration was reached. The built solvated system was minimized and relaxed before the production run of the simulation using OPLS3e force field parameters as the default protocol associated with Desmond. The MD simulation was run using an isothermal isobaric ensemble (NPT) at 300° K temperature and a pressure of 1.013 bar. A 100-nanosecond simulation was run during which 1000 frames were saved to the trajectory. Lastly, using the Simulation Interaction Diagram (SID) tool, MD simulation trajectory was analyzed.^61–66

2.8. Density functional theory (DFT) calculations

The electrical structures of ligands are also linked to their pharmacological activities. It is important to examine the structural behavior of the active compound and to explore how structural orientation, any biological effect part in the structure, and what parameters may strength the biological activeness of the molecule. As a result, single-point energy calculations using DFT were performed in this study to explore the detailed aspects in terms of structure, electronics, and energy states of every atom of compound. Frontier molecular orbitals (FMO) of compounds specify the critical role of charge transfer interactions with the binding site of M. tuberculosis DHFR. The HOMO (highest occupied molecular orbital) and LUMO (lowest unoccupied molecular orbital) are commonly known as FMO and were found to give extremely applicable information about electron density clouds around the molecule. To shed light on the electronic structural properties of virtual screened compounds, both the screened compounds and the control drug trimethoprim were entered into the Jaguar platform in Schrodinger to compute the HOMO and LUMO by incorporating the basis set at the 6-31G* level and hybrid DFT with Becke 3-parameter exchange potential.^67,68

3 Results

3.1. Molecular docking analysis

A total of 1026 ligands were screened in the AutoDock-Vina application and the most desirable candidates were filtered out. Based on the bound conformations and the binding affinity of all the ligands with the protein, suitable ligands are chosen. Given that the lower binding energy indicates a higher affinity, the ligands with the lowest binding energy are considered to be the most suitable drug candidate. The binding energy of the control drug (trimethoprim) was kept at −12 kcal mol⁻¹ to narrow down our dataset. Among all 1026 ligands, 452 exhibited lower binding scores than the control binding energy. The results of the comparative docking from the Autodock-Vina tool are provided in Table S1.

3.2. Drug likeness

Drug-likeness study determines the likelihood of a compound to be the best drug candidate analyzing its structural or physicochemical properties. The SwissADME online tool identifies compounds that meet the criteria for an optimal oral medication using five distinct pharmacokinetic rules (Lipinski, Ghose, Veber, Egan, and Muegge) based on different algorithms. A total of 452 drug-like compounds were examined for drug-likeness properties in our study.

Initially, all of them were refined following Lipinski's rule of five filter criteria, including molecular weight (optimal range: <500), amount of hydrogen bond donors (optimal range: ≤5), amount of hydrogen bond acceptors (optimal range: ≤10), lipophilicity (expressed as Log P, standard range: <5) and molar refractivity (optimal range: 40–130). Only the ligands that fit within these standard ranges were selected. Afterward, they were reanalyzed on the basis of the Ghose, Veber, Egan, and Muegge rules. The number of bioavailability scores, rotatable bonds, topological polar surface area (Standard range: 20–130 Å2), log S all fell within these rules, and ligands having the standard values confirmed zero violation of all these rules. Finally, a total of 248 compounds were chosen for further analyses. The Drug-Likeness properties of the selected ligands are reported in Table S2.

3.3. ADMET (absorption, distribution, metabolism, excretion, and toxicity) analysis

AdmetSAR is an effective platform for the prediction of comparative ADMET (chemical absorption, distribution, metabolism, excretion, and toxicity) profiles of screened ligands. For the molecule to be analyzed, the site uses SMILES format as the input data. The ADME/T test results, including probability score, are reported in Table S3.

In the absorption part, almost all ligands exhibit a high intestinal absorption rate, oral bioavailability, and positive Caco 2 permeability. Only the ligands with the most crucial value in comparison to others were chosen. In terms of distribution, however, not all of them are permeable to the blood–brain barrier (BBB). Following this, when it comes to distribution and excretion, none of the ligands inhibit any CYP450 enzymes, although a few of them act as substrates for CYP3A4, CYP2C9, and CYP2D6, while others have no function as a substrate for any of the CYP450 enzymes. Finally, all the screened ligands showed negative results in carcinogenicity and Ames mutagenicity test. They also have no acute oral toxic or hepatotoxic activity. A total of 30 drug-like molecules were selected following AdmetSAR analysis (Table S4).

3.4. Glide docking and MM-GBSA analysis

All 30 ligands were redocked, and the best candidates were chosen using the XP (extreme precision) mode. These ligands were also subjected to MM-GBSA analysis. Ligands with the highest XP glide and MM-GBSA values were considered to have the highest binding affinity. The glide scores for CHEMBL1770248, CHEMBL161702, and CHEMBL577 were −8.863, −9.62, and −8.813, respectively, whereas the control ligand trimethoprim and methotrexate displayed the lowest binding affinity of −7.08 and −7.797 (Table 1). However, MM-GBSA assessment revealed the comparative free binding energy of these ligands. The MM-GBSA scores of CHEMBL1770248, CHEMBL161702, and CHEMBL577 drug-like compounds were −71.98, −53.09, and −58.91 kcal mol⁻¹, respectively). These three drug-like compounds were further chosen for molecular dynamic simulation. The schematic representation of protein-ligand docked complexes is shown in Fig. 3.

Fig. 3 The non-bonded interaction of the top three docked complexes, where, A, B, and C indicates the interaction diagram of CHEMBL577, CHEMBL161702, and CHEMBL1770248, respectively. The first column depicts the surface view of the protein–ligand complex, with the ligand shown in red. The second column presents the pose view, illustrating ligand orientation and hydrogen bonding (red dashed lines) or hydrophobic contacts (yellow dashed lines) within the active site. The third column displays 2D receptor–ligand interaction maps, highlighting hydrogen bonds (green dashed lines) and hydrophobic interactions (pink dashed lines) with key amino acid residues.

3.5. Molecular dynamics simulations

Following the docking calculations, the best 3 ligands, and control Trimethoprim and methotrexate -Mtb dihydrofolate reductase complexes, MD simulations were carried out using the Desmond Simulation Package for 100 ns. Root–mean–square deviation (RMSD) values, root–mean–square fluctuation (RMSF) values, and protein–ligand contacts were calculated from the MD trajectories. The solvated system of the ligand–protein docked complex was checked for its binding stability based on the RMSD fluctuations during the simulation. The RMSD fluctuation is measured individually for the protein and ligand structures in the trajectory of MD simulations; if it falls within 3 Å, then the complex is considered to be stable.^{61,65,69–71} The RMSD values of the protein–ligand complexes of trimethoprim, methotrexate, CHEMBL577, CHEMBL161702, and CHEMBL1770248 are shown in Fig. 4. The protein–ligand complexes of trimethoprim, methotrexate, CHEMBL577, CHEMBL161702, and CHEMBL1770248 displayed average RMSD values of 2.35 Å, 5.44 Å, 7.93 Å, 3.39 Å, and 4.39 Å, respectively.

Fig. 4 Ligand RMSD profiles over 100 ns molecular dynamics simulations. RMSD values of fit ligand CHEMBL577 (green), CHEMBL161702 (orange), CHEMBL1770248 (gray) and control trimethoprim (gold) and methotrexate (blue) of the protein-ligand docked complexes.

The RMSD of the CHEMBL577 gradually increased for 13 ns to 31 ns time span; after that, the curve showed minor RMSD fluctuation, signifying the stability of CHEMBL577. The RMSD plot of CHEMBL161702 shows minor fluctuation from the initial stimulation period of time to 25 ns, after which a similar RMSD trend was observed with the control complex. In the CHEMBL1770248-complex, ligand RMSD increased substantively for 55 ns and then stabilized RMSD between 1.2 Å and 1.8 Å for the rest of the simulation. The RMSD plot of proteins is displayed in Fig. 5. All proteins initially show RMSD from 1.2 Å to 1.8 Å, after which the RMSD value was less than that of the control protein RMSD for the remainder of the simulation. Lower RMSD values indicate higher stability, whereas larger deviations reflect increased conformational flexibility or displacement from the initial binding pose.

Fig. 5 RMSD values extracted from Cα of the protein–ligand docked complexes, viz. CHEMBL577 (green), CHEMBL161702 (orange), CHEMBL1770248 (gray), trimethoprim (gold), and methotrexate (blue), with respect to 100 ns simulation time.

The RMSF values for Cα atoms of all residues were assessed using 100 ns trajectory data to examine the binding efficiency of lead compounds with Mtb Dihydrofolate reductase. The average RMSF values assessed for dihydrofolate reductase upon binding of trimethoprim, methotrexate, CHEMBL577, CHEMBL161702, and CHEMBL1770248 are 1.5 Å, indicating that all residues fluctuated within the RMSF range of 0.5 to 2.0 (Fig. 6), implying that Mtb dihydrofolate reductase exhibits minimal fluctuation and relative secondary conformational stability upon binding of reported lead compounds. All the RMSF values in the binding pocket are listed in Table 2.

Fig. 6 RMSF analysis of Cα atoms for protein–ligand complexes. RMSF profiles for docked complexes, viz. CHEMBL577 (green), CHEMBL161702 (orange), CHEMBL1770248 (gray), trimethoprim (gold), and methotrexate (blue), show residue-wise flexibility during 100 ns simulation time. Peaks around residues ∼20, ∼40, ∼75, ∼95, and ∼160 indicate higher flexibility, mainly in the loop and terminal regions, with slight variations among ligands.

The average total energies of trimethoprim-protein complexes and methotrexate–protein complexes were 9100 and 9199 kcal mol⁻¹, respectively. For CHEMBL577-protein and CHEMBL161702-protein complexes, the average total energy was below 9100 kcal mol⁻¹. However, the CHEMBL1770248-protein complex showed a slightly higher average energy (9206 kcal mol⁻¹) than the control drug complexes. After 80 ns, the total energy of the selected drug-complexes decreased and went below 9100 kcal mol⁻¹. It showed a slight deviation of total energy among the complexes of desired drug-like molecules and control drugs.

Throughout the 100 ns simulation, protein interactions with the ligand were also monitored. These interactions can be categorized by hydrogen bonds, hydrophobic, ionic, and water bridges, as illustrated in Fig. 7 and 8. The active site amino acid residues of the Mtb-dihydrofolate reductase contributing to binding interactions with trimethoprim were Trp22 (hydrogen bond donor and water mediated H-bond), Ala7 (hydrogen bond acceptor), Gly97 (hydrogen bond donor), and Thr46 (hydrogen bond donor) produced 50%, 32%, 69%, 49%, and 49%, respectively. The amino acid residues Arg45, Gly97, and Gly96 exhibited hydrogen bond contact with the molecule CHEMBL577 with 43%, 66%, and 66%, respectively, interaction in the simulation period. CHEMBL577 also interacted with the hydrophilic residue Arg45 via ionic interactions (Fig. 7A).

Fig. 7 Detailed atomic interaction of three best hits with the protein residues (PDB: 1DG5). The left side of the figure represents protein–ligand contact of the (A) CHEMBL577 (B) CHEMBL161702, and (C) CHEMBL1770248 compounds to the total 100 ns MD simulation time. The right side of the figure depicts the 2D interactions of the ligand–protein complex of 100 ns simulation time.

Fig. 8 Detailed atomic interaction of the control ligand: (A) trimethoprim, and (B) methotrexate, with the protein residues (PDB: 1DG5). The left side of the figure represents protein–ligand contact to the total 100 ns MD simulation time. The right side of the figure depicts the 2D interactions of the ligand–protein complex of 100 ns simulation time.

In the CHEMBL161702-dihydrofolate reductase complex, crucial hydrogen bond interaction was found with Arg45, Thr46, Ser49, Ile94, Gly96, Gly97, and Tyr100 over the course of the simulation. The amino acid residues Gly18, Asp19, Ile20, Arg23, Asp27, Gln28, Thr46, and Ser49 were involved in the H-bond formation with the molecule CHEMBL1770248 in the active site of 1DG5 (Fig. 7). All the amino acid interactions identified upon docking studies of the target molecule were also displayed during the dynamic study. We also studied the interaction pattern of the three selected ligands after simulation (last frame of simulation), as displayed in Fig. S1–S3. This demonstrated that the protein-ligand complex remained stable throughout the simulation period, with only minimal backbone fluctuations in the system. Overall, the simulation indicated more water-mediated hydrogen bond interactions with the lead molecules also found in the MD simulation study.

3.6. Post-dynamic MMGBSA

To execute the post-dynamic MM-GBSA evaluation, the thermal_MMGBSA.py script of the Prime/Desmond module of the Schrodinger suite was applied. The post dynamic-MM-GBSA analysis of free binding energy evaluation was computed with the creation of 1 frame with 1000-step sampling size. The frame was processed and analyzed throughout the post-dynamic-MM-GBSA calculation of 100 ns MD data of control compound trimethoprim, and lead compounds CHEMBL577, CHEMBL161702 and CHEMBL1770248, as revealed by the studies of the dynamics. The post-MD MMGBSA binding free energies revealed notable differences in ligand ranking and stability profiles compared to pre-dynamic calculations. The calculated binding free energy (ΔG) average of the trimethoprim, CHEMBL577, CHEMBL161702, and CHEMBL1770248 was found to be −53.5115 ± 8.01 kcal mol⁻¹, −42.6578 ± 4.54 kcal mol⁻¹, −52.0844 ± 0.79 kcal mol⁻¹, and −68.1447 ± 1.44 kcal mol⁻¹. Compared to the post-dynamic MM-GBSA calculations, the pre-dynamic MM-GBSA evaluation showed a similar binding free energy score (Table 3). Trimethoprim's post-dynamics MMGBSA score (−53.51 ± 8.01 kcal mol⁻¹) was slightly reduced from its pre-dynamics value, suggesting minor destabilization but retained strong binding. CHEMBL577 showed a larger drop (−42.66 ± 4.54 kcal mol⁻¹), indicating weaker dynamic stability. CHEMBL161702 remained highly stable (−52.08 ± 0.79 kcal mol⁻¹), while CHEMBL1770248 had the most favorable score (−68.14 ± 1.44 kcal mol⁻¹), surpassing trimethoprim and indicating superior stability and potential efficacy.

3.7. Density functional theory (DFT) calculations

The HOMO and LUMO orbitals govern how the molecule interacts with other species and help to explain its chemical reactivity and kinetic stability. HOMO and LUMO orbitals play a crucial role in charge transfer between these orbitals during a chemical process. HOMO energy represents the ability to donate an electron (electron-rich orbital) and the LUMO energy is related to accepting an electron (electron deficient orbital). The HOMO energy is proportional to the ionization potential, while the LUMO energy is proportional to the electron affinity. Furthermore, the energy gap (E) between the HOMO and LUMO orbitals is linked to interaction stability. It is a useful tool for determining the most bioactive compound.^64,72 Charge transfer occurs more easily within the molecule if the compound has a smaller energy gap. It is correlated that when the energy gap is wide, the system achieves low reactivity, and when the energy gap is low, the system achieves high reactivity. According to DFT result analysis, the energy gap of compounds CHEMBL577 (−0.2157), CHEMBL161702 (−0.1854), and CHEMBL1770248 (−0.1954) is closer to the control drug trimethoprim (−0.1926), which indicates their similar reactivity with trimethoprim (Fig. 9). In all three investigated compounds—CHEMBL577, CHEMBL161702, and CHEMBL1770248—polar functional groups such as amino (–NH₂), hydroxyl (–OH), and carboxyl (–COOH) were identified as key contributors to hydrogen bonding interactions with the active site residues of the Mtb-dihydrofolate reductase enzyme. These groups acted as both hydrogen bond donors and acceptors, forming stable interactions with amino acid residues such as Arg45, Gly97, Thr46, and Ser49, as observed in MD simulations. The DFT-calculated frontier molecular orbitals (HOMO and LUMO) were predominantly localized around these polar interaction sites. Specifically, the HOMO orbitals were primarily concentrated on the phenyl ring, suggesting that these regions have a higher electron-donating capability and may participate in π–π stacking or hydrogen bonding through lone pair donation. In contrast, the LUMO orbitals were mainly localized around the carboxyl and amino groups, highlighting their potential as electron-accepting regions, which facilitate the formation of hydrogen bonds with electron-rich amino acid residues in the protein's binding pocket.

Fig. 9 Plots of the HOMO, LUMO, and 3D-MESP overlaid onto a surface of constant electron density of CHEMBL577, CHEMBL161702, CHEMBL1770248, and trimethoprim.

The molecular electrostatic potential (MEP) maps further reinforced these observations by identifying electron-rich (red) and electron-deficient (blue) regions on the molecular surface.⁷² The most negative electrostatic potential was observed in CHEMBL577 (−57.46 kcal mol⁻¹) and CHEMBL161702 (−51.30 kcal mol⁻¹), indicating strong regions for electrophilic attack. These negative zones were located near the oxygen atoms of carboxyl and hydroxyl groups, while the most positive potentials were found around the nitrogen atoms of the pyrimidine or aminopyridine rings, suggesting favorable sites for nucleophilic interactions. In conclusion, the DFT-derived electronic parameters complement MD simulation results by confirming the significance of the polar functional groups in mediating ligand–protein interactions. The spatial localization of HOMO and LUMO orbitals around these active regions helps explain the observed stability and binding affinity of the ligand-protein complexes, thereby strengthening the overall mechanistic understanding of their bioactivity.

4 Discussion

Since the introduction of the computer-assisted molecular docking technique, a new era of drug discovery has begun. Because of multi-drug resistance mechanisms, therapeutic options for Mycobacterium tuberculosis are restricted. Consequently, there is an absolute need to look for novel molecular compounds and drugs for M. tuberculosis to fight this life-threatening disease.

Dihydrofolate reductase (DHFR) is a catalytic enzyme responsible for converting dihydrofolate to tetrahydrofolate. This is a crucial step in the synthesis of DNA, RNA, and proteins in M. tuberculosis.²⁰ Compounds that inhibit mycobacterial DHFR (mt-DHFR) might be effective in the treatment of tuberculosis. A range of 1026 ligands were evaluated in this investigation, with trimethoprim and methotrexate acting as the control drug since they all show inhibitory action against dihydrofolate reductase (DHFR).^10,30,31,33 To proceed with the molecular docking interactions, all the ligands were docked separately against DHFR using the PyRx tools AutoDock Vina. This section of the research yielded 452 ligands with the lowest binding energy when compared to trimethoprim and methotrexate.

The idea of drug-likeness, based on the interpretations of the physicochemical and structural attributes of small drug-like compounds, has been widely performed in drug designing to filter out compounds with properties consistent with accepted pharmacokinetic levels.^50,69,73 To screen out the best oral drug candidates, the pharmaceutical industry authorized 5 distinct ruled-based filters to be applied sequentially to all of the evaluated ligands.⁴⁹ The initial Lipinski filter (Pfizer) is based on some basic physicochemical features, and it rigorously values all nitrogen and oxygen as H bond recipients and all nitrogen and oxygen that have a minimum of one hydrogen as H bond providers.⁵⁰ In the lipid bilayer membrane, higher molecular weight is linked to a lower penetration rate, and Log P below 5 is 90% highly probable to be orally bioavailable.⁷⁴ When a substance fits with all the five basic rules, it will have improved pharmacokinetic characteristics and bioavailability in the metabolism.⁷⁵ To show a drug resemblance, all compounds must abide by the rules of Ghose, Veber, and Muegge as well.^51–53 In this experiment, we have carefully chosen only the ligands which showed a null gross violation of all five rules. The bioavailability scores, topological polar surface area (TPSA), and a number of rotatable bonds were also taken into consideration before choosing the final compounds. At the completion of the Swiss-ADME analysis, 248 ligands were picked for further investigation.

The pharmacokinetic profile of a drug-like molecule must be evaluated before discovering a drug to ensure that the drug can survive itself in the body, accomplish its function, and be eliminated effectively, which is termed ADMET (Absorption, Distribution, Metabolism, Elimination, Toxicity) properties. The absorption of a drug compound through the human intestinal cell lining is a fundamental way to determine the bioavailability of a compound after systemic administration.⁷⁶ As a result, the ADMET profile must include the human intestinal absorption (HIA) rate and Caco-2 permeability estimation. Our selected ligands show a high rate of HIA and all of them penetrate the Caco-2 cell line. While the ADMET profiling provided preliminary insights into the pharmacokinetic properties of the compounds; however, these predictions are based solely on computational analyses, and validation requires further in vitro and in vivo assays.

P-glycoprotein is an important transport protein that facilitates the transportation of many drugs, hence it is crucial to see whether the compound is a P-glycoprotein substrate or inhibitor.⁷⁷ The majority of drugs that function in the central nervous system (CNS) require BBB penetration.⁷⁸ Following absorption, the drug must be metabolized by an enzyme in the cytochrome P450 (CYP) family. CYP1A2, CYP2A6, CYP2C9, CYP2C19, CYP2D6, CYP2E1, and CYP3A4 are the principal isoforms of human CYP involved in medication metabolism.⁷⁹ Its catalytic activity fluctuates for different compounds, which is essential in changing the bioavailability and drug–drug interactions.⁸⁰ We focused on five of these CYP enzymes in this research: CYP1A2, CYP3A4, CYP2C9, CYP2C19, and CYP2D6. A compound that serves as a substrate for a CYP450 enzyme can readily be metabolized by that specific enzyme, and when it inhibits that enzyme, that enzyme will no longer be able to metabolize that compound.⁸¹

While designing a drug, various toxicity tests should be performed, including hepatotoxicity, acute oral toxicity, carcinogenicity, and Ames mutagenesis. Hepatotoxicity associated with drugs is increasingly the leading cause of acute liver failure, liver transplantation, and, in the worst-case scenario, death.⁸² The Ames mutagenicity assay examines whether a medication has mutagenic potential, which causes gene mutation by causing genetic harm.⁸³ Carcinogenicity and oral toxicity both are crucial endpoints that must be considered during the drug development process.^84,85 We screened out the ligands that do not have any of these toxicities. Out of 248 ligands, 30 drug-like molecules were chosen based on the ADMET properties.

The chosen ligands were then docked one more time with Maestro Schrodinger Suite. To estimate the ligand-binding affinity and compare the free binding energy of those ligands, the ligands were tested for the Glide XP and the MMGBSA score. The protein-ligand docked constructs were used to conduct the MM-GBSA derived binding free energy estimations.^86,87 CHEMBL577, CHEMBL161702, and CHEMBL1770248 compounds showed the highest Glide XP and the MMGBSA scores.

Molecular dynamics simulation was used to validate the conformational stability of putative drugs following interaction with dihydrofolate reductase. The RMSD of ligands trimethoprim, methotrexate, CHEMBL577, CHEMBL161702, and CHEMBL1770248 showed a minor RMSD fluctuation throughout 100 ns simulation. The binding of the ligands CHEMBL161702 and CHEMBL1770248 did not produce conformational instability within proteins since the measured RMSD change after ligand binding followed the same pattern as the trimethoprim-protein complex. A comprehensive examination of the root mean square fluctuation (RMSF) curve of Cα atoms indicated that our lead drugs maintained tight contact with their binding sites, as seen by their short fluctuation of 0.5 to 2.0 within dihydrofolate reductase. However, because of loop sections on dihydrofolate reductase, RMSF fluctuation was conspicuous. The 100 ns simulation of protein ligand interactions revealed that vGly18, Asp19, Ile20, Arg23, Asp27, Gln28, Thr46, and Ser49 residues were involved in the formation of H-bonds with the molecule CHEMBL1770248 at the active site of 1DG5. The total energy difference of CHEMBL577, CHEMBL161702, and CHEMBL1770248 complexes displayed a similar energy pattern compared to the control drug complexes along the 100 ns simulation. For the CHEMBL161702-dihydrofolate reductase complex, these H-bond-forming residues were Arg45, Thr46, Ser49, Ile94, Gly96, Gly97, and Tyr100, while Arg45, Gly97, and Gly96 amino acid residues form hydrogen bonds with the compound CHEMBL577. The protein-ligand association remained relatively stable throughout the 100 ns simulation, with minor backbone fluctuations in the system. Post-dynamic MMGBSA analysis of control drug trimethoprim and best hits CHEMBL161702, and CHEMBL1770248 also revealed their robust binding free energy, indicating superior stability and binding strength, making them promising candidates for further investigation. On the basis of our findings, dihydrofolate reductase appears to be a viable target for the CHEMBL161702 and CHEMBL1770248 molecules.

Lead compounds CHEMBL1770248, CHEMBL161702, and CHEMBL577 performed well and were considered to be the most efficient inhibitor of DHFR. These compounds CHEMBL577, CHEMBL161702, and CHEMBL1770248 are commercially available and approved, enabling direct procurement for in vitro assays in future studies (https://www.ebi.ac.uk/chembl/). Nevertheless, these ligands need to be studied further in the lab experiment to ensure their potency and inhibitory efficacy in vitro studies. Despite some limitations in terms of the quality of data obtained, distortions, and errors in the existing base of knowledge, computational methods can quickly provide a list of promising drug candidates in a cost-effective and timely manner for an organism like mycobacteria, which is difficult to treat owing to its difficult cultivation and destructive nature.

5 Conclusions

Drug repurposing is a well-known safety approach for the discovery of novel therapeutic agents, with the main benefit of shortening the time and expense of new drugs in a subsequent biomedical study. Tuberculosis caused by Mycobacterium tuberculosis has become a serious public health issue, and discovering a suitable therapy remains a major hurdle. According to the findings of our investigation, an integrated strategy combining pharmacokinetic profiling and molecular modeling techniques, such as molecular docking and molecular dynamics simulations, may be utilized to identify novel therapeutic compounds that precisely target the tuberculosis dihydrofolate reductase enzyme. The ADMET analyses reflected drug-likeness characteristics, whereas the molecular docking data represented binding affinity and hydrogen bond interactions of lead compounds CHEMBL577, CHEMBL161702, and CHEMBL1770248. The molecular dynamics simulations of these lead compounds revealed that the structures were consistent and equilibrated throughout the simulation in the binding site of the dihydrofolate reductase enzyme, and also satisfying results for the binding affinity estimation using the MM-GBSA technique. It should be noted that CHEMBL161702 and CHEMBL1770248 had the highest binding affinity in combination with the enzyme. To conclude, all the repurposed molecules presented might give a comprehensive idea for structure-based drug designing for tuberculosis treatment.

Author contributions

Conceptualization: Sajal Kumar Halder, Arafin Sultana, Ive Sultana. Formal analysis: Sajal Kumar Halder, Arafin Sultana, Md. Oliullah Rafi, Ive Sultana, Fatiha Elma. Investigation: Sajal Kumar Halder, Arafin Sultana, Md. Oliullah Rafi, Ive Sultana, Aparna Shil, Mahbubul Kabir Himel. Methodology: Sajal Kumar Halder, Arafin Sultana, Iqrar Ahmad, Ive Sultana, Md. Oliullah Rafi, Fatiha Elma, Israt Jahan Ananna, Iqrar Ahmad, Harun Patel. Project administration: Mahbubul Kabir Himel and Aparna Shil. Resources: Mahbubul Kabir Himel. Software: Sajal Kumar Halder, Arafin Sultana, Iqrar Ahmad, Ive Sultana, Md. Oliullah Rafi, Fatiha Elma. Supervision: Aparna Shil and Mahbubul Kabir Himel. Writing – original draft: Sajal Kumar Halder, Arafin Sultana, Ive Sultana, Md. Oliullah Rafi, Iqrar Ahmad, Harun Patel. Writing – review & editing: Md. Oliullah Rafi, Aparna Shil, Mahbubul Kabir Himel, Iqrar Ahmad, Harun Patel.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

The datasets supporting the conclusions of this study are included within the article (and its additional files). Supplementary information is available: chemical space visualization (Fig. 1), total energy differences of protein–ligand complexes (Fig. 2), and 100 ns poses of CHEMBL577, CHEMBL161702, and CHEMBL1770248 (Fig. 3–5). Supplementary tables provide binding energies of 1,026 ligands with Mycobacterium tuberculosis DHFR (Table 1), drug-likeness properties (Table 2), admetSAR results of 248 hits (Table 3), and ADMET profiles of the top 30 hits (Table 4). See DOI: https://doi.org/10.1039/d5ra04418a.

References

1. K. Floyd, P. Glaziou, A. Zumla and M. Raviglione, Lancet Respir. Med., 2018, 6, 299–314.
2. S. K. Halder and F. Elma, J. Clin. Tuberc. Other Mycobact. Dis., 2021, 24, 100246.
3. M. Daffé and P. Draper, Adv. Microb. Physiol., 1998, 39, 131–203.
4. B. Gao and R. S. Gupta, Microbiol. Mol. Biol. Rev., 2012, 76, 66–112.
5. R. Singh, S. P. Dwivedi, U. S. Gaharwar, R. Meena, P. Rajamani and T. Prasad, J. Appl. Microbiol., 2020, 128, 1547–1567.
6. L. F. Kuyper, J. M. Garvey, D. P. Baccanari, J. N. Champness, D. K. Stammers and C. R. Beddell, Bioorg. Med. Chem., 1996, 4, 593–602.
7. L. F. Kuyper, D. P. Baccanari, M. L. Jones, R. N. Hunter, R. L. Tansik, S. S. Joyner, C. M. Boytos, S. K. Rudolph, V. Knick, H. R. Wilson, J. M. Caddell, H. S. Friedman, J. C. Comley and J. N. Stables, J. Med. Chem., 1996, 39, 892–903.
8. R. L. Blakley, Adv. Enzymol. Relat. Areas Mol. Biol., 1995, 70, 23–102.
9. P. G. Young, C. A. Smith, P. Metcalf and E. N. Baker, Acta Crystallogr., Sect. D:Biol. Crystallogr., 2008, D64, 745–753.
10. J. Swarbrick, P. Iliades, J. S. Simpson and I. Macreadie, Open Enzyme Inhib. J., 2008, 1, 12–33.
11. B. I. Schweitzer, A. P. Dicker and J. R. Bertino, FASEB J., 1990, 4, 2441–2452.
12. J. M. Scott, Proc. Nutr. Soc., 1999, 58, 441–448.
13. A. D. Hanson and J. F. Gregory 3rd, Annu. Rev. Plant Biol., 2011, 62, 105–125.
14. M. P. Barrett, J. C. Mottram and G. H. Coombs, Trends Microbiol., 1999, 7, 82–88.
15. S. SHAKYA, K. Kasturi and K. Rao, PHARMANEST An, Int. J. Adv. Pharm. Pharmaceut. Sci., 2010, 1, 6–10.
16. C. C. Wang, Parasitology, 1997, 114(Suppl), S31–S44.
17. A. C. Anderson, Drug Discovery Today, 2005, 10, 121–128.
18. A. Rosowsky, J. B. Hynes and S. F. Queener, Antimicrob. Agents Chemother., 1995, 39, 79–86.
19. F. Derouin and C. Chastang, Antimicrob. Agents Chemother., 1989, 33, 1753–1759.
20. S. Hawser, S. Lociuro and K. Islam, Biochem. Pharmacol., 2006, 71, 941–948.
21. A. J. Salter, Rev. Infect. Dis., 1982, 4, 196–236.
22. I. B. Müller and J. E. Hyde, Mol. Biochem. Parasitol., 2013, 188, 63–77.
23. M. Sharma and P. M. S. Chauhan, Future Med. Chem., 2012, 4, 1335–1365.
24. J. A. Kovacs, C. J. Allegra, J. C. Swan, J. C. Drake, J. E. Parrillo, B. A. Chabner and H. Masur, Antimicrob. Agents Chemother., 1988, 32, 430–433.
25. F. R. Sattler, C. J. Allegra, T. D. Verdegem, B. Akil, C. U. Tuazon, C. Hughlett, D. Ogata-Arakaki, J. Feinberg, J. Shelhamer and H. C. Lane, J. Infect. Dis., 1990, 161, 91–96.
26. N. Gonen and Y. G. Assaraf, Drug Resistance Updates, 2012, 15, 183–210.
27. A. L. Jackman and A. H. Calvert, Ann. Oncol., 1995, 6, 871–881.
28. P. Cipriani, P. Ruscitti, F. Carubbi, V. Liakouli and R. Giacomelli, Expert Rev. Clin. Immunol., 2014, 10, 1519–1530.
29. C. Robson, M. A. Meek, J. D. Grunwaldt, P. A. Lambert, S. F. Queener, D. Schmidt and R. J. Griffin, J. Med. Chem., 1997, 40, 3040–3048.
30. E. L. White, L. J. Ross, A. Cunningham and V. Escuyer, FEMS Microbiol. Lett., 2004, 232, 101–105.
31. T.-S. Huang, C. M. Kunin, B.-S. Yan, Y.-S. Chen, S. S.-J. Lee and W. J. Syu, J. Antimicrob. Chemother., 2012, 67, 633–637.
32. A. Sardarian, K. T. Douglas, M. Read, P. F. Sims, J. E. Hyde, P. Chitnumsub, R. Sirawaraporn and W. Sirawaraporn, Org. Biomol. Chem., 2003, 1, 960–964.
33. C. Vilchèze and W. R. J. Jacobs, Antimicrob. Agents Chemother., 2012, 56, 5142–5148.
34. A. B. Gerum, J. E. Ulmer, D. P. Jacobus, N. P. Jensen, D. R. Sherman and C. H. Sibley, Antimicrob. Agents Chemother., 2002, 46, 3362–3369.
35. M. Hekmat-Nejad and P. K. Rathod, Exp. Parasitol., 1997, 87, 222–228.
36. V. Cody, N. Galitsky, J. R. Luft, W. Pangborn, A. Rosowsky and R. L. Blakley, Biochemistry, 1997, 36, 13897–13903.
37. R. Li, R. Sirawaraporn, P. Chitnumsub, W. Sirawaraporn, J. Wooden, F. Athappilly, S. Turley and W. G. Hol, J. Mol. Biol., 2000, 295, 307–323.
38. D. A. Matthews, R. A. Alden, J. T. Bolin, S. T. Freer, R. Hamlin, N. Xuong, J. Kraut, M. Poe, M. Williams and K. Hoogsteen, Science, 1977, 197, 452–455.
39. M. Kumar, R. Vijayakrishnan and G. Subba Rao, Mol. Diversity, 2010, 14, 595–604.
40. G. P. Miller and S. J. Benkovic, Chem. Biol., 1998, 5, R105–R113.
41. H. M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T. N. Bhat, H. Weissig, I. N. Shindyalov and P. E. Bourne, Nucleic Acids Res., 2000, 28, 235–242.
42. A. Gaulton, A. Hersey, M. Nowotka, A. P. Bento, J. Chambers, D. Mendez, P. Mutowo, F. Atkinson, L. J. Bellis, E. Cibrián-Uhalte, M. Davies, N. Dedman, A. Karlsson, M. P. Magariños, J. P. Overington, G. Papadatos, I. Smit and A. R. Leach, Nucleic Acids Res., 2017, 45, D945–D954.
43. S. Dallakyan and A. J. Olson, Methods Mol. Biol., 2015, 1263, 243–250.
44. L. L. C. Schrödinger, The PyMOL Molecular Graphics System, version 1.8, New York, 2015..
45. N. Guex and M. C. Peitsch, Electrophoresis, 1997, 18, 2714–2723.
46. S. K. Halder, M. M. Mim, M. M. H. Alif, J. F. Shathi, N. Alam, A. Shil and M. K. Himel, RSC Adv., 2022, 12, 24319–24338.
47. M. A. Rauf, S. Zubair and A. Azhar, Int. J. Basic Appl. Sci., 2015, 4, 168.
48. G. M. Morris, R. Huey, W. Lindstrom, M. F. Sanner, R. K. Belew, D. S. Goodsell and A. J. Olson, J. Comput. Chem., 2009, 30, 2785–2791.
49. A. Daina, O. Michielin and V. Zoete, Sci. Rep., 2017, 7, 42717.
50. C. A. Lipinski, F. Lombardo, B. W. Dominy and P. J. Feeney, Adv. Drug Delivery Rev., 2001, 46, 3–26.
51. A. K. Ghose, V. N. Viswanadhan and J. J. Wendoloski, J. Comb. Chem., 1999, 1, 55–68.
52. D. F. Veber, S. R. Johnson, H.-Y. Cheng, B. R. Smith, K. W. Ward and K. D. Kopple, J. Med. Chem., 2002, 45, 2615–2623.
53. I. Muegge, S. L. Heald and D. Brittelli, J. Med. Chem., 2001, 44, 1841–1846.
54. H. Yang, C. Lou, L. Sun, J. Li, Y. Cai, Z. Wang, W. Li, G. Liu and Y. Tang, Bioinformatics, 2019, 35, 1067–1069.
55. M. M. A. K. Shawan, S. K. Halder and M. A. Hasan, Bull. Natl. Res. Cent., 2021, 45, 27.
56. Schrödinger, LLC, Maestro, Release S.1, New York, NY, USA, 2018.
57. Schrödinger, LLC, Impact, Release S.1, New York, NY, USA, 2016 (2014–2018)..
58. R. Dash, S. M. Z. Hosen, M. R. Karim, M. S. H. Kabir, M. M. Hossain, M. Junaid, A. Islam, A. Paul and M. A. Khan, J. Appl. Pharm. Sci., 2015, 5, 73–78.
59. L. Xu, W. Jiang, H. Jia, L. Zheng, J. Xing, A. Liu and G. Du, Front. Cell. Infect. Microbiol., 2020, 10, 16.
60. Schrödinger, LLC,Desmond, version 5.9, New York, NY, USA, 2019..
61. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
62. I. Ahmad, D. Kumar and H. Patel, J. Biomol. Struct. Dyn., 2021, 1–13.
63. G. J. Martyna, Phys. Rev. E:Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1994, 50, 3234–3236.
64. I. Ahmad, H. Jadhav, Y. Shinde, V. Jagtap, R. Girase and H. Patel, In Silico Pharmacol., 2021, 9, 23.
65. R. V. Parikh, D. V. Pandya, P. Salaria, R. M. Amarendar, V. K. Vyas, H. G. Bhatt and T. M. Dhameliya, ChemistrySelect, 2024, 9, e202304151.
66. V. N. Savaliya, P. Rastogi, P. D. Savaliya, P. Salaria, R. M. Amarendar, V. K. Vyas, G. Natesan and T. M. Dhameliya, Chem. Select, 2025, 10, e202405449.
67. T. Koopmans, Physica, 1934, 1, 104–113.
68. Z. K. Genc, S. Tekin, S. Sandal, M. Sekerci and M. Genc, Res. Chem. Intermed., 2015, 41, 4477–4488.
69. A. R. Shama and M. L. Savaliya, Bioorg. Chem., 2025, 155, 108150.
70. S. K. Halder, I. Sultana, M. N. Shuvo, A. Shil, M. K. Himel, M. A. Hasan and M. M. A. K. Shawan, BioMed Res. Int., 2023, 2023, 5469258.
71. S. K. Halder, A. R. Sharma, T. A. Arian, S. Saha, A. Shil, M. O. Rafi, S. R. Sarker, M. N. Alam, M. K. Himel, M. A. Hasan and M. M. A. K. Shawan, Coronaviruses, 2024, 5, 85–113.
72. M. S. Ganesan, K. K. Raja, S. Murugesan, B. K. Kumar, G. Rajagopal and S. Thirunavukkarasu, J. Mol. Struct., 2020, 1217, 128360.
73. S. Tian, J. Wang, Y. Li, D. Li, L. Xu and T. Hou, Adv. Drug Delivery Rev., 2015, 86, 2–10.
74. M. P. Pollastri, Curr. Protoc. Pharmacol., 2010, 49, 1–8.
75. X. Chen, H. Li, L. Tian, Q. Li, J. Luo, Y. Zhang and J. Comput, Biol., 2020, 27, 1397–1406.
76. E. V Radchenko, A. S. Dyabina, V. A. Palyulin and N. S. Zefirov, Russ. Chem. Bull., 2016, 65, 576–580.
77. F. Broccatelli, E. Carosati, A. Neri, M. Frosini, L. Goracci, T. I. Oprea and G. Cruciani, J. Med. Chem., 2011, 54, 1740–1751.
78. M. S. Alavijeh, M. Chishty, M. Z. Qaiser and A. M. Palmer, NeuroRx, 2005, 2, 554–571.
79. S. A. Wrighton and J. C. Stevens, Crit. Rev. Toxicol., 1992, 22, 1–21.
80. F. P. Guengerich, Annu. Rev. Pharmacol. Toxicol., 1999, 39, 1–17.
81. P. Glue and R. P. Clement, Cell. Mol. Neurobiol., 1999, 19, 309–323.
82. A. Cheng and S. L. Dixon, J. Comput. Aided Mol. Des., 2003, 17, 811–823.
83. K. Mortelmans and E. Zeiger, Mutat. Res., 2000, 455, 29–60.
84. X. Li, L. Chen, F. Cheng, Z. Wu, H. Bian, C. Xu, W. Li, G. Liu, X. Shen and Y. Tang, J. Chem. Inf. Model., 2014, 54, 1061–1069.
85. X. Li, Z. Du, J. Wang, Z. Wu, W. Li, G. Liu, X. Shen and Y. Tang, Mol. Inf., 2015, 34, 228–235.
86. A. John, M. Sivashanmugam, V. Umashankar and S. K. Natarajan, J. Biomol. Struct. Dyn., 2017, 35, 2155–2168.
87. U. Vetrivel, S. Muralikumar, B. Mahalakshmi, K. Lily Therese, H. N. Madhavan, M. Alameen and I. Thirumudi, Genomics Inform., 2016, 14, 53–61.

---