Gibbs Free Energy Calculation of Mutation in PncA and RpsA Associated With Pyrazinamide Resistance.

A central approach for better understanding the forces involved in maintaining protein structures is to investigate the protein folding and thermodynamic properties. The effect of the folding process is often disturbed in mutated states. To explore the dynamic properties behind mutations, molecular dynamic (MD) simulations have been widely performed, especially in unveiling the mechanism of drug failure behind mutation. When comparing wild type (WT) and mutants (MTs), the structural changes along with solvation free energy (SFE), and Gibbs free energy (GFE) are calculated after the MD simulation, to measure the effect of mutations on protein structure. Pyrazinamide (PZA) is one of the first-line drugs, effective against latent Mycobacterium tuberculosis isolates, affecting the global TB control program 2030. Resistance to this drug emerges due to mutations in pncA and rpsA genes, encoding pyrazinamidase (PZase) and ribosomal protein S1 (RpsA) respectively. The question of how the GFE may be a measure of PZase and RpsA stabilities, has been addressed in the current review. The GFE and SFE of MTs have been compared with WT, which were already found to be PZA-resistant. WT structures attained a more stable state in comparison with MTs. The physiological effect of a mutation in PZase and RpsA may be due to the difference in energies. This difference between WT and MTs, depicted through GFE plots, might be useful in predicting the stability and PZA-resistance behind mutation. This study provides useful information for better management of drug resistance, to control the global TB problem.


