Synthetic in vitro transcriptional oscillators.

The construction of synthetic biochemical circuits from simple components illuminates how complex behaviors can arise in chemistry and builds a foundation for future biological technologies. A simplified analog of genetic regulatory networks, in vitro transcriptional circuits, provides a modular platform for the systematic construction of arbitrary circuits and requires only two essential enzymes, bacteriophage T7 RNA polymerase and Escherichia coli ribonuclease H, to produce and degrade RNA signals. In this study, we design and experimentally demonstrate three transcriptional oscillators in vitro. First, a negative feedback oscillator comprising two switches, regulated by excitatory and inhibitory RNA signals, showed up to five complete cycles. To demonstrate modularity and to explore the design space further, a positive-feedback loop was added that modulates and extends the oscillatory regime. Finally, a three-switch ring oscillator was constructed and analyzed. Mathematical modeling guided the design process, identified experimental conditions likely to yield oscillations, and explained the system's robust response to interference by short degradation products. Synthetic transcriptional oscillators could prove valuable for systematic exploration of biochemical circuit design principles and for controlling nanoscale devices and orchestrating processes within artificial cells.

Introduction

Fundamental goals for synthetic biology are to understand the principles of biological circuitry from an engineering perspective and to establish engineering methods for creating biochemical circuitry to control molecular processes—both in vitro and in vivo (Benner and Sismour, 2005; Adrianantoandro et al, 2006). As the canonical example of non-equilibrium chemical dynamics—and as a design target for synthetic systems—for more than 50 years chemical oscillators have had a central role in shaping our understanding of chemical self-organization and the mechanistic origins of life-like behavior in non-living systems (Goodwin, 1963; Zhabotinsky, 1991; Winfree, 2001). Known chemical oscillators are grouped into five classes: (1) biological oscillators within living cells, such as the circadian rhythm (Panda et al, 2002); (2) biological oscillators reconstituted in vitro (Nakajima et al, 2005; Mori et al, 2007); (3) designed synthetic oscillators engineered into living organisms (Elowitz and Leibler, 2000; Atkinson et al, 2003; Stricker et al, 2008); (4) synthetic chemical oscillators involving small molecule reactions in vitro (Zhabotinsky, 1964; Epstein and Pojman, 1998); and (5) synthetic biochemical oscillators designed from biological components in cell-free in vitro reactions (Wlotzka and McCaskill, 1997; Kim, 2007; Montagne et al, 2011).

The ability to engineer the circuit architecture of synthetic oscillators makes it possible to investigate design principles by exploring the design space (Barkai and Leibler, 2000; Novák and Tyson, 2008; Savageau et al, 2009), to challenge modeling techniques with well-defined biochemical systems of intermediate complexity (Simpson, 2006; Gutenkunst et al, 2007; Cantone et al, 2009), and to orchestrate other molecular processes within natural and artificial chemical systems (Weber and Fussenegger, 2009; Liedl and Simmel, 2005). The third class of oscillators—using in vivo synthetic biology—has already been enormously productive for these reasons. However, a direct comparison between synthetic oscillator designs in vivo remains a challenge because of differences in their regulatory components as well as potential interference with other cellular networks; moreover, in vivo oscillators are not helpful for engineering systems that must avoid using biological organisms within them.

Consequently, the fifth category—using cell-free in vitro synthetic biology—is particularly interesting because reactions such as transcription and translation can be rewired combinatorially using synthetic DNA templates, the resulting systems can be studied and characterized without the complexities and unknowns of living cells, and the creation of artificial chemical systems with complex autonomous dynamics becomes possible (Simpson, 2006; Forster and Church, 2007; Bujara et al, 2010). However, initial attempts to construct a cell-free biochemical oscillator using transcription and reverse transcription (Wlotzka and McCaskill, 1997) were only moderately successful, perhaps because of accumulated sequence mutations. Since then, multistep reaction pathways using transcription and translation (Noireaux et al, 2003), bistable circuit dynamics using RNA transcription and degradation (Kim et al, 2006), and logic gates using multiple enzymes (Takinoue et al, 2008) have been experimentally demonstrated, and highly efficient cell-free platforms for transcription and translation are now available (Shimizu et al, 2001; Jewett et al, 2008). Further, theoretical models developed for these systems are capable of sustained oscillations (Ackermann et al, 1998; Kim et al, 2004; Simpson et al, 2009; Takinoue et al, 2009), suggesting that a variety of biochemical circuit architectures can in principle be synthesized and explored in cell-free reactions.

Here, we make use of a previously proposed class of in vitro biochemical systems, transcriptional circuits, which can be modularly wired into arbitrarily complex networks by changing the regulatory and coding sequence domains of DNA templates. In principle, transcriptional circuits can be wired as continuous time analog neural networks with symmetric or asymmetric weights (Hopfield, 1984), implying that they are a computationally and behaviorally complete circuit architecture (Kim et al, 2004). Individual transcriptional switches have been demonstrated experimentally to exhibit sharp sigmoidal inhibitory regulation, allowing the construction of a two-switch circuit exhibiting bistable dynamics (Kim et al, 2006). Furthermore, a single switch with sharp sigmoidal excitatory regulation has also been demonstrated; it also exhibits bistability when configured to regulate itself (Subsoontorn et al, 2011). Here we combine switches with inhibitory and excitatory regulation to construct and characterize three different oscillator designs with circuit architectures similar to known synthetic and natural oscillators: a two-switch oscillator utilizing both excitatory and inhibitory regulation, loosely analogous to the p53–Mdm2-feedback loop (Bar-Or et al, 2000); the same oscillator augmented with a positive-feedback loop, loosely analogous to a synthetic relaxation oscillator (Atkinson et al, 2003); and a three-switch ring oscillator analogous to the repressilator (Elowitz and Leibler, 2000).

Simplified mathematical models predicted that all three designs could be tuned to reach the oscillatory regime in parameter space and were useful for understanding the oscillators' principles of operation, but were not in quantitative agreement with experiments. More detailed chemical kinetic models proved essential for guiding the exploration of experimental parameters to achieve oscillations, but further extensions to the model were required to account for the build-up of incomplete RNA degradation products that interfered with the designed interactions. To our surprise, the two-switch oscillator was remarkably robust to the accumulation of this interfering RNA ‘waste'. The successful implementation of three oscillator designs underscores the potential of in vitro circuitry, but our experience also serves as a warning of the challenges posed by unintended reactions even after the uncertainties and complexities of living cells have been removed.

Results

Design I: a two-switch negative feedback oscillator

Our in vitro transcriptional switches consist of a synthetic DNA molecule with a regulatory domain, a promoter, and an output domain (Figure 1A, center top and bottom). The mechanism for switch function depends primarily on the Watson-Crick complementarity of key sequence domains, allowing the same switch motif to be implemented with a variety of sequences, thereby enabling the straightforward construction of circuits involving many switches. Specifically, each synthetic switch Swij is controlled by an input signal, RNA species j, and produces an output signal, RNA species i. The response of a switch to its input can exhibit a sharp threshold, which derives from four types of strong DNA and RNA hybridization reactions that we call activation, annihilation, inhibition, and release (Figure 1A, solid line arrows).

Figure 1

