Facing Antitubercular Resistance: Identification of Potential Direct Inhibitors Targeting InhA Enzyme and Generation of 3D-pharmacophore Model by in silico Approach

Introduction

Tuberculosis (TB) is still sounding the alarm of emergency worldwide, according to the World Health Organization’s 2022 global tuberculosis report, a total of 10.6 million people fell ill with tuberculosis in 2021 and 1.6 million people died from it (including 187,000 among HIV positive).¹ TB is the 13th leading cause of death and the second leading infectious killer after COVID-19 (above HIV/AIDS) (World Health Organization, Global tuberculosis report 2021), The WHO recommended treatment for pulmonary TB (HIV-negative patients) is six months of quadruple therapy: two months of four-drug regimen (isoniazid, rifampicin, pyrazinamide, and ethambutol), followed by four months of isoniazid and rifampicin with once-daily or three times weekly dosing.²

Tuberculosis (TB) aggravates the situation with the problem of drug-resistance, multidrug-resistant TB (MDR-TB) and extensively drug-resistant TB (XDR-TB) and new pharmaceutical agents are urgently required to control this situation.

Using genetic methods, Jacobs and his team identified the protein product of the mycobacterial INHA gene as the putative target for isoniazid.³ InhA (Enoyl-[acyl-carrier-protein] reductase, NADH), catalyses the NADH-dependent reduction of the double bond of 2-trans-enoyl-[acyl-carrier protein], an essential step in the fatty acid elongation cycle of the FAS-II pathway.⁴

Since its first discovery in 1952, isoniazid (INH) has been a cornerstone of tuberculosis chemotherapy, showing a potent bactericidal activity against MTB.^5,6 This drug causes acid-fastness loss by targeting InhA protein after its activation by KatG protein due to the inhibition of mycolic acid synthesis, a-alkyl-b-hydroxyl long-chain fatty acids (60–90 carbons long), a major component of the mycobacterial cell wall.^7,8

Knowing that isoniazid itself does not directly inhibit the enzyme, but it must be converted first by catalase peroxidase (KatG) into the active form of the drug^9,10 and then becomes covalently linked to the nicotinamide ring of the NADH and form in return INH-NAD complex, this complex inhibits the action of the enzyme and prevents the formation of the protein responsible for the formation of mycolic acids, with the problem of resistance, mutations detected in the active site of the target protein, make the activation of isoniazid more difficult and unreachable,¹¹ the need to explore INH-NAD analogues is strongly required, which can target directly the InhA enzyme without needing to be activated first.

Figure S1 shows the explanatory pathway of isoniazid activation as well as the role of the InhA protein in mycolic acid synthesis, it also shows the scope of our problem, and the need for direct inhibition of this protein.

Drug discovery is a veritably expensive and time-consuming process, and to discover a new medicine it may take several years, this task has become more accessible using computer-aided drug design (CADD) approaches like a very powerful tool to simulate molecular interactions of large systems. It provides valuable insights on experimental results and the mechanism of action, predict new molecules to synthetize, and most important it can help to make cost effective decisions before starting expensive molecule synthesis. Many compounds were discovered and optimized using computer-aided drug design methods or even gained US FDA approval.¹² Among the techniques used are virtual screening analysis to discover and optimize novel potent molecules by using molecular docking,^13–16 the computational prediction methods in order to predict the impact of mutation on protein properties,¹⁷ toxicity prediction has also been done on the basis of the molecule criteria which show the physicochemical properties and biological activity of the selected drugs and pharmacophore modelling for the classification of the key features of the best compounds found.¹⁸

The aim of this study is to develop direct InhA protein inhibitors based on the understanding of isoniazid mechanism using CADD like a powerful technology, for this purpose, three different approaches were used: molecular modelling including a virtual screening to determine the best leads with the lowest bindings energy against InhA protein that still have a nice binding affinity score with the mutant models generated; computational prediction method in order to predict the impact of all the mutations described in the literature on protein properties; and pharmacophore mapping to discover common features between the best ligands found.

Materials and Methods

Workflow

Figure 1 shows the workflow of the present study.

Figure 1 In silico study workflow.

Dataset Generation and Preparation of Wild-type and Mutated Models of InhA Protein

