Using Simulations and kinetic network models to reveal the dynamics and functions of Riboswitches

Using Simulations and kinetic network models to reveal the dynamics and functions of Riboswitches


Riboswitches, RNA elements found in the untranslated region, regulate gene expression by binding to target metaboloites with exquisite specificity. Binding of metabolites to the conserved aptamer domain allosterically alters the conformation in the downstream expression platform. The fate of gene expression is determined by the changes in the downstream RNA sequence. As the metabolite-dependent cotranscriptional folding and unfolding dynamics of riboswitches is the key determinant of gene expression, it is important to investigate both the thermodynamics and kinetics of riboswitches both in the presence and absence of metabolite. Single molecule force experiments that decipher the free energy landscape of riboswitches from their mechanical responses, theoretical and computational studies have recently shed light on the distinct mechanism of folding dynamics in different classes of riboswitches. Here we first discuss the dynamics of water around riboswitch, highlighting that water dynamics can enhance the fluctuation of nucleic acid structure. To go beyond native state fluctuations we used the Self-Organized Polymer (SOP) model to predict the dynamics of add adenine riboswitch under mechanical forces. In addition to quantitatively predicting the folding landscape of add-riboswitch our simulations also explain the difference in the dynamics between pbuE adenine- and add adenine-riboswitches. In order to probe the function in vivo we use the folding landscape to propose a system level kinetic network model to quantitatively predict how gene expression is regulated for riboswitches that are under kinetic control.


I Introduction and Scope of the Review

The remarkable discovery just over a decade ago that riboswitches, which are RNA elements in the untranslated regions in mRNA, control gene expression by sensing and binding target metabolites with exquisite sensitivity is another example of the versatility of RNA in controlling crucial cellular functions Serganov and Nudler (2013). In the intervening time, considerable insights into their functions have come from a variety of pioneering biochemical and biophysical experiments. In addition, determination of structures of a number of riboswitches has greatly aided in molecular understanding of their functions and in the design of synthetic riboswitches. Typically, riboswitches contain a conserved aptamer domain to which a metabolite binds, producing a substantial conformational change in the downstream expression platform leading to control of gene expression. The variability in the functions of structurally similar aptamer domain is remarkable. The add adenine (A) riboswitch activates translation upon binding the metabolite (purine) whereas the structurally similar pbuE adenine (A) riboswitch controls transcription Mandal et al. (2003). We can classify both these as ON riboswitches, which means that translation or transcription is activated only upon binding of the metabolite. In contrast, OFF riboswitches (for example Flavin Mononucleotide (FMN) binding aptamer) shut down gene expression when the metabolite binds to the aptamer domain. Finally, some of the riboswitches are under kinetic control (FMN riboswitch Wickiser et al. (2005) and pbuE adenine riboswitch Frieda and Block (2012)) whereas others (for example add adenine riboswitch and SAM-III riboswitch) may be under thermodynamic control.

The functions of riboswitches are vastly more complicated than indicated by in vitro studies, which typically focus on limited aspects of their activities. Under cellular conditions the metabolite which binds to the aptamer domain itself is a product of gene expression. Thus, regulation of gene expression involves negative or positive (or a combination) feedback. This implies that a complete understanding of riboswitch function must involve a system level description, which should minimally include the machinery of gene expression, rate of transcription or translation, degradation rates of mRNA, and activation rate of the synthesized metabolite as well as rate of binding (through feedback loop) to the aptamer domain. Many in vitro studies have dissected these multisteps into various components in order to quantify them as fully as possible. In this context, single molecule studies in which response of riboswitches to mechanical force are probed have been particularly insightful Greenleaf et al. (2008); Frieda and Block (2012).

More recently, computational and theoretical studies have been initiated to develop a quantitative description of the folding landscape and dynamics at the single molecule level Lin and Thirumalai (2008); Whitford et al. (2009); Allner et al. (2013); Quarta et al. (2012); Feng et al. (2011). The findings have been combined to produce a framework for describing the function of riboswitches at the system level. Because atomically detailed simulations can only provide limited information of the dynamics in the folded aptamer states it is necessary to develop suitable coarse grained (CG) model for more detailed exploration of the dynamics. The CG models of nucleic acids, first introduced by Hyeon and Thirumalai Hyeon and Thirumalai (2005), are particularly efficacious to deal with riboswitches that undergo large scale conformational fluctuations for functional purposes. The goal of this article is limited to a brief review of the insights molecular simulations have brought to the understanding of the folding landscape of riboswitches using purine-binding and S-adenosylmethionine (SAM) riboswitches as examples. We begin with the description of the fluctuations in the folded state of the small preQ1 riboswitch, which exhibits rich dynamics in the native state that can be probed using atomic detailed simulations in explicit water. The hydration dynamics could have a functional role, which can be resolved by spectroscopic experiments. We then describe the response of three classes of riboswitches to forces and map the entire folding landscape from which we have made testable predictions. The data from these free energy landscapes are used to construct a network model, which provides system level description of the OFF riboswitches (Fig.1).

Ii Hydration dynamics around the folded state: All atom simulations

