Molecular Dynamics Assisted Mechanistic Study of Isoniazid-Resistance against Mycobacterium tuberculosis InhA

Examination of InhA mutants I16T, I21V, I47T, S94A, and I95P showed that direct and water mediated H-bond interactions between NADH and binding site residues reduced drastically. It allowed conformational flexibility to NADH, particularly at the pyrophosphate region, leading to weakening of its binding at dinucleotide binding site. The highly scattered distribution of pyrophosphate dihedral angles and chi1 side chain dihedral angles of corresponding active site residues therein confirmed weak bonding between InhA and NADH. The average direct and water mediated bridged H-bond interactions between NADH and mutants were observed weaker as compared to wild type. Further, estimated NADH binding free energy in mutants supported the observed weakening of InhA-NADH interactions. Similarly, per residue contribution to NADH binding was also found little less as compared to corresponding residues in wild type. This investigation clearly depicted and supported the effect of mutations on NADH binding and can be accounted for isoniazid resistance as suggested by previous biochemical and mutagenic studies. Further, structural analysis of InhA provided the crucial points to enhance the NADH binding affinity towards InhA mutants in the presence of direct InhA inhibitors to combat isoniazid drug resistance. This combination could be a potential alternative for treatment of drug resistant tuberculosis.


Introduction
The clinical isolates of mycobacterium tuberculosis bacilli obtained from the isoniazid (INH)resistant tuberculosis (TB) patients are observed to have mutations in the structural gene of InhA protein [1][2][3][4]. INH react with NADH to form INH-NAD adduct which binds into dinucleotide binding site (DBS) of InhA. INH only blocks the site of hydride transfers by interacting with NADH that causes to inhibit the reduction of fatty acyl substrate. This is the molecular mechanism which makes INH an indirect inhibitor of InhA. However, INH does not contribute in INH-NAD binding into DBS. Thus, the observed mutations at DBS only affect NADH binding which ultimately leads to diminished INH-NAD binding affinity [5][6][7][8]. NADH interacts with DBS through H-bonds with side chain of residues Asp64, Ser20 and Lys165, and main chain of residues Val65, Gly14, Ile194, and Ile21 [7]. In addition to this, Dessen et al. reported water molecules in the DBS of InhA to mediate H-bond interactions between NADH and binding site residues. One of the water molecules is observed near the pyrophosphate moiety of NADH which mediates H-bonds with residues Ser94, Gly14 and Ala22. However, mutation of residue Ser94 into Ala94 causes disruption of bonding pattern which is accounted for considerable changes in NADH conformation. It is consistent with biochemical and genetic studies carried out by Basso et al., where a higher dissociation constant (K d ) for NADH is observed for mutants I16T, I21V, I47T, S94A, and I95P as compared to InhA WT [9]. Oliviera et al. also reported the enhanced K d values for I21V, I47T and S94A with the help of kinetic study [6]. Schroeder et al. also demonstrated the enhanced NADH conformational flexibility in I16T and I21V mutants using MD simulations [10]. These experimental and computational studies emphasize the influence of mutations on NADH binding which leads to INH-resistance. However, mutants I47T, S94A, and I95P are not explored at atomic level for deducing the INH-resistance using computational approaches. In addition to this, all mutant residues are not characterized at thermodynamic ground to compare the reduced binding energy contribution with respect to wild type residues. Many theoretical studies have shown the significance of computational approach in characterizing the mutation and ligand driven structural changes in vital proteins and their complexes [11][12][13]. MD simulation is one of the computational approaches which can be used effectively in the thermodynamic assessment of mutated residues on ligand binding and conformational changes [14][15][16][17]. The reported experimental studies show that mutations do present and plays a critical role in INH-resistance. However, the reported molecular modeling studies do not cover all InhA mutants. In addition, these studies also do not characterize the mutations on energetic ground of NADH binding.
The present study characterizes the mechanism of INH-resistance with respect to structural changes and binding energy of NADH using MD simulations. A deep analysis may substantiate the functional role of mutations in INH-resistance that can provide crucial structural insights to improve the INH binding. In addition to this, these structural insights can be viewed with respect to direct InhA inhibitor (DII) binding site adjacent to the DBS. It may lead us to explore the application of DIIs to enhance the INH binding in InhA mutants.
System preparation for molecular dynamics simulations. Preparation of the topology and coordinate files for the all complexes was carried out using the selected InhA crystal structure complexes. The obtained structures may have missing bond order, connectivity, steric clashes or bad contacts with the neighboring residues. Therefore, the selected structures were corrected for atoms and bonds and energy minimize to potentially relax the structures. It resolved steric hindrance and clashes in the structures. FF99SB force field parameters were set for the protein using the AMBER (Assisted Model Building with Energy Refinement) LeaP module [22]. The parameters missing for the NADH were taken from the amber parameter database, university of Manchester [23,24]. The prepared systems were solvated with TIP3P water model by creating an isometric water box, where distance of the box was set to 10Å from periphery of protein [25]. Molecular systems were neutralized through the AMBER LeaP module by the addition of necessary amount of counter ions (Na + ) to construct the system in electrostatically preferred positions. The whole assembly was then saved as per the requirement of free energy calculations. It involved preparing the parameter and coordinate files for the complex, protein and the ligand without solvation. Further, the prepared topology and coordinate files of solvated complexes were used as input for sander module of the AMBER [26]. The optimization and relaxation of solvent and ions were performed by means of two energy minimization cycles using 1500 and 2000 steps. The initial 1000 steps of each minimization cycle were performed using steepest descent followed by conjugate gradient minimization for rest of the steps. In the first part of minimization, InhA-NADH complex was kept fixed to allow water and ion molecules to move, followed by the minimization of the whole system (water, ions and complex) in the second part. Heating was performed into six steps using a NVT ensemble for 120ps where InhA-NADH complex was restrained with a very small force constant of 5kcal/ mol/Å 2 . The temperature was allowed to increase by 50K in each step till 300K. The system was further equilibrated under constant pressure at 300K for the period of 100ps without restrain on the complex. Final simulations i.e. production phase were performed for 100ns on NPT ensemble at 300K temperature and 1atm pressure. The step size of 2fs was kept for whole simulation study. Langevin thermostat and barostat were used for temperature and pressure coupling. SHAKE algorithm was applied to constrain all bonds containing hydrogen atoms [27]. Non-bonded cutoff was kept on 10Å and long range electrostatic interactions were treated by Particle Mesh Ewald method (PME) with fast Fourier transform grid spacing of approximately 0.1nm [28]. Trajectory snapshots were taken at each 100ps of the production phase, which were used for final analysis. The minimization and equilibration were performed by sander module of Amber11, while production simulation was performed using Pmemd program of AMBER running on NVIDIA Tesla C2050 GPU work station [29]. The production run was considered for the analysis which was carried out using the Ptraj module of the Amber11 and VMD [30,31].
Root-mean-square deviation, B-factor analysis, dihedral and Chi1 angles. The rmsd, Bfactor, dihedral and Chi1 angles were calculated using the ptraj analysis tool in the AMBER program. The rmsd is the measure of the average distance between the atoms (usually the backbone atoms) of superimposed protein structures [32]. The equation below illustrates the rmsd: where δ is the distance between N pairs of equivalent atoms.
To compare the flexibility of the structures, B-factor was used to calculate the mobility of the residues present in the SBL. It was calculated from the mean square fluctuations (msf) using the following equation: B-factor was calculated using the Eq 2 for each residue and plotted against each residue [33]. Dihedral angles (F, C) are calculated for backbone of protein while Chi1 (χ1) dihedral angle referred to side chain dihedral angle of the residue [34]. They are calculated as per the following equations Where, N and C atoms are from NH2 and COO groups respectively, directly bonded to Cα atom of residue while Cβ represents side chain carbon atom bonded to Cα. Xγ represents the carbon/non-carbon atom bonded to Cβ. Free energy calculations. Binding free energy calculations were performed for all complex systems using MM-GBSA (Molecular Mechanics Generalized Born Surface Area) considering single trajectory approach. The various previous study find this approach useful in estimation of binding free energy [35,36]. Bello et al. identified the single trajectory approach useful for characterizing the key residues involved in valproic dehydrogenation using MM-GBSA [37]. Similarly, in another study, Bello et al. explored the binding of four fatty acids with β-lactoglobulin using the single trajectory approach and MM-GBSA method [38]. These cases depict the importance of the single trajectory approach in calculating the binding free energy using MM-GBSA. Binding free energy calculations were performed on production phase considering last 20ns of 100ns run using the Born implicit solvent model of 2 (igb = 2).
The calculations include various energy components consisting of molecular mechanical energy (E MM ) and polar contribution (G GB ) towards solvation energy calculated by generalized Born (GB) method respectively (Eq 7). G SA is the contribution from nonpolar terms towards solvation energy, and TS is the entropic contribution of the inhibitor. E MM was obtained by summing contributions from, electrostatic energy (E ele ), VDW energy (E vdw ), and internal energy including bond, angle, and torsional angle energy (E int ) using the same force field as that of MD simulations (Eq 8).
Further, effective binding free energy was calculated using the frames obtained in every 5ns of 100ns simulations. The average of binding free energy was estimated for each 5 ns and referred as effective binding free energy Energy decomposition analysis. Free energy was decomposed to estimate the contribution of each residue in the ligand binding process and was performed by using MM-GBSA method [39]. Energy of each residue-ligand interaction is given by following equation: Where, ΔG GB is the polar group contribution to the solvation free energy calculated using GB model. ΔG SA is the non-polar group contribution to the solvation free energy calculated using ICOSA method. Similar to the binding free energy, the decomposition energy was also averaged over 200 frames taken at the interval of 100ps over the last 20ns of 100ns production run. In addition, average structures for comparison of trajectories were also obtained using the 200 frames over the last 20ns of 100ns simulations.

