Structural and Dynamic Insights into the W68L, L85P, and T87A Mutations of Mycobacterium tuberculosis Inducing Resistance to Pyrazinamide

Tuberculosis (TB), the most frequent bacterium-mediated infectious disease caused by Mycobacterium tuberculosis, has been known to infect humans since ancient times. Although TB is common worldwide, the most recent report by the WHO (World Health Organization) listed the three countries of India, China, and Russia with 27%, 14%, and 8% of the global burden of TB, respectively. It has been reported that resistance to TB drugs, particularly by the pncA gene to the pyrazinamide drug due to mutations, significantly affects the effective treatment of TB. Understanding the mechanism of drug resistance using computational methods is of great interest to design effective TB treatment, exploring the structural features with these tools. Thus, keeping in view the importance of these methods, we employed state-of-the-art computational methods to study the mechanism of resistance caused by the W68L, L85P, and T87A mutations recently reported in 2021. We employed a molecular docking approach to predict the binding conformation and studied the dynamic properties of each complex using molecular dynamics simulation approaches. Our analysis revealed that compared to the wildtype, these three mutations altered the binding pattern and reduced the binding affinity. Moreover, the structural dynamic features also showed that these mutations significantly reduced the structural stability and packing, particularly by the W68L and L85P mutations. Moreover, principal component analysis, free energy landscape, and the binding free energy results revealed variation in the protein’s motion and the binding energy. The total binding free energy was for the wildtype −9.61 kcal/mol, W68L −7.57 kcal/mol, L85P −6.99 kcal/mol, and T87A −7.77 kcal/mol. Our findings can help to design a structure-based drug against the MDR (multiple drug-resistant) TB.


Introduction
Tuberculosis (TB), the most frequent bacterium-mediated infectious disease caused by Mycobacterium tuberculosis, has been known to infect humans since ancient times [1]. Although TB is common worldwide, the most recent WHO report listed the three countries of India, China, and Russia, with 27%, 14%, and 8% of the global burden of TB, respectively [2]. Tuberculosis (TB) causes more deaths worldwide than HIV/AIDS, remaining a public health concern. Recently in 2019, the WHO also reported an estimated 10.0 million infections of TB in the worldwide population. Drug resistance remains a major obstacle in the successful treatment of TB. It has been reported that among the reported TB cases, 0.5 million TB patients (78% of those diagnosed) reported multidrug-resistant (MDR) TB [2]. In 2019, 3.3% newly reported patients, and 17.7% of patients being treated were reported to have MDR or (RR) rifampicin-resistanct TB [2].
The process of drug resistance is triggered by a spontaneous mutation in the target genes, including the panD, rpoB gene, or the inhA promoter/katG gene. Among the known

Molecular Docking
Before molecular docking, we used Open Babel software v 8 (University College Cork, Co. Cork, Ireland) to prepare and minimize drug molecules by employing the Universal Force Field (UFF). Then for optimization, the conjugate gradient algorithm was employed. The prepared structure of PZA was docked against PZase by defining the active site residues for grid generation using the AutoDock Software (The Scripps Research Institute, La Jolla, CA, USA) [42]. Exhaustiveness was set high (64) for better accuracy. Finally, the resulting top-scoring conformations were used for simulation purposes using AMBER20 software (University of California, California, CA, USA).

Molecular Dynamics Simulation
Investigation of dynamic features of the wildtype and mutant complexes were explored through molecular dynamics simulation by employing the AMBER20 simulation tool (University of California, California, CA, USA) [43]. For protein parameterization FF19SB, the GAFF (general amber force field) was used for drug topologies generation. The Fe ion parameters were also generated and "optimal" point charge (OPC) waterbox (10 Å) and Na + ions were added to solvate to neutralize each system. A total of 7256 water molecules and 15 sodium ions were added to the 1373 total atoms. Iterative cycles were completed, each step had 6000 cycles employing the steepest descent algorithm and then the conjugate gradient algorithm for another round of 3000 cycles. Each complex was heated to 300 K, with equilibration at constant pressure 1 atm and weak restraint, while in the second round of equilibration, no restraint was added. Lastly, a 100 ns production run was completed. The temperature was controlled with a Langevin thermostat [36]. Long-range electrostatic interactions were treated with a particle mesh Ewald algorithm with a cutoff distance of 10.0 Å [44]. Finally, a SHAKE algorithm was recruited to treat a covalent bond [45]. The GPU (graphical processing unit)-reinforced simulation for each system was completed, while CPTRAJ and CPPTRAJ modules were used to investigate the trajectories [46]. The simulations were performed at the College of Applied Medical Sciences, Qassim University using an HP Z800 Xeon workstation (Hewlett-Packard, California, CA, USA) with 32 GB RAM and Asus extreme gaming 1080 Ti GPU. The simulation of each complex was completed in 3 days and post-simulation analyses were performed in 1 day.

