Insights Into Mutations Induced Conformational Changes and Rearrangement of Fe2+ Ion in pncA Gene of Mycobacterium tuberculosis to Decipher the Mechanism of Resistance to Pyrazinamide

Pyrazinamide (PZA) is the first-line drug commonly used in treating Mycobacterium tuberculosis (Mtb) infections and reduces treatment time by 33%. This prodrug is activated and converted to an active form, Pyrazinoic acid (POA), by Pyrazinamidase (PZase) enzyme. Mtb resistance to PZA is the outcome of mutations frequently reported in pncA, rpsA, and panD genes. Among the mentioned genes, pncA mutations contribute to 72–99% of the total resistance to PZA. Thus, considering the vital importance of this gene in PZA resistance, its frequent mutations (D49N, Y64S, W68G, and F94A) were investigated through in-depth computational techniques to put conclusions that might be useful for new scaffolds design or structure optimization to improve the efficacy of the available drugs. Mutants and wild type PZase were used in extensive and long-run molecular dynamics simulations in triplicate to disclose the resistance mechanism induced by the above-mentioned point mutations. Our analysis suggests that these mutations alter the internal dynamics of PZase and hinder the correct orientation of PZA to the enzyme. Consequently, the PZA has a low binding energy score with the mutants compared with the wild type PZase. These mutations were also reported to affect the binding of Fe²⁺ ion and its coordinated residues. Conformational dynamics also revealed that β-strand two is flipped, which is significant in Fe²⁺ binding. MM-GBSA analysis confirmed that these mutations significantly decreased the binding of PZA. In conclusion, these mutations cause conformation alterations and deformities that lead to PZA resistance.

Introduction

Pyrazinamide (PZA), along with isoniazid (INH) and rifampin (RIF) is a very effective and fast therapy against persistent bacilli (Mitchison, 1985; Aggarwal et al., 2018). Pyrazinamide (PZase) encoded by the pncA gene of Mycobacterium tuberculosis (Mtb) transform this prodrug to pyrazinoic acid (POA). POA inhibits the proliferation of latent Mtb at very low pH values (Zhang and Mitchison, 2003; Malik et al., 2019). Studies have shown that resistance is developed against PZA due to mutations in three genes: pncA, panD, and rpsA, among which pncA gene mutations contribute to 72–99% resistance against PZA (Mitchison, 1985; Akhmetova et al., 2015; Miotto et al., 2017).

Mutations in the pncA gene have been mapped both in the coding as well as promoter region (Lemaitre et al., 2001; Miotto et al., 2014; Maningi et al., 2015). Recent investigations indicated that Pzase activity is affected due to mutations in D49A, Y64S, W68G and F94A positions (Miotto et al., 2014). The mentioned mutations have been shown to affect enzyme functionality drastically, and together with other reported mutations, influence protein structure integrity, solubility, function stability, and rate of expression (Petrella et al., 2011). More recently, novel pncA mutations are being described as liable to cause PZA resistance (Tan et al., 2014; Junaid et al., 2018).

The crystallographic structure of apo pyrazinamidase has been reported comprising six β-sheets covered by α-helices. This enzyme has metal and substrate binding sites. Iron (Fe²⁺ ion), histidine (His51, His57, and His71), and aspartate (Asp49) residues are part of the metal-binding site, whereas Asp8, Lys96, and Cys138 make the catalytic triads (Chaturvedi and Shrivastava, 2005; Petrella et al., 2011).

Computational approaches are now in routine to decipher mutations mediated biological mechanisms responsible for neutralizing the action of potent drugs. This atomic-level understanding holds great potential in de nova drug design and as such, speeds up novel drug discovery. In particular, advancements in molecular dynamics simulations allows scientist to analyze protein dynamics in environmental milieu replica of real biological cells (Khan et al., 2018a,b, 2019a, 2020a). It has been noticed that binding of PZA to the PZase enzyme altered protein’s conformation, which is valuable data-keeping their importance in the quest of novel drug design. Likewise, MD simulations made it possible to study conformational variations in the three-dimensional structure of proteins that may arise following mutation(s) in the sequence (Khan et al., 2019c, 2020b,f). Thus MD simulations decrease time, costs and resources by reducing the number of cases for which experimental evaluation is required (Dolatkhah et al., 2017). We also investigated the molecular mechanism behind the resistance caused by D49N, Y64S, W68G, and F94A mutations (Stoffels et al., 2012; Miotto et al., 2014; Wan et al., 2020). Furthermore, extensive post-simulation analyses were employed to get insights into the atomic level with an ultimate objective to design novel chemical structures that can be effectively used in drug-resistant TB infections—with minimum side effects.

Materials and Methods

PZase and PZA Structure Retrieval

The 3D structure of Mtb PZase (accession ID: 3PLI) and PZA (accession ID CID1046) were retrieved from the PDB databank (Rose et al., 2016) and PubChem (Kim et al., 2019) respectively. Water molecules were removed from the protein structure before starting downward analyses. As specific mutant structures of the enzyme were not available, mutations were introduced in the enzyme structure using PYMOL (DeLano, 2002) at particular locations.

Molecular Docking

Energy minimization steps were performed for PZA structure in Open Babel using Universal Force Field (Dallakyan and Olson, 2015). The ligand was optimized with default steepest descent and conjugate gradient algorithms in UCSF Chimera (Goddard et al., 2005). Docking was done in the PatchDock server, where binding conformation clusters were set at RSMD of 4.0 Å (Schneidman-Duhovny et al., 2005). Conformations with the lowest binding score were processed for molecular dynamics simulation using AMBER18 software (Wang et al., 2001; Case et al., 2005).

Impact of Mutations on Protein-Drug Interaction and Stability

Mutations’ effect on protein thermodynamic stability was evaluated using mCSM¹ (Pires et al., 2013). The server utilizes graph-based signatures to predicts structural stability impact caused by mutation. mCSM accepts PDB files as input and a list of mutations to predict their effect on protein stability.

