Progress in constructing biological networks will rely on the development of more advanced components that can be predictably modified to yield optimal system performance. We have engineered an RNA‐based platform, which we call an shRNA switch, that provides for integrated ligand control of RNA interference (RNAi) by modular coupling of an aptamer, competing strand, and small hairpin (sh)RNA stem into a single component that links ligand concentration and target gene expression levels. A combined experimental and mathematical modelling approach identified multiple tuning strategies and moves towards a predictable framework for the forward design of shRNA switches. The utility of our platform is highlighted by the demonstration of fine‐tuning, multi‐input control, and model‐guided design of shRNA switches with an optimized dynamic range. Thus, shRNA switches can serve as an advanced component for the construction of complex biological systems and offer a controlled means of activating RNAi in disease therapeutics.
The engineering of biological systems is critical to providing effective solutions to global challenges including sustainable energy and health and medicine (Keasling, 2008; Savage et al, 2008). Many challenges remain in the longer‐term goal of building biological systems with predictable and reliable behaviour akin to a computer or an automobile (Endy, 2005). One significant challenge is to expand the current collection of biological components that can be combined to construct larger‐scale systems in a way that is amenable to computational design. Introducing properties such as standardization, modularity, and functional programmability into component design will greatly expand the scale at which biological systems can be constructed.
Assembling components out of RNA offers a unique approach that exploits the rich functional properties that have been demonstrated by this biological substrate: specific and tight binding by aptamers to diverse biomolecules comparable with protein antibodies (Gopinath, 2007), simple and predictable base‐pairing interactions that principally dictate overall conformation and activity, and numerous gene regulatory capabilities that rely on RNA's critical placement in the process of gene expression. Researchers have integrated all three properties to construct RNA‐based information processing devices termed riboswitches that control gene expression on the basis of the presence and concentration of a molecular ligand inside the cell (Suess and Weigand, 2008). Several examples exist of engineered riboswitches that have been shown to be amenable to property tuning reflected by predictable shifts in the transfer function linking the concentration of molecular ligand and expression levels of the target genes (Suess et al, 2003; Bayer and Smolke, 2005; Lynch et al, 2007; Win and Smolke, 2007). However, the observed qualitative shifts have never been further investigated to develop computational tools that support the forward engineering of riboswitch function.
We engineered a novel riboswitch termed an shRNA switch that exerts an effect through RNA interference (RNAi), a natural silencing pathway in humans and other eukaryotes that has garnered recent attention as a revolutionary biological research tool and a promising therapeutic strategy (Kim and Rossi, 2008). A composition framework was developed on the basis of the integration of three modular domains, an aptamer, shRNA stem, and competing strand, that utilizes a strand displacement mechanism to establish two distinct conformations with differential processing and subsequent activation of the RNAi machinery (Figure 1A). Ligand binding to the formed aptamer in the inactive conformation reduces processing, resulting in an increase in target gene expression levels in the presence of ligand (Figure 1C and D).
To systematically evaluate transfer function tuning, we developed a simple mathematical model from the proposed mechanism. Modelling results suggested multiple tuning strategies to uniquely tune the transfer function (Figure 3A). Experimental results from modifications to each domain agreed with model predictions and demonstrated the limits of shRNA switch tuning. Furthermore, we examined modularity by independently swapping each domain and demonstrating that proper shRNA switch functionality was maintained. These tuning and domain swapping strategies were combined to yield fine‐tuned multi‐input control of gene expression, a useful capability for the construction of large‐scale biological systems.
Although the initial model served as a qualitative predictor of general sequence modifications, an extended model that provides a quantitative link between an shRNA switch sequence and its in vivo behaviour would facilitate rapid design towards various applications. To this end, we employed RNA secondary structure prediction algorithms (Mathews et al, 2004) to calculate the ΔG separating active and inactive conformations, and used this value in the model (Figure 6A). Upon determining that common methods of calculating ΔG were inaccurate for this system, we devised a novel method that provided a strong correlation between ΔG and switch activity in the absence of ligand (Figure 6C). An extended model that includes empirical constants estimated from the correlation was used to predict the values of ΔG that maximize the dynamic range, an important performance indicator. After testing a series of switches constructed using the deduced design strategies, we identified a switch sequence with the predicted ΔG value that gave a dynamic range approximately equal to the theoretical limit predicted by the model (Figure 6D and E). Our results show that a computational model can be employed with secondary structure prediction algorithms for the forward design of shRNA switches with desired properties.
To explore the generality of our modelling approach, we conducted similar analyses for shRNA switches targeting an endogenous gene. Our devised method again produced the strongest correlation between ΔG and switch activity in the absence of ligand that accurately predicted the relationship between ΔG and dynamic range. The empirical values from these experiments were notably different than the previous values, suggesting that sequence‐specific factors contribute to the correlation serving as the foundation of the extended model. However, these experiments support the generality of shRNA switch design and the utility of model‐guided forward design within our framework.
shRNA switches can be implemented for integrated ligand control of gene silencing through RNA interference in mammalian cells.
Our design framework allows for control by different ligands (through aptamer domain swapping strategies), targeting of heterologous or endogenous genes (through regulator domain swapping strategies), integration into regulatory networks for multi‐input control (through combinatorial implementation strategies), and programming of the resulting transfer function (through coupling composition and modeling frameworks).
RNA switch design strategies based on the integration of distinct domains that carry out separate functions support domain swapping and rational tuning of the transfer function.
By linking designed conformational changes within an RNA switch molecule to competitive hybridization events, secondary structure folding algorithms can be combined with mathematical models to develop computational tools that yield quantitative sequence‐to‐function relationships and provide frameworks for the forward design of RNA switch properties.
To maintain fitness under diverse conditions, biological systems must integrate multiple environmental cues (inputs) to determine appropriate phenotypic outcomes (outputs) over short and long time scales. This relationship, which can be interpreted as an input–output function or transfer function, is specified by the behaviour of the individual system components and their network interactions. The complexity of natural biological systems, reflected by the sheer number of associated components and network interactions, can appear intractable to scientists and engineers seeking to understand and reliably construct biological systems.
Synthetic biological systems that perform information processing operations with specified transfer functions can be constructed through the design of either individual complex components encoding multiple integrated functionalities or simpler components assembled into networks with emergent properties. For example, the ultrasensitive switch behaviour of the mitogen‐activated protein kinase signalling cascade (Huang and Ferrell, 1996) has been replicated with one component (Dueber et al, 2007) or a network of components (Hooshangi et al, 2005). The majority of previously engineered biological systems have employed design strategies focused on the assembly of simpler components into networks (Elowitz and Leibler, 2000; Gardner et al, 2000). However, as biological engineers move towards constructing large‐scale systems with more advanced behaviours, integration of complex components into network design will be critical, especially as current network design strategies do not effectively scale with system complexity (Croft et al, 2003). Furthermore, engineering of complex components and their integration into networks will facilitate the construction of advanced systems with a reduced number of constituent parts. Such design strategies will result in a lower energetic load on the cell and will comply with size limitations associated with packaging synthetic systems for delivery applications (Flotte, 2000; Grieger and Samulski, 2005).
RNA is a rich biological substrate for engineering complex components (Isaacs et al, 2006), where scalable molecular information processing systems have been built in vitro from nucleic acid components (Stojanovic and Stefanovic, 2003; Seelig et al, 2006). Natural biological systems exhibit widespread utilization of regulatory RNAs in larger networks (Shalgi et al, 2007) and integrated functionalities by riboswitches (Sudarsan et al, 2003; Grundy and Henkin, 2006). The latter is an example of a complex RNA component that converts the intracellular concentration of a molecular signal into levels of a target protein. Building on these natural examples, researchers have integrated synthetic regulatory RNAs into larger networks (Deans et al, 2007; Rinaudo et al, 2007) and developed synthetic riboswitch elements (Isaacs et al, 2006; Suess and Weigand, 2008) to yield desired transfer functions. Synthetic riboswitches have been developed through the coupling of regulatory RNA elements to aptamers, sensory elements that bind specific ligands, to achieve in vivo ligand control of transcription (Buskirk et al, 2004), RNA stability (Win and Smolke, 2007), translation (Grate and Wilson, 2001; Suess et al, 2003; Bayer and Smolke, 2005; Lynch et al, 2007; Ogawa and Maeda, 2008; Wieland and Hartig, 2008), splicing (Thompson et al, 2002; Weigand and Suess, 2007), and RNA interference (RNAi) (An et al, 2006). As aptamers can be generated de novo to potentially any molecule (Lee et al, 2004; Gopinath, 2007) and regulatory RNAs can be designed to target practically any gene of interest, the potential exists to construct complex regulatory RNAs that respond to molecular inputs displaying low cytotoxicity, tissue targeting, or endogenous production to impose desired phenotypic outcomes. In addition, the ability to alter or tune the component transfer function, which allows manipulation of component performance towards optimal system performance, has been demonstrated for in vivo (Suess et al, 2003; Isaacs et al, 2004; Bayer and Smolke, 2005; Win and Smolke, 2007) and in vitro systems (Hall et al, 2007) by modifying RNA folding energetics. However, broader implementation of these synthetic regulatory components has been limited, as most examples do not support domain swapping of different sensory and regulatory elements. Furthermore, the design of such regulatory systems has lacked predictive tools for the translation of sequence information into component transfer functions to enable in silico optimization of system behaviour before construction.
We have developed a framework for the construction of shRNA switches that mediate ligand control of RNAi across diverse mammalian cell types. Our platform utilizes a strand displacement strategy, where the functions of ligand binding, RNAi activation (Kim and Rossi, 2008), and translation of the binding interaction into reduced processing by the RNAi machinery are isolated to individual domains, which increases the generality and ease of successful domain swapping and subsequent broad application. In addition, we systematically investigated tunability of the shRNA switch transfer function through a combined experimental and mathematical modelling approach that resulted in the identification of five tuning strategies. Standard RNA folding algorithms (Mathews et al, 2004) were used to establish a quantitative sequence‐to‐function relationship. Our efforts highlight the current limitations of these broadly used algorithms for the design of RNAs that function in vivo and offer a framework for optimizing shRNA switch behaviour in silico. By demonstrating combinatorial tuning strategies, multi‐input control, and model‐guided forward design of shRNA switches with an optimized dynamic range within a specified context, we show that shRNA switches extend the utility of RNAi as a regulatory tool and are valuable components for the construction of complex biological systems.
Design and characterization of a modular shRNA switch platform
We engineered a complex RNAi substrate that encodes a ligand‐controlled gene regulatory function by replacing the loop of a small hairpin (sh)RNA with two domains: an aptamer and a competing strand (Figure 1A). The shRNA switch molecule is designed to adopt distinct ‘active’ and ‘inactive’ conformations due to complementarity between the competing strand and the shRNA stem, similar to previously engineered RNA regulatory systems (Bayer and Smolke, 2005; Lynch et al, 2007; Win and Smolke, 2007). In the active conformation, irreversible processing by the RNAi machinery of the formed shRNA stem results in small interfering (si)RNA production and subsequent RNAi‐mediated silencing of the target gene. In the inactive conformation, base‐pairing by the competing strand disrupts the shRNA stem, which is predicted to inhibit processing by the RNAi machinery (Zeng and Cullen, 2004; Macrae et al, 2006). This base‐pairing coincides with formation of the aptamer domain, such that ligand binding stabilizes the inactive conformation and indirectly reduces siRNA production, thereby linking intracellular ligand concentration to target protein levels through a component transfer function. To decrease the activation energy separating the two conformations, we removed two nucleotides in the passenger strand, thereby mimicking the bulge from the microRNA (miRNA) mir‐30a (Griffiths‐Jones, 2004; Griffiths‐Jones et al, 2006).
The three domains that comprise an shRNA switch perform distinct functions: the shRNA stem encodes the guide strand that activates RNAi‐mediated silencing of the target gene, the aptamer detects the molecular input concentration through a ligand‐binding interaction, and the competing strand translates the binding interaction into a decrease in regulatory activity by affecting processing by the RNAi machinery. On the basis of the action of the competing strand that is complementary to the shRNA stem, the sequences of the shRNA stem and aptamer domains are independent of one another. Therefore, the shRNA stem and aptamer domains can be independently modified without altering the functionality of the opposing domain or requiring sequence reassignment.
We designed an initial shRNA switch (S1) to target EGFP and respond to theophylline by incorporating an EGFP‐targeting guide strand and the theophylline aptamer (Zimmermann et al, 2000) into our switch platform (Figure 1A). We used in‐line probing (Soukup and Breaker, 1999) to assess the structural characteristics of a T7‐transcribed variant similar to S1 (S4t; Figure 1B). In‐line probing provides information on structural changes within the molecule as a result of theophylline binding from the ligand dependence of spontaneous RNA cleavage. Theophylline‐dependent changes in cleavage rates at individual nucleotides were observed in the aptamer domain, competing strand, and the downstream shRNA stem sequence. The results suggest that theophylline binding promotes structural changes in the shRNA switch as expected for dynamic RNAs undergoing ligand‐dependent conformational switching. The apparent dissociation constant (KD) of ∼5 μM, which was determined by quantifying the cleavage products at two positions (Supplementary Figure S1), is an order of magnitude larger than that of the aptamer alone (KD≈0.29 μM) (Zimmermann et al, 2000). The observed increase in KD is in agreement with our proposed mechanism (see below), where only the inactive conformation provides a formed aptamer that can bind ligand. As shRNA switches can occupy both conformations, the apparent affinity will be lower because ligand can only bind the inactive conformation that is transiently present in a fraction of the shRNA switch population.
The functionality of shRNA switches was assessed in mammalian cell culture. We transiently transfected plasmids harbouring S1 and various switch controls transcribed from a U6 promoter into HEK293T cells stably expressing EGFP (Abbas‐Terki et al, 2002). Flow cytometry analysis revealed that S1 elicits intermediate knockdown of EGFP as compared to the original shRNA targeting EGFP (sh) and a scrambled shRNA (neg) (Figure 1C and D), where the observed silencing by S1 can be attributed to activation of RNAi on the basis of antisense inhibition of guide strand activity (Hutvagner et al, 2004; Meister et al, 2004) (Supplementary Figure S2). In the presence of theophylline, GFP levels increased in a dose‐dependent manner for S1 but not for the control shRNAs. The effective concentration to achieve 50% activity (EC50) for S1 of ∼300 μM was much larger than the KD of 5 μM measured in vitro, which can be primarily attributed to a concentration drop in theophylline across the cellular membrane (Koch, 1956) (J Liang, J Michener, C Smolke, unpublished data, 2008). Mutating the aptamer core of S1 (S1′) greatly reduced the observed theophylline dependence without perturbing basal expression levels. We attribute the minor increase in GFP levels at high theophylline concentrations for S1′ to pleiotropic effects of the ligand (An et al, 2006) and potentially to reduced binding capability of the mutated aptamer. Taken together, S1 links theophylline concentration to GFP levels in vivo through a relationship described by a component transfer function (Figure 1D). We obtained qualitatively similar results when shRNA switches targeting EGFP were transiently transfected into other cell lines (Supplementary Figure S3), suggesting that shRNA switches can be broadly applied in different cell lines and types.
Mathematical modelling offers tuning parameters to predictively modulate component transfer functions
Previous switch platforms utilizing strand displacement strategies have demonstrated tuning on the basis of aptamer swapping and modulation of folding energetics (Bayer and Smolke, 2005; Win and Smolke, 2007). We systematically evaluated the tuning capabilities of shRNA switches with the aid of a mathematical model relating ligand concentration and target gene expression levels. Standard model parameters were incorporated to represent each chemical step from our proposed mechanism (Supplementary text 1). We assumed that the two adopted conformations are at thermodynamic equilibrium, that ligand only binds the inactive conformation, and that the active conformation is solely processed to an siRNA with a reduced efficiency as compared to the original shRNA. These assumptions yield the following relationship between relative expression levels of the target gene (f; output) and exogenous ligand concentration (L; input):
where e is the processing efficiency, fshRNA is the relative knockdown achieved by the original shRNA (sh), Kcomp is the partitioning constant between active and inactive conformations (Kcomp=[inactive]/[active]), Kapt is the association constant for binding between ligand and the formed aptamer, and h is the Hill coefficient to account for non‐linearity between siRNA concentration and target expression levels. Although mathematical models have been developed for RNAi (Raab and Stephanopoulos, 2004; Bartlett and Davis, 2006; Malphettes and Fussenegger, 2006), our approach utilizes a minimal parameter set that is experimentally tractable, fully represents RNAi in the context of shRNA switches, and captures the steady‐state behaviour of our system (Supplementary Figure S4). For one shRNA stem sequence and input ligand (fixed fshRNA, h), our model provides three tuning parameters that can be varied to tune the component transfer function: Kcomp, Kapt, and e (Figure 2A–C). Varying Kcomp results in a concomitant and opposing variation in EC50 and basal expression levels, which are independently tuned by Kapt and e, respectively. In addition, as Kcomp approaches zero, basal expression levels approach a lower limit that is dependent on the value of e and is higher than that of the original shRNA (Figure 2D). As each tuning parameter represents individual steps in the proposed mechanism, we examined how modifying the sequence in each domain, specifically the competing strand and aptamer domains, corresponds to parameter variation to identify unique tuning strategies (Figure 3A).
Competing strand tuning strategies enable predictive alteration of the component transfer function
Modifying competing strand base‐pairing interactions is anticipated to reflect changes in Kcomp, as this parameter represents the thermodynamic partitioning between active and inactive conformations. We developed competing strand‐tuning strategies to target modifications to three regions within the competing strand domain: the length of the competing strand on the 3′‐end (Figure 3B and C) or the 5′‐end (Figure 3D and E), or the base‐pairing complementarity (Figure 3F and G). We introduced iterative nucleotide changes under each competing strand tuning strategy and generated component transfer functions as before. Regardless of the selected strategy, each nucleotide change resulted in a shift in the response curve in line with the model prediction for variation in Kcomp. The results suggest that decreasing the extent of base‐pairing interactions between the competing strand and the shRNA stem decreases the stability of, or bias towards, the inactive conformation (lower Kcomp), resulting in lower basal expression levels and a higher EC50. The trend towards higher EC50 is consistent with the order‐of‐magnitude difference between the apparent KD of S4t observed in the in‐line probing experiment and that reported for the aptamer alone (Figure 1B). Thus, sequence modifications to the competing strand that affect the extent of base‐pairing solely map to variation of Kcomp.
Aptamer tuning strategies enable predictive alteration of the component transfer function
Although ligand binding to the formed aptamer directly relates to aptamer affinity, represented by Kapt, sequence changes in the aptamer domain may affect other parameters. To evaluate how sequence modification of the aptamer domain corresponds to parameter variation, we tested two theophylline aptamer variants (S11 and S12) with dissimilar KD values (Zimmermann et al, 2000) and the mutated aptamer (S1′) (Figure 4A and B). Mutating the aptamer core (S1′) without perturbing shRNA switch secondary structure or sequence length resulted in a shift in EC50, whereas decreasing aptamer affinity by decreasing the aptamer stem length (S11 and S12) resulted in a shift in both EC50 and basal expression levels. The shifts in EC50 for S11 and S12 matched the relative KD measured in vitro for the aptamer variants alone (Zimmermann et al, 2000), suggesting that modulating aptamer affinity is reflected by variation in Kapt. However, Kapt affects only EC50, suggesting that either Kcomp or e varies with aptamer size. As the competing strand sequence is preserved for S1, S1′, S11, and S12, we hypothesized that the shift in basal expression levels independent of Kapt (most obvious in comparing the transfer functions of S1 and S11) is solely attributed to the third tuning parameter e (Figure 2C).
To evaluate the relationship between aptamer size and the tuning parameter e, we replaced the theophylline aptamer with the smaller xanthine aptamer (Kiga et al, 1998) or the larger tetracycline aptamer (Berens et al, 2001). As variation of e and Kcomp both affect basal expression levels, sole evaluation of e requires estimation of the lower limit of basal expression levels for vanishingly small values of Kcomp (Figure 2D). To this end, we constructed at least one shRNA switch with each aptamer that strongly prefers the active conformation (low Kcomp) and measured GFP basal expression levels of cells transfected with these constructs (Figure 4C). Assay results indicated that aptamer size strongly correlated with the lower limit of basal expression levels. The results suggest that the tuning parameter e, which is predicted to have a significant effect on the lower limit of basal expression levels, maps to the size of the aptamer domain. Our observations led to the specification of two aptamer tuning strategies: targeted changes in aptamer affinity without changing aptamer size alter Kapt and targeted changes to aptamer size to alter the processing efficiency of the switch (e). Taken together, variation of Kapt and e map to the aptamer domain and depend on the nature of the sequence modification.
We examined whether placement of new aptamers into the aptamer domain imparts new ligand dependence while preserving shRNA switch functionality. Previous RNA‐based regulatory platforms have demonstrated alteration of ligand dependence by the modular incorporation of new aptamers (Bayer and Smolke, 2005; Win and Smolke, 2007) or minimal mutation of the base aptamer (Thompson et al, 2002; Desai and Gallivan, 2004). We evaluated the xanthine aptamer, as it produced low basal expression levels and tightly binds the water‐soluble and non‐cytotoxic small molecule hypoxanthine. Following construction of shRNA switches that incorporate the xanthine aptamer by direct replacement of the aptamer domain, we generated component transfer functions in HEK293T cells stably expressing EGFP. As observed for S1, intermediate basal expression levels of GFP increased in a dose‐dependent manner that was abolished by mutating the aptamer core (Figure 4D and E). Furthermore, the competing strand tuning strategies were preserved as evidenced by the effect of changing the competing strand length on the hypoxanthine response curves (Figure 4F and G). Contrary to model predictions, mutation of the xanthine aptamer resulted in increased basal expression levels, which may be attributed to base‐pairing interactions between the mutated aptamer and the competing strand or shRNA stem sequences, or changes in aptamer folding and stability. However, the shift in basal levels upon mutation of the aptamer sequence is less than that observed for changes in the competing strand, supporting the conclusion that our model serves as a sufficient first approximation. Thus, our shRNA switch design can accommodate different aptamers to alter the identity of the molecular input that regulates gene expression.
Programming transfer functions by combining competing strand or aptamer tuning strategies
The ligand‐regulated behaviour of shRNA switches can be programmed through application of the competing strand and aptamer tuning strategies described above. If these programming strategies could be combined, then a collection of shRNA switches could be constructed that display finely tuned transfer functions and respond to a range of molecular inputs. Such capabilities will be integral to the construction of higher‐order biological networks that display multi‐input control over gene expression.
On the basis of the independence of the competing strand tuning strategies, we examined whether the strategies can be combined to fine‐tune the component transfer function beyond the capabilities of any single strategy. To generate small deviations in the transfer function of a parent shRNA switch, we added compensatory nucleotide changes under each competing strand tuning strategy in a stepwise manner (Figure 5A and B): a point mutation (G68A) within the competing strand to increase complementarity, deletion of two base pairs to decrease the competing strand length at the 5′‐end, and a single insertion at the 3′‐end to increase the competing strand length. Each nucleotide change yielded the expected shift in the transfer function corresponding to the relative stabilization (increased Kcomp) or destabilization (decreased Kcomp) of the inactive conformation. The final switch, S10, displayed a transfer function slightly shifted from that of the parent switch, S4, demonstrating that nucleotide changes following the three competing strand tuning strategies can be combined to yield fine‐tuning of the component transfer function.
In addition to fine‐tuning of a single shRNA switch, combining shRNA switches with different ligand dependencies would contribute to the construction of networks that integrate multiple molecular inputs, as suggested from recent work on siRNA‐based logic evaluator systems (Rinaudo et al, 2007). To evaluate the efficacy of combining shRNA switches, we transfected HEK293T cells stably expressing EGFP with shRNA switches that incorporate the theophylline aptamer (S4) or the xanthine aptamer (X1), where both switches target EGFP and display similar basal expression levels. On the basis of the combined component transfer functions (Materials and methods), we anticipated that the combined regulatory effects of S4 and X1 would require the presence of both hypoxanthine and theophylline to fully inhibit GFP silencing (Figure 5C and D). Relative GFP levels were measured for cells transfected in the presence of theophylline, hypoxanthine, or the combination of the two (Figure 5E). GFP levels were high when individual switches were paired with their cognate ligand or both ligands. Some signal cross talk was observed as evidenced by the lower responsiveness of S4 and X1 to their non‐cognate ligands, and primarily attributed to low aptamer specificity as observed in previous theophylline aptamer studies (Jenison et al, 1994) and not explicitly tested for the xanthine aptamer (Kiga et al, 1998). When the switches were cotransfected, high GFP levels coincided only in the presence of both ligands as expected on the basis of the circuit configuration. On the basis of the results, shRNA switches allow the construction of finely tuned genetic networks that can process multiple inputs.
An in silico framework towards component sequence‐to‐transfer function prediction
The construction of large‐scale biological systems will require the simultaneous optimization of the behaviour of all system components to yield proper network behaviour as suggested for natural (Suel et al, 2007) and synthetic (Gardner et al, 2000; Yokobayashi et al, 2002) systems. Although the transfer functions associated with shRNA switches and other synthetic riboswitches are amenable to physical tuning, a computational framework to effectively navigate through qualitatively functional sequences is necessary for the rapid optimization of switch performance. Folding energetics dictate conformational partitioning and therefore switch performance for a strand displacement mechanism. RNA secondary structure prediction algorithms (Mathews et al, 2004) have the potential to perform accurate in silico prediction of in vivo switch performance, although these algorithms have not been sufficiently tested for in vivo folding dynamics. To investigate the applicability of the secondary structure algorithms to predict in vivo switch behaviour, we sought to develop a sequence‐to‐function relationship for shRNA switches using such algorithms in combination with our model (Figure 6A).
On the basis of our tuning analysis, we identified Kcomp as the sole parameter that reflects partitioning between active and inactive conformations and maps to the competing strand. The free energy difference (ΔG) between conformations is directly related to Kcomp such that transfer function prediction is possible by calculating ΔG from sequence information with the aid of structure prediction algorithms, converting this value into Kcomp, and inserting Kcomp into Equation (1) to quantitatively relate ligand concentration and target gene expression levels. A fully determined model requires values for the remaining model parameters; as these parameters are not currently amenable to calculation in silico, experimental estimation can be conducted with a minimal set of experiments on the basis of our model construction (Supplementary text 1).
We first determined if ΔG values calculated from the algorithm correlate with the measured basal expression levels for shRNA switches with varying competing strand sequences. The implicit assumption is that competing strand alterations affect only conformational partitioning, which can be calculated with the structure prediction algorithms. We evaluated ΔG (ΔGmethod) by separating active and inactive conformations on the basis of the minimal free energy (MFE) and the weighted energies from a partition function (PF) calculation (Supplementary text 2 and Supplementary Figure S5), where both methods are commonly used to evaluate RNA folding in vitro and in vivo. These methods were employed to calculate ΔGmethod for shRNA switches S1–10, which differ only in their competing strand sequence. To measure the correlation strength between ΔGmethod and basal expression levels for either method, we performed a least‐squares fit using a three‐parameter equation of the same form as our model with both data sets. Ideally, the fit relationship between ΔGmethod and measured basal expression levels should align with the same relationship predicted by the model (Figure 6B), where ΔG (ΔGmodel) is related to Kcomp according to Equation (3). Both MFE and PF calculations failed to provide a significant correlation between ΔGmethod and basal expression levels (Supplementary Figure S6), suggesting that these methods are insufficient for accurate prediction of RNA folding dynamics in vivo.
For all competing strand tuning strategies, increasing the stability of the inactive conformation always resulted in an increase in basal expression (Figure 3B–G). The MFE and PF methods did not effectively capture each energetic shift potentially due to the inclusion of binding interactions outside of the major stems. We hypothesized that the interactions outside of the competing strand domain are less prevalent in vivo and are biasing the energetic calculations. To examine this possibility, we devised a third method, the ‘Stems’ method, that accounts only for the energetics of the major stem in each conformation (Supplementary Figure S5). Implementing the ‘Stems’ method resulted in a strong correlation (R2=0.94) between basal expression levels and ΔGmethod (Figure 6C).
Despite the absence of a perfect overlap between the correlation of the ‘Stems’ method and that predicted by our model (Figure 6B and C), the correlation established a significant empirical link between shRNA switch sequence and behaviour in the absence of ligand. This correlation can be assimilated into the model by equating basal expression levels predicted by the fit equation and the model to determine the relationship between ΔGmethod and Kcomp (Supplementary text 2). Doing so yields a predictive component transfer function that is now dependent on the calculated value of ΔGmethod:
where C1−3 are empirical constants from the fit correlation. Our extended model provides a general framework for predicting shRNA switch transfer functions from sequence information, where energetic values produced from structure prediction algorithms are inserted into the model for the prediction of switch behaviour. Although the extended model currently requires parameter fitting to yield the predicted relationship between ligand concentration and target gene expression levels, the framework establishes a starting point for the development of methods that rely on in silico calculations for transfer function prediction from sequence information.
Model‐guided forward design of shRNA switches with optimized transfer functions
To apply the extended model to the forward design of shRNA switches with defined functional properties, we sought to design a theophylline‐regulated shRNA switch displaying a maximized dynamic range (η). η is defined as the ratio of GFP levels at high (3 mM) and low (1 μM) theophylline concentrations. We used our extended model to calculate the range of ΔG values where η is maximized. Model predictions suggest that η is maximized for switches with ΔGmethod∼−3 kcal/mol and that use of the smaller theophylline aptamer (higher e) yields a higher maximum (Figure 6D). To evaluate the predicted landscape, we designed new shRNA switches (S13–25) that include the smaller theophylline aptamer and display ranging ΔG values, generated component transfer functions and calculated η. Plotting ΔGmethod against the measured value of η for all theophylline‐regulated shRNA switches (S1–25; Figure 6E) shows a maximum for switches containing the smaller theophylline aptamer that is higher than that for the switches containing the larger aptamer. Furthermore, both maxima existed at ΔGmethod∼−3 kcal/mol as predicted by the extended model supplemented with the empirical parameter values, and the best switch (S13) was approximately equal to the theoretical maximum of η according to model predictions (ηmax,theor∼5). Flow cytometry data illustrate the improvement in dynamic range (Supplementary Figure S7) for the best shRNA switch (S13; Figure 6G) as compared to the original shRNA switch (S1; Figure 6F).
To examine the generality of shRNA switch design and functionality, we designed a set of shRNA switches targeting the endogenous La protein. Following selection of an shRNA sequence that yielded moderate knockdown of La as assayed by qRT–PCR, six shRNA switches (L1–6) were constructed with the smaller theophylline aptamer and various competing strand sequences covering a range of ΔG values. Each shRNA switch exhibited variable response to 1.5 mM theophylline that was not observed for the base shRNA (Supplementary Figure S8). As observed for the GFP‐targeting shRNA switches, use of the ‘Stems’ method provided a suitable correlation between basal levels and ΔGmethod, whereas the MFE and PF methods did not. Supplying the model with fit values for C1−3 yielded a predicted dynamic range trend that closely matched the experimental data (Figure 7A and B). Interestingly, when the values of fshRNA and e calculated from the La‐targeting base shRNA and an shRNA switch preferentially adopting the active conformation (L6) were combined with the remaining parameter values from the GFP experiments, the resulting trend predicted the same maximum dynamic range with a shifted value of ΔGmethod corresponding to the maximum dynamic range. The results suggest that the shRNA sequence affects calculations with the ‘Stems’ method such that empirical values are specific to individual sequences and experimental conditions. However, the ‘Stems’ method produced a strong correlation such that the model may be implemented in future designs by generating a small set of shRNA switches covering ΔGmethod values of approximately −5 to 0 kcal/mol and measuring basal expression levels. Thus, shRNA switches can be constructed to target different genes, and the model can be used as a tool to guide the forward design of switches displaying optimal behaviour.
A comparison between the framework described here and a recently described ligand‐controlled shRNA system (An et al, 2006) highlights important design strategies to engineer domain swapping and tuning of the transfer function into synthetic riboswitch systems. In the previous design, ligand control of RNAi was achieved through direct coupling of the theophylline aptamer and an shRNA stem. This design inherently limits aptamer swapping, as the aptamer must perform ligand binding coordinated with modulation of Dicer processing, and prevents tuning of the transfer function, as sequence changes that modulate Dicer processing cannot be implemented without a complete loss of ligand responsiveness. In contrast, we propose here a framework on the basis of the coupling of three distinct domains that carry out separate functions necessary to convert ligand binding into modulation of RNAi activity. Our system requires that the aptamer performs one function, ligand binding, and the modulation of RNAi processing is performed by a separate domain, the competing strand. The competing strand permits fine‐tuning of the transfer function and enables modular coupling of the aptamer and shRNA stem domains, as confirmed by independently replacing each domain and demonstrating preservation of functionality.
We developed a model to enhance our understanding of shRNA switch activity and identified five tuning strategies reflected in three model parameters, Kcomp, Kapt, and e, that map specifically to sequence changes in the competing strand or aptamer domains. Our model also established important shRNA switch design guidelines. The first is that basal expression levels are determined by a collection of factors: shRNA potency (fshRNA), shRNA switch processability (e) and prevalence of the active conformation (Kcomp). To achieve a desired basal expression level, all factors must be considered in the switch design. Another guideline originates from the observation that larger aptamers coincided with increased basal expression levels, potentially due to sterically hindering processing by the RNAi machinery. The specific contribution of secondary or tertiary structure to the inhibitory effect is unclear, although further understanding of how the RNAi machinery specifically interacts with the shRNA through crystallographic or mutational studies may shed light on this dependence. Our results suggest that shRNA switch sequence length has an upper limit before compromising activity, where future engineering efforts may focus on alleviating or entirely removing this limitation. Furthermore, if achieving low basal expression levels is critical and a set of aptamers against the same ligand are available, use of smaller aptamers may be preferred even at a cost to aptamer affinity. Such a guideline may even direct library design for the selection of new aptamers by placing an upper limit on the length of the randomized sequence.
We incorporated RNA folding algorithms into our model for in silico prediction of shRNA switch behaviour in vivo. The resulting model yielded a framework for the forward design of shRNA switches with specified functional properties. This was achieved by linking RNA secondary structure prediction algorithms, which convert sequence information into energetic values, to our model, which converts energetic values into switch behaviour, to provide an empirical sequence–function relationship. The specific method used to calculate the free energy difference (ΔGmethod) between active and inactive conformations deviated from commonly used methods (MFE and PF calculations) on the basis of observations from the experimental tuning trends. Our alternative method may provide a better correlation with experimental results by focusing the prediction of Kcomp to the region of the switch in which the competing strand binding events are occurring, ignoring energetic contributions from other regions of the switch molecule that may not be relevant to the in vivo conformational switching process. Our analysis moves towards direct sequence‐to‐function relationships and suggests that commonly used methods for predicting RNA structure and behaviour should be carefully evaluated when applied to in vivo environments. RNA folding in vivo is a complex process, and algorithms that account for folding kinetics (Danilova et al, 2006) and ulterior structural formation (Parisien and Major, 2008), such as pseudoknots or non‐canonical base‐pairing interactions, may increase the accuracy of the model as well as provide insight into sequences that deviate from model predictions (Figures 6E and 7B). For practical application, newer algorithms will need to be more fully developed to offer the same functionalities as existing algorithms, such as the ability to rapidly scan suboptimal structures, calculate the energetics of multiple RNA strands, and perform a partition function calculation. Although the PF method did not produce a strong correlation using existing algorithms, non‐canonical base‐pairing interactions may have been an important factor that will be accounted for with newer algorithms.
On the basis of the demonstrated modularity and tunability of our platform, shRNA switches can be implemented towards various applications. The required dynamic regulatory range of a given application will be one of the main considerations in utilizing shRNA switches, as the switches are practically limited to an ∼10‐fold induction ratio on the basis of the maximum achievable knockdown with an endogenously expressed shRNA. However, many endogenous non‐coding RNAs, including miRNAs, exhibit similar restrictions on dynamic regulatory range and have important functions in diverse biological processes, suggesting that this limited dynamic range is not absolutely restrictive to the utility of shRNA switches as dynamic gene regulatory components.
As one potential application, shRNA switches can be applied to disease therapy by sensing intracellular disease markers and inducing apoptosis or cell cycle arrest only in the affected cells as suggested previously (Rinaudo et al, 2007). When a context‐dependent concentration threshold divides diseased and normal cells, tunability is essential to reduce the likelihood of false positives or negatives. In addition, shRNA switches can be integrated into synthetic genetic circuits to generate advanced control schemes in biological systems. Such systems often exhibit complex dependencies on the dynamics of component interactions, and tuning of component behaviour is often necessary to achieve optimal system performance. Through the fine‐tuning strategies and model‐guided forward design tools described here, shRNA switches may be used to address challenges faced in biological network design and serve as complex regulatory components in synthetic biology (Endy, 2005; Keasling, 2008; Savage et al, 2008).
Materials and methods
All shRNAs were cloned into pSilencer 2.1‐U6 puro (Ambion). The original shRNA present in pSilencer was used as a scrambled shRNA control. The pSilencer backbone was modified to co‐express DsRed‐Express in 293T cells by cloning the SV40 origin of replication, CMV IE promoter, and DsRed‐Express into the NsiI/MfeI restriction sites. The original XhoI site present in the backbone was removed by XhoI cleavage, extension with the large Klenow fragment (New England Biolabs), and ligation. To clone the shRNA switches, the original shRNA followed by a 6‐nucleotide (nt) string of T's was cloned into BamHI/HindIII directly downstream of the U6 promoter. The original shRNA was converted into an shRNA switch by cloning the remaining sequence (Supplementary Table I) into XhoI/XbaI contained within the shRNA loop region. This allowed cloning in parallel of multiple shRNA switches that are comprised of the same shRNA region. All cloning steps involved annealing of 5′‐phosphorylated synthetic oligonucleotides (Integrated DNA Technologies) and ligation into the backbone vector. All restriction enzymes and T4 DNA ligase were purchased from NEB. All constructs were sequence‐verified (Laragen Inc.), where sequences are provided in Supplementary Table I.
Preparation of RNAs
S4t was transcribed in vitro from an annealed template containing the T7 promoter (5′‐TTCTAATACGACTCACTATAGGG‐3′, where G is the first transcribed nucleotide) using the Ampliscribe T7 transcription kit (Epicentre) according to the manufacturer's instructions. Following transcription and DNase treatment, unincorporated NTPs were removed using a NucAway clean‐up column (Ambion). 5′‐phosphates were subsequently removed using Antarctic phosphatase (NEB). Dephosphorylated RNA was then gel‐purified on a 6% denaturing polyacrylamide gel and quantified using an ND‐1000 spectrophotometer (NanoDrop). RNAs were 5′‐radiolabeled using T4 PNK (NEB) and [γ‐32P]‐ATP, purified using a NucAway clean‐up column, and gel‐extracted on a 6% denaturing polyacrylamide gel.
In‐line probing was conducted as described previously (Soukup and Breaker, 1999). After heating at 70 °C for 2 min followed by slow cooling to room temperature, 5′‐radiolabeled RNAs (0.2 pmol) were incubated for 40 h at 25 °C in varying amounts of theophylline with 50 mM Tris–HCl pH 8.5 and 20 mM MgCl2. Reactions were terminated by adding an equal volume of loading buffer (10 M urea, 1.5 mM EDTA). The alkaline hydrolysis ladder was generated by incubating RNA in 50 mM NaHCO3/Na2Co3 pH 9.2 and 1 mM EDTA for 6 min at 95 °C. The G‐specific cleavage ladder was generated by incubating RNA in 1 U RNase T1 (Ambion) with 20 mM sodium citrate pH 5.0, 1 mM EDTA, 7 M urea, and 3 μg yeast RNA for 25 min at 25 °C. RNAs were resolved on an 8% denaturing polyacrylamide gel, dried for 90 min at 70 °C, then visualized on an FX phosphorimager (Bio‐Rad). Band quantification was performed using the Quantity One software package (Bio‐Rad). To account for well‐loading variability, quantified band intensities were normalized to an adjacent band of similar intensity showing negligible theophylline dependence.
Cell culture and transfection
All cells were maintained at 37 °C in a 5% CO2‐humidified incubator. HEK293T, HEK293, HeLa, and HEK293T tTA‐d2EGFP cells were maintained in minimal essential medium alpha media (Invitrogen) supplemented with 10% fetal bovine serum (FBS) (Invitrogen), whereas MDA‐MB‐231 cells were maintained in RPMI 1640 with glutamine (Invitrogen) supplemented with 10% FBS. Cells were transfected 1 day after seeding using Fugene 6 (Roche) according to the manufacturer's instructions, followed by the immediate addition of ligand. HEK293T tTA‐d2EGFP were transfected with shRNA vector (250 ng), whereas cells lacking endogenous GFP were cotransfected with shRNA vector (250 ng) and pcDNA3.1(+) (Invitrogen) harbouring the d2EGFP gene (25 ng) (Clontech). One day post‐transfection, the media and ligand were replaced. Transfected cells were collected 3 days post‐transfection for flow cytometry analysis.
Cell fluorescence analysis
Three days post‐transfection, cells were trypsinized and subjected to flow cytometry analysis using the Cell Lab Quanta SC MPL (Beckman Coulter). Cells were first gated twice for (1) viability as assessed by electronic volume versus side scatter and (2) green fluorescence above autofluorescence to remove a non‐fluorescent subpopulation. Cells were then gated for either low or high DsRed‐Express fluorescence, representing untransfected or transfected cells, respectively. To minimize well‐to‐well variability, the median green fluorescence value of transfected cells were divided by that of untransfected cells in the same well and reported as GFP (%). For cells cotransfected with shRNA and GFP plasmids, GFP (%) is the relative GFP levels when normalized to mean red fluorescence followed by normalization to cells transfected with the scrambled shRNA. See Supplementary Figure S9 for representative plots and the corresponding gates for transfected and untransfected cells.
Modelling and RNA energetic calculations
Calculation of RNA free energy and partition functions were performed using RNAStructure (Mathews et al, 2004). Kcomp and the energy difference between inactive and active conformations are related by the following expression:
where NA is Avogadro's number, kB is the Boltzmann constant, and T is temperature (K). See Supplementary texts 1 and 2 for a full description of the model derivation, methods for calculating folding energetics, and prediction of the transfer function for a given shRNA switch sequence. Equation fits to measure the correlation strength between ΔGmethod and basal expression levels were performed by least‐squares analysis using the following expression that has the same mathematical form as Equation (1):
where C1–3 are the fit constants. Supplementary Table II contains energetic values calculated under each method along with experimentally determined expression levels.
To model the multi‐input system attained by cotransfecting equimolar concentrations of plasmids harbouring S4 and X1, the knockdown achieved by the original shRNA was halved to reflect the 50% decrease in delivered plasmid DNA. The ligand‐dependent contributions to decreased expression levels were combined into a single expression to reflect the additive nature of shRNA levels mediating knockdown of the target gene:
where the subscripts S4 and X1 designate parameter values generated by a fit to the corresponding individual component transfer functions, and Lth and Lxan represent the exogenous theophylline or hypoxanthine concentration, respectively.
The authors declare competing financial interests in the form of a pending patent application.
We thank G Soukup and D Endy for critical reading and comments on this manuscript, P Aebischer for providing the HEK293T tTA‐d2EGFP cells, and S Culler for providing assistance with qRT–PCR. This work was supported by the Caltech Joseph Jacobs Institute for Molecular Engineering for Medicine (grant to CDS), the Defense Advanced Research Projects Agency (grant to CDS), the National Institutes of Health (training grant to TSB; fellowship to KGH), and the Department of Defense (fellowship to CLB).
Supplementary Information [msb200862-sup-0001.pdf]
This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.
- Copyright © 2008 EMBO and Nature Publishing Group