INTRODUCTION
Evolution may have optimized proteins to perform proper functions, native to the host organism, in different environmental conditions. Pharmaceutical industries desire changes in the thermodynamic properties of a protein (Liszka et al., 2012;Gapsys et al., 2016) to enhance the thermal stability, improving the protein-protein interactions. These desired changes are oftenly accomplished by mutations, and the free-energy changes are predicted to gain the desired properties. However, natural mutation in drug target may cause a resistance to the therapeutic drugs. Such mutations pose a great threat to the treatment of major infectious diseases. Understanding the forces like thermodynamic properties and protein folding involved in maintaining the protein structures, is of central interest when working on drug resistance. The folding process is most often affected by mutations (Carra and Privalov, 1996). To explore the dynamic properties behind mutations, molecular dynamic (MD) simulations have been widely performed, and have been especially useful in unveiling the mechanism of drug failure behind mutation (Carter Childers and Daggett, 2017;Dong et al., 2018;Hashemzadeh et al., 2019;Kaushik et al., 2019). MD simulation studies of ligand-protein interactions are a widely applied approach for explaining the mechanisms of drug resistance behind mutations (Aggarwal et al., 2017;Carter Childers and Daggett, 2017;Bera et al., 2018;Liu et al., 2018;Pandey et al., 2018;Ishima et al., 2019). During in vivo analysis, the crystal structure is analyzed for drug resistance. However, it can be formed based on some experimental conditions where none of the protein-drug complexes provide the mechanism of resistance, and none of the structures can be attained by X-ray. Investigating the insight mechanism at molecular level, MD simulation has got a certain advantage over experimental approaches of exploring drug resistance behind mutations (Liu and Yao, 2010;Khalaf and Mansoori, 2018;Liu et al., 2018;Meng et al., 2018;Mehmood et al., 2019). Furthermore, the dynamics and residues level analysis could be performed which was difficult to achieve through experimental approaches (Hou et al., 2008;Xue et al., 2012;Ding et al., 2013;Khan et al., 2018).
The effect of mutations on a protein complex is experimentally performed by different methods including isothermal titration calorimetry (ITC) (Ghai et al., 2012), surface plasmon resonance (Masi et al., 2010), Fluorescence resonance energy transfer (FRET) (Phillip et al., 2012), and some other procedures as described earlier (Kastritis and Bonvin, 2013). However, all these techniques are considered to be time consuming as well as costly. The mechanism of resistance behind mutation is of key interest where free energy is commonly altered. To estimate changes in the thermodynamics of wild types and mutant proteins, MDbased free energy calculations allow a precise measurement of changes (Aldeghi et al., 2019). Gibbs free energy (GFE) or free enthalpy (Greiner et al., 1995;Matthews, 2000;Li et al., 2014;Rietman et al., 2016) can be used to estimate the maximum level at which the process is reversible, performed through a thermodynamic system. The GFE is the non-expansion work, calculated from a thermodynamically closed system where this maximum can be achieved individually in an entirely reversible procedure. The reversible transformation of a system is going to decrease in GFE, from initial state to a final state, equal to the work done by the system to its surroundings, minus the work of the pressure forces (Matthews, 2000).
The internal motion of the system is measured using Principal Component Analysis (PCA), which is performed on the mass-weighted cartesian coordinates, and the long dynamics are able to recognize low modes in proteins (Jencks, 1981;Rajendran et al., 2018). In a long trajectory, PCA reduces the complicated motion (Novotny, 1991;Zídek et al., 1999;Datar et al., 2006). In a comparative analysis of two sets of proteins, a transformed set of variables z1, z2. . ., zp called principal components (PCs) where the PC1 and PC2 are the first two components, give the trajectories on the primary two principal components of motion (Verma et al., 2008;Martis et al., 2015).
Binding free energy calculations yield either absolute free energies (Molecular mechanics generalized Born surface area and Molecular mechanics Poisson-Boltzmann surface area) or relative free energies (Alchemical method) (Michel and Essex, 2010;Chodera et al., 2011;Mobley and Klimovich, 2012). Alchemical free energy calculations work by introducing a series of intermediate unphysical states spanning between the desired end states. Molecular docking combined with MD simulations followed by Molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) analysis is an efficient approach for Free energy calculation. The results of MM/PBSA are in reasonable agreement with previous experiments Kollman, 2000, 2001; and less computationally demanding than alchemical free energy methods. These two methods have been widely applied in biomolecules such as protein folding, protein-ligand binding, protein-protein interaction, etc. (Hou, 2010;Xu et al., 2013;Chen et al., 2015Chen et al., , 2019Sun et al., 2018;Wang et al., 2019). MM-PBSA and Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) have been the two most efficient methods to rapidly evaluate binding ability and to compute binding free energies (Hou, 2010;Sun et al., 2014a).
In previous studies, we have investigated the PZA drug sensitivity testing and then sequencing to find mutations in pncA and rpsA genes associated with PZA-resistance ,a, Khan et al., 2019c. MD simulation of some MTs in comparison with WT have been investigated as the cause behind resistance (Junaid et al., 2018;Khan et al., 2018Khan et al., ,a,c, 2019bRehman et al., 2019). In the current paper, we aimed to reanalyze the free energy differences, predicted via MM/GBSA and MM/PBSA, of WT and MTs that may be applied as a measure of stability in the binding affinity of drug and targets.

Mutants Selection in pncA
The primary cause behind PZA resistance have been associated with mutations in the pncA gene. The majority of studies have been conducted to investigate the drug resistance mechanism behind mutation by analyzing the root mean square deviation, root mean square fluctuation, and motion of MTs and WT PZase. However, the comparison of free energy as a mechanism of changes that occur behind a mutation is required to be investigated for better understanding of PZA-resistance. Here we selected N11K, P69T, D126N, L19R, R140H, and E144K to analyze the effect of mutations on free energy by comparing the MTs and WT (Junaid et al., 2018;Khan et al., 2018,a). A threedimensional structure (PDB ID 3pl1) was retrieved from the Brookhaven Raster Display (BRAD) protein data bank (PDB) (Berman et al., 2000). Using the mutate_Model script of Modeller (Webb and Sali, 2016) and PYMOL (DeLano, 2002), mutants were created at specific locations.

Protein-Ligand Interaction
Protein and ligand structures were prepared as described in earlier studies (Aggarwal et al., 2017;Friesner et al., 2004) using MOE. Incorrect hydrogen atoms were corrected and selenomethionine were changed into methionine. Protein-drug interactions were examined in MOE as a flexible docking. WT and MTs structure were subjected to MD simulations in apo and complex with the drug.

