Novel Tunable Spatio-Temporal Patterns From a Simple Genetic Oscillator Circuit

Multicellularity, the coordinated collective behavior of cell populations, gives rise to the emergence of self-organized phenomena at many different spatio-temporal scales. At the genetic scale, oscillators are ubiquitous in regulation of multicellular systems, including during their development and regeneration. Synthetic biologists have successfully created simple synthetic genetic circuits that produce oscillations in single cells. Studying and engineering synthetic oscillators in a multicellular chassis can therefore give us valuable insights into how simple genetic circuits can encode complex multicellular behaviors at different scales. Here we develop a study of the coupling between the repressilator synthetic genetic ring oscillator and constraints on cell growth in colonies. We show in silico how mechanical constraints generate characteristic patterns of growth rate inhomogeneity in growing cell colonies. Next, we develop a simple one-dimensional model which predicts that coupling the repressilator to this pattern of growth rate via protein dilution generates traveling waves of gene expression. We show that the dynamics of these spatio-temporal patterns are determined by two parameters; the protein degradation and maximum expression rates of the repressors. We derive simple relations between these parameters and the key characteristics of the traveling wave patterns: firstly, wave speed is determined by protein degradation and secondly, wavelength is determined by maximum gene expression rate. Our analytical predictions and numerical results were in close quantitative agreement with detailed individual based simulations of growing cell colonies. Confirming published experimental results we also found that static ring patterns occur when protein stability is high. Our results show that this pattern can be induced simply by growth rate dilution and does not require transition to stationary phase as previously suggested. Our method generalizes easily to other genetic circuit architectures thus providing a framework for multi-scale rational design of spatio-temporal patterns from genetic circuits. We use this method to generate testable predictions for the synthetic biology design-build-test-learn cycle.

Multicellularity, the coordinated collective behavior of cell populations, gives rise to the emergence of self-organized phenomena at many different spatio-temporal scales. At the genetic scale, oscillators are ubiquitous in regulation of multicellular systems, including during their development and regeneration. Synthetic biologists have successfully created simple synthetic genetic circuits that produce oscillations in single cells. Studying and engineering synthetic oscillators in a multicellular chassis can therefore give us valuable insights into how simple genetic circuits can encode complex multicellular behaviors at different scales. Here we develop a study of the coupling between the repressilator synthetic genetic ring oscillator and constraints on cell growth in colonies. We show in silico how mechanical constraints generate characteristic patterns of growth rate inhomogeneity in growing cell colonies. Next, we develop a simple one-dimensional model which predicts that coupling the repressilator to this pattern of growth rate via protein dilution generates traveling waves of gene expression. We show that the dynamics of these spatio-temporal patterns are determined by two parameters; the protein degradation and maximum expression rates of the repressors. We derive simple relations between these parameters and the key characteristics of the traveling wave patterns: firstly, wave speed is determined by protein degradation and secondly, wavelength is determined by maximum gene expression rate. Our analytical predictions and numerical results were in close quantitative agreement with detailed individual based simulations of growing cell colonies. Confirming published experimental results we also found that static ring patterns occur when protein stability is high. Our results show that this pattern can be induced simply by growth rate dilution and does not require transition to stationary phase as previously suggested. Our method generalizes easily to other genetic circuit architectures thus providing a framework for multi-scale rational design of spatio-temporal patterns from genetic circuits. We use this method to generate testable predictions for the synthetic biology design-build-test-learn cycle.

