Design Principles for Riboswitch Function

Scientific and technological advances that enable the tuning of integrated regulatory components to match network and system requirements are critical to reliably control the function of biological systems. RNA provides a promising building block for the construction of tunable regulatory components based on its rich regulatory capacity and our current understanding of the sequence–function relationship. One prominent example of RNA-based regulatory components is riboswitches, genetic elements that mediate ligand control of gene expression through diverse regulatory mechanisms. While characterization of natural and synthetic riboswitches has revealed that riboswitch function can be modulated through sequence alteration, no quantitative frameworks exist to investigate or guide riboswitch tuning. Here, we combined mathematical modeling and experimental approaches to investigate the relationship between riboswitch function and performance. Model results demonstrated that the competition between reversible and irreversible rate constants dictates performance for different regulatory mechanisms. We also found that practical system restrictions, such as an upper limit on ligand concentration, can significantly alter the requirements for riboswitch performance, necessitating alternative tuning strategies. Previous experimental data for natural and synthetic riboswitches as well as experiments conducted in this work support model predictions. From our results, we developed a set of general design principles for synthetic riboswitches. Our results also provide a foundation from which to investigate how natural riboswitches are tuned to meet systems-level regulatory demands.


Introduction
The breadth of function exhibited by biological systems provides a foundation from which to develop solutions to global challenges including sustainability, renewable energy production, material synthesis, and medical advancement. Underlying these systems-level functions are regulatory components that evaluate molecular information in the extracellular and intracellular environments and ultimately translate that information into phenotypic responses over varying time scales. The properties of individual regulatory components and genetic networks composed of these components are tuned to control critical functions, including survival in fluctuating environments [1,2], minimization of energy expenditure in metabolism [3,4], developmental fate assignment [5], and proper information transmission through regulatory cascades [6][7][8]. To more effectively approach the reliable construction of synthetic biological systems, it is critical to advance our understanding of the degree to which individual component properties are tuned in natural systems, the underlying mechanisms that support tuning of biological components, and the effect of tuned components on resulting systems-level functions.
Riboswitches are RNA-based regulatory components that mediate ligand control of gene expression. Natural riboswitches have been identified in all three kingdoms of life [9] and primarily function by sensing a variety of essential cofactors, amino acids, and nucleotides and regulating the expression levels of proteins in associated metabolic pathways [10]. Riboswitches typically exploit three properties of RNA to translate changes in ligand concentration to changes in the expression of a target protein: specific and high affinity ligand binding by aptamer sequences, formation of distinct functional conformations primarily dictated by base-pairing interactions, and diverse gene expression regulatory mechanisms based on the central location of mRNA in the process of gene expression. With the exception of the glmS ribozyme [11,12], natural riboswitches function through a general mechanism in which the RNA molecule can primarily adopt two conformations and ligand binding to the formed aptamer in one conformation biases partitioning toward the ligand-bound conformation. Each conformation is associated with differential regulatory activities such that increasing ligand concentrations either increase (ON behavior) or decrease (OFF behavior) gene expression depending on which conformation contains the formed aptamer. Synthetic riboswitches have been constructed based on this functional mechanism to expand on the regulatory potential exhibited by natural riboswitches [13,14]. There has been significant interest in engineering riboswitches as tailored ligandresponsive genetic control elements by integrating aptamers selected in vitro [15] against diverse molecular ligands appropriate for different applications.
Natural and synthetic riboswitches have been demonstrated to be highly tunable regulatory components. Targeted nucleotide changes in synthetic riboswitches can shift the response curve [16][17][18][19][20]. Studies of natural riboswitches functioning through transcrip-tional termination found that the time lag between transcription of the aptamer and the terminator stem can tune the effective ligand concentration at which a half-maximal response is achieved (EC 50 ) [21,22]. These previous studies explored tuning in limited contexts by focusing on one aspect of riboswitch function (EC 50 ) for one type of regulatory mechanism. However, advancing the characterization or design of new riboswitches requires a general quantitative framework that applies broadly to different regulatory mechanisms. Due to the link between RNA secondary structure and function and the relative ease with which RNA molecules can be modeled, riboswitches present an interesting class of regulatory components through which researchers can examine links between physical composition, tuned component response properties, and resulting systems-level behavior.
In this study, we employed mathematical modeling to explore how the dynamics of riboswitch function dictate its performance, where performance is evaluated based on the response curve quantitatively linking ligand concentration and protein levels. To draw general conclusions regarding riboswitch performance, we considered three representative regulatory mechanisms: transcriptional termination, translational repression, and mRNA destabilization. Parameter space for all three mechanisms was surveyed in order to understand the relationship between model parameters and riboswitch performance. Our results show that the competition between reversible and mechanism-specific irreversible rate constants primarily dictates riboswitch performance and response curve tuning properties. Complete dominance of irreversible rate constants renders a riboswitch non-functional, although ligand binding during transcription can rescue riboswitch performance. Our results also demonstrate that placing an upper limit on the ligand concentration alters the observed tuning properties such that a maximum dynamic range exists for intermediate conformational partitioning. Model predictions are supported by published experimental data and new data obtained through the modification of a synthetic riboswitch. We provide a set of design principles for the construction of synthetic riboswitches based on our modeling results. In addition, our results lend insights into the inherent flexibility and potential biological relevance of tuning of natural riboswitches.

