Tuberculosis is the foundation of the 2^nd major reason of human fatalities just after HIV which leads to ‘AIDS’ (World Health Organization, 2023). The persistent spread of anti-microbial resistant strains has reduced the number of usable anti-tubercular drugs (Soltanet al., 2021). Genetic research has denoted that resistance is the outcome of spontaneous point mutations, deletions and insertions in genes that are involved in drug inactivation (Alcocket al., 2020). Drug resistance to antimicrobials are divided into 4 categories (Reygaert, 2018). There are 3 virulence factors like number of enzyme coding genes, superfamilies of enzymes and active drug effiux systems which are involved in antimicrobial resistant activities of M. tuberculosis strains (Coleet al., 1998). The cell wall decreases permeability of antimicrobial agents (Bentley, 1997). The virulence factors encoded by a number of genes of pathogens, produce disease in host cells (Ruiz-Bacaet al., 2021) and these genes are related to virulence associated life style (Echeverria-Valenciaet al., 2018). These virulence factors rely upon its transmission across cell membrane of bacteria (Comaset al., 2010). Moreover, various transporter systems are involved in drug resistance (Poseyet al., 2006). The ‘KatG’ which is the most important amid all protective enzymes, detoxifies the toxic elements and neutralizes redox or oxidative stress (Trivediet al., 2012). MTB strains having mutations at codon 315 in ‘KatG’ gene persist in virulent form (Pymet al., 2002). Moreover, ‘KatG’ protein is actively connected with pathogenicity (DeVito and Morris, 2003). The in silico approaches, are designed on the basis of wet-lab based results and literature reports denoting its successful application (Barhet al., 2011). Now a days, bioinformatic approaches, have achieved remarkable importance in the area of potential drug target identification (Boeckmannet al., 2003) and development of promising drug compounds as they have brought down the costs and time needed for wet-lab investigation (Boeckmannet al., 2003; Linet al., 2020).

The present study has, been oriented to exploit various computational approaches like newly designed DNN (Mukul and Sarkar, 2022) including CodonW and Genetic Algorithms (G.A) for biological data extraction of various MDR-strains of M. tuberculosis. The subtractive genomics approaches, structure-based approaches, molecular dynamics simulation including ADMET properties analysis have been undertaken.

The present study has been conducted to identify promising candidate for KatG inhibition in M. tb_DKC2 strain by employing advanced computational methods.

The whole genome sequence of M. tb_DKC2 were retrieved from NCBI GenBank and submitted to CodonW followed by GA with Random Forest for data extraction. The DNN trained with genomic data of MTB strains was executed on dataset to predict essential proteins corresponding to essential genes followed by subtractive genomics methods for drug target identification.

The katG protein was predicted as drug target by retrieval of essential proteins from NCBI GenBank and it was submitted to SWISS MODEL (Waterhouseet al., 2018) for homology model. PROCHECK, PSIPRED, PROSA and ERRAT Computational tools were used for validation of katG protein.

The active binding site with co-crystalized ligand (HEME) of template protein was identified by Biovia Discovery Studio Visualizer (BIOVIAiovia, 2017). A ‘Grid Box’ was constructed on the same active binding site using 86.4388, 37.3258 and 50.228 co-ordinates of center pointing with a size of 18.3818Å, 15.9234Å and 22.0152Å in X, Y and Z directions, respectively. The ‘Grid Box’ with same parameters was formed on KatG modelled protein through virtual screening using UCSF Chimera v 1.12 user interface to predict active binding site and the same was confirmed by using DoGSiteScorer (Volkameret al., 2012).

The naturally occurring ligand ‘HEME’ found within the active binding site of template protein (2CCD) was retrieved from RCSB PDB for validation of docking parameters of Autodock Vina by redocking of template protein and the validated Autodock Vina parameters was used during docking of modelled protein with ligands of ‘ZINC’ database.