INTRODUCTION
Multicellularity and collective cell behavior exemplify the emergence of complex patterns and structures across scales in living systems. When cells interact they can generate higher order patterns of gene expression (differentiation) as well as patterns of mechanical stresses and strains (Chan et al., 2017;Vining and Mooney, 2017). This process takes place in natural phenomena such as embryonic development, tumor formation, wound healing, among others (Velardo et al., 2004;Aboobaker et al., 2005;Khain and Sander, 2006;Gjorevski and Nelson, 2010;Santos-Moreno and Schaerli, 2019). Understanding how these patterns are generated and maintained will enable applications in tissue engineering and regenerative medicine. However, natural emergent multicellular phenomena present numerous unknown processes that pose difficulties for understanding the fundamental mechanisms underlying pattern formation.
Synthetic biology applies design principles to generate combinations of genetic parts that perform a given function, for example oscillation, which at the same time helps us to understand the complexity inherent to natural systems. The prototypical engineering process is the design-build-test-learn cycle, which is an iterative process relying heavily on models of genetic circuit function. A variety of genetic circuits have been designed, analyzed, simulated, and then implemented in this way. These synthetic circuits simplify biological systems reproducing a specific function (Xie and Fussenegger, 2018) such as toggle switches (Gardner et al., 2000;Yeung et al., 2017), oscillators (Elowitz and Leibler, 2000;Stricker et al., 2008;Danino et al., 2010;Potvin-Trottier et al., 2016), logic gates (Tamsir et al., 2011;Nielsen et al., 2016;Green et al., 2017;Kim et al., 2018), and arithmetic operators (Wong et al., 2015;Ausländer et al., 2018).
While these circuits are often studied as dynamical systems in single cells or well mixed populations, the function of genetic circuits has also been studied in cell colonies (Luo et al., 2019;Santos-Moreno and Schaerli, 2019) through the engineering of patterns of gene expression such as symmetry breaking (Nuñez et al., 2017), Turing patterns (Karig et al., 2018), fractal patterns (Rudge et al., 2013), tissue like structures (Toda et al., 2018;Healy and Deans, 2019), among others. These emergent spatio-temporal patterns depend on mechanical constraints on cells, which are the result of cell-cell and cell-substrate interactions. Thus, synthetic gene circuits can be engineered to generate higher order spatio-temporal patterns when coupled to mechanical constraints.
We focus here on the repressilator (Elowitz and Leibler, 2000), a gene network that encodes a ring oscillator topology consisting of three repressors, where repressor 1 inhibits repressor 2, repressor 2 inhibits repressor 3, and repressor 3 inhibits repressor 1 (Figure 1). In the original realization of this circuit topology (Elowitz and Leibler, 2000) the circuit was subject to significant effects of noise and oscillations quickly became desynchronized. Recently, the circuit was revisited by Potvin-Trottier et al. (2016) with microfluidics systems that allowed them to observe single cells oscillating synchronously in chemostatic conditions for long periods of time. In this work sources of noise were reduced in several ways. Firstly, the fluorescent reporters were integrated into the same low-copy plasmid as the repressilator reducing the standard deviation in amplitude greatly. They also removed the degradation tags and used a protease deletion strain ( clpXP) as the chassis to remove noise from enzymatic queuing (Cookson et al., 2011;Steiner et al., 2016). They also increased the effective repression threshold with a high-copy titration "sponge" plasmid that sequesters a proportion of the TetR repressor, since low repression thresholds imply sensitivity to noisy repressor expression levels. These modifications allowed regular and sustained synchronous oscillations that peaked around every 14 generations. The circuit oscillated for approximately 18 periods before accumulating half a period of drift, demonstrating that cells remained synchronized for more than 250 generations. Strikingly, Potvin-Trottier et al. (2016) were able to observe whole flasks of liquid bacterial culture oscillate synchronously, and bacterial colonies form coherent ring patterns at macroscopic scale. These findings show that the repressilator can be effectively isolated from noise, function in a robust and synchronous fashion, and is capable of forming spatial patterns.
Models are essential in the design process, they allow engineers to screen the parameter space looking for possible functional constructions (Endy and Brent, 2001;De Jong, 2002). Synthetic biology has gone from intracellular dynamic models using ODEs (Elowitz and Leibler, 2000;Gardner et al., 2000), and SSA (Potvin-Trottier et al., 2016;Karig et al., 2018) to sophisticated collective behavior models based on individual agents (Rudge et al., 2012;Gorochowski, 2016) and integrated circuit-host models (Sickle et al., 2020). Using cellular scale individual-based models (IBMs) gives rich information about the emergent collective properties of cell populations due to the interactions between themselves and their environment. These models track cell growth and gene expression in ways analogous to experiments performed in controlled environmental conditions with specified properties such as viscosity, chemical concentrations, etc. This makes them an essential tool in the analysis and design of emergent properties of genetic circuits operating in multicellular chassis. However these models are complex and require significant computation time, highlighting the need for simple tractable mathematical and computational methods.
In this study we uncover novel spatio-temporal patterns of gene expression generated by the repressilator in growing cell colonies, and establish a simple method for their design. Since it is generalizable, this work provides a quantitative framework for multi-scale rational design of spatio-temporal patterns from genetic circuits. We provide testable predictions for the synthetic biology design-build-test-learn cycle for engineering repressilator spatio-temporal pattern.