It is known that unlike proteins there are many are several low free energy excitations (alternate structures) that a folded RNA can access. Consequently, dynamical fluctuations of the folded states are critical for RNA to execute their biological functions. There is growing evidence that hydration plays a key role in triggering conformational fluctuations in RNA. First, RNA can access low-lying excitation states via local melting of bases Jacob et al. (1987); Dethoff et al. (2012). Recent NMR studies suggest that a potential pre-melting of the hydration shell is required for the base pair disruption in response to elevation in temperature Rinnenthal et al. (2010); Nikolova and Al-Hashimi (2010). Second, the versatile functional capacity of RNAs can be attributed to their ability to access alternative conformations Zhang et al. (2006). Local conformational fluctuations from a few nanosecond dynamics enable RNA to explore a heterogeneous conformational ensemble, giving them the capacity to recognize and bind a diverse set of ligands. Third, binding of metabolites to riboswitches to control gene expression Montange and Batey (2008) may also be linked to local fluctuations in specific regions of co-transcriptional folded UTR regions of mRNA. In all of these examples, hydration of RNA is likely to play important role.

In order to illustrate the importance of hydration we performed atomically detailed simulations of PreQ riboswitch. We showed that water dynamics is spatially heterogeneous with metastable functionally relevant states whose dynamics and spans many orders of magnitude. This behavior is reminiscent of glassy behavior Thirumalai et al. (1989). The glassy behavior of water molecules may indicate that RNA molecule is able to access low-lying free energy states around the putative folded state. We identified distinct classes of water molecules near the RNA surface, which can be classified as “bulk”, “surface”, “cleft”, and “buried” water in the order of increasing water hydrogen bond relaxation time Yoon et al. (2014). In this section we review the molecular details of hydration around various regions of RNA and discuss how water dynamics gives rise to local structural fluctuations by using atomic simulations of PreQ-riboswitches Yoon et al. (2013, 2014).

Water hydrogen bond kinetics around nucleotides. The time- and ensemble-averaged auto-correlation function (see definition in Ref.Yoon et al. (2014); Luzar and Chandler (1996)) are used to quantify the structure and dynamics of water molecules near the surface of RNA. The relaxation kinetics of the water HB at T=310 K around three different nucleotide groups (B: base, R: ribose, P: phosphate) is fitted to a multiexponential function where with and denotes B, P, or R (Fig.2a, left). The time constant ranges from ps to ps, but 90 % of kinetics is described by the dynamics of ps (Fig.2a, left). The average lifetime of the water HB at each nucleotide group () reveals that water molecules exhibit the slowest relaxation dynamics near bases instead of phosphate groups as = 193 ps, = 81 ps, = 289 ps, and hence . It is worth pointing out that these time scales are much longer than found in proteins, and far exceed by a few orders of magnitude hydrogen bond dynamics in bulk water.

Excess monovalent counter-ions are distributed around RNA to neutralize the negative charges on the phospho-diester backbone of nucleic acids. The auto-correlation functions computed for Na ions bound to P, R, and B (Figure 2a, right) show that the time constant of ion relaxation is a few orders of magnitude greater than the water hydrogen bonds, = 294 ns, = 63 ns, = 9.2 ns. The order of lifetime differs from that of water HB as . In contrast to water HB, monovalent counterions have the slowest dynamics near the phosphate group. Most importantly, while binding or release of a Na ion to or from the surface of RNA certainly perturbs the water environment Song et al. (2014), the time scale separation between water and counterion dynamics ensures that the hydration dynamics around RNA occurs essentially in a static ionic environment.

Heterogeneity of Water Dynamics on the RNA Surface. The time scale of hydrated water varies many orders of magnitude depending on its location on the surface of RNA. Calculations of electrostatic potential on the solvent accessible surface confirm Yoon et al. (2014, 2013) that the charge distribution on RNA surface is indeed not uniform but heterogeneous. Multi-exponential function with different weights () and well-separated time constants () is needed to quantitate the relaxation dynamics of water molecules around four selected nucleotides of preQ-riboswitch, 24U, 29A, 33C, and 35A (Fig.2c). The rich dynamics reflects the heterogeneity and justifies the interpretation that there are distinct class of water molecules, which can be divided into multiple classes such as “bulk”, “surface”, “cleft”, and “buried” water Yoon et al. (2014). At high temperatures, the population of fast, bulk water-like dynamics is dominant, but as the temperature decreases, the population of slow dynamics grows. The average lifetime of water molecules near RNA is at least 1–2 orders of magnitude slower than that of bulk water over the broad range of temperatures (Fig.2c).

Water-induced fluctuations of base-pair dynamics. Dynamic feature of water that induce local conformational fluctuation of RNA is captured by probing the base pair dynamics along with surface water Yoon et al. (2013). The space made of base stacks and base pairings is generally dry and hydrophobic, and thus devoid of any water molecules. However,in base pairs located at the end of stacks, it is possible to observe an enhanced fluctuation of base pair. Figure 2c shows the dynamics of base pairs A3-U24 located in the 5’- and 3’-end in preQ riboswitch in aqueous solution. Remarkably, when the time series of water density around H3 of U24 and breathing dynamics of the base pair are compared, the change in water density always precedes the change in base-pair distance. The water densities calculated in the first and second solvation shell around H3 of AU24 show that water population starts to increase before the base pair disruption; the decrease of water population always precedes the event of base-pair formation. Thus, we conclude that the dynamics of water hydration and dehydration induces the breathing dynamics of base pairs. The spontaneous fluctuations in base pair opening induced by water are important in protein-DNA interactions as well and may be responsible for transcription initiation by RNA polymerase.

Iii Stability of isolated helices control the folding landscapes of purine riboswitches

