The MSM reveals a jump-from-cavity mechanism of the PPi release

To fully sample the PPi release process, we first conducted steered MD (SMD) simulations along five potential pathways to pull the PPi group from inside around the active site toward the outside into the bulk water. To eliminate the bias introduced from the SMD simulation, the conformations sampled along the pulling pathways were grouped into tens of clusters. From each cluster, we then performed nanosecond MD simulations starting from centered conformations (see Methods). A total of ~100 MD trajectories were collected with a sum of ~1,000,000 snapshots aligned, so that the RMSDs of the PPi group were calculated to track the PPi location. Alternatively (see Supporting Information and S1A Fig), one can also measure the distance between the PPi group and a Mg²⁺ ion (MgA) kept inside the active site.

We firstly mapped the sampled population density along the PPi RMSD coordinate. As shown in Fig 2, one can identify three population states on the profile, centered at ~ 3 Å, 6 Å, and 10 Å along the RMSD coordinate, and labeled as S1a, S1b, and S2, respectively. For each state, we provide one representative snapshot. In S1a, the PPi group is grabbed by three positively charged residues: Lys631, Arg627, and Lys472, as well as by two negatively charged residues: Asp537 and Asp812, through a Mg²⁺ ion (MgB) bound with PPi. These two aspartic acids are highly conserved in RNA polymerases, and play essential roles in catalysis [48,49]. Separation of PPi from these two aspartic acids is necessary for the PPi release, and this separation becomes visible in S1b. We also projected the sampled conformations along the PPi RMSD and the distance between the PPi-bound MgB and one oxygen of Asp537 (see S1B Fig). While in S2, the PPi group is not only separated far from Asp537 and Asp812, but also loses the contact with Lys631. In addition, the PPi group in S2 forms marginally stable interactions with Lys472 and Arg627, which locate at the exit to the bulk water (see Fig 2). We also conducted convergence tests for the sampled population distributions (see S2 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. The population density profile for conformations sampled along the simulated PPi release pathways.

The profile is shown on a potential-like surface along the PPi RMSD coordinate. Three population states are identified: S1a, S1b, and S2, with representative snapshots shown for each state. The color scheme in the structure is the same as that in Fig 1. The convergence tests of the profile is provided in SI S2 Fig.

https://doi.org/10.1371/journal.pcbi.1004624.g002

We next constructed the MSM based on all the conformations sampled from the short MD simulations by firstly clustering these conformations into hundreds of microstates. We then further lumped these microstates into three macrostates in order to visualize key metastable states (see Methods and Fig 3A and 3B). The three states match well with the above mapping onto S1a (~53%), S1b (~39%), and S2 (~8%). Hence, one can see that the total ‘inside’ conformations (S1a and S1b, together as S1) dominate over the ‘outside’ one (S2). Correspondingly, the PPi release from the inside to the outside is not energetically favored. Instead, the release is an activation process, and the rate-limiting step is the ‘jump’ transition from the inside to the outside. The mean first passage times (MFPTs) of the macrostate transitions were estimated ~ 0.2 μs with a standard error < 0.01 μs for S1a → S1b, and ~ 1.0 μs with a standard error < 0.1 μs for S1b → S2 (the values vary slightly according to different numbers of microstates in the MSM construction, see Methods and SI), while the upper bounds of the MFPTs for S1a → S1b and S1b → S2 were estimated at ~9 μs and ~ 20 μs, respectively (see Methods).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. The MSM construction for the PPi release in the T7 RNAP elongation.

(A) A structural view of the PPi distribution (shown in spheres) sampled from the MD simulations. Each sphere represents the center of the mass of the PPi group. The distribution was generated from randomly chosen MD snapshots for states S1a (in purple sphere), S1b (in green), and S2 (in orange). (B) The three-state macrostate representation generated from the MSM construction. The equilibrium populations from the MSM are indicated. (C) Distance analyses between several key residues and the PPi group in the three-state representation.

https://doi.org/10.1371/journal.pcbi.1004624.g003

Next, for each macrostate, we analyzed the key residue interactions affecting the PPi release, by measuring the distance between each key residue and the PPi group. As shown in Fig 3C, three pairs of the distances, MgA-PPi, Asp812-MgB and Asp537-MgB, consistently increase from S1a to S1b and to S2, revealing the three-state character. In comparison, the distance Lys631-PPi increases and the corresponding association loosens only upon transition to S2, which demonstrate instead the two-state behavior. Interestingly, the distance Lys472-PPi gradually decreases so that the Lys472 association with PPi becomes comparatively stable in S2. In addition, the distance Arg627-PPi does not change obviously and the Arg627 association with PPi keeps moderate all along.

The microsecond simulations suggest that the PPi release is not tightly coupled to the O-helix opening but is assisted by Lys472 swing

As noted above, the MSM is constructed based on extensive short MD simulations for nanoseconds each, which are individually too short to account for continuous slow motions. In order to examine better the slow motions in the PPi release, we conducted microsecond long MD simulations to the following systems: (i) The unperturbed PPi-bound product complex of T7 RNAP; (ii) A modified product complex with the stabilizing charges to PPi inside (Asp537, Asp812, and Lys631) turned off (off-charge); and (iii) Another modified product complex with PPi removed (no PPi). All simulations started from some S1a configurations. Later in the simulation, (i) was kept well in S1a, (ii) moved slightly toward S1b (see S1A and S1B Fig). Although no PPi release event was yet captured, we observed bending and some opening motions of the O-helix. Essentially, we also identified significant side chain swings of Lys472.

The O-helix is supposed to rotate to open (~ 22°) in the post-translocated configuration [50], while in the PPi-bound product complex and prior to the RNAP translocation, the O-helix is in a closed configuration [7]. The closed to open transition of the O-helix was suggested to couple the PPi release to the polymerase translocation in the PS model [7]. Accordingly, in our simulation, we wanted to know whether and how the O-helix switches from closed to open in response to the PPi release. To monitor opening motions of the O-helix, we recorded the overall rotation angles of the O-helix (see Fig 4), as well as the rotation angles measured from the N- and C-terminus of the O-helix (see S3 Fig), respectively. In general, we found that the overall rotation of the O-helix is smaller than that measured from the two terminus, which indicates occasional bending of the O-helix. Our analyses show that in the unperturbed simulation (i), the O-helix fluctuates near the closed configuration (< 10°, see Fig 4A and 4B). Occasionally, the N-term opens much more (up >15°, see S3B Fig) than the C-term, so that the helix bends significantly. The O-helix bending, however, imposed limited effects on the PPi group, which kept stable in the active site. Interestingly, in the off-charge simulation (ii), we noticed substantial openings of the O-helix (up > 15°, see Fig 4A and 4C), in particular, from the N-terminus (up > 20°, see S3C Fig). Nevertheless, PPi was still kept inside without significant movements, even though the PPi-Lys631 and PPi-Asp537/Asp812 distances largely increase (see S1D Fig). Hence, the PPi group does not seem to respond to the increased O-helix opening within our microsecond simulation. Moreover, we noticed that some non-charged local residues, including Gly538, Cys540, and Ser541, also stabilize PPi through forming hydrogen bonds or coordination with the MgB ion (see S4 Fig). These local interactions were present to stabilize PPi inside around the active site, in both the unperturbed simulation (i) and the off-charge simulation (ii). Finally, in the simulation (iii) with PPi removed, the O-helix gradually opened in the beginning of the simulation (up ~ 12°). The opening, however, was not sustained, as the rotation angle reduced (after ~ 300 ns) to that similar to the unperturbed simulation in the presence of PPi (see Fig 4A and 4D). Therefore, the O-helix could not open in response to the removal of PPi.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. The O-helix rotation detected in microsecond MD simulations.

(A) The rotation angles of the O-helix were measured from simulation (i), (ii) and (iii), and are shown in blue, green, and red, respectively. (B) The molecular snapshot view of the O-helix (green) taken from the unperturbed simulation (i). The initial closed form of the O-helix is also shown (cyan). (C and D) The molecular snapshots taken from the off-charge simulation (ii) and the no PPi simulation (iii), respectively, with the O-helix shown as in (B).

https://doi.org/10.1371/journal.pcbi.1004624.g004

On the other hand, we found dramatic swing motions of Lys472 that could substantially aid the PPi release. Lys472 is located on a flexible loop (residue 469 to 475) on the protein surface, which ensures that the Lys472 side chain can rotate freely. In Fig 5A, we show two major configurations of the Lys472 side chain, sampled from both the short SMD simulations and the microsecond long simulations (i and iii). In one configuration, the Lys472 side chain contacts on PPi and points inside toward the active site; in the other configuration, the side chain swings outward and points toward solvent. Indeed, in the first ~500 ns of the unperturbed long simulation (i), the Lys472 side chain kept close to PPi (in the active site) but occasionally, it swung away from PPi to point outward, at a frequency about once per 100 ns (see Fig 5B). In the latter part of the unperturbed simulation, the Lys472 side chain settled down and pointed only to the active site, grabbing PPi all along. In the off-charge simulation (ii), the Lys472 grabbed on PPi most of time, with only a couple of times, it swung away from PPi (not shown). In the absence of PPi in the simulation (iii), however, the Lys472 side chain swung away and relaxed to point outward, in response to the removal of PPi.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. The side chain swing of Lys472 assists the PPi release.

(A) The molecular view of Lys472 side chain that swings from pointing inside toward PPi in the active site (S1) to pointing outside to the solvent (S2). The two major configurations of the side chain were sampled from the long simulation (i and iii) and from the SMD simulations. (B) The Lys472-MgA distance measured from the long simulation (i) and (iii), with and without PPi, respectively. (C) The 2D density map (-ln P) generated from many short simulations used for the MSM construction. The map is depicted along the distance between Lys472 (NZ) and MgA and along the PPi RMSD. (D) The Lys472-MgA distance change along with the PPi-MgA distance change, obtained from three SMD simulations for nanoseconds, ran from different initial conditions/directions to pull PPi out.

https://doi.org/10.1371/journal.pcbi.1004624.g005

Notably, the two configurations of the Ly472 side chain correspond well to the inside and outside states S1 (S1a and S1b) and S2 from the MSM, respectively. As shown in Fig 5C, we projected all the PPi conformations from the short MD simulations onto two reaction coordinates: One is the PPi RMSD, and the other one is the distance between the Lys472-NZ atom and MgA. It is clearly seen that the swing motion of the Lys472 side chain correlates well with the PPi release. That is, the Lys472 points towards the active site in S1 while it swings out to assist PPi to release toward S2. This high correlation was also observed in our SMD simulations. As shown in Fig 5D, in three PPi pulling events (with different initial conditions or pulling directions), we found that the Lys472 side chain swung along with PPi as PPi was pulled from the inside to the outside (see SMD S1 Movie in SI). Altogether with the above observations, we suggest that the side chain swing of Lys472 greatly promotes the PPi transition from S1 to S2, or say, the PPi release.

Conclusion

In this computational work, we examined structural dynamics of the PPi release process in the T7 RNAP transcription elongation at an atomistic resolution. We employed a large number of short MD simulations at nanoseconds to construct the MSM, combined with a small number of long simulations at microseconds to detect further slow essential motions. In contrast to the facilitated charge hopping mechanisms of the PPi release found in multi-subunit RNAPs, the PPi release in T7 RNAP is dominated by a jump-from-cavity activation process. Charge interactions from Lys631 on the O-helix and from two highly conserved aspartic acids Asp527 and 812 stabilize the PPi-MgB group inside around the active site to prevent the release. In addition, a group of non-charged residues interact locally with the PPi group to further hinder the activation process. Remarkably, a highly fluctuating lysine residue Lys472 grabs on PPi and pulls it out through significant side chain swings. The lysine and a similar arginine residue playing the same role can be found in other phage or mitochondria RNAPs, in structurally similar DNAPs, and even in the multi-subunit RNAPs. Hence, the flexible lysine/arginine side chain swing appears to be a common module to assist the PPi release, while the assisted jump-from-cavity mechanism of the PPi release may apply in general to a wide group of single subunit polymerases. Experimental mutations can be designed on those key residues to test their specific impacts. For example, turning off charge interactions individually from Lys631 and Lys472 is expected to accelerate and hinder the PPi release, respectively. Furthermore, the opening of the O-helix on the fingers domain does not appear to be tightly coupled to the PPi release according to our MD simulations, nor could the activated PPi release process energetically support the O-helix opening or further the translocation. Hence, current study favors the Brownian ratchet mechano-chemical coupling scenario over the power stroke one in the T7 RNAP elongation, consistent with single molecule studies revealing the Brownian ratchet nature of the RNAPs, no matter for single or multi-subunit species.

1. MD simulation setup for the PPi-bound T7 RNAP complex

The crystal structure of the PPi-bound complex of T7 RNAP (PDB: 1S77) [7] was solvated with TIP3P water in a cubic box, and the minimum distance from the protein to the wall was set to 10Å. The RESP charges of the PPi group were calculated in AMBER package [54], after using the RED V software via HF/6-31G* method for the structure optimization [55]. The additional structural parameters were derived from the general Amber force field (GAFF) [56,57] (see SI S1 and S2 Tables). To neutralize the system and keep the ionic concentration of 0.1 M, 116 Na⁺ ions and 68 Cl^- ions were added. There were totally ~113,000 atoms in the final system. The MD simulations were performed using the GROMACS-4.6.5 software package [58]. The AMBER99sb force field with PARMBSC0 nucleic acid parameters was used to describe the model [59,60]. The cut off value for van der Waals (vdW) and short-range electrostatic interactions was set to 10 Å. Long-range electrostatic interactions were treated with the Particle-Mesh Ewald (PME) summation method [61]. The time-step was 2 fs and the neighbor list was updated every five steps. Firstly, the system was minimized using the steepest decent method for 50,000 steps, then 100 ps MD simulation under canonical ensemble was carried out. The temperature was set at 310 K using the velocity-rescaling thermostat [62]. After that, 100 ps MD simulation was run under NTP ensemble at 1 bar and 310 K, using the Parrinello-Rahman barostat [63,64] and the velocity rescaling thermostat, respectively. In the MD simulation, position restraints were imposed on the heavy atoms at the beginning of the simulation (for ~ 1 ns), and were then removed in the later equilibration. The above procedures were conducted to prepare equilibrated initial structures for the SMD simulations, to perform short MD simulations for the MSM construction, and to perform the microsecond long MD simulations.

2. Generating initial PPi release pathways using the Steered MD (SMD) simulations

We implemented the SMD simulations [65] to generate the initial PPi release pathways, as that were conducted previously [25,26]. To consider all possible PPi release pathways, the PPi-MgB group was pulled out of the active site along five directions (see Supporting Information and S5 Fig for details). The external force was applied to the center of mass of the PPi group with a force constant 10 kJ/mol/Ų and the pulling rate was 0.01 Å/ps. During the SMD simulation, 100 kJ/mol/Å force restraints in the x, y, and z directions were imposed on protein Cα atoms and nucleic acid heavy atoms to maintain protein complex stability under the external force. For each pulling direction, three independent 3-ns SMD simulations were conducted (fifteen SMD simulations in total).

3. Seeding unbiased MD simulations to construct the Markov State Model (MSM)

The conformations obtained in the SMD simulations were grouped into 20 clusters using the K-center clustering algorithm [66]. In the clustering, the distance between a pair of conformations was set to the RMSD value of the PPi group (two phosphorus atoms and the bridged oxygen atom). The RMSD was computed by firstly aligning each conformation to one particular reference conformation where the PPi group locates in the active site right after the catalysis, according to the C_α atoms of the O-helix (627 to 640) on the finger domain. Then 3–5 conformations were randomly selected from each cluster (a total of 100 conformations) to conduct the following unbiased MD simulations. Each MD simulation was run for ~ 20 ns and the snapshots were saved every 2 ps. Altogether, an aggregation of ~2 μs simulations with 1000,000 conformations were obtained.

4. Construction the MSM for the PPi release

The Markov State Models (MSMs) describe molecular systems by partitioning the high-dimensional conformation space into a number of discrete metastable states [27–31], such that the conformational transitions within a certain metastable state are relatively fast comparing to the inter-state transitions. The disparate timescales allow constructing a markovian process where the probability of the transition from state i to state j (T_ij) after a certain lag time Δt is only dependent on the current state i and not any former states. The T_ij for any two states i and j can be estimated from extensive short MD simulations with each state well equilibrated. Then by building a transition probability matrix T, we can describe the stochastic transitions between different discrete states and propagate the markovian dynamics to a given long timescale: where P(nΔt) is a vector of state populations at time nΔt and T is the transition probability matrix. Here we adopted a two-step procedure to construct the MSM for the PPi release: 1) Clustering the MD conformations into hundreds of microstates based on their geometric differences. 2) Lumping the microstates into several macrostates in order to visualize the key intermediate states during the PPi release [66]. The detailed procedures are described below.