Molecular Dynamics Simulation (MD)
MD simulation was performed on all the MTs and WT using the Amber14 package (Salomon-Ferrer et al., 2013;Sun et al., 2014a,b) with the ff14SB force field. The TIP3P water model was used to solvate each system and counterion were added to neutralize the system (Jorgensen et al., 1983). The neutralized systems were minimized with the steepest descent minimization step (6000 cycles) and conjugate gradient (3000 cycles) followed by heating upto 300K. The systems were equilibrated at 1 atm and 300 K. For control of the temperature, the Langevin thermostat was turned on. For Long-range electrostatic interactions, the Particle Mesh Ewald algorithm was used (Darden et al., 1993;Essmann et al., 1995) and the treatment of the covalent bonds was performed with the SHAKE algorithm (Ryckaert et al., 1977). The production step of MD simulation was performed with pmemd code 30 (Götz et al., 2012). The cpptraj package in Amber 14 was used to analyze the trajectories.

Principal Component Analysis and Gibbs Free Energy Calculation
The high fluctuations in residues of protein were captured through principal component analysis (PCA) (Amadei et al., 1993) while variation in GFE values has been accounted for in the calculation of stability level in proteins molecules to perform proper function (Rajendran et al., 2018;Martis and Coutinho, 2019;Sohaib Shahzan et al., 2019). GFE calculation is a useful process for understanding the thermodynamic properties of antibody-antigen complex formation and proteinsproteins interactions (Jencks, 1981;Novotny, 1991;Zídek et al., 1999). Using a cpptraj package, the covariance matrix was calculated considering only the Cα coordinates followed by the diagonalization to calculate the eigenvectors and eigenvalues. PCA was calculated from the trajectory, containing 5000 snapshots. PC1 and PC2, the first two components, were used for the plotting. The binding free energy was calculated as described in previous studies (Sun et al., 2014a,b;Martis and Coutinho, 2019).
The binding events involved many interactions (Wang and Wade, 2001;Datar et al., 2006;Verma et al., 2008;Martis et al., 2015), therefore the classical binding free energy equation may be written as follows: In equation 2, Gs ol : solvation energy, G conf : conformational energy, G int energy due to interaction with residues in the vicinity, G motion : energy of motions (translational, rotational, and vibrational).
The free energy landscape (FEL) was developed using g_sham module to capture the lowest energy stable state. The deep valleys on a plot show the stable state while the boundaries between deep valleys represent the intermediate conformations (Hoang et al., 2004). The first two principal components were used to calculate the FEL based on the equation: PC1 and PC2 are reaction coordinates, KB symbolizes the Boltzmann constant, and P (PC1, PC2) illustrate the probability distribution of the system along the first two principal components.
The changes in enthalpy ( H), standard free energy ( G), and entropy ( S) are calculated using the following equation (Basu, 2010;Gautam and Chattopadhyaya, 2016); Where, H = Enthalpy, T = temperature in Kelvin, S = entropy, G = Gibbs Free Energy. In the current review we analyzed the GFE of MTs and WT PZase and RpsA that might be useful to measure the resistance among drug target proteins for better management of drug resistance.

Solvation Free Energies of Wild Type and Mutants
The solvation free energy is the product of the atomic solvation parameter and the accessibility of the atom to the solvent. This method estimates the relative stability of protein conformations, and estimates the free energy of proteins binding to ligands (Eisenberg and McLachlan, 1986). The stability and fluctuation of protein are measured through the solvation. Protein and solvent interactions at atomic level is quantified by solvation free energy (SFE). Free energy of protein hydration (solvation) is carried out with explicit solvent and all-atom treatment (Weber and Asthagiri, 2012;Kokubo et al., 2013;Matubayasi, 2017). Here we calculated the Solvation Free Energy ( Gsolv) of WT and MTs PZase and RpsA to find the effect of mutation on the proteins solvation free energy.

