Synchrony and pattern formation of coupled genetic oscillators on a chip of artificial cells.

Understanding how biochemical networks lead to large-scale nonequilibrium self-organization and pattern formation in life is a major challenge, with important implications for the design of programmable synthetic systems. Here, we assembled cell-free genetic oscillators in a spatially distributed system of on-chip DNA compartments as artificial cells, and measured reaction-diffusion dynamics at the single-cell level up to the multicell scale. Using a cell-free gene network we programmed molecular interactions that control the frequency of oscillations, population variability, and dynamical stability. We observed frequency entrainment, synchronized oscillatory reactions and pattern formation in space, as manifestation of collective behavior. The transition to synchrony occurs as the local coupling between compartments strengthens. Spatiotemporal oscillations are induced either by a concentration gradient of a diffusible signal, or by spontaneous symmetry breaking close to a transition from oscillatory to nonoscillatory dynamics. This work offers design principles for programmable biochemical reactions with potential applications to autonomous sensing, distributed computing, and biomedical diagnostics.

Synchrony and pattern formation are manifestation of nonlinear reaction dynamics in discrete or continuous systems (1, 2). A population of independent oscillators reduces to a globally synchronized oscillation when the coupling between them is strong (3). Pattern formation can result from spontaneous symmetry breaking as in Turing patterns (4) and Belousov–Zhabotinsky reaction (5), or induced by spatially varying signals in morphogenesis (6), and wave-front coupled to gene-expression oscillations (7). Inanimate closed chemical systems exhibit collective modes transiently toward a spatially homogenous chemical equilibrium (8). In living systems, however, spatial self-organization stems from nonequilibrium internal cellular dynamics of biochemical networks, combined with molecular interactions between cells, all of which are difficult to isolate and control (9).

Synthetic gene networks have recently been engineered to program and reconstitute oscillatory behavior in single cells (10, 11), as well as synchrony and pattern formation in populations (12–14). Cell-free systems provide another level of simplicity and control, offering a means to design reactions and overcome the inherent entanglement of processes in living systems. Recent examples of oscillatory dynamics include gene-expression (15–17) and transcription-only (18) networks, as well as purified enzyme networks (19, 20). First steps toward spatial patterns in cell-free systems were demonstrated in gels (21), protein surface waves (22), DNA enzymatic reactions (23, 24), and gene-expression networks (25). Here, we assembled a chip of DNA compartments (15, 25) to program a one-dimensional system of up to 15 coupled oscillators driven by a gene-expression network, and revealed mechanisms leading to synchrony and pattern formation.

Results and Discussion

Compartmentalization and Design of Oscillations in Cell-Free Gene-Expression Reaction.

Gene constructs encoding an oscillatory network were immobilized as a DNA brush on the surface of a 2D compartment carved in silicon to a height of and radius , connected by a capillary of length and width to a flow channel feeding in a cell-free transcription–translation reaction based on Escherichia coli extract (15, 26) (Fig. 1 and ). The combination of the thin layer of compartments and capillaries with the deep and wide flow channel creates a scenario in which transport into and out of the compartment occurs solely by diffusion. The junction between channel and capillary fixes a zero-concentration boundary condition for newly synthesized molecules, thereby creating a source-sink dynamic with an effective lifetime of expressed proteins (15), , where we use a typical diffusion constant, . We designed the gene network based on general principles of biochemical oscillators, which include negative feedback, nonlinearity, and time delay (27). The network consisted of a nonlinear activator-repressor loop with E. coli transcription factor and lambda phage repressor (15) (Fig. 1 and ). To achieve oscillations in a wide parameter range we expressed constitutively two delay elements: an anti- inhibitor to sequester the activator (25), , and a protease complex, to degrade the repressor by targeting its ssrA degradation tag (28, 29). These high-affinity elements create a nonlinear threshold of their target activation at low concentrations, thereby acting as an effective delay of both repression and activation (30). The dynamics of the network was reported by EGFP regulated by either activator or repressor promoter.

Fig. 1.