Growth Rate Variation in Growing Microcolonies
We consider the case of Escherichia coli, the cellular chassis for which the repressilator (Figure 1) was designed, growing on a viscous substrate such as a hydrogel or PDMS (polydimethylsiloxane) and supplied with fresh nutrients via microfluidic channels. Growth is constrained by forces between cells and between cells and the substrate due to viscous drag (Rudge et al., 2012). The cells in such a system are constrained to a monolayer and form a quasi-two-dimensional array of extending rod shapes (Farrell et al., 2013;Grant et al., 2014). We used an individual-based model (Rudge et al., 2012) to characterize the distribution of cell growth rates across a two dimensional cell monolayer growing in such conditions over time. We simulated the growth of microcolonies from single cells to populations of approximately 60,000 cells with a radius of approximately 200 cell diameters. Figures 2A-C show the development of a colony from approximately 5,000 to 50,000 cells. The distribution of growth rates across the colony has a clear radially symmetric pattern with a maximum at the edge of the colony ( Figure S1). This is as expected (Vicsek et al., 1990;Smith et al., 2017) since the cells at the edge of the colony are relatively unconstrained. Thus individual bacteria inside microcolonies perceive a different biophysical environment depending on their spatial position. Taking radial averages of growth rate on growing colonies over a range of time points we see the same exponential decay relative to the colony edge ( Figure 2E), leading to a simple model for the cell growth rate as a function of the radial position of the cell with respect to the colony edge r(t),μ where r 0 (8.23 ± 1.69 cell diameters) is the characteristic length scale of the radial variation in growth rate, and we have normalized by µ 0 -the maximal unconstrained growth rate of the cells. These results suggest that growth rate time dynamics are determined by radial distance from the colony edge. After a short transient, the colony edge moves with constant velocity v front = 5.00 so that R max increases linearly ( Figure 2F). Assuming a two-dimensional densely packed monolayer with random cell orientations, growth is isotropic and expansion is equal in all directions. The area expansion rate approximates the growth rate and is given by the divergence of the velocity field, where A is the cell area and v is the velocity. Since growth is isotropic we may decompose the expansion rate equally into its radial and perpendicular components. Considering the velocity v in the radial direction r, and velocity w in the perpendicular direction s, and expanding the divergence term, Equation (2) gives, Hence, with v = dr/dt, and considering only the radial direction, we can rescale time as t → tµ 0 and radial distance as r → r/r 0 to obtain, d dr Integrating by r and t results in, where τ = −2log(exp(r(0)) − 1), and r(0) is the initial radial position of the cell. Equations (5)- (7) give us valuable insights into the system behavior ( Figure 3). The velocity in the radial direction (away from the colony edge) is a sigmoidal logistic function and saturates to a velocity of v = 1/2 as r increases ( Figure 3A). This gives the front velocity as v front = 1/2 in rescaled units. Correspondingly the radial position relative to the colony edge r tends toward linear increase at velocity v = 1/2 as the cell becomes effectively stationary relative to the colony center ( Figure 3B). In real units this means that the front velocity is v front = r 0 /2, where r 0 is the length scale of growth rate variation (Equation 1). This is consistent with our individual based simulations (Figure 2F), in which v front = 5.00 and r 0 = 8.23 ± 1.69. The growth rate is also sigmoidal and tends from maximum at the colony edge to zero as the cell moves away from the growing front of the colony ( Figure 3C).
The critical time τ , depending on the initial cell position, is the time at which the growth rate and velocity are at their half maximum values, and the cell is at radial position r = log(2) (Figure 3, dashed lines). At this time the cell switches from a high FIGURE 2 | Average growth rate decays exponentially with distance relative to the edge. (A-C) Are growing colonies at 5,000, 25,000 and 50,000 cells, respectively. (D) Is a zoom from the white square region in (C). Single cells are colored according to their growth rate, which ranges between 0 and 1. (E) Average growth rate µ at different distances from the colony edge, blue line is an exponential function (e −r/r0 ) with r0 = 8.23 fitted to the binned data (blue dots). Inset shows the log scale of the average growth rate µ at different distances from the edge. (F) Colony radius (R max ) over time (dashed red line), linear fit to data at t>10 giving v front = 5.00 the velocity of the colony edge (blue line).
growth, low velocity regime (remaining close to the colony edge), to a high velocity, low growth regime (remaining stationary while the colony edge propagates). The time τ for this switch to occur is greater for cells closer to the colony edge, that is they remain in the fast growing regime for longer. Thus, cells effectively experience a switch in growth rate and velocity at their critical time τ , which depends on their initial radial position in the colony.

Dynamic Growth Rate Dependent Mathematical Model of the Repressilator
Here we develop a simple mathematical model coupling the repressilator genetic circuit to growth rate variation via simple dilution of proteins. The repressilator can be considered as an abstract genetic circuit topology. We consider an implementation of this topology following the design modifications made by Potvin-Trottier et al. (2016), which essentially isolated the circuit from noise and allowed sustained and synchronous oscillations over time scales up to 250 generations. Stochastic simulations performed with relevant parameters reproduced this behavior, showing essentially deterministic behavior ( Figure S2), therefore we may use a simpler differential equation model to track the repressor concentration of each cell over time. A simple two-step model of the balanced repressilator genetic circuit (Figure 1), a type of genetic ring oscillator, can be formulated as follows, where i is one of the three genes, j is its corresponding repressor gene, m i , p i are the mRNA and protein concentrations respectively, a is the constitutive transcription rate, b is the leaky or repressed transcription rate, µ is the instantaneous growth rate of the cell or population of cells, γ is the protein degradation rate, and δ is the mRNA degradation rate.
Order of magnitude estimates of these parameters are given in Supplementary Material. Since mRNAs are typically short lived (see Supplementary Material), we may assume quasi-steady state concentrations and the system is, Rescaling protein concentration as p j → p j /K j and time by t → tµ 0 with µ 0 the maximal growth rate, assuming that basal Frontiers in Bioengineering and Biotechnology | www.frontiersin.org FIGURE 3 | Analytical solutions for cell velocity (A) and position (B) relative to the colony edge, and growth rate (C), based on exponential decay of growth rate with distance from colony edge. At critical time τ the system switches between two regimes: high growth/low velocity, and low growth/high velocity.
expression is negligible, and combining with (Equation 7), where α = ca/δµ 0 K (order of magnitude 10 4 , see Supplementary Material) is the steady state maximal gene expression rate and the constantγ = γ /µ 0 is the protein degradation rate as a fraction of the maximal growth rate (order of magnitude 1). This model depends on three dimensionless parameters α,γ , and n.
In this model we assume that the dominant effect of growth rate variation is by dilution of proteins. While there is evidence for growth rate dependencies of transcription and translation rates and plasmid copy number (Neubauer et al., 2003;Klumpp et al., 2009;Klumpp, 2011), all of which affect the parameters of the model, these effects have only been observed due to different biochemical environments. In spatially constrained cell populations the shape of the growth profileμ (t) depends both upon the biochemical environment and mechanical constraints (Andersen and von Meyenburg, 1980;Matsushita and Fujikawa, 1990;Tuson et al., 2012;Farrell et al., 2013;Rudge et al., 2013;Smith et al., 2017;Winkle et al., 2018). At the typical bacterial microcolony scale the biochemical environment is essentially uniform in space due to the fast diffusion of small molecules like sugars and aminoacids (Matson and Characklis, 1976;Fraleigh and Bungay, 1986;Guélon et al., 2012). Using microfluidic devices cells can be maintained in constant fresh media allowing observation of the long term dynamics of genetic circuits (Danino et al., 2010;Long et al., 2013;Potvin-Trottier et al., 2016). Under these conditions then the predominant factors determining growth rate are physical forces and constraints.