Principal Component Analysis
Using the CPPTRAJ package (University of California, California, CA, USA), the PCA of each complex was performed to map the high amplitude oscillating motions [47]. The covariance matrix was computed and diagonalized to estimate the eigenvectors in each complex and their respective eigenvalues. The whole trajectory of each complex was used to capture high amplitude fluctuations. PC1 and PC2 were computed and plotted against each other to demonstrate the motions.

Free Energy Landscape (FEL)
The free energy landscape (FEL) approach was used to extract the lowest energy state stable conformation shown by the contour. The intermediate conformations are represented by the separating boundaries between the subspaces [48]. The g_sham module integrated in Gromacs was employed to calculate the FEL. The first two PCs were then used to calculate the FEL based on the following equation:

Binding Affinity Calculations
The total binding energy of the wildtype and W68L, L85P, and T87A mutant complexes with the PZA drug was calculated through MMPBSA.py script (University of Florida, Gainesville, FL, USA) [49][50][51][52][53][54]. For free energy computation from each trajectory, 5000 structures were subjected to calculation of TBE based on the following equation: The subsequent equation was used to compute the value of each component:

Results and Discussion
The apo structure of PZase using PDB ID 3PL1 was retrieved from PDB and the mutations, i.e., W68L, L85P, and T87A were generated using the Chimera rotamers tool. Figure 1A represents the 3D structure of wildtype PZase and Fe 2+ binding. The modeled 3D structures of the mutants are shown in Figure 1B-D, where the mutations are encircled in red. Prior to docking conformation predictions, the structures were minimized and any bad contact was addressed. Superimposition of the wildtype on the three mutants revealed variations in the structures. A superimposition of the wildtype on W68L revealed an RMSD (root mean square deviation) of 0.127 Å, the wildtype and L85A demonstrated an RMSD of 0.115 Å, while the T87A reported an RMSD of 0.111 Å respectively. The docking of the PZA drug was against the active site residues His71, Lys96, Ala134, His137, and Cys138 residues. Two essential residues, i.e., His137 and Cys138 are involved in direct hydrogen bonding interactions with PZase. The docking results revealed that PZA binds to the wildtype more robustly than the three mutants. The docking score for the wildtype was −6.07 kcal/mol; for W68L, it was −5.65 kcal/mol; for L85P, it was −5.73 kcal/mol; while for T87A the docking score was −5.72 kcal/mol. Figure 1E-H shows the variations in the binding pattern of PZA to the wildtype and the three mutants. The wildtype established three essential interactions with ile133, Ala134, and Cys138, which have also previously been reported to interact and produce inhibitory properties. On the other hand, the W68P-PZA complex has three hydrogen bonds (Lys96 and Ala102) and three salt bridges (His57 and His71). This shows that the binding of PZA was altered by the leucine substitution at position 68 and lost the essential interactions with Cys138. Moreover, the binding of L85P also revealed significant variations. With three salt bridges (His57, His71, and Ala134) and three hydrogen bonds (Asp8, Lys96, and Ala102) the docking score of L85P remained relatively higher than W68L. Many salt bridges were established by this mutant; however, due to a lower energy contribution the binding was altered. Moreover, the T87A complexes also established several interactions including hydrogen bonds by Asp8 and Ala102, which were actively involved in the hydrogen bonding. The residues His57, His71, and Ala102 also established salt bridges in this complex. Hence, the docking and binding analysis confirmed that these particular substitutions W68L, L85P, and T87A affected the binding of the PZA. Similar findings of altered binding networks by other mutations by different studies have also been reported previously [55][56][57]. Moreover, we also calculated the binding cavity volume and surface area of all the structures including the wildtype and mutants. The results revealed that the wildtype had a volume of 141.229, the W68L mutant had 136.730, the L85P mutant reported 135.669, while the T87A mutant had 138.945. The findings are consistent with the previous reports, which also confirmed a reduction in the binding cavity volume for the drug-resistant mutations in PZA [56,58]. The top scoring conformations shown in Figure 1 were further used to check the dynamic features altered by these mutations and to demonstrate the binding variation using MM/GBSA. due to a lower energy contribution the binding was altered. Moreover, the T87A complexes also established several interactions including hydrogen bonds by Asp8 and Ala102, which were actively involved in the hydrogen bonding. The residues His57, His71, and Ala102 also established salt bridges in this complex. Hence, the docking and binding analysis confirmed that these particular substitutions W68L, L85P, and T87A affected the binding of the PZA. Similar findings of altered binding networks by other mutations by different studies have also been reported previously [55][56][57]. Moreover, we also calculated the binding cavity volume and surface area of all the structures including the wildtype and mutants. The results revealed that the wildtype had a volume of 141.229, the W68L mutant had 136.730, the L85P mutant reported 135.669, while the T87A mutant had 138.945. The findings are consistent with the previous reports, which also confirmed a reduction in the binding cavity volume for the drug-resistant mutations in PZA [56,58]. The top scoring conformations shown in Figure 1 were further used to check the dynamic features altered by these mutations and to demonstrate the binding variation using MM/GBSA.