A key event in the function of riboswitches is the conformational change in the aptamer domain leading to the formation of the terminator with the downstream expression platform (Fig.1a) or sequestration of the ribosome binding site upon ligand binding (Fig.1b). In order to assess the time scale in which such conformational change takes place and how it competes with ligand binding it is first important to quantitatively map the folding landscapes of riboswitches. From such landscapes the time scales for the conformational change in the switching region in the aptamer can be estimated Hyeon et al. (2008); Lin and Thirumalai (2008).

In a pioneering experiment Block and coworkers used single molecule pulling experiments to map the folding landscape as a function of the extension of the RNA. Purine (guanine and adenine) riboswitches are remarkably selective in their affinity for ligands and carry out markedly different functions despite the structural similarity of their aptamers. For the pbuE adenine (A) riboswitch, whose response to force was first probed in the LOT experiments, ligand binding activates the gene expression when an antiterminator is formed. In the absence of adenine, part of the aptamer region is involved in the formation of a terminator stem with the expression platform resulting in transcription termination. The add A-riboswitch activates the gene expression by forming a translational activator upon ligand binding. In the absence of adenine, the riboswitch adopts the structure with a translational repressor stem in the downstream region. At the heels of the first single molecule studies, we reported the entire folding landscape and calculated the time scale for switching of helix that engages in hairpin formation with the downstream sequence using the self-organized polymer (SOP) model Hyeon et al. (2006); Hyeon and Thirumalai (2007); Lin and Thirumalai (2008). As we show below comparison of the landscapes of these two riboswitches underscores the importance of the stability of the isolated helices in the assembly and rupture of the folded straucture.

Structures of purine riboswitch aptamers are characterized by a three-way junction consisting of P1, P2 and P3 helices, which are further stabilized by tertiary interactions in the folded state (Fig.3a). For pbuE A-riboswitch, binding of metabolite (adenine) activates the gene expression by enabling the riboswitch to form an antiterminator. Without adenine, the molecule forms a terminator stem with the expression platform, resulting in transcription termination. On the other hand, the add A-riboswitch uses adenine to regulate the process of translation. Recent single molecule experiments Greenleaf et al. (2008); Neupane et al. (2011) and our simulation studies Lin and Thirumalai (2008); Lin et al. (2014) have shown that, despite the marked structural similarity, these two aptamers have different folding landscapes, thus providing a fingerprint of their function.

Single molecule optical tweezer experiments have been used to directly observe the hierarchical folding of both pbuE A- and add A-riboswitch aptamers Greenleaf et al. (2008); Neupane et al. (2011). Here, we summarize force ()-triggered unfolding and refolding of the A-riboswitch aptamer theoretically using Brownian dynamics simulations of the SOP model Hyeon et al. (2006); Hyeon and Thirumalai (2007). The crystal structure of add A-riboswitch (PDB id: 1Y26 (U17 to A79)) is available while that of pbuE A-riboswitch is not. However, since the sequence similarity between add-A and pbuE A-riboswitch is unusually high, we modeled the atomic structure of pbuE A-riboswitch by substituting the sequences of pbuE A-riboswitch into the crystal structure of add A-riboswitch and produced an ensemble of pbuE A-riboswitch structures via conformational sampling with molecular dynamics simulations Lin et al. (2014).

In the absence of adenine our simulation show that force-induced unfoldings of both pbuE-A and add A-riboswitches occur in three distinct steps. Force extension curves of riboswitch generated under constant loading condition ( pN/s) reveal three distinct steps for both RS. Investigating the loss of secondary and tertiary contacts during the unfolding process, we found that the order of unfolding events differs qualitatively in add A-riboswitch and pbuE A-riboswitch. In add A-riboswitch, the unfolding occurred in the order of P1P2/P3P3U. The order of forced unfolding of pbuE A-riboswitch is P1P2/P3P2U, where P2/P3 denotes the disruption of kissing loop interaction between P2 and P3 due to force. In the absence of adenine thermal fluctuations transiently disrupt this kissing-loop interaction, which is consistent with the observation that stable P2/P3 tertiary interactions require adenine.

The presence of adenine in the binding pocket in the triple-helix junction of add A-riboswitch changes the force-response of RS completely: (i) The unfolding force increases from pN to pN, the value of which is comparable to the one found in experiments for the pbuE A-riboswitch aptamer Greenleaf et al. (2008); and (ii) the unfolding of RS occurs in all-or-none fashion without intermediate unfolding steps. After the complete unfolding, when refolding of the add A-riboswitch is initiated by reducing the force, we find that the refolding pathway follows the reverse order of unfolding pathway as UP3P2P2/P3P1. Refolding of P3 preceding that of P2 implies that P3 is more stable than P2, which is consistent with the implication from the stability of each helix ( kcal/mol kcal/mol) calculated using the Vienna RNA package Hofacker (2003).

Remarkably, despite the structural similarity between pbuE-A and add A-riboswitch aptamers, experiments show that P2 in pbuE unfolds at the last moment, which implies that P2 is the first structural element to refold upon force quench (or reduction). In agreement with the experiments, our results also imply that P2 ought to be more stable than P3 in the pbuE A-riboswitch aptamer, and Vienna RNA package indeed predicts that the stability of P2 is lower than that of P3 by 2 kcal/mol ( kcal/mol kcal/mol) (Fig.3). The difference of the stability in the two RS aptamers explains the reversed order of the folding of P2 and P3 in the pbuE A-riboswitch aptamer. The order of unfolding of the helices, which is in accord with single molecule pulling experiments, is determined by the relative stabilities of the individual helices. Our results show that the stability of isolated helices determines the order of assembly and response to force in these non-coding regions. Thus, the folding landscape is determined by the local stability of the structural elements, a finding that also holds good in the thermal refolding of a number of RNA pseudoknots Cho et al. (2009).