The wild-type X-ray structure of InhA in complex with the isonicotinic inhibitor-ACYL-NADH (PDB ID: 1ZID) with a resolution of 2.70 Å was extracted from the Protein Data Bank,¹⁹ and the attached ligand was removed. Fifteen known InhA mutations were collected from literature. These mutations are: D148G, I16T, I21V, I47T, I95P, I194T, K8N, K165A/Q, S94A, T266A/D, V78A, and Y158A/S.^4,20–25 Then, a mutant model was generated for each of them using “Rotamer” function of UCSF Chimera tool V1.16.²⁶ Finally, the protein structures (including mutated models) were prepared using AutoDock Tools – V.1.5.7,²⁷ by removing water molecules and adding polar hydrogens and Kollman charges then saved under pdbqt format.

Ligands Preparation and Filtration

To search for similar INH-NAD molecules, a search by similarity was carried out using PubChem, 1,000 molecules were selected with a Tanimoto threshold 90%, and then downloaded under sdf format with 2D conformation, the conversion to 3D, the elimination of duplicate and then the conversion into a pdb file format was carried out by graphical interface of OpenBabel.²⁸

The toxicity properties of compounds were evaluated using two software products, STopTox^29,30 and ProTox-II.³¹ In this step, 892 molecules with a toxicity class-4/5 with a predicted LD50 (mg/kg) ranged from 2,000 to 5,000 for Class V and more than 5,000 for class VI by ProTox-II and noncardiotox by STopTox were kept.

Drug likeness score predictions of the best compounds were evaluated using DruLiTo.³² The 823 compounds were selected for molecular docking preparation, all the ligands were prepared and contain atom types supported by AutoDock tools plus extra records that specify rotatable bonds and then saved in pdbqt format.

Mutation Effect Prediction

To study the impact of mutation on protein function, the site directed mutator (SDM)³³ and DynaMut³⁴ online were used to predict the mutations impact on protein stability and flexibility respectively.

Calculation of RMSD of Re-docking Test

As a positive control, and to perform re-docking calculation from Protein Data Bank, the X-ray crystallographic complex protein-ligand (PDB ID: 1ZID). The re-docking test was evaluated using BALLView software.³⁵

Docking and Scoring

After selecting the active site residues of the target protein using AutoDock Tools V.1.5.7,³⁶ allowed the preparation of the grid box, with a spacing of 1 Å for each macromolecular protein and center coordinates fixed at X= −8.411 Å, Y=38.576 Å, Z=11.141 Å and a size of X=20 Å, Y= 20 Å, Z=20 Å. The grid settings file was saved from the grid menu output option.

In order to choose ligands with lower affinity scores, Autodock Vina was used to assess the binding affinity and bound conformation of those selected ligands. This tool is designed for protein-ligand docking, using multiple CPUs at a time, making it faster and more accurate and uses the Lamarckian genetic algorithm and semi-empirical free energy force field that generates free binding energy of 10 ligand poses after docking.³⁷ The 3D protein-ligand interactions were visualized using PyMOL V.2.5.2.³⁸

3D-pharmacophore Model Building

To build a pharmacophore profile of the three lead compounds, pharmacophore query methodology was used, this method is based on searching 3D distances between features like HBD/HBA hydrogen bond donors/acceptors, PI positive ionizable, ARO aromatic ring, Hyd hydrophobic centres.³⁹ For this purpose, the Molecular Operating Environment (MOE) tool was used. An MOE database was created using the builder option, the energy was minimized and then the alignment was launched to generate the best pharmacophoric model with a pharmacophoric map with the first common features.

Results

Mutation Effect Prediction

In order to start the second molecular docking, a detailed study of the mutational profile of InhA protein was carried out.

Fifteen mutations were collected in the literature and the first filtration was carried out by PROVEAN by which it was possible to identify the neutral mutations from those deleterious (Figure 2) and also their effect on the solvent accessibility of the protein (Figure 3A). Subsequently the ten putative deleterious mutations were chosen to study their effect on the flexibility and stability of the protein using SDM and DynaMut (Figure 3B).

Figure 2 The distribution of the 15 mutations according to the function, deleterious or neutral.

Figure 3 (A) The effect of the deleterious mutations on solvent accessibility surface area (SASA) between the wild-type protein (in blue) and the mutanttype protein (in orange). (B) The effect of the deleterious mutations on the flexibility (in orange) and stability (in blue) of the InhA protein.

Re-docking Score

The re-docking test was evaluated by calculating the RMSD between the predicted pose and the co-crystallographic ligand. The score of RMSD between the two complexes show a value lower than 2 Å with a value of 0.38 Å.

Docking and Scoring

To filter out the best compounds from molecular docking, 823 compounds were docked with the wild-type InhA protein, and only 34 compounds were chosen with a binding energy to InhA protein less than −10.7 kcal/mol, since the reference ligands INH-NADH binds to InhA protein with an affinity of −10.7 kcal/mol, and NADH alone binds with an affinity of −7.5 kcal/mol. The 34 compounds selected from the first docking that showed better affinity than the reference were docked again separately with the 10 mutated model generated.

