In silico identification of novel chemical compounds with antituberculosis activity for the inhibition of InhA and EthR proteins from Mycobacterium tuberculosis

Tuberculosis (TB) continuously poses a major public health concern around the globe, with a mounting death toll of approximately 1.4 million in 2019. Reduced bioavailability, elevated toxicity, increased side effects, and resistance of multiple first-line and second-line TB medications, including isoniazid, ethionamide necessitate studies of new drugs. The method of computational biology and bioinformatics approach allows virtual screening of a large number of drugs, reduces growing side effects of medications, and predicts potential drug resistance over time. In this study, we have analyzed fifty small molecules with antituberculosis properties using in silico approach including molecular docking, drug-likeness assessment, ADMET (absorption, distribution, metabolism, excretion, toxicity) profile evaluation, P450 site of metabolism prediction, and molecular dynamics simulation. Among those fifty compounds, 3-[3-(4-Fluorophenyl)-1,2,4-oxadiazol-5-yl]-N-(2-methylphenyl) piperidine-1-carboxamide (C22) and 5-(4-Ethyl-phenyl)-2-(1H-tetrazol-5-ylmethyl)-2H-tetrazole (C29) were found to pass the two-step molecular docking, P450 site of metabolism prediction and pharmacokinetics analysis successfully. Their binding stability for target proteins has been evaluated through root mean square deviation and root mean square fluctuation, Radius of gyration analysis from 10 ns Molecular Dynamics Simulation (MDS). Our identified drugs (C22 and C29) performed better than the control drugs (Isoniazid, Ethionamide) regarding binding affinity and molecular stability with the regulatory proteins (InhA, EthR) of Mycobacterium tuberculosis. The study proposed these compounds as effective therapeutic agents for Tuberculosis drug discovery, but further in vitro and in vivo testing are needed to substantiate their potential as novel drugs and modes of action.