Stability Evaluation of the Wildtype and Mutant Complexes
Assessment of structural stability in a dynamic environment always reveals important information regarding the binding pattern and impact of different mutations. To demonstrate the impact of the aforementioned mutations on the binding complexes, we also calculated the RMSD for each complex. Our analysis revealed that the wildtype complex remained more stable than the mutant complexes. It can be seen that the RMSD of the wildtype exhibited a rigid behavior and no significant deviation was experienced by the system. The complex demonstrated stable dynamics with an average RMSD of 1.2 Å . Comparative assessment of wildtype and W68L revealed significant differences in the stability of the mutant complex. The W68L complex faced significant structural perturbation during the first 35 ns and, then, although the RMSD increased gradually, no significant

Stability Evaluation of the Wildtype and Mutant Complexes
Assessment of structural stability in a dynamic environment always reveals important information regarding the binding pattern and impact of different mutations. To demonstrate the impact of the aforementioned mutations on the binding complexes, we also calculated the RMSD for each complex. Our analysis revealed that the wildtype complex remained more stable than the mutant complexes. It can be seen that the RMSD of the wildtype exhibited a rigid behavior and no significant deviation was experienced by the system. The complex demonstrated stable dynamics with an average RMSD of 1.2 Å. Comparative assessment of wildtype and W68L revealed significant differences in the stability of the mutant complex. The W68L complex faced significant structural perturbation during the first 35 ns and, then, although the RMSD increased gradually, no significant deviation was experienced. The average RMSD for W68L was reported to be 2.5 Å. Moreover, the L85P complex initially demonstrated rigid and stable dynamics until 40 ns but then the RMSD deviated by a larger proportion and significant convergence was experienced by this complex. The average RMSD for L85P was reported to be 2.4 Å. On the other hand, the T87A and wildtype RMSD remained comparable. The complex did not converge, and the average RMSD was reported to be 1.1 Å. The findings are consistent with previous studies, which also reported that drug resistant mutations destabilized the structures in a dynamic environment, thus implying that these mutations produce destabilizing effects and, consequently, pose resistance to PZA [55][56][57]. The RMSD graphs of the wildtype and mutant complexes are shown in Figure 2. We also calculated the RMSDs of the ligands only, which confirmed the stability of the ligands during the simulation. Binding and unbinding of the ligands during the simulation resulted in an increase or decrease in the RMSD, particularly in the mutant complexes, as can be seen at different time intervals. The RMSDs of the ligands are shown in Figure 3. Moreover, the secondary structural elements were also calculated, which revealed further variations induced by the mutations. In the case of the wildtype, 32.43% of the residues were alpha-helix, 19.46% extended strand, and 8.11% beta sheets, while 40% were random coiled. In the case of the W68L mutant, 34.59% of the residues were alpha-helix, 17.30% extended strand, and 7.57% beta sheets, while 40.54% were random coiled. Consequently, the percentage of beta sheets was reduced, which is directly associated with protein stability. In the case of L85P and T87A, 35.14% and 37.84% of the residues were alpha-helix, respectively, 17.30% and 18.38% were extended strand, and 7.57% and 17.30% were beta sheets, while 38.92% and 34.93% were random coiled. It can be seen that the number of beta sheets was increased in the T87A complex, which in turn increased the stability, as can be seen in the RMSD plot, where the RMSD graph is comparable with the wildtype. deviation was experienced. The average RMSD for W68L was reported to be 2.5 Å . Moreover, the L85P complex initially demonstrated rigid and stable dynamics until 40 ns but then the RMSD deviated by a larger proportion and significant convergence was experienced by this complex. The average RMSD for L85P was reported to be 2.4 Å . On the other hand, the T87A and wildtype RMSD remained comparable. The complex did not converge, and the average RMSD was reported to be 1.1 Å . The findings are consistent with previous studies, which also reported that drug resistant mutations destabilized the structures in a dynamic environment, thus implying that these mutations produce destabilizing effects and, consequently, pose resistance to PZA [55][56][57]. The RMSD graphs of the wildtype and mutant complexes are shown in Figure 2. We also calculated the RMSDs of the ligands only, which confirmed the stability of the ligands during the simulation. Binding and unbinding of the ligands during the simulation resulted in an increase or decrease in the RMSD, particularly in the mutant complexes, as can be seen at different time intervals. The RMSDs of the ligands are shown in Figure 3. Moreover, the secondary structural elements were also calculated, which revealed further variations induced by the mutations. In the case of the wildtype, 32.43% of the residues were alpha-helix, 19.46% extended strand, and 8.11% beta sheets, while 40% were random coiled. In the case of the W68L mutant, 34.59% of the residues were alpha-helix, 17.30% extended strand, and 7.57% beta sheets, while 40.54% were random coiled. Consequently, the percentage of beta sheets was reduced, which is directly associated with protein stability. In the case of L85P and T87A, 35.14% and 37.84% of the residues were alpha-helix, respectively, 17.30% and 18.38% were extended strand, and 7.57% and 17.30% were beta sheets, while 38.92% and 34.93% were random coiled. It can be seen that the number of beta sheets was increased in the T87A complex, which in turn increased the stability, as can be seen in the RMSD plot, where the RMSD graph is comparable with the wildtype.