Three ligands were preferred over all compounds with their lowest binding affinity score toward the wild-type InhA protein and the 10 mutant models. Comparative docking results are depicted in Figure 4.

Figure 4 The three best leads with lowest affinity binding to the wild-type protein and the the mutant models.

Table 1 shows the reference ligand INH-NAD, NADH and the three molecule candidates after the second docking against the mutant models with lowest binding affinity than INH-NAD ligand with their essential information: molecular formula, binding affinity energy, chemical structure, PubChem CID, molecular weight, toxicity test, and finally predicted fragment contribution in toxicity test.

Table 1 The essential information about the reference ligand INH-NAD, natural ligand NADH and the Three best ligands including molecular formula, binding affinity score, hydrogen bond CID in PubChem, molecular weight, and toxicity test

Figure 5 illustrates the 3D protein-ligand interactions visualized using PyMol. Figure 5A and B shows the interaction between the reference ligand INH-NAD and NADH alone against the InhA protein respectively. While Figure 5C–E depicts interactions between ligands (CID: 25176455, 91754233 and 101256471) and InhA protein accordingly.

Figure 5 Visualization of 3D interactions against InhA protein. (A) Interaction between the reference ligand Inh-NAD and the inhA protein. (B) Interactions between NADH and InhA protein. (C) Interaction between ligand (CID:25176455) and inhA protein. (D) Interaction between ligand (CID:91754233) and InhA protein. (E) Interaction between ligand (CID:101256471) and InhA protein.

3D-pharmacophore Modelling

Using the Molecular Operating Environment (MOE) tool, a pharmacophoric map of the best three compounds was generated.

Common features with inter-features distances of these compounds are shown in Figure S2.

Figure 6 shows a simplified view of the pharmacophoric map with the aligned ligands (Figure 6A) and without ligands (Figure 6B). The 3D-pharmacophoric model consists of a four hydrophobic (Hyd) center localized in the middle of the structure connected on both side to two Acc and ML features with a distance of 8.71Å between F11 and F1 and 5.35 Å between F2 and F8. The structure is also characterized by an aromatic centre (Aro) F6 leading to another aromatic/acceptor/metal ligator F16 (ML/Aro/Acc). The aromatic regions play an important role in binding to the enzyme and the HBA/HBD group help the ligand to bind tightly to the target.

Figure 6 The pharmacophore model generated after the flexible alignment of the best three compounds. (A) shows the different features explored in the model structure with the aligned ligands. (B) shows a simplified view of the model without ligands.

Discussion

In this study, computer prediction methods have been adopted to explore new inhibitors of Mycobacterium tuberculosis directly targeting the InhA protein in order to overcome the problem of resistance to treatments, linked to mutations in this bacterium. First, a prediction of the impact of different single mutations on the properties of the InhA protein was studied in order to reveal the deleterious mutations, then, a collection of 1,000 inhibitors similar to the INH-NAD complex and a study of toxicity were carried out. The interaction of the 823 nontoxic ligands with the wild-type InhA target was analyzed by docking in order to select those which present a better affinity than the INH-NAD complex and which are of the order of 34, including the interaction with the models of the InhA protein mutated by deleterious mutations was studied by a second docking. This second docking highlighted three compounds with a better affinity than that of the INH-NAD complex which were used for the generation of the pharmacophore model.

After having extracted all the mutations of the InhA protein described in the literature, 10 mutations among 15 were identified as deleterious. Analysis of the effect of each deleterious mutation on the protein’s flexibility, stability, and surface area solvent accessibility (SASA) revealed that all the mutations described as deleterious D148G, I16T, I21V, I47T, I95P, I194T, K8N, K165A/Q, S94A, T266A/D, V78A and Y158A/S would have a direct effect on function by altering the various properties studied, thus modifying the folding of the protein, the protein-protein or protein-ligand interaction, confirmed by several studies.^17,40,41 This is confirmed by our results, which show a direct relationship between SASA and stability, especially in these mutations I95P, I16T, D148G and I194T, this relationship was also mentioned in a study that stability is the fundamental property improving biomolecular function, and regulation.⁴²

Several studies have indicated that many diseases, including drug resistance-related mutations, are located in solvent-accessible sites, suggesting that analysis of these mutations could also shed light on the mechanisms underlying the diseases.^43,44