Structural features of InhA-NADH complex
InhA catalyzes the fatty acyl substrate with the help of NADH co-factor. It is comprised of two nucleotide moieties coupled together through the phosphate groups making a pyrophosphate bridge. The nucleotide moieties of NADH consist of nitrogen base (one moiety consists of nicotinamide while other one consists of adenine), ribose sugar and phosphate group, while nitrogen base with ribose sugar is known as nucleoside ( Fig 1A). The co-factor NADH contains ten single rotatable bonds (dihedral angle: D1-D10), and 8 of them are present in the pyrophosphate bridge contributing to its flexibility. NADH binds to the DBS through a network of eleven H-bond interactions with the active site residues. DBS is comprised of hydrophobic and polar residues. The residues Gly14, Asp64, Val65 and Gly96 are observed to make H-bond interactions with the adenine ring and adenine ribose moiety, while residues Ile95, Lys165 and Ile194 are found to interact with nicotinamide and ribose moiety of NADH. The residues Ser20 and Ile21 are observed to form H-bonds with the phosphate moieties of adenine and nicotinamide nucleotide respectively. In addition to these interactions, residue Phe41 is observed to make π-π interaction with the adenine ring of NADH. All mutants show H-bond interactions similar to the WT ( Table 1). The details of H-bond donor and acceptor atoms with the distances are given in Table 1.  A distance cut-off of 3.30Å is considered for the H-bond with minimum cut-off of 140°for H-bond angle.
It is important to note that the residues that undergo mutation do not involve in direct Hbond interactions to NADH, except the Ile21 which does not show mutational effect as Val21 forms H-bond similar to WT. It may be due to the backbone atom of Ile21, which is involved in the H-bond interaction with NADH while the mutation only altered the side chain of the residue without affecting the backbone conformation. Further, these interactions are also aided by additional water mediated bridged H-bond interactions with NADH (Table 2). S94A is observed to make additional bridged H-bond interactions with NADH with the help of crystal waters WAT271 and WAT. In contrast to this, I16T, I47T, S94A, I95P and WT are not observed to form bridged H-bond interactions with these water molecules through Arg43 and, Ile95 and Lys165 through WAT396 and WAT287 respectively.
In addition to this, I21V, I47T and S94A are not observed to make bridged H-bond interactions with WAT380 through Thr196. I21V is observed with maximum bridged H-bond interactions while I47T and S94A are observed with minimum bridged H-bond interactions. Overall, a minimum 19 interactions are observed between the active site residues of DBS and NADH in the studied structures including direct and bridged H-bond interactions with residues. Among all complexes, mutant I21V is observed to have 23 interactions between the residues and NADH. Moreover, it is clear from the analysis that the maximum interactions are found near the pyrophosphate moiety of NADH. These interactions are supposed to impart rigidity to NADH conformation. Also, it is necessary to stabilize the flexibility of pyrophosphate to make a stable InhA-NADH complex.
Apart from the InhA-NADH interactions analysis, the distance between nearest atoms of mutated residues and pyrophosphate moiety of NADH is estimated. Comparison of distances shows that residues Ile16, Ile21, Ser94 and Ile95 are located within the radius of 5Å around the pyrophosphate bridge while residue Ile47 is located within a distance >6Å from pyrophosphate bridge. It suggests that the residues present in the vicinity of pyrophosphate may directly influence InhA-NADH interaction due the mutation. However, Ile47 will not probably affect InhA-NADH binding due to the mutational effect. In addition to location from NADH,