Determination of Structural Compactness as Rg (Radius of Gyration)
Determination of structural compactness in a dynamic environment demonstrates important binding and unbinding events and the impact of different mutations [59,60]. To demonstrate the binding and unbinding events, we calculated the radius of gyration (Rg) for each complex as a function of time. The Rg results of the wildtype and mutants are strongly in uniformity with the RMSD results. The Rg(s) of the wildtype and T87A demonstrated similar structural compactness. The average Rg for the wildtype and T87A was reported to be 15.52 Å and 15.50 Å, respectively. Moreover, the L85P and W68L experienced significant variations and perturbations in the Rg during the simulation. The Rg of the W68L gradually increased over the simulation time and reached 16.0 Å. During the first 35 ns, the average Rg value was 15.8 Å and then increased during the last 65 ns up to 16.1 Å. The increment in structural compactness after 90 ns was due to the reduction in the motion of the lid and flap region of the protein structure, which stabilized at that particular time point, hence experiencing reduction in the structural compactness. On the other hand, the Rg pattern of T87A also revealed a consistent graph as did the RMSD. Initially, the structure remained well packed until 45 ns, with a minor deviation, and then, a significant deviation until 100 ns was reported. The average Rg for the T87A complex was reported to be 15.8 Å. The Rg graphs of the wildtype and the mutant complexes are shown in Figure 4.