Molecular Dynamics Simulation

AMBER18 package was used to perform extensive MD simulations. This was done to investigate the stability of the PZA at the active site of both normal and mutant PZase. Parameters of protein were generated through ff14SB force field, and ligand preparation was done via Amber general force field (GAFF) (Wang et al., 2001; Case et al., 2005). MD simulations were performed for all five systems, including one wild (WT) and four mutants (D49N, Y64S, W68G, and F94A). Each system is solvated in TIP3P water box. Counter ions were added to each system to get charge neutralization. Afterward, two-step energy minimization procedure was adopted; (i) steepest decent minimizations of 6,000 cycles and (ii) conjugate gradient minimization of 3,000 cycles was applied on each system to remove steric clashes and allow system relaxation. Complexes were then heated to 300 K for 0.2 ns, followed by systems equilibration for 2 ns at 300 K. Temperature hold was achieved via Langevin thermostat (Zwanzig, 1973). For all systems, MD simulations production run was completed on GPU supported PMEMD code for 100 ns, and each simulation was repeated three times. Long-range electrostatic interactions (Darden et al., 1993; Essmann et al., 1995; Toukmaji et al., 2000) were detected with the particle mesh Ewald method using a cutoff distance of 10.0 Å. SHAKE method was applied for covalent bond treatment (Kräutler et al., 2001) (Salomon-Ferrer et al., 2013). CPPTRAJ and PTRAJ (Roe and Cheatham, 2013) packages in AMBER18 were considered for trajectories analysis.

Principal Component Analysis

Principal component analysis was utilized to measure structural fluctuations within the protein of all used complexes (Amadei et al., 1993). CPPTRAJ package calculated the covariance matrix based on Cα coordinates. Eigenvectors and eigenvalues estimation was performed by diagonalizing the covariance matrix, and these values indicate motion direction and fluctuation, respectively. In total, 5000 frames from each system MD trajectories were used to get PCA calculations. The plotting performed on PC1 and PC2 was used for motion monitoring. The lowest energy stable state was determined by the free energy landscape (FEL) and is indicated by deep valleys on the plot, whereas the intermediate state is shown by boundaries between deep valleys (Hoang et al., 2004). In this study, FEL calculations based on PCI and PC2 were obtained by the following equation:

$\mathrm{\Delta }G(\mathrm{PC1},\mathrm{PC2})=-K\text{B}T\text{ln}P(\mathrm{PC1},\mathrm{PC2})$

Where KB indicates Boltzmann constant, PC1 and PC2 were used to estimate the reaction coordinates, and probability distribution P of the system is shown along PC1 and PC2.

Binding Affinity Estimation

PZA binding free energy with PZase (native and mutants) was estimated through MMPBSA.py script of AMBER over 500 snapshots of simulation trajectories (Miller et al., 2012; Mishra and Koča, 2018). The equation given below is used for binding free energy calculations

$\mathrm{\Delta }{G}_{\text{bind}}=\mathrm{\Delta }{G}_{\text{complex}}-[\mathrm{\Delta }{G}_{\text{receptor}}+\mathrm{\Delta }{G}_{\text{ligand}}]$

where ΔG_bind, ΔGcomplex, ΔGreceptor, and ΔGligand indicate net binding free energy, binding free energy of the complex, protein, and ligand, respectively. The following equation was used to calculate the value of each component:

$G=G_{\text{bond}}+G_{\text{ele}}+G_{\text{vdW}}+G_{\text{pol}}+G_{\text{npol}}$

where the energy of bonds, electrostatic, van der Waals interactions, the polar and non-polar contributions are shown by the Gbond, Gele, GvdW, Gpol, and Gnpol, respectively. Whereas Gpol and Gnpol were calculated by the generalized Born (GB) implicit solvent method with SASA.

Results

Mutant PZase Structural Modeling and Docking With PZA

The PZase apo structure (available as crystal structure) with ID: 3PL1 was retrieved from the protein databank and subjected to mutagenesis module in the PyMOL software where D49N, Y64S, W68G, and F94A mutants were created. Before molecular docking, all the structures were minimized by removing bad contacts from newly mutated residues as well as other residues. Following the minimization process, the docking process was completed blindly. Docking results suggested that our docking protocol is reliable, as indicated by the involvement of similar residues in interaction, as reported by a previous study (Junaid et al., 2018, 2020). Two residues such as His137 and Cys138, were reported to be involved in hydrogen bonding interactions with the oxygen of PZA. In the present study, similar results were obtained. The docking score of all complexes, including wild type and mutants, are tabulated in Table 1. The more negative binding energy implies better PZase-PZA intermolecular complementarity and higher binding affinity in contrast to the positive binding energy. The complex structure of wild type PZase and the PZA and its interaction pattern are given in Figure 1A. The docking score of PZA with both wild type and mutants PZase is in the following order: wild type (−5.21 kcal/mol), D49N (−4.75 kcal/mol), Y64S (−4.1 kcal/mol), W68G (−4.51 kcal/mol) and F94A (−4.18 kcal/mol). This data suggests that the PZA drug has a higher binding affinity for the wild type PZase enzyme in contrast to the mutants. Among the mutants, the lowest binding affinity of the PZA drug was noticed for Y64S and F94A. There is a high possibility that the mutations alter the active pocket conformation and thus not allowing proper PZA binding. The binding interaction pattern of each complex is given in Figures 1–E. MD simulations were performed on top scorer conformations to ascertain the effect of the mutation on the PZase structure as well as its binding with PZA.

TABLE 1

Table 1. Molecular docking scores of the wild and mutant complexes. The mCSM predicted stability changes upon mutation. All the energies are given in kcal/mol.

FIGURE 1

Figure 1. 3D structure of the PZase along with the PZA drug and the Fe²⁺ metal shown in the circle. The figure also shows the binding of Fe²⁺ ion and PZA drug to the wild type (A), D49N (B), Y64S (C), W68G (D), and F94S (E).

Dynamics Characterization of Wild Type and Mutant Complexes

