Computational approaches: atom-based 3D-QSAR, molecular docking, ADME-Tox, MD simulation and DFT to find novel multi-targeted anti-tubercular agents

Abstract

Tuberculosis (TB) has become the biggest threat to human society because of the rapid rise in resistance to the causative bacteria Mycobacterium tuberculosis (MTB) against the available anti-tubercular drugs. There is an urgent need to design new multi-targeted anti-tubercular agents to overcome the resistance species of MTB through computational design tools. With this aim in mind, we performed a combination of atom-based three-dimensional quantitative structure–activity relationship (3D-QSAR), six-point pharmacophore (AHHRRR), and molecular docking analysis on a series of fifty-eight anti-tubercular agents. The created QSAR model had a R² value of 0.9521, a Q² value of 0.8589, and a Pearson r-factor of 0.8988, all of which are statistically significant. This means that the model was effective at making predictions. We performed the molecular docking study for the data set of compounds with the two important anti-tubercular target proteins, Enoyl acyl carrier protein reductase (InhA) (PDBID: 2NSD) and Decaprenyl phosphoryl-β-D-Ribose 20-epimerase (DprE1) (PDBID: 4FDO). We used the similarity search principle to do virtual screening on 237 compounds from the PubChem database in order to find strong anti-tubercular agents that act against multiple targets. The screened compound, MK3, showed the highest docking score of −9.2 and −8.3 kJ/mol towards both the target proteins InhA and DprE1, which were picked for a 100 ns molecular-dynamic simulation study using GROMACS. The data showed that the compound MK3 was thermodynamically stable and effectively bound to both target proteins in their active binding pockets without much movement. The analysis of the highest occupied molecular orbital (HOMO), lowest unoccupied molecular orbital (LUMO), and energy gap predicts the molecular reactivity and stability of the identified molecule. Based on the result of the above studies, the proposed compound MK3 can be successfully used for the development of a novel multi-targeted anti-tubercular agent with high binding affinity and favourable ADME-T properties.

Introduction

Tuberculosis (TB) is one of the oldest, contagious, fatal, and pervasive respiratory infections caused by the gram-positive bacteria Mycobacterium tuberculosis (MTB) [1, 2]. In recent years, during the COVID-19 pandemic, TB has re-emerged as a major world health problem that causes severe impairment in patients who require long-term treatment [3]. The World Health Organization's (WHO) report on TB indicates a significant increase in infection and death rates during this pandemic, primarily due to an increase in the frequency of multiple drug-resistant TB (MDR-TB) and extremely drug-resistant TB (XDR-TB) cases. This increase is attributed to patients' non-adherence and non-compliance with the available drug regimen [4, 5]. However, the emergence and spread of resistance to the currently available chemotherapeutic agents pose a growing risk to the global population, given the increasingly favourable conditions for the bacteria. These conditions include the HIV epidemic, other co-morbidities, such as type 2 diabetes, and low-quality living conditions in underdeveloped and economically backward countries. This underscores the urgent need for the development of drugs with shorter treatment times, simpler regimens, increased potency, and multi-targeted anti-tubercular agents that can combat the drug-resistant forms of this disease [6,7,8]. To achieve this goal, we used a computer-based drug design approach that aims to identify potential drug candidates and targets against MTB drug-resistant strains [9]. In this study, we used different computational approaches such as atom-based three-dimensional quantitative structure active relationship (3D-QSAR), pharmacophore modelling, molecular docking, pharmacokinetic, dynamic, toxicity study, and molecular dynamic simulation study to identify potential multi-targeted drug candidates used to treat drug-resistant tuberculosis [10,11,12].

In recent decades, many promiscuous drug targets for anti-tubercular action were reported, but two targets, Enoyl acyl carrier protein reductase (InhA) and Decaprenyl phosphoryl-β-D-Ribose 20-epimerase (DprE1), are considered the most clinically reproducible, effective, and highly vulnerable targets for treatment against MTB, MDR-TB, and XDR-TB [13,14,15]. The present work employs a computational approach to identify new and effective antagonists for these two crucial druggable targets in TB treatment.

The NADH-dependent enoyl-ACP reductase (InhA) enzyme is clinically validated as the target of the frontline anti-TB drug isoniazid (INH) and second-line drug ethionamide (ETA), encoded by the gene InhA of MTB [14]. The enzyme InhA catalyses the biosynthesis of mycolic acid, which is the central constituent of mycobacterial cell walls (Fig. 1). Mycolic acid production occurs through the fatty acid synthase (FAS) pathway. Two enzyme systems, fatty acid synthase I (FAS I) and fatty acid synthase II (FAS II), comprise it [16, 17]. FAS I produces short-chain fatty acids, while the FAS II pathway elongates these chains [18]. The X-ray structure of InhA reveals that each subunit has several α-helices and β-strands that contain NADH binding sites. For the last part of FAS II, the InhA enzyme reduce the double bond in the fatty acyl-ACP (acyl carrier protein) to make the saturated fatty acyl-ACP. This helps with the last part of the fatty acid elongation process [19, 20]. Therefore, compounds that can directly inhibit InhA without any activation disrupt the biosynthesis of mycolic acid in the mycobacterium and ultimately lead to the organism's death [21]. Hence, InhA inhibitors have a very promising opportunity towards the treatment of MTB, MDR-TB, and XDR-TB [22].

Fig. 1

Biosynthetic role of InhA and DprE1 enzyme for cell wall synthesis of Mycobacterium tuberculosis

Researchers have reported Decaprenyl phosphoryl-D-ribose 20-epimerase (DprE1) as a potential drug target for the treatment of tuberculosis (TB). The heteromeric protein DprE1 is an essential component for growth and survival of mycobacterium (Fig. 1). The DprE1 enzyme-mediated redox reaction synthesises polysaccharide arabinogalactan, the composition of the mycobacterial cell wall. During the reaction, the oxidase enzyme DprE1 carried out the conversion of decaprenylphosphoryl-d-ribose (DPR) to decaprenylphosphoryl-d-arabinose (DPA) by epimerization via an intermediate decaprenylphosphoryl-2-keto-β-derythro-pentofuranose (DPX) [23,24,25]. Inhibition of DprE1 disrupts the synthesis of arabinogalactan, weakening the bacterial cell wall and making the bacteria more susceptible to the chemotherapeutic agents used for the treatment of MTB, MDR-TB, and XDR-TB [26].