Based on the stability hypothesis as the determining factor of the RNA folding landscape, we make an interesting prediction for pulling experiments in a mutant of the add A-riboswitch. One of the major differences that contributes the different stability of P2 helix in the two purine riboswitches is that the P2 of add A-riboswitch has one G-U and two G-C base pairs, whereas P2 in the pbuE A-riboswitch has three G-C base pairs (see Fig.3a and Fig.4a). At the level of stability of secondary structure, a point mutation of U28C in the add A-riboswitch, which leads to three G-C base pairs in P2, would increase the secondary structure stability of P2 to 7.3 kcal/mol. Thus, the U28C mutation stabilizes the add A-riboswitch P2 by 1.1 kcal/mol lower than P3. The folding landscape of the U28C add A-riboswitch would be qualitatively similar to the WT pbuE A-riboswitch. As a consequence, we predict that U28C would reverse the order of unfolding of the add A-riboswitch.

It is noteworthy that the number of contacts between P2 and P3 hairpin loops with and without adenine is almost identical, which suggests that adenine binding does not affect the interactions between P2 and P3 hairpin loops. Rather, the adenine binding stabilizes the triple-helix junction and makes unfolding of P1 more difficult. Hence, the rate limiting step in the fully folded aptamer is the formation of P1.

The free energy profiles (), obtained from the distribution of the molecular extension at force () using , make explicit the hierarchical characteristics of RNA assembly and disassembly. In the absence of adenine, the position of first transition barrier from the folded state ( nm) is nm. Thus, the unfolding transition over the first barrier amounts to the unzipping of three base pairs in P1 helix. In the presence of adenine, the position of the first barrier shifts to nm, implying that five base pairs of P1 in direct interaction with adenine should be disrupted at the first TS. Thus, adenine binding makes the first unfolding barrier the rate limiting step. We also find that, at pN, binding of adenine stabilizes the folded state by 6 and increases the energy barrier for leaving the folded state by 2 . These results support the hypothesis that the ligand binding stabilizes P1, a feature that is common to most riboswitches.

The binding of adenine stabilizes the folded basin of attraction, which slows down the unfolding transition by two orders of magnitude. The slow unfolding rate, which make the conformational sampling difficult, is in agreement with that observed for FMN riboswitches Wickiser et al. (2005), suggesting that functions of riboswitches are kinetically, rather than thermodynamically, controlled (see below for further discussion). The result is crucial for transcription of the complete riboswitch.

Iv Folding landscapes of SAM riboswitch. Is SAM riboswitch under thermodynamic control?

Upon binding S-adenosylmethionine (SAM) the riboswitch undergoes undergoes an allosteric transition to control translation. There are at least five distinct classes of riboswitches that bind SAM or its derivative S-adenosylhomocysteine (SAH). One of them, the SAM-III riboswitch Fuchs et al. (2006) in the metK gene (encodes SAM synthetase) from Lactobacillales species, inhibits translation by sequestering the Shine-Dalgarno (SD) sequence (Fig.5a) Fuchs et al. (2006). When SAM is bound, the somewhat atypical SD sequence (GGGGG shown in the blue shaded area in Fig.5a) is sequestered by base pairing with the anti-Shine-Dalgarno (ASD) sequence, thus hindering the binding of the 30S ribosomal subunits to mRNA. In the absence of SAM, the sequence outside the binding domain enables the riboswitch to adopt an alternative folding pattern, in which the SD sequence is exposed and free to engage the ribosomal subunit Lu et al. (2011). Thus, SAM-III is an OFF switch for translation, in contrast to add A-riboswitch, which is an ON switch. Translation control is determined by competition between the SD-ASD pairing and loading of ribosomal subunit onto the SD sequence. In order to function as a switch, the SD sequence has to be exposed for ribosome recognition, which implies that at least part of the riboswitch structure accommodating SAM III has to unfold (Fig.5). These considerations prompted us to quantitatively determine the folding landscape and the rates of conformational transitions between the ON and OFF states of the SAM-III riboswitch Anthony et al. (2012).

In order to understand if SAM-III functions under thermodynamic control we calculated the force-extension () or FEC as well as constant -dependent free energy profiles as a function of for SAM-III. The FEC (black curve in Fig.5b) shows that there are two intermediate states, with extension nm and nm in the absence of SAM. Helices P1 and P4 rupture in a single step at 9 pN and P2 unravels at pN. Helix P3 unfolds fully only when 12pN. Interestingly, when SAM is bound, the riboswitch unfolds in an apparent all-or-none manner at pN (Fig.5b). The distribution (blue curve in Fig. 5A), shows the presence of two intermediate states ahead of global unfolding. The nm peak corresponds to rupture of P1 and P4 helices. P3 unfolds in the later stages creating a peak at nm. The hierarchical unfolding pathway of SAM-III riboswitch is , where means helices P1 and P4 are ruptured, and P2 represents additional unfolding of P2. The observed order of unfolding is also reflected in the rupture of contacts, a more microscopic representation of unfolding dynamics (Fig. 5c). The intermediate states in Fig.5b can be traced to the breaking of contacts within the helices. Just as in Fig. 5b the order of contact rupture corresponds to the order in which the helices unfold in the FEC in Fig.5b.