Protein Dilution Enables the Repressilator to Produce Traveling Waves in Growing Microcolonies
The model presented above (Equations 5-7 and 11) describes the trajectories of cells as they move in the radial direction and change their protein concentrations over time. Assuming that the motion and growth of cells is not affected by the expression of repressor genes, this model describes the mean behavior of cells starting from some initial radial position. It is obvious from these equations that in the absence of growth the system can only produce plane wave, homogeneously synchronized oscillations. However, in the presence of growth we have an explicit relation between cell position and protein expression rate, enabling spatio-temporal pattern formation.
Since growth dilutes proteins the effective degradation rate of the repressors isγ +μ(t). The effect of the sigmoidal growth rate switch on the repressilator is therefore to decrease the effective degradation rate of the repressors fromγ + 1 toγ as cells move out of the growing regime (Equation 7). Potvin-Trottier et al. (2016) showed that increasing the degradation rate of repressors by appending a degradation tag reduced the period of oscillations T. This decrease was by approximately a factor of two at 37 • C, with less effect at lower temperatures likely due to decrease in protease activity (Purcell et al., 2012). We confirmed this result using our model by numerically integrating Equation (11) at fixed effective degradation ratesγ +μ (Figure 4). The frequency of oscillations f = 1/T was proportional to the degradation rate, with a slope depending on α, the maximum gene expression rate. Hence in colonies, as cells move away from the edge due to mechanical constraints the effective repressor degradation rate decreases and the period of their oscillations increases.
After the critical time τ the period of oscillations increases as the cell switches from high growth rate and low velocity to low growth rate and high velocity (Equations 5-7). This means that there is effectively an interior region oscillating with long period T int and an exterior region oscillating with short period T ext . The phase offset between peaks of the two signals after some time t is, FIGURE 4 | Oscillation frequency f (or period T) depends on the effective degradation rate of repressorsγ +μ. For a given maximum gene expression rate α the frequency is proportional to the effective degradation rate, and increasing α decreases the frequency.
When the phase offset T = T ext we have in phase oscillations. The time required for a cell to achieve this phase offset is the time spent in the high growth regime τ , plus the time spent in the low growth regime, hence, is the time at which the cell is in phase with the edge of the colony, the wave source. At time t * the distance from the edge r(t * ) of this cell can be obtained from Equation (6), and since this is the peak-to-peak distance it gives the wavelength λ. Assuming that exp(t * ) ≫ 1, The wave propagation velocity is v p = λ/T ext −v front = λ/T ext − 1 2 , wherev front is the velocity of the colony edge in normalized distance units, hence, Equations (14) and (15) show that when T ext < T int the system generates traveling waves with finite wavelength and wave speed. When T int = T ext , there is no effect of mechanical constraint on oscillation period and we find that v p = ∞ and λ = ∞, meaning that the system forms homogeneous plane waves with the whole colony oscillating synchronously. As T int → ∞ the interior does not oscillate and we find that v p → 0 and λ → T ext /2 and we thus have static rings of gene expression with spatial wavelength T ext /2 =v front T ext . This is the case when protein degradation is negligible (γ = 0), a condition under which the repressilator has been shown to form static rings in growing colonies (Potvin-Trottier et al., 2016). Thus we have shown analytically that growth rate heterogeneity induces the repressilator to form either static rings or traveling waves in growing cell colonies, depending on the degradation rate of the repressors.