The process of developing new molecules using the virtual screening workflow has a crucial significance due to the addition of artificial intelligence (AI) and machine learning (ML) [27]. Identifying hit molecules through computational drug discovery has proved to be a meaningful methodology in recent years [28]. Structure-based similarity search and screening, one of the various approaches to drug design and discovery, has become a routine concept in the design and discovery of new chemotherapy molecules [29]. Similarity search is based on the concept that the two molecules having structural similarity share similar properties and biological action [30]. Thus, finding molecules that are similar to a known active molecule is one of the keys to drug discovery. Drug discovery based on similarity searches improves the odds of researchers finding more active molecules at the lowest cost and with the highest probability of success [30, 31]. Today's involvement of various in silico modules of computer-aided drug design (CADD), such as 3D-QSAR, molecular docking, ADME-T prediction, and simulation studies, has become increasingly important and significantly aids in identifying the most effective drug compounds for specific disease targets [32]. In this regard, we have conducted a study on all 58 2-nitroimidazooxazines derivative anti-tubercular agents, utilising CADD techniques to detect and identify highly effective multi-targeted drug candidates. These candidates are expected to form more stable chemical bonds with the two most potential protein targets, InhA and DprE1, of mycobacterium for the treatment of tuberculosis.

During the initial stage of our research, we utilized atom-based 3D-QSAR and ligand-based pharmacophore hypotheses to pinpoint the characteristics that contribute to the biological activity of the compounds in the data set, specifically those with anti-tubercular properties. We then conducted a molecular docking study on the ligands to determine their intermolecular interaction with the amino acid residues at the active site of the two target proteins, InhA and DprE1.

In the second phase, we conducted virtual screening of the PubChem database, using the best docked compound from the series as a reference compound to identify structurally similar compounds. We then screened the selected compounds using their docking results with the two target proteins, InhA and DprE1. Using the docking results, we finally subjected the screened compounds to ADME-Tox and drug likeliness studies, applying the Lipinski rule of five. The work has concluded with a molecular dynamic simulation study and density functional theory analysis to investigate the stability and reactivity of the identified ligand within the protein–ligand complex against InhA and DprE1 proteins.

Material and methods

Data set of ligands

A set of 58 2-nitroimidazooxazines derivatives was taken from the previously published literature for the present study, which are sharing the same activity and assay procedure with significant variations in their structure and potency [33]. The observed potencies of the compounds in the data set have IC⁵⁰ values ranging from 0.035 to 2.8 µm, which were further converted to pIC50 by using the mathematical formula given as Eq. 1:

$${\text{pIC}}^{{{5}0}} = \, - {\text{log 1}}0^{{({\text{IC5}}0)}}$$

(1)

To generate the 3D-QSAR models, the dataset of 58 compounds was randomly divided into a training set of 41 compounds and a test set of 17 compounds, as presented in Table S1. We generated the models using the training set of compounds and validated the developed models using the test set of compounds [34].

Preparation of ligands and alignment

Molecules selected for the study were constructed using the Schrodinger suite's Chem Sketch and then subjected to geometrical optimisation using the Ligprep module. After energy minimisation, low-energy 3D structures were obtained for each ligand. We aligned the ligands using the flexible ligand alignment option of Maestro software [34]. This is a crucial step in producing precise and accurate 3D-QSAR models [35]. We aligned all the data set ligands so they superimpose on each other, facilitating the study and observation of structural entity variations and their relationships (Fig. 2).

Fig. 2

Alignment structure of ligands

Pharmacophore modelling

Pharmacophore hypothesis modelling is commonly the spatial arrangement of different chemical features similar to two or more active ligands, which explains the interaction involved in binding ligands with the target protein [36]. We divided the ligands of the series into active and inactive according to their activity threshold value to generate a common pharmacophore hypothesis [37]. We kept the activity threshold values for active and inactive ligands at 6 and 5, respectively. We generated the pharmacophore model using a dataset of ligands with pIC50 distribution ranges from 5.553 to 7.523. The PHASE module of Schrodinger Maestro software was used to generate a pharmacophore model that provides a standard set of six pharmacophoric features like hydrogen bond acceptor (A), hydrogen bond donor (D), hydrophobic group (H), aromatic ring (R), negatively ionisable (N), and positive ionisable (P) groups that affect the ligand-target interaction [38, 39]. PHASE generates models based on the active ligands superimposed on features associated with the hypothesis.

A six-point common pharmacophore hypothesis was identified from all the active ligands having identical sets of features with very similar spatial arrangements and keeping a minimum intensities distance of 2.0 Å. The best common pharmacophore hypothesis was selected depending on the survival score. The high scoring hypothesis was used to create QSAR models.

Building of QSAR models

PHASE modules of software have two types of molecular alignment: the first is pharmacophore-based alignment, and the second is atom-based alignment [36]. The pharmacophore-based model falls short in elucidating the ligand's features and comprehensive molecular structure analysis, essential for the ligand's stearic interaction with target proteins [34]. The atom-based QSAR models study the entire molecular structure of the ligands, making them more useful in explaining structure–activity relationships. During the generation of atom-based 3D-QSAR models, the structural features of each atom are treated as van der Waals spheres [35, 40]. The atoms are treated as hydrogen bond donors-D (hydrogen bonded to elements like N, O, P, and S), hydrophobic or nonpolar-H (C, Cl, Br, F, and I), negative ionic groups-N (atoms of negative charge), positive ionic groups-P (atoms of positive charge), electron withdrawing, including hydrogen bond acceptors-W (non-ionic atoms like N, O), and miscellaneous-X (other types of atoms) as per simple internal rules [41, 42]. The study maps the ligand's features to a 3D cubic grid space. We achieved the generation of QSAR models by setting all the parameters to default and the PLS factor to 8. We generate atom-based 3D-QSAR models by assigning 70% and 30% of ligands to the training and test sets, respectively. We developed the models by considering descriptors as independent variables and biological activity as dependent variables.

