Pyrophosphate ion (PPi) release during transcription elongation is a signature step in each nucleotide addition cycle. The kinetics and energetics of the process as well as how it proceeds with substantial conformational changes of the polymerase complex determine the mechano-chemical coupling mechanism of the transcription elongation. Here we investigated detailed dynamics of the PPi release process in a single-subunit RNA polymerase (RNAP) from bacteriophage T7, implementing all-atom molecular dynamics (MD) simulations. We obtained a jump-from-cavity kinetic model of the PPi release utilizing extensive nanosecond MD simulations. We found that the PPi release in T7 RNAP is initiated by the PPi dissociation from two catalytic aspartic acids, followed by a comparatively slow jump-from-cavity activation process. Combining with a number of microsecond long MD simulations, we also found that the activation process is hindered by charged residue associations as well as by local steric and hydrogen bond interactions. On the other hand, the activation is greatly assisted by a highly flexible lysine residue Lys472 that swings its side chain to pull PPi out. The mechanism can apply in general to single subunit RNA and DNA polymerases with similar molecular structures and conserved key residues. Remarkably, the flexible lysine or arginine residue appears to be a universal module that assists the PPi release even in multi-subunit RNAPs with charge facilitated hopping mechanisms. We also noticed that the PPi release is not tightly coupled to opening motions of an O-helix on the fingers domain of T7 RNAP according to the microsecond MD simulations. Our study thus supports the Brownian ratchet scenario of the mechano-chemical coupling in the transcription elongation of the single-subunit polymerase.
Transcription elongation is a continuous process of producing messenger RNA as an RNA polymerase (RNAP) moves along double stranded (ds) DNA, copying information from the DNA template to synthesize an RNA strand [1–3]. For each nucleotide addition cycle (NAC) of the transcription elongation, a nucleoside triphosphate (NTP) is recruited into the active site of RNAP. Then phosphoryl transfer reaction happens which leads to the addition of the nucleoside monophosphate (NMP) to the existing RNA strand, followed by the release of a pyrophosphate ion (PPi). The PPi release hence serves as a nice signal for detecting each NAC, for example, in real-time DNA sequencing [4,5]. Although the PPi release in RNAP can be detected, how it happens in structural dynamics detail remains elusive from existing experimental measurements. Since the polymerase works as a molecular motor driven by chemical free energy from each NAC, how the PPi release after catalysis couples with essential conformation changes and translocation of the polymerase determines the mechano-chemical coupling nature of the motor system [1,6–9]. In studying the RNAP structures from bacteriophage T7, it was suggested that the PPi release directly drives the polymerase translocation [7]. Such a mechanism is called the ‘power stroke’ (PS), as the product release reaction couples tightly to the mechanical movement of the protein [10,11]. On the other hand, single molecule measurements implementing optical-tweezer forces on RNAPs, for single or multi-subunit ones, consistently suggested a ‘Brownian ratchet’ (BR) mechanism [12–14], in which the translocation proceeds in Brownian motions, thermally activated without being energetically coupled to the preceding PPi release. As such, it is necessary to clarify whether T7 RNAP makes an exception to the BR type RNAP motors. Besides, studies that are able to zoom into structural detail and to characterize kinetic or energetics of the process are highly demanded for understanding the basics of the transcription engine. T7 RNAP is a prototypical single subunit polymerase that works self-sufficiently without transcription factors. It maintains high promoter specificity and drives transcription strongly [15]. T7 RNAP has been widely utilized in protein over-expression systems as well as in synthetic gene circuits [16,17]. It is also characterized well from structural and kinetic studies at single molecule level [7,14,18–22] as well as from biochemical ensemble measurements [6,23,24]. The PS mechanism was suggested based on the high-resolution structural studies [7]. It is observed that a flexible fingers domain with five alpha helices including an O-helix adopts two major conformations, the closed and open forms in each NAC. In the PS mechanism, a transition from the closed to open form was supposed to couple the PPi release to the polymerase translocation with a significant free energy decrease. Nevertheless, in the single molecule force measurements of T7 RNAP, the elongation rates measured at variable NTP concentrations and forces were fitted to a kinetic scheme, yielding only a very small free energy bias upon the translocation (~ 1kBT) [14,20], supporting the BR rather than the PS scenario. Furthermore, transient state kinetics measurements on the T7 RNAP transcription have also been conducted [23]. The study suggested that an isomerization from the open to the closed form prior to the catalysis is rate limiting, while the steps afterwards including the PPi release and translocation are fast. The study accordingly provides a kinetic framework to investigate the structure-function detail of the transcription elongation in T7 RNAP. The dissociation constant of PPi was determined at ~ 1. 2 mM [23]. Recently, computational work combining extensive all-atom molecular dynamics (MD) simulations with the Markov state model (MSM) has been conducted to study the PPi release in multi-subunit RNAPs from bacteria to yeast [25,26]. The advantage of implementing the MSM is to efficiently group high-dimensional data and extract kinetic information from a large number of short MD simulations sampled widely in the conformational space [27–31]. One big challenge the conventional MD simulations face with is to reach physiologically relevant timescales, especially for the system as large as the RNAPs (system size of a hundred thousand to a million atoms, solvated in explicit water). The MSM had been implemented successfully to describe relatively long dynamics in protein [32,33] and RNA folding [34], receptor-ligand binding [35–37], etc. . The MSM constructed for the multi-subunit RNAPs suggested that the PPi release happens at around a microsecond [25,26], faster than the polymerase translocation (estimated to be tens of microseconds at least) [38]. In the multi-subunit RNAPs, there is an important structure called trigger loop (TL) [39,40] that folds (closes) and unfolds (opens) upon the NTP binding and PPi release, respectively, analogous to the O-helix closing and opening in the single-subunit T7 RNAP. In the yeast RNAPII, it was noticed that the TL increases its fluctuation after catalysis to aid the PPi release, which occurs through a hopping mode whereby four meta-stable states were identified [26]. In comparison, a relatively simple two-state model of the PPi release was discovered in bacterial RNAP [25]. The TL only slightly opens in the yeast RNAP II upon the PPi release, while the TL does not open in a same way in the bacterial RNAP [25,26]. Computational studies on single-subunit polymerases have also been conducted recently [41–47]. The structurally similar single-subunit polymerases include not only RNAPs from bacteriophage and human mitochondria, but also a group of DNA polymerases (DNAPs). The MD study on the translocation mechanism of DNAP I [44] showed that the PPi release facilitates the opening transition of the fingers domain. The authors therefore suggested that the PPi release triggers the translocation by facilitating the domain opening transition of the polymerase [44]. More recently, structured-based kinetic modeling and all-atom MD simulations on T7 RNAP have been conducted [42,43]. The kinetic modeling was built upon single molecule force measurements on the T7 RNAP transcription elongation[14,20], supporting the BR mechanism [43]. The ensued MD study focused on nucleotide selection prior to the catalysis for transcription fidelity control [42]. Further computational studies on the post-catalysis coordination from the PPi release to the translocation are highly desired, therefore, for a complete description of the elongation cycle and for resolving the conflicting views on the mechano-chemical coupling of the transcription machine. In the single-subunit T7 RNAP, the active site is buried inside the polymerase within a distance of ~15Å away from the bulk water. Right after the catalysis, the negatively charged PPi sits inside a cavity, surrounded by three positively charged residues (Lys631, Arg627, and Lys472) from the outside, as well as by two negatively charged residue (Asp537 and Asp812; through a magnesium ion bound with PPi) from the inside (see Fig 1). Accordingly, the releasing pathway of the PPi group from the active site is blocked by the side chains of the above residues. The structural features of the PPi release channel appear quite different from that of the multi-subunit RNAPs, which have a fairly long pore region with a length of ~30 Å. A comparative study of the PPi release from the viral T7 RNAP with that from the bacteria and high organisms would then also help to elucidate how variable structural designs accommodate a convergent bio-molecular function, achieved through enzyme evolution from different species. In this work, we started from an atomistic structure of the T7 RNAP elongation complex in the PPi-bound product state and conducted extensive MD simulations. The MSM was then constructed to extract kinetic information from the grouped snapshots sampled along the PPi release pathways. Comparing to the four-state and two-state PPi hopping model proposed previously for the multi-subunit RNAPs, we investigated the PPi release mechanism in T7 RNAP and built a corresponding model, concerning (1) how many metabstable states exist during the PPi release, (2) is the energetics going up or down along the release path, and (3) what is the rate-limiting transition in the full process. We obtained a jump-from-cavity kinetic model of the PPi release from the extensive simulations, with two relatively stable states inside the cavity, and one marginally stable state at the exit. We also found that the PPi jump from the cavity is rate limiting over a microsecond. In addition, we performed microsecond long MD simulations in both unperturbed and modified conditions. With that we identified which slow motions essentially couple with or impact on the PPi release. Overall, we demonstrated that the PPi release in T7 RNAP is a jump-from-cavity activation process closely assisted by a flexible lysine residue, without being tightly coupled to the opening of the O-helix to further support the translocation. In the end, we compare T7 RNAP with other RNA or DNA polymerases on conserved structural features that are key to 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 Mg2+ 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 Mg2+ 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). 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). 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. 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. 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. 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.
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.
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. 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/Å2 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). 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. 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 (Tij) after a certain lag time Δt is only dependent on the current state i and not any former states. The Tij 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: P (nΔt) =[T (Δt) ]nP (0) 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: τk=−τlnμk (τ) 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: fij=∑kPik× (tik+fkj) Where fij is the MFPT for the transition from state i to state j, Pik is the transition probability from state i to state k under a certain lag time tik (here tik is the markovian time and equals to 3 ns). Under the boundary condition of fjj = 0, the fij 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Å). 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.
RNA polymerases (RNAPs) are usually recognized to work as Brownian ratchet machines, with mechanical movements and chemical reactions loosely coupled during transcription elongation. Nevertheless, structural studies on the single-subunit T7 RNAP had suggested an alternative power stroke scenario, which requires the PPi product release reaction tightly couples with or drives the polymerase translocation. To resolve the conflicting views in the mechano-chemical coupling scenario, we conducted atomistic molecular dynamics (MD) simulations to investigate the PPi release mechanism of T7 RNAP. Using a large number of nanosecond short simulations, we constructed the Markov state model (MSM) and found that the PPi release undergoes a jump-from-cavity activation process, unlikely to further support the translocation. The MSM implementations to multi-subunit RNAPs previously demonstrated instead the charge facilitated hopping mechanisms of the PPi release. To further explore essential slow motions during the PPi release, we also performed microsecond long MD simulations. The PPi release does not appear to be tightly coupled to the rotational opening of an O-helix that is directly tied to the polymerase translocation. Hence, the study again disfavors the power stroke scenario. Remarkably, we discovered a key residue Lys472 that is able to assist PPi to jump out of the cavity. The mechanism can be general to a group of polymerases sharing conserved structural features with T7 RNAP, while the lysine or arginine module to assist the PPi release can be universal to both the single and multi-subunit RNAPs, even though the overall structural features and the PPi release mechanisms are significantly different.
A Jump-from-Cavity Pyrophosphate Ion Release Assisted by a Key Lysine Residue in T7 RNA Polymerase Transcription Elongation