SN.
Residues 30Å is considered with a minimum cut-off of 140°for H-bond angle. *atom of NADH mutated residues are also analyzed with respect to their location in the secondary structure of the InhA. Mutations located in the secondary structure may induce some distortional changes in protein topology. The mutations I16T, S94A and I95P are observed in the loop regions, and thus, are not supposed to affect the secondary structures of the InhA. I21V is found in the helix region, however, it is observed to restore the bonding pattern similar to corresponding secondary structure in WT. Similarly, I47T is also found in the helix region and observed to restore the bonding pattern similar to WT. Thus, the observed mutations are not observed to bring the distortional changes in the secondary structure of InhA.

Comparative molecular dynamics simulations
Root mean square deviation, B-factor and effective energy. To obtain an estimate of the quality and convergence of the MD trajectory, the backbone root mean square deviation (rmsd) of InhA is calculated for the WT and MTs (Fig 2A). After a rapid increase during the first 40ns, the backbone rmsd of WT stabilized with an average rmsd of 0.72Å. The backbone rmsd of the MTs is also stabilized after 40ns of the simulation; however, it shows a higher average rmsd of~1.2Å as compared to the WT. It indicates that enhanced flexibility of the mutant structures which may have occurred due to the distorted H-bonding pattern between InhA and NADH. In addition to this, the mutational effect is also expected for the disturbed InhA-NADH interactions leading to an increase in the fluctuation of binding residues. The order of the average rmsd for WT and MTs is observed as: 44Å). This order demonstrates that InhA-NADH bonding pattern is drastically disturbed resulting in high fluctuations in backbone of MTs as compared to WT. Overall, these observations lead us to infer that there is a weakening of InhA-NADH interactions in the MTs due to mutations at the DBS. However, the average rmsd of the InhA backbone is observed below 0.8Å for the WT while it remains near to~1.5Å for the MTs. It shows the convergence of trajectory and supports its use for further analysis.
The large differences between the backbone rmsd of MTs and WT can be rationalized on the basis of higher flexibility of loop regions estimated as B-factor. The major contribution to rmsd is shown by the residues 101-107 and 194-210 corresponding to A-loop and substrate binding loop (SBL) respectively (Fig 2B). It is important to note that these loops are the part of substrate binding pocket (SBP) and play a crucial role in substrate binding. Thus, the large humps are observed for SBL and A-loop in Fig 2B which is consistent with the necessity of these regions to undergo conformational changes upon binding of fatty acyl substrate. Apart The rmsd pattern provides sufficient information about time when the system reached equilibrium. However, as suggested by Bello et al., the effective binding free energies can also provide information about system convergences to support rmsd based convergence [40]. Fig 3 shows that all six systems shows fluctuation in average effective binding energy before 60ns. However, after 60ns of simulation trajectory shows a smooth convergence into a stable energy curve. This shows that all complexes achieve a stable effective energies after 60ns.
Finally, these observations indicate the enhanced flexibility of binding residues in DBS as a result of disturbed InhA-NADH interactions. It will lead to conformational changes in NADH and binding residues particularly mutated residues which can be analyzed by measuring dihedral angle for NADH and side chains of residues.
NADH dihedral angle distributions in WT and MTs. The variations in the rmsd of WT and MTs indicate a probability of weakening the InhA-NADH interactions due to the influence of the mutation at the DBS. It leads to the enhanced flexibility of NADH as well as binding site residues. The generated hypothesis from the earlier discussion is also consistent with the reported experimental study where mutations associated with the DBS lead to an increase in the NADH dissociation constant [20].
In addition to this, side chains of the mutated residues are small as compared to that in WT which may lead to a small gap between NADH and the residues in MTs structures. This small space may be a sufficient space for the conformational flexibility of NADH. A large gap between interacting atoms will also lead to weak interactions. In case of S94A, a large difference is observed in the length of side chain and polar nature of residues which may lead to a little large gap between interacting atoms of mutated residue and NADH. Consequently, it will increase the flexibility of NADH conformation. Similarly, in case of I95P, small side chain of mutated residues may also lead to more conformational freedom to NADH. I16T also undergoes the reduced side chain length and change in polarity which is supposed to affect InhA-NADH interactions. Although, residue Ile47 is located about >6Å away from the NADH binding site, change of Ile47 into Thr47 leads to change in polarity which may affect the microenvironment around NADH. Further, it may lead to enhanced conformational flexibility of NADH. Unlike other mutations, I21V does not undergo a change in amino acid polarity, however, reduced side chain length due to mutation is supposed to diminish the InhA-NADH interactions and enhance NADH conformational flexibility. Thus, it can be presumed from the discussions that a difference in side chains and polarity of residues upon mutation may lead to the conformational flexibility of the NADH. Conversely, these changes will also enhance the mobility of mutated residues. These conformational changes indicate the weakening of the InhA-NADH interactions and can be measured on the basis of dihedral angle calculations for NADH and mutated residues. In case of I16T, comparison of NADH conformation shows a rmsd difference of 1.28Å from average NADH conformation from WT structure. The major differences are observed at the pyrophosphate region shown in blue colored circle in Fig 5 which makes H-bond interactions with Ser20 and Ile21. In addition to this, conformational differences are also observed for adenine and nicotinamide moiety due to variations in the pyrophosphate conformations. Similarly, NADH conformation in I95P shows a rmsd difference of 1.27Å from NADH conformation in WT. The large differences at pyrophosphate regions are observed in I17T and S94A where a rmsd difference of~1.60 is estimated with respect to NADH conformation in WT. In case of I21V, a minimum rmsd difference of 0.90 is found with respect to NADH conformation in WT which depicts little better NADH binding in I21V as compared to other MTs. It clearly shows that the distribution of dihedral angles D5, D6, D7, D8 and D9 is highly scattered, which suggests sufficient conformational flexibility of the corresponding fragment in NADH i.e. pyrophosphate moiety. Moreover, enhanced flexibility of NADH is expected to diminish InhA-NADH interactions which will also increase flexibility of binding site residues.
Dihedral angle distribution of mutated and active site residues in WT and MTs. The enhanced flexibility of NADH conformation indicates weakening of NADH interactions with active site residues. Therefore, corresponding changes are also expected to occur in the side chain dihedral angle (χ1) of mutated residues involved in interactions with pyrophosphate moiety. The χ1 dihedral angle is estimated for the mutated residues as their side chains directly interact with pyrophosphate moiety of NADH. Fig 6 is shows the dihedral angle distribution of mutated residues involved in interactions with NADH. All the estimated dihedral angles are observed to be highly scattered across a window of 200°to -200°. The side chain distributions of mutated residues show large deviations from the corresponding WT residues. Only, Thr47 show side chain distribution similar to Ile47 throughout the simulation. On the other hand, Thr16 shows side chain distribution similar to Ile16 till 60ns of simulations. Similarly, Val21 also shows side chain distribution corresponding to Ile21 from 20ns to 35ns. Among the mutated residues, Ala94 shows large deviation from the corresponding residue Ser94. The smaller side chain of Ala94 allows more conformational flexibility to side chain due to diminished interactions with NADH. Unlike Ser94, side chain of the mutated residue Pro95 shows little deviation from the Ile95. Thus, it is clear from the discussion that the mutated residues with smaller side chain show deviation in side chain dihedral angle distribution as compared to the corresponding residues in WT. The space created due to side chain difference of WT and MT residues can be accounted for reduced InhA-NADH interactions which increases the conformational flexibility of NADH.
In addition to mutated residues, other binding residues around pyrophosphate moiety of NADH are also expected to face conformational changes due to distorted InhA-NADH interactions due to mutations. The χ1 dihedral angle distribution for the residue Ser20 is found scattered within a window of -100°to 100°which shows inconsistent interactions with NADH ( Fig  7). In contrast to this, dihedral angle distribution for the residues Asp64, Lys165 and Thr196 is observed to be consistent. The χ1 angle distribution for the residue Asp64 is bifurcated into the two rows at about 150°and -150°. However, the H-bond with NADH is still supposed to be maintained as the rotation of side chain by this range of angle will restore the H-bond with another oxygen atom of carboxyl group.
The χ1 dihedral angle for Lys165 is observed consistent within the deviation of~50°for all the mutants except S94A. Similarly, the χ1 dihedral angle for Thr196 is also observed very consistent for the WT and MTs within the deviation of~50°except for I21V and S94A. Moreover, as per the comparison of χ1 among the studied residues, Asp64 is found to be more consistent than the rest of the residues. Overall, differences in dihedral angles of active site residues and H-bond interactions between InhA-NADH in WT and MTs. The conformational changes in the side chain of binding site residues, with progress of the simulations, indicate weakening interactions with NADH. Further, to validate the deduced hypothesis, H-bond occupancy is estimated which is defined as the number of frames with H-bond for particular atom pair to the total number of frame produced in the simulation. H-bond occupancy obtained for residue-NADH interaction with a minimum cut-off of 40% is shown in Table 3 Only, six residues are observed to makes interactions with NADH in the studied complexes as compared to nine residues in crystal structures. The residues Asp64, Ile194 and Thr196 are  Table 3. Summary of H-bond interactions between NADH and InhA.

