|Home | About | Journals | Submit | Contact Us | Français|
We present results from molecular dynamics simulations and free energy calculations of the twister ribozyme at different stages along the reaction path in order to gain insight into its mechanism. Results, together with recent biochemical experiments, provide support for a mechanism involving general acid catalysis by a conserved adenine residue in the active site. Although adenine has been previously implicated as a general acid acting through the N1 position in other ribozymes such as the hairpin and VS ribozymes, in the twister ribozyme there may be - a twist. Biochemical experiments suggest that general acid catalysis may occur through the N3 position, which has never before been implicated in this role; however, there currently lacks a detailed structural model for the active state of the twister ribozyme in solution that is consistent with these and other experiments. Simulations in a crystalline environment reported here are consistent with X-ray crystallographic data, and suggest that crystal packing contacts trap the RNA in an inactive conformation with U-1 in an extruded state that is incompatible with an in-line attack to the scissile phosphate. Simulations in solution, on the other hand, reveal this region to be dynamic and able to adopt a conformation where U-1 is stacked with G33. In this state, the nucleophile is in line with the scissile phosphate, and the N1 position of G33 and N3 position of A1 are poised to act as general base and acid, respectively, as supported by mutational experiments. Free energy calculations further predict the electrostatic environment causes a shift of the microscopic pKa at the N3 position of A1 toward neutrality by approximately 5 pKa units. These results offer a unified interpretation of a broad range of currently available experimental data that points to a novel mode of general acid catalysis through the N3 position of an adenine nucleobase, thus expanding the repertoire of known mechanistic strategies employed by small nucleolytic ribozymes.
The twister ribozyme is a recently discovered member of the group of small nucleolytic ribozymes that catalyze the site-specific cleavage of the RNA phosphodiester backbone. This group of self-cleaving ribozymes serves as an important testbed in the study of the mechanisms of RNA catalysis,1–6 and has impacted in the design of new biomedical technology7–12 and theories of evolution.13,14
The twister ribozyme was originally identified from bioinformatics, and subsequently has been studied extensively by X-ray crystallography,15–17 mutagenesis,15,16,18 and other biochemical experiments.15,16 The twister ribozyme is of special interest in that experiments suggest that metal ions, although important for folding, do not play a direct role in catalysis.18 This implies that catalysis derives mainly from the electrostatic engineering of the active site and participation of nucleobase functional groups, which are generally considered as fairly chemically inert when compared to amino acids.3,19 In the case of the twister ribozyme, mutational studies have implicated that the N3 position of a conserved adenine residue, A1, is critical for activity, and may serve, in protonated form, as a general acid in the reaction. If this supposition is true, it would represent an intriguing new twist on general acid catalysis, and therefore be important in furthering our understanding of the mechanistic strategies employed by RNA enzymes.
Crystallographic data15–17 are suggestive in that the N3 position of A1 is in reasonably close proximity to the O5′ leaving group to potentially act as a general acid. However, none of the currently available crystal structures15–17 are likely representative of the active state of the twister ribozyme, obfuscating their functional interpretation. Thus, there remain significant gaps in our understanding of the twister ribozyme based on currently available experimental data that prevent definitive conclusions to be drawn about its mechanism. Specifically, there currently does not exist a structural (or dynamical) model of the active state of the twister ribozyme in solution that is consistent with the current body of available experimental data. As a result, there remain several key questions as to: 1) what is the origin of the inactive state observed crystallographically15,17? 2) what conformational events are required to form an active in-line conformation in solution, and what is the probability of finding the ribozyme in this conformation? 3) what residues likely act as the general base and acid? and 4) how do electrostatics shift the microscopic pKa values of implicated general acid and base residues, and help to stabilize the transition state?
In this paper, we utilize molecular dynamics simulations in a crystalline environment and in solution, along with free energy calculations, to gain insight into these key questions. Our results provide a unified molecular-level interpretation of a broad range of experimental data15–18 that sheds new light onto the mechanism of the twister ribozyme.
Crystallographic data provides valuable information about ribozyme structure and function; however, it remains an open question as to the degree to which a static picture of a deactivated ribozyme in a crystalline environment may have direct catalytic relevance for the active ribozyme in solution.20 In the case of the twister ribozyme, with each of the three available crystal structures (PDB ID: 40JI,15 4RGE,16 and 4QJD/4QJH17), U-1 both lacks the O2′ nucleophile and is seen to be extruded from the active site. As a consequence, these structures do not reflect a conformation where the nucleophile is both in line (so as to be able to attack the scissile phosphate) and directly interacting, through hydrogen bonding, with a residue with the potential to act as the general base. These features are generally considered to be a requirement for catalysis.21 This observation then begs the question: is the structure in the crystal catalytically relevant? In order to address these questions, we performed molecular dynamics simulations of the twister ribozyme, departing from the high resolution structure15 (PDBID: 4OJI), both in a crystalline environment and in solution (see Methods). We first set out to ascertain the consistency of our crystal simulations with available crystallographic data, and second, to characterize the dynamics in the region of the active site and, in particular, the local geometry associated with the orientation of U-1.
Overall, the crystal simulations were in very close agreement with available experimental data. The average structure from the crystal simulation was almost identical to the crystal structure (Figure 1), particularly in the active site where the heavy-atom root-mean-square deviation (RMSD) was 0.44 Å. Furthermore, the B-values determined from the simulated mean-square positional fluctuations (Figure 2), represent the degree of structural variation, and agree closely with those determined from crystallographic data (linear correlation coefficient of 0.90). Additional details on the observed structural fluctuations within the crystal simulation can be found in Figure S2. Taken together, this provides credence for the reliability of the simulations, and also underscores the supposition that the region of residues 17 and 18 are highly disordered as interpreted by the lack of electron density in the 2Fo-Fc map.15 The crystal simulations predict that U-1 is stabilized by a crystal packing contact with G14 in a symmetry related monomer, where it remains trapped for the duration of the simulation and does not sample configurations that allow the 2′ position to come in line with the scissile phosphate (Figure S3). Together with the relatively modest experimental and predicted fluctuations in this region, this suggests U-1 does not sample catalytically relevant conformations in the crystal.
Over the course of a 120 ns unrestrained simulation of the twister ribozyme in solution, we identified three distinct conformational states of U-1 (Figure S4) that we denote as "Extruded", "Stacked" and "Triple". In the "Extruded" state, U-1 is completely solvent exposed in an orientation similar to that of the crystal (but with greater dynamical fluctuations since it is not stabilized by crystal packing contacts). In the "Stacked" state, U-1 stacks under G33 while forming a hydrogen bond with a non-bridging oxygen of A34, and in the "Triple" state the additional hydrogen bonds characteristic of a WC/H/WC (U-1/A34/A19) base triple are observed. We found these different states could be distinguished by a single coordinate, the distance between U-1:N3 - A34:N7, henceforth denoted as Dstack: "Extruded" (Dstack > 7.0 Å), "Stacked" (3.75 Å< Dstack < 7.0 Å) and "Triple" (Dstack < 3.75 Å). The average structures from a set of three 75 ns simulations restrained to each of the aforementioned states are shown in Figure 3 (with a similar set of structures for G33 deprotonated shown in Figure S5) Additionally, this coordinate can be used to assess the occupancy of each state in solution, and estimate the kinetics of transitions between these states. Figure 4 shows the free energy profile along the stacking coordinate for simulations with G33, which has been implicated as a general base, in its neutral form and deprotonated at the N1 position as would be required for general base catalysis.
In both cases, simulations predict that U-1 rotation from the extruded state to the stacked state (the global free energy minimum) is essentially unhindered in solution. While the triple state is higher in free energy than the stacked state by less than 2.5 kcal/mol with G33 in its neutral form, it is considerably higher (greater than 8 kcal/mol) for G33 in the deprotonated form. Further, there is a considerable forward (~3.5 kcal/mol) and reverse (6-11.5 kcal/mol) barriers to transitions between the triple and stacked states. Overall, the stacked state is predicted to be long-lived and the dominant state in solution. As will be discussed below, this has important implications for catalysis.
To explore the impact of the different stacking states on the active site structure where the nucleophile adopts a catalytically active in-line conformation, we performed sets of simulations restrained to each of these states, considering G33 both in the neutral (protonated at N1) and anionic (deprotonated at N1) forms (see Methods).
Figure 5 shows a profile for the probability of finding the nucleophile in an active in-line attack conformation for each of the stacking states with G33 in its ground state neutral form. Whereas active in-line attack conformations are not observed in the extruded state, and are partially observed (51.0% occupancy) in the triple state, they are dominantly clustered in the stacked state. This clustering is preserved (96.4% occupancy) in simulations with G33 activated (deprotonated) for base catalysis (Figure 5: G33N1− Stacked), and further analysis reveals this state is marked by preservation of key hydrogen bonds between the endocyclic and exocyclic nitrogens of G33 with the nucleophile 2′OH and non-bridging oxygen of the scissile phosphate, respectively (Figure 5). This provides strong computational support that the twister ribozyme forms in-line attack conformations while in the dominant stacked state in solution, and upon deprotonation of G33, forms a hydrogen bond network that facilitates nucleophile activation with the N1 position of G33 acting as a general base catalyst.
With G33 poised to act as a base (G33N1−), the sampling of in-line conformations significantly increases for both the "Triple" and "Extruded" states, with occupancies of 93.1% and 57.5%, respectively (Figure S6). However, the sampling of in-line conformations is only one of the bases for developing a model of the twister ribozyme’s mechanism. With this in mind, the "Triple" state is disfavored as the most relevant active state, firstly, because of the higher free energy associated with sampling that state relative to the "Stacked" or "Extruded" states (Figure 4) and secondly, because of the loss of the hydrogen bond between G33N2 and the non-briding phosphoryl oxygen of A1, which aids in neutralizing the negative charge on the phosphate (Figure 3c and Figure S5c). In contrast, with G33 deprotonated, the "Extruded" state does exhibit the important hydrogen bonding interactions between G33 and A1. However, the "Stacked" state is both lower in free energy and has a significantly higher probability of the nucleophile being positioned in line, when compared to the "Extruded" state.
In addition to the six restrained simulations exploring the impact of U-1 stacking on the sampling of in-line conformations and the potential role of G33 as both a general base and hydrogen bond donor to the scissile phosphate, one additional restrained simulation was run to examine the interaction between the proposed general acid (A1:N3) and the 5′-oxygen leaving group. As with the other simulations of U-1 in the "Stacked" state, the probability of finding the ribozyme in an in-line conformation (Figure 6a) is extremely high (99.1% occupancy). Additionally, the hydrogen bonding between G33 and either the nucleophile and the pro-R non-briding oxygen are observed in a majority of simulation snapshots. As seen in Figure 6b, the distribution of distances between A1:N3 and A1:O5′ is bimodal, indicating that there are two distinct states being sampled. The first state (Figure 6c), is strikingly similar to the G33N1− "Stacked" state (Figure S5b) where the general acid is neutral, and there is no direct interaction between the presumed acid and the leaving group. However, a small fluctuation in the position of A1 (Figure 6d) allows for N3 to donate a hydrogen bond to the 5′-oxygen, enhancing it as a leaving group. This data is consistent with the proposed active state model, where U-1 stacks beneath the general base G33 and A1:N3 is poised to act as the general acid.
Simulations of a transition state mimic (TSM) allow the prediction of the active site architecture and hydrogen bonding network at a critical point along the reaction pathway. Simulations of the TSM were performed with the twister ribozyme in a protonation state that would support general acid/base catalysis in accordance with implicated residues identified through mutagenesis and biochemical studies.15 Specifically, both the N3 position of A1 (the presumed general acid) and the N1 position of G33 (the presumed general base) were protonated, representing a transition state that follows from the active state (designated as AH+EB−) where G33 was negatively charged and has accepted a proton from the 2′-oxygen nucleophile. Under these assumptions, examination of the TSM simulations serves as a predictive test; if the integrity of the active site and catalytic requirements of the hydrogen bond network is preserved, this provides support for the roles of A1 and G33 as general acid and base, respectively. On the other hand, if the TSM simulations are not stable, this creates strong evidence against the working hypothesis.
Overall, simulations of the TSM in the stacked state preserve the key hydrogen bonds necessary for catalysis, while providing important electrostatic stabilization of the negatively charged transition state. Figure 7 shows that G33 maintains hydrogen bonding with the pro-R oxygen of the scissile phosphate, while at the same time, A1 donates a hydrogen bond via the protonated N3 position to the O5′ leaving group. Results of the TSM simulations in the triple state (Figure S7b) are generally similar with those of the stacked state. Whereas, comparison of the TSM simulations in the extruded state (Figure S7a) are markedly different. In particular, there is a notable rotation of the backbone alpha and gamma torsion angles for the scissile phosphate. The alpha/gamma torsion angles flip from g−/g− in the "Extruded" state to t/g+ and g+/t for the "Triple" and "Stacked" states, respectively. The rearrangement of the backbone to t/g+ or g+/t allows for the 5′ leaving group to be positioned in close proximity to A1. These results suggest U-1 must stack underneath G33 in order to preserve the necessary hydrogen bond network to promote general acid/base catalysis by A1/G33, consistent with the "Stacked" state as the preferred working model of catalysis for the twister ribozyme.
Although the present simulations, along with crystallographic data,15 indicate that A1:N3 is positioned to donate a hydrogen bond to the O5′ leaving group, the microscopic pKa of the N3 site for adenosine in solution has been predicted to be quite low (1.5 ±0.3),22 raising possible objections to its effectiveness as an acid catalyst near neutral pH.17 This poses a question as to whether the electrostatic environment in the active site might cause the pKa of A1:N3 to shift toward neutrality, and if so, to what extent. The observation of such a large shift in the pKa of a nucleic acid fuctional group would not be unprecedented in either experimental23,24 or computational studies.25 Therefore, we have performed free energy calculations to predict the apparent pKa values for the proposed general acid (A1:N3) and base (G33:N1) in order to address this key question.
In keeping with the microscopic kinetic model for general-acid base catalysis26 (Figure 8), the free energy (and thus the pKa) for each leg of the cycle was calculated using a series of thermodynamic integration calculations (see Methods). From these free energy calculations, predicted fractions of each microstate could be determined from the partition function Q,
with the fraction of the active state, AH+EB−, determined as
The observed rate, which is pH-dependent, is presumed to be directly proportional to the active fraction according to the relation:
where kcl is a proportionality constant, representing the pH independent rate of cleavage. Furthermore, the simulated active fraction (AH+EB−) can be fit to a two state apparent pKa model27 that assumes the pKa values of the general acid and base are uncorrelated. This model has the three free parameters kcl, pKa,A, and pKa,B, with the latter two parameters representing the apparent pKa of the acid and base participating in catalysis, respectively:
By fitting the function as described in Eq. (4) to the simulated active fraction (i.e., in the same manner that experimental activity-pH curves are analyzed and interpreted), one obtains a predicted pH-activity curve and apparent pKa values that can be compared with experiment (Figure 9). The predicted apparent pKa values of 6.5±0.4 and 9.0±0.4 for A1:N3 and G33:N1 acting as the general acid and base respectively, are very close to the experimentally observed values of 6.9 and 9.5, respectively, and the overall shape of the simulated and experimental pH-activity profiles are very similar.15 These calculations support the hypothesis that the local environment within the active site of the twister ribozyme significantly shifts the pKa of A1:N3 (by ≈5 pKa units) towards neutrality, to provide a new mode of general acid-base catalysis.
The simulations presented here provide a uniform interpretation of a large body of experimental data in terms of a molecular-level description of catalytic mechanism for the twister ribozyme. The results depict a mechanism whereby U-1, which is trapped in an extruded state in the crystal due to crystal packing contacts, adopts a stable stacked state with G33 in solution. In this stacked state, the active site is able both to adopt an active in-line attack conformation of the nucleophile with the scissile phosphate and to position nucleotide residues to act as general acid and base catalysts. The transition state mimic simulations support a model whereby G33:N1 maintains hydrogen bonding with the O2′ position (which has a partial covalent bond to phosphorus) and A1:N3 donates a hydrogen bond to the O5′ position of the leaving group. Finally, free energy calculations suggest that the active site provides electrostatic stabilization of the transition state, and shifts the pKa of A1:N3 toward neutrality by almost 5 pKa units. Together, these results support the hypothesis that the twister ribozyme employs a new mode of general acid-base catalysis involving the N3 position of a protonated adenine, adding to the growing body of knowledge about the diversity of strategies employed by catalytic RNAs.
In keeping with the methodology set forth in the literature,28–31 we completed a 120 ns simulation of the ribozyme in the crystalline environment in an effort to both validate our computational models and provide the foundation for investigating twister in solution. The simulated crystal was generated by applying the appropriate symmetry operations to the A conformation of the X-ray structure deposited as PDB ID: 4OJI,15 such that the unit cell contained 12 monomers arranged in the space group P 65 2 2. The system was then solvated with TIP4P-Ew32 water, and equilibrated over 20 ns at constant temperature and pressure conditions to ensure that the system volume was within 0.3% of the experimental value. In a continued effort to replicate the experimental conditions, ammonium acetate and magnesium chloride were added to the system to match the experimental bulk concentrations of 0.1M and 0.02M, respectively. Following approximately 20 ns of equilibration, the simulation was carried out using the GPU accelerated version of PMEMD1433,34 at 280K and 1 atm using the AMBER FF14SB force field35 with ion parameters36 for use with the TIP4P-Ew water model. The average monomer structure from this simulation was generated by first reversing the symmetry operations for each of the 12 monomers in the unit cell in each 10 ps snapshot, then averaging over the coordinates of all monomers as well as over all snapshots in the last 100 ns of the 120 ns trajectory. Analysis of the data was conducted using both the XtalAnalyze30 module in AMBER14 and in-house scripts.
For each of the solution phase simulations, the ribozyme was solvated in a truncated octahedral box of TIP4P-Ew water with 12 Å between the solute and box, additionally including the crystallographically observed Mg2+ as well as 140 mM NaCl. All of the following solution phase simulations were carried out at 300K and 1 atm, using the AMBER FF14SB force field. The simulated annealing and solute equilibration procedure applied here has been used for other ribozymes and is detailed elsewhere.37 Details on restraints applied and data collection specific to the various simulations completed is discussed below.
Two free energy profiles were generated via umbrella sampling, by first minimizing representative structures of the "Stacked" state (Figure 3b or Figure S5b) while restraining the U-1:N3 - A34:N7 distance to 5.5 Å. This distance restraint was applied using a harmonic potential with a force constant of 50 kcal/mol-Å2. A short equilibration step of 1 ns was performed under constant temperature and pressure conditions, while continuing to enforce the above restraint. The structures used for each of the umbrella sampling windows were spaced every 0.5 Å from 3.0 to 11.0 Å and were equilibrated in serial from the previously equilibrated, adjacent window (starting from the initial window at 5.5 Å). The distance distribution data for each of the 17 windows, used for generating the free energy profile, was gathered every picosecond from the last 10 ns of 15 ns trajectories. For the data collection stage of the simulations, the restraint weight was reduced to 10 kcal/mol-Å2 to provide adequate overlap between adjacent sampling windows. One additional window at 3.75 Å was needed for the G33N1− profile to improve sampling at the top of the barrier between the "Stacked" and "Triple" states. The free energy profiles were built using vFEP38 and are presented in Figure 4. In order to estimate the statistical error of these calculations, 100 random subsamples of 10% of the data were collected. The 95% confidence intervals calculated from this distribution of subsamples are essentially the width of the lines in the plots presented in Figure 4.
During the equilibration stage (20 ns) restraints enforcing both the stacking state of U-1 and an "active" in-line conformation were maintained. The additional restraints supporting the "active" in-line conformation were slowly relaxed in the same fashion as described in the literature,37 and then removed for the data production stage. Keeping the restraints on Dstack ensured that throughout the course of the simulations statistics could be gathered on a single stacking state of U-1, since the barriers are sufficiently small to allow for transitions between the states. The nucleophile attack angle (θinl) defined as the angle between U-1:O2′, A1:P, and A1:O5′, and distance between the nucleophile and the phosphate (Dinl).21 This data was collected every 10 ps over the course of each of seven 75 ns restrained simulations, and clustered based on the hydrogen bonding patterning involving G33, the 2′-hydroxyl nucleophile, and the pro-R phosphate oxygen as seen in Figure 5 and Figure 6a.
Initially, restraints enforcing the characteristic hydrogen bonding observed for each of the U-1 stacking states were applied. Once minimized, these structures were equilibrated, removing the restraints, using the same procedure as referenced previously and then simulated at constant temperature and pressure for 120 ns. Data was collected from the last 100 ns of the simulations.
In order to determine the pH dependent probabilities of each of the four microstates in Figure 8, a series of thermodynamic integration calculations were completed. First, as a reference, the free energy of deprotonating the N3 position of adenosine and the N1 position of guanosine was calculated and assigned pKa values of 1.5±0.322 and 9.4±0.2,39 respectively. The model nucleotides were simulated as methyl capped trinucleotides, including the residues flanking the proposed general acid and base in twister (e.g. A1 was modeled by the trinucleotide: CH3-pU-pA-pA-CH3). Each model was simulated for a total of 5 ns (1 ns per value of lambda), with the data being collected over the final 750 ps of each trajectory. The free energy for each leg of the thermodynamic cycle presented as Figure 8 was calculated for the full ribozyme in the U-1 "Stacked" conformation, from 500 ps simulations at each of nine values for λ. A set of restraints reinforcing the hydrogen bonding interactions between the A2-G33 base pair and from residues A1 to C16:OP1 and C17:OP2 were applied to maintain the integrity of the active site throughout the simulation, in addition to restraining Dstack to between 3.5 Å and 7.0 Å. These restraints allowed for a more stable calculation by restricting the sampling to just the conformational basin (U-1 Stacked) predicted to be catalytically active. The difference in free energy between the models and the ribozyme microstates was then used to calculate the shifted (at 300 K) microscopic pKa values relative to the experimentally determined reference values for the nucleotides in solution:
From this data, the macroscopic "apparent pKa" values could then be estimated, by fitting a two state apparent pKa model27 to the active fraction, AH+EB−, as described in the text.
It should be noted that the thermodynamic cycle imposes the contraint that the legs of the cycle sum to zero, however, the free energy calculations performed here were done so independent of one another. As a result the calculated cycle did not close and instead summed to −0.53 ± 0.19 kcal/mol. In order to close the cycle, the residual energy (−0.53 kcal/mol) was distributed evenly across the four microstates and the apparent pKa model was fit to the active fraction generated from the "closed cycle" pKa values. The data for these free energy calculations and the resulting pKa values used in the fitting procedure are presented in Table S1.
We thank Dr. David M.J. Lilley and Dr. Timothy J. Wilson for discussion. Additionally, the authors are grateful for financial support provided by the National Institutes of Health (GM62248 to DY). This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number OCI-1053575 (project number TG-MCB110101 to DY), and also supported in part by the Blue Waters sustained petascale computing project (NSF OCI 07-25070 and PRAC OCI-1036208).
Additional details on residue numbering, crystal simulation, molecular dynamics simulations of the impact of U-1 stacking on the active site, and data used to predict apparent pKa values.
Competing Financial interests
The authors declare no competing financial interests.