Mutations were found to confer instability in enzyme structure as predicted by mCSM web server and also by RMSD plots from a triplicate run of 100 ns MD simulations. mCSM server predicts the impact of each substitution by forecasting the change in conformational energy. As given in Table 1, it can be easily pointed that the given mutations induced greater instability compared to the wild type and hence classified as highly destabilizing. Among the four mutants, it was observed that W68G has a profound destabilizing effect on the PZase enzyme with ΔG of −3.14 kcal/mol. This was followed by F94A mutation that contributes to enzyme destabilizing change of −2.94 kcal/mol (Miotto et al., 2017). Among others, the predicted destabilizing energy change for Y64S is −2.2 kcal/mol whereas, for D49N, the energy change is 2.0 kcal/mol. These findings are in line with the docking score of the systems and together, both analysis demonstrated the mutations are responsible for the change of PZase active pocket conformation, thus destabilize the binding network of the PZA drug, as can be seen in Figure 1.

For the stability assessment of each system, Cα atoms root-mean-square deviation (RMSD) was calculated based on simulated trajectory. The WT system reached an equilibrium state up to 60 ns, followed by a minor RMSD increased up to a maximum of 1.5 Å. Later, the RMSD continued over 1.5 Å with insignificant fluctuation (Figure 2). The D49N mutant system reached an equilibrium state of 2 Å in the first 20 ns, and then the RMSD fluctuated high throughout the simulation time due to system instability compared to WT. The Y64S system, like the WT gained equilibrium in the 50 ns and remained stable with slight fluctuations in the RMSD. The W68G system is in stable conformation till 30 ns with RMSD of 2.2 Å, then retained with RMSD at 1.5 Å and fluctuating slightly from the WT for the rest of time. The F94A system gains equilibrium in the first 10 ns and afterward showing minor fluctuations up to 2 Å. This unstable dynamics behavior of the mutants supports the enzyme conformation changes upon mutations to show resistance against PZA. Further inspection of Cα-RMSD rise for mutants compared to the WT showed that the D49N, Y64S, W68G, and F94A might weaken the active site residues interactions with the PZA. The RMSD of the mutant complexes is comparable with the wild type in terms of RMSD value, but the destability justifies that the different convergences at different intervals faced by the mutant structures but not in the wild type. This explanation of the wild and mutant complexes elucidates that due to small protein, the systems have reached the equilibrium point earlier. Furthermore, it can also be seen that the wild type reached the stability at 1.0Å; however, the other systems gained the equilibrium at ∼1.5Å, which shows the mutations induced structural perturbation in mutant complexes. The RMSD results for the other two replicates are given in Supplementary Figures 1, 2.

FIGURE 2

Figure 2. RMSD of wild and mutants’ complexes. RMSD of each mutant is superimposed on to the RMSD of the wild type. The X-axis shows the simulation time in nanoseconds, while the y-axis shows the RMSD in Angstrom.

Local fluctuations due to mutations were examined through Cα, root-mean-square fluctuation (RMSF). Residues fluctuation was noted significantly in the mutant systems compared to the WT. WT system fluctuates at the N-terminus. The D49N mutant system reveals several point fluctuations as compared to WT and other mutations. The RMSF high fluctuation from the WT discloses that the mutations greatly affect the binding of the drug to the active site of the protein. The flexibility of the mutants may justify the binding differences, which can be better revealed by exploring the binding affinity differences. In the case of the mutants, the specific fluctuations at the site of the mutation can be easily distinguished. The RMSF of all the systems is given in Figure 3. The RMSF results for the other two replicates are shown in Supplementary Figures 3, 4.

FIGURE 3

Figure 3. RMSF of the wild and mutants’ complexes. The RMSD of each mutant is superimposed over RMSF of the wild type. The x-axis shows residues number, while the y-axis shows RMSF in Angstrom. Shadowed regions depict enzyme amino acids stretch highly affected by the mutation.

Figure 4 presents broader distance distributions in mutants in contrast to WT, indicating more conformation dynamics in the former systems. As three residues: His51, His57, and His71 form a catalytic triad, it is important to understand the effect of these substitutions on the triad dynamics. It can be seen that the wild type, Y64S, and F94A showed a similar pattern of dynamics, while the D49N and W68G possess different triad distance network dynamics. Distinct changes of His57 is due to the loop harboring this residue. Fe²⁺ ion disturbance may reduce PZase activity and may explain the resistance phenomenon of these mutations. This effect was also confirmed by calculating the distance between PZA and the PZase. Supplementary Figure 5 shows that the distance between the wild type and the PZA is conserved, and the average distance reported was 8Å. However, this distance significantly fluctuates in the case of D49N, W68G, and F94A. While in the case of Y64S, the distance between the PZA and the receptor molecule remained somewhat similar to the wild type. Thus, these results also confirm that mutations have induced structural destabilization and favor PZA unbinding due to their weak attachment.

FIGURE 4

Figure 4. Distances between the Fe²⁺ ion and its coordinating four residues. Furthermore, within each figure inside, there is a legend that shows the distance between Asp49 [O], His51[N], His57[N], His71[N] and Fe²⁺ ion. Each residue from the metal coordinates is differently colored.

Furthermore, we also calculated Rg to estimate the compactness of each system. The calculated Rg for each complex is given in Supplementary Figures 6–8. The results show that D49N is less compact than the wild type. For D49N initially, the higher Rg was observed, which then decreased; however, similar pattern of increasing and decreasing was experienced until 100 ns. In the case of Y64S, the Rg pattern was comparable with the wild type, but at 40–60 ns the Rg converged and a similar pattern was also observed between 95–100 ns.

Similarly, W68G systems were significantly affected. The Rg value significantly increased, and the average Rg was reported to be 15.6Å. The results of F94A and Y64S are comparable. No significant convergence was observed; however, at different intervals, the Rg increased.

Dimensionality Reduction and Clustering the Protein Motions