Schematics for in vitro transcriptional oscillators. (A) Reaction diagram for the two-switch negative-feedback oscillator (Design I). On the top left is a block diagram, wherein arrowheads indicate activation or production and circular ends indicate inhibition. The block diagram corresponds to the detailed diagram as highlighted by gray shading: T21A1 (ON-state switch Sw21) and T21 (OFF-state switch Sw21) are summarized by the Sw21 block; RNA inhibitor rI2, together with its threshold, DNA activator A2, and their complex, A2rI2, is summarized by the rI2 block; and similarly for the Sw12 and rA1 blocks. The sequence domains are color coded to indicate identical or complementary sequences; for the switch templates, the dark blue sequence domain inside the rectangle indicates the T7 RNAP promoter sequence with arrows pointing in the direction of transcription, with transcription domains indicated by light blue dashed circles. For fluorescence monitoring, OFF-state switches are labeled with fluorophores, T21 with Texas Red (red circle) and T12 with TAMRA (green circle), and both activators A1 and A2 are labeled with Iowa Black RQ quenchers (black circle). Four types of hybridization reactions are indicated by arrows: activation (magenta), inhibition (orange), annihilation (brown), and release (blue). Hybridization reactions are reversible; the arrowhead corresponds to the thermodynamically favorable direction, and the reverse reactions are expected to be so slow as to be negligible (Wetmur, 1991; Yurke and Mills, 2003; Zhang and Winfree, 2009). Transcription by RNAP is shown as black dashed arrows and degradation of RNA within RNA-DNA hybrids by RNase H is shown as black dotted arrows. DNA and RNA sequences of single-stranded species and complexes for all three oscillators are shown in Supplementary Figure S1. (B) List of hybridization and enzyme reactions. See Supplementary information section 1.4 for details. (C) Theoretical end states of hybridization reactions in the absence of enzymes. As the input RNA inhibitor rI2 concentration increases, initially the free DNA activator A2 is consumed without affecting switch state. When all free A2 is consumed (i.e., [rI2]=[A2tot]–[T12tot]), rI2 displaces A2 from the T12A2 complex in stoichiometric amounts until all A2 is consumed (i.e., [rI2]=[A2tot]), resulting in a piecewise linear graph. Similarly, the response of switch Sw21 to rA1 input is a piecewise linear graph. See Supplementary information section 1.1 for details. (D) Block diagrams for the amplified negative-feedback oscillator (Design II) and the three-switch ring oscillator (Design III). See Supplementary Figures S2 and S3 for detailed reaction diagrams.

