- 1School of Pharmacy, Lanzhou University, Lanzhou, China
- 2State Key Laboratory of Applied Organic Chemistry and Department of Chemistry, Lanzhou University, Lanzhou, China
Tuberculosis is still one of the top 10 causes of deaths worldwide, especially with the emergence of multidrug-resistant tuberculosis. Rifampicin, as the most effective first-line antituberculosis drug, also develops resistance due to the mutation on Mycobacterium tuberculosis (Mtb) RNA polymerase. Among these mutations, three mutations at position 451 (H451D, H451Y, H451R) are associated with high-level resistance to rifampicin. However, the resistance mechanism of Mtb to rifampicin is still unclear. In this work, to explore the resistance mechanism of Mtb to rifampicin due to H451D/Y/R mutations, we combined the molecular dynamics simulation, molecular mechanics generalized-Born surface area calculation, dynamic network analysis, and residue interactions network analysis to compare the interaction change of rifampicin with wild-type RNA polymerase and three mutants. The results of molecular mechanics generalized-Born surface area calculations indicate that the binding free energy of rifampicin with three mutants decreases. In addition, the dynamic network analysis and residue interaction network analysis show that when H451 was mutated, the interactions of residue 451 with its adjacent residues such as Q438, F439, M440, D441, and S447 disappeared or weakened, increasing the flexibility of binding pocket. At the same time, the disappearance of hydrogen bonds between R613 and rifampicin caused by the flipping of R613 is another important reason for the reduction of binding ability of rifampicin in H451D/Y mutants. In H451R mutant, the mutation causes the binding pocket change too much so that the position of rifampicin has a large movement in the binding pocket. In this study, the resistance mechanism of rifampicin at the atomic level is proposed. The proposed drug-resistance mechanism will provide the valuable guidance for the design of antituberculosis drugs.
Introduction
Tuberculosis (TB), an infectious disease caused by Mycobacterium tuberculosis (Mtb), is the leading cause from a single infectious agent worldwide. Mtb, a pathogenic bacterium species of the family Mycobacteriaceae, can attack the lung of people and spread in the population through the droplets from the throat of an active TB infection patient (Mishra and Surolia, 2018). Millions of people continue to fall sick with TB each year. The 2018 WHO Global Tuberculosis Report (Organization, 2018) estimated 10.0 million new cases of TB and 1.6 million deaths in 2017.
Although it spreads widely, TB is preventable and curable. The 6-month short-course regimen with a combination of four anti-TB drugs (rifampicin, isoniazid, pyrazinamide, and ethambutol) has been used as a standard treatment for active, drug-susceptible TB patients over the past decades (Service and Council, 1981; Chang et al., 2018; Seid et al., 2018; Tiberi et al., 2018). Furthermore, the risk of relapse is generally below 5% reported among drug-susceptible TB patients after treated with standard 6-month regimens in clinical trials (Chang et al., 2006). However, the multidrug-resistant tuberculosis (MDR-TB), which at least resists to both rifampicin and isoniazid, had been emerged in the early 1990s due to multiple factors (He et al., 2008; Sandgren et al., 2009; Ahmad et al., 2018).
Rifampicin, one of the most effective anti-TB drugs, has been used as the first-line treatment in drug-susceptible TB patients, and it is also effective against initial isoniazid resistance (Mitchison and Nunn, 1986). Unfortunately, rifampicin resistance in Mtb arises due to the residues' mutations on its molecular target, Mycobacterium tuberculosis RNA polymerase (Mtb-RNAP). More than 95% of the rifampin-resistant strains have mutations in a small region defined “rifampicin resistance-determining region” in Mtb-RNAP (Morlock et al., 2000; Zaw et al., 2018). The most common mutation in rifampicin resistance-determining region are S456, H451, and D441, corresponding to S531, H526, and D516 in Escherichia coli, respectively. Studies have shown that 70% of rifampicin-resistant clinical isolates have point mutation in two residues (S456 and H451) (Morlock et al., 2000), and H451 is most usually substituted for Asp (D), Tyr (Y), and Arg (R) (Telenti et al., 1993; Caws et al., 2006; Ma et al., 2006; Chikaonda et al., 2017; Wu and Hilliker, 2017). As early as 1995, the in vitro activity experiment of rifampicin by Bodmer et al. (1995) had been demonstrated that H451D/Y/R mutations could cause high-level resistance to rifampicin. After more than 20 years, the level and frequency of resistance to rifampicin are also increasing.
In 2017, there is about 558,000 new cases of rifampicin-resistant tuberculosis (RR-TB), of which 82% are MDR-TB and about 230,000 deaths from MDR/RR-TB (Organization, 2018). Currently, although MDR/RR-TB can be cured with the second-line drugs (e.g., fluoroquinolone and an injectable aminoglycoside), poor efficiency, high toxicity, and expensive price of these drugs make it still difficult for many MDR-TB patients. In some cases, more severe extensively drug-resistant TB may occur, and it will not respond to the most effective second-line anti-TB drugs (Sotgiu et al., 2015; Jeon, 2017; Tiberi et al., 2018). Obviously, the development of new anti-TB drugs is urgent, and exploring the resistance mechanism of rifampicin is of great significance for the discovery of effective drugs.
In this work, in order to uncover the resistance mechanism of Mtb to rifampicin due to the mutation of Mtb-RNAP at position 451, three independent molecular dynamics (MD) simulations for the wild-type Mtb-RNAP and H451D/Y/R mutants were carried out. Based on the obtained trajectories, the molecular mechanics generalized-Born surface area (MM-GBSA) method was applied to calculate the binding free energy of rifampicin with Mtb-RNAP. Furthermore, dynamic network analysis combined with residue interaction network (RIN) analysis was used to show the detailed changes of interactions among the residues surrounding the binding pocket. With the structural and energy analysis, a possible rifampicin-resistant mechanism was also proposed. Compared with the traditional experimental method, MD simulations can show the intuitive and dynamics interaction change process between rifampicin and Mtb-RNAP due to the point mutation. Together with the energy analysis and the dynamics network analysis, the present study show the essential reason of Mtb-RNAP resistant to rifampicin, which can provide the useful guidance for the further drug design against drug resistance.
Materials and Methods
Systems Preparation
The initial atomic coordinate of the wild-type Mtb-RNAP with rifampicin was obtained from Protein Data Bank (Protein Data Bank ID: 5UHB). The crystal structure of Mtb-RNAP reported by Lin et al. (2017) reveals that Mtb-RNAP is composed of six chains, for the A, B chains encoded by the rpoA gene, the C chain encoded by the rpoB gene (Miller et al., 1994), and the D, E, F chains encoded by rpoC, rpoZ, and rpoD, respectively. Rifampicin binds at the active site of the C chain (shown in Figure 1) and inhibits the DNA-directed RNA synthesis of Mtb (McClure and Cech, 1978; Campbell et al., 2001; Somoskovi et al., 2001). Considering that the speed to simulate the whole Mtb-RNAP (~3,826 residues) is too slow, only the C chain complexed with rifampicin was extracted and used as the initial structure of simulations. Furthermore, the deletion of other chains will make the residues of the interface between the two chains unstable, which is inconsistent with that in the multimer. Thus, to simulate the state of interface in the multimer, some relatively flexible and far from the active site amino acid residues were deleted. The three-dimensional structures of three mutants (H451D/Y/R) were obtained by mutating H451 residue in wild type.
Figure 1. (A) The structural overview for Mycobacterium tuberculosis RNA polymerase bound with rifampicin (the protein is shown as cartoon, rifampicin is shown as a stick, the binding pocket of rifampicin is shown as surface). (B) The molecular scheme of the rifampicin.
To generate the force field parameters for the ligand, the Gaussian 09 program (Frisch et al., 2009) was used to optimize the structure of rifampicin and calculate the electrostatic potential at the Hartree–Fock level with 6-31G* basis set. Then, the restraint electrostatic potential protocol (Bayly et al., 1993; Cieplak et al., 1995; Fox and Kollman, 1998) was employed to fit the atomic partial charges. The general amber force field (gaff) (Wang et al., 2004) generated by the antechamber program in the Amber14 package (Case et al., 2014) was applied to describe the ligand. The standard ff99SB force field (Hornak et al., 2006) was used to describe the protein. Then, the LEaP module was applied to add all missing hydrogen atoms and a certain amount of sodium counter-ions to neutralize the unbalanced charges and maintain the systems electro-neutrality. Finally, a rectangular periodic water box of TIP3P (Jorgensen et al., 1983) was added to each system with the water molecules extended 10-Å distance around the complex. The size of the periodic boundary box is 92.6 × 114.4 × 114.9 Å. The whole system has a total of ~100,000 atoms per periodic cell.
Molecular Dynamics Simulations
All MD simulations were performed with Amber14 package (Case et al., 2014). The process of energy minimization, heating, and equilibration was carried out with the Particle Mesh Ewald Molecular Dynamics module. Initially, the energy minimization of each solvated complex includes three steps. For each step, energy minimization was carried out by the steepest descent method for the first 2,500 steps and conjugated gradient method for the subsequent 2,500 steps. In the first step, all the atoms of the complex were restrained with a force constant of 2.0 kcal/(mol·Å2) to only minimize the solvent and ion molecules. After that, the protein backbone atoms were fixed with a restraint force of 2.0 kcal/(mol·Å2) in the second step. Finally, all atoms in the system were minimized without any restraint. After energy minimization, all systems were heated up from 0 to 310 K in the canonical (NVT) ensemble over 100 ps by restraining the protein and the ligand with a 2.0 kcal/(mol·Å2) force constant and using a Langevin thermostat with a coupling coefficient of 2.0/ps. After heating, five steps MD pre-equilibration at 310 K were performed in the NPT ensemble by restraining all the atoms of the complex with decreasing restraints from 2.0 to 1.5, to 1.0, to 0.5, to 0.1. Then, 50-ns equilibration MD simulation without any restraints was performed to eliminate collisions between atoms. Finally, 500-ns production MD simulation was carried out without any restraints on each system in the NPT ensemble at the temperature of 310.0 K and pressure of 1 atm. During the simulations, all the bonds involving hydrogen atoms were restrained with SHAKE algorithm (Ryckaert et al., 1977) to avoid too fast vibration of the hydrogen atoms. In addition, periodic boundary conditions were employed, and the long-range Coulombic interactions were treated using the Particle Mesh Ewald (Darden et al., 1993). The time step was set to 2 fs.
Molecular Mechanics Generalized-Born Surface Area Calculation
To uncover the effects of mutations on the binding affinity of rifampicin to Mtb-RNAP, the MM-GBSA method was applied to estimate the binding free energy of each complex, which has been successfully used in a lot of researches (Pan et al., 2011; Yang et al., 2012). Here, we extracted 1,000 snapshots at 100-ps interval from the last 100-ns trajectory for each system. The binding free energy was calculated from the equation:
where Gcomplex, Greceptor, and Gligand are the free energy of complex, protein, and ligand, respectively. The free energy for each molecular species was calculated based on an average over the extracted snapshots. Each of them can be estimate with the following equations:
where EMM is the gas-phase energy calculated using the Amber ff03 molecular mechanics force field. Eint is the internal energy, including the energy of bond (Ebond), angle (Eangle), and torsion (Etorsion). Eele and Evdw are the Coulomb and van der Waals energy, respectively. Gsol is the solvation free energy and can be decomposed into polar solvation free energy (GGB) and non-polar solvation free energy (GSA). GGB was calculated by solving the GB equation and the dielectric constants for solute as well as solvent were set to 1.0 and 80.0, respectively (Rocchia et al., 2001). GSA was estimated by the solvent accessible surface area determined using a water probe radius of 1.4 Å. The surface tension constant γ was set to 0.0072 kcal/(mol·Å2), and the non-polar contribution to the solvation free energy term β was set to 0 (Sitkoff et al., 1994). T and S are the temperature and the total solute entropy, respectively. The entropy contributions can be estimated by normal mode analysis (Pearlman et al., 1995). However, here, we did not calculate the entropy contributions since our aim is not to obtain the absolute Gibbs energy but to identify the key residues of binding pocket and the detailed interaction features. In addition, previous studies have proven that it is sufficient to compare the binding ability of receptors and ligands based on the values of enthalpy changes (ΔHbind) (Aruksakunwong et al., 2006; Xue et al., 2012).
Moreover, in order to identify the key residues responsible for the binding of rifampicin, the MM-GBSA binding free energy decomposition process was used to decompose the interaction energy to each residue by considering molecular mechanics and solvation energy without considering the contribution of entropy.
Dynamic Network Analysis
Dynamic network analysis, as an effective method to extract information from the obtained molecular dynamics trajectories, has been successfully applied in protein misfolding (Zhou et al., 2019) and protein–protein interaction analysis (Sethi et al., 2009; Bai et al., 2014). Here, in order to observe the dynamic changes of the residues interaction network, 2,000 snapshots were extracted from the last 50-ns trajectory for each system. In the network, one node represents one residue, and the position of each node is defined at the center of Cα atom of residue. The edge represents the interactions of two residues. Furthermore, the edge weight (Wij) between two nodes (i, j) was defined with the following equation:
where Cij represents the pairwise correlations, which is calculated by Carma program (Glykos, 2006), a plugin in VMD (Humphrey et al., 1996). Finally, the NetworkView module in VMD was used to visualize the residue interaction network.
Residue Interaction Network Analysis
RIN uses a network diagram to simplify the inter-residue interaction, which considers the residues as nodes and physico–chemical interactions such as covalent and non-covalent bonds as edges. RIN method has been successfully used to analyze the effects of mutations on drug resistance (Xue et al., 2012, 2013, 2014). In this work, the residue interaction network generator 2.0 (RING-2.0) (Piovesan et al., 2016) software was applied to generate the network for the representative structures. The calculation process of RING-2.0 is described as follows: (i) the calculation of the secondary structure elements by incorporating the DSSP algorithm (Kabsch and Sander, 1983); (ii) hydrogen atom placement based on geometric criteria; (iii) hydrogen bond calculation; and (iv) the calculation of van der Waals interactions. Moreover, Cytoscape (Shannon et al., 2003) and the plugin RINalyzer (Doncheva et al., 2011) were used to visualize the residue interaction network.
Results
H451D/Y/R Mutations Increased the Flexibility of the Active Pocket
Firstly, the root-mean-square-deviations (RMSD) value for the protein backbone atoms, the active pocket, and the heavy atoms of rifampicin relative to the initial structure were calculated to monitor the equilibrium of each system. As shown in Figure 2 and Supplementary Figures 1–3, three parallel MD simulations have similar fluctuations, suggesting each parallel trajectory can produce reproducible results. Thus, the following analysis was based on one of three parallel MD simulations. As can be seen from the RMSDs of wild-type Mtb-RNAP and three mutants, each system achieves equilibrium after 100 ns. Therefore, the last 100-ns trajectory for each system was used for the following structural and energetic analysis. Additionally, from the monitoring of the RMSD value of the heavy atoms of ligand, we can justify roughly if the ligand can bind to the target stably. From Figure 2C and Supplementary Figures 1–3 the RMSDs of rifampicin in mutants were larger than that in wild type, indicating that rifampicin had a large fluctuation in mutants.
Figure 2. Root-mean-square-deviations (RMSDs) for the wild-type system of three independent molecular dynamics simulations: Sim1 (black), Sim2 (green), Sim3 (orange): (A) RMSDs for the backbone atoms of protein versus time. (B) RMSDs for the backbone atoms of active pocket vs. time. (C) RMSDs for the heavy atoms of rifampicin vs. time. (D) Root-mean-square-fluctuation for the backbone atoms of active pocket vs. residue number of wild type (black), H451D (red), H451Y (green), and H451R (purple) in Sim1.
The root-mean-square-fluctuation (RMSF) of each residue was calculated based on the last 100-ns trajectory for each system, and the corresponding results were shown in Figure 2D. It can be seen from Figure 2D that the RMSF values of the residues have similar trends for all systems. However, H451D/Y/R mutations will increase the RMSF values relative to the wild-type Mtb-RNAP, which indicates that these three mutations increased the flexibility of the binding pocket and weakened the interaction of rifampicin with Mtb-RNAP.
H451D/Y/R Mutations Weaken the Binding Ability of Mycobacterium tuberculosis RNA Polymerase With Rifampicin
To explore the effects of three mutations on the binding of Mtb-RNAP with rifampicin, the binding free energy calculation was performed based on MM-GBSA method. As shown in Table 1, the enthalpy changes (ΔHbind) of the wild-type Mtb-RNAP and H451D/Y/R mutants with rifampicin are −43.89, −31.20, −35.55, and −28.58 kcal/mol, respectively. As expected, the binding affinity of Mtb-RNAP to rifampicin reduced obviously due to H451D/Y/R mutations.
Table 1. The calculated binding free energy and the detailed contribution of different energy terms (kcal/mol).
By assessing the contributions of individual energy terms, we found that the non-polar interactions (sum of van der Waals interaction energy ΔEvdw and non-polar interaction energy in solvation free energy ΔGSA) are the driving force for the binding of rifampicin. However, the energy contributions of Evdw and ΔGSA decrease in the mutants. Relative to non-polar interactions, the polar interaction (sum of electrostatic interaction energy ΔEele and polar interaction energy in solvation free energy ΔGGB) seems like unfavorable for the binding of rifampicin and more apparent in the mutants (15.22, 19.21, 20.70, and 23.15 kcal/mol for the wild type, H451D, H451Y, and H451R mutants, respectively). Although the contributions of intermolecular electrostatic interactions (ΔEele) are very favorable, their contributions are compensated by the large desolvation penalties.
Key Residues Responsible for the Reduced Binding Ability of H451D/Y/R Mutants
By decomposing the binding free energy of the wild-type Mtb-RNAP with rifampicin, 10 key residues such as Q438, F439, D441, R454, P489, E490, N493, I497, R613, and Q614 (Figure 3A) with energy contributions over 2 kcal/mol are identified. It can be seen that relative to the wild-type Mtb-RNAP, the contributions of Q438, F439, D441, E490, N493, R613, and Q614 have an obvious decrease in H451D mutant. For H451Y mutant, the reduced energy contribution of residues Q438, D441, R454, N493, R613, and Q614 should be responsible for the reduced binding affinity of rifampicin to Mtb-RNAP. Moreover, the energy reduction in H451R mutant is more obvious relative to that in H451D/Y mutants. Actually, the profile of each residue's energy contribution in three mutants shares some similar features. For example, in all three mutants, the residues Q438, D441, N493, R613, and Q614 have obvious reduced contribution for the binding of rifampicin, suggesting that the drug-resistance mechanisms due to three mutations have some similarities.
Figure 3. (A) The total energy contributions of the key residues for binding of rifampicin: wild type (black), H451D (red), H451Y (green), H451R (purple). (B) The electrostatic energy contributions of the key residues for binding of rifampicin. (C) The van der Waals energy contributions of the key residues for binding of rifampicin.
To further explore the origin of the reduced residues' contribution, we compared the electrostatic contribution and van der Waals contribution of the key residues (shown in Figures 3B,C, respectively). From Figure 3B, we can obtain when H451 is mutated to D451, the electrostatic energy contribution of D451 (4.80 kca/mol) is detrimental for rifampicin binding. Moreover, the contributions of electrostatic energy of Q438, F439, E490, and R613 decrease after H451D mutation. For H451Y mutant, the electrostatic energy contribution of three residues (Q438, R454, and R613) also have an obvious reduction. In H451R mutant, the contributions of electrostatic interactions of almost all key residues reduced. In addition, Figure 3C shows that the reduction of energy contribution of D441, P489, E490, N493, and Q614 is mainly from the loss of the van der Waals interaction contribution in H451D/Y/R mutants.
The Dynamic Network Analysis and Residue Interaction Network Analysis Reveal That H451D/Y/R Mutations Weaken the Interaction of Mutated Residue With Its Adjacent Residues
To investigate how H451D/Y/R mutations change the binding pocket and further lead to the resistance of Mtb-RNAP to rifampicin, the dynamic network analysis was further carried out based on the 2,000 snapshots extracted from the equilibrium phase for each system. The obtained results are shown in Figure 4. The strength of the total interactions (including van der Waals interaction and hydrogen bond interaction and so on) is indicated by the edge thickness. Moreover, the type of interactions is shown by the two-dimensional RIN interactively based on the representative structure of each system (Figure 5). Here, the representative structures were extracted by clustering analysis, and the conformation with the lowest RMSD to the cluster center was selected.
Figure 4. The picture of dynamic network analysis for the active pocket of rifampicin binding: (A) WT; (B) H451D; (C) H451Y; (D) H451R; The purple spheres represent the residues and the sticks represent the total interactions. The strength of interactions of two residues is indicted by the thickness degree of stick, the thicker of the stick, and the stronger interaction of two residues.
Figure 5. RIN for the active pocket of rifampicin binding (A) wild type; (B) H451D; (C) H451Y; (D) H451R. The pink octagons represent residues, and the yellow octagons represent the mutated residues; the edges represent van der Waals (blue) and hydrogen bond (red) interactions and salt bridge interaction (cyan).
In the wild-type Mtb-RNAP, H451 residue can be identified as the central node with six first neighbors in the interaction network according to Figures 4A, 5A. H451 not only can form strong van der Waals interactions with R454, Q438, F439, M440, and D441 but also form hydrogen bonds with D441, S447. Moreover, S447 acts as a key node in the network; in addition to forming hydrogen bond with H451, it also can form strong interactions with D441 (van der Waals interaction), R613 (hydrogen bond interaction), and Q614. I497, as another key node in the other side of the network, can form van der Waals interactions with R454, P489, and N493. Finally, R454 as a joint node can connect the two sides of the networks together by forming the van der Waals interactions with I497, H451, and the hydrogen bond with Q438. In addition, the weak “triangular shape” interaction network formed by R454, R613, and P489 residues that make the binding pocket more compact and coherent.
However, in H451D mutant, it is evident that the mutation causes the disappearance of the interactions between D451 and Q438, F439, M440, D441, and S447 mentioned previously (Figures 4B, 5B). For H451Y/R mutants, the corresponding interactions are also obviously weakened (Figures 4C,D). Overall, in H451D/Y/R mutants, the “triangular shape” interaction network is broken. Based on these results, it can be concluded that the mutations on 451 indeed reduced the interaction connection of the residues in the binding pocket and then result in the active pocket more flexible and open.
The Comparison of Binding Modes of Rifampicin With the Wild-Type Mycobacterium tuberculosis RNA Polymerase and Three Mutants
To show the detailed rifampicin-resistance mechanism to Mtb-RNAP H451D/Y/R mutants intuitively, further structural analysis was performed. The representative structures are depicted in Figure 6. As shown in Figures 3C, 6A, the van der Waals interactions of Q438, D441, H451, R454, P489, E490, N493, I497, and Q614 with rifampicin are pivotal for the binding of rifampicin to the wild-type Mtb-RNAP. In addition, the calculation of hydrogen bond occupancy was carried out to monitor the formation of hydrogen bonds between rifampicin and Mtb-RNAP over the whole MD simulations. Based on the result in Table 2, we can see that the O8, O11, O9/O2, and O12 atoms of rifampicin can form stable hydrogen bonds with the side chains of Q438, F439, and R613 with high occupancy rate. The formation of these hydrogen bonds makes F439, R613, and Q438 have large electrostatic energy contributions (Figure 3B).
Figure 6. Three-dimensional representation for the binding mode of rifampicin with wild-type Mycobacterium tuberculosis RNA polymerase and three mutants based on the obtained representative structure: (A) wild type; (B) H451D; (C) H451Y; (D) H451R. The protein is shown as a cyan cartoon; rifampicin and the key residues are represented as hotpink and cyan carbon sticks, respectively; the intermolecular hydrogen bonds are indicated as dashed yellow lines.
Table 2. The occupancy (%) of hydrogen bonds between rifampicin and Mycobacterium tuberculosis RNA polymerase.
By comparing the binding modes of wild-type and H451D mutant (Figures 6A,B), the position of rifampicin has a clear movement. Combined with the results of RIN analysis (Figures 5A,B), H451D mutation causes the interactions of residue 451 between with some residues of binding pocket disappear, further causing the conformations of some amino acids (such as D441, P489, E490, N493, and Q614) that changed a lot. Finally, the residues D441, P489, E490, N493, and Q614 are away from rifampicin, causing the weakened van der Waals interaction of these residues. From Table 2, some hydrogen bonds between rifampicin and Q438, F439, and R613 disappeared in H451D mutant, which causes the reduction of the electrostatic energy contribution (Figure 3B). Moreover, the H451D mutation leads to the electrostatic repulsion between carboxyl group of D451 and O10 atom of rifampicin and is unfavorable for the binding of rifampicin.
For H451Y mutant, the binding mode of rifampicin with Mtb-RNAP is more similar to that in the wild-type Mtb-RNAP, and the position of rifampicin just has a slight movement (Figures 6A,C). Despite this, the mutation still causes some residues (such as N493, I497, and Q614) away from rifampicin, reducing the nonpolar contribution of these residues. In addition, the result in Table 2 shows that the disappearance of the hydrogen bond between R613 and rifampicin is responsible for the reduction of electrostatic contributions of R613. The reduction of hydrogen bond occupancy rate formed between Q438 and rifampicin causes the loss of electrostatic contribution of Q438.
From the energetic and structural analysis discussed previously, it can be seen that the flipping of R613 greatly affects the binding of rifampicin in H451D/Y mutants. Therefore, it is worthy to explore how the mutation of H451 affects the flipping of R613. For this aim, we superimpose the binding pocket of the wild-type Mtb-RNAP with those of H451D/Y mutants. Figure 7 shows the superimposition results. From Figure 7A, it can be seen that the position of H451 is on H1. When H451 is mutated to D451, no hydrogen bond can be formed between D451 and S447 (Figure 5B), which causes the helix structure of S447 to transform into disordered loops. At the same time, the fluctuation of this loop further causes the disappearance of the hydrogen bond between R613 (located on H2) and S447 (located on H1). Therefore, the R613 residue flips with H2 due to the weaker interaction between H1 and H2. However, the mutation of H451Y enhances the van der Waals interaction of Y451 and S447 (Figure 5C), which increases the distance of S447 and R613 and further to interfere the formation of hydrogen bond. Moreover, the rotation of H2 makes the flipping of R613 more obvious. The flipping of R613 directly causes the disappearance of the hydrogen bonds with rifampicin, further leading that the electrostatic contribution of R163 reduced.
Figure 7. (A) The superimposition of wild-type Mycobacterium tuberculosis RNA polymerase and H451D mutant. (B) The superimposition of wild-type Mycobacterium tuberculosis RNA polymerase and H451Y mutant; the protein is shown as cartoon with cylindrical helices colored with cyan (wild type) and magenta (mutant); rifampicin and R613, S447, H451, D451, and Y451 are shown as stick colored with cyan (wild type) and magenta (mutant); the intermolecular hydrogen bond between S447 and R613 is indicated as dashed yellow lines.
For H451R mutant, there was a large conformational transition of the binding pocket and rifampicin as shown in Figure 6D. Moreover, R451 has a steric clash with rifampicin and so that rifampicin moves out of the binding pocket. Thereby, the movement of rifampicin causes the energy contribution of some key residues such as R613, E490, S447, D441, M440, and F439 reduce (Figure 3). In addition, some hydrogen bonds between rifampicin and F439/R613 also disappear (Table 2). In order to observe the changes of the binding pocket more intuitively, the surface map for the binding pocket of the wild-type and H451R mutant is depicted in Figure 8. It can be seen that when H451 is replaced by R451, the long side chain of R451 made the binding pocket smaller and cannot accommodate the relatively rigid rifampicin, which leads to the movement of rifampicin.
Figure 8. The binding pocket of rifampicin: (A) wild type; (B) H451R mutant; the pocket is shown as surface, and rifampicin is shown as stick colored with cyan.
In summary, there are three main binding sites between rifampicin and the active pocket in wild type (Figure 9A). The polar pocket formed by the residues N610, R613, Q614, etc. (S1) can act as the hydrogen bond acceptor to form the hydrogen bond with the oxygen atom of rifampicin. The other polar pocket (S2) consists of residues Q438, M440, D441, etc., which also forms stable hydrogen bonds with rifampicin. In addition, the hydrophobic pocket composed of residues such as P489, N493, and I497 just accommodates the naphthalene ring of rifampicin. However, H451D/Y/R mutants changed the initial binding mode for rifampicin resulting from the change of the side chain size of 451. The N atom from imidazole ring of H451 can form the hydrogen bond interaction with rifampicin in wild-type RNA polymerase, which is interfered by these mutations. Furthermore, the mutations in 451 position impaired the interaction between the residue 451 and other key residues in the active pocket, interfering with the specific orientation of the key residues' side chains, thereby increasing the flexibility of the key residues such as R613. As a result, the binding affinity of rifampicin reduced, and the important interactions between rifampicin and active pocket were disturbed. One of the most apparent changes from the results of binding free energy decomposition was that the reduction of energy contribution from R613, which was mainly because the disappearance of the hydrogen bond between R613 and O atom at the C15 position of rifampicin (Figure 9B). Therefore, according to the obtained mechanism, to overcome the drug resistance induced by H451D/Y/R, one possible strategy is to enhance the interaction of the inhibitor with R613 by replacing the carbonyl group at the C15 position in rifampicin with a longer and negatively charged group (R). Such group may recover the hydrogen bond interaction between R613 and inhibitor, which could stabilize the binding of rifampicin.
Figure 9. (A) Three main binding sites between rifampicin and the active pocket. Rifampicin is shown as stick colored with magenta, and the residues are shown as lines colored with cyan. The yellow dotted lines represent the hydrogen bond interactions. (B) Structural optimization of rifampicin. R represents a longer negatively charged group.
Conclusion
In this study, MD simulations together with MM-GBSA, dynamic network analysis, and RIN analysis were carried out on the complexes of rifampicin with the wild-type Mtb-RNAP and H451D/Y/R mutants to explore the resistant mechanism of rifampicin. The results from the MM-GBSA calculations are well-consistent with the experimental result. The reduced binding affinity for the studied mutants mainly comes from the loss of van der Waals contribution of D441, P489, E490, N493, I497, and Q614 and electrostatic contribution of Q438, F439, and R613 in mutants. The electrostatic energy contribution of R613 decreases obviously with the disappearance of the hydrogen bonds between R613 and rifampicin, which caused by the conformation flipping of R613 in H451D/Y mutants. The binding modes and dynamic network analysis show that the weakened interactions among D/Y/R451 with Q438, F439, M440, D441, and S447 increase the flexibility of the binding pocket, thereby reducing the binding affinity of rifampicin to Mtb-RNAP. In addition, these mutations caused the key hydrogen bond interactions between residue 451 and rifampicin disappear. Finally, the position of rifampicin had a clear movement, which changed the stable binding mode of rifampicin. We firstly proposed the atomic level resistance mechanism of Mtb to rifampicin due to H451D/Y/R mutations on Mtb-RNAP. In addition, we also proposed some guidance for the alteration of the rifampicin and the development of new drugs in the future. Though the hypothesis is still unvalidated on the traditional experimental, it can provide some theoretical underpinnings for the design of new anti-TB drugs to some extent, or a particular aspect.
Data Availability Statement
All datasets generated for this study are included in the article/Supplementary Material.
Author Contributions
QZ and HuL designed the research and wrote the manuscript. QZ, XA, HoL, SW, and TX performed the dynamic simulations and analyzed the data.
Funding
This work was supported by the National Nature Science Foundation of China (Grant No. 21675070).
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/fchem.2019.00819/full#supplementary-material
References
Ahmad, S., Mokaddas, E., and Al-Mutairi, N. M. (2018). Prevalence of tuberculosis and multidrug resistant tuberculosis in the Middle East region. Exp. Rev. Anti. Infect. Ther. 16, 709–721. doi: 10.1080/14787210.2018.1519393
Aruksakunwong, O., Wolschann, P., Hannongbua, S., and Sompornpisut, P. (2006). Molecular dynamic and free energy studies of primary resistance mutations in HIV-1 protease-ritonavir complexes. J. Chem. Inf. Model. 46, 2085–2092. doi: 10.1021/ci060090c
Bai, Q., Shen, Y., Jin, N., Liu, H., and Yao, X. (2014). Molecular modeling study on the dynamical structural features of human smoothened receptor and binding mechanism of antagonist LY2940680 by metadynamics simulation and free energy calculation. Biochim. Biophys. Acta 1840, 2128–2138. doi: 10.1016/j.bbagen.2014.03.010
Bayly, C. I., Cieplak, P., Cornell, W., and Kollman, P. A. (1993). A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys. Chem. 97, 10269–10280. doi: 10.1021/j100142a004
Bodmer, T., Zürcher, G., Imboden, P., and Telenti, A. (1995). Mutation position and type of substitution in the β-subunit of the RNA polymerase influence in-vitro activity of rifamycins in rifampicin-resistant Mycobacterium tuberculosis. J. Antimicrob. Chemother. 35, 345–348. doi: 10.1093/jac/35.2.345
Campbell, E. A., Korzheva, N., Mustaev, A., Murakami, K., Nair, S., Goldfarb, A., et al. (2001). Structural mechanism for rifampicin inhibition of bacterial RNA polymerase. Cell 104, 901–912. doi: 10.1016/S0092-8674(01)00286-0
Case, D. A., Babin, V., Berryman, J. T., Betz, R. M., Cai, Q., Cerutti, D. S., et al. (2014). AMBER 14. San Francisco, CA: University of California.
Caws, M., Duy, P. M., Tho, D. Q., Lan, N. T., Hoa, D. V., and Farrar, J. (2006). Mutations prevalent among rifampin-and isoniazid-resistant Mycobacterium tuberculosis isolates from a hospital in Vietnam. J. Clin. Microbiol. 44, 2333–2337. doi: 10.1128/JCM.00330-06
Chang, K. C., Leung, C. C., Yew, W. W., Chan, S. L., and Tam, C. M. (2006). Dosing schedules of 6-month regimens and relapse for pulmonary tuberculosis. Am. J. Respir. Crit. Care Med. 174, 1153–1158. doi: 10.1164/rccm.200605-637OC
Chang, K. C., Nuermberger, E., Sotgiu, G., and Leung, C. C. (2018). New drugs and regimens for tuberculosis. Respirology 23, 978–990. doi: 10.1111/resp.13345
Chikaonda, T., Ketseoglou, I., Nguluwe, N., Krysiak, R., Thengolose, I., Nyakwawa, F., et al. (2017). Molecular characterisation of rifampicin-resistant Mycobacterium tuberculosis strains from Malawi. Afr. J. Lab. Med. 6:463. doi: 10.4102/ajlm.v6i2.463
Cieplak, P., Cornell, W. D., Bayly, C., and Kollman, P. A. (1995). Application of the multimolecule and multiconformational RESP methodology to biopolymers: charge derivation for DNA, RNA, and proteins. J. Comput. Chem. 16, 1357–1377. doi: 10.1002/jcc.540161106
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
Doncheva, N. T., Klein, K., Domingues, F. S., and Albrecht, M. (2011). Analyzing and visualizing residue networks of protein structures. Trends Biochem. Sci. 36, 179–182. doi: 10.1016/j.tibs.2011.01.002
Fox, T., and Kollman, P. A. (1998). Application of the RESP methodology in the parametrization of organic solvents. J. Phys. Chem. 102, 8070–8079. doi: 10.1021/jp9717655
Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., et al. (2009). Gaussian 09. Wallingford, CT: Gaussian Inc.
Glykos, N. M. (2006). Software news and updates carma: a molecular dynamics analysis program. J. Comput. Chem. 27, 1765–1768. doi: 10.1002/jcc.20482
He, G. X., Zhao, Y. L., Jiang, G. L., Liu, Y. H., Xia, H., Wang, S. F., et al. (2008). Prevalence of tuberculosis drug resistance in 10 provinces of China. BMC Infect. Dis. 8:166. doi: 10.1186/1471-2334-8-166
Hornak, V., Abel, R., Okur, A., Strockbine, B., Roitberg, A., and Simmerling, C. (2006). Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins 65, 712–725. doi: 10.1002/prot.21123
Humphrey, W., Dalke, A., and Schulten, K. (1996). VMD: visual molecular dynamics. J. Mol. Graph. 14, 33–38. doi: 10.1016/0263-7855(96)00018-5
Jeon, D. (2017). WHO treatment guidelines for drug-resistant tuberculosis, 2016 update: applicability in South Korea. Tuberc. Respir. Dis. 80, 336–343. doi: 10.4046/trd.2017.0049
Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. (1983). Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935. doi: 10.1063/1.445869
Kabsch, W., and Sander, C. (1983). Dictionary of protein secondary structure: pattern recognition of hydrogen-bonded and geometrical features. Biopolymers 22, 2577–2637. doi: 10.1002/bip.360221211
Lin, W., Mandal, S., Degen, D., Liu, Y., Ebright, Y. W., Li, S., et al. (2017). Structural basis of Mycobacterium tuberculosis transcription and transcription inhibition. Mol. Cell 66:e168. doi: 10.1016/j.molcel.2017.03.001
Ma, X., Wang, H., Deng, Y., Liu, Z., Xu, Y., Pan, X., et al. (2006). rpoB Gene mutations and molecular characterization of rifampin-resistant Mycobacterium tuberculosis isolates from Shandong Province, China. J. Clin. Microbiol. 44, 3409–3412. doi: 10.1128/JCM.00515-06
McClure, W. R., and Cech, C. L. (1978). On the mechanism of rifampicin inhibition of RNA synthesis. J. Biol. Chem. 253, 8949–8956.
Miller, L. P., Crawford, J. T., and Shinnick, T. M. (1994). The rpoB gene of Mycobacterium tuberculosis. Antimicrob. Agents Chemother. 38, 805–811. doi: 10.1128/AAC.38.4.805
Mishra, A., and Surolia, A. (2018). Mycobacterium tuberculosis: surviving and indulging in an unwelcoming host. IUBMB Life 70, 917–925. doi: 10.1002/iub.1882
Mitchison, D. A., and Nunn, A. J. (1986). Influence of initial drug resistance on the response to short-course chemotherapy of pulmonary tuberculosis. Am. Rev. Respir. Dis. 133, 423–430.
Morlock, G. P., Plikaytis, B. B., and Crawford, J. T. (2000). Characterization of spontaneous, in vitro-selected, rifampin-resistant mutants of Mycobacterium tuberculosisstrain H37Rv. Antimicrob. Agents Chemother. 44, 3298–3301. doi: 10.1128/AAC.44.12.3298-3301.2000
Pan, D., Sun, H., Bai, C., Shen, Y., Jin, N., Liu, H., et al. (2011). Prediction of zanamivir efficiency over the possible 2009 Influenza A (H1N1) mutants by multiple molecular dynamics simulations and free energy calculations. J. Mol. Model. 17, 2465–2473. doi: 10.1007/s00894-010-0929-8
Pearlman, D. A., Case, D. A., Caldwell, J. W., Ross, W. S., Cheatham, I. I. I., DeBolt, T. E., et al. (1995). AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules. Comput. Phys. Commun. 91, 1–41. doi: 10.1016/0010-4655(95)00041-D
Piovesan, D., Minervini, G., and Tosatto, S. C. (2016). The RING 2.0 web server for high quality residue interaction networks. Nucleic Acids Res. 44, W367–W374. doi: 10.1093/nar/gkw315
Rocchia, W., Alexov, E., and Honig, B. (2001). Extending the applicability of the nonlinear Poisson– Boltzmann equation: multiple dielectric constants and multivalent ions. J. Phys. Chem. B 105, 6507–6514. doi: 10.1021/jp010454y
Ryckaert, J. P., Ciccotti, G., and Berendsen, H. J. C. (1977). Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23, 327–341. doi: 10.1016/0021-9991(77)90098-5
Sandgren, A., Strong, M., Muthukrishnan, P., Weiner, B. K., Church, G. M., and Murray, M. B. (2009). Tuberculosis drug resistance mutation database. PLoS Med. 6:e1000002. doi: 10.1371/journal.pmed.1000002
Seid, M. A., Ayalew, M. B., Muche, E. A., Gebreyohannes, E. A., and Abegaz, T. M. (2018). Drug-susceptible tuberculosis treatment success and associated factors in Ethiopia from 2005 to 2017: a systematic review and meta-analysis. BMJ Open 8:e022111. doi: 10.1136/bmjopen-2018-022111
Service, H. K. C., and Council, B. M. R. (1981). Controlled trial of four thrice-weekly regimens and a daily regimen all given for 6 months for pulmonary tuberculosis. Lancet 1, 171–174.
Sethi, A., Eargle, J., Black, A. A., and Luthey-Schulten, Z. (2009). Dynamical networks in tRNA: protein complexes. Proc. Natl. Acad. Sci. U.S.A. 106, 6620–6625. doi: 10.1073/pnas.0810961106
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303
Sitkoff, D., Sharp, K. A., and Honig, B. (1994). Accurate calculation of hydration free energies using macroscopic solvent models. J. Phys. Chem. 98, 1978–1988. doi: 10.1021/j100058a043
Somoskovi, A., Parsons, L. M., and Salfinger, M. (2001). The molecular basis of resistance to isoniazid, rifampin, and pyrazinamide in Mycobacterium tuberculosis. Respir. Res. 2, 164–168. doi: 10.1186/rr54
Sotgiu, G., D'Ambrosio, L., Centis, R., Mura, I., Castiglia, P., Spanevello, A., et al. (2015). The multidrug-resistant tuberculosis threat: old problems and new solutions. J. Thorac. Dis. 7, E354–E360. doi: 10.3978/j.issn.2072-1439.2015.09.21
Telenti, A., Imboden, P., Marchesi, F., Lowrie, D., Cole, S., Colston, M. J., et al. (1993). Detection of rifampicin-resistance mutations in Mycobacterium tuberculosis. Lancet 341, 647–650. doi: 10.1016/0140-6736(93)90417-F
Tiberi, S., du Plessis, N., Walzl, G., Vjecha, M. J., Rao, M., Ntoumi, F., et al. (2018). Tuberculosis: progress and advances in development of new drugs, treatment regimens, and host-directed therapies. Lancet Infect. Dis. 18, e183–e198. doi: 10.1016/S1473-3099(18)30110-5
Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A., and Case, D. A. (2004). Development and testing of a general amber force field. J. Comput. Chem. 25, 1157–1174. doi: 10.1002/jcc.20035
Wu, E. Y., and Hilliker, A. K. (2017). Identification of rifampicin resistance mutations in Escherichia coli, including an unusual deletion mutation. J. Mol. Microbiol. Biotechnol. 27, 356–362. doi: 10.1159/000484246
Xue, W., Ban, Y., Liu, H., and Yao, X. (2014). Computational study on the drug resistance mechanism against HCV NS3/4A protease inhibitors vaniprevir and MK-5172 by the combination use of molecular dynamics simulation, residue interaction network, and substrate envelope analysis. J. Chem. Inf. Model. 54, 621–633. doi: 10.1021/ci400060j
Xue, W., Jin, X., Ning, L., Wang, M., Liu, H., and Yao, X. (2013). Exploring the molecular mechanism of cross-resistance to HIV-1 integrase strand transfer inhibitors by molecular dynamics simulation and residue interaction network analysis. J. Chem. Inf. Model. 53, 210–222. doi: 10.1021/ci300541c
Xue, W., Liu, H., and Yao, X. (2012). Molecular mechanism of HIV-1 integrase-vDNA interactions and strand transfer inhibitor action: a molecular modeling perspective. J. Comput. Chem. 33, 527–536. doi: 10.1002/jcc.22887
Yang, Y., Shen, Y., Li, S., Jin, N., Liu, H., and Yao, X. (2012). Molecular dynamics and free energy studies on Aurora kinase A and its mutant bound with MLN8054: insight into molecular mechanism of subtype selectivity. Mol. Biosyst. 8, 3049–3060. doi: 10.1039/c2mb25217a
Zaw, M. T., Emran, N. A., and Lin, Z. (2018). Mutations inside rifampicin-resistance determining region of rpoB gene associated with rifampicin-resistance in Mycobacterium tuberculosis. J. Infect. Public Health 11, 605–610. doi: 10.1016/j.jiph.2018.04.005
Keywords: tuberculosis, rifampicin, drug resistance, molecular dynamics simulation, molecular mechanics generalized-Born surface area, dynamic network analysis
Citation: Zhang Q, An X, Liu H, Wang S, Xiao T and Liu H (2019) Uncovering the Resistance Mechanism of Mycobacterium tuberculosis to Rifampicin Due to RNA Polymerase H451D/Y/R Mutations From Computational Perspective. Front. Chem. 7:819. doi: 10.3389/fchem.2019.00819
Received: 10 July 2019; Accepted: 12 November 2019;
Published: 03 December 2019.
Edited by:
Zunnan Huang, Guangdong Medical University, ChinaReviewed by:
Bing He, University of Michigan, United StatesGuohui Li, Dalian Institute of Chemical Physics (CAS), China
Copyright © 2019 Zhang, An, Liu, Wang, Xiao and Liu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Huanxiang Liu, hxliu@lzu.edu.cn