Novel Spatio-Temporal Patterns Emerging From Repressilator Dynamics
To test the predicted spatio-temporal patterns we integrated Equation (11) from a range of initial cell radial positions to construct the kymograph p i (R, t), where R = t/2 − r is the rescaled distance from the center of the colony. A kymograph (Figures 5A-C) represents the spatial dynamics of a one-dimensional system, such as ours, evolving over time.
Each point in the kymograph represents the state of the system at a given time (x-axis) and position (y-axis) with a color. By taking radial averages the kymograph fully characterizes radially symmetric spatio-temporal patterns with the vertical axis representing distance from the center of the colony. Here we reflect the kymograph to represent the symmetric structure of the pattern (Figures 5A-C), whereby the growth of the colony can be seen as two linearly expanding borders forming a triangular shape. The slope of this border is the front speed v front = 1/2. The color represents the radially averaged repressor protein concentrations (red, green, and blue) at each position at each time point, normalized to their maximum values. The corresponding predicted colony pattern is shown inset in Figure 5C. Stripes in the kymograph represent rings of gene expression. Horizontal stripes show static rings since their radial position does not change ( Figure 5A). Vertical stripes would represent in-phase homogeneous oscillations since they do not vary in space ( Figure 5B). Diagonal stripes represent traveling waves, moving rings of repressor expression, since they vary in both space and time ( Figure 5C). Hence we confirm our theoretical prediction of traveling wave patterns. The spatio-temporal dynamics of the system are described by two parameters that can be extracted from the kymographs. The slope of each stripe gives the wave speed v p , and the vertical peakto-peak distance gives the spatial wavelength λ = (v p +v front )T = (v p + 1 2 )T, where T is the period of oscillations at the colony edge (the wave front) andv front is the front velocity in normalized distance units (Figure 5C).

Tuning the Repressilator to Control Spatio-Temporal Pattern Formation
In order to quantitatively test our theoretical predictions and to characterize the design space of these traveling wave patterns we scanned the parameter space within physiologically relevant ranges. We measured the wave speed v p and wavelength λ of the system for 625 combinations of α andγ spanning four orders of magnitude to construct the phase space (Figures 5D,E). We see that wavelength was predominantly determined by α, while wave speed depended onγ . Traveling waves were observed for all values of α but clearly require non-zero protein degradation rateγ ( Figure 5F).
We observed that wave speed was proportional toγ , with v p ≈γ /2 (linear fit v p = 0.535γ + 0.0214) . Static rings (v p = 0) form whenγ = 0. Asγ → ∞, as is the case at growth arrest (µ 0 = 0), we saw earlier that the system can only form plane waves with all parts of the colony oscillating in phase, which corresponds to v p = ∞. Wavelength was predominantly but weakly affected by maximum gene expression rate α (Figure 5G) following approximately α ≈ 10 λ−1 [from the linear fit λ = 1.01 log 10 (α) + 0.998 ]. These results are consistent with our theoretical predictions based on the oscillation period from Equations (14) to (15) (Figures 5F,G triangles).
We demonstrate the tuning parameters α andγ above using kymographs in Figure 6. With no protein degradation (γ = 0) the system produces static rings of gene expression following the phase of each repressor (Figures 6A,B). In the kymograph this spatio-temporal pattern is observed as horizontal stripes of consecutive red, green, and blue, representing the three repressors. This confirms the observation of fixed ring patterns in colonies hosting a repressilator with stable repressor proteins (Potvin-Trottier et al., 2016). The static ring patterns observed are therefore a special case of the more general traveling wave solution with velocity v p = 0. These traveling waves are induced and modulated by protein degradation. In the intermediate case when protein degradation and growth are similar (Figures 6C,D,F) we see the clear emergence of a traveling wave solution. This spatio-temporal pattern is characterized by diagonal stripes in the kymograph, with steeper sloped lines indicating higher velocity of the waves (Figure 5A). At lower protein degradation rate (γ = 0.3) we see traveling waves with lower velocity (Figures 6C,D). Hence protein degradation rate tunes the speed of the traveling waves. Changing α however does not affect the speed of the waves (Figures 6A-F) but does change the wavelength of the traveling waves resulting in more spatial rings at lower α values. We note that increasing α also stabilizes the oscillations as found by Osella and Lagomarsino (2013) and Potvin-Trottier et al. (2016) (Figure 6 panels below each kymograph).