1) Clustering the MD conformations into various numbers of microstates. Based on the geometric differences calculated from the RMSD of the PPi group (see the previous section), we clustered all the 1,000,000 MD conformations into 200, 300, 400 or 500 microstates using the K-center clustering algorithm implemented in the MSMBuilder software [67,68]. For each cluster size, we then constructed MSM and calculated the implied timescales at different lag times using: where μ_k is the eigenvalue of the transition probability matrix T obtained in a certain lag time τ. Each implied time scale, representing an average transition time between two subsets of states, is one of the indicators to validate if the model is markovian [29]. For a sufficient long lag time τ, further increase of the lag time will not change the implied time scale when the model becomes markovian. As shown in S6A Fig, the implied timescales for different cluster sizes demonstrate similar kinetic behaviors for the slowest dynamics. Our result thus suggests that the MSM is robust to choices of different numbers of microstates, we therefore finally used one of the models, 200-state model, to further lump the microstates into macrostates and calculated the associated thermodynamic and kinetic properties. In the 200-state model, the average RMSD within each cluster was determined to be ~1.2 Å, indicating that the conformations within each cluster are geometrically similar to each other so that the kinetics within each cluster is relatively fast. The implied time scale curves level off after the lag time of 3 ns, indicating that the model becomes markovian. We therefore used the lag time of 3 ns as the markovian time to extract the thermodynamic and kinetic properties.