Using simulations at constant force we also calculated the free energy profiles using where is the probability distribution of . At pN and in the absence of SAM the riboswitch is in the folded basin of attraction (Fig. 5d). The free energy profile in Fig. 5d also shows that binding of SAM consolidates the formation of helix P1 and P4, further stabilizing the folded state. At pN, SAM binding stabilizes the folded state by , and increases the energy barrier for leaving the folded state by . The distance from the folded state to the first barrier in the absence of SAM is nm, which indicates unzipping of 2.5 base pairs, assuming a contour length increase of 0.4 nm/nt. In the presence of SAM, the position of the first barrier shifts to nm, implying that 4 base pairs of P1 next to the nucleotide G48 (Fig.5a) that has direct contacts with SAM are ruptured at the transition state. Thus, disruption of contacts with SAM becomes the key barrier in the first unfolding step, and must be an important step in translational regulation.

V Is SAM riboswitch under thermodynamic control?

As shown in Fig.6a, the kinetic processes in riboswitches that control transcription are determined by a number of time scales. In the transcription process, the ability to function as an efficient switch depends on an interplay of the time scales: (i) metabolite binding rate (), (ii) the folding times of the aptamer (), (iii) the time scales to switch and adopt alternate conformations with the downstream expression platform (), and (iv) the rate of transcription. In “OFF” riboswitches that shut down gene expression upon metabolite binding, a decision to terminate transcription has to be made before the terminator is synthesized, which puts bounds on the metabolite concentration, and the aptamer folding rate (). For simplicity, can serve as a simple criterion to determine whether the co-transcriptional folding of riboswitches is under thermodynamic or kinetic control. In the limit transcript synthesis is faster than the equilibration time of the riboswitch conformation. For typical values of these parameters in both FMN and pbuE A-riboswitch, efficient function mandates that the riboswitches be under kinetic control, which implies that the “OFF” and “ON” states of riboswitch are not in equilibrium.

In contrast, the function of SAM-III, which controls translation, is different. The major time scales that control the function of SAM-III RS, and those that regulate translation in general, are (i) bimolecular binding rate of SAM to RS (), (ii) dissociation rate of SAM from the riboswitch complex (), (iii) the rate of mRNA degradation (). Thus, the only clear physical bound on the function of SAM-III is that binding of metabolite should occur multiple times before the mRNA degradation, which leads to , where is the concentration of SAM. Typical values of , and requires that nM. It is worth pointing out that our estimates of folding and unfolding times based on simulations at low forces, and other time scales are all much less than sec, which sets the longest time for translational control. Hence, the multiple transitions between the OFF and ON states can occur before mRNA is degraded, which gives additional credence to the argument that the function of SAM-III is under thermodynamic control. There is a caveat to this conclusion. It is known that in bacteria transcription and translation are coupled, which is likely to complicate our arguments. In order to provide a complicate description we require a network model that includes transcription-translation coupling. Because is small it is still possible that the SAM riboswitch could be under thermodynamic control.

Vi Kinetic Network Model of gene regulation and the role of negative feedback in control of transcription

As stated earlier, gene expression is mediated by binding of metabolites to the conserved aptamer domain, which triggers an allosteric reaction in the downstream expression platform. However, the target metabolites are usually the products or their derivatives of the downstream gene that the riboswitches control. Hence, metabolite binding to riboswitches serves as a feedback signal to control RNA transcription or translation initiation. The feedback through metabolite binding is naturally designed to be a fundamental network motif for riboswitches. In ON-riboswitches, metabolite binding thus stabilizes the aptamer structure during transcription and prevents the formation of the terminator stem before transcription is completed (pbuE A-riboswitch) or the formation of translation repressor stem before translation is initiated (add A-riboswitch). Whereas in OFF-riboswitch, metabolite binding shuts down the gene expression by promoting the formation of terminator stem (see Fig.6). In order to understand the in vivo riboswitch we developed a kinetic network model taking into account the interplay between the speed of RNA transcription, folding kinetics of the nascent RNA transcript, and the kinetics of metabolite binding to the nascent RNA transcript, and the role of feedback arising from interactions between synthesized metabolities and the transcript. The effects of speed of RNA transcription and metabolite binding kinetics have also been investigated experimentally in vitro in an insightful study involving the flavin mononucleotide (FMN) riboswitches (Wickiser et al., 2005). They argued that FMN riboswitch is kinetically driven implying that the riboswitch does not reach thermodynamic equilibrium with FMN before a decision between continued transcription and transcription termination needs to be made. The mathematical solution of the kinetic network model, which uses as partial input the rates of switching obtained from the folding landscapes, show that in general riboswitches that control transcription are under kinetic control Lin and Thirumalai (2012). A brief summary is presented here.