Synchrony and pattern formation in an array of DNA compartments. (A) Overlay image of expressed GFP (488 nm) and fluorescently labeled DNA patterns (white square, 647 nm) in a circular compartment carved in silicon, connected by a diffusive capillary to a feeding channel flowing a cell-free reaction mix. (Scale bar, ) (Network diagram) Activator–repressor network with activator and repressor , tagged with an ssrA degradation tag. The protease complex is synthesized and assembled in the compartment under a promoter and degrades the -ssrA protein, controlling the delay in repression. Finally, a promotor expresses binding to the , which sequesters its activity to control delay of the activator. The reporter gene is either under the regulation of or promoters. (B) Dynamics of 15 different isolated oscillators with varying gene composition. (C) Overlay image of fluorescently labeled DNA and of GFP expressed in three oscillators coupled in an array. Distance between compartments , compartment capillary length L = 200 μm, and , the capillary length between connecting capillary and the main flow channel. (D) Protein expression profile in an array of coupled DNA compartments originating from a single DNA source constitutively expressing GFP under promoter. Data are fitted to an exponential profile (solid line) with compartments. (Scale bar, ) (E) Dynamics of the 15 oscillators in a coupled array. (F) Space–time images of GFP in an array of identical oscillators with and without an activator source at the first compartment.

The network architecture provided flexibility in exploring the effect of gene ratios and composition on the dynamics, resulting in a wide range of stable and unstable oscillatory behavior, with period from 2.4 to 4.5 h and amplitude up to protein observed for 15–18 h (). Fig. 1 shows the dynamics of 15 different oscillators in separate compartments (). To couple the oscillators we used a chip in which compartments are connected laterally, apart, through thin capillaries, such that signals emanating from a compartment diffuse to neighboring ones, with concentrations decaying exponentially (25), , with decay length to 5 compartments (Fig. 1 ). Strikingly, when these 15 oscillators (Fig. 1) were coupled by diffusive transport of newly synthesized proteins, their frequency and phase synchronized (Fig. 1), thereby creating long-range order on a scale of the system size, mm. The hierarchy in amplitudes of the uncoupled and coupled oscillators was conserved.

The synchrony and long-range order can be reasoned by considering the classic Kuramoto model for oscillators that are all mutually coupled through their phases. Within this model there exists a critical coupling above which oscillators synchronize (3, 31), , with coupling strength and frequency variance . Experimentally, we varied the distance between the main feeding channel and the connecting capillary, at fixed (Fig. 1). This length scale controls the typical timescale for diffusion between compartments, . We therefore estimate the coupling strength as , which increases with for ; for . The discontinuous jump in occurs because the compartments are completely isolated by the main channel when . We note that the genetic oscillators are locally coupled; hence, long-range synchrony is less expected than Kuramoto oscillators. Nonetheless, we find that with coupling strength , and variations in the frequency, , implying that synchrony is consistent with the Kuramoto model.

Having demonstrated long-range order in the synchrony of a coupled array of different oscillators, we next sought to explore the emergence of pattern formation. Inspired by morphogenesis, in which identical cells respond to a concentration gradient and develop patterns of expression, we assembled an array of coupled identical oscillators subject to a symmetry-breaking signal. We immobilized in the first compartment the gene coding for the activator controlled by a constitutive promoter that is not influenced by the oscillatory network. The remaining 14 compartments were encoded by identical oscillators. The activator source diffused into the array of oscillators, locally increasing the concentration of activator along an exponentially decaying profile. The gradient of activator in the array disrupted the synchrony of the oscillators and created dynamic patterns, typically changing over a timescale of (Fig. 1). Most notably, we observed states with spatial oscillations in which neighboring compartments exhibited anticorrelated patterns that change in time, as measured for and (). While the oscillator phases varied in space, their period remained nearly constant independently of the distance from the gradient source (). Space–time plots of these dynamics show a checkerboard pattern under the action of the activator gradient, but a spatially uniform pattern of synchronized oscillators without the gradient (Fig. 1).

Properties of the Biochemical Oscillator at the Single-Compartment Level.

To further understand the nature of synchrony and pattern formation, we varied network parameters and investigated the dynamics of isolated oscillators. We measured the oscillator frequency and amplitude as a function of the activator concentration by changing the fraction of its gene in the DNA brush, , while keeping the rest of the components of the oscillator at constant stoichiometric ratio (Fig. 2 and ). Similarly, we varied the gene fractions of the repressor [R] and the protease complex beyond its basal level in the cell-free reaction (26). The oscillator frequency decreased continuously by a factor of 2.6 from for . Increasing [R], at a midrange activator , also resulted in a decrease of frequency, but with a weaker effect. In contrast, the frequency increased with by a factor of 1.5. Thus, enhancing the negative feedback either by increasing activator or repressor, or by decreasing the protease delay element, slows down the oscillator. These data are consistent with a numerical solution of the network dynamics model (). Notably, by removing the respective genes and degradation tags we verified that each delay element was sufficient and necessary to enable weak amplitude oscillations, while their combination resulted in stable pronounced oscillations (Fig. 2 and ). We conclude that negative feedback retards the oscillator frequency and that delay elements stabilize oscillations in parameter space.