The modelled protein prepared by computing Gasteiger charges and adding polar hydrogens, was used to perform virtual screening against compounds (ligands) obtained from ‘ZINC’ database using Autodock Vina with validated docking parameters. Eventually, best possible ZINC compounds were obtained.

The complexed ligands obtained from ZINC database, were extracted using UCSF Chimera v 1.12 (Buttet al., 2020) while 3D modelled protein defined as ‘Receptor’, was prepared by computing Gasteiger charges and adding polar hydrogens using UCSF Chimera v 1.12 (Buttet al., 2020).

Autodock Vina v 1.12 (Tanchuket al., 2016; Gaillard, 2018) was executed for molecular docking of KatG (PDB ID: 2CCD) with its ‘HEME’ ligand. An iterated Local Search global optimization algorithm based on BFGS method was executed in docking study and UCSF Chimera v 1.12 was used for various orientations and conformations of ‘HEME’ ligand binding in active site of template. The lowest binding (-16.8 kcal/mol) score was selected which indicated stable binding of ‘HEME’ ligand. The modelled protein was then docked with prepared ZINC ligands using same parameters which were used for AutoDock validation. Altogether 3 hit compounds were identified as drug like candidates and ranked as per their binding energies to verify the ligand binding site. The post docking analysis of 3 hit compounds was performed using BIOVIA Discovery Studio Visualizer v 4.5 (BIOVIAiovia, 2017).

The prediction of ADMET properties of 3 selected compounds (ZINC-65407173, ZINC-09849203 and ZINC-09349296) was performed using SWISS ADME (Dainaet al., 2017) and pkCSM (Pireset al., 2015) computational tools (Table S2 and S3) while OSIRIS property explorer (Sander, 2001; Oladejoet al., 2023) was executed to ascertain drug-likeness and other properties of the same drug like compounds (Table 2).

Both protein-HEME and protein-ZINC65407173 complexes were subjected to MD Simulation using GROMACS 2021.4 software (Hesset al., 2008). The pdb2gmx module of GROMACS and SWISS-PARAM (Zoeteet al., 2011) were used to built topology of protein and ligands, respectively. CHARMM v 27 (Bjelkmaret al., 2010) was assigned to protein and ligands for intramolecular interactions. The long range coulombic electrostatic force was treated by PME algorithm. The neutralization of charges was performed by adding Na⁺ and Cl^– ions. The solvation of both complexes was performed using TIP3P water cube model. The steepest descent method was used to minimize energy. The system was set up maintaining NVT and simulation was performed for 100ps at constant 300k using modified Berendsen Thermostat algorithm. NPT ensemble was performed to stabilize pressure at 1 atm for 100ps employing Parrinello-Rahman and Berendsen algorithms. Finally, simulation was performed at both constant pressure and temperature of 1 bar and 310K, respectively employing time step of 2fs. H-bonds were constrained with LINCS algorithm and 60ns simulation was generated by recording coordinate trajectories at every 10ps. The graphs were plotted using XMGrace software (Turner, 2005). The trajectory files were analysed by HeroMDAnalysis tool (Rawatet al., 2021) to calculate change in Gibb’s free energy of binding for protein-ligand complexes.

Figure 1:
A. Predicted binding of Catalase-peroxidase: Displayed by Pro100(3.270Å), Gly273(3.316 Å), Gly269(3.417Å), His270(2.452Å), Lys274(3.228Å), Trp321(3.384Å), Thr275(3.262Å), Leu265(3.215Å), His276(2.501Å), Val230(3.363Å), Thr314(3.259Å), Thr315(2.905Å), residues. B. The superimposed docked complex of pre-docked ‘HEME’ (Sky Blue) over post-docked (Khaki) depicting the validation of docking study in terms of RMSD value indicated as 0.525Å.

The most active binding site of KatG was selected using ‘HEME’ ligand in order to perceive functions of protein and the active binding site was further confirmed using ‘DoGSiteScorer’ (Figure 1). RMSD value>2.0 Å justified suitable method for re-creation of binding pose (Shoichetet al., 2002).