Growing Cell Colonies Generate Traveling Waves in Quantitative Agreement With Predictions
To test if these predictions hold in constrained growing microcolonies of cells we used our individual based biophysical model of bacterial cell growth and division. We grew colonies from a single cell up to 60,000 cells, tracking each cell's motion and protein expression levels according to Equation (11) without the growth rate term (dilution was computed by the biophysical model). The results show, as predicted, the formation of symmetrical rings relative to the center of the colony ( Figure 7A, Supplementary Material, Video 1 and Video 2).
In order to test the dependency of wavelength and wave speed on protein degradation and maximal expression rate, we simulated a range ofγ and α (kymographs in Figures 7B-E). Our findings matched with the prediction of the 1D model; no waves were formed forγ = 0, wave speed increased whenγ was increased, and wavelength increased when α was increased. The spatio-temporal dynamics of each repressor is regulated by protein dilution (γ ), which moves the system from fixed rings (Figure 7B) to an oscillatory behavior which gives rise to traveling waves with different wavelengths controlled by maximum expression levels and speeds controlled by protein degradation (Figures 7C-E).
Since we tracked every cell as they grow, replicate, and express proteins (Figures 7B-E right column) we could follow the dynamics of individual cells as they move through the colony, changing their growth rate depending on their mechanical environment (Figures 7B-E middle column). We selected central and peripheral cells for each of the colonies to study the most restricted and the most unconstrained cells. For colonies with traveling waves the constrained non-growing central cells oscillated with constant frequency and phase. Cells starting at the periphery of the colony however experience changes in growth rate as they move away from the colony edge that cause a sharp decrease in frequency and a resulting phase offset with respect to the central cell. We found that cells in the periphery exhibit higher frequency oscillations compared to central cells, and that difference is increased when increasingγ (Figures 7B-E central  column). This is consistent with our theoretical prediction that growth rate reduction increases the period of oscillations, causing a phase offset between the interior and peripheral regions of the colony.
The wavelength and wave speed obtained from growing microcolonies was closely correlated to the predictions of our simple model (Pearson's correlation coefficient 0.983 for wavelength and 0.999 for wave speed, Figure 8). The length scale of wave speed and wavelength is set by r 0 , hence fitting a linear model between the predicted and simulated speed and wavelengths gives an estimate of r 0 for the growing colony. From the wave speed values we obtained r 0 = 11.6 and from wavelength r 0 = 9.01, which is in close agreement with that estimated from the growth rate distribution of colonies (Figure 2). Our one dimensional model predicts that the front velocity of the colony should be v front = r 0 /2. From wave speed we obtain a value of v front = 5.80 and from wavelength v front = 4.51, which is extremely close to the estimated value of v front = 5.00 for our individual based model (Figure 2).

Manipulating Mechanical Growth Constraints to Control Pattern Formation
Microfabricated cell culture devices and microfluidics provide fine control over the mechanical as well as biochemical conditions in which cells grow. As well as providing fresh nutrients via flow, maintaining cells in steady state, these devices provide techniques to physically constrain cell growth and therefore another mode of design of spatio-temporal patterns induced by growth rate heterogeneity as we have demonstrated. Commonly microfluidic devices are designed to constrain cells to a monolayer, while allowing loading of seed cells into a chamber or channel (Figures 9A,B). We imposed two such constraints on our colonies. In a long thin channel (400 × 20 × 1 cell diameters) cells form one dimensional traveling waves directed along the channel axis ( Figure 9A, See Supplementary Material, Video 4, Video 5, Video 6). When constrained to a chamber (80 × 80 × 1 cell diameters) we observed the emergence of traveling waves during unconstrained growth (Figure 9B, t+4, See Supplementary Material, Video 3). These waves were sustained over long time periods after the cells became constrained and stopped growing altogether (Figure 9B, t + 8, t + 12). Growth is necessary to form traveling waves but the established phase offsets between different radial positions are maintained after growth arrest, continuing to produce traveling waves. The history of the shape of wavefront is therefore retained in the pattern. Finally, we demonstrate that growth rate heterogeneity in three dimensional colonies also generates traveling waves as layers (Figure 9C), showing that the spatio-temporal pattern is not specific to monolayers (see Supplementary Material, Video 7 and Video 8).  . The length scale of wave speed and wavelength is set by r 0 and the slope of these plots gives estimates for r 0 as 11.6 (wave speed) and 9.01 (wavelength) cell diameters. Speed and wavelength of traveling waves were closely correlated between the models, with Pearson's correlation coefficient 0.999 for wave speed, and 0.983 for wavelength.

DISCUSSION
Here we demonstrated how biophysical constraints on growth can induce spatio-temporal pattern formation from a simple genetic circuit. By coupling the repressilator (Potvin-Trottier et al., 2016) to a heterogeneous growth rate pattern via protein dilution we generated emergent traveling waves of gene expression. These traveling waves can be characterized by two properties; the wavelength and the wave speed. These properties are determined by two simple parameters that are feasible to control experimentally; the protein degradation rate, which controls the wave speed, and the maximal protein expression rate, which controls the wavelength. Our results make quantitative and qualitative predictions about the spatio-temporal patterns produced by the repressilator in growing cell colonies.
Our analysis predicts that traveling waves will be observed if the ratio of protein degradation to growth rateγ = γ /µ 0 , is sufficiently high for the waves to form in the time of the experiment. Forγ = 0 we predict the formation of static rings of gene expression as observed in experiments (Potvin-Trottier et al., 2016), however we show here that this pattern could be generated purely by protein dilution and does not require cells to transition into stationary phase. We show that increasingγ , which means increasing protein degradation rate or decreasing growth rate, will increase the speed of the waves. This could be achieved by choosing one of several protein degradation tag sequences that target the proteins for proteolysis (Purcell et al., 2012). Further, our model suggests that increasing maximum protein expression rate α, for example by choosing a more efficient ribosome binding site (RBS) (Salis et al., 2009) will increase the wavelength of the pattern. We parameterize simple empirical models for the effects of each of these genetic design modifications; log 10 (α) = λ − 1 and v p = γ /2 . Thus, we have effectively generated a quantitative datasheet for the repressilator gene circuit topology operating in a simple microcolony chassis.
We derive a simple model of coupling genetic circuits to growth rate via protein dilution, and show that it accurately predicts traveling wave patterns in growing cell colonies, their speed, and wavelength. The model also accurately predicts the colony front velocity. The mathematical and computational approach outlined here is not specific to the repressilator nor to bacterial colonies and could make predictions about spatial patterns produced by other circuit topologies and chassis. Here we did not consider gene circuits that affect growth rate, for example by regulation of metabolism, which may produce more complex spatial patterns (Nuñez et al., 2017), however it could be included in our framework leading to a more complex set of coupled differential equations. Thus, this analysis approach implements the rational design of spatio-temporal patterns of gene expression, enabling the design stage of the design-build-test-learn cycle.
Oscillators are important in regulation of multicellular systems and many studies have reproduced oscillations in synthetic genetic circuits by assembling different devices combining modular parts (Liu et al., 2015;Niederholtmeyer et al., 2015;Perez-Carrasco et al., 2018;Riglar et al., 2019). Studying and engineering synthetic oscillators can direct us to understand complex multicellular behaviors at multiple scales, in particular here the emergence of traveling waves of gene expression in populations of cells. There are a wide range of phenomena in which a key element to a developmental process is the appearance of a traveling wave of chemical concentration, mechanical deformation (Espeso et al., 2016), electrical or other type of signal (Murray, 2002). Two examples are the chemical and mechanical waves which propagate on the surface of many vertebrate eggs (Deneke and Di Talia, 2018). A developing embryo presents a large number of wave like events that appear after fertilization (Kimmel et al., 1995). Thus, one importance of this work is that we were able to rationally design and manipulate in silico genetic circuits to recapitulate such patterns with tunable wavelength and wave speed.
Noise is known to affect oscillators in various ways including stochastic coherence which makes the oscillations more consistent (Hilborn and Erwin, 2008), and may therefore Traveling waves start at the sides and merge in the center (white arrows). (B) Growing colony in a square simulated microfluidic chamber. Before the colony reaches the constraints it grows with radial symmetry and initiates traveling waves (t, t+4). At t+12 the colony has used all the space and is constrained by chamber, growth arrest occurs, but the traveling waves continue (white arrows). (C) Spherical colonies grow when not constrained to a plane, forming traveling layers of gene expression. Image shows a cutaway of half the colony. (D) Transversal section of (C).
stabilize spatio-temporal patterns to stochastic fluctuations in gene expression. We do not consider the role of noise in this study because at the parameter values we explore, the stochastic behavior approximates the continuous model, with regular and sustained synchronous oscillations (Woods et al., 2016) (Figure S2). We also note that since (Potvin-Trottier et al., 2016) observed synchronous long-term oscillations that form ring patterns in growing colonies, cells must be synchronized on average over long length and time scales, and so noise is not likely to be important in the pattern formation process described here. However, it would be interesting to consider the role of noise in the generation of spatio-temporal patterns (Sagués et al., 2007;Zhou et al., 2008) due to lower gene copy or other circuit properties (Vilar et al., 2002;Lestas et al., 2010).
We reason that the traveling waves described here are generated by phase and frequency changes induced by reduction in growth rate as cells become more distant from the edge of the colony, but maintained by protein degradation. In the absence of protein degradation no traveling waves but simple static rings form (Potvin-Trottier et al., 2016). The phase differences are locked in as growth rate decays to zero, such that even after total growth arrest the traveling waves continue (Figure 9B). The scale of the waves, their speed and their wavelength are determined by r 0 , the characteristic length scale of decay in growth rate. However, the radius of the colony also scales with r 0 and so the overall pattern is in a sense scale invariant, showing the same relative wavelength and speed for any exponential growth rate profile.
Our results show that the speed of traveling waves in growing bacterial colonies is approximately 10 cell diameters per doubling time (approximately 10µm per hour for E. coli) toward the colony center, but the colony border grows at only around 4 cell diameters per doubling. Hence gene expression information can be transferred faster via a traveling wave than by the physical transmission of cells. The ability to tune the wavelength λ and wave speed v p of these patterns could enable design of novel cell-cell communication systems based on oscillatory signals. Further, coupling the oscillator to production of pulses of diffusing chemicals such as acyl-homoserine lactones (AHLs) could be used to enhance information transmission (Hopfield, 1974;Mangan et al., 2003). We note also that in a sense the traveling wave pattern, its speed and wavelength encode the history of the shape of the wavefront as the colony expands, which may be useful for example for information storage.
A fundamental result of this work is to demonstrate that mechanical constraint gives rise to higher order gene expression patterns in cell colonies, and provide such a system for analysis. There are a vast number of experimental conditions which could be created to induce different spatio-temporal patterns in such microcolonies. Microfluidics has shown to be of particular help to control mechanical constraints (Ruprecht et al., 2017), substrate stiffness (Wang et al., 2018), nutrients (Alnahhas et al., 2019), chemical inducers (Danino et al., 2010), cell-cell signaling (Alnahhas et al., 2019), and pattern formation (Kantsler et al., 2020). As we showed in Figure 9, controlling biophysical constraints using different channel layouts and mechanical properties of the substrate could produce different patterns of growth rate that give rise to structures that mimic different stages of the development of organisms (Johnson et al., 2017;Toda et al., 2018). A simple example is the one dimensional channel (Figure 9A) which mimics in a simple way an embryo growing along its axis and sending back waves of gene expression from the front of the cell population.
In summary we report here novel traveling wave spatiotemporal patterns resulting from the growth rate dependent dynamics of a repressilator genetic oscillator circuit. We developed an analytical framework to predict the spatio-temporal behavior of such genetic circuits in growing colonies. This framework allows multi-scale rational design of spatio-temporal patterns from genetic circuits and makes testable predictions for the synthetic biology design-build-test-learn cycle.

Individual Based Model
We grew colonies from 1 up to 60,000 cells, simulated using CellModeller (Rudge et al., 2012) with parameters Ŵ = 10 and t = 0.05. Ŵ = γ cell /γ s is the ratio of cell stiffness to substrate stiffness, which we estimate order of magnitude 10 (see Supplementary Material). Briefly, CellModeller simulates cells as rods extending along their axis but otherwise rigid. As cells expand the resulting constraint energy is minimized to find the new arrangement of cells. Cells experience viscous drag with the substrate (γ s ) and along their length axis (γ cell ), and divide when they reach a target length set to l 0 = 3.5 cell diameters. At division the cell is divided into two equal sized rods, which are randomly perturbed slightly in their axis orientation. Cells were constrained to lie in a plane, except in Figure 9C in which cells grew in three dimensions.
Colonies were grown for approximately 48 doubling times and the final radius of the colonies was approximately ∼230 cell diameters. Since these simulations correspond to E. Coli cells, these units represents approximately ∼230µm. Simulations were stored in a file for each time step. This file contains information about the state of each cell present in the colony, including the position, protein concentrations, growth rate, volume, length, among other variables.

Colony Growth Analysis
Growth rate µ and radial position R were obtained for each cell from 3 growing colonies from 5,000 up to 60,000 cells. At each time point the colony radius R max (t) was calculated and divided into n bins of size r = 5 cell diameters from the edge b 0 to the center b n . The growth rates µ of all cells with r ∈ [b n , b n+1 ) were averaged to get µ . An exponential of the form e −r(t)k was fitted to µ at each time, obtaining k at colonies with different R max . A linear model was fitted to the colony radius R max (t) when t > 10 to compute the front velocity v front using numpy.polyfit.

Kymograph Construction for 1D Model
To obtain values of p i (r, t) we integrate Equation (11). Starting from some initial colony radius R 0 we initialize p i (r, 0) for r a regularly spaced lattice on (0, R 0 ). We use p 2 (r, 0) = 5 for all r, that is homogeneous expression of only p 2 . At each step of an explicit Euler integration scheme we find the new cell positions and construct a new regularly spaced lattice in (0, R(t)) = (0, R 0 + t/2) and interpolate p i (r, t + t) onto this grid before taking the next integration step. The algorithm is as follows: 1. Initialize r(0) as a regular lattice on (0, R 0 ) and p * i (r, 0) to some values. 2. Compute p i (r, t + t) by an explicit Euler step of Equation (11) , and r(t + t) using Equation (6). 3. Compute a new regular mesh r ′ (t + t) on (0, R(t + t)). 4. Interpolate the protein concentrations to get p i (r ′ , t + t). 5. Set t → t + t, and r(t + t) → r ′ (t + t), and repeat from step 2.
At the end of this procedure we have constructed a set of samples p i (r, t) which we then interpolate to form the kymograph.

Dynamical Simulations of Gene Expression
Using the file stored for each simulation in the IBM, we have a representation of the biophysical model decoupled from the genetic circuit. Using these data we performed simulations of the gene expression model derived in Equation (11). In order to keep updating the state of the cells, which is affected by cell division, we constructed a graph of parent-child relations. Thus, we integrate Equation (11) forward using explicit Euler method between each state of the biophysical model. One assumption made is that when a cell divides the children inherit the value of the protein concentration his parent. We assume the number of proteins divides equally between the two cells, as does the volume of the cell, keeping protein concentration constant. Resultant simulations then serialized to a JSON file. These files were later used to perform analysis and create visual representations.

Kymograph Construction for Individual Based Simulations
Using the JSON file obtained in the temporal simulation of gene expression with the biophysical model, we generated positions and growth rates of cells. Then we binned cells according to their radial position using bin size R = 5. Finally we take the average protein concentration in each bin and repeat for all time steps to get p i (R, t).

Wave Speed Estimation
First we take each row of the kymograph and identify the radial peaks (scipy.signal.find_peaks) in each protein concentration for each time step. Next the peaks are paired with the nearest peak in the previous time step, and the average distance between them used to calculate the wave speed as v = x /δt, where x is the mean peak to peak distance and δt is the simulation time step.

Wavelength Estimation
In order to estimate the wavelength λ of the traveling waves we note that λ = (v p + v front )T, where v p is the wave speed, T is the period of oscillations, and v front = 1 2 is the velocity of the colony edge or wavefront. To estimate the oscillation period we take the leading edge of the colony and compute the peaks in its time varying protein concentration p i (r = 0, t). Then as above we estimate the period as the mean of the peak to peak times so that T = t . The wave speed is taken from the calculations described above, and the resulting estimate for wavelength is λ = (v p + 1 2 )T = (v p + 1 2 ) t .

Analytical Estimates of Wavelength and Wave Speed
We used Equations (14)-(15) to estimate the wavelength and wave speed of traveling waves that the repressilator would produce with an exponential growth rate profile (Equation 1). First we numerically integrated Equation (11) with fixed growth rateμ. For each combination of parameters we simulated oscillations at the colony edge (μ = 1) and the colony interior (μ = 0), and measured the periods T ext and T int as described above. These values were then substituted into Equations (14)-(15) to compute the estimated wavelength and wave speed.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.github. com/synbiouc/spatialoscillator.