Computational approaches to identify novel inhibitors for the drug‐ resistant Mycobacterium tuberculosis DprE1 enzyme

Mycobacterium tuberculosis causes tuberculosis (TB), which is a common but life‐debilitating disease. The continued development of resistance to frontline anti‐TB drugs such as isoniazid and rifampicin threatens the efficacy of currently available treatment procedures. This highlights the need to explore diverse approaches essential for drug development against multi‐drug‐resistant strains of tuberculosis. Drug development relies on the findings associated with novel protein targets, which play a crucial role in the disease life cycle. DprE1, an enzyme that plays a critical role in the cell wall synthesis of M. tuberculosis , has been recognized as a promising target for drug development. In the present study, based on previous experimental findings, seven mutant models of DprE1 involved in DprE1 resistance are predicted using homology modeling. Further, potential inhibitors are selected based on their efficacy and IC 50 values. Shortlisted inhibitors are docked with the wild‐type and mutant structures of DprE1. The deduced inhibitor molecule (ZINC5) is found to possess high potential as a lead inhibitor for all the models of DprE1. It can be used to circumvent drug resistance in the current treatment regime.


Introduction
The bacillus Mycobacterium tuberculosis affects millions of people worldwide, mainly by causing tuberculosis (TB).Among the deaths that occurred due to infectious diseases, TB was one of the leading causes (Kim et al. 2020).With advancements in science and technology, the number of TB cases worldwide has been decreasing in the last two decades (Kyu et al. 2018).However, in 2018, the World Health Organization estimated 10 million new cases of TB and 1.45 million deaths (Liu et al. 2020), with a re ported number of 1,300,000 TB deaths in HIVnegative patients and 300,000 TB deaths in HIVpositive patients (Hariguchi et al. 2020).In spite of advancements in the treatment of the disease, there is persistence in the num ber of TB cases worldwide.It is a result of the signif icant emergence of multidrugresistant (MDR) and ex tensively drugresistant (XDR) strains of M. tuberculosis (Wilsey et al. 2013).Novel target identification and novel drug development are the two potential approaches that scientists are focusing on to achieve effective treatment for the disease.The novel drug development approach is used to identify suitable agents against a known target with the help of computeraided drug discovery (CADD).This pathway includes methods like quantitative structure activity relationship (QSAR), virtual screening, molecular docking, molecular dynamics simulation, etc.In the past several years, multiple druggable targets of M. tuberculo sis have been identified and tested for clinical use (Maharaj et al. 2015).
Indonesian Journal of Biotechnology 28(3), 2023, 180-190 followed by reductase DprE2 catalyse the epimerization of decaprenylphosphorylβᴅribofuranose (DPR) to DPA.Inhibition of DprE1 halts the synthesis of DPA, which subsequently reduces the formation of arabinan, an important cell wall component (Oh et al. 2021).There fore, M. tuberculosis cell wall synthesis is prevented due to the depletion of DPA, a necessary precursor for the synthesis of arabinan (Maharaj et al. 2015; Wilsey et al. 2013).Hence, DprE1 is one of the most studied protein targets for drug development against M. tuberculosis.However, the polymorphism and mutation in DprE1 confound the current treatment with the problem of drug resistance (Neres et al. 2015; Foo et al. 2016).Thus, it is essential to develop a novel DprE1 inhibitor for its potential use against MDR and XDR M. tuberculosis.
In the present study, we employed in silico approaches to discover novel inhibitors for the M. tuberculosis target protein DprE1, incorporating recently identified mutant residues.Subsequently, virtual screening has been per formed to identify potential compounds with good bind ing affinity against DprE1 target.The identified novel in hibitor molecule will be acting as potential inhibitor and overcome the drug resistance in M. tuberculosis.