The highest simple score (0.63) predicted by DoGSiteScorer, confirmed the most active binding site of target protein amid 26 druggable binding sites for the purpose of docking analysis (Uddin and Rafi, 2017). Molecular docking analysis of 2632 druglike compounds was performed to ascertain potential lead compound against KatG (Figure 2). Three (3) best hit compounds viz. ZINC65407173, ZINC09849203 and ZINC09349296 with binding energy of -11.3, -10.3 and -10.9 Kcal/mol., respectively (Table 1) were identified by Autodock Vina 1.1.2.

Hit compounds retained substantial interactions like H-bonding, Van der Waals interactions, Pi-Pi stacking, Pi-Sigma interactions (Figure 3) which may have important role in inhibition of KatG (Oladejoet al., 2023). Highest confidence score of 0.997, was generated for ZINC65407173 (Table S1) The docking interactions was more in respect of ZINC65407173 RMSD value of 0.525Å as output of re-docking validated the virtual screening of hit compounds, ADME profiling (Table S2) and Toxicity properties (Table S3) of hit compounds were analysed.

The toxicity risks, drug likeness, topological polar surface area (Å), drug score and ClogP values of 3 drug like compounds (Table 2) were investigated using ‘OSIRIS property explorer tool’ (Sander, 2001). MD simulation of two docked complexes viz HEME-KatG and ZINC65407173-KatG as well as virtual drug screening were performed to identify potent inhibitor which may be used as antitubercular compound. The graphs of Figure 4 represented molecules A) HEME and B) ZINC65407173 bound inside the active site of modelled protein KatG which was displayed in ‘Cartoon’ and ligands were shown in ‘CPK’ (Corey-Paulin-Koltun) representation (Corey and Pauling, 1953; Koltun, 1965). The equilibrium ‘RMSD’ value for both protein complexes stood around 0.4nm. The effect of ‘HEME’ and ZINC65407173 ligands conformational changes on KatG was observed from the fluctuation in residues as shown in ‘RMSD’, ‘RMSF’, ‘Radius of gyration’ (R_g), ‘SASA’ and ‘Free energy of solvation’ plots (Figure 5).

The SASA plot displayed similar results (320 nm² as average value) for packing distribution of residues in respect of both HEME-protein and ZINC65407173-protein complexes and it showed insignificant conformational changes of complexes.

The average ‘R_g’ value of ZINC65407173 and HEME were recorded as 2.8 nm and 2.9 nm, respectively.

KatG-HEME and KatG-ZINC65407173 complexes were also analysed based on H-bonds % occupancy (Figure 6) during MD simulation. In respect of ZINC65407173 ligand, amino acid residues namely TRP82 and HIS245 were found common with ‘HEME’ active site after screening of (%) occupancy. The interaction of ZINC65407173 with TRP82 and HIS245 residues were stable for 3.2 and 64.60% duration of simulation, respectively while the interaction of ‘HEME’ ligand with TRP82 and HIS245 residues were for 2.59 and 24.13% duration, respectively. The interaction of ‘HEME’ ligand with Lys249 residue in the binding loop between protein and HEME was most stable for 70.18% duration of simulation. The other amino acid residues involved in stable binding included TRP82, LYS188, HIS245, HIS251, THR289, SER290 and LYS532.

Figure 2:
A. Virtual screening of 2632 compounds. B. Identified leads like compounds. C. proposed 3 lead compounds.

The increased average ΔG binding free energy (Table 3 and Figure 5F) in case of ‘HEME’ and ZINC65407173 were reported as -40.13 Kcal/mol and -22.64 Kcal/mol, respectively which denoted their ability to interact better with modelled protein i.e. KatG protein rather than the initial docked pose.

Figure 3:
Docked compounds incorporating modelled protein: A. For ZINC65407173, B. For ZINC09849203. C. For ZINC09349296.

Figure 4:
Graphical representation of molecules A. HEME and B. ZINC65407173 bound inside active site of modelled protein where protein is shown in cartoon representation and ligand is shown in CPK representation.