Kinetic modeling of riboswitch function
We started with a detailed molecular description of riboswitch function ( Figure S1) that accounts for folding and ligand binding during discrete steps of transcription. Three regulatory mechanisms were considered: translational repression, transcriptional termination, and mRNA destabilization. Translational repression occurs through ribosome binding site (RBS) sequestration in a double-stranded secondary structure that prevents ribosome recruitment. Transcriptional termination occurs through a rhoindependent mechanism such that hairpin formation directly upstream of a polyuridine stretch induces dissociation of the transcript from the template and the polymerase. We also considered the regulatory mechanism of a recently-described synthetic riboswitch that undergoes ribozyme self-cleavage [19], thereby initiating mRNA destabilization [23]. In these examples, two inter-converting conformations (A/B) are associated with differential protein levels subject to the specified regulatory mechanism. Ligand binding to the formed aptamer harbored in B promotes conformational stabilization, thereby increasing the prevalence of B.
We assigned a rate constant to each mechanistic step in the models to yield a detailed relationship between ligand concentration (L) and protein levels (P). In all models, transcriptional initiation produces a partial-length riboswitch in either conformation A (k fA ) or B (k fB ) to reflect transcriptional folding. Transcription is broken into discrete steps that represent different sequence lengths. Each step determines the extent of conformational switching (k 1 , k 1 9), the ability to bind and release ligand (k 2 , k 2 9), and the rate of progression to the next step (k E ). For transcriptional termination, riboswitches effectively choose between termination (k TA , k TB ) and extension (k MA , k MB ) after transcription of the terminator stem. To ensure that both conformations make the decision with the same frequency, we set the sum of termination and extension rate constants for each conformation equal to a single parameter k M : Following transcription of the full-length riboswitch for translational repression and mRNA destabilization or extension through the terminator stem for transcriptional termination, the transcript can be translated into protein (k PA , k PB ) or undergo degradation (k dMA , k dMB ). A single constant is assigned when the rate constants are equal between conformations (k P , k dM ). Values for the rate constants can vary widely depending on the organism and regulatory mechanism (Table 1). Therefore, we explored how each rate constant contributes to riboswitch performance. Riboswitch performance was evaluated under steady-state conditions for both ON and OFF behaviors by calculating a collection of performance descriptors that define the response curve ( Figure 1A): EC 50 , dynamic range (g) defined as the difference between high and low protein levels, basal protein levels (P(L = 0)), and ligand-saturating protein levels (P(LR')). While the dynamic range can be alternatively defined as the ratio of high and low protein levels, we selected the difference definition based on

Author Summary
Riboswitches are RNA-based components that integrate ligand binding and gene regulation to dynamically respond to molecular signals within cells. Natural riboswitches are employed to regulate metabolism and other cellular processes, while synthetic riboswitches have been constructed to expand the sensory and regulatory capabilities exhibited in nature. Characterization studies have revealed that sequence modifications can tune properties of the riboswitch response curve, which links ligand concentration to expression levels. Tunability is critical when matching component properties to the regulatory demands of biological systems; however, the characterization of riboswitch tuning strategies is complicated by the integration of numerous regulatory mechanisms and various processes, such as RNA folding and turnover, that impact riboswitch performance. To develop a generalized framework for examining quantitative aspects of riboswitch tuning, we modeled the kinetics of riboswitch function operating under common regulatory mechanisms. Our results reveal that riboswitch performance is primarily dictated by the competition between reversible and mechanism-specific irreversible rate constants. We demonstrate that practical system restrictions can significantly alter the requirements for riboswitch performance, necessitating a variety of tuning strategies. We developed design principles to guide the construction of synthetic riboswitches and a quantitative framework from which to investigate how natural riboswitches are tuned to meet systems-level regulatory demands.
the mathematical symmetry between the equations representing ON and OFF behaviors (Text S1).
Transcription can be considered as a discrete multistep process ( Figure S1). The conformations that can form at each step depend on the ordering of elements along the riboswitch sequence, such as the relative location of the aptamer or gene regulatory elements. Matching the number of steps and parameter values to particular sequence configurations becomes burdensome and restricts the elucidation of general principles. Therefore, we simplified the transcription process by assuming that synthesized transcripts appear in either conformation A or B and are immediately subject to conformational partitioning, ligand binding, and the regulatory mechanism ( Figure 1B-D). As a result, the outcome of the transcription process is reflected by biased folding into either conformation A (k fA ) or conformation B (k fB ). This simplification excludes ligand binding during transcription, which has been demonstrated for natural riboswitches functioning through transcriptional termination [21,22]. Therefore, we separately accounted for ligand binding during transcription in our analyses.

Competition between reversible and irreversible rate constants suggests three operating regimes
We first derived expressions for the performance descriptorsdynamic range, EC 50 , and basal and ligand-saturating levels-for riboswitches functioning through each regulatory mechanism (Text S2). From these derivations two common parameters emerged: the partitioning constant (K 1 = k 1 9/k 1 ) and the ligand association constant (K 2 = k 2 /k 2 9). K 1 reflects the relative stability of conformation A and is present in all performance descriptor expressions. K 2 reflects the affinity between the aptamer and its cognate ligand and is only present in the expression for EC 50 .
For all regulatory mechanisms, K 1 and K 2 reflect reversible conformational switching and ligand association, the core processes of riboswitch function. These processes are opposed by irreversible events that deplete the abundance of both conformations: mRNA degradation for translational repression and mRNA destabilization, and the riboswitch's decision to terminate or extend for transcriptional termination. The ratio between the rate constants for irreversible and reversible events is prevalent in all expressions for the performance descriptors (Text S2). This ratio is encapsulated in two reduced parameters c 1 and c 2 : Here, k iA and k iB represent the irreversible rate constants for conformation A or B, respectively, for translational repression (k dM ), transcriptional termination (k M ), and mRNA degradation (k dMA and k dMB ). Three operating regimes can be defined based on the ratio of reversible and irreversible rate constants within c 1 and c 2 . The first regime occurs when both reversible rate constants dominate (c 1 , c 2 converge to one), the second begins when either of the reversible rate constants is balanced with the associated irreversible rate constant (either c 1 or c 2 is less than one), and the third begins when the irreversible rate constant k iA dominates over k 1 (c 1 converge to zero). Each regime is generally determined by the competition between reversible and irreversible rate constants. We next evaluated the tuning properties of each regime for all regulatory mechanisms.

Riboswitches display similar tuning trends in the thermodynamically-driven regime
For dominating reversible rate constants (c 1 = c 2 = 1), a riboswitch molecule can sample both conformations and bind and unbind ligand many times before the irreversible event occurs. We define this regime as 'thermodynamically-driven' in accord with previous uses of this term in the study of natural riboswitches [21,24,25], since energetics dictate the prevalence of each conformation.
In the thermodynamically-driven regime, riboswitch function is captured for the three regulatory mechanisms by a general molecular description ( Figure 2A). The associated response curve is captured by a single equation that includes the partitioning constant (K 1 ), the aptamer association constant (K 2 ), mRNA degradation rate constants (k dMA , k dMB ), and representative regulatory activities of conformations A (K A ) and B (K B ): The values of K A and K B depend on the selected regulatory mechanism and are provided in Text S2. Parameter variation has a unique effect on the response curve for both ON ( Figure 2B and 2C) and OFF behaviors ( Figure S2). Increasing K 1 stabilizes conformation A, resulting in more riboswitch molecules adopting this conformation. Since conformation A has lower regulatory activity for ON behavior (K A ,K B ), basal levels decrease. Concomitantly, EC 50 increases as higher ligand concentrations are required to offset the decreased abundance of conformation B. Increasing K 2 reduces EC 50 as expected when aptamer affinity is modulated. However, K 2 has no effect on dynamic range and ligand-saturating levels since we assumed sufficient ligand can be added to saturate the response curve. Previous mutational studies of two synthetic riboswitches [16,17] support these model predictions. However, these studies Rates for endogenous mRNA degradation (norm) and ribozyme cleavage (rib) are separately described. a Reflects the time to reach the termination stem following transcription initiation and is dependent on the rate of polymerase extension, pausing, the length of the transcribed sequence, and nucleotide concentration. b Includes mRNA or protein degradation and dilution due to cell division. examined trans-acting mechanisms, calling into question whether model insights apply to cis-acting mechanisms. Finally, rate constants distinct from the core processes of riboswitch function such as transcription initiation (k f ) and protein decay (k dP ) affect both basal levels and dynamic range by modulating the steadystate mRNA and protein levels.
Stabilizing conformation A (increasing K 1 ) improves the dynamic range to an upper limit set by the regulatory activities (K A , K B ) and mRNA degradation rate constants (k dMA , k dMB ) associated with each conformation ( Figure 2D). While all four parameters affect the dynamic range, k dMA and k dMB also impact the dependence of the dynamic range on conformational partitioning ( Figure S3). This latter effect results from the dominant influence of the larger mRNA degradation rate on steady-state mRNA levels, which can be countered by biasing partitioning toward the more stable conformation. Therefore, when conformation A degrades faster (higher k dMA , ON behavior), less partitioning toward conformation A (lower K 1 ) is required to separate basal and ligand-saturating levels, whereas more partitioning toward conformation A (higher K 1 ) is required when conformation B degrades faster (higher k dMB , OFF behavior). As a result, thermodynamically-driven riboswitches functioning through mRNA destabilization require more (OFF behavior) or less (ON behavior) partitioning toward conformation A to achieve a larger dynamic range. In contrast, riboswitches functioning through translational repression and transcriptional termination display similar trends in dynamic range as a function of conformational partitioning for ON and OFF behaviors as the degradation rate constants are the same for each conformation.
EC 50 is also dependent on the value of K 1 according to the following relationship: All riboswitches can reversibly switch between conformations A and B that display different regulatory activities or different rates of degradation. Conformation B contains a formed aptamer that can reversibly bind ligand. Models assume negligible ligand binding during transcription. Green arrows designate mRNA synthesis with biased transcriptional folding, red arrows designate species degradation, and blue arrows designate translation that is proportional to mRNA levels. Under transcriptional termination (C), riboswitches effectively choose between termination to form a truncated product (T) and extension to form the full-length mRNA (M). To ensure both conformations make the decision with the same frequency, we designated a rate constant k M equal to the sum of the rate constants for extension and termination for either conformation (k M = k MA +k TA = k MB +k TB ). doi:10.1371/journal.pcbi.1000363.g001 EC 50 approaches the lower limit set by the aptamer dissociation constant when riboswitches principally adopt conformation B (low K 1 ). Therefore, although stabilizing conformation A (increasing K 1 ) can improve the dynamic range, excessive stabilization can be detrimental due to the increase in EC 50 . As a result, tuning strategies based on increasing K 1 require higher ligand concentrations to access the improved dynamic range. The ratio of the mRNA degradation rate constants in the expression for EC 50 offsets the modified dependence of the dynamic range on K 1 for riboswitches functioning through mRNA destabilization. Therefore, riboswitches functioning through any of the regulatory mechanisms exhibit the same trade-off between EC 50 and dynamic range.

Riboswitches display expanded tunability with reduced performance in the kinetically-driven regime
The second regime begins when either of the irreversible rate constants balances the associated reversible rate constant (either c 1 or c 2 is between zero and one). We call this regime the 'kinetically-driven' regime in accord with uses of this term in the study of natural riboswitches [21,22], where performance is driven by kinetics over energetics. In this regime, riboswitch molecules have fewer opportunities to sample both conformations and bind and release ligand before the irreversible event occurs, where the number of opportunities is governed by the competition between reversible and irreversible rate constants. Since c 1 is coupled to K 1 and k fA while c 2 is coupled to K 2 , both c 1 and c 2 are anticipated to have a significant impact on the response curve and impart several tuning properties distinct to this regime. We initially use riboswitches functioning through transcriptional termination to highlight two of these tuning properties ( Figure 3A-C).
First, irreversible rate constants modulate all performance descriptors, often at a cost to riboswitch performance. As the rate constant for terminator stem formation (k M ) increases, riboswitch molecules become trapped in a given conformation after transcriptional folding or conformational switching as reflected in c 1 . This effect reduces the dynamic range ( Figure 3A) and shifts basal and ligand-saturating levels according to the extent of transcriptional folding ( Figure 3B). The reduction in dynamic range can be offset by increasing the overall mRNA and protein abundance through modulation of the rates of transcription (k f ), translation (k P ), and mRNA (k dM ) and protein (k dP ) degradation. However, such changes also increase basal levels. c 1 and c 2 both impact EC 50 according to the following relationship: Since c 1 and c 2 reflect the ratios of k M /k 1 and k M /k 2 9, respectively, the relationship between EC 50 and k M depends on both conformational switching (k 1 ) and ligand release (k 2 9) ( Figure 3C). Increasing k M , which reduces both c 1 and c 2 , can increase or decrease EC 50 based on the opposing contributions of c 1 and c 2 . Reducing c 1 decreases EC 50 by restricting the time available to switch between conformations, while reducing c 2 increases EC 50 by decreasing the half-life of the ligand-aptamer complex (BL). Therefore, the relative values of k 1 and k 2 9 must be known to predict the effect of modulating the irreversible rate constant on EC 50 . Second, biased transcriptional folding can modulate the relationship between irreversible rate constants and the dynamic Figure 2. Thermodynamically-driven riboswitches display similar tuning properties. (A) Simplified molecular description that captures the three riboswitch regulatory mechanisms in the thermodynamically-driven regime. K 1 is the conformational partitioning constant (k 1 9/k 1 ) and K 2 is the aptamer association constant (k 2 /k 2 9). Coloring is the same as in Figure 1B-D. K A and K B reflect the regulatory activity of conformations A and B, respectively, and are specific to each regulatory mechanism: translational repression (k PA , k PB ), transcriptional termination (k P ?k MA /k M , k P ?k MB /k M ), and mRNA destabilization (k P , k P ). (B) K 1 affects both basal levels and EC 50 . (C) K 2 only affects EC 50 . Parameter values for red response curves: K 1 = 10; K 2 = 1/mM; K A = 10 23 /s; K B = 10 22 /s; k f = 10 211 M/s; k dP = 10 23 /s; k dMA = k dMB = 10 23 /s. (D) Biased conformational partitioning toward conformation B maximizes the dynamic range (g) at the cost of increased EC 50 . doi:10.1371/journal.pcbi.1000363.g002 range ( Figure 3A). When transcriptional folding is biased toward conformation A (k fB /k f = 0), riboswitch molecules must have sufficient time to switch between conformations to maintain activity. Therefore, the dynamic range declines as k M approaches and surpasses k 1 . In contrast, when transcriptional folding is biased toward conformation B (k fB /k f = 1), riboswitch molecules must switch to conformation A before the irreversible event occurs. In this case, k M must exceed the sum 2k 1 +k 1 9 to reduce the dynamic range. As a result, biasing transcriptional folding toward conformation B in the kinetically-driven regime increases the dynamic range.
A third tuning property is associated with riboswitches functioning through mRNA degradation. The rate constants for mRNA degradation (k dMA , k dMB ) are responsible for both the irreversible event and the steady-state basal and ligand-saturating levels, resulting in complex tuning properties ( Figure S4). Increasing either k dMA (ON behavior) or k dMB (OFF behavior) initially improves the dynamic range by separating the steady-state basal and ligand-saturating levels. However, the impact of larger values of k dMA and k dMB depends on riboswitch behavior. For ON behavior, if a riboswitch predominantly folds into conformation A during transcription (k fB /k f = 0), then values of k dMA in excess of k 1 diminish the dynamic range as conformation A is degraded before it can switch conformations. However, if a riboswitch predominantly folds into conformation B during transcription (k fB /k f = 1), then the dynamic range plateaus as each molecule either binds ligand or irreversibly switches to conformation A. In contrast, transcriptional folding has a negligible impact on the relationship between the dynamic range and k dMB for OFF behavior, since molecules that adopt conformation A will switch to conformation B before undergoing degradation. Furthermore, as observed for the thermodynamically-driven regime, more partitioning toward conformation A (higher K 1 ) is required to counteract the influence of k dMB on mRNA steady-state levels. Increasing k dMB eventually dominates basal levels when partitioning is maintained, leading to a loss in the dynamic range ( Figure S3).
As such, a tailored design approach is required to account for the difference between ON and OFF behavior for kineticallydriven riboswitches functioning through mRNA destabilization. Transcriptional folding is a key tuning parameter for ON behavior and should be the predominant focus before tuning the degradation rate of conformation A. In contrast, transcriptional folding can be largely ignored for OFF behavior, and the degradation rate of conformation B must be properly tuned to optimize the dynamic range.

Rescuing riboswitch performance in the non-functional regime
Higher irreversible rate constants require increased ligand concentrations to achieve a diminishing change in protein levels. We define the 'non-functional' regime as one in which riboswitches are effectively trapped in the conformation formed during transcriptional folding (c 1 = 0). In this regime, ligand has a negligible effect on performance. The fast time scales of terminator stem formation and mRNA cleavage may drive riboswitches functioning through these regulatory mechanisms into this regime. Our analysis of the kinetically-driven regime revealed that performance can be preserved by biasing transcriptional folding toward conformation B and ensuring that k 1 9 exceeds the irreversible rate constant k iA . However, these approaches do not alleviate the increased EC 50 caused by the reduced half-life of the ligand-aptamer complex (BL) when c 2 approaches 0. As a potential solution, studies of natural riboswitches have suggested that ligand binding during transcription can preserve EC 50 [21,22]. Therefore, we examined the effect of ligand binding during transcription under the assumption that conformation B is solely available (k fB /k f = 1) prior to polymerase extension (k E ) through the gene regulatory element responsible for the irreversible event.
We examined ligand binding during transcription for riboswitches functioning through transcriptional termination ( Figure 4A). We assumed that terminator stem formation (k M ) occurs much faster than ligand release (k 2 9) and the progression from conformation A to B (k 1 ) to limit consideration to nonfunctional riboswitches. Under these assumptions, the dynamic range is dependent on the ratio of read-through efficiencies for conformations A (k MA /k M ) and B (k MB /k M ), the progression from conformation B to A (k 1 9), and the rate of terminator stem formation (k M ). The dynamic range is maximized when conformational progression occurs much faster than terminator stem formation ( Figure 4B) as predicted from our analysis of the kinetically-driven regime ( Figure 3A). An in vitro study of the ribD FMN riboswitch operating through transcriptional termination yielded a reduced dynamic range when removing the polymerase pause site in the terminator sequence, increasing the nucleotide concentration, and withholding the NusA protein responsible for increasing polymerase residence time at pause sites [22]. These manipulations are expected to reflect an increase in k M and thus support our model predictions. If increasing k 1 9 above k M maximizes the dynamic range, riboswitches operating in this regime are expected to display strong stabilization of conformation A reflecting rapid progression from conformation B. In support of this claim, full-length riboswitches operating under transcriptional termination strongly prefer the aptamer-disrupted conformation and exhibit negligible ligand binding affinity [22,25,26].
EC 50 tuning properties are strikingly different for riboswitches in which ligand binding during transcription allows for improved performance than those for thermodynamically-driven riboswitches. EC 50 depends on model parameters in Figure 4A according to the following relationship: Both ligand release (k 2 9) and the time necessary to transcribe the sequences required for the formation of conformation A (k E ) have a significant impact on the value of EC 50 ( Figure 4C). Interestingly, tuning of k E decouples EC 50 and basal levels such that EC 50 can equal the aptamer dissociation constant (k 2 /k 2 9) without impacting the dynamic range. In contrast, the EC 50 of a thermodynamically-driven riboswitch approaches the aptamer dissociation constant as conformation B is stabilized, resulting in a concomitant decrease in the dynamic range ( Figure 2D). A previous theoretical study of the pbuE adenine riboswitch using experimentally measured kinetic rates also concluded that modulating polymerase extension time can tune EC 50 when the extension time is not significantly slower than ligand release [21].

Restricting the ligand concentration upper limit alters observed tuning properties
In our analyses thus far, we assumed that the maximum ligand concentration always saturates the response curve. However, studies of synthetic riboswitches have demonstrated that the response curve may not be saturated by the accessible upper limit in ligand concentration ( Figure 5A) due to various system properties including aptamer affinity, ligand solubility, permeability of the ligand across the cell membrane, and cytotoxicity of the ligand [16,17,19,[27][28][29]. Furthermore, natural riboswitches may regularly function in response to physiologically-relevant changes in metabolite concentrations that are much smaller than the ,1000-fold range necessary to access the full riboswitch response curve. To assess the effect of establishing an upper limit to the ligand concentration, we evaluated the response curve descriptors for a maximum ligand concentration of L'. An apparent EC 50 (EC 50 APP ) was calculated according to protein levels at L = 0 and L'.
Restricting L' alters the dependence of the dynamic range ( Figure 5B) and the apparent EC 50 ( Figure 5C) on model parameters as illustrated for riboswitches operating in the thermodynamically-driven regime. L' acts as a system restriction that prevents access to the full response curve such that increasing K 1 shifts the actual EC 50 beyond L', thereby reducing the maximum dynamic range that can be achieved. This behavior was recently observed for a trans-acting synthetic riboswitch operating under a limited ligand concentration range [17], supporting model predictions. Reflecting this behavior, the apparent EC 50 has the following dependence: where the apparent EC 50 converges to L'/2 as expected for a linear response when L' is below the EC 50 for an unbounded ligand concentration range ( Figure 5D). Our modeling results demonstrate that restricting the ligand concentration upper limit reduces riboswitch performance and establishes a unique relationship between dynamic range and conformational partitioning. In addition to serving as a design constraint for synthetic riboswitches, natural riboswitches may inherently operate under defined limits in ligand concentration. Future experiments may focus on measuring the physiologically-relevant metabolite concentration range experienced by natural riboswitches to examine what section of the response curve is utilized.

Application of tuning strategies to a synthetic riboswitch supports model predictions
To begin evaluating how the predicted tuning trends apply to both natural and synthetic riboswitches, we physically manipulated a recently-described synthetic riboswitch functioning through translational repression that up-regulates gene expression (ON behavior) in the presence of theophylline [18] (Figure 6A). Under the naming convention from Figure 1B, conformation A comprises a base-paired structure between the aptamer and RBS, while conformation B includes a formed aptamer and a single-stranded RBS. This riboswitch was selected because it closely resembles natural riboswitches functioning through translational repression, experimental data suggest that this riboswitch operates in the thermodynamically-driven regime [18], the ligand concentration upper limit does not saturate the response curve [28], and the demonstration that different sequences yield different response curves suggests riboswitch tuning [18]. A theophylline concentration of 1 mM was used as an upper limit, as exceeding this concentration inhibited cell growth. In studies performed by Lynch and coworkers, sequences associated with desirable response curves were identified by randomization of the RBS and screening for variants with low basal activity and a large activity increase in the presence of theophylline. Since the randomized sequence was located in a region responsible for conformational partitioning and translation, mutations most likely reflect simultaneous modulation of K A , K B , and K 1 . We therefore sought to introduce directed mutations to solely modulate individual model parameters and test model predictions for a thermodynamically-driven riboswitch with a ligand concentration upper limit that prevents response curve saturation. We examined two model predictions that could not be supported with available experimental data for cis-acting riboswitches: (1) solely modulating conformational partitioning (K 1 ) affects both EC 50 and basal levels ( Figure 2B), and (2) the dynamic range can be optimized by modulating K 1 when the ligand concentration upper limit does not saturate the response curve ( Figure 5B). We modulated K 1 by introducing systematic mutations into the aptamer stem while preserving the RBS sequence (m1-4; Figure 6A). Mutant sequences were ordered with increasing K 1 based on the energetic difference between conformations predicted by the RNA folding algorithm Mfold [17]. The mutations were not anticipated to significantly affect aptamer affinity (K 2 ) [30,31] or translational efficiency for either conformation (K A , K B ). Additional mutants were examined that are predicted to entirely favor either conformation A (mA) or conformation B (mB) to establish the regulatory activity of either conformation. Riboswitch performance was evaluated by measuring b-Galactosidase levels over a range of theophylline concentrations. The introduced mutations altered the response curve in agreement with model predictions (Figure 6B-D). Protein levels in the presence and absence of theophylline correlated with the relative stability of conformation A. Furthermore, complete stabilization of conformation A (mA) and conformation B (mB) established respective lower and upper limits for the observed expression levels. As predicted for a non-saturating value of L', an intermediate conformational partitioning value optimized the dynamic range to a value that was below the maximum dynamic range (g max = 15,600 MU) ( Figure 6B), and EC 50 approached 0.5 mM (L'/2) for increased stabilization of conformation A ( Figure 6C and 6D). Dynamic range optimization is clearly observed when evaluating the ratio of high and low protein levels, which is predicted to display the same qualitative tuning behavior ( Figure S5). The data agree with our model predictions for K 1 modulation in the thermodynamically-driven regime under conditions where the ligand concentration upper limit does not saturate the response curve, although we cannot rule out the possibility that stabilization of conformation A inadvertently drove the riboswitch into the kinetically-driven regime. The introduction of the aptamer sequence to the regulatory element decreased the regulatory activity of conformation B as observed when comparing protein levels for mB and a construct harboring only the RBS and aptamer basal stem (empty; Figure 6B). Our previous construction and characterization of a trans-acting synthetic riboswitch functioning through RNA interference [17] also showed submaximum dynamic range optimization when the ligand concentration was limiting and compromised activity of the regulatory element due to introduction of the aptamer element of the riboswitch. Thus, the results support the extension of our model predictions to synthetic riboswitches. In addition, our modeling results may have direct implications for the performance and tuning of natural riboswitches based on the similarity between the synthetic riboswitch examined here and natural riboswitches operating under translational repression.

Discussion
While our modeling efforts focused on translational repression, transcriptional termination, and mRNA destabilization, the predicted tuning trends generally apply to riboswitches utilizing other regulatory mechanisms. For example, riboswitches that function through alternative splicing mechanisms [32] can be modeled using the same approach applied to transcriptional termination-based riboswitches, where the splicing event occurs during a time interval in which riboswitch conformation determines the final exon composition. The identity and relative value of the irreversible rate constant for each conformation are important in determining the tuning properties associated with different regulatory mechanisms. The varying effects of irreversible rate constants on riboswitch performance are highlighted by the different tuning properties for riboswitches functioning through transcriptional termination and mRNA destabilization in the kinetically-driven regime.
The competition between reversible and irreversible rate constants establishes three operating regimes with distinct tuning properties. Therefore, measuring the reversible and irreversible rate constants is critical when predicting the impact of parameter modulation on the response curve. While well-established methods allow measurement of the rates of mRNA degradation and ligand binding and release, measuring the rates of RNA folding and conformational inter-conversion is currently an active area of research. New technologies are emerging that allow the measurement of kinetic folding rates: site-specific incorporation of aminopurines [25,33], single-molecule force experiments [34][35][36], and single-molecule fluorescence resonance energy transfer [37]. Studies of natural and synthetic riboswitches that apply these approaches may yield a comprehensive understanding of the relationship between riboswitch function and performance [38].
An alternative approach to measuring conformational switching relies on parameter predictions with RNA folding algorithms. Most algorithms calculate the free energy of individual conformations and can be used to estimate the value of K 1 for a riboswitch sequence [39,40]. Algorithms have also been developed that provide estimates of the rate constants for conformational switching (k 1 , k 1 9) [41]. By employing these algorithms, sequences can be rapidly screened in silico to identify riboswitches with tuned conformational partitioning according to model predictions. Mutations that impact other parameters, such as mutations to the RBS sequence that affect regulatory activity, can also be screened in silico to evaluate the impact on secondary structure and conformational partitioning. However, these algorithms are often inaccurate when predicting RNA folding in vivo, requiring modified approaches [17] or the development of more advanced algorithms [40].

Design principles for synthetic riboswitches
Synthetic riboswitches can be divided into two categories based on the intended application: inducible regulators and autonomous regulators. The applicable category depends on the identity and source of the detected ligand and requires distinct approaches to riboswitch design. We provide the following design principles assembled from our modeling results to guide the design of synthetic riboswitches as inducible or autonomous regulatory systems.
The desired properties of inducible regulatory systems include large dynamic ranges, low basal expression levels, and titratable control over expression levels. Selecting an effective regulatory mechanism is critical since numerous factors reduce the dynamic range, such as conformational partitioning, dominating irreversible rates, upper limits to ligand concentration, and reduced gene regulatory efficiencies from the incorporation of other riboswitch elements [17,19]. A design that is biased toward forming the disrupted-aptamer conformation (high K 1 ) will generally increase the dynamic range, although such strategies require higher ligand concentrations to modulate protein levels. The rates of events separate from core riboswitch processes, such as transcription, translation, and protein decay, can be modulated to increase the dynamic range difference at the expense of increased basal levels.
The selected regulatory mechanism will likely dictate the values of the irreversible rate constants and thus the operating regime. In support of this, studies on natural riboswitches have suggested a consistent pairing between translational repression and the thermodynamically-driven regime [25] and transcriptional termination and the non-functional regime with ligand binding during transcription [21,22,25,26,33]. Therefore, the design of inducible regulatory systems may rely on the tuning properties associated with each regime. While thermodynamically-driven riboswitches generally provide for the largest dynamic range, kinetically-driven and non-functional riboswitches can be designed to perform similarly using insights from our modeling efforts. In general, placing the aptamer toward the 59 end of the riboswitch sequence will preserve the dynamic range by biasing transcriptional folding toward conformation B. The exception is OFF-behaving riboswitches acting through mRNA destabilization, which are insensitive to biased transcriptional folding ( Figure S4). In addition, introducing pause sites and ensuring rapid conformational switching from the aptamer-formed conformation (k 1 9) will allow kinetically-driven and non-functional riboswitches to exploit ligand binding during transcription, thereby decreasing the amount of ligand required to induce gene expression.
In many practical applications, system restrictions will limit the accessible range of the response curve ( Figure 7A-C). Such limitations need to be addressed through parameter tuning in order to access the appropriate section of the response curve. For most biological systems, a predominant restriction is a limit to the maximum ligand concentration. In situations where the maximum ligand concentration does not saturate the response curve, designs for thermodynamically-driven riboswitches should be based on intermediate conformational partitioning values (K 1 ) that achieve a suboptimal maximum dynamic range. An alternative strategy is the design of non-functional riboswitches that bind ligand during transcription, which can respond to ligand at lower concentrations without sacrificing the dynamic range.
Genes often exist in regulatory networks that dictate cellular phenotype such that complex relationships exist between the expression levels of individual genes and systems-level functions. To effectively regulate these genes with synthetic riboswitches, a variety of tuning strategies must be employed to tune the response curve to operate within system restrictions. The properties of the regulated gene, its integration into a biological network, and the ultimate systems-level functions must be considered. One can envision an application-specific regulatory niche that defines the acceptable ranges of basal and maximum-ligand expression levels for proper system performance ( Figure 7B). For example, the properties of a given system may require that the riboswitch be tuned to minimize basal expression over maximizing dynamic range, such as when the regulated enzyme exhibits high activity or cytotoxicity.
The engineering of synthetic riboswitches that act as autonomous regulatory systems presents an even greater design challenge. Here, the upper and lower ligand concentrations that the system fluctuates between establish the accessible section of the response curve such that the regulatory niche is further restricted ( Figure 7C). For example, riboswitches responsive to an endogenous central metabolite will likely be operating under a defined concentration range characteristic of the organism and the environment. In this case, the response curve must be tuned to place the desired expression levels at the limits of this defined concentration range by modulating the appropriate performance descriptors. Depending on system restrictions, proper tuning of riboswitches acting as autonomous control systems may require minimization of basal levels, operation across higher expression levels, or maximization of the change in expression levels.
Many parameters can potentially be modulated to tune the response curve. However, current practical considerations favor the modulation of a subset of these parameters in the laboratory. As one example, a given riboswitch may require a higher EC 50 value to meet the performance requirements. Aptamer affinity (K 2 ), conformational partitioning (K 1 ), and the irreversible rates associated with the gene regulatory mechanism can be modulated to increase EC 50 . However, rational modulation of aptamer affinity is restrictive since most mutations effectively abolish ligand binding, while the method and ease of modulating irreversible rates depend on the regulatory mechanism. Modulating conformational partitioning is an attractive approach since simple base-pairing interactions principally establish each conformation. However, conformational partitioning also impacts basal levels and the dynamic range, such that other parameters may need to be modulated to compensate for any undesired changes. Thus, the effective design of synthetic riboswitches requires knowledge of the relationship between riboswitch sequence and model parameters and may require the coordinated modulation of multiple parameters to meet application-specific performance requirements.
The relationship between riboswitch sequence and model parameters depends in part on the composition framework used in the riboswitch design. A synthetic riboswitch can be designed such that parameters map to individual domains [16,17,19] or multiple domains [18,42,43]. Each design strategy offers distinct advantages depending on whether rational design or random screening is used to select riboswitch sequences. Individual domain mapping strategies allow for insulated control over each parameter and domain swapping without requiring redesign of the riboswitch, thereby presenting significant advantages for rational design approaches. Multiple domain mapping strategies may be more desirable for random screening approaches, where assigning multiple parameters to a single sequence domain can reduce the number of randomized nucleotides required to sufficiently sample parameter space.

Evolutionary implications for tuning in natural riboswitches
Natural riboswitches primarily serve as key autonomous regulators of diverse metabolic processes [10]. Recent characterization of eleven known S-adenosylmethionine riboswitches in Figure 7. The accessibility of the riboswitch response curve depends on application category and associated system restrictions. Categories include an inducible regulatory system with (A) no ligand limitations or (B) a ligand concentration upper limit, and (C) an autonomous regulatory system with defined lower and upper limits for the ligand concentration. The accessible dynamic range (g) for each response curve depends on the system restrictions. The properties of other components in the network will dictate which riboswitch design best meets performance requirements. For example, under the autonomous regulatory system (C) the red curve may be more appropriate for the regulation of cytotoxic genes, the orange curve may be more appropriate for the regulation of enzymes with low catalytic activity, and the blue curve may be more appropriate for regulatory networks that require a large change in protein levels. doi:10.1371/journal.pcbi.1000363.g007 Bacillus subtilis demonstrated that these riboswitches exhibit a diverse range of values for basal expression levels, EC 50 , and dynamic range [44], suggesting that natural riboswitches are finely tuned to match their occupied regulatory niche. However, this study is the only one to date to characterize the response curves of multiple natural riboswitches responsive to the same ligand. Two questions emerge from these observations and our modeling results that underlie the biological utilization of natural riboswitches as dynamic regulators of metabolism: (1) how finely tuned are natural riboswitches to their regulatory niche, and (2) what sequence modifications are associated with response curve tuning?
Understanding the extent to which natural riboswitches are tuned to their regulatory niches will provide insights into riboswitch utilization and the underlying principles of genetic regulatory control. Similar to the tuning of synthetic riboswitches to match their intended regulatory niche, investigating the extent and biological relevance of natural riboswitch tuning requires knowledge of the functional properties of the regulated genes and their contribution to cellular fitness. Furthermore, the typical ligand concentration range encountered in the intracellular environment designates the operational section of the response curve, such that determining this range is critical to advancing our understanding of natural riboswitch tuning within regulatory niches.
The composition of a natural riboswitch dictates the relationship between its sequence and model parameters. One way to gain insights into this relationship is investigating sequence deviations between natural riboswitches in the same organism or different organisms that recognize the same ligand and employ the same regulatory mechanism. Using the response curve as a basis of comparison, these mutations may be neutral or shift the response curve in line with modulation of single or multiple parameters. Identifying which parameters are modulated will provide insights into how accessible each parameter is to random point mutations and how evolution effectively tunes the response curve through parameter modulation. Advances in our understanding of the biological utilization of natural riboswitches will enable researchers to better define regulatory niches in a biological system and more effectively design synthetic riboswitches to match these niches. Beyond riboswitch design and implementation, insights into the fine-tuning of natural regulatory components and networks will enable the construction of biological networks that reliably control systems-level functions.

Mathematical modeling
All modeling assumptions and methods are fully described in Text S2. Briefly, time-dependent differential equations were generated using mass action kinetics to describe each mechanistic step in the simplified molecular descriptions of riboswitch function for translational repression, transcriptional termination, and mRNA degradation. The resulting equations were simplified by assuming steady-state conditions. Relevant tuning properties were identified based on the impact of model parameters on the response curve descriptors, including dynamic range (g) defined as the difference between high and low protein levels, ligand concentration to induce a half-maximal response (EC 50 ), basal protein levels (P(L = 0)), and maximum-ligand protein levels (P(LRL' or ')).

Plasmid construction
pSAL8.3 served as the base plasmid for all experimental studies [18]. A theophylline-dependent synthetic riboswitch functioning through translational repression resides between the upstream P tac1 promoter and the downstream Tn10-b-Galactosidase fusion gene. Mutant sequences were cloned into the unique KpnI and HindIII restriction sites located directly upstream of the riboswitch and approximately 200 nucleotides into the fusion gene coding region. Primers harboring mutant sequences (Table S1) and a 59 KpnI site were used to PCR amplify the 59 untranslated region extending through the HindIII restriction site. The resulting PCR product was digested with KpnI/HindIII, ligated into pSAL8.3 digested with the same restriction enzymes, and transformed into Escherichia coli strain DH10B. Assembled plasmid constructs were verified by sequencing (Laragen). All molecular biology reagents and enzymes were obtained from New England Biolabs. b-Galactosidase activity assay b-Galactosidase assays were conducted using E. coli DH10B cells harboring the pSAL8.3 plasmid mutants based on modifications to previously described protocols [18,45]. Cells harboring each construct were grown overnight in Luria-Bertani (LB) broth supplemented with 50 mg/ml ampicillin. Overnight cultures were back-diluted into three separate wells containing 500 ml LB broth with 50 mg/ml ampicillin and the appropriate concentration of theophylline and grown at 37uC for 3 hrs with shaking at 210 RPM. Approximately 3 ml of the overnight culture was added to each well. Following the 3-hr incubation with shaking, optical density was recorded by transferring 175 ml into a 96-well microplate with a mClear bottom (Greiner) and measuring on a Safire fluorescence plate reader (Tecan). Cells were lysed by mixing 20 ml of culture with 80 ml permeabilization solution (100 mM Na 2 HPO 4 , 20 mM KCl, 2 mM MgSO 4 , 0.6 mg/ml CTAB, 0.4 mg/ml sodium deoxycholate, and 5.4 ml/ml bmercaptoethanol) and mixed at room temperature for approximately 10 min. In a fresh 96-well microplate, 25 ml of the lysed culture was mixed with 150 ml substrate solution (60 mM Na 2 HPO 4 , 40 mM NaH 2 PO 4 , 1 mg/ml ONPG, and 5.4 ml/ml b-mercaptoethanol). ONPG hydrolysis was stopped with the addition of 75 ml 1 M Na 2 CO 3 when a faint yellow color was observed. Absorbance at 420 nm was then measured on the fluorescence plate reader and protein levels were calculated in Miller Units (MU):

MU~1000
: where t is in minutes and absorbance values reflect the difference between each sample and blank media. The MU value of cells carrying a blank plasmid was also subtracted from each sample measurement. Figure S1 Detailed molecular description of riboswitch function. Molecular descriptions shown for riboswitches functioning through (A) translational repression, (B) transcriptional termination, and (C) mRNA destabilization. All riboswitches can reversibly switch between conformations A and B that display different regulatory activities or different rates of degradation. Conformation B contains a formed aptamer that can reversibly bind ligand (L). Transcription is represented as two discrete steps designated by the subscripts I and II, although the model can be altered to include more or less steps. Riboswitches at each step may switch between conformations or reversibly bind ligand depending on the transcribed sequence. The lag between steps is captured by the rate constant k E . Once the full riboswitch sequence is transcribed, the riboswitch is susceptible to the regulatory mechanism that controls protein (P) production. Green arrows designate mRNA synthesis with biased transcriptional folding, red arrows designate species degradation, and blue arrows designate translation that is proportional to mRNA levels. Under transcriptional termination (B), riboswitches effectively choose between termination to form a truncated product (T) and extension to form the full-length mRNA (M). To ensure both conformations make the decision with the same frequency, we designated a rate constant k M equal to the sum of the rate constants for extension and termination for either conformation (k M = k MA +k TA = k MB +k TB ). We assumed that the riboswitch sequence immediately prior to termination or extension cannot undergo degradation. Found at: doi:10.1371/journal.pcbi.1000363.s001 (3.60 MB TIF) Figure S2 Thermodynamically-driven riboswitches exhibiting OFF behavior display similar tuning properties to riboswitches exhibiting ON behavior. K 1 is the conformational partitioning constant (k 1 9/k 1 ) and K 2 is the aptamer association constant (k 2 / k 2 9). Biased transcriptional folding significantly affects riboswitches displaying ON behavior. Riboswitches displaying OFF behavior show a negligible dependence on transcriptional folding (C, inset) for the selected parameter values. Parameter values: k 1 = 5*10 23 ; k 1 9 = 2*10 21 ; k 2 = 10 6 /M*s; k 2 9 = 10 23 /s; k P = 10 23 /s; k f = 10 211 M/s; k dP = 10 23 /s; k dMA = 10 24 /s for OFF behavior; k dMB = 10 24 /s for ON behavior.

Supporting Information
Found at: doi:10.1371/journal.pcbi.1000363.s004 (1.78 MB TIF) Figure S5 The dynamic range difference and ratio exhibit qualitatively similar tuning properties. Model predictions for the dynamic range difference (g D , A) and ratio (g R , B) when subjected to a ligand concentration upper limit (L'). In the absence of a ligand concentration upper limit, the dynamic range converges on a maximum (g max ). The optimum value of the conformational partitioning constant (K 1 ) is higher for the dynamic range ratio as the ratio favors lower basal levels. Increasing the aptamer association constant (K 2 ) or L' improve the suboptimal dynamic range maximum. Parameter values for the red curves are identical to those reported in Figure 5, and notation is identical to that used in Figure 5B. (C) b-Galactosidase assay results from Figure 6B, where the dynamic range is calculated as the ratio of b-Galactosidase levels in the presence (filled circle) and absence (empty circle) of 1 mM theophylline. The positive control construct (empty) harbors only the RBS and aptamer basal stem. A slight increase in b-Galactosidase activity was observed in the presence of theophylline for the control construct. The experimental data follow the general trends predicted from the model, including the higher optimum K 1 value for the dynamic range ratio as compared to the dynamic range difference. b-Galactosidase levels are reported in Miller Units (MU). Data represent independent measurements of triplicate samples, where the standard error was below 5% of each mean value.