Efficient function of RS, implying a large dynamic range (quantified by response of the RS to varying metabolite concentration) without compromising the requirement to suppress transcription or translation, is determined by a balance between the transcription speed (), the folding and unfolding rates of the apatmer ( and ), and the binding and unbinding rates of the metabolite ( and , where is the metabolite concentration). In order to capture the physics behind the dynamics it is necessary to consider kinetic network model describing the coupling between aptamer dynamics and transcription. In Fig.6 demonstrating the kinetic network model for the transcription regulation by OFF-riboswitches, the upstream of the protein-coding gene consists of sequences involving the transcriptions of aptamer (B), antiterminator (B) and terminator (B) of the riboswitch. The transcription initiation is followed by elongation, folding of the RNA transcript, and metabolite binding. RNA polymerase first transcribes the aptamer (B), and moves on to the synthesis of the RNA transcript for anti terminator B at a rate of , and terminator sequence at , resulting in the production of the regulatory region of RNA. is the transcript with the sequence of the protein-coding region starting to be transcribed, and eventually grows to , the full protein-coding region transcribed, with a rate of . During the process of transcription elongation, each of the transcript states, B and B, can form states with the aptamer domain folded (B and B) with a folding rate of and , respectively. The folded aptamers bound with metabolite (M) are BM and BM with binding rate constant and , respectively. The transcripts in state B and BM can further elongate until the terminator sequence is transcribed with their expression platform forming a transcription terminator stem and dissociate from the DNA template with a rate of , forming B and BM. The fraction of transcription termination, , is determined from the amount of the terminated transcripts (in green block) relative to nonterminated transcripts (in blue block). The activated metabolte (M), produced from protein and activated by the enzyme () encoded by the gene OF, can bind to the folded aptamer and can abort transcription, which imposes a negative feedback on the transcription process.

For a riboswitch to function with a large dynamic range, transcription levels should change significantly as the [M] increases from a low to high value. (i) In the high [M], RNA transcript in the aptamer folded state binds a metabolite with . In FNM-riboswitch, small value results in the formation of a terminator stem, which subsequently terminates transcription. (ii) In the low [M] limit, the aptamer folded state is mostly unbound and can remain folded until transcription termination or can fold to the antiterminator state, enabling the synthesis of full RNA transcript. The levels of transcription termination are thus controlled by the transition rates between the aptamer folded and unfolded states (, , , ). Equilibrium between B and B can be reached only if the rate of transcription is much slower than the rates of folding/unfolding and metabolite binding. By varying the rate of transcription, which can be experimentally realized by adding transcription factors such as NusA Zhou et al. (2011), it may be possible to drive the cotranscriptional folding of riboswitch from kinetic to thermodynamic control. However, for realistic values of the various rates in Fig.6 we predict that transcription in vivo is under kinetic control.

In the presence of a negative feedback loop the concentration of target metabolites is also regulated by gene expression. Under nominal operating conditions () binding of target metabolites, products of the downstream gene that riboswitches regulate, significantly suppresses the expression of proteins. Negative feedback suppresses the protein level by about half relative to the case without feedback. In vivo, the presence of RNA binding proteins, such as NusA Zhou et al. (2011), may increase the pausing times, thus effectively reducing the transcription rates. Thus, the repression of the protein level by the riboswitch through metabolite binding may be up to 10-fold. Faster RNA folding and unfolding rates than those we obtained may also increase the suppression by negative feedback and broaden the range of transcription rates over which maximal suppression occurs. These predictions are amenable to experimental test.

In response to changes in the active operon level, the negative feedback speeds up the response time of expression and modestly reduces the percentage change in the protein level relative to change in the operon level. The steady-state level of expression for autoregulation varies as a square root of the DNA concentration. Adaptive biological systems may minimize the variation in gene expression to keep the systems functioning normally even when the environments change drastically. One may need to consider more complex networks than the single autoregulation in the transcription network to find near perfect adaptation to the environmental change Ma et al. (2009).