Residue
Residue atom NADH atom H-bond occupancy (%)   WT  I16T  I21V  I47T  S94A  I95P   Ser20  HG  O2  68  -56  ---Ile21  NH  O9  42  -58  ---Asp64  OD1  H65  -41  46  40  -40   OD2  H65  ------OD1/OD2  H65  49  ---56  -Val65  NH  N4  ---47  -45   Ile194  NH  O1  --40  ---NH  O14  62  62  -51  53  55   CO  H86  19  -32  ---Thr196  HG1  O8  98  89  79  -45  87 doi:10.1371/journal.pone.0144635.t003 observed to make the H-bond interactions with NADH with H-bond occupancy > 40%. However, residues Ser20 and Ile21 are observed to make H-bond interactions with NADH in WT and I21V only. Despite the mutation, Val21 makes H-bond interactions with NADH similar to Ile21. It shows that the change in the side chain of residue due to mutation does not bring change in the interaction as the backbone participates in H-bond interactions in these residues. In contrast to Ser20 and Ile21, Val65 is observed to make H-bond interactions with NADH in I47T and I95P only. In context to high flexibility of pyrophosphate moiety, it is interesting to note that residues Ser20 and Ile21 are observed to maintain H-bond interaction with NADH in WT and I21V. It probably indicates the better stability of the pyrophosphate in WT and I21V as compared to other structures. Further, the detailed analysis of H-bonding within the NADH revealed the two intra-molecular H-bond interactions within NADH in WT and I21V. These interactions are observed between O2 and O5, and O8 and N7 atoms of NADH which can be accounted for the acquired rigidity of the pyrophosphate moiety. However, these residues do not make Hbond interactions with pyrophosphate moiety in I16T, I47T, S94A and I95P. It leads to comparatively better flexibility of NADH conformation in these structures as compared to WT and I21V.
It also indicates the comparatively better binding affinity of NADH in WT and I21V as compared to other complexes. Similarly, average H-bond analysis throughout the trajectory also indicates presence of~5 H-bonds in WT (Fig 8). It is followed by the I21V which shows about 4 H-bonds throughout the simulations. On the other hand, I95P shows about 3 H-bond interactions while I16T and S94A show 3 H-bonds interactions between binding site residues and NADH. Among all studied complexes, Only I47T shows about~2 H-bond interactions with NADH. These results are consistent with H-bond occupancy analysis which also shows better H-bond interactions in WT and I21V as compared to other complexes. Therefore, comparatively good binding energy of NADH is expected in WT. Further, it is expected to continue on with following sequence for NADH binding energy: I21V, I95P, I16T, S94A and I47T.
Water mediated H-bond interactions in WT and MTs. The water mediated H-bonds are also monitored in the WT and MTs throughout the 100ns simulation. Table 4 shows the occupancy of water mediated and direct H-bond to the NADH. The six crystal water molecules (WAT271, WAT272, WAT273, WAT274, WAT275 and WAT283) are observed to interact with NADH in WT and MTs. It is interesting to note that only the WAT271 is observed to maintain direct H-bond interactions in the WT and MTs except the I47T.
In addition to this, WAT271 also maintains the water mediated H-bond interactions through residues Gly14, Ala22 and Ser94 in WT and I21V. However, it shows interactions with residues Gly14 and Ala22 in I16T, S94A and I95P. WAT272 and WAT274 are also observed to make the H-bond with NADH in WT, S94A and I95P, and WT, I16T and S94A respectively. It is interesting to note that most of the water molecules are observed with inconsistent interactions with NADH. Moreover, with the progress of the simulations to 100ns, some new water molecules like WAT273, WAT274, WAT283 and WAT6120 are observed in the DBS. These water molecules are also found to have water mediated interactions with NADH. The analysis of the water mediated H-bond interactions demonstrates the significance of water molecules in stabilizing the NADH binding conformation in WT and MTs. It also shows the importance of WAT271 as a crucial water molecule involved in the InhA-NADH interaction with H-bond occupancy !55% in WT and MTs. However, with the progress of simulations, direct and water mediated interactions with NADH are drastically reduced which can be demonstrated by average H-bond analysis between NADH and water molecules.  These observations clearly depict that the overall interactions with NADH are reduced due to mutations. WT is fairly observed to make maximum interactions with NADH including water mediated and direct H-bond interactions. It is followed by I21V which shows about 6 interactions with NADH. Similarly, I95P also show about~6 H-bond interactions with NADH. It is followed by I16T and S94A with 5 interactions while I47T show minimum~3 interactions with NADH. Therefore, binding affinity of NADH towards WT and MTs is also expected in the similar sequence which can be verified by calculating binding free energy.
Binding free energy calculations for NADH. The comparison of the H-bond occupancy generally demonstrates the difference in the percentage of particular H-bond interaction among the trajectories produced in the MD simulations. However, it does not illustrate the stability of these interactions till the end of the simulations . Figs 8 and 9 show that interaction of NADH with residues gradually decreasing till the end of the simulations. Therefore, average of binding free energy for last 20ns of simulations can give a better estimate of binding affinity of NADH towards WT and MTs ( Table 5). As per the earlier discussion on interaction of NADH with residues, WT is observed with higher binding free energy (ΔG total , -68.53Kcal/mol) for NADH. It is followed by I21V (ΔG total , -54.73Kcal/mol). It is consistent with observed H-bond interactions with NADH where WT is found with about~7 H-bond integrations while about 6 H-bond interactions are estimated for I21V. Despite the same number of H-bond interactions, I95P shows comparatively less binding free energy (ΔG total , -50.60 Kcal/mol) than I21V. Extending the comparison to I16T and S94A, only 5 H-bond interactions are observed with NADH which leads to a binding free energy of about -45kcal/mol. I47T shows the lowest binding free energy for NADH (ΔG total , -37.78 Kcal/mol) among the studied complexes which is consistent with observed interactions. Thus, the following order is obtained by comparing the binding free energy for NADH: WT > I21V > I95P > I16T > S94A > I47T. It is consistent with the observed interactions pattern. In addition to this, a high binding free energy for NADH is also in agreement with reported low dissociation constant for NADH in WT (Table 5). In contrast to this, higher dissociation constants are reported for NADH in MTs which further verifies the estimated binding affinity of NADH against MTs. It simply implies that the observed results are reasonably in agreement with the reported experimental results. In addition, it depicts the usefulness of the methodology to find the molecular details for the drug resistance phenomenon. Further, per residue energy decomposition analysis can provide a quantitative difference in NADH binding with respect to mutated residues.
Per residue decomposition of binding free energy. To get a deeper insight into the influence of mutations, interactions are analyzed on the basis of free energy decomposition to find the per residue contribution in the NADH binding ( Table 6). The majority of the active site residues in WT show highest contribution to the NADH binding as compared to MTs. I16T, I21V, S94A and I95P to show the influence of mutation as the binding energies of the mutated residues are observed to decrease by -3.06, -1.56, -2.17 and -2.01kcal/mol respectively. In case of I16T and S94A, a low binding energy contribution by mutated residues can be accounted due to a change in the non-polar residue into polar and vice versa. It is supposed to disturb the NADH binding due to change in the nature of interacting residues. However, in case of I21V, similar interactions are maintained by the mutated residue Val21 and the small difference in energy is accounted due to reduced length of side chain only. Similarly, I95P shows the large difference in the side chain of mutated residue which is accounted for reduced interactions with NADH. However, I47T does not show the direct effect of mutation as compared to the respective residues in the WT as it is located 6Å away from the NADH binding site. Therefore, mutated residue in I47T does not show its effect on NADH binding. Moreover, apart from the mutated residues, it is important to note that the binding energy of most of the residues is also observed to decrease in the MTs except for a few residues. It shows that a single mutation is sufficient to disturb InhA-NADH binding balance which is reflected as reduced binding energy contribution of other WT residues in the MTs.