In silico identification of novel chemical compounds with anti-TB potential for the inhibition of InhA and EthR from Mycobacterium tuberculosis

Tuberculosis (TB) continuously pose a major public health concern around the globe, with a mounting death toll of approximately 1.4 million in 2019. The reduced bioavailability, increased toxicity and resistance of several first-line and second-line anti-TB drugs such as isoniazid, ethionamide have necessitated the search for new medications. In this research, we have identified several novel chemical compounds with anti-TB properties using various computational tools like molecular docking analysis, drug-likeness evaluation, ADMET profiling, P450 site of metabolism prediction and molecular dynamics simulation study. This study involves fifty drug-like compounds with antibacterial activity that inhibit InhA and EthR involved in the synthesis of one of the major lipid components, mycolic acid, which is crucial for the viability of Mycobacterium tuberculosis. Among these 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 filtering analysis successfully. Their binding stability for target proteins have been evaluated through RMSD, RMSF, Radius of gyration analysis from 10 ns Molecular Dynamics Simulation (MDS) run. Our identified drugs could be a capable therapeutic for Tuberculosis drug discovery, having said that more in vitro and in vivo testing is required to justify their potential as novel drug and mode of action.


INTRODUCTION
Tuberculosis (TB) is a transmissible disease brought about by Mycobacterium tuberculosis (Mtb), which brings severe damage to lungs (pulmonary TB), along with other regions of the body, i.e. brain, spinal cord, lymph node, abdomen, bones/joints, intestinal and genitourinary system (extra pulmonary TB) [1,2,3]. Active TB symptoms depends on the severity of spread, generally some of clinical features are: persistent coughing more than 3 weeks, bloody sputum (hemoptysis), chest pain, fever, fatigue/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 [4]. 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 [2,7].
Although TB is a treatable infection, the massive emergence of resistance to antibiotics makes it a global threat. The classes of resistance that are becoming more and more prevalent includes (i) RR-TB, resistance to the first line drug rifampin only (ii) MDR-TB or multidrug resistance TB defines as the resistance towards at least two first line drugs, isoniazid and rifampin (iii) XDR-TB or extensively drug resistance TB defines as MDR-TB combined with resistance to a minimum of one drug from every of the two classes of second-line drugs (fluoroquinolones and injectables) and (iv) TDR-TB or totally drug resistance defines as the resistance to more drugs than strains categorized as XDR-TB [5,6].
It is undeniable that, the world needs new drugs to improve current treatment of TB to get around the escalation of drug-resistance. Therefore, our in-silico approach seeks to offset the burden of existing drug resistance scenarios of through exploration of new drugs performance. As a part of the experiment, a total fifty drug-like compounds were collected to assess their effectiveness in the treatment of TB and clinical management. Our experiment proceeded through evaluating binding affinity of fifty ligands utilizing two different docking tools DockThor and Autodock vina individually against the InhA (PDB ID: 3FNG) and EthR (PDB ID: 3G1M) receptor proteins. Later on, leading four compounds were analyzed for predicting probable metabolic sites of CYP. The study was further validated by assessing drug-likeness properties and ADMET prediction so as to examine their pharmacokinetics and pharmacodynamics properties. Then the best compounds with receptor complex were subjected to molecular dynamics simulations to test the thermodynamic properties.

InhA (Enoyl acyl carrier protein reductase), a pivotal enzyme used in Mycolic acid biosynthesis pathway:
Mycobacteria are distinctive in having two Fatty acid synthase (FAS) systems which ends up synthesizing the mycobacterial mycolic acid (MA) precursors [8]: the "eukaryotic-type" multifunctional FAS-I for de novo synthesis of long-chain fatty acyl-CoA and their subsequent elongation by the "bacterial-type" FAS-II [9,10]. Malonyl-CoA is produced by the carboxylation of acetyl-CoA catalyzed by acetyl-CoA carboxylase (ACC) and integrated into the developing acylic chain during the repeated cycle of reactions of fatty acids synthase I and II (FAS I/FAS II) [11]. The mero-mycolic chain precursors are synthesized in four stages of iterative cycles with four FAS-II enzymes catalyzing acyl carrier protein (ACP) derivatives in each elongation cycle [8]. The transformation of malonyl-CoA to malonyl-ACP by malonyl-CoA:ACP transacylase (mtFabD) facilitates its integration into the FAS-II enzymatic system [12]. Both FAS systems are interconnected by a CoA-dependent linker enzyme β-ketoacyl-ACP synthase III (mtFabH) which catalyzes Claisen condensation of acyl-CoA precursors produced by FAS-I and malonyl-ACP to form β-ketoacyl-ACPs [13]. Next, β-ketoacyl-ACP reductase (MabA) performs CoA, later on the synthesis of malonyl-ACP occurs by malonyl-CoA:AcpM transacylase (mtFabD). FAS II is triggered by the β-ketoacyl-ACP synthase (mtFabH), which condenses the acyl-CoA to malonyl-ACP to generate a β-ketoacyl-ACP. Following the condensation stage, the β-ketoacyl-ACP passes through a series of keto-reduction, dehydration, and enoyl-reduction catalyzed by the β-ketoacyl-ACP reductase (MabA), the β-hydroxyacyl-ACP dehydratase complex (HadABC) and the enoyl-ACP reductase (InhA), respectively. Further elongation cycles are initiated either by KasA or KasB β-ketoacyl ACP synthase. The last condensation carried out by polyketide synthase (Pks13) delivers trehalose mycolic β-ketoester and then reduced to trehalose monomycolate (TMM) [ADAPTED FROM 19].