The effect of negative feedback accounting for the binding of metabolites, which themselves are the product of genes that are being regulated. Our previous work showed that because of the interplay of a number of time scales determining the riboswitch function at the system there are many scenarios that can emerge, which can be encapsulated in terms of a dynamic phase diagram. An example dynamic phase diagram (for a full discussion see Lin and Thirumalai (2012)) in terms of the transcription rates , , , illustrates the complexity of the transcription process. The interplay between folding of RNA transcripts, transcription, and metabolite binding regulate the expression of , which can be quantified using the production of the protein, , on the transcription rates and the effective binding rate . The dynamic phase diagram in Fig. 4B, calculated by varying both () and with the equilibrium binding constant of the metabolite to the aptamer fixed to nM, a value that is appropriate for FMN Wickiser et al. (2005). We expect that after the aptamer sequence is transcribed, the formation of the aptamer structure is the key step in regulating transcription termination. Thus, regulation of should be controlled by the folding rate ,the effective metabolite binding rate, and for regulation of . Figure 4B shows three regimes for the dependence of on and . In regime I, , the folding rate is slow relative to transcription to the next stage (Fig. 1), which implies that the aptamer structure does not form on the time scale set by transcription. The dominant flux is from to , which leads to high probability of fully transcribed RNA downstream because of the low transition rate from to . The metabolite binding has little effect on protein expression in this regime, particularly for large , and hence the protein is highly expressed. In regime II, , the aptamer has enough time to fold but metabolite binding is slow. The dominant flux is , leading to formation of antiterminator stem () or transcription termination (). The expression level of protein is thus mainly determined by and , and the protein production is partially suppressed in this regime. In regime III, and , the aptamer has sufficient time to both fold and bind metabolite, the dominant pathway is , leading to transcription termination. The protein production is highly suppressed in this regime. The arrow shows that for parameters that are appropriate for FMN riboswitch (see Tables 1 and Table 2 in Lin and Thirumalai (2012)( fall on the interface of regime I and regime III. The metabolite binding fails to reach thermodynamic equilibrium due to low dissociation constant. However, the effective binding rate is high because the steady state concentration of metabolites ( 25 M) is in large excess over RNA transcripts. Thus, the riboswitch is kinetically driven under this condition even when feedback is included.

Vii Concluding Remarks

Based on our previous works, we have provided broad overview, from atomic scale to systems level, of how the complex dynamics of riboswtiches emerge depending many inter-related rates. At the atomic scale, dynamics of the surface water around riboswitches plays critical role in inducing the local fluctuation in the riboswtich. At the level of single riboswitch, we have shown that explicit simulations of riboswitches, in conjunction with single molecule experiment, is a powerful tool to understand the conformational dynamics of riboswtich both with and without metabolites. At the systems level, in which the minimal model of cellular environment is considered, the dynamics of riboswitches in isolation are modulated by a number of factors, which is made explicit by the kinetic network model of FMN riboswitch. Our collective works show that combination of theory, experiments, and simulations are needed to understand the function of riboswitches under cellular conditions.

Riboswitches also provide novel ways to engineer biological circuits to control gene expression by binding small molecules. As found in tandem riboswitches Breaker (2008); Sudarsan et al. (2006), multiple riboswitches can be engineered to control a single gene with greater regulatory complexity or increase in the dynamic range of gene control. Synthetic riboswitches have been successfully used to control the chemotaxis of bacteria Topp and Gallivan (2007). Our study provides a physical basis for not only analyzing future experiments but also in anticipating their outcomes.

Acknowledgements: This work was supported by grants from the National Institutes of Health (GM 089685) and the National Science Foundation (CHE 13-61946).

Figure 1: Principles of transcription and translation regulation by OFF-riboswitches. (A) Metabolite binding to the aptamer domain results in the formation of the terminator hairpin exposing a stretch of Uracil nucleotides. The polymerase disengages from the gene, thus terminating transcription. (B) Binding of the metabolite encrypts the Ribosome Binding Site (RBS) preventing the ribosome to initiate translation. The figure adapted from Ref.Kim and Breaker (2008).
Figure 2: a. The data for the dynamics of water HBs (left) and Na (right) with base, ribose, and base (left) are described by the appropriate auto-correlation functions. b. Temperature dependence of water HB relaxation time of bulk water and around four nucleotides, 24U, 29A, 33C, and 35A or PreQ. Arrhenius fits are made to K for the four nucleotide and bulk water. Alternatively, the relaxation dynamics of bulk water HB can be fit using the Vogel-Fulcher-Tamman (VFT) equation over the full range of temperatures (dashed line) Yoon et al. (2014). c. Water dynamics induced fluctuations of a base pair. The status of base pair (N1–N3), quantified by calculating logistic function as well as distance , shows apparent fluctuations, whose time scale is 10 ns. Figures adapted from Refs.Yoon et al. (2013, 2014)
Figure 3: Force-induced dynamics of add A-riboswitch (RS). a. Structure of the conserved domain of purine riboswitch containg a three way junction. On the left is the secondary structure map and on the right the three dimensional structure is shown. b. Force-extension curves (FECs) obtained by pulling the RS at loading rate of 960 pN/s in the presence (left) and absence (right) of metabolite. The FEC in red on the right panel was obtained during the refolding of the RS while the exerted force is reduced. c. Free energy profile with (red) and without (blue) the metabolite. d. Force-dependent transition rates. The data points are directly from simulation; The lines were obtained by calculating mean first passage time using with a force-independent diffusion constant, which was calibrated by equating the theoretical and simulation rates. Figure adapted from Ref.Lin and Thirumalai (2008)
Figure 4: pbuE adenine-riboswitch. a. Secondary structure map (on the left) and the three dimensional structure on the right. b. Free energy profile calculated from simulations at , 13 pN. Figure adapted from Ref.Lin et al. (2014)
Figure 5: Dynamics of SAM-III riboswitch under force. a. Structure of SAM-III RS. The blue shaded area on the left indicates the Shine-Dalgarno sequence recognized by the ribosome. b. Simulated force-extension curve of SAM-III riboswitch in the absence of metabolite (black) produced at pN/s. The distribution of molecular extension () during the pulling simulation is shown in blue at the bottom. FEC in red was produced in the presence of metabolite at the binding pocket. c. Average number of contacts in each helix from P1 to P4. d. Free energy profile at zero force calculated from streching simulation with and without metabolite (SAM) in the binding pocket. e. Transition rates between F and P2/P3 states at varying forces. The data points are from explicit simulations. The lines were obtained from mean first passage time calculation on . Figure adapted from Ref.Lin and Thirumalai (2013)
Figure 6: a. Kinetic network model for transcription regulated by “OFF”-riboswitches. b. Dependence of protein production on the network parameters with “negative feedback” ( , , , , , , , , , , , , , , , and . (Top) Protein levels (color coded) as functions of and ). The dependence of on and is categorized into three regimes. Points on the dashed line separating regime II and regime III satisfy . The major pathway in the transcription process in each regime is shown on the right. The arrow indicates the data point from the value of and . (Bottom) as functions of and . Points on the dashed line separating regime I and II/III satisfy . The data point corresponding to the arrow results from using the value of and . Figure adapted from Ref.Lin and Thirumalai (2012)


  1. A. Serganov and E. Nudler, Cell 152, 17 (2013).
  2. M. Mandal, B. Boese, J. E. Barrick, W. C. Winkler, and R. R. Breaker, Cell 113, 577 (2003).
  3. J. K. Wickiser, W. C. Winkler, R. R. Breaker, and D. M. Crothers, Molecular cell 18, 49 (2005).
  4. K. L. Frieda and S. M. Block, Science 338, 397 (2012).
  5. W. J. Greenleaf, K. L. Frieda, D. A. N. Foster, M. T. Woodside, and S. M. Block, Science 319, 630 (2008).
  6. J. Lin and D. Thirumalai, J. Am. Chem. Soc. 130, 14080 (2008).
  7. O. Allner, L. Nilsson, and A. Villa, RNA 19, 916 (2013).
  8. G. Quarta, K. Sin, and T. Schlick, PLOS Comp. Biol. 8, e1002368 (2012).
  9. J. Feng, N. G. Walter, and C. L. Brooks, III, J. Amer. Chem. Soc. 133, 4196 (2011).
  10. P. C. Whitford, A. Schug, J. Saunders, S. P. Hennelly, J. N. Onuchic, and K. Y. Sanbonmatsu, Biophys. J. 96, L7 (2009).
  11. C. Hyeon and D. Thirumalai, Proc. Natl. Acad. Sci. U.S.A. 102, 6789 (2005).
  12. W. F. Jacob, M. Santer, and A. E. Dahlberg, Proc. Natl. Acad. Sci. U.S.A. 84, 4757 (1987).
  13. E. A. Dethoff, J. Chugh, A. M. Mustoe, and H. M. Al-Hashimi, Nature 482, 322 (2012).
  14. J. Rinnenthal, B. Klinkert, F. Narberhaus, and H. Schwalbe, Nucleic Acids Research 38, 3834 (2010).
  15. E. N. Nikolova and H. M. Al-Hashimi, RNA 16, 1687 (2010).
  16. Q. Zhang, X. Sun, E. D. Watt, and H. M. Al-Hashimi, Science 311, 653 (2006).
  17. R. K. Montange and R. Batey, Annu. Rev Biophys. 37, 117 (2008).
  18. D. Thirumalai, R. D. Mountain, and T. R. Kirkpatrick, Phys. Rev. A. 39, 3563 (1989).
  19. J. Yoon, J.-C. Lin, C. Hyeon, and D. Thirumalai, J. Phys. Chem. B 118, 7910 (2014).
  20. J. Yoon, D. Thirumalai, and C. Hyeon, J. Am. Chem. Soc. 135, 12112 (2013).
  21. A. Luzar and D. Chandler, Phys. Rev. Lett. 76, 928 (1996).
  22. J. Song, J. Franck, P. Pincus, M. W. Kim, and S. Han, J. Am. Chem. Soc. 136, 2642 (2014).
  23. C. Hyeon, G. Morrison, and D. Thirumalai, Proc. Natl. Acad. Sci. U.S.A. 105, 9604 (2008).
  24. C. Hyeon, R. I. Dima, and D. Thirumalai, Structure 14, 1633 (2006).
  25. C. Hyeon and D. Thirumalai, Biophys. J. 92, 731 (2007).
  26. K. Neupane, H. Y. Daniel, A. N. Foster, F. Wang, and M. T. Woodside, Nucleic Acids Res. 39, 7677 (2011).
  27. J.-C. Lin, C. Hyeon, and D. Thirumalai, Phys. Chem. Chem. Phys. 16, 6376 (2014).
  28. I. Hofacker, Nucleic Acids Res. 31, 3429 (2003), ISSN 0305-1048.
  29. S. Cho, D. Pincus, and D. Thirumalai, Proc. Natl. Acad. Sci. U.S.A. 106, 17349 (2009).
  30. R. T. Fuchs, F. J. Grundy, and T. M. Henkin, Nature Struct. Mol. Biol. 13, 226 (2006).
  31. C. Lu, A. M. Smith, F. Ding, A. Chowdhury, T. M. Henkin, and A. Ke, J. Mol. Biol. 409, 786 (2011).
  32. P. C. Anthony, C. F. Perez, C. García-García, and S. M. Block, Proc. Natl. Acad. Sci. U.S.A. 109, 1485 (2012).
  33. J.-C. Lin and D. Thirumalai, Biophys. J. 103, 2320 (2012).
  34. J. Zhou, K. S. Ha, A. La Porta, R. Landick, and S. M. Block, Molecular cell 44, 635 (2011).
  35. W. Ma, A. Trusina, H. El-Samad, W. A. Lim, and C. Tang, Cell 138, 760 (2009).
  36. R. R. Breaker, Science 319, 1795 (2008).
  37. N. Sudarsan, M. C. Hammond, K. F. Block, R. Welz, J. E. Barrick, A. Roth, and R. R. Breaker, Science 314, 300 (2006).
  38. S. Topp and J. P. Gallivan, J. Am. Chem. Soc. 129, 6807 (2007).
  39. J. N. Kim and R. R. Breaker, Biology of the Cell 100, 1 (2008).
  40. J.-C. Lin and D. Thirumalai, J. Am. Chem. Soc. 135, 16641 (2013).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description