To understand the protein motion and cluster the related structural frames, PCA was performed. PCA is a mathematical method that transforms several correlated variables into smaller uncorrelated variables called principal components. To comprehensively understand the impact of the substitution on the protein motion initially, the eigenvectors were calculated and presented in Figure 5.

FIGURE 5

Figure 5. Fractions of the first ten eigenvectors. Using the MD trajectory, the fraction of motions is calculated and given in percentage against the eigenvector numbers.

As given in Figure 5, the first three eigenvectors showed significant variations while rest of the eigenvectors showed localized fluctuations. It was reported that the wild type contributed 41% variance by the first three eigenvectors to the total motion. For D49N, Y64S, W68G, and F94A, variance contribution by the first three eigenvectors is 55, 41, 63, and 32%, respectively. These results, particularly the D49N and W68G mutations, are significantly in uniformity with the RMSD, RMSF, and Rg results because these two mutations significantly affected the overall dynamics of the proteins and PZA binding.

We further plotted the principal components (PC1 and PC2) to cluster the trajectories motion for a perusable understanding. The conformational transition from one to another is represented in different colors (red to blue). Given in Figure 6, each dot represents a single frame from the trajectory. The mutant complexes variable phase space as compared to the WT. Together, all these results indicate that mutations significantly affect the structure that has led to the resistance against PZA drug.

FIGURE 6

Figure 6. Principal component analysis of all systems, including the Wild type and the four mutants. The first two principal components (PC1 and PC2) are used to project motion in the space phase at 300 K.

Destabilization of Fe²⁺ Ion by Mutations Induced Conformational Changes

Three histidine residues and one asparagine residue coordinate the Fe²⁺ ion. Li/Merz ion parameters for divalent Fe²⁺ ion was used to generate the topology. Mutations induced by Fe²⁺ destabilization during the simulation were determined by using the free energy landscape. It was found that Fe²⁺ is greatly influenced by the mutations. As given in RMSD and RMSF that the stability of each system is differentially affected, while the residual flexibility also showed variations. As presented in Figure 7, in the wild type structure, the Fe²⁺ did not move out significantly, but other regions showed little dynamic differences. The lowest energy conformation was attained at 92 ns. The only metastable state was extracted for wild type PZase is given below, which shows that the protein conformation is not altered during simulation.

FIGURE 7

Figure 7. Structural rearrangement of Fe²⁺ and the other regions in the protein given above (WT and D49N mutant). The lowest energy conformation from the wild type (92 ns) for the D49N (70 ns) was extracted and compared with the native state. The circle represents the lowest energy conformation.

On the other hand, as presented in Figure 7, the mutant system D49N showed destabilization of the Fe²⁺ ion. The structural coordinates extracted from the simulation trajectory at 70 ns represent the lowest conformation. In the case of D49N, the β-sheet two is significantly affected by transforming conformation.

Y64S has no significant effect on the enzyme and has the lowest energy conformation state attained at 72 ns. As given in RMSD and RMSF, the structural dynamics are not significantly affected by the Y64S mutation. All the analysis performed for Y64S in the manuscript discover consistent results and found Y64S as a comparatively less-lethal mutation than others. On the other hand, as reported above, W68G was significantly involved in structural destabilization and Fe²⁺ rearrangements. Along with the Fe²⁺ replacement and distortion of the coordination, the β-sheet 2 also flipped and thus causes a displacement of Asp49 residues that forms Fe²⁺ coordination along with the three histidine residues. The lowest conformational state of the W68G was extracted (5ns) after attaining the equilibrium (Figure 8).

FIGURE 8

Figure 8. Structural rearrangement of Fe²⁺ and the other regions in the protein given above (Y64S, W68G, and F94A mutants). The lowest energy conformation from each trajectory was extracted and compared with the native state. The circle represents the lowest energy conformation.

Mutation Diminishes the Binding Affinity of PZA

The MM-GBSA approach was employed to assess the binding affinity of WT and mutated receptors and ligand [1,2]. The last 10 ns trajectory, 500 snapshots, were used as input to estimate dominant forces between the protein and ligand interactions. The total binding free energies ΔGbind of WT and mutants (WT/-8.13, D49N/-5.93, Y64S/-4.88, W68G/-4.02, and F94A/-4.03) were calculated in kcal/mol (Table 2). The total energies of mutants compared to the WT indicates that these mutations drop the binding strength of the PZA. The vdW, Elec, and ΔPS energies contribution to the binding energies of the mutants compared to the WT were significantly low. It explores that the mutated proteins have weak binding to PZA. Mutations that are not involved in the direct interaction with the PZA affect orientation coordination of active site residues involved in direct contact with the PZA.

TABLE 2

Table 2. shows the binding affinity comparison between the wild type and mutant systems.

Discussion

Different studies have revealed that the administration of PZA, along with RIF and INH, is efficacious in treating Mtb infections (Gu et al., 2016; Khan et al., 2018c). Mtb resistance to these drugs renders front-line therapy ineffective, and as a consequence TB patient are exposed to a higher dose of the drugs. This leads to strong side effects on the patients and lower survival chances (Akhmetova et al., 2015; Khan et al., 2018d). Mitchison (1985), Miotto et al. (2014, 2017), Shi et al. (2014), Junaid et al. (2018), and Khan et al. (2020f). Since 1972, PZA was used as an active drug against the Mtb by targeting panD gene and had played key role in clearing persistent Mtb. Mutations in pncA, resulting in a loss of function of PZase, represent the primary molecular mechanism for PZA resistance in clinical strains. Pyrazinoic acid (POA) binds to the pncA active site and any conformational changes efflux the drug from the active site and this compromise the activity of the drug (Yadon et al., 2017). Additionally, the POA binding pocket is relatively small so any conformational changes result in unviability of the drug (Hewlett et al., 1995). It is clear that how the conformational changes affect the PZA binding, so, in this study, we selected D49N, Y64S, W68G, and F94A mutations at the Fe²⁺ binding site and PZA binding of PZase for their possible role in resistance to PZA (Miotto et al., 2017). We determined how the conformational changes may affect the binding of PZA and hinder the treatment of Mtb and eventually lengthen the eradication of tuberculosis.