Determination of Structural Compactness as Rg (Radius of Gyration)
Determination of structural compactness in a dynamic environment demonstrates important binding and unbinding events and the impact of different mutations [59,60]. To demonstrate the binding and unbinding events, we calculated the radius of gyration (Rg) for each complex as a function of time. The Rg results of the wildtype and mutants are strongly in uniformity with the RMSD results. The Rg(s) of the wildtype and T87A demonstrated similar structural compactness. The average Rg for the wildtype and T87A was reported to be 15.52 Å and 15.50 Å , respectively. Moreover, the L85P and W68L experienced significant variations and perturbations in the Rg during the simulation. The Rg of the W68L gradually increased over the simulation time and reached 16.0 Å . During the first 35 ns, the average Rg value was 15.8 Å and then increased during the last 65 ns up to 16.1 Å . The increment in structural compactness after 90 ns was due to the reduction in the motion of the lid and flap region of the protein structure, which stabilized at that particular time point, hence experiencing reduction in the structural compactness. On the other hand, the Rg pattern of T87A also revealed a consistent graph as did the RMSD. Initially, the structure remained well packed until 45 ns, with a minor deviation, and then, a significant deviation until 100 ns was reported. The average Rg for the T87A complex was reported to be 15.8 Å . The Rg graphs of the wildtype and the mutant complexes are shown in Figure 4.

Residual Flexibility (RMSF) Estimation
Estimation of residual flexibility is an important parameter to reveal the conformational flexibility and strength conferred by a particular residue in the binding [61]. Herein, to calculate the residual flexibility we used the root mean square fluctuation (RMSF) approach. The results, as shown in Figure 4, revealed that the flexibility of the wildtype and

Residual Flexibility (RMSF) Estimation
Estimation of residual flexibility is an important parameter to reveal the conformational flexibility and strength conferred by a particular residue in the binding [61]. Herein, to calculate the residual flexibility we used the root mean square fluctuation (RMSF) approach. The results, as shown in Figure 4, revealed that the flexibility of the wildtype and mutant complexes remained comparable; however, variations in different regions particularly 15−25, 30-45, 55-75, and 90-105 demonstrated higher fluctuation in the wildtype and W68L, specifically. The other variants demonstrated higher fluctuation between 30 and 45 aa. The mutation W68L in the loop region resulted in higher fluctuation of the loop, thus rendering higher fluctuation in this particular mutation. In conclusion, the mutations affected the residual flexibility of the protein structure differently and caused functional variations. The RMSF graph of all the complexes is shown in Figure 5.

Residual Flexibility (RMSF) Estimation
Estimation of residual flexibility is an important parameter to reveal the conformational flexibility and strength conferred by a particular residue in the binding [61]. Herein, to calculate the residual flexibility we used the root mean square fluctuation (RMSF) approach. The results, as shown in Figure 4, revealed that the flexibility of the wildtype and mutant complexes remained comparable; however, variations in different regions particularly 15−25, 30-45, 55-75, and 90-105 demonstrated higher fluctuation in the wildtype and W68L, specifically. The other variants demonstrated higher fluctuation between 30 and 45 aa. The mutation W68L in the loop region resulted in higher fluctuation of the loop, thus rendering higher fluctuation in this particular mutation. In conclusion, the mutations affected the residual flexibility of the protein structure differently and caused functional variations. The RMSF graph of all the complexes is shown in Figure 5.

Hydrogen Bonding Analysis
Analysis of the hydrogen bonds of the protein-ligand trajectories is essential to understand the binding variations produced by epigenetics or mutations [62]. To quantify the total number of hydrogen bonds and reveal the average, we subjected the simulation

Hydrogen Bonding Analysis
Analysis of the hydrogen bonds of the protein-ligand trajectories is essential to understand the binding variations produced by epigenetics or mutations [62]. To quantify the total number of hydrogen bonds and reveal the average, we subjected the simulation trajectories to hydrogen bonding calculations as a function of time. Our results revealed that the wildtype and W68L, L85P, and T87A had a different total number of hydrogen bonds, as shown in Figure 5. In each complex, the average hydrogen bonds were estimated to be 93 (wildtype), 90 (W68L), 82 (L85P), and 88 in the T87A mutant. This shows that the mutations destabilized the binding and, thus, caused resistance to PZA. Similar findings were also reported by previous studies, thus justifying our results [55][56][57]. The hydrogen bonds graph of all the complexes is shown in Figure 6.