In order to further validate our MSM, we conducted the Chapman-Kolmogorov test [69]: we firstly selected the six most populated microstates, and for each state, we predicted its self-transition probability based on the transition probability matrix T generated at the markovian time of 3 ns. A good agreement is reached between the predicted values and the ones that are directly obtained from the MD simulations (see S6B Fig).

2) Visualizing key intermediates by lumping microstates into macrostates. We coarse-grained the microstates into a few macrostates in order to visualize the key intermediate states, which in turn, can help to elucidate specific mechanisms of the PPi release. Based on the 200-state model, we lumped the 200 microstates into three macrostates using the PCCA+ algorithm [70] implemented in the MSMBuilder. The equilibrium populations of the three macrostates and the MFPT of each macrostate-level transition were then calculated.

3) Calculating the Mean First Passage Time (MFPT). The MFPT for each pair of the states can be calculated using the following formula: Where f_ij is the MFPT for the transition from state i to state j, P_ik is the transition probability from state i to state k under a certain lag time t_ik (here t_ik is the markovian time and equals to 3 ns). Under the boundary condition of f_jj = 0, the f_ij can be determined by solving a set of linear equations. Based on the 200-state transition probability matrix constructed at the lag time of 3ns, we generated ten 10 ms long Monte Carlo (MC) trajectories. Then from each MC trajectory, we collected all the MFPT values for both the S1a → S1b and S1b → S2 transitions and plotted the corresponding occurrence counts (see S7 Fig for the results of one of the MC trajectories). Next, we truncated those very small and very large MFPT values that correspond to the fastest and slowest dynamics, respectively. In such a way, we obtained a high-quality linear regression (see S7 Fig) according to the Poisson distribution. Upon that, we finally calculated the average of the MFPT for each macrostate transition. The average MFPTs were then obtained by averaging over the ten MC trajectories and the corresponding standard errors of the averages were calculated.