In this regard, in silico techniques such as MD simulations were utilized to study the said mutations role in PZase resistance to PZA. These methods are widely used to understand the mechanism of resistance and any binding perturbation caused by mutations (Khan et al., 2020c,e, 2021a,b). It unveils conformational changes of proteins caused by any intrinsic mutations or ligand binding. This information is vital for devising novel strategies to combat drug resistance strains (Khan et al., 2019c, 2020b). Initial investigation of our selected mutations revealed that these mutations had altered the binding affinity of the PZA drug which shows that these mutations have clear role in resistance. Further characterization using biophysical tools revealed that RMSD and RMSF values suggest smaller fluctuations in wild type and higher fluctuations in mutant types. This probably suggest that wild type is highly stable, whereas more fluctuations in mutant type during the course of simulation suggest that the selected mutations are classified as highly destabilizing, and these findings are in line with previous experimental studies conducted on native and mutant (Q10P, D12A, G97D, R123P, T76P, G150A, H71R, W68R, W68G, and K96R). They reported that the mutations causes structural flexibility and thus weaken the drug binding (Khan et al., 2019b). The RMSF high fluctuation in mutant as compared to the WT discloses that the mutations have profound effect on the binding of the drug to the active site of the protein. A previous study carried out by Muhammad et al. also concluded that mutations in the PZA enzyme affect the binding orientation of PZA drug by shortening active pocket volume (Junaid et al., 2018, 2020; Khan et al., 2019c, 2020d). Findings of the current study may also suggest that said mutations affect the binding pocket, due to which the binding pocket volume as a whole is disturbed. Any distortion in the functional cavity volume might alter the binding affinity of PZA. This supports the previous study carried out by Vats et al. that mutations at the active pocket decrease the optimum affinity of the drug (Junaid et al., 2018, 2020). The residual flexibility also showed that each mutation displays a different frequency of fluctuations. Conformational dynamics, such as principal component analysis and free energy landscape, which are handy techniques reported by other studies, explored that the binding of Fe²⁺ is significantly affected. The main four residues coordinating the metal ion are disturbed during the simulation. In current study we observed different Rg pattern for the wild type and mutants. In case of wild type, initially Rg value increased and then it remains flat, whereas various patterns of increase/decrease were observed for all the mutants. These patterns suggest that the internal dynamics of each system is impacted by the mutation and eventually contributed to the PZA resistance. This notion is also supported by previous study (Jamal et al., 2020; Karmakar et al., 2020). The lowest energy minima conformation from each trajectory was extracted and compared with the native state that reported significant variations in Fe²⁺ binding and β–stands 2 specifically also confirmed by published literature (Khan et al., 2019c; Junaid et al., 2020). Analogous results have been reported that demonstrate that Fe²⁺ position is affected by catalytic and non-catalytic residues mutation. Furthermore, The Gibbs free energy to estimate the impact of the said substitutions on the binding of PZA. It was witnessed that mutations have significantly reduced the binding affinity of PZA and D49N and W68G being the major which contribute significantly to PZA resistance (Rehman et al., 2019). The study of dynamic behavior provided highly adequate knowledge on the PZase mutation that affected its structure, as well as perspectives into how conformational differences influence protein-ligand interactions which would aid the development of structure-based drug designing against the PZA target of Mtb.

In conclusion, we performed extensive MD simulations in triplets to explore the impact of D49N, Y64S, W68G, and F94A mutations on the PZase resistance to PZA. Our analysis revealed that these mutations affect stability, internal structural dynamics, and the binding energy of PZA. Our study further suggests that the stabilization of Fe²⁺ and β–stand 2 were affected. Hence, there is a dire need to design more potent drugs that would potently inhibit Mtb.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

AN, AK, SU, SSA, and SA conceptualized the study. AK, SS, SY, SHA, and MA did the analysis. AK, MA, and SSA draft the manuscript. AN, AK, SHA, and D-QW revised and finalized the manuscript. LA performed the revision and writing improvement. D-QW is an academic supervisor and supervised the study.

Funding

D-QW is supported by grants from the Key Research Area Grant 2016YFA0501703 of the Ministry of Science and Technology of China, the National Science Foundation of China (Grant Nos. 32070662, 61832019, and 32030063), the Science and Technology Commission of Shanghai Municipality (Grant No. 19430750600), the Natural Science Foundation of Henan Province (162300410060), as well as SJTU JiRLMDS Joint Research Fund and Joint Research Funds for Medical and Engineering and Scientific Research at Shanghai Jiao Tong University (YG2017ZD14). The computations were partially performed at the Pengcheng Lab. and the Center for High-Performance Computing, Shanghai Jiao Tong University.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2021.633365/full#supplementary-material

Footnotes

1. ^ http://biosig.unimelb.edu.au/mcsm/

References

Aggarwal, M., Singh, A., Grover, S., Pandey, B., Kumari, A., and Grover, A. (2018). Role of pncA gene mutations W68R and W68G in pyrazinamide resistance. J. Cell. Biochem. 119, 2567–2578. doi: 10.1002/jcb.26420

PubMed Abstract | CrossRef Full Text | Google Scholar

Akhmetova, A., Kozhamkulov, U., Bismilda, V., Chingissova, L., Abildaev, T., Dymova, M., et al. (2015). Mutations in the pncA and rpsA genes among 77 Mycobacterium tuberculosis isolates in Kazakhstan. Int. J. Tuberc. Lung Dis. 19, 179–184. doi: 10.5588/ijtld.14.0305

PubMed Abstract | CrossRef Full Text | Google Scholar

Amadei, A., Linssen, A. B., and Berendsen, H. J. (1993). Essential dynamics of proteins. Proteins 17, 412–425. doi: 10.1002/prot.340170408

PubMed Abstract | CrossRef Full Text | Google Scholar