Introduction
Tuberculosis is among the earliest known infectious diseases [1]. Since the introduction of Mycobacterium tuberculosis as the causative agent of tuberculosis by Robert Koch in 1882, it has grown and evolved as one of the fatal pathogens in the history of infectious disease in humans [2]. Mycobacterium tuberculosis is a distinct acid-fast, slowly developing, aerobic bacteria with a size of 0.  µm that belongs to the Mycobacteriaceae family [3]. The most regularly studied strain H37Rv of M. tuberculosis has a genome composed of 4.4 Mb that encodes 4018 genes with an average G + C content of 65.6% [4].
Tuberculosis (TB) is a communicable disease, transmitted by a person suffering from active TB with the release of small nuclei respiratory droplets (ranges between 0.65 to greater than 7.0 μm) containing viable airborne bacteria [5,6]. Once the bacilli reach the lung alveoli, either the mycobacteria are phagocytized by mature alveolar macrophages, or the active mycobacteria start reproducing within the macrophage, which then causes their lysis [7]. This local inflammation attracts monocytes from neighboring blood vessels and leads to granuloma formation by encapsulating bacilli for controlled progression (LTBI) [8]. In case of failure, the infection starts to spread through both lymphatic channels and blood circulation which can cause serious damage to the lungs (pulmonary TB), along with other regions of the body, i.e., brain, spinal cord, lymph node, abdomen, bones/joints, intestinal and genitourinary system [7,[9][10][11]. Active TB symptoms depend on the severity of the spread; generally, some of the clinical features are persistent coughing for more than 3 weeks, bloody sputum (hemoptysis), chest pain, fever, fatigu e/weakness, weight loss, anorexia, breathlessness, etc. The person with latent tuberculosis infection (LTBI) presents no symptoms and cannot transmit Mtb to others, but may develop active TB [12].
One of the major barriers to therapeutic interference against Mtb is its unusual cell envelope structure that has a high lipid content with unique long-chain highly hydrophobic α-alkyl-β-hydroxy fatty acids, mycolic acids (MAs) [13,14]. MA functions to maintain the structural integrity and viability of the mycobacterium by acting as an extraordinarily efficient permeability barrier and contributes to the intrinsic resistance towards host bactericidal agents or different classes of antibiotics [15][16][17][18]. It also protects against oxidative stress, elicits differentiation of macrophages into foamy macrophages [19,20]. Therefore, in our in-silico approach, we chose to target a pivotal enzyme, enoyl-acyl carrier protein reductase (InhA) involved in the synthesis of MA [21] and EthR, a negative transcriptional regulator which reduces the expression of EthA enzyme [22]. The inhA gene encodes NADHdependent enoyl-acyl carrier protein (ACP)-reductase enzyme which is associated with long chain fatty acids (mycolic acids) synthesis and reduction phase of fatty acid production. Mutation of the inhA gene facilitates resistance to several first-line drugs such as isoniazid, making it an appropriate target for drug discovery [23]. EthR belongs to the tetR/CamR family of transcriptional repressor that negatively regulates the synthsis of EthA enzyme which is involved in the activation of thiocarbamide-bearing drugs ethionamide. Low activation of ethionamide by the EthA enzyme increases the risk of Mycobacterium tuberculosis becoming resistant to this type of drugs [22,24]. This mechanism implies its significance as therapeutic target for resistance to these drugs.
WHO recommended a standard 6-month administrative treatment for active, pulmonary, and pharmaco-sensitive TB with different combinations of four first-line medications: isoniazid, rifampin, ethambutol, and pyrazinamide as bases for treatment [10,25]. Unfortunately, the complex long-term administrative duration leads to non-compliance of treatment leading to failure and relapse which causes the bacterium to become resistant over time to anti-TB medications [26]. The antibioticresistant strains of Mycobacterium tuberculosis are sabotaging conventional tuberculosis therapy due to resistance to a wide range of anti-TB medications. Primary and secondary resistance contributes to the transfer of resistant strains to newer hosts and the development of drug resistance to two or more drugs respectively [27]. Resistance in this bacterium includes overexpressed drug targets due to mutation of Fig. 1. Complete methodology of this study in a concise flowchart.

S.K. Halder and F. Elma
repressor or promoter region of targets, alteration of drug targets by enzymatic process, cancelation of mechanism of prodrugs such as Ethionamide that induces inactivation of first-line drugs like Isoniazid, direct inactivation of drugs within the bacterial cell by enzymatic process, drug efflux process through which drugs are transferred outside the bacterial cell [28]. The onset of MDR-TB or multidrug resistance TB exacerbated the treatment of TB because of its poor efficacy, toxicity also the extra cost of about 20 months of extended administration of second-line drugs [29]. Furthermore, in 2016 an estimated 6.2% incidence of extensively drug-resistant TB (XDR-TB) among people with MDR-TB cases around the world has made the treatment more complex with limited efficacy of existing drugs [30,31].
Although TB is a treatable infection, the massive emergence of resistance to antibiotics makes it a global threat. It is undeniable that the world needs shorter and simpler yet affordable, effective, tolerable, and safe new drug regimens to improve the current treatment of TB to get around the escalation of drug resistance. Our in-silico approach seeks to offset the burden of existing drug resistance scenarios through the exploration of new drug performance. As a part of the experiment, a total of fifty drug-like compounds were collected to assess their effectiveness in the treatment of TB and clinical management. Our experiment proceeded through evaluating the binding affinity of fifty ligands utilizing two different docking tools DockThor server and Autodock vina software individually against the InhA (PDB ID: 3FNG) and EthR (PDB ID: 3G1M) receptor proteins. Identification of metabolic sites and pharmacokinetics and pharmacodynamics properties of tested compounds contributed to confirm their drug-like properties. Molecular dynamics simulations were conducted to determine the molecular stability of the best drug-like compounds with respective receptor complexes.

Initial molecular docking in DockThor program
The crystal structures of two receptors i.e., Crystal structure of InhA bound to triclosan derivative (PDB ID: 3FNG) and EthR from Mycobacterium tuberculosis in complex with compound BDM31381 (PDB ID: 3G1M), were taken from Protein Data Bank (https://www.rcsb.org/ search). Fifty compounds were collected with anti-TB properties from both ACD: Antibacterial Chemotherapeutics Database (http://amdr. amu.ac.in/acd/index.jsp) and PubChem (https://pubchem.ncbi.nlm. nih.gov/) servers. Isoniazid (INH) and Ethionamide (ETH) were used as control drugs against InhA and EthR proteins respectively. In the beginning, the PDB structures were modified using PyMOL tools (PyMOL) by clearing the water molecules from the structure [32] and then minimizing the structure employing Swiss-PdbViewer [33]. The preliminary docking of two receptor protein and fifty ligands were carried out by DockThor docking server (https://www.dockthor.lncc. br/v2/). The algorithm of the program is based on flexible-ligand and rigid-receptor grid-based system [34].

Autodock-vina binding affinity prediction
Following the initial affinity prediction, desired ligands with lower affinity scores than control drugs were chosen to evaluate the bound     [36] was exploited to determine ligand binding sites of proteins and calculate grid box size.

Drug-likeness properties analysis and ADMET prediction
Determination of drug-likeness property of drug-like compounds is one of the vital steps of drug discovery. Prediction of these properties was completed utilizing Lipinski's rule of five [38], Ghose's rule [39], Veber's rule [40], Muegge's rule [41], TPSA, and No of rotatable bonds. The calculation was carried out using SwissADME online tool (http ://www.swissadme.ch/index.php) [42].

Visualization and interaction analysis
Visualization of the non-bonded relationship of 2D and 3D conformation of protein-ligand docked complexes was done by BIOVIA Discovery Studio 4.1 Visualizer [45]. This tool was utilized to get the total number of hydrogen bonds and interacting amino acids of proteins with respective ligands.

Molecular dynamics simulation
The study of molecular dynamics simulation (MDS) is a thermodynamics-based operation that helps to study the dynamic perturbation found in the protein-ligand complexes. In our experiment, to ensure the stability of protein-ligand complex, we subjected the best ligands screened from previous steps to the molecular dynamics simulation (MDS) study with their respective proteins. We simulated the docking complexes using the NAMD_2.14bNAMD_2.14b2_Win64multicore-CUDA version [46] implying CHARMM 36 force field [47] and TIP3P water model. A multi-step time algorithm was used, with an integration time step of 2 femto seconds. Visual molecular dynamics (VMD) [48] was used to generate psf files of protein-ligand complexes, water box and for neutralizing the system with sodium (Na+) and chloride (Cl-) ions. Ligand topology and parameter files were generated using the CHARMM-GUI web service [49]. The simulation was run for 10 ns, where the system was minimized for 1000 steps. Langevin thermostat was used to maintain a constant temperature of 310 k. Periodic boundary conditions were applied surrounding the system. The complete workflow of our methodology was summarized in Fig. 1.

Molecular docking analysis
To filter out the best ligands from molecular docking analysis, two separate docking tools were used with a different algorithm. Initial docking results of the DockThor server were further scanned by Autodock-vina software. Detailed information of all the collected druglike compounds was provided in Supplementary material.
Total fifty ligands were primarily docked individually against the InhA and EthR receptors via the DockThor server. Against target proteins, InhA and EthR, all the compounds showed better affinity than control 1 (Isoniazid) and control 2 (Ethionamide). The binding affinity score between InhA, EthR receptors, and fifty ligands, as well as control drugs, is depicted in Table 1.
After initial docking from the DockThor server, all the ligands were considered for computing binding energy and bound conformation with the assistance of Autodock Vina software. Two ligands were preferred over all ligands considering their lowest binding score towards InhA protein. Likewise, two ligands with the lowest binding score towards EthR protein were chosen for further analysis. Comparative Docking results of the Autodock-vina tool are listed in Table 2.
The relative ADMET profiles of screened ligands are described in Table 5. The evaluated pharmacokinetic data showed that all of the selected molecules had a high intestinal absorption rate and oral bioavailability. Each of the chosen molecules were capable of pervading Caco2 cell lines. C22 acted as P-glycoprotein substrate, P-glycoprotein inhibitor, and substrate for CYP3A4 cytochrome.
None of the ligands acted as CYP2D6 substrate or CYP2D6 inhibitor. Only C29 showed an inhibitory effect on CYP2C9 cytochrome. C22 acted as an inhibitor for CYP3A4, CYP2C19, CYP1A2 cytochrome but acted as a substrate for CYP2C9 cytochrome. C21 inhibited CYP3A4, CYP2C19, CYP1A2 cytochrome isoform but represented as a substrate for CYP3A4 isoform.
All selected ligands were non-carcinogenic. C22 and C29 were nonhepatotoxic. The best ligands and control drugs with respective receptors (C22 and InhA complex, Isoniazid and InhA complex, C29 and EthR complex, Ethionamide and EthR complex) were depicted in Figs. 2 and 3.  Table 4 The Drug-Likeness properties of the best four ligands.

Molecular Dynamics Simulation:
In our study, Molecular Dynamics Study was performed to assess the conformational stability of protein-ligand complexes as well as InhA and EthR receptors. The overall conformational similarity between drugs (control, C22, and C29) with target proteins was compared assessing Root Mean Square deviation (RMSD). As represented in Fig. 4  (A), RMSD of Free InhA protein c3ontinued stability between 1 ns and 5 ns timescale at about 1.4 Å and 5 ns (nanosecond) to 9 ns at about 1.7 Å. (Isoniazid + InhA) complex displayed steadiness in RMSD between 1 ns and 5 ns timescale at around 1.5 Å and 8 ns to till the end of run at around 1.8 Å. (C22 + InhA) complex held stable backbone stability from 3 ns to 7 ns, with RMSD around 2 Å. Free EthR protein, (Ethionamide + EthR) protein-ligand complex and (C29 + EthR) complex showed a backbone RMSD under 2 Å, which suggested minor structural changes. Free EthR protein remained stable between 1 ns and 6 ns, exhibiting RMSD around 1.2 Å. (Ethionamide + EthR) complex also showed stability from 1 ns to 6 ns timeline, with an RMSD backbone of 1.2 Å, then exhibiting small increased RMSD. On the other hand, RMSD of (C29 + EthR) complex persisted stable from 1 ns to 5 ns timeline at about 1.5 Å, then after slight fluctuation, remained stable from 8 ns to the end of runon average 1.4 Å.
Root Mean Square Fluctuation (RMSF) analysis per residue for backbone atoms was conducted to assess changes in the conformation of Cα backbone of the systems. From Fig. 5(C), Free InhA, (Isoniazid + InhA), (C22 + InhA) complexes mostly had RMSF from 0.4 to 2 Å that indicated close conformational contact between protein and ligands. Nevertheless, the higher fluctuation of RMSF between 205 and 211 residues confirmed the presence of loop within this region. Our MDS study showed that the RMSF of Free EthR, (EthR + Ethionamide), (EthR + C29) complexes fluctuated mostly between 0.5 and 1.5 Å, indicating close contact between the active pocket of receptors and drugs. However, higher fluctuation from 73 to 75 and 170 to 172 amino acid residues indicated that the free protein and its complexes were within the loop regions.
Radius of Gyration with time was calculated to assess the change of compactness after ligand binding with receptors. From Fig. 6(E), the Rg of Free InhA was reported between 17.9 and 18.2 Å; (InhA + Isoniazid)

Discussion
In our study, virtual screening of high-affinity ligands was accomplished by two stages of molecular docking analysis. As for the initial docking evaluation, fifty ligands were docked separately against the InhA and EthR receptor proteins through the DockThor server [34]. Both InhA and EthR proteins are associated in vital function in which InhA is involved in the synthesis of mycolic acids, essential part of mycobacterial cell wall and EthR, a repressor of monooxygenase EthA, brings about resistance to Ethionamide [50,51]. For both InhA and EthR proteins, both drug-like compounds had a lower affinity score than their respective control drugs. Subsequently, Autodock Vina software was used for calculating binding energy and bound conformation in which two ligands for InhA protein and two ligands for EthR protein were selected based on their lower binding score [35].
Cytochrome P450 (CYPs) enzymes are metabolic enzymes, responsible for the biotransformation of ~90% FDA certified drugs [52]. The oxidative metabolism of drugs in phase I is achieved by the CYP system [53,54]. Nine of the isozymes under the CYP system scanned for the prediction of the metabolically vulnerable points using RS-WebPredictor tool [37]. This tool helps to predict the regioselectivity of isozymespecific CYP mediated xenobiotic metabolism on any set of usersubmitted molecules [55]. C22, C38, C21, C29 drugs displayed several sites of metabolism for CYP1A2, 2A6, 2B6, 2C8, 2C19, 2E1, 3A4, 2C9, 2D6 isoforms that suggested a satisfactory result.
Human Intestinal Absorption (HIA) is a pharmacokinetic process that determines the effectivity of intestinal absorption or bioavailability of a drug upon oral administration, an anticipated route of drug administration [59]. Our four ligands showed a high rate of intestinal absorption and oral bioavailability. Caco-2 cell line, a prominent substitute to human intestinal epithelium (mucosa) is one of the in-vitro models to determine in vivo human intestinal absorption of drug molecules due to their morphological and functional resemblances with human enterocytes [60]. All of our potential drug candidates showed the ability to penetrate Caco-2 cell line. C22 acted as P-gp substrates as well as inhibitors. On the other hand, C38, C21, C29 were neither P-gp substrates nor P-gp inhibitors.
Since the occurrence of hepatotoxicity caused by anti-TB is one of the frequent reasons behind the termination of anti-TB products, it is important to evaluate the possible hepatotoxicity of novel drugs [61]. As well as a carcinogenicity test is necessary to classify a tumorigenic potentiality of drugs to assess the relevant risk in humans [62]. Only two ligands, C22 and C29 passed two of the criteria: hepatotoxicity, and carcinogenicity as they showed negative results.
To validate the conformational stability of our proposed drug upon binding with receptors, Molecular Dynamics Simulation was performed. Root Mean Square deviation (RMSD) analysis showed that binding of both control drugs and proposed drug C29 did not induce structural instability to the proteins because the reported RMSD change after ligand binding remained under 2 Å for both receptors. A thorough study of Root Mean Square Fluctuation (RMSF) curve revealed our tested compounds kept close contact with their active sites which was evident from their small range fluctuation under 1.5 Å and 2 Å for EthR and InhA complexes respectively. Though, because of loop regions on receptor proteins, greater fluctuation of RMSF was seen [63]. Radius of Gyration analysis showed that our proposed drug C29 and C22 caused less fluctuation when bound to their respective receptors relative to control drugs. This confirmed C29 and C22 did not cause instability to receptors. We compared our proposed drugs C29 and C22 with control drugs (Ethionamide and Isoniazid) through Molecular Dynamics Simulation (MDS), in which C29, C22 protein complexes represented stability throughout 10 ns simulation.
Despite the fact that drug repurposing earned tremendous success, there are some caveats to the in-silico methodology. Like, one of the downsides of molecular docking is the provision of appropriate scoring functions and algorithms, which may otherwise jeopardize molecular screening [64]. Although our entire approach was based on computational tools, the findings were validated by dynamic molecular simulation and pharmacokinetic profiling of drug-like compounds.
Our study aims to identify the linkage between tuberculosis and already established antibacterial drugs applying several bioinformatics tools. This methodology can evade the experimental hurdles of screening thousands of ligands for tuberculosis. Integration of in silico and experimental approaches led to discovery of Lidocaine, Methotrexate, Mifepristone and Zidovudine as a potent therapeutic against Arrhythmia, Arrhythmia, Cushing's syndrome and HIV (human immunodeficiency virus) respectively [65,66]. To prove the authenticity of our findings, in vivo research involving animal model is required. The curated dataset can be a great deal of interest to those researchers working in this field to find novel anti-TB drugs.

Conclusion
The number of people dying annually of TB is growing rapidly due to multiple drug resistance scenarios around the world. This situation demands newer anti-TB drugs to deal with the crisis. Drug repurposing is an easier and cheaper option to look for novel candidates as anti-TB drugs using different computational tools. The main objective of this study is to a find novel inhibitor against fast-mutating anti-TB drug targets. Our work incorporates pharmacophore analysis, ADMET profiling, two-step molecular docking, followed by 10 ns Molecular Dynamics Simulation. Drugs with greater binding affinity than the control drugs are considered for determining Drug-likeness and ADMET analysis to evaluate their bioavailability and toxicity. Two of our screened compounds: 3-[3-(4-Fluorophenyl)-1,2,4-oxadiazol-5-yl]-N-(2-methylphenyl) piperidine-1-carboxamide and 5-(4-Ethyl-phenyl)-2-(1H-tetrazol-5-ylmethyl)-2H-tetrazole showed promising results with higher binding affinity with respective receptors and standard pharmacophoric properties. Molecular Dynamics Simulation study including RMSD, RMSF, Rg analysis confirmed their binding stability with respective proteins throughout the simulation timeline. Our present work could be productive in discovering potential therapeutics against multiple drug resistant tuberculosis, stating that substantial in vitro and in vivo experiments are needed to prove our hypothetical study.