We also estimated the upper bound of the macrostate MFPT by obtaining the weighted sum of the microstate MFPTs without truncation. When the Markovian property is well maintained at the macrostate level, the estimation should provide a similar value as the average MFPT calculated above. However, in current case, as indicated above, the macrostate Markovian property is not well maintained without truncation (see S7A and S7B Fig). Consequently, the estimation here is dominated by the large microstate MFPTs, or say, by slow microstate transitions. In this way, the estimation gives an upper bound value of the MFPT for S1a → S1b at ~ 9 μs, and that for S1b → S2 at ~ 20 μs.

Besides, partitioning the continuous phase space into a number of discrete states can result in losing the dynamic information within each discrete state. Therefore, using the Markov process to describe the stochastic motions of the biological system is an approximation to the continuous dynamics. We thus should take the errors caused by the discretization process into account when building the MSM. In order to evaluate the effects of the discretization errors on the MFPT, we calculated the MFPT values for other three models with different discretization resolutions (300, 400 and 500 microstates, respectively) by lumping the corresponding microstates into three macrostates. From S8 Fig, we can see that the MFPT values have very limited changes for different models, varying from 0.16 μs to 0.26 μs for the S1a → S1b transition and from 0.94 μs to 1.24 μs for S1b → S2 transition, suggesting that the cluster size does not impact on the MFPT values significantly. It is also noted that although increasing number of microstates leads to smaller discretization errors, it also increases the statistic errors due to fewer transition counts between different states. We thus regard the 200-state model a reasonable one to build the MSM and extract the thermodynamic and kinetic properties, since the conformational space of each microstate is small enough to ensure a fast transition kinetics (the average RMSD comparing to the center conformation within each cluster was ~1.2Å).

5. Identifying slow modes and implied time scales using the trajectory mapping (TM) method

In order to confirm that the MSM construction above is robust using distance metrics other than the PPi RMSD, we utilized the trajectory mapping (TM) method [51] to find slow modes in the PPi release process. A total of 45 distances between oxygen or phosphorous atoms of PPi and Cα atoms of residues near the PPi release channel were used to identify configurations. After mapping the simulated trajectories, these distances are linearly combined to the principle components of the trajectory mapped vectors, which approximately correspond to slow variables in the PPi release process. As shown in S9 Fig, the first principle component (corresponding to the slowest variable) is closely correlated to the distance MgA-PPi (or Lys472-MgA), indicating the consistence between the MSM using the PPi RMSD and the TM analysis using the general distances. The implied time scale plots according to the first and first five principle components are also provided, showing similar results as that from the MSM.