Fig. 2.

Oscillations at a single-cell level. (A) Oscillator frequency as a function of the gene fraction of activator , repressor , and protease delay element in the brush. (B) Different oscillatory dynamics observed for combinations of both delay element, the inhibitor , and protease (Deg). Note: Degradation was eliminated from the circuit by removing the ssrA tag from the repressor. (C) Oscillations as a function of time with the activated gene (orange) and repressed gene (blue) as a reporter. (D) Distribution of period and width of the oscillations for the activated and repressed genes with . Each histogram contains 50 isolated oscillators. Variation in width was normalized separately to the mean of each peak.

The ability to replicate oscillators on the chip in isolated compartments enabled us to study the variability inherent to the network topology. We therefore assembled 50 oscillator replicas, and measured separately the repressor and activator expression dynamics (Fig. 2). Each oscillator dynamics is characterized by a period and a peak width . The distribution of the activated promoter exhibited variability in with an SD, , and similarly for , with (Fig. 2 and ). For the repressed promoter we find , and a narrower distribution of the peak width, . Since and , we conclude that period variability primarily stems from the activated gene. The reduced variability in the peak width of repressed promoter is consistent with previous observations showing that negative feedback reduces noise in gene circuits (32).

Entrainment of Oscillations in a Pair of Coupled Compartments.

The coupling of nonlinear oscillators leads to frequency selection that is less predictable than linear oscillators, in which the frequencies are often linear combinations of the natural ones. Cases where slow oscillator dominates the dynamics are less prevalent (33), as observed in certain cases of circadian clocks (34, 35). This motivated us to study the synchrony of a pair of nonlinear genetic oscillators by changing the coupling strength and network parameters (Fig. 3). We designed a dual-compartment geometry to create coupling with minimal asymmetry due to residual pressure difference. The two compartments were connected to a single point by a widened entry capillary at the junction of the feeding channel. The coupling strength was varied by changing the position of an auxiliary capillary, such that the amplitude of a single oscillator was reduced in the adjacent compartment, where corresponds to 20% of the source value, corresponds to 40%, and corresponds to 100% (). The dynamics of each oscillator A in a pair was measured when its coupled oscillator B was either present or absent from the neighboring compartment.

Fig. 3.

Entrainment and synchrony in coupled compartments. (A) Overlay image of two oscillator gene networks patterned in a coupled pair of compartments. (Scale bar, .) (B) Pairs of coupled oscillators in their different configurations, A and B: (i) uncoupled–defined by a natural period difference ; (ii) coupled–synchronized; coupled and uncoupled with a period difference (iii) ; (iv) . (C) as a function of , measured for three coupling length s values as denoted. (D) Array of 10 different oscillators (A–J) patterned in 15 compartments interconnected by a diffusive capillary of and varying . (Scale bar, .) (E) Space–time plot of oscillators A–J at different coupling strength. Blue (green) color represents low (high) protein concentration in arbitrary units. (F) Synchrony measure of coupled oscillators as a function of geometry, as defined in . (G) Spatial correlations of protein concentration between oscillators separated by a distance averaged over time and space. Correlations are measured up to a distance of smaller than array size. (H) Fitted correlation length as a function of .

In Fig. 3 we show a pair of different oscillators measured both coupled and uncoupled, with a natural period difference of . We define the effect of coupling as the difference between the natural period of the oscillator and the period when it is coupled, . We find that the coupled pair synchronized with the slower oscillator entraining the faster one. The selected pair period was identical to the natural slower period within error of (). Remarkably, this result holds for every one of 10 different oscillators measured at three different couplings, in a range of natural period difference, . Whenever oscillator A was slower, it maintained its natural period, , and entrained oscillator B, , and vice versa (Fig. 3). Because the network is an effective negative-feedback loop, coupling of two oscillators is inhibitory; thereby, all interactions are expected to slow down the dynamics, just as the negative feedback slows down the period of an isolated oscillator (Fig. 2). Furthermore, we observed that at low frequencies the oscillator is characterized by high repressor amplitudes (34, 36), which could further explain the dominance of the slow oscillator (). The entrainment of two coupled oscillators is captured in the numerical model with the slow one dominating the dynamics, yet the slow oscillator entrains the fast one to a period that is up to 10 min different from its natural period ().

