Aqueous Molecular Dynamics Simulations of the M. tuberculosis Enoyl-ACP Reductase-NADH System and Its Complex with a Substrate Mimic or Diphenyl Ethers Inhibitors

Molecular dynamics (MD) simulations of 12 aqueous systems of the NADH-dependent enoyl-ACP reductase from Mycobacterium tuberculosis (InhA) were carried out for up to 20–40 ns using the GROMACS 4.5 package. Simulations of the holoenzyme, holoenzyme-substrate, and 10 holoenzyme-inhibitor complexes were conducted in order to gain more insight about the secondary structure motifs of the InhA substrate-binding pocket. We monitored the lifetime of the main intermolecular interactions: hydrogen bonds and hydrophobic contacts. Our MD simulations demonstrate the importance of evaluating the conformational changes that occur close to the active site of the enzyme-cofactor complex before and after binding of the ligand and the influence of the water molecules. Moreover, the protein-inhibitor total steric (ELJ) and electrostatic (EC) interaction energies, related to Gly96 and Tyr158, are able to explain 80% of the biological response variance according to the best linear equation, pKi = 7.772 − 0.1885 × Gly96 + 0.0517 × Tyr158 (R2 = 0.80; n = 10), where interactions with Gly96, mainly electrostatic, increase the biological response, while those with Tyr158 decrease. These results will help to understand the structure-activity relationships and to design new and more potent anti-TB drugs.