This threshold mechanism is analogous to biological threshold mechanisms such as ‘inhibitor ultrasensitivity' (Ferrell, 1996) and ‘molecular titration' (Buchler and Louis, 2008). It is sufficiently sharp to allow arbitrary analog or digital circuits to be created in principle (Kim et al, 2004). The OFF state of the switch consists of a double-stranded (ds) DNA template (‘T') with a partially single-stranded (ss) and thus incomplete T7 RNA polymerase (RNAP) promoter region. The switch is turned ON by the binding of an ssDNA activator (‘A') that completes the RNAP promoter region (reaction 1, activation). The resulting template (‘T·A') has a nicked promoter but still transcribes well (Kim et al, 2006). An inhibitor strand (either ssRNA, ‘rI', or ssDNA, ‘dI') can bind to a complementary free-floating activator (either ssDNA, ‘A', or ssRNA, ‘rA'), resulting in a functionally inert activator–inhibitor complex such as ‘A·dI' (reaction 2, annihilation). Thus, a population of inhibitor strands in solution can absorb any free activator strands in solution, leaving either only inhibitor strands or only activator strands, depending on which population was initially larger. In addition, because the activator contains a ‘toehold', an ss overhang on T·A, an inhibitor can bind and initiate a strand displacement reaction (Yurke and Mills, 2003) that removes A from the template. Thus, the switch can be turned OFF upon addition of an inhibitor strand, if it exceeds the amount of free activator in solution (reaction 3, inhibition). The DNA inhibitor dI, in turn, has a toehold region on A·dI that allows RNA activator rA to bind and initiate a toehold-mediated strand displacement reaction that releases A to activate target switches (reaction 4, release). Together, these reactions (Figure 1B, top) ensure that there is always a fast kinetic pathway to the end state of the hybridization reactions.

These four hybridization reactions can be assembled to create either an inhibitable switch (Figure 1A, right and bottom) with a threshold set by the total concentration of its DNA activator strand (Figure 1C, bottom) or an activatable switch (Figure 1A, left and top) with a threshold set by its DNA inhibitor strand concentration (Figure 1C, top). (In an activatable switch, the DNA activator strand concentration should be roughly comparable with the template concentration; it should be at least as high, so that all the templates can be turned ON, but it need not be higher, as excess activator merely disables a stoichiometric amount of inhibitor. Following Subsoontorn et al (2011), activatable switches were designed to use indirect activation because direct complementation of the missing promoter region by an ssRNA activator would have resulted in a weaker activation than that provided by an ssDNA activator (McGinness and Joyce, 2002). Conceptually, species are separated into those involving a single signal species (Figure 1A, left and right)—which we call a signal block—and those involving the transcription of one signal species regulated by a second (Figure 1A, top and bottom)—a switch block. The thresholding reactions occur in the signal blocks, which must each function either exclusively as an activator or exclusively as an inhibitor.

For a molecular implementation, the lengths and sequences of each domain must be designed well to obtain the desired hybridization reactions. This entails ensuring that (1) for each reaction, the resulting complex is thermodynamically more favorable than the starting complex (Wetmur, 1991), (2) strand displacement reactions are initiated by sufficiently strong toeholds (Yurke and Mills, 2003; Zhang and Winfree, 2009), and (3) there is minimal unintended secondary structure or spurious binding between unrelated strands. (Sequences for this work were largely derived from those in Kim et al (2006) and Subsoontorn et al (2011). See Supplementary Figure S1 for sequences of DNA complexes for all three oscillator designs.)

Using these design motifs for switch and signal blocks, networks with arbitrary connectivity can be constructed modularly. Given independent, well-designed sequences for each signal block, a switch block connecting any chosen pair may be created simply by using the upstream signal block's sequence in the regulatory domain and using the downstream signal block's sequence in the output domain. Thus, while Figure 1A shows a specific network wiring, other block diagrams can be systematically translated into DNA implementations by replacing each block with its DNA reaction motif exactly as in this figure. For instance, consider the molecular implementation of a switch block, Sw11, acting as a positive-feedback loop (Figure 1D, Design II). This switch block is implemented as a DNA template using the signal sequence rA1 for both the activatable regulatory domain and the output domain. Similarly, consider the molecular translation of the block diagram for a ring oscillator (Figure 1D, Design III). The inhibitable regulatory motif (as illustrated for rI2 regulating Sw12) can be applied to create inhibitable switch templates for Sw31 and Sw23 to meet the design specifications of the block diagram. For example, the molecular realization of Sw23 is simply a DNA template using the signal sequence rI3 for the inhibitable regulatory domain and the signal sequence rI2 for the output domain. (See Supplementary Figures S2 and S3 for detailed reaction diagrams.)

In a typical reaction network, RNA outputs will be produced by RNAP from upstream switches using NTP as fuel, and these outputs will serve as inputs for downstream switches. At the same time, the degradation of RNA signals by Escherichia coli ribonuclease H (RNase H) releases the DNA strands from the functionally inert DNA–RNA hybrid state, thereby undoing the regulatory effect of the RNA inputs and allowing the downstream switch to revert to its native state (Figure 1B, bottom). Consequently, the dynamics of RNA signal concentrations result from the balance (or imbalance) of transcription and degradation processes.

Using the inhibitable and activatable transcriptional switch motifs, we constructed a two-switch negative feedback oscillator with the connectivity shown in Figure 1A (inset). A total of seven DNA strands are used, in addition to the RNAP and RNase H in the reaction buffer. RNA activator rA1 activates the production of RNA inhibitor rI2 by modulating switch Sw21, whereas RNA inhibitor rI2, in turn, inhibits the production of RNA activator rA1 by modulating switch Sw12. The fact that such a negative-feedback loop can lead to temporal oscillations can be seen from a mathematical model of transcriptional networks. Consider a four-dimensional model, termed the ‘simple model', the dynamic variables of which are the two RNA concentrations and the concentrations of the two ON-state switches. We assume that the production rates of the two RNA species are solely determined by the ON-state switch concentrations and their degradation rates depend only on their own concentrations, as follows:

where k is the apparent first-order rate constant assuming a fixed RNAP concentration in excess, and similarly k is the apparent first-order rate constant for RNase H. Changes in ON-state switch concentrations are in turn governed, through hybridization reactions (1–4), by the amounts of RNA activator and inhibitor species present. The steady-state switch response to RNA inputs can be reasonably approximated by Hill functions (Kim et al, 2006), and for simplicity we take that approach here:

where n and m are the apparent Hill exponents, τ is an apparent relaxation time for the hybridization reactions, K is the threshold set by DNA activator A2, K is the threshold set by DNA inhibitor dI1, and [Tijtot] is the sum of concentrations of all molecular species containing Tij—i.e., the sum of ON-state switch concentration [TijAj] and OFF-state switch concentration [Tij]. As shown in Figure 1C, reasonable approximations for the thresholds are given by and . Apparent Hill exponents are determined by the relative levels of the switch template and the DNA activator or inhibitor strands that set the switch's threshold; reasonable approximations are and . Apparent Hill exponents between 3 and 6 were experimentally measured in Kim et al (2006).

Depending on the parameters, two types of solutions are possible: the system may converge toward a stable steady state, or the steady state may become unstable, leading to sustained limit cycle oscillations. To simplify analysis of the design space for this oscillator, we first non-dimensionalize the equations exactly:

where x=[rA1]/K is the ratio of RNA activator to its threshold, y=[rI2]/K is the ratio of RNA inhibitor to its threshold, u=[T12A2]/[T12tot] is the fraction of Sw12 in the ON state, v=[T21A1]/[T21tot] is the fraction of Sw21 in the ON state, s=t/τ rescales time; in addition to these non-dimensional dynamic variables, there are new non-dimensional parameters , and . Up to linear scaling, system behavior depends only on α, β, γ, n, and m.

To visualize the oscillator behavior, we project onto the x–y plane. (We did not use the u–v plane, which would be more directly comparable to experimental fluorescence data, because u and v quickly saturate as x and y move away from their respective thresholds.) The steady-state solution can then be found easily by plotting the nullclines dx/ds=0 and dy/ds=0 with u and v determined by solving du/ds=0 and dv/ds=0. Figure 2, (a) and (b), show nullclines, steady state, and sample trajectories for parameters that yield oscillating and non-oscillating behaviors, respectively. The stability of a steady state can be determined by computing eigenvalues of linearized reaction equations—and in this case, an unstable steady state leads to limit cycle oscillations (see Supplementary information section 1.1). Interestingly, if γ≫1 or γ≪1, i.e., hybridization is either much faster or much slower than RNase H degradation, then the system is essentially two-dimensional and cannot oscillate (see Supplementary information section 1.1 and Novák and Tyson, 2008). Thus, experimentally we must aim for γ≈1, and further treatment of the simple model will assume this. We are left with a four-dimensional parameter space; the oscillatory regime can be identified numerically and visualized in two-dimensional cross-sections, e.g., assuming α=β and n=m (Figure 2, center, and Supplementary information section 1.1).

Figure 2

Simple model characterization and experimental results of the two-switch negative-feedback oscillator (Design I). The phase diagram of dynamic behaviors shows two domains with respect to changes in α and β (unitless parameters proportional to [T12tot] and [T21tot]) versus n and m (apparent Hill exponents of the activation and inhibition functions). (a, b) Phase portraits using unitless dynamic variables proportional to [rA1] and [rI2], x and y. The nullclines for x (violet lines) and y (black lines) are superimposed with temporal trajectories (cyan lines). There are fixed points where x and y nullclines intersect; stable (gray circle) and unstable (red circle) fixed points are marked. The system exhibits (a) limit cycle oscillation or (b) damped oscillation to a monostable steady state. Parameters are [T12tot]=[T21tot]=[A1tot]=100 nM, K=K=1 μM, n=m=5, τ=500 s, γ=1, k=0.002/s, and respectively k=0.04/s or k=0.025/s. In experiments, normalized fluorescence time courses measured the percent OFF-state switch Sw12 (TAMRA-labeled T12, green) and the percent OFF-state switch Sw21 (Texas Red-labeled T21, red). Two examples are shown: reaction #8 with stable oscillations and reaction #13 with strongly damped oscillations. Experimental parameters were, respectively, ([T21tot], [A1tot], [dI1tot],[T12tot], [A2tot], [RNAPtot], [RNaseHtot])=(250, 250, 700, 120, 350, 92, 8.3) nM and (250, 250, 1000, 80, 500, 125, 15.0) nM. Reaction #13 used higher [dI1tot], higher [A2tot], and lower [T12tot], compared with reaction #8. Thus, reaction #13 is mapped to a higher n and a lower α than reaction #8, and lay within the monostable domain rather than the oscillating domain. Experimental results (Supplementary Figure S4) were mapped to the phase diagram with respect to a reference oscillation, as discussed in Supplementary information section 1.1, and shown as ‘stable' (circles, damping coefficient <0.15/h), ‘damped' (dots, damping coefficient between 0.15/h and 0.5/h), and ‘strongly damped or too slow to measure' (crosses, damping coefficient >0.5/h) oscillations.

While insightful for basic phenomena, this simple model neglects important details of the experimental system, such as asymmetries between switches and the complexities of enzyme and hybridization reactions, and therefore is insufficient for quantitative modeling or for selection of experimental conditions. Specifically, as reaction rate constants are not under our control (unless we redesign the molecules), we can only adjust system parameters by varying the concentrations of enzymes and DNA species; to identify experimental conditions that yield oscillations, we need an accurate model that works directly with these concentration parameters. Therefore, we developed a mathematical model, termed the ‘detailed model', derived from the reaction mechanisms shown in Figure 1, using the Michaelis—Menten model for enzyme reactions (see Supplementary information section 1.4).

Recognizing uncertainties both in the model and rate constants from previous studies (Kim et al, 2006; Subsoontorn et al, 2011), we selected experimental DNA and enzyme concentrations based on the robustness with which they elicited oscillatory behavior in the detailed model. To do so, we used a straightforward random sampling technique inspired by stochastic high-dimensional sensitivity analysis (Feng et al, 2004). In brief, we considered rate and concentration parameters distributed log-uniformly over a preassigned ‘reasonable' range derived from previous studies (Kim et al, 2006; Subsoontorn et al, 2011), then estimated the marginal probability of oscillation conditioned on a particular parameter value (e.g., [A2tot]=100 nM) by random sampling of other parameters and numerical simulation of the detailed model. This analysis revealed clear trends for some parameters—e.g., [A1tot] should not be too high, [dI1tot] must be high (the higher the better), whereas [A2tot], [RNAPtot], and [RNaseHtot] each have optimal values. (These trends make sense; e.g., the high activation threshold set by dI1 introduces a delay in the switching response, and consequently in the oscillator dynamics, which is known to help achieve stable oscillations in similar simple negative feedback circuits (Novák and Tyson, 2008; Stricker et al, 2008).) Initial experimental parameters were chosen by hand to roughly maximize their estimated marginal probability, while also ensuring that, taken together, the simulations predicted reasonably fast limit cycle oscillations (see Supplementary information section 1.8 for details).

Using DNA and enzyme concentrations predicted to produce oscillations, experiments were run at 37°C and monitored real-time by fluorescence: the ON-state switches have low fluorescence because of quenching (Marras et al, 2002), whereas the OFF-state switches have high fluorescence (Figure 1A). Early experiments produced damped oscillations (e.g., Figure 2, bottom right, and Supplementary Figure S4) and further empirical optimization identified concentration parameters that gave multiple complete oscillation cycles (e.g., Figure 2, top right, and Supplementary Figure S4).

To examine our understanding of the phase diagram of the system, we mapped the final 37 experiments (Supplementary Figure S4 and Table S5) to the non-dimensionalized parameters for the simple model and plotted them with an indication of how strongly damped they were (see Supplementary information section 3 for calculation of damping coefficients). Because of oversimplifications in the simple model, rather than using the formulas to map experimental conditions to non-dimensionalized parameters directly, we assumed only that their overall scaling laws may hold, but the linear coefficients would have to be adjusted empirically. As an easy way to do this, our empirical mapping placed a reference experiment in the oscillating regime and then placed the other experiments relative to it on the basis of the formula's scaling predictions (see Supplementary information section 1.1 for details). Therefore, the exact positions of experimental data points are somewhat arbitrary. Nevertheless, the qualitative agreement between model and experiment over a considerable range of conditions is encouraging.

The fully optimized system revealed five complete oscillation cycles with a nearly 50% amplitude swing (Figure 3A) until, after roughly 20 h, the production rate could no longer be sustained. (The limited lifetime of the batch reaction could be due to, e.g., exhaustion of NTP fuel or buffer capacity, build-up of waste products, degradation of enzyme functionality, and so on.) Aliquots were taken from this reaction and run on a denaturing gel to measure RNA concentrations, and on a non-denaturing gel to obtain an independent measure of OFF switch concentrations (Figure 3B and C, Supplementary Figures S5 and S6). T12 showed little temporal variation in either fluorescence or gel measurements, whereas T21 showed consistent (but slowing) oscillation in both measurements. Consistent with the high percentage of Sw12 in the OFF state, RNA activator rA1 levels were not detectable. However, to our surprise, rather than oscillations with constant amplitude and constant mean as predicted by the detailed model, the denaturing gel measurements revealed that the RNA inhibitor rI2 concentration builds up after each cycle.

Figure 3

Experimental characterization of a two-switch negative feedback oscillator (Design I, reaction #37). (A) The fluorescence time courses report the OFF-state switch Sw12 (TAMRA-labeled T12, green) and the OFF-state switch Sw21 (Texas Red-labeled T21, red); corresponding extended model fits are shown as lines in lighter shades. The non-denaturing gel measurement of the OFF-state Sw21 (black squares, from C below) showed reasonable agreement with fluorescence results. (B) The total rI2 concentration measurement (blue circles, from C below) showed five complete oscillation cycles; extended model fits are shown as lines (light blue). The concentration of incomplete degradation products was estimated from the band of ≈35 nucleotides in the denaturing gel (magenta). The rA1 concentration was not measured because most bands were not significantly above background. (C) The concentrations of RNA signals and incomplete degradation products were measured in a denaturing gel (top and middle) and the OFF-state switch concentrations were measured in a non-denaturing gel (bottom). (See Supplementary Figures S5 and S6 for complete gels and Methods for normalization procedures.)

One hypothesis is that the short fragments of rI2 generated by RNase H degradation, which also build up roughly linearly over time (Figure 3B and C, and Supplementary Figure S5), may interfere with proper hybridization of rI2 signals to their regulatory target, Sw12. These predicted degradation products would encompass the 8-base toehold binding sequence of rI2 because RNase H cannot process several bases on the 5′ side of the RNA strand on an RNA/DNA hybrid substrate (Lima and Crooke, 1997). Thus, the short fragment can block the (otherwise freely available) toehold region essential for fast kinetics, effectively requiring a higher signal concentration to achieve equivalent kinetics for the inhibition reaction. Intriguingly, rather than destroying the oscillations, the build-up of interfering waste products merely elicits a compensatory shift in mean RNA concentrations and a slight slowdown of the oscillations. To show that rI2 dynamics can be explained by the process described above, we augmented the detailed model to take account of short degradation products; this ‘extended model' was then used for fitting model parameters to results from all 37 experiments simultaneously (as described in Supplementary information section 1.7 and Supplementary Tables S1 and S2).

While the extended model did not capture the experimental dynamics exactly, it was consistent with qualitative features of the experiments (e.g., baselines, amplitude, frequency, and damping) for most experimental conditions (Supplementary Figures S4 and S7). For example, both the experiments and the extended model showed a trade-off, preventing oscillations with a simultaneous high frequency and high amplitude. The dynamics in the absence of interference could be inferred by reducing the hybridization rate of short degradation products in simulation (Supplementary Figure S8), indicating that the interference from incomplete degradation product introduces further delay in oscillator dynamics and qualitatively changes its behavior. Furthermore, in experiments starting from four distinct RNA input combinations that initiate the switch states at various locations in the phase plane, the trajectories converged toward a damped limit cycle oscillation, the mean level of which slowly drifts; extended model simulations suggest that, while the damping can be explained by degradation waste interfering with switching, the drift is consistent with interference with the fluorescence readout mechanism (Supplementary Figure S9).

Despite the introduction of short degradation products, the extended model fits could not capture quantitative aspects for many reaction conditions. Poor fits can be expected for the initial part of fluorescence measurement because of a ‘burst phase' in enzyme kinetics (Jia and Patel, 1997) and for the final part of fluorescence measurement because of buffer exhaustion, NTP depletion, build-up of waste products, product inhibition (Arnold et al, 2001), and degradation of enzyme functionality—none of which were modeled. Rigorous mechanistic modeling that explicitly considers these factors may improve quantitative accuracy for our in vitro transcription systems. However, a further difficulty is that experimental results are sensitive to uncharacterized differences between enzyme batches; we were forced to fit some model parameters independently for the two enzyme batches used in the final 37 experiments with the Design I oscillator. Although there were quantitative differences, comparable oscillations were observed using both batches; we include both data sets here to illustrate expected experimental variation and to acknowledge difficulties with mechanistic modeling. Although a single enzyme batch was used for the Design II and Design III oscillators, we must keep in mind that model parameters are empirical and likely conflate many unmodeled effects.

Design II: an amplified negative feedback oscillator

Understanding systems-level behaviors in biochemical circuits relies heavily on comparative studies. The flexible architecture of our synthetic transcriptional network allows us to synthesize and characterize alternative oscillator designs. As an example, consider that incorporating a positive-feedback loop in addition to negative feedback (Tsai et al, 2008) has been identified as an important element for robust oscillation (Barkai and Leibler, 2000) in the cell cycle (Pomerening et al, 2005) and in engineered circuits in vivo (Atkinson et al, 2003; Stricker et al, 2008). To explore this principle in our in vitro system, we added a positive-feedback loop to the two-node oscillator. In addition to Sw21 and Sw12, the new autoregulatory switch Sw11 allows the RNA activator rA1 to promote its own production (Figure 4, top center). This modification requires only a single additional DNA strand, as the bottom strand of T11 is identical to that of T12, bringing the total to eight (Supplementary Figure S2). Nevertheless, the effect on system dynamics can be profound. Exploring qualitative behaviors using the ‘simple model' modified to include Sw11, again we calculated a phase diagram in terms of non-dimensionalized parameters (Figure 4, bottom center, and Supplementary information section 1.2). When the parameter proportional to [T11tot], δ, was zero, the system is identical to a Design I oscillator, and only two phases (oscillating and monostable) are observed. However, for higher δ, several new qualitative behaviors appear, including a bistable regime, a new ‘high' monostable regime, and three oscillator variants near the intersection of the four main regimes—e.g., a regime with both a stable steady state and a limit cycle attractor (Figure 4(a)–(f) and Supplementary information section 1.2). Near this intersection area, but still within the main regimes, behaviors similar to an excitable medium (Figure 4(d)) or to a relaxation oscillator (Figure 4(e)) can be found. Although the parameter choice used in Figure 4(e) did not give rise to characteristic relaxation oscillator behavior, the nullclines exhibit the characteristic cubic nonlinearity that leads to conditional bistability of one dynamic variable and gives rise to a hysteresis loop (Hasty et al, 2001; Pomerening et al, 2003). More clear relaxation oscillator behavior was found for slightly higher α, β, and n (Supplementary information section 1.2).

Figure 4

Simple model characterization and experimental results for the amplified negative feedback oscillator (Design II). A detailed reaction diagram is shown in Supplementary Figure S2. The phase diagram shows several domains of different dynamics with respect to changes in α and δ (unitless parameters proportional to [T12tot] and [T11tot]). (a–f) Phase portraits using unitless dynamic variables proportional to [rA1] and [rI2], x and y. The nullclines for x (violet lines) and y (black lines) are superimposed with temporal trajectories (cyan lines). Stable (gray circles) and unstable (red circles) fixed points are shown. The system exhibits several kinds of different dynamics: (a) damped oscillation to a low monostable steady state, (b) limit cycle oscillation, (c) damped oscillation to a high monostable steady state, (d) an ‘excitable' monostable steady state, (e) a ‘relaxation' oscillation, (f) bistable steady states, as well as more complex behaviors in the regions between (e) and (f). In experiments, normalized fluorescence time courses measured the percent OFF-state switch Sw12 (TAMRA-labeled T12, green) and the percent OFF-state switch Sw21 (Texas Red-labeled T21, red). The self-activating switch Sw11 was not labeled and hence not observed by fluorescence. The concentration of Sw11, [T11tot], was increased from 0 nM (reaction #5) to 30 nM (reaction #6) and 90 nM (reaction #8). Accordingly, a damped oscillation (reaction #5) became a stable oscillation (reaction #6) and then a much slower oscillation (reaction #8). Experimental reaction #8 maintained the Sw21 ON-state for the initial 5 h, consistent with the simple model prediction of monostable high behavior; yet it did not maintain an ON-state throughout the experiment, possibly because the extensive transcription caused early exhaustion of the buffer. Experimental results (Supplementary Figure S10) were mapped to the phase diagram as described in Supplementary information section 1.2 and shown as ‘stable' (circles, damping coefficient <0.15/h), or ‘strongly damped or too slow to measure' (crosses, damping coefficient >0.5/h) oscillations.

The monostable low and bistable regimes were previously demonstrated experimentally using just Sw11 (i.e., α=β=0) in Subsoontorn et al (2011). Here, we investigated the effect of Sw11 on oscillation. Consistent with the qualitative predictions of the simple model, starting from conditions that are strongly damped, the positive feedback element can reinforce the excitatory signal of rA1 such that sustained oscillation is achieved despite the weak activation from Sw12 because of its low concentration (Figure 4, right). Stronger positive feedback, achieved by increasing the concentration of Sw11, increased both the oscillation period and amplitude, until the system no longer appeared to oscillate, presumably entering the second monostable regime. Equations for the detailed and extended models for Design II, extended model fitting parameters, experimental conditions, and all eight experimental trajectories along with extended model fits can be found in Supplementary information sections 1.5 and 1.7, Supplementary Tables S3 and S6, and Figure S10.

While we did not experimentally demonstrate the full behavioral richness of Design II, this study illustrates the potential of in vitro transcriptional circuitry. In particular, it illustrates the ease with which the strengths of network connections can be tuned to explore a range of dynamic behaviors. Note, however, that because the thresholds are set within the signal blocks rather than within switch blocks, the threshold for self-activation by Sw11 must be the same as the threshold for Sw21. This and other intrinsic limits make it difficult to directly implement the exact dynamics of many previously studied positive/negative-feedback oscillators.

Design III: a three-switch ring oscillator

Our original goal for this project was to recreate in vitro a simplified version of the repressilator, a synthetic genetic regulatory network configured as a three-node ring oscillator operating within individual E. coli cells (Elowitz and Leibler, 2000). However, our original attempt suffered from unintentional secondary structure in some transcripts that degraded switch performance, prompting us to look for a simpler oscillator design. After our success with the Design I and Design II oscillators, we returned to the ring oscillator design.

Starting with the Design I system, we replaced the excitatory connection by a chain of two inhibitory connections, reusing components where possible and redesigning new DNA sequences to avoid undesirable secondary structure. Switch Sw12 is identical to the one in Designs I and II, but its transcript rA1 is rechristened rI1 and used in an inhibitory signal block. Transcript rI2 remains in an inhibitory signal block and a separately designed RNA inhibitor rI3 completes the cycle (Figure 5, inset). Using the modular domains within RNA signals, Sw31 and Sw23 were designed in a straightforward manner (Supplementary Figure S1). The resulting system has nine DNA strands (Supplementary Figure S3).

Figure 5

Experimental characterization of a three-switch ring oscillator (Design III, reaction #27). A detailed reaction diagram is shown in Supplementary Figure S3. (A) Fluorescence time course of the OFF-state switch Sw12 (TAMRA, green) showed four oscillation cycles; the extended model fit is shown as a line in lighter shades. (B) Gel data (circles) and extended model fits (lines) of RNA inhibitor concentrations (rI1: magenta, rI2: blue, rI3: brown) showed periodic changes in the correct order dictated by their regulatory connections. Note that the minima of rI2 (blue) correspond in time with the minima of T12 (green traces in A), as expected. (C) Denaturing gel measurement of RNA inhibitor concentrations. To account for some overlap between bands, the three inhibitor bands were fit as three Gaussian peaks with fixed center locations and widths (see Methods for details). The black lines are drawn as a guide for the eye (see Supplementary Figure S11 for complete gels and lane profiles).

To identify promising experimental conditions, again we developed a ‘detailed model' of the enzyme and hybridization reactions ignoring interfering degradation products (Supplementary information section 1.6) and used a random-sampling technique (Supplementary information section 1.8). Unlike the two-node oscillator model, strong sensitivity was not observed for any single DNA concentration parameter. Instead, to achieve sustained oscillations, the behavior of the three switches had to be approximately matched. Experimentally, however, Sw12 was slow to turn off in response to rI2 (data not shown), possibly because of fluorophore-quencher interaction stabilizing the ON state of the only fluorophore-labeled switch (Marras et al, 2002; Moreira et al, 2005). (In preliminary work, we observed mild to severe differences between fluorophore-labeled switches and their unlabeled twins, depending on the fluorophore used.) Balancing the strengths of the inhibitory connections was achieved instead by adjusting switch and activator concentrations.

After empirical optimization, the resulting oscillations (Figure 5A) were still slower than those of the two-node oscillator; therefore, typically only four cycles were observed before the batch reaction ran down. Although we monitored only one switch state by fluorescence, gel studies confirmed that all three RNA signals showed periodic changes in the order dictated by their regulatory connections (Figure 5B and C). Again, we developed an ‘extended model' that included the effects of degradation waste products (Supplementary information section 1.7) and fit the model parameters to 27 experimental trajectories, obtaining adequate fits for many experimental conditions (Figure 5, Supplementary Tables S4, S7 and Figure S12).

When initiated with no externally supplied RNA inhibitors (as in Figure 5), the first oscillation amplitude was small, but the following cycles showed increasing amplitudes and periods (Figure 5B). This ‘spiraling out' behavior was reproduced by the extended model; however, to our surprise we could not attribute the behavior to interference from waste products; the detailed model also exhibited the same qualitative behavior. Attributing the behavior instead to saturation of the degradation machinery, we developed a ‘simple model' incorporating this feature (Figure 6, top left, and Supplementary information section 1.3), which can exhibit monostable behavior, limit cycle oscillation (Figure 6 (a)), or spiraling out behavior (Figure 6 (b)) depending on parameters.

Figure 6

Simple model characterization of the three-switch ring oscillator (Design III). This model always has a single fixed point that can be calculated from the intersection of the production and degradation functions (upper left), here expressed for non-dimensionalized dynamic variables x=[rIi]/K and x=[rIj]/K, where (i, j)∈{(1,2),(2,3),(3,1)}. The phase diagram (lower left) shows several domains of different dynamic behaviors with respect to changes in α (the unitless parameter proportional to [Tijtot] and inversely proportional to v, the maximum degradation rate) versus n (the apparent Hill exponent of the inhibition functions). (a) For α=1, the maximum production rate of x is the same as the maximum degradation rate of x (top left) and the simulation generates a limit cycle oscillation (top center) with attracting phase plane dynamics (top right). (b) For α=2 with the degradation speed v halved, the maximum production rate of x is twice the maximum degradation rate of x (top left) and the simulation generates an oscillation with increasing amplitude (bottom center) and the phase plane trajectory spirals out (bottom right) because of saturation of the degradation machinery. This is similar to what we observed for reaction #27 as measured by gel (Figure 5 and Supplementary Figure S11). Experimental results (Supplementary Figure S12) are mapped to the phase diagram as described in Supplementary information section 1.3 and shown as ‘stable' (circles, damping coefficient <0.15/h), ‘damped' (dots, damping coefficient between 0.15/h and 0.5/h), and ‘strongly damped or too slow to measure' (crosses, damping coefficient >0.5/h) oscillations.

Specifically, the following simplified oscillator model was used:

where (i, j)∈{(1,2),(2,3),(3,1)}. Parameters for the simulation were chosen to be similar to those of the extended model of the Design III oscillator as follows: [Tijtot]=100 nM, k=0.01/s, K=0.333 μM, n=5, K=0.1 μM, v=1 nM/s (for a) or v=0.5 nM/s (for b). Note that and K represent the Michaelis constant of RNase H. No hybridization delay was modeled, unlike the simple models for Designs I and II, because it was unnecessary for oscillations. Further, saturation of production by RNAP was not modeled because it is not relevant to the behavioral distinction illustrated here. Incidentally, note that saturation of the production and degradation machinery also presumably occurs in the Design I and II oscillators, but it was omitted from the simple models for those systems in order to simplify the presentation. Further note that (unlike the detailed and extended models) this simple model treats each species independently and does not account for competition between substrates for saturated enzymes.

Again, we computed a phase diagram for the simple model (Figure 6, bottom left, and Supplementary information section 1.3) in terms of non-dimensionalized parameters , , and n. For β=0.3 and n>2, we found that system behavior depends almost exclusively on α, which can be considered an effective strength of transcription relative to degradation. When transcription is more than twice as strong as the maximal degradation, spiraling-out behavior is observed. An attempt as before to qualitatively map our 27 experimental trajectories onto the phase diagram revealed that our best (spiraling out) oscillations indeed had the highest α, but we did not observe a clear ‘limit cycle oscillation' regime between the ‘spiraling out' regime and the ‘monostable' regime.

With this, we conclude our experimental investigation of in vitro transcriptional oscillators.

Discussion

Proposing theoretical models of synthetic biochemical oscillators is easier than implementing them successfully in the laboratory. Indeed, a previous attempt to create an in vitro biochemical oscillator using reverse transcriptase and T7 RNA polymerase yielded excellent oscillatory behavior in theory (Ackermann et al, 1998) and every designed reaction step was shown to occur (Wlotzka and McCaskill, 1997). Yet, sustained oscillations were not reported, perhaps because of accumulated sequence mutations, as both DNA and RNA were synthesized and degraded. In contrast, in our transcriptional networks, only RNA is synthesized and degraded, whereas DNA molecules are preserved, thus grounding the dynamics. Furthermore, our transcriptional switches are capable of producing a threshold response with a high input/output gain and low background (Kim et al, 2006), which is essential for driving system dynamics from the stable regime into the oscillatory regime.

With the successful implementation of three oscillator designs here, together with previous implementations of bistable dynamics (Kim et al, 2006; Subsoontorn et al, 2011), transcriptional networks are beginning to live up to their promise (Kim et al, 2004) as a powerful architecture for systematic and modular construction of dynamic circuitry in vitro. The modularity and flexibility of our oscillator architecture were highlighted by the fact that a single additional DNA strand was required to transform the Design I to the Design II oscillator. A possible extension would be the addition of one more DNA strand for a switch Sw22 (which possesses the regulatory domain of Sw12 and the output domain of Sw21) that would add a self-regulating negative-feedback loop producing an analog to the in vivo oscillator reported by Stricker et al (2008) (cf., Smolen et al, 1998), the most robust synthetic oscillator demonstrated so far (Purcell et al, 2010). In principle, not only can transcriptional networks be wired up by sequence design to form arbitrary neural network or boolean logic circuits, but also—as illustrated here—even simple circuits can be tuned to elicit a wide range of behaviors simply by adjusting the concentrations of switch and threshold molecules.

However, several significant difficulties must be considered. First, to emphasize that all observed dynamics are due to chemistry and not due to the experimental apparatus, we have exclusively run closed batch reactions. Consequently, the reactions have a limited lifetime as sources of material and energy (such as NTPs) are used up, enzymes stop working, and waste products accumulate (cf., Klungsøyr et al, 1968). If these constraints were to be relaxed, e.g., by using a chemostat (Atkinson et al, 2003), a dialysis bag (Madin et al, 2000), or vesicles (Noireaux and Libchaber, 2004), we anticipate that sustained limit cycle oscillations could be observed for all three oscillator designs investigated here.

Furthermore, we must recognize that waste products, from incomplete degradation and abortive or promiscuous transcription, form a second independent difficulty. As seen here, waste products can interfere with intended reactions and dynamically change phenomenological rate constants, making systems behavior difficult to model and control. This ‘noise' in the form of unknown partial degradation and transcription products is quite different in character from the more commonly studied noise of stochastic reactions and uncertainties in rate constants. Whereas in vivo systems will dilute waste products during exponential growth phase, our in vitro systems can be limited by ribonuclease saturation, complicating the elimination of waste products. Although identifying more effective degradation pathways would be of considerable help, the deeper issue is to understand how to design circuits whose behavior is robust to interference from a spectrum of unknown molecular species. We saw hints of such robustness in the difference between the Design I oscillator, which responded to accumulation of waste by a compensatory increase in RNA signal concentration while maintaining oscillations, and the Design III oscillator, which responded with damped oscillations unless it was in the ‘spiraling out' regime.

The combination of waste products and other unknown species and interactions creates a third difficulty, that of predictive modeling. This problem is particularly apparent in the asymmetry of switch components. For example, in the Design I oscillator, the simple model assumes symmetric switch characteristics and produces comparable-amplitude oscillation of both switches, as does the detailed model when provided with symmetric parameters. Random sampling of the detailed model using symmetric parameter ranges suggested initial experimental conditions that did, in fact, lead to relatively large amplitudes for both switches, although only as damped oscillations (e.g., Supplementary Figure S4, reaction #1). Nevertheless, asymmetric switch behaviors became more pronounced while optimizing experimental conditions to achieve sustained oscillation. This could be attributable to the influence of various waste species, to unintended secondary structure in designed single-stranded species, to unintended binding between strands, or to interactions between fluorophores, quenchers, and DNA (Marras et al, 2002; Moreira et al, 2005), any of which could differentially slow down only certain hybridization or strand displacement reactions. Such asymmetric behaviors can be captured in the detailed and extended models a posteriori by the phenomenological fitting parameters; yet they make predictive modeling and design difficult.

These difficulties go against the original spirit of this work, which was to simplify the in vivo repressilator by constructing an analogous in vitro oscillator that could be more rigorously characterized and modeled, thanks to the knowledge of exactly what went into the test tube and to the absence of a myriad of unknown species and reactions within the living cell. Yet, in our hands, the in vitro systems were still too complex for satisfactory quantitative and predictive modeling: ‘simple models' were useful for theoretically characterizing the phase diagrams of possible behaviors, but parameters could only be roughly mapped to experimental conditions; ‘detailed models' were effective for identifying promising experimental conditions and gave us confidence in the major reaction pathways but failed to account for several striking features observed experimentally; and ‘extended models' could qualitatively account for these features but still failed to quantitatively match some experimental traces even when fit a posteriori. With up to 18 dynamic variables and 33 system parameters, the latter models inevitably yield fits with several ill-constrained and phenomenological parameters, implying that sophisticated Bayesian methods would be required to make effective predictions (Brown and Sethna, 2003; Gutenkunst et al, 2007). Thus, our experience is somewhat in contrast to the optimistic perspectives of Rosenfeld et al (2007), Canton et al (2008), and Lucks et al (2008) for predictive synthetic biology in vivo—perhaps because the in vitro setting allows us to more easily appreciate the scope of the difficulties, or perhaps because living cells are effective at creating an insulated and reliable local environment for circuit function.

In any case, it is natural to expect that in vitro and in vivo biochemical networks are governed by many similar design principles. Our in vitro oscillators exhibit several design principles previously observed in vivo. (1) Introducing delay in a simple negative-feedback loop can help achieve stable oscillation (Novák and Tyson, 2008; Stricker et al, 2008), as obtained by high activation and inhibition thresholds in the Design I oscillator. (2) The addition of a positive-feedback self-loop to a negative-feedback oscillator provides access to rich dynamics and improved tunability (Tsai et al, 2008), as illustrated in the simple model phase diagram for the Design II oscillator. (3) Oscillations in biochemical ring oscillators (such as the repressilator) are sensitive to parameter asymmetry among individual components (Tuttle et al, 2005), which crippled our original design for the three-switch ring oscillator (Design III) and necessitated adjustments of switch and threshold concentrations even with the redesigned sequences reported here. (4) The saturation of degradation machinery and the management of waste products can have a surprisingly large role for establishing logical regulation at steady state (Guet et al, 2002; Kim and Tidor, 2003) and for improving the robustness of oscillators (Fung et al, 2005; Wong et al, 2007). (It is possible that these issues could be part of the explanation for the lack of robustness and the arrest of oscillation upon entry of E. coli into stationary phase for the repressilator (Elowitz and Leibler, 2000).)

All this begs the question of whether even simpler or more robust biochemical oscillators can be designed for in vitro studies. In biology, minimal circuitry required for oscillation can be quite simple, as exemplified by the circadian clock of cyanobacteria that requires only three proteins and ATP fuel to function when reconstituted in vitro (Nakajima et al, 2005). Indeed, recent work has shown the de novo design and implementation of an in vitro biochemical oscillator even simpler than any of ours (Montagne et al, 2011). However, an alternative view—which we take—is that increasing complexity is inevitable if we wish to create and explore increasingly interesting and powerful information-based chemical systems; hence, we will have to learn how to deal with it one way or another, and cell-free in vitro systems offer a valuable training ground (Simpson, 2006).

Natural next steps would be to improve our characterization and modeling of our oscillators by independent measurements of each reaction (Ronen et al, 2002); to implement more robust oscillator designs (Purcell et al, 2010); to couple them to other chemical processes such as DNA nanomachines (Dittmer and Simmel, 2004) and thereby control temporal self-organization; to prepare them as two-dimensional reaction-diffusion media (Grzybowski, 2009) and thereby control spatial self-organization; to express them within cell-like volumes to examine stochastic behavior (Shahrezaei and Swain, 2008); and to integrate them with other transcriptional circuitry to provide embedded controllers within prototype artificial cells (Szostak et al, 2001) such as water-in-oil emulsions (Griffiths and Tawfik, 2006) and vesicles (Noireaux and Libchaber, 2004).

Materials and methods

DNA oligonucleotides and enzymes

The sequence of all DNA molecules and expected RNA transcript sequences were chosen to minimize the occurrence of alternative secondary structures, checked by the Vienna group's DNA and RNA folding program (Flamm et al, 2000). Sequences can be found in Supplementary information section 4 and Supplementary Figure S1. All DNA oligonucleotides were synthesized and PAGE or HPLC purified by Integrated DNA Technologies. T21-nt was labeled with a Texas Red fluorophore at the 5′ end, T12-nt was labeled with a TAMRA fluorophore at the 5′ end, A1 and A2 were labeled with Iowa Black-RQ quenchers at the 3′ end. The T7 RNA polymerase (enzyme mix), transcription buffer, and NTP as part of the T7 Megashortscript kit and E. coli RNase H were purchased from Ambion. Note that according to the manufacturer's patent (# 5,256,555), the enzyme mix contains pyrophosphatase to extend the lifetime of the transcription reaction; as pyrophosphatase is involved in regulating the power supply for our transcriptional circuits and is not directly involved in the dynamics, we neglect this enzyme in our models and do not call it an ‘essential enzyme' for the circuit dynamics.

Transcription

Switch templates (T-nt and T-t strands) were annealed with 10% (v/v) 10 × transcription buffer from 90 to 20°C over 1 h at five times the final concentrations used. To the annealed templates, we added 15% (v/v) 75 mM each NTP (1.5 times the suggested amount of the kit), 8% (v/v) 10 × transcription buffer, 5% (v/v) 300 mM MgCl2, and water to fill up the volume. The NTP concentration was increased in an attempt to extend the lifetime of the batch reaction and the extra Mg2+ was added to balance the negatively charged NTPs. Transcription reactions for spectrofluorometer experiments were carried out at a total volume of 60 μl in a quartz cuvette. Once fluorescence signals were collected for 5 min, the DNA activators (A1 and A2) from high concentration stocks (50 μM) were added and mixed in the cuvette. After 10 min, the DNA inhibitor dI1 was added from 50 μM stock and mixed. After 10 min, purified RNA signals were added from 30 μM stocks and mixed, for those reaction conditions that used initial RNA signals. Fifteen minutes after adding all ssDNA strands and RNA signals, 4–11% (v/v) RNAP was added and mixed, 1 min after which 0.35–2.1% (v/v) RNase H was added and mixed. The fluorescence data used in the text and Supplementary information reveal the fluorescence traces after enzyme addition is complete. The nominal concentrations quoted by the manufacturer (RNAP stock=1.25 μM and RNase H stock=1.25 μM) were used as the enzyme stock concentrations. The reaction conditions and estimated final enzyme concentrations are listed in Supplementary Tables S5 to S7. Samples were taken at 30 min intervals (or as indicated in Figures 3 and 5) for gel studies and the transcription reactions were stopped by the addition of denaturing dye (80% formamide, 10 mM EDTA, 0.01 g XCFF).

In an attempt to obtain better oscillations by avoiding the large initial transient, and to explore system dynamics from different initial conditions, some experiments were started with RNA present before the addition of enzymes. For the purification of RNA signals for those experiments and for gel controls, full-length template side strands (the complement of T-nt rather than T-t) were used to prepare fully duplex DNA templates for rA1 (which is rI1 in the Design III oscillator) and rI2. The transcription reactions were prepared as a total volume of 60 μl with 0.2 μM fully duplex DNA templates. The transcription conditions were the same as above, except that 20% (v/v) RNAP was used and RNase H was omitted. After a 6 h incubation at 37°C, the reaction mixture was treated with 2.5 μl DNase I for 30 min to remove DNA templates, and stopped by the addition of denaturing dye. The reaction mixture was run on 8% denaturing gel; RNA bands were excised and eluted from gel by the crush-and-soak method, ethanol precipitated, and resuspended in water.

Data acquisition

For spectrofluorometer experiments, excitation and emission for TAMRA-labeled T12 were at 559 nm and 580 nm, whereas excitation and emission for Texas Red-labeled T21 were at 597 nm and 615 nm. Fluorescence was recorded every minute using a SPEX Fluorolog-3 spectrofluorometer (Jobin Yvon) and converted to switch activity (percent of switch in the OFF state) by normalizing against minimum fluorescence (measured before the addition of enzymes with excess quencher-labeled activators) and maximum fluorescence (measured at the end of reaction after addition of excess DNA inhibitors to displace activators).

Denaturing polyacrylamide gels (8% 19:1 acrylamide:bis and 7 M urea in TBE buffer) were allowed to run for 50 min with 10 V/cm at 65°C in TBE buffer (100 mM Tris, 90 mM boric acid, 1 mM EDTA). The non-denaturing gels (10% 19:1 acrylamide/bis in TAE/Mg buffer) were allowed to run for 100 min with 13 V/cm at 35°C in TAE buffer containing 12.5 mM Mg2+ (40 mM Tris-Acetate, 1 mM EDTA, 12.5 mM Mg-Acetate, pH 8.0). The 10-base dsDNA ladder (Invitrogen) was used in the marker lane and the gels were stained with SYBR gold (Molecular Probes) for quantitation. The gel data were collected on a Molecular imager FX gel scanner (BioRad) and measured using the rectangular box tool and the lane profile tool of the Quantity One software package (BioRad). The rectangular box tool counts total fluorescence within a box drawn around the gel band of interest, with automatic subtraction of mean band densities of a background box, whereas the lane profile tool counts fluorescence at each pixel on a lane. For Design I, a background correction was performed with respect to regions between T21-nt (T23-nt for Design III) and T12-t bands, except for RNAs in control lanes, for which a background correction was performed within the control lanes. Background correction was performed for each gel separately.

The rI2 concentrations in denaturing gels for the Design I oscillator were measured using the rectangular box tool, normalizing with respect to rI2 (62 bases) band densities in control lanes that contained 1 μM of purified rA1 (67 bases) and rI2 (62 bases). In many cases, the rA1 (67 bases) bands were not significantly above background or were obscured by the tails of rI2 bands when analyzed using the lane profile tool. Quantifying the rA1 concentrations in the presence of high concentrations of rI2 would lead to an overestimation of rA1 concentrations because of significant contributions from the tails of rI2 bands; therefore, we chose not to quantify the concentrations of rA1 bands. Waste w35 (∼35 bases) was normalized to the rI2 signal level in control lanes; because staining is roughly linear in length, this underestimates the concentration, and because there are many other waste lengths that we did not measure we report waste concentrations in ‘arbitrary units' (a.u.; Supplementary Figure S5). The OFF-state switch concentrations in non-denaturing gels were measured using the rectangular box tool with respect to band densities in control lanes containing 1 μM each of OFF-state switches T12 and T21 (Supplementary Figure S6).

The RNA concentrations in denaturing gels for the Design III oscillator required more careful methods for measurement because the three bands overlapped somewhat. As shown in Supplementary Figure S11, alignment of lane profiles based on three template strands (T12-t, T31-t, and T23-t) can assign the RNA peaks to rI1, rI2, and rI3. Because rI3 (64 bases) migrates close to rI2 (62 bases) and cannot be unambiguously distinguished using the rectangular box tool, we chose to model the three RNA bands within lane profiles as Gaussian peaks and fit their heights as a measure of RNA concentrations. The 1 μM of rI1 and rI2 bands in control lanes were well described by two five-pixel-wide Gaussian peaks. Using the alignment results to assign the center locations of rI1 and rI3 to be 15 pixels apart and those of rI1 and rI2 to be 20 pixels apart, the heights of all three Gaussian peaks were simultaneously fit by custom MATLAB code using the fminsearch function. For normalization, the fitted heights of rI1 and rI2 in control lanes with the same center locations and widths were used. Because staining is roughly linear in length, SYBR gold fluorescence from rI3 was assumed to be the mean of those from rI1 and rI2.

Model simulation

The kinetic simulations and parameter fittings were implemented in MATLAB. Differential equations were solved using the ode23s function. For parameter fitting, we used a cost function using least-squared errors of fluorescence trajectories, gel results, and characteristics of oscillation (oscillation amplitude, frequency, and damping coefficient) between simulation results and experiments. The cost function was minimized using the fmincon function. Details of the experimental parameters and model equations are in Supplementary information. MATLAB and SBML files are available as Supplementary information. SBML files have been submitted to the BioModels database (http://www.ebi.ac.uk/biomodels/) with accession numbers MODEL1012090000 to MODEL1012090006.