Validation of the developed models

The developed QSAR models were used to predict the biological activities of new compounds; hence, to check the robustness of the generated atom-based 3D-QSAR models, both internal and external validation was performed [43]. The data set was divided randomly into training and test sets containing 41 and 17 compounds, respectively. Atom-based 3D-QSAR models were generated for the training set of compounds, and external validation was performed for the test set of compounds to check its predictiveness. The developed models were validated by considering statistical parameters like squared correlation coefficient (R²), cross-validated correlation coefficient (Q²) for the test set, standard deviation of regression, variance ratio (F), Pearson’s correlation coefficient (Pearson-r), root mean square error (RMSE), and significance level of variance ration (P). The predictive ability of the QSAR models for both training and test sets was analyzed based on the regression coefficient value (R²) and crossed validation coefficient (Q²) value [34, 35, 44].

Molecular docking study

The molecular docking simulation study is a computational approach that helps to find ligands that can effectively fit geometrically and energetically into the binding pockets of the target proteins. It also aids in predicting the types of energy interactions between ligands and target proteins [45]. In the present study, molecular docking was performed by PyRx (Autodock Vina) tools version 0.8 programs [46,47,48]. The docking poses with the least interaction energy were analysed and visualised by using Discovery Studio visualizer.

Protein preparation

The whole data sets of compounds were docked into the active site of the two most druggable targets of anti-tubercular action, NADH-dependent enoyl-ACP reductase (InhA) and Decaprenyl phosphoryl-β-D-Ribose 20-epimerase (DprE1). The X-ray diffraction-based, 3D crystallography structures of InhA and DprE1 having PDB ID 2NSD and 4FDO with good resolutions of 1.9 and 2.4 Å were retrieved from the RCSB protein data bank (www.rcsb.org). Further optimisation of the protein structure was done by using Biovia Discovery Studio. The missing hydrogen atoms and residues were added. We removed all the water molecules not involved in binding and co-crystallising ligands and performed energy minimization. The final 3D structure of the target proteins was evaluated using Biopredicta modules. The obtained Ramachandran plot (Fig. S1) revealed more residues in the most favoured regions, suggesting that the proteins are suitable for molecular docking studies.

Protein–Ligand docking

The protein–ligand docking study of the chosen protein–ligand complex was performed by using the virtual screening software interface PyRx (Autodock Vina) tools version v0.8. During docking analyses, protein structures were kept rigid and ligands were kept flexible [49]. The exhaustiveness was set at 8. The program performed energy minimization with the Universal Force Field (UFF) after uploading the chosen target proteins and ligands. The software's Open Babel tool then saved both ligands and protein structures in the '.pdbqt' format. We generated a grid box around the active binding site. We adjusted the grid box's size and coordinates by tracking the box's boundary line. PyRx employs the Lamarckian genetic algorithm as its conformational search algorithm. The present work employed semi-flexible docking as the docking method. The software displayed the binding energy with different conformers after docking and saved it in the '.csv' format. Autodock Vina splits the docking results into individual conformers. Next, we analysed the docking output files using Discovery Studio Visualiser (47) to study the interactions between the ligands and the amino acid present at the active site of target proteins. We loaded each conformer and the protein into Discovery Studio Visualiser and observed the interactions. We selected the best conformer based on the docking score and better non-covalent bond interaction.

Virtual screening

Virtual screening is an in silico, cost-effective, and high-speed technique that functions as a computational analogue of high-throughput screening (HTS). It involves the computational screening of chemical compounds from large libraries such as ZINC, PubChem, ChEMBL, and ChEBI, among others, for bioactive molecules [50, 51]. Researchers greatly benefit from this approach as it reduces the number of candidate molecules they need to test to manageable levels, thereby avoiding costly experiments testing thousands of compounds. Different approaches for the virtual screening of compounds are, first, the parallel approach, in which both ligand-based and structure-based are run independently and the best candidate compounds selected separately from both are considered for biological evaluation [52]. Secondly, the hybrid method, which combines both ligand-based and structure-based techniques into a standalone method, involves two approaches: (a) interaction-based methods and (b) a combination of molecular similarity and docking techniques [53]. Third, the reverse sequential approach includes structure-based virtual screening followed by 2D similarity searching, using the best hit molecule as a reference molecule. In this approach, the first docking of ligands on the target protein was performed to identify the active compound, and then the libraries of ligands for 2D similarity search with the initial active compound were explored [54].

In the current study, we selected compound number 56 as the most active hit molecule for 2D similarity-based virtual screening, based on the docking results of ligands from the data set with both target proteins. We then used this compound as a reference to identify 2D similar ligands from the PubChem database, using a similarity percentage of 70%. Based on the similarity search, approximately 237 ligands were identified, which were then screened using a docking study, drug-likeness study, and ADME-Tox study. The docking procedure was validated by re-docking the co-crystal ligand against the respective drug target proteins.

Pharmacokinetic and drug likeness prediction