Synchrony of Oscillators Set by Intercompartment Geometry.

To address the question of just how synchrony emerges in an ensemble of coupled oscillators (Fig. 1), we varied the coupling strength between neighboring compartments and measured the collective dynamics in an array of 15 compartments. Array size was chosen to be bigger than the largest decay length in the system compartments. The coupling was varied using a capillary connecting the compartments, whose distance from the feeding channel sets the concentration gradient (Fig. 3): For large the decay length was maximal and compartments are strongly coupled, whereas for compartments were isolated. Space–time plots of the dynamics show a gradual transition to synchrony as the coupling parameter increases (Fig. 3 and ). The degree of synchrony was defined as the normalized time-averaged fluctuations of the concentrations of all oscillators (37) (): varying from a random ensemble, to perfect synchronization, . As expected, synchrony in the coupled array increased gradually with coupling strength reaching a highly synchronized state, for (Fig. 3). To further characterize the synchrony in the ensemble we computed the spatial correlations of protein concentration between every two oscillators separated by a distance , averaged over time and compartment location, . We find that correlations decay exponentially, , with distance proportional to the geometrical coupling, S (Fig. 3 ). The decay of correlations is consistent with local coupling between neighboring oscillators, and a gradual transition to long-range order in the limit of very strong coupling.

Mechanisms for Pattern Formation in an Array of Locally Coupled Oscillators.

We next studied morphogen-induced patterns in the coupled array of identical oscillators, elaborating the results presented in Fig. 1. Without a gradient source, the dynamics was homogeneous in space, as reflected by straight lines in the space–time plots (Fig. 4 , I). Oscillators in this configuration exhibit a period variation of 10–15%, corresponding to their location along the array (). This variation is likely due to the boundaries of the array, which alter the lifetime and steady-state concentrations of the compartments at the edges, and to residual flow along the connecting capillaries, creating a small asymmetry in concentrations along the array. In the presence of a morphogenetic source, the spatial symmetry was broken, resulting in inhomogeneous expression dynamics. We used two independent signals to induce patterns: the activator and the inhibitor delay element (Fig. 4 , II and III and ). The morphogen signal was constitutively expressed throughout the duration of the experiment from a source located at the first compartment. With either signal, the expression dynamics of the remaining 14 compartments initiated by a single synchronized pulse, , followed by an intermediate state of constant expression levels up to . Consistent with the observation that high concentration of activator inhibits oscillatory dynamics (Fig. 2), we find that activator morphogen resulted in a prolonged period of low expression, depending on the distance from the source. Compartments close to the source light up later than farther ones, while those outside the gradient range oscillate from the start. This trend resulted in a front of low to high expression propagating toward the source, after which we observed antiphasing checkerboard patterns in the space–time plot.

Fig. 4.

Mechanisms for pattern formation in an array of coupled oscillators. Space–time plot of an array of 14 identical coupled oscillators with (A) no external “morphogen” gradient; (B) a morphogen source of activator protein ; (C) a morphogen source of an inhibitor delay element . Sources were located at the first compartment. Blue (green) color represents low (high) protein concentration. (D) Dynamics of two oscillators located at adjacent compartments along the array (I) with no source, (II) with an activator source, (III) with a delay element source. (E) Distribution of phase difference between adjacent couples of identical oscillators along an array with a source of for (red) before the gradient was established, and for (white). (F) Spatial correlations averaged over time between couples of adjacent compartments without a source and with an activator source. (G) Velocity of backward propagation measured for six different oscillators under the influence of a gradient of activator as a function of the oscillation period (compartment per hour). (H) Transition to nonoscillatory regimes measured in isolated compartments, and in coupled compartments. The transition occurs at . (I) Enlarged space–time plot of spontaneous pattern formation at the transition. Dynamics obtained with as reporter.