Principle Component Analysis (PCA) and Free Energy Landscape (FEL)
To identify the collective motions of the wildtype and mutant complexes, we explored the dynamic behavior using PCA. It is a statistical method for dimensionality reduction without compromising vital information [63]. Plentiful fluctuations were revealed by the first three eigenvectors, while localized fluctuations were demonstrated by the other eigenvectors in each complex, as shown in Figure 6. In the case of the wildtype-PZA complex, the first 3 eigenvectors were 53%, the W68L-PZA complex demonstrated 63%, L85P displayed 44%, and T87A 57% of the overall observed motion. This behavior may explain the structural rearrangement due to the mutations; thus, it may empirically be proposed that the PZA interaction may destabilize the protein by increasing the dynamics of the active regions to a higher level. The current findings also support the previous reports with higher magnitude of motion in the PZA complexes and directions identical to the wildtype [55][56][57][58]. Moreover, we also explored the principal components, i.e., PC1 and PC2 (principle components) of the simulation trajectories to project each frame in a phase space (Figure 7). The interpretation of the pink to purple displays the conformational transformations during the simulation. These conformational states were found closer to each other in the three mutants and were considered as an energetically unstable conformational state by attaining an unstable conformational state purple color.
trajectories to hydrogen bonding calculations as a function of time. Our results revealed that the wildtype and W68L, L85P, and T87A had a different total number of hydrogen bonds, as shown in Figure 5. In each complex, the average hydrogen bonds were estimated to be 93 (wildtype), 90 (W68L), 82 (L85P), and 88 in the T87A mutant. This shows that the mutations destabilized the binding and, thus, caused resistance to PZA. Similar findings were also reported by previous studies, thus justifying our results [55][56][57]. The hydrogen bonds graph of all the complexes is shown in Figure 6. Figure 6. Hydrogen bonding analysis of the wildtype and W68L, L85P, and T87A mutants. The wild type structure is shown as black, W68L (red), L85P (magenta), and T87A (purple).

Principle Component Analysis (PCA) and Free Energy Landscape (FEL)
To identify the collective motions of the wildtype and mutant complexes, we explored the dynamic behavior using PCA. It is a statistical method for dimensionality reduction without compromising vital information [63]. Plentiful fluctuations were revealed by the first three eigenvectors, while localized fluctuations were demonstrated by the other eigenvectors in each complex, as shown in Figure 6. In the case of the wildtype-PZA complex, the first 3 eigenvectors were 53%, the W68L-PZA complex demonstrated 63%, L85P displayed 44%, and T87A 57% of the overall observed motion. This behavior may explain the structural rearrangement due to the mutations; thus, it may empirically be proposed that the PZA interaction may destabilize the protein by increasing the dynamics of the active regions to a higher level. The current findings also support the previous reports with higher magnitude of motion in the PZA complexes and directions identical to the wildtype [55][56][57][58]. Moreover, we also explored the principal components, i.e., PC1 and PC2 (principle components) of the simulation trajectories to project each frame in a phase space (Figure 7). The interpretation of the pink to purple displays the conformational transformations during the simulation. These conformational states were found closer to each other in the three mutants and were considered as an energetically unstable conformational state by attaining an unstable conformational state purple color.

Free Energy Landscape (FEL) Analysis
In addition, a free energy landscape (FEL) was created to link structural characteristics to thermodynamic attributes. FEL was used to map the minimum energy conformation of the complexes within the examined time scale, and then to connect the structural shifts between these minima based on the likelihood of provided data points of the MD trajectories. Figure 8 presents the FEL of the wildtype, W68L, L85P, and T87A where all reached only one energy minima, thus explaining global conformational variations adjusted by the mutant complexes in response to mutations.