In addition to the optimum binding affinities of the lead molecules with the target protein, the potency of the hit molecules is another driving factor in the drug development process. To become therapeutically successful and effective, the identified hits must possess high biological actions with low toxicity [55]. Experimenting with ADMET (A: Absorption, D: Distribution, M: Metabolism, E: Excretion, T: Toxicity) properties on small molecules is expensive and time-consuming. Therefore, the computational evaluation of pharmacokinetic (PK) and toxicity profiles of small molecules has proven to be an effective and crucial element in evaluating small molecules as drug candidates during the initial stages of drug development [56, 57]. Nowadays, the study of ADME-Tox properties has become an essential field of drug discovery, which significantly reduces the clinical failure of lead compounds. We made ADME-Tox predictions for all the ligands selected from virtual screening using ADMET Lab 2.0, a user-friendly, freely available web server (https://admetmesh.scbdd.com) [24, 32]. The properties assessed during the study are partition coefficient, aqueous solubility, % of oral absorption, plasma protein binding, skin permeability, blood–brain barrier, plasma protein binding, metabolism, and elimination. Additionally, we studied various toxicity aspects, including the maximum tolerated human dosage, hepatotoxicity, skin reactivity, mutagenicity, and hERG inhibitor. For drug-likeness analysis, we also examined the number of rotatable bonds, molecular weight, number of hydrogen bond donors, number of hydrogen bond acceptors, and topological polar surface area. We further subjected the lead compounds to estimating their drug-like properties using the Lipinski rule of five [58].

Molecular dynamic simulation

The molecular dynamics (MD) simulation study is very important for finding new small compounds that could be used as biological drug targets because it shows how ligands and proteins interact at the atomic level [59]. MD simulations help the researchers study the conformational changes, binding events, and structural stability of both protein targets and ligands. Molecular dynamics studies bridge the gap between structural information and the dynamic behaviour of target proteins, which aids in the rational design of potential drug candidates by providing a deeper understanding of their binding mechanisms and interactions with the various target proteins [60, 61]. The simulation study and generation of trajectory files were performed by GROningen MAchine for Chemical Simulations (GROMAC) software [62]. The best docking conformations of selected ligands with both the target proteins of PDB ID 2NSD and 4FDO were selected for the MD simulation study. The CHARMM27 force field and simple point charge (SPC) water solvation models were selected for study. A cubic boundary box and the counter ion Na⁺ Cl⁻ of concentration 0.15 M were added to neutralise the system. The energy minimisation was performed by selecting the steepest descent algorithm as an EM integrator with 5000 steps. Simulation was conducted under the equilibration parameters NPT and NVT at 300 k, 1 bar pressure, and a thermostat relaxation time of 100 ps. Leap frog was selected as simulator, and 100 ns under simulation time were executed using the mdrun program in GROMAC. Trajectories files were generated for analysing various dynamic parameters such as root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (RoG), solvent accessible surface area (SASA), binding free energy estimate (MM-PBSA), and H-bonds [63, 64].

Density functional theory (DFT) analysis

We conducted a density functional theory (DFT) analysis to examine the electronic properties of the best-identified hit from virtual screening [65]. We used the Gaussian 09 software tool to perform geometry optimisation and total energy calculations, utilising the Becke-3-Lee–Yang–Parr (B3LYP) function with the standard 6–311 + + G (d,p) basis set. Visualisation of the structure and the analysis of the outputs were carried out with Gauss-View software [66]. Frontier Molecular Orbital (FMO) studies can predict compounds' chemical reactivity and identify their stability. We calculated the energies of the highest occupied molecular orbital (HOMO), the lowest unoccupied molecular orbital (LUMO), and the gap between them to study the chemical stability of the identified molecule [67]. We found out the new inhibitor's chemical potential (μ), chemical hardness (η), chemical softness (S), and global electrophilicity (ω). We derived the 'index' mathematically using Eq. (2), taking into account the frontier molecular orbitals LUMO and HOMO. The chemical hardness (η), chemical softness (S), and global electrophilicity (ω) were computed using the expressions (3), (4), and (5), respectively [68, 69].

$$\mu \, = \, \left( {{\text{E}}_{{{\text{HOMO}}}} + {\text{E}}_{{{\text{LUMO}}}} } \right)/{2}$$

(2)

$$\eta \, = {\text{ E}}_{{{\text{LUMO}}}} - {\text{E}}_{{{\text{HOMO}}}}$$

(3)

$${\text{S }} = { 1}/\eta$$

(4)

$$\omega \, = \, \mu {2}/{2}\eta$$

(5)

Result and discussion

Pharmacophore model design

During pharmacophore model generation, the data set of 58 compounds was divided into active and inactive sets of compounds. The PHASE module of the Schrodinger software was used to generate six features (A, D, H, R, N, and P) based on 3D pharmacophoric models. The developed models help to predict biological activity by prognosis of the features necessary for ligand binding to target proteins. Using twenty-two active compounds from the "pharmaset," we are creating models that share common pharmacophoric features with these active sets of compounds. The scoring and ranking of generated pharmacophoric models were done to identify the best hypothesis. The scoring algorithm considered the site point configuration, vector magnitude, selectivity, and activity with overall energies. Based on its scoring, Table 1 selects the six-point pharmacophore hypothesis as the finest. This hypothesis includes features such as one H-bond acceptor (A), two hydrophobic groups (HH), and three aromatic rings (RRR), denoted as AHHRRR.

Table 1 Score of multiple Pharmacophore hypothesis AHHRRR

The top pharmacophore model with good predictive power for both active and inactive ligands was associated with the six-point hypothesis AHHRRR. Further, the predictability of a well-pharmacophoric hypothesis was confirmed by considering its survival score and adjusted score. The hypothesis, AHHRRR_1, has the highest survival score of 5.222 and an adjusted score of 3.81, making it the best hypothesis for predicting the structural features required by both active and inactive ligands to bind with their target protein and perform therapeutic action. The image of distances and angles between the pharmacophoric sites (Fig. 3a) and hypothesis images for active ligand (compound no. 3) and inactive ligand (compound no. 18) are represented in Fig. 3b, c, respectively. The pharmacophoric features (A) map the hydrogen bond acceptor to the etheric 'O' atom between the imidazooxazine and methyl biphenyl ring, while the two hydrophobic groups (HH) map to the -CF3 group attached to the 4th position of the benzene ring and the oxazine ring of the fused imidazooxazine ring. Among three aromatic rings (RRR) features, the first was present on the imidazole ring of the bicyclo imidazooxazine ring, and the second and third were visible on biphenyl rings attached to the imidazooxazine ring. The hypothesis suggests that the identified pharmacophoric features are crucial for the effective binding of ligands with target proteins, thereby demonstrating anti-tubercular action.

Fig. 3

a Common Pharmacophoric hypothesis (AHHRRR_1). The Pharmacophoric feature A: H-bond acceptor; appear as light pink sphere with two arrows, H: hydrophobic group; appear as green spheres, R: aromatic rings appear as orange torus in the plane of the ring. b Pharmacophoric hypothesis (AHHRRR_1) for active ligand 3. c Pharmacophoric hypothesis (AHHRRR_1) for inactive ligand 18

Generation of atom based 3D-QSAR model

We generated models relying on the alignment of the ligands in 3-dimensional space using the atom-based 3D-QSAR study. We randomly divided the data set of 2-nitroimidazooxazines, which contained 58 compounds, into 41 training sets and 17 test set molecules. We developed the atom-based 3D-QSAR models using the PHASE modules of Schrodinger software. The PHASE algorithm has the advantage of producing a 3D contour map based on favourable and unfavourable regions. The present study developed atom-based QSAR models for training sets, considering partial least squares (PLS) factor 8, and further validated them using compounds from the test set.

Analysis of developed QSAR models

The predictivity of the developed atom-based 3D-QSAR models with eight PLS factors was validated internally and externally for both training and test set compounds. The statistical parameters, squared correlation coefficient (R²), cross-validated correlation coefficient (Q²), standard deviation of regression, variance ratio (F), Pearson’s correlation coefficient (Pearson-r), root mean square error (RMSE), and significance level of variance ratio (P), were used to evaluate the quality of the QSAR models. The summary of the statistical data of all the developed atom-based QSAR models is listed in Table 2.

Table 2 Summary of atom based 3D- QSAR results

The PLS factor 8 model has the smallest standard deviation (SD), which is 0.1424. The squared correlation coefficient (R2) for the training set is 0.9525, and the cross-validated correlation coefficient (Q2) for the test set compounds is 0.8589. This shows that the model is effective at predicting the test set of compounds. The built QSAR model exhibits good precision and is suitable for further analysis and study, as indicated by the higher values of F (80.2), Pearson r (0.8988), and other statistical parameters falling within the acceptance range. Figure 4a, b present the linear scattered plots of actual versus predicted pIC50 for training and test sets, demonstrating the predictive power of the generated QSAR model.

Fig. 4

a Linear scattered plots of actual Vs predicted pIC⁵⁰ for training set. b Linear scattered plots of actual Vs predicted pIC⁵⁰ for test set

We developed the QSAR model based on features of the atoms attached to the core ring system, including hydrophobicity or nonpolarity, positive and negative ionic interactions, the electron-withdrawing effect, and other interactions. Table 3 tabulates the contribution of atom type fraction to the developed atom-based 3D-QSAR models. The atom type fraction contribution result indicates that the presence of hydrophobic or non-polar substitutions and electron withdrawing groups significantly contributes to the anti-tubercular activity, while the presence of positive and negative ionic interaction groups has a minor role in the anti-TB activity.

Table 3 Atom type fraction contribution of atom based 3D- QSAR models

We visualised the developed atom-based 3D-QSAR models in PHASE and conducted a study to correlate activity with various atomic contributions, using coloured cubes for both training and test set compounds. The developed QSAR models allowed us to find different atomic contributions to anti-tubercular activity, such as the presence of hydrophobic or non-polar groups, electron-withdrawing groups, and positive and negative ionic groups. This method used atom types and their occupancy positions in a grid of cubes to predict properties and visualise the regions that are favourable and unfavourable for anti-tubercular activity. Figure 5a–d for a training set compound (18) and Fig. 6a–d for a test set compound (41) display the maps generated for different atomic contributions in atom-based 3D QSAR. In these contour pictorial presentations of hydrophobic or non-polar interaction, the magenta colour cube shown is unfavourable, and the green colour cube is orable. For negative and positive ionic interaction, the yellow cube contributes positively, while the red and purple cube contribute negatively. Lastly, the electron-withdrawing map indicates that the green colour cube is favourable for the bioactivity of the ligand, while the red colour cube is avorable. The contribution map made for the atom-based 3D-QSAR study shows the structural details that are needed for ligands to interact with their target proteins. These maps also help us figure out which atoms or groups are attached to the core ring system and need a certain physiochemical property to make ligands more effective at fighting tuberculosis.

Fig. 5

a–d Atom based 3D QSAR visualization map of various atomic contribution for a training set compound: a Hydrophobic or non-polar b Negative ionic c Positive ionic d Electron withdrawing

Fig. 6

a–d Atom based 3D QSAR visualization map of various atomic contribution for a test set compound: Hydrophobic or non-polar b Negative ionic c Positive ionic d Electron withdrawing

Molecular docking study

All the 58 compounds of the data set were docked into the binding pockets of the two most effective and potential drug target proteins for anti-tubercular action, NADH-dependent enoyl-ACP reductase (InhA) and Decaprenyl phosphoryl-β-D-Ribose 20-epimerase (DprE1) having PDB ID 2NSD and 4FDO, respectively, by PyRx Tools version v0.8 using Autodock Vina. Table S2 lists the drug-binding scores in kcal/mol for all of the compounds. Afterwards, we selected compound 56, which had the highest interaction energy of −8.2 kcal/mol and −9.6 kcal/mol with both the target proteins InhA and DprE1, for analysis. We then used this compound to retrieve compounds with structural similarity up to 70% from the PubChem database. The selected compounds from the database were further screened by using docking, ADME analysis, and molecular simulation (MD) studies.

Structural similarity based virtual screening

To find the effective multi-targeted anti-tubercular agent, compounds from the PubChem database were screened based on structure similarity by taking compound no. 56 as a reference compound. About 237 ligands were identified, and their structure was retrieved from the database for further screening. These compounds were then docked into the grid pockets of both target proteins having PDB ID 2NSD and 4FDO. Based on the significant docking scores of >−6.5 kcal/mol for both targets, only nine compounds were selected for screening of their drug-like properties and ADMET predictions. For validation of the docking study, the co-crystallised ligand of the receptor was extracted and re-docked into the binding pockets of the respective target proteins. The result of the docking study of the top-ranked compounds has been reported in Table 4.

Table 4 Docking result of top identified hits and co-crystalize ligands

After analysing the compounds that were screened, it was found that CHEMBL566642 (MK3) had the best docking score, scoring −9.2 and −8.3 kcal/mol in the binding pockets of both of the chosen druggable targets for anti-tubercular activity. Fig. S2 presents the 3D docking poses of reference compound (56), screened compound (MK3), and co-crystallise ligands for the protein targets. Table 5 reports the 2D docking interaction results of the reference compound (56), the screened compound (MK3), and the co-crystallise ligands for both receptors. The 2D docking poses of the reference compound (56), the screened compound (MK3), and the co-crystallised ligand with InhA and DprE1 proteins are given in Fig. 7a–f, respectively. Upon examining the docking features of the identified hit (MK3) with target proteins InhA (2NSD), it was found that it has formed one H-bond with Ile194, П- П stacked interaction with Phe149, Ile194, and П-alkyl interaction with Ala157 and Ile215 residues at the active site of protein (Fig. 7b) and formed two H-bonds with Asn135, Asn144, carbon-hydrogen bond with Thr225, Glu190, Gly140, П- П stacked interaction with His415, and П-alkyl interaction with Tyr226 and Ala139 residues present at the active site of target protein DprE1(4FDO) (Fig. 7e).

Table 5 Docking interactions result of the reference, screened and co-crystallize ligands with amino acid residues of target proteins

Fig. 7

a–f Docking interactions with InhA (PDB ID: 2NSD) a Docking interaction diagram of Reference ligand (56) b Docking interaction diagram of Identified hit (MK3) c Docking interaction diagram of Co-crystallize ligand; Docking interactions with DprE1 (PDB ID: 4FDO) d Docking interaction diagram of Reference ligand (56) e Docking interaction diagram of Identified hit (MK3) f Docking interaction diagram of Co-crystallize ligand. H-bond shown as bold green line, light green indicates carbon hydrogen bond, purple colour bond is П- σ interaction, dark pink bond is П- П stacked interaction, light pink indicates П-alkyl interaction with amino acid residues

In silico drug likeness and pharmacokinetic (ADME-T) analysis

After the virtual screening for structural similarity, the nine best hits were used to make drug-likeness and ADME-T predictions using the open web server ADME-T Lab 2.0. The Lipinski rule of five violations will be used to analyse the druglikeness property. We also analysed other ADME-T properties such as water solubility, pharmacokinetics, and ligand toxicity. Drug-likeness studies qualitatively measure the chance of a molecule turning into an oral drug in terms of its bioavailability. Table 6 summarises the drug-likeness and Rule of Five prediction properties for the top nine compounds.

Table 6 Predictions of drug-likeness and rule of five for the top screened compounds

The results demonstrated that all nine screened compounds exhibited good drug-likeness, with zero violations of rules. It was important to find the right values for the molecular weight (≤ 500), LogP (≤ 5), number of hydrogen bond acceptors (0–12), number of hydrogen bond donors (< 05), number of rotatable bonds (0–11), and topological polar surface area (< 140Å2) because they all affected how well the drug could be absorbed by the body when taken by mouth. All the compounds demonstrate an excellent synthetic accessibility score of less than 06, indicating the complexity of their molecular structure and ring system and their ease of synthesis. Advanced knowledge of pharmacokinetic and toxicity study results is helpful in designing potential drug candidates with less toxicity. All the screened compounds were evaluated for their drug-like behaviour through analysis of pharmacokinetic properties and toxicity studies. Table 7 lists the results. For all the identified compounds, the LogS value is between −4.5 and 0.5 log mol/l, indicating good aqueous solubility, which is important for estimating the absorption and distribution of drugs within the body (Fig. S3). For all the compounds that were tested, the predicted values for plasma protein binding (PPB) and blood–brain barrier (BBB) are between 0.0 and 0.3 for PPB and 80 to 90% for BBB. The expected values for drug-induced liver injury (DILI) and human hepatotoxicity (H-HT) are within the acceptable range. This means that the compounds were toxic to the liver at high doses. The Ames mutagenicity and skin sensitisation study's acceptance value indicates that these compounds are safe from carcinogenicity and inflammatory skin reactions. All of the results for the pharmacokinetic properties are well within the acceptable range for use in humans. This shows that they could be used as new multi-targeted anti-tubercular drugs.

Table 7 Pharmacokinetic (ADME) and Toxicity prediction results for the top screened compounds

Prediction of anti-tuberculosis sensitivity

Further, the screened compounds were investigated to predict their minimum inhibitory concentration (MIC) against eight different Mycobacterium species by employing an online mycoCSM server [70]. Only Mycobacterium tuberculosis (MTB) MIC values were extracted and analysed with marketed standards (isoniazid and rifampicin). Table 8 displays the predicted MIC values calculated by mycoCSM. The result indicates that the MIC value of the hit molecule MK3 (−6.181 µM) was close to the MIC value of rifampicin (−6.130 µM) and higher than that of isoniazid (−4.942 µM). The compound MK3 with a lower MIC value requires less to inhibit the growth phase of the organisms, indicating its potential as an anti-tubercular agent for further study [70].

Table 8 Anti-TB activity prediction of screened hits through online server mycoCSM

Molecular dynamic simulation study

The MD simulation study is a computational technique that informs alterations in protein structure and behavior that occur throughout the simulation period. MD simulation can also be helpful in studying protein dynamics, folding, stability, and interaction of proteins with ligands. The present study conducted MD simulations to confirm the stability of the InhA-MK3 and DprE1-MK3 complexes in physiological environments, a feat that molecular docking could not accomplish. Based on the molecular docking scores, we selected the best-screened compound (MK3) for MD simulation analysis, along with the reference compound (56) to co-crystallise ligands. We ran MD simulations using Gromacs software, setting the best dock poses of the hit molecule with target proteins. We estimated the stability of the MK3 binding complex with both target proteins InhA and DprE1 by evaluating the plots of RMSD, RMSF, RoG, SASA, H-bonds, and the binding free energy estimate (MM-PBSA).

RMSD (Root Mean Square Deviation)

RMSD found the average distance between atom positions in the simulated structure and the initial reference structure. This shows how far the molecular dynamics (MD) simulated structure has changed from its original shape. A system with a lower RMSD value has less structural drift and is therefore more stable. This picture (Fig. 8a) shows the RMSD plot of the reference compound (56), the identified compound (MK3), and the co-crystallise ligand with the InhA protein. It shows that the complex form is stable. The InhA-56 complex initially stabilised, then transitioned into instability between 30 and 50 ns before stabilising again from 70 ns to the end. In the InhA-MK3 complex, we observed steady confirmations at the beginning, followed by unsteady confirmations from 30 to 65 ns, and then a return to linearity until 100 ns. The InhA-cocrystallize ligand complex initially displayed fluctuations and maintained a consistent RMSD value between 30 and 100 ns. The RMSD plot for the DprE1 protein is shown in Fig. 8b. It shows the reference compound (56), the identified compound (MK3), and the co-crystallise ligand. This plot shows no significant deviation from the unbound protein. The plot for all three compounds shows stability throughout the entire 100 ns simulation time. The simulated analysis conducted against the proteins InhA and DprE1 reveals that the screened compound MK3 has RMSD values that do not surpass 0.35 nm. It follows that the complex formed with the target proteins InhA and DprE1 has no major modifications and is stable during a 100 ns simulation analysis, as indicated by the identified compound MK3 having an RMSD less than 0.5 nm [71]. As a result, it offers a good chance to continue developing and getting better as a strong anti-tubercular drug.

Fig. 8

a RMSD trajectory plot of InhA protein (Apo protein) with reference compound (56), identified hit (MK3) and co-crystallize ligand. b RMSD trajectory plot of DprE1 protein (Apo protein) with reference compound (56), identified hit (MK3) and co-crystallize ligand

RMSF (Root Mean Square Fluctuation)

During MD simulation, root-mean-square fluctuation (RMSF) was assessed to analyse the impact of lead compounds binding on the flexible portion of the targeted protein. The RMSF result also estimates each residue's variations around its average location. Higher RMSF values indicate that the residues are more flexible. RMSF can assist in identifying flexible and stiff protein regions. The found hit molecule MK3 interacts with both the InhA and DprE1 proteins and stays stable in their residues during the simulation study. The RMSF analysis in Fig. 9a, b shows that the hit compound was more stable than the reference compound (56) and the co-crystallise ligand complex with the proteins InhA and DprE1. However, several residues such as Arg45, Phe109, Arg153, Ile202, Gly205, Trp249, and Leu269 for the InhA complex and Thr8, Arg41, Phe267, Arg304, pro329, Phe362, Arg372, and Lys398 for the DprE1 complex are highly flexible and showed significant RMSF. The RMSF analysis shows that the hit molecule is much more stable than 56 and forms co-crystallised ligands for both target proteins. The presence of significant RMSF residues outside the active site of target proteins suggests that any conformational changes undergone may not have a substantial impact on the binding ability of MK3 into the active site of the target proteins of InhA and DprE1 [72].

Fig. 9

a RMSF plot of InhA protein with reference compound (56), identified hit (MK3) and co-crystallize ligand. b RMSF plot of DprE1 protein with reference compound (56), identified hit (MK3) and co-crystallize ligand

RoG (Radius of Gyration)

The Radius of Gyration (RoG) measures the dispersion of a protein's mass around its centre of mass, which helps to identify the protein structure's expansion and compactness. A reduction in RoG indicates a folded or compact structure of the protein, while an increase in RoG value indicates an unfolded structure. We analysed the complexes of both target proteins InhA and DprE1 with the reference compound (56), the identified compound (MK3), and the co-crystallise ligand for RoG, and presented the results as Fig. 10a, b, respectively. When MK3 complexed with InhA and DprE1, it showed lower RoG values of 1.82 and 2.16 nm compared to the complexes of InhA and DprE1 proteins with reference compound 56 and co-crystallise ligands, which had RoG values of 1.84, 1.86, and 2.2, 2.22, respectively. This means that the structures of InhA-MK3 and DprE1-MK3 are more compact than the structures of InhA and DprE1 when they mixed with 56 and co-crystallise ligands. This RoG study result shows that the Rg trajectory stayed stable throughout the whole simulation study for MK3. This means that it had a compact structure while interacting with both target proteins, which is a good sign for RoG stability [71].

Fig. 10

a RoG plot of InhA protein with reference compound (56), identified hit (MK3) and co-crystallize ligand. b RoG plot of DprE1 protein with reference compound (56), identified hit (MK3) and co-crystallize ligand

SASA (Solvent Accessible Surface Area)

The surface area of protein that is accessible to the solvent is referred to as its "solvent accessible surface area" (SASA). Changes in SASA may be a sign of protein-ligand interactions, folding, ligand binding, or conformational changes. SASA is frequently used to examine the dynamics and stability of proteins. For the SASA analysis, a complex was made up of the reference compound (56), the identified compound (MK3), and the co-crystallise ligand with the InhA and DprE1 proteins (Fig. 11a, b). According to the results, 135.561 and 180.221 nm² were predicted as the average SASA values of screened compound MK3 with InhA and DprE1 proteins, whereas reference compounds and co-crystallise ligands had 136.281, 135.531, and 192.215, 189.287 nm², respectively. The study showed a very slight deviation during the simulation due to minor structural changes during complex formation. The MK3 complex with InhA and DprE1 confirms that this compound has acceptable stability [71].

Fig. 11

a SASA plot of InhA protein with reference compound (56), identified hit (MK3) and co-crystallize ligand. b SASA plot of DprE1 protein with reference compound (56), identified hit (MK3) and co-crystallize ligand

The SASA measurements further supported the results from the RMSD, RMSF, and RoG studies by providing additional information on the stability of MK3 in interaction with InhA and DprE1 target proteins for anti-tubercular action.

H-bonds analysis (Hydrogen Bonds Analysis)

H-bond analysis during MD simulation illustrates the degree of interaction between ligand and protein and the stability of the protein–ligand complex throughout the simulation. Figure 12a, b present the H-bond analysis of the complex of the reference compound (56), the identified hit (MK3), and the co-crystallisation of ligand with target proteins InhA and DprE1. During the simulation study, the screened compound MK3 formed an average of 2–4 and 2–5 H-bonds within the target sites of InhA and DprE1, respectively, according to the H-bonding interaction plot. The graph clearly shows that both complexes conserved the H-bonds formed between the ligand and amino acids of the target protein during the 100 ns simulation. Finally, the H-bond analysis confirmed what other structural studies had found by showing that MK3 is stable when it comes into contact with both InhA and DprE1 proteins [71].

Fig. 12

a H-bond analysis plots for reference compound (56), identified hit (MK3) and co-crystallize ligand with InhA. b H-bond analysis plots for reference compound (56), identified hit (MK3) and co-crystallize ligand with DprE1

Binding free energy calculation (MM-PBSA)

Calculating the binding free energy allows for the assessment of the protein–ligand complex's energetic stability and the prediction of their binding strength. We use the gmx_mmpbsa tool of GROMACS to calculate the binding free energy of the ligand–protein complex. This gmx_mmpbsa applies the Molecular Mechanic-Poisson-Boltzmann Surface Area (MM-PBSA) method for binding energy calculation. We calculated the binding free energy for the docked complexes of target proteins InhA and DprE1 with compounds 56 and MK3, and co-crystallised ligands using expression 6.

$$\Delta {\text{G}}_{{{\text{bind}}}} = {\text{ G}}_{{{\text{complex}}}} - \, \left( {{\text{G}}_{{{\text{protein}}}} + {\text{ G}}_{{{\text{ligand}}}} } \right)$$

(6)

where G_complex is the energy of the protein–ligand complex, G_protein and G_ligand are the energy of protein and ligand in aqueous solvent, respectively. Other energies like van der Waals energy (E_vdw), electrostatic energy (E_elec), polar solvation energy (G_polar), non-polar solvation energy (G_nonpolar) were also calculated and shown in Table 9.

Table 9 MM-PBSA calculations for the complex of DprE1 and InhA protein with compound 56, MK3 and co-crystallize ligands

The binding free energy calculation indicates that Evdw, Eelec, and Gnon-polar have a significant contribution to the total binding energy because of their negative values, whereas Gpolar has no contribution because of its positive value. Compound MK3 is found to have a very good binding energy of −42.27 kJ/mol and −32.39 kJ/mol with the target proteins InhA and DprE1, respectively, which suggests that this compound has strong and effective interaction within the active binding pockets of these proteins.

The simulation study's results collectively indicate that the screened ligand MK3 possesses strong affinity and energy stability, indicating its potential as an InhA and DprE1 inhibitor. These promising findings underscore the potential of compound MK3 as a multi-targeted anti-tubercular agent.

DFT analysis

The Gaussian 09 software tool performed the DFT analysis for compound MK3, analysing its electronic properties and chemical reactivity using the Becke-3-Lee–Yang–Parr (B3LYP) function and the standard 6–311 + + G (d,p) basis set. Visualisation of the structure and the analysis of the outputs were carried out with Gauss-View software. FMO analysis of the identified compound was performed by measuring the energy of the highest occupied molecular orbital (HOMO), lowest unoccupied molecular orbital (LUMO), and the gap between them. Measurements of HOMO and LUMO indicate the compound's nucleophilicity and electrophilicity. The energy gap between HOMO and LUMO checks the reactivity of the compound; a compound with a small energy gap allows the molecule to be more reactive but less stable, while a compound with a higher difference margin is less active but more stable. Along with HOMO, a HOMO and LUMO energies were used to figure out the new inhibitor MK3's bioactivity. Global reactivity indices like chemical potential (μ), chemical hardness (η), chemical softness (S), and global electrophilicity (ω) were also calculated. Study are represented in Table 10.

Table 10 Global indices of the screened compound MK3

Furthermore, the Frontier Molecular Orbitals (FMOs) of the screened molecule are given in Fig. 13. The derived values of ΔEgap and global indices indicate significant chemical and bio-reactivity of the identified molecule MK3 [69].

Fig. 13

The geometries of the HOMO and LUMO orbitals, along with the value of ΔE_gap of the compound MK3

Conclusion

Tuberculosis is regarded as one of the fatal infections caused by Mycobacterium tuberculosis in humans, which triggers mobility and morbidity throughout the world because of the upbringing of drug resistance cases. To cut down on the time needed to treat MTB strains that are resistant and the cost of making new drug candidates, this study combined a number of computer-aided drug design (CADD) methods to create new, effective molecules that target two possible anti-tubercular drug targets, InhA and DprE1. These molecules were found through virtual screening using similar structure-based drug discovery methods. We used atom-based 3D-QSAR analysis, ADMET profiling, molecular docking, MD simulation, and DFT analysis to learn more about the title molecule and see if it could be used to fight drug-resistant MTB. We used predictive validated atom-based 3D-QSAR and pharmacophore hypothesis models, through a rigorous assessment process, to identify molecular descriptors that influence the enhancement of anti-tubercular activity. The proposed compounds, MK1-MK9, confirmed the Lipinski rule of five. We saw that compound MK3 had the highest binding energy of −9.2 and −8.3 kcal/mol during the molecular docking process. This meant that it was a lead-like drug. In addition, it had the best interaction with the residues in the binding pockets of the target receptors, InhA and DprE1. The ADME-T results showed good absorption, no penetration into the brain, and no toxicity for the newly identified molecule. Furthermore, the 100 ns MD simulation results support the molecular docking results and indicate that the compound, MK3, demonstrated stable interactions with effective RMSD, RMSF, RoG, SASA, and H-bond formation within the active pockets of InhA and DprE1 proteins. In MM-PBSA analysis, the compound MK3 has a binding energy of −42.27 kJ/mol and −32.39 kJ/mol towards InhA and DprE1 proteins, suggesting this molecule has the strongest and most effective binding within these proteins' active pockets. Again, the result of DFT analysis suggests that the molecule MK3 tends to exhibit more active anti-tubercular action because it has a smaller ΔEgap of 0.14806 between HOMO and LUMO. After validating all the theoretical results, we found that the identified molecule MK3 has a significant potential to function as a new inhibitor of the two most druggable targets, InhA and DprE1, for the treatment of tuberculosis. The outcome of this computational study suggests that the compound MK3 could potentially develop into a potent multi-targeted anti-tubercular agent.