Initial Molecular Docking in DockThor program
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. This web-based tool spontaneously produces the topology files (i.e., atom types and partial charges) for the protein, ligand, and cofactors according to the MMFF94S force field [29]. 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 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 drug 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 [30] and then minimizing the structure employing Swiss-PdbViewer [31].
Following the initial affinity prediction, desired ligands with lower affinity score than control drugs were chosen to evaluate bound conformation and binding energy. Autodock vina software (AutoDock 4.2) was used to evaluate bound conformation and binding energy of those selected ligands [32]. This tool uses Lamarckian Genetic Algorithm to evaluate binding energy of ligands with receptor proteins. Exhaustiveness defines how many times the system will repeat the computation on binding site of protein and this was kept 10 for the analysis. CASTp server [33] was used to identify ligand binding sites of proteins and prepare grid box.

Visualization and interaction analysis
The 2D and 3D visualization of non-bonded interactions of protein-ligand docked complexes were done by BIOVIA Discovery Studio 4.1 Visualizer [34]. This tool was utilized to get the total number of hydrogen bonds and interacting amino acids of proteins.

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. These properties were predicted utilizing Lipinski's rule of five [36], Ghose's rule [37] Veber's rule [38] Muegge's rule [39], TPSA and No of rotatable bonds. The calculation was carried out using SwissADME online tool (http://www.swissadme.ch/index.php) [40].
Absorption, distribution, metabolism, excretion and toxicity, all these chemical properties were also crucial determinant of drug like compounds. ADMET properties of each of four compounds were predicted using another online tool, admetSAR (http://lmmd.ecust.edu.cn/admetsar2/) [41].
In each of the prediction study, canonical smiles of the compounds were used from PubChem database (https://pubchem.ncbi.nlm.nih.gov/).

Molecular Dynamics Simulation
The study of molecular dynamics simulation (MDS) is a thermodynamics-based operation which helps to study the dynamic perturbation found in 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_Win64-multicore-CUDA version [42] implying CHARMM 36 force field [43] and TIP3P water model. A multi-step time algorithm was used, with an integration time step of 2 femto second. Visual molecular dynamics (VMD) [44] was used to generate psf files of protein-ligand complexes, water box and for neutralizing the system with sodium (Na+) and a chloride (Cl-) ions. Ligand topology and parameter files were generated using CHARMM-GUI web service [45]. 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 310k. Periodic boundary conditions was applied surrounding the system. Finally, the results were visualized and analyzed in VMD.

Molecular docking analysis
To filter out the best ligands from docking analysis, two separate docking tools were used with different algorithm. Initially used DockThor (server-based docking tool) results were further scanned through Autodock-vina software. Detailed information of all the collected drug-like compounds were provided in Supplementary material.
Total fifty ligands were primarily docked individually against the InhA (PDB ID: 3fng) and EthR (PDB ID: 3g1m) receptors protein via the DockThor server. Against both 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.    Table 2.

P450 Site of Metabolism (SOM) prediction
Cytochromes P450 plays an important role in metabolizing external materials like drugs. The  Table 3.

3.3.Drug-likeness and ADMET analysis
The concept of drug-likeness has appeared as an approach that can screen the high-affinity ligands with acceptable ADME (absorption, distribution, metabolism, excretion) properties. In our study,  (Table 4). 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 chosen molecule was capable of pervading Caco2 cell lines.
C22 acted as P-glycoprotein substrate, P-glycoprotein inhibitor and substrate for CYP3A4 cytochrome.

3.4.Molecular Dynamics Simulation:
In our study, Molecular Dynamics Study was performed to assess conformational stability of protein-ligand complexes as well as InhA and EthR receptors. The overall conformational similarity between control drugs and the drugs being tested with target proteins was compared assessing Root Mean Square deviation (RMSD). As represented in Figure 5

DISCUSSION:
In our study, the screening of high affinity ligands was accomplished by two stages of molecular docking analysis. As an initial docking evaluation, fifty ligands were primarily docked separately against the InhA and EthR receptors protein via the DockThor server [29]. Both for InhA and EthR proteins, all the compounds showed lower affinity score than their respective controls.
Subsequently Autodock Vina software was employed for computing binding energy and bound conformation in which two ligands for InhA protein and two ligands for EthR protein were chosen considering their lower binding score [32].

F
Drug-likeness is an approach that qualitatively evaluate the probability for a molecule to qualify as an oral drug with respect to bioavailability [46,47]. To analyze the drug-likeness property of our compounds we used SwissADME, a free web-based tool [40]. Lipinski's Rule of 5, a rule of thumb was established to evaluate 'drugability' of new chemical entities having certain pharmacological or biological activity [48]. In the drug discovery and optimization approach, the Rule of 5 reveals that a drug compound should follow the criteria with the acceptable range having According to Ghose filter, to qualify as a drug, it should have log P value between -0.4 and 5.6, molecular weight from 160 to 480, molar refractivity between 40 and 130, total atomic count between 20 and 70 [37]. Veber rule states that compounds which meet two criteria of: 10 or fewer rotatable bonds and polar surface area equal or less than 140 Å are predicted to have good oral bioavailability [38]. The Egan rule considers candidate drug to have good oral bioavailability with −1.0 ≤logP ≤5.8 and TPSA ≤130 Å2 [49]. Muegge rule utilizes pharmacophore point filter to distinguish between drug-like and nondrug-like compounds which is based on simple structural rules. To pass the filtering, candidate drug should obtain two to seven pharmacophore points [39].
In addition to following the Lipiniski rule, they followed the Ghose, Veber, Egan and Muegge rules as well. The aqueous solubility of a compound expressed as log S determines the oral bioavailability resulting in affecting its absorption and distribution in the ADME/T properties [50].
The "Topological Polar Surface Area" (TPSA) is an approach that enables the calculation of the Polar Surface Area which is correlated with passive molecular transport through membranes, thus used to predict the transport properties of medicines such as human intestinal absorption, permeability to Caco-2 monolayer and permeation of blood-brain barrier [51]. Bioavailability results define the probability compounds to have minimum value of 10% oral bioavailability in rat or measurable Caco-2 permeability [52]. Considering the standard values, our four compounds pass through the Log S, bioavailability score, TPSA (Å²), Number of rotatable bonds value. The drug-likeness property experiment showed that our four compounds possess drug-like properties against tuberculosis.
Cytochrome P450 (CYPs) enzymes are metabolic enzymes, responsible for the biotransformation of ∼90% FDA certified drugs [53]. The oxidative metabolism of drugs in phase I is achieved by the CYP system [54,55]. Nine of the isozymes under the CYP system scanned for the prediction of the metabolically vulnerable points using RS-WebPredictor tool [56]. It is a publicly available, web-based tool for predicting the regioselectivity of isozyme specific CYP mediated xenobiotic metabolism on any set of user-submitted molecules [57]. 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.
In silico analysis of ADMET properties, absorption, distribution, metabolism, elimination and toxicity can help in determining the efficacy, pharmacological profile, mode of administration and safety of drugs. [58]. 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 probable drug candidates showed ability of penetrating Caco-2 cell line. P-glycoprotein (P-gp) is an efflux, ATP powered membrane-bound transport protein pump that is widely distributed throughout the body. It is a unique defensive barrier network that protects tissues from toxic xenobiotics by blocking entry to the cytosol and pumping its substrates to the exterior and also determines the uptake and distribution of pharmacotherapeutic drugs [61]. 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. Blood-Brain Barrier, a unique property possessed by the blood vessels in the Central Nervous System (CNS) tightly controls the transport of molecules, ions, nutrients and cells from blood to the brain and vice versa [62]. Consideration of BBB is an utmost priority for those drugs whose primary target is brain cells. Out of four ligands, only C38 was not capable of penetrating BBB.
As the incidence of anti-TB induced hepatotoxicity is one of the frequent causes behind the termination of anti-TB medication, assessing the potential hepatotoxicity of novel drugs is important [63]. As well as carcinogenicity test is necessary to classify a tumorigenic potentiality of drugs to assess the relevant risk in humans [64]. Only two ligands, C22 and C29 passed two of the criteria: -hepatotoxicity, and carcinogenicity as they showed negative result.
To validate 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 [65]. Our Radius of Gyration (Rg) analysis disclosed our proposed drug C29 and C22 induced fewer fluctuation of Rg upon binding with their respective receptors compared to the control drugs. This confirmed C29 and C22 did not cause instability to receptors. We compared our proposed drugs C29 and C22 with

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 candidate 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 pharmacore 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