Mechanistic analysis elucidating the relationship between Lys96 mutation in Mycobacterium tuberculosis pyrazinamidase enzyme and pyrazinamide susceptibility
Continuous emergence of new mutations associated with drug resistance in Mycobacterium tuberculosis is posing a great danger to the conventional therapies. The methods need to be revised
according to the current situation, which requires understanding of the mechanisms
of resistance. In order to identify one such mechanism we have tried to analyze the
mutation K96R occurring in the active site of Mycobacterium tuberculosis. Various analyses have been done in order to understand the mechanism at molecular
level.
Preparation of mutant protein
Since the structure of mutant protein was not available in the database, in silico mutation was created at position 96 by replacing lysine with arginine. In order to
allow the protein to rearrange the backbone and the side chains of the mutated residue
and the nearby amino acids so as to reach the conformational state with no steric
clashes and minimum internal energy, it was subjected to MD simulation for 15 ns.
RMSD of each frame in reference to the first frame was plotted against the simulation
time to comment on the equilibrium state of the mutated protein. Initially an RMSD
of 4Ã… was observed for the first 6 ns after which an almost stable trajectory was
obtained. The coordinates of all the frames within the most stable time frame, i.e.
from 6 to 15 ns, were averaged to compute a structure representing the equilibrated
state of the mutated protein. This representative average structure was used for further
study.
Binding cavity analysis
Binding pocket analysis is an important measure to identify the impact of mutation
on the cavity. The pocket volume for both native and wild type were calculated using
CASTp server. The structures obtained after molecular dynamics were used for the calculation
of the pocket volume. The volume of the cavity for native PncA protein was 551.9 Ã…3 and that of mutated form was 1314.2 Ã…3. A considerable difference can be seen in the two volumes (Figure 1). The increase in volume affects the binding of PZA within active cavity thereby
decreasing the stability of the ligand. In order to further evaluate the behavior
of the mutant, we performed docking and molecular dynamics simulation studies.
Figure 1. Ribbon representation of the superimposition of structures of the native and mutant
M. tuberculosis PncA protein.
The change in the cavity size was further validated using MOLE 2.0. A significant
increase in the active site volume was observed for the mutant protein with this program
as well.
Interaction analysis between receptor and ligand molecule
Previous studies on PncA have shown the significance of van der Waals interactions
in binding to PZA. Binding energy calculations provide information on both van der
Waals and electrostatic energy of the complex. The value of van der Waals energy and
electrostatic energy for native complex calculated using Glide was found to be -17.16
kcal/mol and -10.27 kcal/mol respectively. The glide energy or the total ligand-receptor
energy was found to be -27.43 kcal/mol for the native PncA and PZA complex. On the
contrary, there was a significant difference in the values for the complex with mutated
protein. The values for van der Waals, electrostatic and glide energy were -10.86,
-13 and -23.43 kcal/mol respectively (Table 1). As indicated by the total energy, PZA in complex with mutated PncA was less stable
in comparison to PZA docked with native PncA. Hydrogen and hydrophobic interactions
also play a major role in making the binding stable. The interaction plots for these
two types of interactions were generated using Ligplot 48].
Table 1. Scores from docking analysis of native and mutant protein
In native complex, an oxygen atom of PZA was found making one hydrogen bond with the
nitrogen of Ala 134 of PncA and the distance was 3.19Ã… (Figure 2A). The ligand also made a significant number of hydrophobic bonds with Asp 8, Phe
13, leu1 9, Val 21, Ile 133, His 137, Cys 138 and Val 163 of PncA (Figure 3A). On the other hand the mutant protein was found to be making two hydrogen bonds.
An oxygen atom of the ligand was making hydrogen bonds with nitrogen of Ala 134 and
Cys 138 of mutated protein and the bond lengths were 3.11 Ã… and 3.08 Ã… respectively
(Figure 2C). However, the number of residues involved in making hydrophobic contacts was less
in the case of mutated protein-PZA complex. Amino acids of mutated PncA participating
in hydrophobic interactions included Asp 8, Gln 10, Ile 133 and His 137 (Figure 3C). We observed a difference in the interaction pattern of the two proteins with PZA.
Cys 138, one of the key active residues was involved in hydrophobic interaction in
native protein complex while in mutant complex it was forming a hydrogen bond with
the ligand. Over all, PZA had a lower binding affinity towards the mutated form of
PncA as indicated by the docking score of -3.77 in comparison to the wild type protein,
with which it had a glide docking score of -4.025. To further validate the binding
affinity two online servers, PatchDock and SwissDock were also used. The results obtained
indicated a lower affinity of PZA for the mutated protein in comparison to the wild
type. It was in coherence with the results obtained using Glide.
Figure 2. Depiction of hydrogen bond interactions in A) native protein bound to PZA after docking,
B) native protein bound with PZA after MD simulations, C) docked mutant PncA-PZA complex,
D) mutant protein bound to PZA post MD simulations.
Figure 3. Hydrophobic interaction pattern in A) native protein bound to PZA after docking, B)
native protein bound with PZA after MD simulations, C) docked mutant PncA-PZA complex,
D) mutant protein bound to PZA post MD simulations.
All the three factors, namely, volume of binding pocket, docking score and interaction
energy, typically used here to explain the reason for resistance to drug supported
the idea of weak binding of the ligand with the protein in mutated form, which was
hampering its activation. However, increased number of hydrogen bonds in mutant complex
was suspicious. Therefore, in order to study the conformational changes occurring
in the complexes in dynamic environment similar to in vivo conditions, we performed molecular dynamics simulation studies. Calculation and analysis
of parameters like RMSD, radius of gyration and RMSF was done to understand the behavior
of binding of the ligand with the enzyme.
Investigation of flexibility of native and mutant PncA enzyme
PZA bound native and mutant proteins were simulated in a SPC water box for around
15 ns. RMSD trajectory was analyzed for both the complexed systems. As shown in Figure
4, the native unbound and native protein complex, both were quite stable. The native
unbound display low fluctuations and remain stable throughout the simulations. The
native bound also did not deviate much from its initial docked conformation. However,
in the mutant structure, the protein complex showed deviation up to 3Ã… in the first
half phase of the simulation run and then gave a stable RMSD trajectory (Figure 4).
Figure 4. Graph showing the RMSD of backbone of native unbound protein, native bound to PZA,
mutant protein in unbound form and mutant protein in complex with PZA.
Some changes were observed in the interaction pattern of the native and mutant complexes
post molecular dynamics studies. In native complex, the number of hydrogen bonds increased
to two. The oxygen atom of PZA was making two hydrogen bonds with nitrogen of Ala
134 and Cys 138 with bond lengths of 3.33Ã… and 2.89Ã… respectively (Figure 2B). The hydrophobic interactions however reduced. PZA was now making hydrophobic contacts
with Asp 8, Phe 13, Val 21, Ile 133 and His 137 of native PncA (Figure 3B). On the other hand, mutant complex was now making only one hydrogen bond with Ala
134 (Figure 2D). The bond length for the corresponding hydrogen bond was observed to be 3.06Ã…. Seven
amino acid residues of mutant proteins, namely Asp 8, Gln 10, Phe 13, Asp 49, Ile
133, His 137 and Cys 138 were now interacting hydrophobically with PZA (Figure 3D).
The radius of gyration is defined as the mass-weighted root mean-square distance of
a cluster of atoms from their common center of mass. In other words, the radius of
gyration of a protein is a measure of its compactness. A protein will maintain a steady
value of Rg if it is stably folded, and the value will fluctuate if the protein unfolds. The radius
of gyration of all the frames during the simulation run was plotted against time and
the data was analyzed. The curve for the native complex was stable throughout the
simulation period and fluctuations were negligible. This steady curve can be attributed
to high stability and compactness of the protein. The plot for the mutant protein
complex was however more fluctuating and the value for Rg reached as high as 17 Ã… (Figure 5). This fluctuation was prevalent during the entire simulation period. Hence it can
be inferred that the protein was losing its compactness due to the change in its conformation.
Figure 5. Radius of gyration of native and mutant protein in complex with PZA.
Root mean square fluctuation (RMSF) with respect to each residue in both native and
mutant complex was then calculated to explain the loosening of the mutant protein
structure. The fluctuations were widespread in the mutant protein, while in native
sequence small fluctuations were observed in certain regions. The fluctuations in
certain residues of native protein were conferring flexibility to the protein and
keeping it compact at the same time. While in case of mutant protein, high level of
fluctuations depicting high mobility of the residues was resulting in a loosely packed
protein. As illustrated in Figure 6, fluctuations were highest in the loop region i.e. from 51 to 71 amino acid residues.
This result was in accordance with information obtained from superposition shown in
Figure 1. The loop region of the mutated protein moved away from the active cleft, which was
also the reason behind increase in the volume of catalytic cavity. According to a
study, residues 51 to 71 compose an important loop region or a flap which is involved
in holding the ligand in position after binding. Thus this loosening of the loop region
was affecting the binding of PZA with the active residues.
Figure 6. Root mean square fluctuation of all the residues in native protein bound to PZA, mutant
protein in unbound form and mutant protein in complex with PZA.
The RMSF of catalytically important binding residues, Asp 8, Phe13, IIe133, Ala134
and Cys138 of native and mutant structures were also analyzed. Table 2 lists the RMSF values of all these residues in unbound mutant, mutant and native
protein bound with PZA. The fluctuation in binding residues was very small in case
of native protein while these residues were highly flexible in case of mutant protein
in both bound and unbound form.
Table 2. RMSF values of binding cavity residues
All these analyses brought the discussion to the point that in mutated form the flap
region of the protein becomes more flexible and moves away from the active site. This
results in a much bigger active cavity. Once the ligand, PZA comes and binds to the
active residues, flap is not able to hold it in bound state thereby preventing the
activation of the prodrug.
Previously studies have been performed on PncA mutations namely, D8G, S104R and C138Y.
The analysis showed that the binding cavity volume decreased and the cavity became
rigid. This did not permit proper alignment of PZA inside the cavity and hampered
hydrogen bond formation. Our study included a novel mutation. The mutation being present
in the flap region, affected the volume of the cavity in contrast to the previous
mutations. The mutation led to a substantial increase in the binding cavity. This
prohibited the enzyme from holding the drug properly and therefore PZA could not take
its active form 24].