To further study loss of synchrony due to the activator gradient we analyzed the distribution of nearest-neighbor phase difference and time-average correlation in expression as a function of distance from morphogen, (). At early time, , the phase difference was sharply distributed, whereas for the distribution broadened, reflecting loss of synchrony (Fig. 4). Close to the source, neighboring oscillators were anticorrelated, with a transition to synchrony far from the source, compartments (Fig. 4). A similar set of coupled oscillators showed that backward propagation is a robust feature of the activator morphogen, yet checkerboard patterns do not always occupy the entire space–time plot (). Interestingly, the front propagation velocity scaled inversely with the natural oscillator period, (compartments per hour) (Fig. 4). The slowing down of the propagation velocity for longer periods is in line with the enhancement of negative feedback by the activator. In contrast, when replacing the activator by a source of the delay element , which negates the inhibition, we observed an inverted pattern with backward propagation of a transition from high to low expression. When the activator or sources were coupled to the network feedback, by placing them under regulation of the repressed promoter of the oscillator, the resulting patterns reverted to a nearly synchronized state, implying entrainment of the source ().

Finally, we addressed the question whether patterns could emerge without a morphogen source. A mechanism for symmetry breaking by fluctuations could be of importance in biological processes, for example, as suggested is cases of early development (38). We first varied the amount of activator gene fraction in isolated compartments, , arbitrarily choosing [XP] = 0, using only basal levels of endogenously found in the cell-free reaction. We observed a transition from nonoscillatory to oscillatory dynamics at (Fig. 4). We next coupled an array of identical oscillators and observed their dynamics near the transition. In the oscillatory and nonoscillatory regimes the dynamics was synchronized and homogeneous in space with a variation of 15% in period in the oscillating regimes (). However, at the transition, we observed the emergence of spontaneous oscillatory patterns at one edge that slowly penetrated into the array (Fig. 4). The stripe pattern initiated at the edge of the array due to local asymmetry caused by the boundaries of the array and residual flow in the system (). This pattern likely stems from fluctuations near the transition, similarly to the effect of fluctuations near a transition from a monostable to bistable network dynamics in DNA compartments (25).

Summary

Pattern formation and synchrony by coupled biochemical oscillators have been challenging to study in cell-free gene-expression systems, essentially due to the experimental difficulty to assemble spatially distributed reactions that communicate by diffusible signals. The DNA compartment enables steady-state reaction conditions, complex oscillatory dynamics, and diffusion-based communication (15, 25), offering a means to control parameters difficult to access in living systems. The oscillator network used here is based on an activator–repressor − network, but with additional elements that create nonlinear activation threshold by enzymatic repressor degradation and by activator sequestration using the . To this end, the protease complex is encoded in the DNA brush, synthesized, and assembled to target degradation of the repressor.

At the single-compartment level the oscillatory dynamics is controlled by three important parameters: lifetime of the reaction in the compartment, network topology, and gene concentration ratio in the network. We found that nonlinear activation thresholds of activator and repressor are essential to stabilize the oscillatory dynamics. Altering the gene ratio of different regulatory elements in the network revealed that an increase in negative feedback, either by increasing activator or repressor, or by decreasing degradation, slows down the oscillation period. We further assembled multiple identical single-oscillator replicas, and observed variability in oscillation width and period. We find that width variability is higher for activated than repressed genes, whereas the variability in period is of the same order for both genes. We therefore deduce that period variability primarily stems from the regulation of the activated gene. These properties can serve as design principle for implementation of synthetic networks using different transcriptional regulatory elements.

Control over spatial distribution and frequency of the oscillators allowed us to study frequency selection in pairs of coupled oscillators. We found that oscillators coupled by a diffusion of all network elements are entrained to the frequency of the slow oscillator. This result is captured in numerical solutions of the reaction equations and agrees with the observation that negative feedback in the circuit slows down the period.

An additional advantage of the coupled DNA compartment platform is the geometrical control of interaction strength, and population variability both in amplitude and frequency. It is therefore simple to measure a transition from an uncoupled system with different oscillation dynamics to a coupled state, in which oscillations synchronize to a common frequency. Finally, we sought a mechanism for pattern formation in a system that favors synchronization and homogeneous spatial solutions. We introduced symmetry breaking, either spatially by localized concentration gradient, or by spontaneous symmetry breaking close to a transition. These systems reveal a rich spatiotemporal scenario, suggesting plausible mechanisms for pattern formation in developmental biological processes.

To conclude, the construction of complex biochemical systems in vitro provides unique access to the understanding of molecular interactions involved in gene regulation. Our work exemplifies programmable gene expression outside a living organism, from the gene, to a scale of a compartment and to multicompartment collective behavior.