PCA and Entropy
Intra-protein information is transmitted over distances via allosteric processes. This ubiquitous protein process allows for protein function changes due to ligand binding events. Understanding protein allostery is essential in protein functions. Allostery in the protein has been inspected using a rigid residue scan method along with configurational entropy calculation and PCA. Based on a covariance correlation analysis of simulations, the contributions from individual residues to whole-protein dynamics have been systematically assessed and the entropic contributions of individual residues to whole-protein dynamics were also evaluated. When individual residues are held rigid, the FIGURE 1 | Gibbs free energy of PZase and mutants (N11K, P69T and D126N) in apo and complex state with PZA. The Gibbs free energy landscape for wild and mutated proteins in their apo and wild states is depicted along with their value bars against PC1 and PC2. Noticeable differences can be observed. The red color represents the high energy state, yellow and green low and blue represents the lowest stable state. variations of overall protein entropy favor the rigidity/flexibility equilibrium in protein structure. Further, the change of entropic contribution from each residue has been linked to the intrinsic differences among all the residues. These findings provide a systematic approach to dig out the contribution of individual residue's internal motion to overall protein dynamics and allostery (Bhakat et al., 2014;Kalescky et al., 2016).

Gibbs Free Energy Comparison Between Wild Type and Mutant in PncA
Geographically distinct and novel mutations have been detected in our recent studies ,b, Khan et al., 2019c after the drug susceptibility testing followed by pncA and rpsA sequencing of PZA resistance Mycobacterium tuberculosis isolates. Changes in values of GFE might be important in calculating the stability of proteins' confirmation. In order to explore the protein conformational shift from WT to mutant, the GFE for the first two principal components (PC1 and PC2) has been calculated. The energy landscape of both the apo and complex states of WT, and three mutants, N11K, P69T and D126N have been shown in Figure 1. The minimum energy area is indicated by the blue color. WT protein shows a clear large global energy minima basin (in blue), whereas the MTs reveal several different energy minima states. The blue areas depict more stability while more blue areas indicate transitions in the protein conformation followed by the thermodynamically more favorable state. The WT shows low energy state as compared to the MTs. The result demonstrates that native PZase has a more stable cluster as compared to the MTs that might be involved in low binding affinity with PZA, causing resistance (Junaid et al., 2018;Yang et al., 2018). Calculating the GFE in case of PZA resistance might be a useful way to analyze the MTs stability and also aid in alternative drug discovery.
The differences in GFE values of WT and MTs PZase, L19R, R140H, and E144K showed that mutations may alter the stability (Figure 2) which could be a measure to evaluate the PZA resistance.  ,b, Khan et al., 2019bRehman et al., 2019). The MtRpsA CTD is the POA binding site. All these MTs at 100 and 50 ns of MD simulations showed significant change in the structure and activity of RpsA (Figures 3-5).

Gibbs Free Energy Comparison Between Wild Type and Mutants RpsA
A native RpsA structure has the minimum GFE ,  exhibiting  significant  variations  when  compared  with MTs, D342N, D343N, A344P and I351F FIGURE 3 | Gibbs free energy (GFE) of WT and MTS, D342N, D343N, A344P, and I351F in apo and complex states. Wild type has a significant GFE difference to MTs as indicated by the color of the GFE plot. WT exhibited a more stable state as compared to mutants. POA resistance might be due the GFE states altering the affinity of RpsA (Khan et al., 2019b).
FIGURE 4 | Comparison of Gibbs free energy of MTs and wild type RpsA. WT exhibited a significant difference in GFE as indicated by the peak color of GFE plot.
( Figure 3). The color (red) in the plot is more prevalent in mutants, and seems less stable when compared with WT. The differences in GFE values of other MTs RpsA (S324F, E325K, and G341R) has been shown (Figure 4), revealing that they may have altered the stability of MTs RpsA. WT attained a significant value in comparison with MTs ( Figure 5). The peak color in both states of the native is seems to be more stable, indicating the importance of GFE calculation when measuring the effect of mutation on proteins dynamics characteristics.
In a more insightful study of WT and MTs, T370P and W403G RpsA, the comparison shows a significant difference not only in GFE states but also a difference in the loop structure ( Figure 5). The loop structure in MTs is seemed to be more open in both apo and complex with POA. These changes may cause POA resistance, resulting in weak or no binding with RpsA. Further, the high energy state might be involved to exhibit a more open loop residue. However, further confirmation through experimental approaches will enhance the understanding of low, medium, and high levels of POA resistance.
The differences in GFE values may have effects on the binding affinity and the stability calculation, resulting in weak interactions or loss of interactions with POA. In a farther site mutation, WT exhibited a significant difference in GFE in comparison with mutants, T370P and W403G (Figure 5). Mutations in C-terminal site of RpsA might be involved in the alteration of GFE, resulting in a loss of binding affinity with the drug.