Apparently, one can clearly see that the change in amino family makes the protein more accessible to solvent, as shown by the I16T and I194T mutations, a switch from isoleucine to tyrosine, a nonpolar aliphatic amino acid to an amino acid polar and uncharged; a clear change of family, polarity and charge loss causes an increase in solvent accessibility with a difference of +16.7% (from 54.3% to 71.0%) and 10.3% more (from 5.3% to 15.6%), respectively. I194T also causes a decrease in stability at −2.51 kcal/mol.

Regarding the I95P mutation, there was no change in the family of amino acids, although our work shows that this mutation would cause a decrease in stability at −3.45 kcal/mol.K and that it would make the protein more accessible to solvent with a difference of 7.2% (from 20.4% to 27.6%) as clearly shown in (Figure 3A). This is confirmed by a study by Kumar and Sobhia as a mutation influencing the binding energy by decreasing it by −2.01 kcal/mol and believed to disrupt the INH-NADH bond due to the change in the nature of the interacting residues, and strongly reducing the direct and water-mediated interactions of the H bonds between NADH and binding site residues allowing conformational flexibility to NADH, particularly at the pyrophosphate region, and causing a weakening of its binding to the dinucleotide binding site.⁴⁵

For the D148G mutation, a clear difference was noticed in SASA between the wild-type and the mutant-type protein, a change from aspartic acid to glycine amino acid; a polar and negatively charged amino acid to a polar and uncharged amino acid, a change of family and a loss of charge made the protein more accessible to the solvent from 0.1% to 11.5% with a difference of 11.4%. This mutation is described in a study by Hartkoorn et al as indeed responsible for pyridomycin and isoniazid resistance.²²

Changes in amino acid size, amino acid charge and hydrogen bond number and changes in SASA have been shown to cause the loss of stability as well as aberrant folding and aggregation of the protein.⁴²

To explore effective inhibitors to each of the deleterious mutations studied that would prove to be involved in resistance to treatment, we analyzed the interaction of the 34 compounds by docking for the models of these mutations, which showed a better affinity than the INH-NAD complex to wild-type InhA.

Three leads showed a lower binding affinity better than the reference ligand INH-NAD using AutoDock Vina with the following CID 101256471, 91754233, and 25176455, the first one (C25H30N7O16P2+) binds with a binding affinity of −12.2 kcal/mol with the wild-type InhA protein with a number of eight hydrogen bonds; one bond of 2.7 Å with the ALA-191, five bonds with three water molecules with a distance ranged from 2.3 Å to 3.5 Å respectively, and two bonds of 1.9 Å and 2.6 Å with GLY (14–192) respectively, the second compound (C27H34N8O15P2) binds with the wild-type protein with an energy of −12.0 kcal/mol, and a number of eight hydrogen bonds; four hydrogen bonds with four water molecules with a distance ranged from 2 Å to 3.5 Å, and three bonds of 2.5 Å, 2.6 Å and 2.2 Å with ILE-21, LEU-63 and LYS-165, respectively. The last compound (C23H28N7O14P2+) binds with the wild-type protein with an energy of and −13 kcal/mol and a number 10 hydrogen bond; with three water with molecules with four hydrogen bonds ranged from 2.1 Å to 3.3 Å, two bonds of 2.0 Å and 2.9 Å with GLY-192 and two bonds of 2.1 Å and 2.1 Å with GLY (14–96), respectively, one bond of 2.5 Å with TYR-158 and one bond of 2.7 Å with SER-20 as demonstrated in Figures 4 and 5.

Based on the structure of these three compounds, the pharmacophoric map showed that the hydrophobic Hyd groups and the aromatic Aro regions play a vital role in their binding to the hA protein. HBA groups enhance the interaction between ligands and the active site of the target^46,47 and should be considered as crucial structural features.⁴⁸ We hope that our contribution can offer a suitable pharmacophore for further studies.

Conclusion

The number of people dying yearly because of tuberculosis is growing rapidly due to the emergence of multidrug-resistant tuberculosis, extensively drug-resistant tuberculosis and totally drug-resistant tuberculosis, which makes the discovery of a new drug an emergency to confront this infectious disease. InhA protein is one of the most effective targets for eradicating M. tuberculosis as long as this protein is a key player in the biosynthesis of mycolic acid cell wall.

In this sense we were able to propose in silico, through this study, three new direct inhibitors of InhA, for both for wild-type and resistance causing mutant-type, nontoxic showing better results than INH-NAD adduct, from which a pharmacophore map was generated to fdetermine common features between the best ligands found. This model can open up a new template for further studies to discover novel powerful anti-tuberculosis drugs.