Introduction
Tuberculosis (TB) is an infecto-contagious disease caused by Mycobacterium tuberculosis (MTB), mainly affecting lungs, but it can also infect others vital organs, such as central nervous, genitourinary, and osteoarticular systems [1][2][3]. In 2013, according to the World Health Organization (WHO), about 1.5 million people die every year from TB, and this disease is the second cause of death worldwide [4]. Since the 1990s, the WHO recommends the DOTS (directly observed treatment, short-course) strategy that includes a chemotherapy regimen combining four first-line drugs (i.e., isoniazid, rifampicin, ethambutol, and pyrazinamide) for at least six to nine months [5,6]. Despite this scheme having high rates of a successful cure, many patients discontinue treatment due to serious side effects, such as hepatotoxicity, which can be fatal [4,7].
The main consequence of the TB treatment abandonment is the emergence of multi-drug resistant (MDR) strains, showing resistance to at least isoniazid and rifampicin, and extensively-drug resistant (XDR) strains, with the same profile of MDR strains plus resistance to a fluoroquinolone (e.g., ofloxacin or ciprofloxacin) and at least one injectable drug (e.g., amikacin, kanamycin or capreomycin) [8,9]. As a result, chemotherapy required against MDR-TB and XDR-TB has several disadvantages, e.g., it is less effective, more expensive, and more toxic in comparison to the first-line drugs [9].
Thus, due to the prolonged treatment, toxicity, and emergence of resistant strains to the anti-TB drugs, it is necessary to develop new, more potent and effective drugs. Therefore, the enoyl-acyl carrier protein (ACP) reductase (ENR) is an attractive target to design new drugs, because it is an essential enzyme of the type II fatty acid synthesis (FAS-II) [10]. The ENR from M. tuberculosis (InhA) catalyzes the reduction of trans-2-enoyl-ACPs that is dependent of the NADH cofactor ( Figure 1). In fact, the anti-tuberculosis drug isoniazid acts as a pro-drug, activated by oxidation catalyzed by M. tuberculosis catalase-peroxidase (KatG). Hence, this product and cofactor (NADH or NAD+) react to form an adduct that inhibits InhA, disrupting the biosynthesis of mycolic acids (FAS-II), the main components of the mycobacterial cell wall, thus causing cell death [11]. The required activation of isoniazid is responsible for the emergence of isoniazid-resistant strains related to KatG mutations [12][13][14]. This fact led to efforts to identify direct InhA inhibitors that do not require activation through KatG [15][16][17][18]. In this context, diphenyl ethers are a promising class because they can directly inhibit InhA. Triclosan (5-chloro-2-(2,4-dichlorophenoxy)phenol, TCL (Figure 1), an important representative of this class, is an antimicrobial (antibacterial and antifungal) agent found in toothpaste and deodorants, which inhibits ENRs in pathogenic organisms, such as Escherichia coli and Staphylococcus aureus [19][20][21].
Recently, 3D structures of several diphenyl ethers inhibitors, including TCL and TCU ( Figure 1 and Table 1), in ternary complexes with InhA and the oxidized cofactor form (NAD+), solved by X-ray diffraction and available in the Protein Data Bank (PDB; http://www.rcsb.org/pdb/) [22], allowed researchers to describe the main H-bonding and hydrophobic enzyme-inhibitor interactions in the substrate binding pocket [23,24]. Importantly, these data come from crystallization experiments using only the oxidized form of the cofactor, where the enzyme-cofactor-inhibitor ratio is 1:5:200 [25].  50 ) was converted to the inhibition constant (K i ) [29], which was converted to pK i (i.e., -LogK i ). n.d = not determined.
In the IC50 assays, fixed concentrations of InhA, the substrate mimic (DD-CoA, Figure 1), and reduced cofactor are used to evaluate the conversion rate from NADH to NAD+ [25]. For example, TCL shows IC50 = 1000 nM [23], but this inhibitor has a greater affinity for the enzyme-NAD+ complex with NAD+ pre-incubation, according to kinetic studies with ENR from E. coli and B. napus (FabI) [30], and M. tuberculosis (InhA) [31]. TCU shows IC50 = 5.3 and 50.3 nM at InhA concentrations equal to 10 and 100 nM, respectively [25], in a kinetic study considering only the oxidized cofactor step.
New derivatives of this class have been synthesized and evaluated [16], but none showed IC50 better than that reported for TCU. It is noteworthy that these studies aim to find compounds which dissociate slowly from the InhA-NAD+ complex generated after catalysis [16,25,28]. This proposed inhibition mechanism is related to α-helix-6 motion, in the substrate-binding pocket, via a slow conformational conversion from closed to open states. However, in these studies, it was not possible to determine clearly the influence of structural changes in the diaryl ethers class with this motion. Furthermore, some authors argue that maintaining the α-helix-6 structure is directly related with the inhibitor residence time in the InhA-cofactor-inhibitor complex, influencing the biological response [23,25,28,[32][33][34]. Recent isothermal titration calorimetry and thermal melting studies with other classes of InhA inhibitors, such as methyl-thiazol [35] and NITD-564 [36], have shown the importance in evaluating whether inhibition occurs with the enzyme apo (free InhA) or holoenzyme (InhA-NADH or InhA-NAD+) forms. Both inhibitors, methyl-thiazol (IC50 = 3 nM) [35] and NITD-564 (IC50 = 590 nM), bind preferentially to InhA-NADH [36].
Moreover, homologous enzymes to InhA, such as FabI from E. coli, B. natus, and P. falciparum, show marked differences [10,25,26,[37][38][39] The hypothesis that diphenyl ethers bind to the enzyme-NAD+ complex generated after catalysis, which is consistent with kinetic data from in vitro assays, does not rule out the hypothesis that, in vivo, the inhibitor binds to the enzyme-NADH complex before catalysis, since there is a gap between in vitro assays and the actual biological environment. Therefore, in the current work, molecular dynamics (MD) simulations in aqueous solvent of the holoenzyme (InhA-cofactor), holoenzyme-substrate, and 10 holoenzyme-inhibitor systems were performed considering the cofactor reduced form, in order to gain more insight about the solvent influence on the H-bond and hydrophobic interactions, and the dynamic behavior of the secondary structures that compose the binding site. The results obtained from this MD study could help to design new, more potent and effective InhA inhibitors in order to improve the pharmacological treatment against TB.

Results and Discussion
In order to evaluate the main structural changes that lead to the inhibition of the NADH-dependent enoyl-ACP reductase enzyme from Mycobacterium tuberculosis (InhA), we carried out up to 20 or 40 ns of molecular dynamics (MD) simulations of 12 aqueous protein systems, using the GROMACS 4.5 package, by investigating one binary (holoenzyme) and 11 ternary (holoenzyme-ligand) complexes. The ligands are THT (a substrate mimic) and 10 potent InhA inhibitors from the diphenyl ether class, including TCL (triclosan), which were divided into two groups: TCL derivatives (TCL, JPL, JPM, JPJ, and 8PC) and alkyl diphenyl ethers (5PP, 8PS, TCU, JUS, and 1TN). The starting point for the MD simulations are the corresponding X-ray crystal structures of the binary (InhA-cofactor) and ternary (InhA-cofactor-ligand) complexes available in the Protein Data Bank (PDB). The holoenzyme (InhA-NADH) and holoenzyme-substrate (InhA-NADH-THT) complexes were used as reference for the evaluation of the holoenzyme-inhibitor complexes. The MD study was conducted in order to gain more insight into the enzyme dynamic behavior before and after binding to the ligand (substrate or inhibitor), focusing on the secondary structure motifs close to the substrate binding pocket of the InhA active site. Therefore, in order to understand the structure-activity relationships (SAR) into the design of new and more potent anti-TB drugs, it was monitored the lifetime of the main intermolecular interactions in these systems: hydrogen bonds (including water-bridge H-bonds) and hydrophobic contacts between protein and ligands.
The structural stability of both protein complexes, 2AQ8 (InhA-NADH) and 1BVR (InhA-NADH-THT), was monitored ( Figure 3). In 2AQ8 ( Figure 3A), the InhA (black line) and cofactor (green line) atoms fluctuations tend to reach a plateau before 10 ns, but only after 20 ns of simulation the protein and NADH RMSD values reached the established criteria (rectangle). Analyzing the conformational change of NADH in the last 5 ns of simulation, the main variations occur in the pentose and phosphate groups. In 1BVR ( Figure 3B), the protein reach a plateau at approximately 10 ns of simulation; however, according to the established criteria, the RMSD fluctuation reached a value of less than 1.5 Å faster than 2AQ8. Analyzing the last 5 ns, the cofactor heavy atoms showed RMSD values below 0.5 Å, with the major variations occurring in the nicotinamide group located near the substrate mimic. Regarding THT, the average RMSD is 1.6 Å, but the irregularities shown in Figure 3B are due to large conformational freedom of the linear alkyl group, since the alkene group has little deviation of atomic positions in the last 5 ns of simulation. In order to obtain more structural details about the dynamic behavior of the protein, the local structure flexibility was monitored by means of the spatial RMSD (i.e., the RMSF of the α-C atoms throughout the last 5 ns of simulation), aiming to evaluate the local atoms fluctuation mainly close to the substrate binding pocket (Figure 4). Therefore, the spatial RMSD of 2AQ8 and 1BVR are depicted as a tubular backbone model, where less flexible regions (low RMSD values) are represented by thinner tubes, while more flexible regions (high RMSD values) are represented by wider tubes.
The analysis of the spatial RMSD of 2AQ8 ( Figure 4A) shows little fluctuation of residues belonging to the substrate binding pocket. However, AH-6 suffered partial denaturation changing to a 310-helix conformation (Leu197-Ser200), while AH-7 suffered partial unfolding (Glu209-Ile215). In 1BVR ( Figure 4B), THT seems to increase the mobility of LP-6 and AH-7, mainly comprising Val203 to Gly221, without changing the secondary structure of AH-6. Increased mobility in these two domains may be related to large fluctuation of the linear alkyl group of THT due to its high conformational freedom. Interestingly, studies with EcFabI-NADH-TCL (PDB ID: 1QSG) [40] and FtuFabI-NAD-TCL (PDB ID: 3NRC) [41] show a type of lid, above the nicotinamide ring, which prevents access to the solvent. In the case of InhA, this lid would correspond to the intersection region of LP-6 and AH-7. In order to characterize different conformational states (e.g., open and closed) [33,34], related to a lid protecting access to the substrate binding pocket, interatomic distances between the pairs of residues Phe97/Ala198 (LLP-4/AH-6) and Ile105/Ala206 (MLP-4/LP-6) were measured in the average structures of 2AQ8 and 1BVR, which are shown superimposed in Figure 4C. The LLP-4/HA-6 and MLP-4/LP-6 distances ( Figure 4D and Table 2) are smaller in InhA-NADH (9.71 and 7.19 Å, respectively) than in InhA-NADH-THT (14.95 and 10.06 Å, respectively). These data indicate that InhA adopts a closed conformational state in the InhA-NADH complex, where a lip protects the substrate-binding pocket ( Figure 4A), and an open conformational state in the InhA-NADH-THT complex, where the substrate is exposed to the solvent due to the large distance between MLP-4 and LP-6. Together, these results point to the importance of evaluating conformational changes [42][43][44] and structural disorder [45,46] of the enzyme in order to understand the binding affinity and dissociation of the ligand.

Comparative Analysis of the MD Simulations of the Holoenzyme-Inhibitor Complexes
Several experimental [16,23,25,32] and theoretical [33,34,47,48] studies of InhA in complex with inhibitors from the diphenyl ether class use as a reference an X-ray crystal structure containing the cofactor in its oxidized form (NAD+). Most of the InhA-cofactor-inhibitor complexes available in PDB were obtained by crystallization assays from an enzyme solution (saturated with an inhibitor) containing a molar excess of the cofactor in its oxidized form (NAD+) [49]. However, kinetic studies of InhA with TCL demonstrate that the interaction also occurs with NADH, and the TCL binding to InhA (wild-type) is uncompetitive with respect to both NADH and DD-CoA [50]. Therefore, all MD simulations in our work were carried out considering the cofactor in its reduced form (NADH) in order to evaluate the possibility of these diphenyl ether inhibitors bind to the InhA-NADH complex. Moreover, to prevent misinterpretations on the actual potency of these inhibitors, due to variations in the enzyme concentration on the IC50 assays, these values were normalized by converting to pKi [29].
In the temporal RMSD graphs of the holoenzyme-inhibitor complexes of the most potent (TCU, pKi = 9.86 M) and least potent (TCL, pKi = 6.16 M) diphenyl ethers of this series ( Figure 5), NADH maintains a constant and low fluctuation level during the entire simulation time, reaching a stability plateau, in both cases, near 10 ns. In the last 5 ns of the 2X23 and 2B35 complexes simulation (rectangle, Figure 5), the protein backbone has an average RMSD equal to 1.32 Å, while TCU and TCL have an average RMSD equal to 0.85 and 0.30 Å, respectively. In InhA-NADH-TCU (2X23), the greater atomic fluctuation occurs in the alkyl group, while for InhA-NADH-TCL (2B35), the greater fluctuation corresponds to the phenol ring.  According to the spatial RMSD graphic of 2X23 ( Figure 6), containing the most potent alkyl diphenyl ether derivative (TCU), which is also the most potent of the entire series, there is low mobility of the active site residues (rectangle), where the secondary structures of AH-6 and AH-7 are maintained. The LP-6/MLP-4 and AH-6/LLP-4 distances ( Table 2, Figures 4D and 6) observed in 2X23 (6.16 and 9.12 Å, respectively) are even lower than those found in 2AQ8 (InhA-NADH; 7.19 and 9.71 Å, respectively), likewise corresponding to a closed state of the enzyme. Due to spatial proximity between MLP-4 and LP-6, the inhibitor should remain protected from the solvent. However, this cavity must have some stability due to H-bond interactions between Ile105(N) and Ala206(O) with 79.7% of occupancy.
The LP-6/MLP-4 and AH-6/LLP-4 distances ( Figure 4D and Table 2), which are observed in 4OIM and 4OXY (Figure 6), corresponds to an open conformation, where the lid opening (LP-6) leaves JUS and 1TN exposed to solvent, thus reducing their affinities to the holoenzyme. In 4OIM (Figure 6), the secondary structure of AH-6 is maintained from Met199 to Ile202, but there is a 310-helix from Gln100 to Gly102 in MLP-4 due to H-bonds comprising Met103(N)-Gln100(O) and Gly102(N)-Pro99(O), with occupancies up to 60.6% and 96.4%, respectively. In 4OXY (Figure 6), AH-6 loses its secondary structure due to an increased movement of Leu197-Gly212, associated with the increase of the LLP-4/AH-6 distance, thus exposing the inhibitor to the solvent. In 2B36 and 2B37, whose corresponding inhibitors (5PP and 8PS) differ only in the length of the 5-alkyl phenol chain (pentyl and octyl, respectively), the MLP-4/LP-6 and AH-6/LLP-4 distances are smaller (but close) than those found in 2AQ8 (Table 2), i.e., the enzyme is in its closed conformation. Nevertheless, there is an increased mobility at AH-6 and LP-6 with partial unfolding of AH-6 (Ser200-Ala201), while in 2B37, even with the low mobility of AH-6 and LP-6, there is also a partial unfolding of AH-6 with the formation of a 310-helix from Gln100 to Gly102 in the MLP-4 region ( Figure 6).  According to the spatial RMSD graphic of 2B35 (InhA-NADH-TCL, Figure 7), containing TCL, which is the least potent compound of this subgroup (and of the entire series), there is an increased mobility of residues from ULP-4, AH-6, and LP-6, with a 310-helix formation in MLP-4. In addition, there is a H-bond interaction between Gln100(NE2) and Ile202(O) from MLP-4 and AH-6, respectively, with an occupancy of 72.8%. This H-bond along with the increased mobility of LP-6 could be responsible for the overall AH-6 unfolding. Furthermore, in 2B35, the binding site adopts a closed conformation ( Figure 4D) as in holoenzyme, according to AH-6/LLP-4 and LP-6/MLP-4 distances ( Table 2). In fact, these results are consistent with those described for EcFabI-NADH-TCL [40] and FtuFabI-NAD-TCL [41], where the lid closing (due to the movement of AH-6) maintains the inhibitor inside the binding site.
In 3FNE (InhA-NADH-8PC) ( Figure 4D and Table 2), the LP-6/MLP-4 and AH-6/LLP-4 distances indicate that the enzyme assumes an open conformation due to the high mobility of ULP-4, AH-6, and LP-6, but maintaining the secondary structure of AH-6, while MLP-4 assumes a 310-helix fold (Figure 7). In 3FNF and 3FNG, whose corresponding inhibitors (JPM and JPL) differ only in the 5-position substituent of the phenol ring (benzyl and cyclohexylmethyl, respectively), ULP-4, AH-6, and LP-6 have lower mobility than in 3FNE (Figure 7). In addition, the AH-6/LLP-4 and LP-6/MLP-4 distances are even lower than those found in 2AQ8 ( Figure 4D and Table 2), indicating that despite the substituent volume, the enzyme maintains a closed conformation. Finally, in 3FNH, containing JPJ (which differs from JPM and JPL by the presence of the 5-phenetyl substituent), we observed regions of low-fluctuation from the substrate binding pocket except for UPL-4 ( Figure 7). In addition, the AH-6/LLP-4 and LP-6/MLP-4 distances are similar to those found in 3FNF and 3FNG ( Figure 4D and Table 2), indicating again that, despite the substituent volume, the enzyme maintains a closed conformation.

Analysis of Holoenzima-Inhibitor Interactions by Hydrogen Bond and Hydrophobic Contact
The analysis of the H-bond and hydrophobic contact interactions occurring between the holoenzyme (InhA-NADH) and inhibitors is crucial to understanding the structural changes observed and discussed above. According to literature [23,24], Tyr158 (LP-5) is described as responsible for H-bond interaction in the substrate binding pocket, while eleven residues distributed in LP-5 (Phe149, Met155, Pro156, and Ala157), AH-5 (Met161), AH-6 (Ala198, Met199, Ile202, and Val203), and AH-7 (Leu218 and Trp222) are described as responsible for hydrophobic contacts.
Thus, in order to obtain more insight about the binding mode of these inhibitors and to identify residues that actually can contribute to explain the biological activity profile of this series of compounds, it was monitored the lifetime (i.e., occupancy) of the main intermolecular interactions in these systems during the MD simulations, i.e., hydrogen bonds (H-bonds) between holoenzyme and inhibitors, including indirect H-bonds mediated by water molecules (water-bridge), and hydrophobic contacts. Table 3 shows the monitored interactions (H-bond and hydrophobic contact) between all diphenyl ethers inhibitors and InhA-NADH during the last 5 ns of simulations. Figures 8 and 9 show the close view of the substrate binding pocket with the alkyl diphenyl ethers subgroup (TCU, JUS, 1TN, 5PP, and 8PS) and TCL and its derivatives (8PC, JPM, JPL, and JPJ), respectively.      Considering the alkyl diphenyl ethers subgroup, in InhA-NADH-TCU (closed conformation) (2X23, Table 3, Figure 8A), there is an H-bond between the phenol hydroxyl group of TCU and the Gly96 backbone (LLP-4) with occupancy of 26.5%. In addition, this inhibitor also makes hydrophobic contacts of high and low occupancies with residues from MLP-4 (Met103, 92%) and AH-6 (Ala198, 25.2% and Ile202, 25.6%), respectively. In InhA-NADH-JUS (open conformation) (4OIM, Table 3, Figure 8B), there are H-bonds mediated by a water molecule (water-bridge H-bond) between the phenol hydroxyl group of JUS and Tyr158 (UAH-5), which are of high occupancy (HOH-inhibitor, 89.2%; HOH-Tyr158, 98.2%). This inhibitor also makes hydrophobic contacts of high and low occupancies with residues from UAH-5 (Tyr158, 71.1%) and AH-6 (Ile202, 2.2%), respectively. In InhA-NADH-1TN (open conformation) (4OXY, Table 3, Figure 8C), there is an H-bond between Tyr158 and the phenol hydroxyl group of 1TN with high occupancy (91.4%). The inhibitor also makes hydrophobic contacts of high and medium occupancies with residues from UAH-5 (Met161, 82.4%) and AH-6 (Val203, 37%), respectively. Comparing 4OXY and 4OIM, the absence of a solvent-inhibitor interaction may be related to maintaining the lid closed, even with the AH-6 motion, thus preventing solvent-inhibitor interaction.

Steric and Electrostatic Interactions of the Holoenzyme-Inhibitor Complexes
The correlation of the biological response (pKi) with the steric (ELJ), electrostatic (EC) and total (ELJ + EC) interaction energies from the MD simulations was done by the PLS/GFA analysis using the WOLF program [51]. In the best linear equation (pKi = 7.772 − 0.1885 × Gly96 + 0.0517 × Tyr158; R 2 = 0.80; n = 10; Figure 10) the selected descriptors, Gly96 and Tyr158, correspond to the total energy (ELJ + EC). Descriptor contribution to the activity (pKi) should be analyzed by considering the energy value and its coefficient signals (positive and negative). As the Gly96 descriptor has negative signals of value and coefficient, it increases the pKi, i.e., this descriptor contributes increasing the potency. However, Tyr158 descriptor decreases the potency, since it has negative value and positive coefficient. Simple correlation analysis between pKi and independent variables (Supplemental Material) shows high correlation (R ≥ |0.80|) between Gly96 descriptor and activity (R = −0.85) ( Figure 10A), mainly due to electrostatic term (Gly96C, R = −0.82), since there is not significant correlation between pKi and steric term (Gly96LJ, R = −0.38). According to Table 3, only TCU interacts with Gly96 for 26.5% of the time by H-bond, thus, other types of electrostatic interaction are more important in explaining the high correlation. In the case of Tyr158 descriptor (Figure 10B), the correlation with pKi is not as high (R = 0.64) and is not related separately to any of these terms: steric (Tyr158LJ, R = 0.20) and electrostatic (Tyr158C, R = 0.39). According to Table 3, Tyr158 makes both types of interactions: H-bond (e.g., 73.6% of the time with JPL) and hydrophobic (e.g., 71.1% of the time with JUS). Despite not being selected in the best equation, Phe97 descriptor also shows some correlation with activity (R = −0.60), mainly due to steric contribution (Phe97LJ, R = −0.65), but not to electrostatic term (Phe97C, R = −0.24).
In the cross correlation matrix of the descriptors ( Figure 10B), it is possible to identify pairs of descriptors that may be contributing redundantly in the explanation of the biological response. This matrix has been calculated only for the descriptors that have R ≥ |0.60| with activity, i.e., Gly96 (EC and ELJ + EC), Phe97 (ELJ and ELJ + EC) and Tyr158 (ELJ + EC). The best equation descriptors (Gly96 and Tyr158) are poorly correlated (R = −0.48), indicating that each descriptor provides almost unique information to the model. The Gly96/Phe97 (R = 0.66) and Gly96/Phe97LJ (R = 0.65) descriptors are not highly correlated, but have some correlation, which can be explained by the spatial proximity of these residues. On the other hand, highly correlated descriptors, such as Phe97/Phe97LJ (R = 0.97) and Gly96/Gly96C (R = 0.93), correspond to the same amino acid. Figure 10 shows the holoenzyme in complex with the most ( Figure 10C) and least ( Figure 10D) active compounds, highlighting only the amino acids with R ≥ |0.60| with pKi. The proximity of TCU with Gly96 ( Figure 10C) causes conformational change of the cofactor adenine group, favoring a pi-stacking interaction between Phe97 and the purine ring, while for TCL ( Figure 10D), this conformational change is not observed.

Final Considerations
The analysis of the intermolecular interactions of the holoenzyme-inhibitor complexes allows us to establish a qualitative relationship between the maintenance of secondary structures that compose the substrate-binding pocket, particularly AH-6 (α-helix-6), and the relative potency of inhibitors. In general, the most potent compounds (high pKi values) interact through H-bonds directly with Gly96 backbone (TCU, pKi = 9.86 M) or the lateral chain of Lys165 (8PS, pKi = 8.47 M) or Tyr158 (JPM, pKi = 7.93 M and 5PP, pKi = 7.91 M. On the other hand, the less potent compounds (low pKi values) interact by H-bond only with NADH cofactor (TCL, pKi = 6.16 M) or they also interact with enzyme amino acids, but only through water-bridge H-bonds (JUS, pKi = 6.77 M). This behavior is in agreement with literature, where TCL binds weakly to holoenzyme when the cofactor is reduced, since no significant interaction was identified [30,31,40].
Regarding to the enzyme conformational state, in the holoenzyme-inhibitor complex, the open conformation is favored when the inhibitor makes hydrophobic contacts of high occupancy with residues from AH-5 (Tyr158 or Met161). In addition, water-bridge H-bonds only occurs when the mobility of LP-6 and AH-7 intercession increases, which favors an open conformation, consequently, there is a decrease in the biological response. Furthermore, when inhibitor makes hydrophobic contacts with only one residue from AH-6 (Ala198, Ile202, or Val203) with occupancy of 30% or more, it induces the unfolding of this secondary structure, while hydrophobic contacts with two residues or more favor the maintenance of the AH-6 secondary structure. Overall, the most potent compounds are those in which the enzyme assumes a closed conformational state, similar to holoenzyme, as is the case of the four most potent inhibitors: TCU (pKi = 9.86 M), 8PS (pKi = 8.47 M), JPM (pKi = 7.93 M), and 5PP (pKi = 7.91 M). On the other hand, the less potent compounds are usually those in which the enzyme has an open conformational state (somewhat similar to that adopted in the enzyme-cofactor-substrate complex), as is the case of the two less potent inhibitors: TCL (pKi = 6.16 M) and JUS (pKi = 6.77 M).
Regarding the steric (ELJ) and electrostatic (EC) interactions between inhibitors and residues, according to the linear equation shown in Figure 10 (pKi = 7.772 − 0.1885 × Gly96 + 0.0517 × Tyr158; R 2 = 0.80; n = 10), the total energies (ELJ + EC) related to Gly96 and Tyr158 are able to explain 80% of the biological response variance. Interactions with Gly96, especially electrostatic, contribute to increasing the pKi value, while both types of interactions with Tyr158 decrease the pKi value.

Structure of Protein Complexes Selection
The X-ray crystal structures of InhA in complex with cofactor (NADH/NAD+), a substrate mimic (THT), and 10 inhibitors (diphenyl ethers derivatives) were retrieved from the Protein Data Bank (PDB), using the following PDB IDs  [24]; 4OIM (InhA-NAD-JUS; R = 1.85Å) [16]; and 4OXY (InhA-NAD-1TN; R = 2.35Å) [28]. As the PDB files 2B35 and 2B37 had missing atomic coordinates of amino acid residues located in the active site, 2X23 was used for the molecular docking simulations of the corresponding TCL and 8PS inhibitors.
The 3D structures of the ligands (cofactor, substrate, and inhibitors) were extracted from their respective complexes with InhA according to the following PDB ligand codes: NAI (NADH, reduced cofactor form, 1,4-dihydronicotinamide-adenine-dinucleotide), THT (a substrate mimic, trans- In all calculations, it was used the cofactor in its reduced form (NADH). The 3D structures of two diphenyl ethers inhibitors (1TN and JPL) have been corrected due to disagreements with experimental data published in the original articles [16,24]. The n-pentyl group of 1TN was corrected as n-hexyl group, while the 1,5-dien-1-ylmethyl group of JPL was corrected as cyclohexylmethyl group.

Biological Data
The biological activity (measured as the half-maximal inhibitory concentration, IC50) of the diphenyl ethers derivatives was compiled from various papers and it was determined at different concentrations of the enzyme (InhA), cofactor (NADH), and C12 fatty acyl substrate (2-trans-docenoyl-coenzyme A, DD-CoA) [16,[23][24][25]28]. This concentration variation, mainly of enzyme, may lead to misinterpretations of the actual inhibitor potency, as in case of TCU having two IC50 values (5.3 and 50.3 nM) depending on the InhA concentration (10 and 100 nM, respectively) [25]. Therefore, to avoid comparisons of IC50 values measured at relatively different experimental conditions, these data were normalized by their conversion to the inhibition constant (Ki) according to Equation 1 [29]. This is because the diphenyl ethers derivatives are classified as uncompetitive inhibitors that binds slowly in the same cavity of the substrate [50].
In Eqation 1, Ki is the inhibition constant; IC50 is the half-maximal inhibitory concentration; E and S are the concentrations of enzyme (E) and substrate (S), respectively; and Km is the Michaelis-Menten constant for each E-S system. In the InhA inhibition assays using DD-CoA as substrate and NADH as cofactor, Km is equal to 29 µM [23]; and the InhA (E) and DD-CoA (S) concentration values varied between 1-100 µM and 25-300 µM, respectively [16,[23][24][25]52]. Finally, to allow a linear and direct comparison of the biological activity of these derivatives, the Ki value was converted to its corresponding negative logarithm (pKi = −LogKi) expressed in Molar (M) unit (Table 1).

Molecular Docking Simulation
Before the molecular dynamics (MD) simulations, molecular docking was performed for TCL, 8PS, 1TN, and JPL in order to obtain a starting atomic coordinate of them with holoenzyme (InhA-NADH). Firstly, the 3D structures of these compounds were submitted to the default conformational analysis process of the Spartan program (version 10, Wavefunction, Inc., Irvane, CA, USA) [53], using the RM1 semi-empirical method [54]. The most stable conformation of each compound was submitted to a single-point calculation by the DFT/B3LYP method, using the 6-31+G(d) base function available on Spartan, in order to derive the partial atomic charges.
These compounds were docked within the InhA-NADH complex (PDB ID: 2X23), excluding TCU inhibitor and all water molecules, using the Molegro Virtual Docker (MVD) program [55], according to the previously works published by our group [56,57]. The MVD automatic preparation module was used to correct the atoms types and bond orders, to add the hydrogen atoms and to assign the default atomic partial charges in the protein structure.
Potential binding sites (cavities) were detected using the grid-based cavity prediction algorithm. The population size, maximum interactions, scaling factor, and crossover rate were set to 150, 1500, 0.50, and 0.90, respectively. For each complex, we performed 100 independent runs with the MolDock optimizer algorithm, returning five solutions (poses, i.e., conformation and orientation) for each run. The MolDock score function with a grid resolution of 0.30 Å was used to precompute score grids for rapid dock evaluation. Guided differential evolution and a force-field-based docking scoring function were used to search for the binding orientation and conformation of each candidate molecule. The best pose of each inhibitor was selected for the subsequent MD simulations.
To guarantee the reliability of the docking protocol previously described, the TCU inhibitor was removed and docked back into the active site of InhA-NADH. The root-mean-square deviation (RMSD) calculated between the best pose and the original coordinates was less than 2 Å, an acceptable value as suggested by different authors in literature [58][59][60][61][62].
The ligands (cofactor, substrate, and inhibitors) topology were created using the AnteChamber PYthon Parcer InterfacE (ACPYPE) tool [67], and the ACPYPE missing parameters were generated by the MKTOP software [68]. The partial atomic charges were assigned according to the AM1-BCC parameters [69], through the ANTECHAMBER package [70,71]. The holoenzyme and holoenzyme-ligand complexes were inserted and centered into the cubic periodic box of spc216 water molecules and it was considered the SPC/E water model [72], thus each complex was neutralized by six Na + counter-ions. Each complex was submitted to a preliminary optimization using the steepest descent algorithm with position restrained (PR) of the ligands (cofactor, substrate, and inhibitors) and convergence criteria of 1000.00 kJ·mol −1 ·nm −1 , followed by steepest descent without PR, and conjugate gradients until an energy of 100.00 kJ·mol −1 ·nm −1 . Then, the minimized complexes were submitted to 1000 ps of MD, at 300 K, considering NVT and NPT state with PR to the entire system, except to water molecules, using the V-rescale thermostat [73], and the Parrinello-Rahman scheme [74] for pressure coupling. The PR was applied to guarantee the distribution of the solvent molecules around the protein. All bonds involving hydrogen atoms in the complexes were frozen by the LINCS scheme [75], the long-range electrostatic interactions were treated using the PME algorithm (Particle-Mesh Ewald) [76,77], and a cut-off value of 1 nm was applied for the van der Waals and Coulomb interactions. Afterward, it was carried out 20 to 40 ns of MD simulations without any restriction, using 2 fs of integration time and a cut-off of 10Å for long-distance interactions.
The criteria used to finalize the MD simulations followed a rigorous protocol to determine if the system is indeed in an equilibrium state. Thus, from 10 ns, we calculated the RMSD value during every 5 ns (e.g., from 10 to 15 ns, from 15 to 20 ns, etc.), verifying if the RMSD variation was lower than 1.5 and 1.0 Å for the protein backbone and ligand, respectively. In case of InhA-NADH-THT, we establish the RMSD variation should be less than 2 Å, due to the high degree conformational freedom of the substrate mimic.

Molecular Dynamics Analysis
The root mean standard deviation (RMSD) calculation, cluster analysis, mean structure extraction, interatomic distance measurements, and hydrogen bond analysis were performed using diverse modules available in GROMACS 4.5 [63][64][65]. The H-bonds (D-H…A) were computed considering the cutoff distance between donor (D) and acceptor (A) atoms until 0.35nm and the cutoff H-D-A angle until 30°. The H-bonds frequency was calculated using hbmap2grace package [78]. The cutoff value for distance (nm) between non-hydrogen atoms of ligand, cofactor and residues was until 0.40nm. The occupancy percentage (%) for H-bonds and interatomic distances were calculated only for the last 5 ns of simulation. The 3D structure was rendering using VMD program [79]. The Grace program was used to plot the RMSD graphics [80]. The 3D root mean square fluctuation (RMSF) or sausage representation was created for the last 5 ns, using the MOLMOL package [81].

Steric and Electrostatic Interactions Calculations of the Holoenzyme-Inhibitor Complexes
The steric (Lennard-Jones, LJ) and electrostatic (Coulomb, C) interaction energies (ELJ and EC) between the inhibitors and holoenzyme were calculated only for the binding site residues, using the average energy of the 100 conformations of the most populated cluster from the last 5 ns of each simulation. In order to estimate the influence of these descriptors (independent variables) on the biological response (pKi) (dependent variable), we considered the steric (ELJ), electrostatic (EC) and total (ELJ + EC) average energies by residue for each inhibitor. The 192 descriptors resulting were submitted to a combined Partial Least Squares (PLS) and Genetic Function Approximation (GFA) analysis available in the WOLF program [51]. In the PLS/GFA calculations, it was considered three principal components, population of 600-1000 equations, crossover operations of 50,000-500,000, mutation rate of 100%, and smoothing factor of 0.1-1.0.
The variance (σ 2 ) of the pKi values shows a small dispersion (σ 2 = 1.02), which means that descriptors with small variance do not contribute with relevant information to the energy-activity relationship. Thus, we discarded descriptors with σ 2 < 0.30. Then, we submitted the 71 descriptors resulting to the WOLF program, in order to generate equations that may explain the biological response. Due to the limited number of compounds (n = 10), we selected only equations containing two descriptors. In order to select descriptors correlated with pKi and intercorrelated, we calculated the simple correlation (R) by the cross-correlation matrix of the descriptors and pKi (Supplementary Material).

Conclusions
The MD simulations showed the importance of the H-bond and hydrophobic contact interactions between holoenzyme and ligands, responsible for structural changes and the steric and electrostatic interaction energies to identify important residues correlated with the biological response of diaryl ethers derivatives. Analysis of 2AQ8 and 1BVR revealed that the open and closed conformational states of InhA are related to the AH-6 motion, while the solvent-inhibitor interaction is related to the mobility of LP-6/AH-7 intersection that works as a lid. The enzyme assumes an open conformation mainly when the inhibitor makes: (i) H-bonds with the cofactor, (ii) indirect H-bonds (water-bridge H-bonds) with enzyme residues as a consequence of the lid opening, and (iii) hydrophobic contacts (>30% of occupancy) with AH-5. All these interactions are related to a decrease in the inhibitor affinity by the holoenzyme. Moreover, the protein-inhibitor total steric (ELJ) and electrostatic (EC) interaction energies, related to Gly96 and Tyr158, are able to explain 80% of the biological response variance according to the best linear equation, pKi = 7.772 − 0.1885 × Gly96 + 0.0517 × Tyr158 (R 2 = 0.80; n = 10), where interactions with Gly96, mainly electrostatic, increase the biological response, while those with Tyr158 decrease. Taken together, our results demonstrate the importance in studying the enzyme-cofactor-inhibitor dynamic behavior, considering the structural changes, solvent influence, and interaction energies for understanding the biological activity. These results will help to understand the structure-activity relationships and to design new and more potent anti-TB drugs.