Solvation Free Energies of Wild Type and Mutant PZase and RpsA
Hydrogen bonding is an important part of molecular interactions where the solvent is water. Free energy of protein hydration (solvation) is carried out with explicit solvent and all-atom treatment (Weber and Asthagiri, 2012;Kokubo et al., 2013;Matubayasi, 2017). The solvation free energy is the product of the atomic solvation parameter and the accessibility of the atom to the solvent. This method estimates the relative stability of protein conformations, and free energy of proteins binding to ligands (Eisenberg and McLachlan, 1986). Protein and solvent interactions at an atomic level is quantified by solvation free energy (SFE). The Solvation Free Energy ( Gsolv) of WT and MTs have been given ( Table 2). Interestingly the MTs PZase exhibited a decreased level of SFE than WT except P69T (−681.60) and D126N (−385.29). Similarly during the investigation of protein kinetics and thermodynamics, the MTs, Y91Q exhibited a decreased SFE and hydrophobicity compared to WT. This may be due to the more exposed and solvated hydrophilic side chains in the R1-region in acylphosphatase . Another important property while studying the protein thermodynamics is the energy of solvation (ES), recorded when dissolving a solute in a solvent. A positive and negative SE represents endothermic and exothermic processes respectively. This process of solvation is thermodynamically favored only when the overall GFE of the solution is decreased, as compared to the GFE of separated solvent and solute. A negative value is obtained when the change in enthalpy minus the change in entropy is multiplied by the absolute temperature or GFE of the system decreases. All the MTs exhibited lower ES than WT except P69T and D126N (Table 2). Similarly MTs RpsA attained a much lower ES than WT except E325K. Entropy of WT and MTs has been found in significant variation, a measure of a system's thermal energy per unit temperature, unavailable for useful work. Molecular disorder, or randomness, of a system may also be measured through entropy (Chong and Ham, 2012;Caro et al., 2017;Verteramo et al., 2019).
Overall, SFE changes by point mutation in PZase and RpsA casusing PZA-resistance during TB treatment regime. The SFE is commonly influenced by the hydrophilic residues. In a previous study the SFE of Y91Q has been found lower by 25.1 kcal/mol than WT acylphosphatase, indicating that Y91Q is less hydrophobic . Two mutations (R1s40H, E144K) that have been detected in α-helix of PZase exhibited the lowest SFE and SE as shown (Table 2 and Figures 6, 7). All the MTs RpsA attained lower SFE and SE than WT except E325K and G341R. Further, the solvation entropy of all the MTs is higher than WT (−3664.43 kcal/mol) except E325K (−3682.23 kcal/mol) and G341R (−3705.82 kcal/mol). The standard deviation of total free energy has also been given in Table 1

CONCLUSION
In the current analysis GFE along with SFE and SE of WT and MTs exhibited a significant difference which might be useful in predicting the drug resistance level behind mutations in PZase and RpsA. Molecular dynamics simulations, binding free energy, and PCA clearly show the impact of mutations on thr thermodynamics of proteins. These findings depict that mutations affect the overall enzyme's conformational landscape and distort the atomic interaction network. The GFE differences provide rapid potential, key for further designing of novel inhibitors to combat MTB resistant strains. The physiological effect of mutations in drug targets might be due to the energy differences. Evolutionary pressures might have maintained a protein folding integrity and stability while mutations may have decreased and posed severe consequences in disturbing bonds of intrinsic energy. The level of resistance might be analyzed through further experimental analysis and alternative drug discovery for better achieving the goals of the global TB eradication program 2030.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.