Homology modelling of wild-type and mutant DprE1
The protein sequence of DprE1 (UniProt ID P9WJF0) was retrieved from the UniProt protein sequence database (https://www.uniprot.org,The UniProt Consortium (2021)).There are two reviewed entries composed of 461 AA residues.Further, both sequences are 100% identical.Therefore, we selected one entry, P9WJF0, for our study.
The template search was performed using BLASTp against the PDB database (https://www.rcsb.org).There are a total of 27 structures of DprE1 available in the PDB with different resolutions.However, all crystal structures are not full length and have unmodelled regions.From the obtained results, the structure of 4P8C.pdb(Neres et al. 2015) was selected for further studies as it shows 100% identity with the query sequence.Anyway, this structure lacks coordinates for amino acid regions 1-7, 268-274, and 323-326.The initial regions 1-6 are not involved in the interaction with the ligand.However, the 268-274 and 323-326 regions near the catalytic groove of DprE1 might be involved in the interaction with the ligand molecule (Neres et al. 2015).Therefore, these regions were modelled by the loop modelling option of Discovery Studio version 3.5 (Barage andSonawane 2014; Barage et al. 2017).Finally, 20 homology models of DprE1 were generated using the MODELLER 10.2 program based on the obtained template structure.Among them, a single model was selected on the basis of the discrete optimized protein energy (DOPE) score (Webb and Sali 2021).
The structural validation and characterization of the predicted models were performed using the PROCHECK and ERRAT tools available on the SAVES v6.0 server (https://saves.mbi.ucla.edu).Further, the ProSAweb server was also used to measure the protein quality (Laskowski et al. 1993; Wiederstein andSippl 2007).Subsequently, the selected model was used to generate the mutant models.Based on previous experimental studies (Makarov et al. 2009; Neres et al. 2015; Foo et al. 2016) (Makarov et al. 2009; Neres et al. 2015; Foo et al. 2016).All of the gener ated mutant and wildtype models were subjected to energy minimization using the CHARMM force field (http://mackerell.umaryland.edu/charmm_ff.shtml).Subsequently, all models were used for docking studies.

Preparation of ligand dataset
All known DprE1 inhibitors were retrieved from previous experimental reports.The inhibitors PBTZ169, BTZ043, Azaindole and BTO were shortlisted based on their lowest IC 50 values (Table 1).These shortlisted inhibitors were used for similarity searches against the ZINC database by setting the Tanimoto index to 90 (Sterling and Irwin 2015).From the obtained ligands, only purchasable lig ands were extracted in (.sdf) format.In total, 46 ligand data sets were prepared, which include 42 similar and 4 known inhibitors, as listed in Table 2.All retrieved ligand molecules were loaded in PyRx through Open Babel and minimized using the MMFF94 force field (Dallakyan and Olson 2015).The minimized compounds were converted to an autodock compatible format (.pdbqt).The protona tion state and Gasteiger charges were assigned to all ligand atoms.

Target preparation and virtual screening
Molecular docking studies were carried out using AutoDock 4.2, implemented in PyRx (Dallakyan and Ol son 2015).The docking protocol has been validated using a redocking experiment in which the crystallographic pose of Y22 was separated from DprE1 (4P8C.pdb)and redocked into the active site of DprE1.The predicted binding mode of Y22 was compared with the crystallo graphic pose of Y22 using RMSD calculations (Barage et al. 2017; Meshram et al. 2020).The present docking protocol successfully reproduced the crystallographic position, and the same has been used for the unknown ligand molecules.The model of DprE1WT was loaded in PyRx, and the charges and protonation states were  assigned to the protein atoms.It has been reported that Cys387, Val365, Tyr60, His132, Lys418, Gln334 and Asn385 residues are involved in catalysis and the interaction of their substrates (Bhutani et al. 2015), hence they are considered flexible residues for docking.The rest of the residues were treated as rigid and converted to autodock macromolecules.Grid preparations were performed using the AutoGrid tool of PyRx, where the grid box was designed in such a way that it could encompass all the flexible residues, with dimensions of 46 × 49 × 50.The grid spacing was set to 0.375 Å.
In this study, we performed 100 docking runs using the Lamarckian Genetic Algorithm (LGA) for each ligand.For Lamarckian GA, the parameters were set as 10 runs of each GA, the number of individuals in the population was set to 150, the maximum number of generations was set to 27,000, the rate of mutation was set to 0.02, the rate of crossover was set to 0.8, and other parameters were used as default values.Based on an Autodock scoring function, which was used to score the docked poses, the best confor mation of each ligand was given the highest rank that also showed the least binding free energy.Adequate spaces of 0.25 Å for translation and 5 Å for rotation were chosen with a maximum energy evaluation of 2.5 × 106.All other docking parameters were set to their default values.Simi larly, the abovementioned docking protocol was used for all the mutant models considered in this study.Confor mational clustering was performed on docked conformers with an RMSD tolerance of 1 Å.The top three most pop ulated clusters were selected from each docking round of the mutant.The single lowest energy conformation was selected from the highest populated cluster and used for further interaction analysis.Docking studies were done to focus on hydrogen bond interactions formed between the ligand and surrounding active site residues of DpRE1, along with van der Waals, hydrophobic and other interactions that may contribute to the stability of the protein-ligand complexes.

Structural modelling of DprE1 wild-type and its mutants
We have predicted the fulllength structure of DprE1 with missing residues using homology modelling with the MODELLER 10.2 program, as described in Section 2.1.The stereochemical quality of the resultant structure was analysed using various tools on the SAVES and ProSA web server.
The Ramachandran plot of the model showed 94% residues in the favoured region and 4% residues in the al lowed region, as shown in Figure 1c.ERRAT used to de termine the quality of the model (Figure 1a) showed an overall quality factor of 91.49, which signifies an excel lent quality of the predicted model (Colovos and Yeates 1993).The ProSAWeb Zscore of the model is 10.9 (Fig ure 1c), which indicates good model quality (Wiederstein and Sippl 2007).Also, the predicted knowledgebased en ergy used to evaluate the local model quality was found to be in the acceptable range.Overall, with a good qual ity score, the predicted homology model predominantly consists of sheets, helices, and coils with wellvalidated geometry.The superimposed structure of the predicted DprE1 model with the template model is shown in Fig ure 2, and the RMSD Cα value is 0.303 Å.This predicted homology model was used for the generation of mutants, which was followed by docking.Information related to the antibiotic resistant mutants of DprE1 was retrieved from the literature (Makarov et al. 2009; Neres et al. 2015; Foo et al. 2016).The seven mutant residues discussed in Sec tion 2.1 (MutC387A, MutC387S, MutC387T, MutC387G, MutC387N, MutG17C and MutL368P) were individually incorporated in the predicted model at corresponding po sitions using Discovery Studio 3.5.Furthermore, a sin gle model was prepared by incorporating all the mutant residues (MutAll).All the mutant models, along with DprE1WT and MutAll were used for the docking stud ies.

Virtual screening of compounds against predicted DprE1 models
To characterise the binding mode of a known inhibitor, initially, the Y22 molecule cocrystalized with DprE1 (4P8C.pdb)was separated and redocked to check whether Autodock is able to reproduce the binding mode.The obtained RMSD value, 0.72 Å indicates the good relia bility and consistency of the docking protocol (Maharaj et al. 2015; Barage et al. 2017; Meshram et al. 2020).Further, we analysed the molecular interactions of Y22 with DprE1, as depicted in Supplementary Material Figure S1.Also, the binding mode of the identified novel ligand molecules needs to be evaluated in comparison with the known ligands.Thus, to characterize the binding mode of the identified ligand molecules, molecular docking was used as described in Section 2.3.In total, 46 ligands were used for docking with the predicted DprE1WT model and other mutant structures.The docked conformation for each ligand was subjected to conformational cluster ing with a RMSD tolerance of 1 Å.

Selection of top leads as DprE1 inhibitors
We docked each of the 46 ligand molecules with DprE1 WT, MutAll, and the seven mutant models of DprE1.
The number of highest populated clusters (HPC) and their  corresponding binding energies (BE) were used for short listing the leads.Conformational clustering of 46 ligand molecules docked with the DprE1WT structure showed three highest populated clusters with 20, 16 and 15 con formations for the ligands ZINC5, ZINC8 and PBTZ169, respectively (See Supplementary Material Table S1 for ZINC IDs).The ligands ZINC5, ZINC8 and PBTZ169 corresponding to the highest populated clusters showed the lowest binding energies of -10.38, -10.25 and -9.08 kcal/mol, respectively.The molecular interactions of these ligand molecules with respect to the receptor were analysed.Similarly, the docked conformations of Mut All and all other mutant structures were clustered based on an RMSD tolerance of 1 Å.From each docked com plex, three lowest energy conformations were extracted from the highest populated clusters.Several ligand con formations were observed for their number of highest pop ulated clusters and binding energies, as listed in Table 2. Furthermore, from all ligand molecules docked with all DprE1 structures, the top three docked conformations were analysed and evaluated for their binding mode in the catalytic groove of DprE1.Supplementary Material Table S1 enlists the ZINC IDs and their short forms as used in this study.We identified all the ZINC ligand molecules based on the similarity search with known in hibitors.So, the ZINC5 ligand has substructure similar ity with the known inhibitor Benzothiazinone compound, which is used as an antitubercular drug (Karoli et al. 2012).Interestingly, among the top ligand molecules bound with   the nine DprE1 models, ZINC5 was found to be common in eight complexes, while PBTZ169 was common in four.
Conversely, all other top three ligands were not observed in common with any docked complexes.Although the lig ands ZINC5 and PBTZ169 were found to be common in multiple docked conformations, the top three ligands of all docked complexes were evaluated for their molecular interactions and the binding mode in the catalytic groove of DprE1 (Supplementary Material Figure S1).Our quest was to understand the substructure specificities of these ligands with the active site of DprE1 and to study the im portant residues involved in interactions with the ligand molecules.
The docked conformations of the ligands with the Gly17Cys17 mutated model were subjected to confor mational clustering analysis with a RMSD cutoff of 1 Å.The top 3 ligands, ZINC10, ZINC39 and ZINC5, showed the highest populated clusters with 35, 35 and 34 conformations, respectively.Similarly, the top 3 lig ands for the Leu368Pro368 mutated model were ZINC5, ZINC28 and PBTZ169, which showed the highest popu lated clusters with 58, 53 and 51 conformations, respec tively.Since ZINC5 showed both better binding energy and a populated number of clusters, it was considered to be the best ligand.Similarly, for the Cys387Ala387 mutated model, ZINC5, ZINC28 and ZINC42 were ob served to have the highest populated clusters with 55, 44 and 42 conformations, respectively, where ZINC5 exhib ited the best binding energy of -11.86 kcal/mol.Addi tionally, ZINC36, ZINC5 and PBTZ169 had the highest populated clusters of 31, 28 and 27 conformations, re spectively, for the Cys387Ser387 mutated model; for this model too, ZINC5 was selected owing to its better bind ing energy of -10.69 kcal/mol.For the Cys387Thr387 model, BTZ043, PBTZ169 and ZINC5 showed the high est populated clusters with 64, 56 and 43 conformations, respectively.The binding energies of PBTZ169 (-10.36 kcal/mol) and ZINC5 (-10.12 kcal/mol) were found to be better as compared to the BTZ043 ligand.In the case of the Cys387Gly387 model, ZINC5 showed the highest populated clusters and lowest binding energy of 54 and -11.12 kcal/mol, respectively, in comparison to ZINC28 and ZINC17.With respect to the Cys387Asn387 mutated model, ZINC36, ZINC10 and BTZ043 exhibit the highest populated clusters of 48, 48 and 42 conformations, respec tively.When all three residues, Gly17Cys, Leu368Pro, and Cys387Ala, were mutated (MutAll), ZINC5 showed the lowest binding energy of -11.63 kcal/mol with the highest populated cluster of 35 conformations.Other top ligands were ZINC7 and ZINC17, with binding energies of -9.09 and -8.69 kcal/mol, respectively.

Interaction analysis of topmost leads
Molecular interactions of the top 3 ligand molecules were analysed to understand the substructure specificity and complementarity with the active site of MutAll DprE1.All the bound conformations of ZINC5 were observed in the active site groove of the MutAll DprE1 model, as de picted in Figure 3.In this section, we focus on the topmost lead identified in each protein model.The detailed inter action analysis of all protein models is provided in Sup plementary Material.
WildType: The top ligand molecule ZINC5 inter action with wildtype analysis revealed hydrogen bond ing interactions with Thr118, His132 and Gln336 residues with the oxygen atoms.The residues Leu317, Leu363, Val365, Ile131, Cys129 and Ala417 show alkyl and pi alkyl interactions with the ZINC5 molecule (Figure 4).Tyr60 forms pistacking interactions with the aromatic ring of ZINC5.Overall, an excellent interaction network of DprE1 residues has been observed with ZINC5 (Figure 5).
Cys387Ala387: ZINC5 exhibits the best interactions with the C387A mutant model.Tyr415 shows a pisulfur interaction with the sulfur atom and also shows pipi T shaped interactions with the aromatic rings of ZINC5 (Fig ure 4).Also, the His132 residue exhibits Hbonding inter actions with the oxygen atom of ZINC5.Ala417, Cys129, Val121 and Arg58 form alkyl and pialkyl interactions with ZINC5 (Figure 5).
Cys387Gly387: The interactions of ZINC5 with mu tated model C387G have the best HPC and binding en ergy.The fluorine atom of ZINC5 makes halogen interac tions with the residues Asn385, Ile386, Phe366 and Val365 (Figure 4).The residues Ile131, Pro116, Ala417, Lys367 and Arg58 also show alkyl and pialkyl interactions with the ZINC5 ligand molecule.Also, Val121 shows pisigma interactions with the aromatic ring of the ZINC5 ligand (Figure 5).
Cys387Asn387: The topmost ligand interacting with the C387N mutant model is BTZ043.BTZ043 shows alkyl and pialkyl interactions with Pro116, Ile131, Ala417, Tyr415 and Lys367 residues (Figure 4 and Supplemen tary Material Figure S4).Val365, Asn385 and Phe366 show halogen interactions with the fluorine atoms of the BTZ043 ligand molecule.Cys129, Lys418 and His132 make hydrogen bonding interactions with the oxygen atom of BTZ043 (Figure 5).
Cys387Ser387: ZINC5 is the top ligand against the MutC387S model.Tyr415 shows a pisulfur interaction with the sulfur atom ZINC5.The fluorine atom exhibits Hbonding interactions with the Ser387 residue of ZINC5 (Figure 4).The fluorine atom of ZINC5 also makes halo gen interactions with residues Asn385 and Val365.The residues Val365, His132, Ala417, Ile131, Pro116, and Cys129 interact by alkyl and pialkyl interactions with ZINC5.Also, Val121 and Arg58 show alkyl and pialkyl interactions, and Cys129 shows a pisulfur interaction with the aromatic ring of ZINC5 (Figure 5).
Leu368Pro368: The fluorine atom of ZINC5 makes Halogen interactions with the Gln336 residue.The oxygen atom of ZINC5 makes Hbonding interactions with His132 residues (Figure 4).The residues Ile131, Ala417, Arg58, Val121 and Cys129 show alkyl and pialkyl interactions with the ZINC5 ligand molecule.Also, Val121 and Arg58 show alkyl and pialkyl interactions, and Cys129 shows pisulfur interactions with the aromatic rings of ZINC5.Try415 residue shows a pipi Tshaped interaction with ZINC5 (Figure 5).
Gly17Cys17: The top ligand, ZINC5, shows a pi pi Tshaped interaction with Tyr415 of the MutG17C model.His132 residue makes Hbonding interactions with the NH 2 group of ZINC5.The oxygen atom of ZINC5 makes Hbonding interactions with the Lys418 residue.The fluorine atom of ZINC5 makes Halogen interactions with residues Asn385 and Val365.The residues Phe369, Lys367, Ile131, Pro116 and Ala417 show alkyl and pi alkyl interactions with ZINC5 (Figure 4).Also, Val121 shows pisigma interaction, Arg58 shows alkyl and pi alkyl interaction, and Cys129 shows pisulfur interaction with the aromatic rings of ZINC5 (Figure 5).
MutAll (Gly17Cys, Leu368Pro, Cys387Ala): The top ligand complexed with the MutAll model is ZINC5.Val121, Arg58, Ala417 and Pro116 react with the carbon rings of ZINC5 via alkyl and pialkyl interactions.Cys129 shows pisulfur and pialkyl interactions with ZINC5 (Fig ure 4).Along with alkyl and pialkyl reactions and pi sigma reactions with ZINC5, the Tyr415 residue also ex hibits pisulfur interactions with the sulfur atom.His132 and Lys418 represent hydrogen bonding interactions with the oxygen atom of ZINC5.Also, His132 makes halogen interactions with the fluorine atom and forms alkyl and pi alkyl interactions with the carbon atoms of ZINC5.Ad ditionally, Gln336 and Asn385 residues exhibit halogen interactions with the fluorine atoms of the ZINC5 ligand (Figure 5).
Overall, molecular interactions were analysed for the top three ligands with the wildtype and all mutant struc tures of DprE1.We show that ZINC5 has better inter action chemistry than Y22 with wild and mutant DprE1 (Supplementary Material Figure S1 and Figure 4).Inter estingly, the ZINC5 molecule was observed to efficiently accommodate in the active site of wildtype as well as all mutant structures of DprE1 (Figure 3).Furthermore, the interactions of DprE1 residues with the ZINC5 ligand in all the docked complexes were found to be conserved.In addition, the binding free energy of the ZINC5 ligand was found to be the lowest in most of the docked com plexes.Whereas, compounds ZINC28 and ZINC17 has been found to be potential inhibitor for some mutant mod els.

Conclusions
In the present study, 3dimensional structures of the DprE1 molecule were predicted using homology modelling and clinically relevant mutations were incorporated into the predicted structure at corresponding positions.The mutant models of DprE1 were subjected to energy minimization and were further used for docking studies.It was observed that the mutant structures did not exhibit a noticeable en ergy difference in comparison to the wildtype structure of DprE1, which indicates that the mutations did not af fect the overall conformation and activity of the enzyme.The identified novel ligand molecules, along with known inhibitors, were screened against the wildtype and each mutant structure of DprE1.The results obtained from the docking studies revealed that the ZINC5 molecule pos sessed a good binding affinity towards the wildtype and all mutant structures of DprE1, except for MutC387N.The ZINC5 molecule was observed to have efficiently accom modated itself in the catalytic groove of DprE1.It depicted conserved hydrophobic and hydrophilic interactions with the key residues of DprE1 in all the docked conforma tions.Overall, the ZINC5 molecule can function as a po tential inhibitor candidate against DprE1 inhibitor, which can overcome the resistance of M. tuberculosis towards contemporary antibiotic strategies.Further, experimen tal validation essential to test the efficacy of compound against DprE1 for its potential use as lead molecule.

FIGURE 3
FIGURE 3 Superimposition of all docked ZINC5 conformations with Mut-All protein model.

TABLE 1
Highest populated clusters (HPC) and binding energies (BE) of the top ligands in nine models of DprE1.

TABLE 2
Known inhibitor structures of DprE1 used for similarity search against ZINC database.

TABLE 3
Receptor-ligand interactions of the multiple drugs against the modelled DpRE1 targets.