Case, D. A., Cheatham, T. E. III, Darden, T., Gohlke, H., Luo, R., Merz, K. M. Jr., et al. (2005). The Amber biomolecular simulation programs. J. Comput. Chem. 26, 1668–1688. doi: 10.1002/jcc.20290

PubMed Abstract | CrossRef Full Text | Google Scholar

Chaturvedi, U. C., and Shrivastava, R. (2005). Interaction of viral proteins with metal ions: role in maintaining the structure and functions of viruses. FEMS Immunol. Med. Microbiol. 43, 105–114. doi: 10.1016/j.femsim.2004.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Dallakyan, S., and Olson, A. J. (2015). Small-molecule library screening by docking with PyRx. Methods Mol. Biol. 1263, 243–250. doi: 10.1007/978-1-4939-2269-7_19

CrossRef Full Text | Google Scholar

Darden, T., York, D., and Pedersen, L. (1993). Particle mesh Ewald: an N⋅ log (N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092. doi: 10.1063/1.464397

CrossRef Full Text | Google Scholar

DeLano, W. L. (2002). Pymol: an open-source molecular graphics tool. CCP4 Newsletter Protein Crystallogr. 40, 82–92.

Google Scholar

Dolatkhah, Z., Javanshir, S., Sadr, A. S., Hosseini, J., and Sardari, S. (2017). Synthesis, molecular docking, molecular dynamics studies, and biological evaluation of 4 h-chromone-1, 2, 3, 4-tetrahydropyrimidine-5-carboxylate derivatives as potential antileukemic agents. J. Chem. Inf. Model. 57, 1246–1257. doi: 10.1021/acs.jcim.6b00138

PubMed Abstract | CrossRef Full Text | Google Scholar

Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H., and Pedersen, L. G. (1995). A smooth particle mesh Ewald method. J. Chem. Phys. 103, 8577–8593. doi: 10.1063/1.470117

CrossRef Full Text | Google Scholar

Goddard, T. D., Huang, C. C., and Ferrin, T. E. (2005). Software extensions to UCSF chimera for interactive visualization of large molecular assemblies. Structure 13, 473–482. doi: 10.1016/j.str.2005.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, Y., Yu, X., Jiang, G., Wang, X., Ma, Y., Li, Y., et al. (2016). Pyrazinamide resistance among multidrug-resistant tuberculosis clinical isolates in a national referral center of China and its correlations with pncA, rpsA, and panD gene mutations. Diagn. Microbiol. Infect. Dis. 84, 207–211. doi: 10.1016/j.diagmicrobio.2015.10.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Hewlett, D., Horn, D. L., and Alfalla, C. (1995). Drug-resistant tuberculosis: inconsistent results of pyrazinamide susceptibility testing. JAMA 273, 916–917. doi: 10.1001/jama.1995.03520360030022

CrossRef Full Text | Google Scholar

Hoang, T. X., Trovato, A., Seno, F., Banavar, J. R., and Maritan, A. (2004). Geometry and symmetry presculpt the free-energy landscape of proteins. Proc. Natl. Acad. Sci. U. S. A. 101, 7960–7964. doi: 10.1073/pnas.0402525101

PubMed Abstract | CrossRef Full Text | Google Scholar

Jamal, S., Khubaib, M., Gangwar, R., Grover, S., Grover, A., and Hasnain, S. E. (2020). Artificial intelligence and machine learning based prediction of resistant and susceptible mutations in Mycobacterium tuberculosis. Sci. Rep. 10:5487. doi: 10.1038/s41598-020-62368-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Junaid, M., Khan, M. T., Malik, S. I., and Wei, D.-Q. (2018). Insights into the mechanisms of the pyrazinamide resistance of three pyrazinamidase mutants N11K, P69T, and D126N. J. Chem. Inf. Model. 59, 498–508. doi: 10.1021/acs.jcim.8b00525

PubMed Abstract | CrossRef Full Text | Google Scholar

Junaid, M., Li, C.-D., Li, J., Khan, A., Ali, S. S., Jamal, S. B., et al. (2020). Structural insights of catalytic mechanism in mutant pyrazinamidase of Mycobacterium tuberculosis. J. Biomol. Struct. Dyn. doi: 10.1080/07391102.2020.1761879 [Epub ahead of print]

PubMed Abstract | CrossRef Full Text | Google Scholar

Karmakar, M., Rodrigues, C. H. M., Horan, K., Denholm, J. T., and Ascher, D. B. (2020). Structure guided prediction of Pyrazinamide resistance mutations in pncA. Sci. Rep. 10:1875. doi: 10.1038/s41598-020-58635-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, A., Ali, S. S., Khan, M. T., Saleem, S., Ali, A., Suleman, M., et al. (2020a). Combined drug repurposing and virtual screening strategies with molecular dynamics simulation identified potent inhibitors for SARS-CoV-2 main protease (3CLpro). J. Biomol. Struct. Dyn. doi: 10.1080/07391102.2020.1779128 [Epub ahead of print]

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, A., Heng, W., Wang, Y., Qiu, J., Wei, X., Peng, S., et al. (2021a). In silico and in vitro evaluation of kaempferol as a potential inhibitor of the SARS-CoV-2 main protease (3CLpro). Phytother. Res. doi: 10.1002/ptr.6998 [Epub ahead of print]

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, A., Junaid, M., Kaushik, A. C., Ali, A., Ali, S. S., Mehmood, A., et al. (2018a). Computational identification, characterization and validation of potential antigenic peptide vaccines from hrHPVs E6 proteins using immunoinformatics and computational systems biology approaches. PLoS One 13:e0196484. doi: 10.1371/journal.pone.0196484

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, A., Junaid, M., Li, C.-D., Saleem, S., Humayun, F., Shamas, S., et al. (2020b). Dynamics insights into the gain of flexibility by Helix-12 in ESR1 as a mechanism of resistance to drugs in breast cancer cell lines. Front. Mol. Biosci. 6:159. doi: 10.3389/fmolb.2019.00159

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, A., Kaushik, A. C., Ali, S. S., Ahmad, N., and Wei, D.-Q. (2019a). Deep-learning-based target screening and similarity search for the predicted inhibitors of the pathways in Parkinson’s disease. RSC Adv. 9, 10326–10339. doi: 10.1039/c9ra01007f

CrossRef Full Text | Google Scholar

Khan, A., Khan, M., Saleem, S., Babar, Z., Ali, A., Khan, A. A., et al. (2020c). Phylogenetic analysis and structural perspectives of RNA-dependent RNA-polymerase inhibition from SARs-CoV-2 with natural products. Interdiscip. Sci. 12, 335–348. doi: 10.1007/s12539-020-00381-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, A., Khan, M. T., Saleem, S., Junaid, M., Ali, A., Ali, S. S., et al. (2020d). Structural insights into the mechanism of RNA recognition by the N-terminal RNA-binding domain of the SARS-CoV-2 nucleocapsid phosphoprotein. Comput. Struct. Biotechnol. J. 18, 2174–2184. doi: 10.1016/j.csbj.2020.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, A., Rehman, Z., Hashmi, H. F., Khan, A. A., Junaid, M., Sayaf, A. M., et al. (2020e). An integrated systems biology and network-based approaches to identify novel biomarkers in breast cancer cell lines using gene expression data. Interdiscip. Sci. 12, 155–168. doi: 10.1007/s12539-020-00360-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, A., Saleem, S., Idrees, M., Ali, S. S., Junaid, M., Kaushik, A. C., et al. (2018b). Allosteric ligands for the pharmacologically important Flavivirus target (NS5) from ZINC database based on pharmacophoric points, free energy calculations and dynamics correlation. J. Mol. Graph. Model. 82, 37–47. doi: 10.1016/j.jmgm.2018.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, A., Zia, T., Suleman, M., Khan, T., Ali, S. S., Abbasi, A. A., et al. (2021b). Higher infectivity of the SARS-CoV-2 new variants is associated with K417N/T, E484K, and N501Y mutants: an insight from structural data. J. Cell. Physiol. doi: 10.1002/jcp.30367 [Epub ahead of print].

CrossRef Full Text | PubMed Abstract | Google Scholar

Khan, M., Malik, S., Bhatti, A., Ali, S., Khan, A., Zeb, M., et al. (2018c). Pyrazinamide-resistant mycobacterium tuberculosis isolates from Khyber Pakhtunkhwa and rpsA mutations. J. Biol. Regul. Homeost. Agents 32, 705–709.

Google Scholar

Khan, M. T., Ali, S., Ali, A., Khan, A., Kaushak, A. C., Irfan, M., et al. (2020f). Insight into the PZA Resistance Whole Genome of Mycobacterium Tuberculosis Isolates from Khyber Pakhtunkhwa, Pakistan.

Google Scholar

Khan, M. T., Junaid, M., Mao, X., Wang, Y., Hussain, A., Malik, S. I., et al. (2019b). Pyrazinamide resistance and mutations L19R, R140H, and E144K in Pyrazinamidase of Mycobacterium tuberculosis. J. Cell. Biochem. 120, 7154–7166. doi: 10.1002/jcb.27989

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, M. T., Khan, A., Rehman, A. U., Wang, Y., Akhtar, K., Malik, S. I., et al. (2019c). Structural and free energy landscape of novel mutations in ribosomal protein S1 (rpsA) associated with pyrazinamide resistance. Sci. Rep. 9:7482.

Google Scholar

Khan, M. T., Malik, S. I., Ali, S., Sheed Khan, A., Nadeem, T., Zeb, M. T., et al. (2018d). Prevalence of pyrazinamide resistance in Khyber Pakhtunkhwa, Pakistan. Microb. Drug Resist. 24, 1417–1421. doi: 10.1089/mdr.2017.0234

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, S., Chen, J., Cheng, T., Gindulyte, A., He, J., He, S., et al. (2019). PubChem 2019 update: improved access to chemical data. Nucleic Acids Res. 47, D1102–D1109.

Google Scholar

Kräutler, V., Van Gunsteren, W. F., and Hünenberger, P. H. (2001). A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations. J. Comput. Chem. 22, 501–508. doi: 10.1002/1096-987x(20010415)22:5<501::aid-jcc1021>3.0.co;2-v

CrossRef Full Text | Google Scholar

Lemaitre, N., Callebaut, I., Frenois, F., Jarlier, V., and Sougakoff, W. (2001). Study of the structure–activity relationships for the pyrazinamidase (PncA) from Mycobacterium tuberculosis. Biochem. J. 353, 453–458. doi: 10.1042/bj3530453

CrossRef Full Text | Google Scholar

Malik, S. I., Ali, S., Masood, N., Nadeem, T., Khan, A. S., and Afzal, M. T. (2019). Pyrazinamide resistance and mutations in pncA among isolates of Mycobacterium tuberculosis from Khyber Pakhtunkhwa. Pakistan. BMC Infect. Dis. 19:116.

Google Scholar

Maningi, N. E., Daum, L. T., Rodriguez, J. D., Mphahlele, M., Peters, R. P., Fischer, G. W., et al. (2015). Improved detection by next-generation sequencing of pyrazinamide resistance in Mycobacterium tuberculosis isolates. J. Clin. Microbiol. 53, 3779–3783. doi: 10.1128/jcm.01179-15

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, B. R. III, McGee, T. D. Jr., Swails, J. M., Homeyer, N., Gohlke, H., and Roitberg, A. E. (2012). MMPBSA. py: an efficient program for end-state free energy calculations. J. Chem. Theory Comput. 8, 3314–3321. doi: 10.1021/ct300418h

PubMed Abstract | CrossRef Full Text | Google Scholar

Miotto, P., Cabibbe, A. M., Feuerriegel, S., Casali, N., Drobniewski, F., Rodionova, Y., et al. (2014). Mycobacterium tuberculosis pyrazinamide resistance determinants: a multicenter study. mBio 5:e01819–14.

Google Scholar

Miotto, P., Tessema, B., Tagliani, E., Chindelevitch, L., Starks, A. M., Emerson, C., et al. (2017). A standardised method for interpreting the association between mutations and phenotypic drug resistance in Mycobacterium tuberculosis. Eur. Respir. J. 50:1701354. doi: 10.1183/13993003.01354-2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Mishra, S. K., and Koča, J. (2018). Assessing the performance of MM/PBSA, MM/GBSA, and QM–MM/GBSA approaches on protein/carbohydrate complexes: effect of implicit solvent models, QM methods, and entropic contributions. J. Phys. Chem. B 122, 8113–8121. doi: 10.1021/acs.jpcb.8b03655

PubMed Abstract | CrossRef Full Text | Google Scholar

Mitchison, D. (1985). The action of antituberculosis drugs in short-course chemotherapy. Tubercle 66, 219–225. doi: 10.1016/0041-3879(85)90040-6

CrossRef Full Text | Google Scholar

Petrella, S., Gelus-Ziental, N., Maudry, A., Laurans, C., Boudjelloul, R., and Sougakoff, W. (2011). Crystal structure of the pyrazinamidase of Mycobacterium tuberculosis: insights into natural and acquired resistance to pyrazinamide. PLoS One 6:e15785. doi: 10.1371/journal.pone.0015785

PubMed Abstract | CrossRef Full Text | Google Scholar

Pires, D. E., Ascher, D. B., and Blundell, T. L. (2013). mCSM: predicting the effects of mutations in proteins using graph-based signatures. Bioinformatics 30, 335–342. doi: 10.1093/bioinformatics/btt691

PubMed Abstract | CrossRef Full Text | Google Scholar

Rehman, A. U., Khan, M. T., Liu, H., Wadood, A., Malik, S. I., and Chen, H.-F. (2019). Exploring the pyrazinamide drug resistance mechanism of clinical mutants T370P and W403G in ribosomal protein S1 of Mycobacterium tuberculosis. J. Chem. Inf. Model. 59, 1584–1597. doi: 10.1021/acs.jcim.8b00956

PubMed Abstract | CrossRef Full Text | Google Scholar

Roe, D. R., and Cheatham, T. E. III (2013). PTRAJ and CPPTRAJ: software for processing and analysis of molecular dynamics trajectory data. J. Chem. Theory Comput. 9, 3084–3095. doi: 10.1021/ct400341p

PubMed Abstract | CrossRef Full Text | Google Scholar

Rose, P. W., Prlić, A., Altunkaya, A., Bi, C., Bradley, A. R., Christie, C. H., et al. (2016). The RCSB protein data bank: integrative view of protein, gene and 3D structural information. Nucleic Acids Res. 45, D271–D281.

Google Scholar

Salomon-Ferrer, R., Götz, A. W., Poole, D., Grand, S. Le, and Walker, R. C. (2013). Routine microsecond molecular dynamics simulations with AMBER on GPUs. 2. Explicit solvent particle mesh Ewald. J. Chem. Theory Comput. 9, 3878–3888. doi: 10.1021/ct400314y

PubMed Abstract | CrossRef Full Text | Google Scholar

Schneidman-Duhovny, D., Inbar, Y., Nussinov, R., and Wolfson, H. J. (2005). PatchDock and SymmDock: servers for rigid and symmetric docking. Nucleic Acids Res. 33(Suppl._2), W363–W367.

Google Scholar

Shi, W., Chen, J., Feng, J., Cui, P., Zhang, S., Weng, X., et al. (2014). Aspartate decarboxylase (PanD) as a new target of pyrazinamide in Mycobacterium tuberculosis. Emerg. Microbes Infect. 3:e58.

Google Scholar

Stoffels, K., Mathys, V., Fauville-Dufaux, M., Wintjens, R., and Bifani, P. (2012). Systematic analysis of pyrazinamide-resistant spontaneous mutants and clinical isolates of Mycobacterium tuberculosis. Antimicrob. Agents Chemother. 56, 5186–5193. doi: 10.1128/aac.05385-11

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, Y., Hu, Z., Zhang, T., Cai, X., Kuang, H., Liu, Y., et al. (2014). Role of pncA and rpsA gene sequencing in detection of pyrazinamide resistance in Mycobacterium tuberculosis isolates from southern China. J. Clin. Microbiol. 52, 291–297. doi: 10.1128/jcm.01903-13

PubMed Abstract | CrossRef Full Text | Google Scholar

Toukmaji, A., Sagui, C., Board, J., and Darden, T. (2000). Efficient particle-mesh Ewald based approach to fixed and induced dipolar interactions. J. Chem. Phys. 113, 10913–10927. doi: 10.1063/1.1324708

CrossRef Full Text | Google Scholar

Wan, L., Liu, H., Li, M., Jiang, Y., Zhao, X., Liu, Z., et al. (2020). Genomic analysis identifies mutations concerning drug-resistance and beijing genotype in multidrug-resistant mycobacterium tuberculosis isolated from China. Front. Microbiol. 11:1444. doi: 10.3389/fmicb.2020.01444

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Wang, W., Kollman, P. A., and Case, D. A. (2001). Antechamber: an accessory software package for molecular mechanical calculations. J. Am. Chem. Soc. 222:U403.

Google Scholar

Yadon, A. N., Maharaj, K., Adamson, J. H., Lai, Y.-P., Sacchettini, J. C., Ioerger, T. R., et al. (2017). A comprehensive characterization of PncA polymorphisms that confer resistance to pyrazinamide. Nat. Commun. 8:588. doi: 10.1038/s41467-017-00721-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., and Mitchison, D. (2003). The curious characteristics of pyrazinamide: a review. Int. J. Tuberc. Lung Dis. 7, 6–21.

Google Scholar

Zwanzig, R. (1973). Nonlinear generalized Langevin equations. J. Stat. Phys. 9, 215–220. doi: 10.1007/bf01008729

CrossRef Full Text | Google Scholar