Figure 5:
Graphical representation of plots: A. Protein RMSD, B. Ligand RMSD, C. Protein RMSF, D. Solvent Accessible Surface Area, E. Radius of Gyration and F. Free Energy of Solvation.

Figure 6:
Percent (%) Occupancies of h-bond protein-ligand contacts: A. Modelled protein-HEME complex, B. Modelled protein-ZINC65407173 complex.

In modern era, computational approaches have achieved appreciable attention in the field of identification of drug target and drug candidates (Boeckmannet al., 2003; Linet al., 2020). In the present study, 2682 compounds of ZINC library were screened against KatG protein to identify promising drug candidates. Eventually, ZINC-65407173, ZINC-09849203 and ZINC-09349296 were shortlisted as drug-like molecules based on docking analysis. ADMET (Absorption, Distribution, Metabolism, Excretion and Toxicity) analysis for 3 drug-like molecules were performed. MD simulation was executed to investigate stability of binding conformation and mode of residual interactions between ZINC-65407173 compound and KatG at active site. The predicted results were discussed in this section. It was revealed from investigation (Table 2, Tables S2 and S3) that ZINC-65407173 compound having molecular weight<500 Da and solubility value of -2.182 is soluble in water and it might reach to active site. The Blood-Brain-Barrier (BBB) value (-0.837>-1) justified that ZINC65407173 was poorly permeable and it can’t cross BBB (Tamai and Tsuji, 2000). CYP4501A2 metabolizing enzymes (Di, 2014) is not inhibited by this compound. The CaCO-2 permeability value of 0.18, % value of HIA (65.379), skin permeability value (-2.822) and Biochemical PAINS with no false positive result for ZINC-65407173 compound indicated that this compound has good oral doping (Rasulev and Casañola-Martin, 2018), good intestinal absorption (Pireset al., 2015), high permeability (Potts and Guy, 1992) and no inhibiting property against Renal Organic Cation transporter protein. Compound ZINC-65407173 possesses ‘ds’ value of 0.904, MRTD value of 0.187 log (<0.477) (Stampferet al., 2019), Minnow toxicity (Chenet al., 2020), T. pyriformis toxicity (Yoshioka et al., 1985), Oral rat acute toxicity and hepatotoxicity but this compound is Ames negative (Vijayet al., 2018) and does not show skin sensitization. Lowest CLogP and Druglikeness values denote that ZINC-65407173 is suitable drug-like compound within no-risk range. TPSA value of 113.3 Å (<140 Å) indicates that this compound does not possess poor cell permeability (Gregson and Plowe, 2005). The stability and conformation of protein and ligand complexes were verified employing RMSD, RMSF, R_g, SASA and % occupancies of H-bonds in protein-ligand complexes during simulation. The RMSD value of 0.4 nm denoted that both HEME-KatG and ZINC65407173-KatG complexes were stable. RMSF showed similar residual fluctuation (<0.5 nm) of lower magnitude for both complexes. SASA plot did not show any occurrence of unfolding or folding incident. H-bonds % occupancy also indicated stable interactions between ‘ZINC-65407173’ compound and ‘HIS-245’ residues for 64.60% duration of simulation while interactions of ‘HEME’ with ‘HIS245’ residues were 24.13% duration. The ‘R_g’ measurement denoted a similarity in protein stability and compactness between KatG-HEME and KatG-ZINC-65407173 (Srinivasanet al., 2021). The good binding affinity, minimal deviation in simulation and decent ADMET properties including drug score with no toxicity risk observed for ZINC-65407173 indicated that compound ZINC-65407173 (C₁₆H₁₇N₇O₃) is a potential candidate for KatG inhibition.

The wet lab studies and clinical trials coupled with this in silico approaches are required to be performed to evaluate the reported activities of Compound No.: ZINC65407173 as promising lead against various strains of M. tuberculosis including M. tb_DKC2 strain.