Free Energy Landscape (FEL) Analysis
In addition, a free energy landscape (FEL) was created to link structural characteristics to thermodynamic attributes. FEL was used to map the minimum energy conformation of the complexes within the examined time scale, and then to connect the structural shifts between these minima based on the likelihood of provided data points of the MD trajectories. Figure 8 presents the FEL of the wildtype, W68L, L85P, and T87A where all reached only one energy minima, thus explaining global conformational variations adjusted by the mutant complexes in response to mutations.

Free Energy Landscape (FEL) Analysis
In addition, a free energy landscape (FEL) was created to link structural characteristics to thermodynamic attributes. FEL was used to map the minimum energy conformation of the complexes within the examined time scale, and then to connect the structural shifts between these minima based on the likelihood of provided data points of the MD trajectories. Figure 8 presents the FEL of the wildtype, W68L, L85P, and T87A where all reached only one energy minima, thus explaining global conformational variations adjusted by the mutant complexes in response to mutations.

Free Energy Calculations
Using an MM/GBSA (molecular mechanics generalized Born and surface area continuum solvation) approach to re-evaluate the TBE of a protein and ligand is a common approach [1,2,64]. The TBE (total binding energy) was calculated for wildtype, W68L, L85P, and T87A-PZA complexes. The ΔGbind of wild type, W68L, L85P, and T87A-PZA complexes was reported to be −9.61 kcal/mol, −7.57 kcal/mol, −6.99 kcal/mol, and −7.77 kcal/mol, respectively (Table 1). When comparing the overall energies of mutants to the WT, it is clear that these mutations reduce the PZA's binding strength. In comparison to the WT, the contributions of vdW (Van der Waal), Electrostatic, and PS (polar solvation) energies to the binding energies of the mutants were substantially lower. It was discovered that the mutant proteins have a low affinity for PZA. The configuration of the active

Free Energy Calculations
Using an MM/GBSA (molecular mechanics generalized Born and surface area continuum solvation) approach to re-evaluate the TBE of a protein and ligand is a common approach [1,2,64]. The TBE (total binding energy) was calculated for wildtype, W68L, L85P, and T87A-PZA complexes. The ∆G bind of wild type, W68L, L85P, and T87A-PZA complexes was reported to be −9.61 kcal/mol, −7.57 kcal/mol, −6.99 kcal/mol, and −7.77 kcal/mol, respectively (Table 1). When comparing the overall energies of mutants to the WT, it is clear that these mutations reduce the PZA's binding strength. In comparison to the WT, the contributions of vdW (Van der Waal), Electrostatic, and PS (polar solvation) energies to the binding energies of the mutants were substantially lower. It was discovered that the mutant proteins have a low affinity for PZA. The configuration of the active site residues in direct contact with the PZA is affected by mutations that are not engaged in direct interaction with the PZA. The current findings strongly corroborate the previous findings where reduction in the binding free energy was reported to be induced by different mutations [55][56][57][58].

Conclusions
Forecasting the effect of mutations on the structure of a protein's function and binding offers prodigious potential for formulating therapeutic approaches. A thorough knowledge of biological processes, the relationship between genetic patterns and phenotypic characteristics, and how these systems interact with the host can aid in the development of new and successful treatments. As proteins and their biological functionality are controlled by genes, it is essential to understand that both genomics and proteomics patterns hold significant importance in therapeutics research. Thus, keeping in view the importance of these methods, we have employed advanced computational methods to study the mechanism of resistance caused by the W68L, L85P, and T87A mutations recently reported in 2021. We employed a molecular docking approach to predict the binding conformation and studied the dynamic properties of each complex using molecular dynamics simulation approaches. Our analysis revealed that compared to the wildtype, these three mutations altered the binding pattern and reduced the binding affinity. Moreover, principal component analysis, free energy landscape, and the binding free energy results revealed variation in the protein's motion and the binding energy. The total binding free energy for the wildtype was −9.61 kcal/mol, W68L −7.57 kcal/mol, L85P −6.99 kcal/mol, and T87A −7.77 kcal/mol. Our findings could help to design structure-based drugs against the MDR TB. Using these structural features, drugs which can induce structural stability and particularly the Fe ion stability could effectively overcome the drug resistance problem. This study offers a profound understanding on the mechanism of resistance caused by the reported mutations in tuberculosis.

Data Availability Statement:
The data presented in this study are available within the article.