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

1. 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.

FIGURE 1

Figure 1. Repressilator genetic oscillator circuit. (A) Network representation of the repressilator, in the schematic 1 represses 2, 2 represses 3, and 3 represses 1. (B) Genetic circuit diagram of a plasmid encoding the repressilator. (C) Oscillating repressor concentrations over time computed by solving Equation (11) with fixed growth rate. CDS 1, 2, and 3 are the complementary DNA sequences coding for each repressor, and Ori is the origin of replication.

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.

2. Results

2.1. 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),

$\begin{array}{ll}\bar{\mu }(t)=e^{-r(t)/r_0},& (1)\end{array}$

where r₀ (8.23 ± 1.69 cell diameters) is the characteristic length scale of the radial variation in growth rate, and we have normalized by μ₀—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).

FIGURE 2

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).

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,

$\begin{array}{ll}\nabla ·\text{v}=\frac{1}{A}\frac{dA}{dt},& (2)\end{array}$

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,

$\begin{array}{ll}\frac{dv}{dr}+\frac{dw}{ds}=2\frac{dv}{dr}=\mu (r).& (3)\end{array}$

Hence, with v = dr/dt, and considering only the radial direction, we can rescale time as t → tμ₀ and radial distance as r → r/r₀ to obtain,

$\begin{array}{ll}\frac{d}{dr}\left(\frac{dr}{dt}\right)=\frac{1}{2}{e}^{-r}.& (4)\end{array}$

Integrating by r and t results in,

$\begin{array}{ll}v\left(t\right)=\frac{1}{2}{\left(1+\text{exp}\left(\frac{\tau -t}{2}\right)\right)}^{-1},& (5)\end{array}$

$\begin{array}{ll}r\left(t\right)=\text{log}\left(1+\text{exp}\left(\frac{t-\tau }{2}\right)\right),& (6)\end{array}$

$\begin{array}{ll}\bar{\mu }\left(t\right)={\left(1+\text{exp}\left(\frac{t-\tau }{2}\right)\right)}^{-1},& (7)\end{array}$

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₀/2, where r₀ 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₀ = 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).

FIGURE 3

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.

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 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.

2.2. 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,

$\begin{array}{ll}\frac{d{m}_i}{dt}=\frac{a+b{(p_j/K_j)}^n}{1+{(p_j/K_j)}^n}-\delta m_i,& (8)\end{array}$

$\begin{array}{ll}\frac{d{p}_i}{dt}=c{m}_i-\gamma p_i-\mu (t)p_i,& (9)\end{array}$

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,

$\begin{array}{ll}\frac{d{p}_i}{dt}=\frac{c}{\delta }\frac{a+b{(p_j/K_j)}^n}{1+{(p_j/K_j)}^n}-\gamma p_i-\mu (t)p_i.& (10)\end{array}$

Rescaling protein concentration as p_j → p_j/K_j and time by t → tμ₀ with μ₀ the maximal growth rate, assuming that basal expression is negligible, and combining with (Equation 7),

$\begin{array}{ll}\frac{d{p}_i}{dt}=\frac{\alpha }{1+p_j^n}-\overline{\gamma }{p}_i-p_i{\left(1+\exp \left(\frac{t-\tau }{2}\right)\right)}^{-1},& (11)\end{array}$

where $\alpha =ca/\delta {\mu }_{0}K$ (order of magnitude 10⁴, see Supplementary Material) is the steady state maximal gene expression rate and the constant $\overline{\gamma }=\gamma /{\mu }_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 α, $\bar{\gamma }$, 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 $\bar{\mu }(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.

2.3. 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 $\bar{\gamma }+\bar{\mu }(t)$. The effect of the sigmoidal growth rate switch on the repressilator is therefore to decrease the effective degradation rate of the repressors from $\bar{\gamma }+1$ to $\bar{\gamma }$ 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 $\bar{\gamma }+\bar{\mu }$ (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.

FIGURE 4

Figure 4. Oscillation frequency f (or period T) depends on the effective degradation rate of repressors $\bar{\gamma }+\bar{\mu }$. For a given maximum gene expression rate α the frequency is proportional to the effective degradation rate, and increasing α decreases the frequency.

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,

$\begin{array}{ll}\Delta T=\left(\frac{T_{int}-T_{ext}}{T_{int}}\right)\Delta t.& (12)\end{array}$

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,

$\begin{array}{ll}{t}^{*}=\tau +\Delta t=\tau +\frac{T_{int}{T}_{ext}}{T_{int}-T_{ext}},& (13)\end{array}$

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,

$\begin{array}{ll}\lambda =\frac{T_{int}{T}_{ext}}{2(T_{int}-T_{ext})}.& (14)\end{array}$

The wave propagation velocity is $v_p=\lambda /T_{ext}-{\bar{v}}_{front}=\lambda /T_{ext}-\frac{1}{2}$, where ${\bar{v}}_{front}$ is the velocity of the colony edge in normalized distance units, hence,

$\begin{array}{ll}{v}_p=\frac{T_{ext}}{2(T_{int}-T_{ext})}.& (15)\end{array}$

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={\bar{v}}_{front}{T}_{ext}$. This is the case when protein degradation is negligible ($\bar{\gamma }=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.

2.4. 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 ${\bar{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.

FIGURE 5

Figure 5. Parameter dependence characterization for 1D model. (A–C) Kymographs of growing colonies show protein concentrations as a function of radial position over time, with no protein degradation (A, $\bar{\gamma }=0$, α = 10⁴) we obtain static rings, with no growth rate dilution (B, $\bar{\gamma }=1.5$, α = 10⁴) we see plane waves, and with growth rate dilution and protein degradation (C, $\bar{\gamma }=1.5$, α = 10⁴) we have traveling waves. White dashed lines show wave trajectories. The distance between two trajectories is the wavelength λ, and the slope is the wave speed v_p. Inset in (C) is the whole colony at the end of the experiment. (D,E). Heatmaps of wave speed v_p and wavelength λ, respectively over a range of α and $\bar{\gamma }$. (F). Effect of $\bar{\gamma }$ on wave speed v_p at different α. (G). Effect of α on wavelength λ at different $\bar{\gamma }$. Triangles in (F,G) show analytical estimates for α = 100 and $\bar{\gamma }=1$, respectively.

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 peak-to-peak distance gives the spatial wavelength $\lambda =(v_p+{\bar{v}}_{front})T=(v_p+\frac{1}{2})T$, where T is the period of oscillations at the colony edge (the wave front) and ${\bar{v}}_{front}$ is the front velocity in normalized distance units (Figure 5C).

2.5. 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 $\bar{\gamma }$ 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 $\bar{\gamma }$. Traveling waves were observed for all values of α but clearly require non-zero protein degradation rate $\bar{\gamma }$ (Figure 5F).

We observed that wave speed was proportional to $\bar{\gamma }$, with $v_p\approx \bar{\gamma }/2$ (linear fit $v_p=0.535\bar{\gamma }+0.0214$) . Static rings (v_p = 0) form when $\bar{\gamma }=0$. As $\bar{\gamma }\to \infty$, as is the case at growth arrest (μ₀ = 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₁₀(α)+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 $\bar{\gamma }$ above using kymographs in Figure 6. With no protein degradation ($\bar{\gamma }=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 ($\bar{\gamma }=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).

FIGURE 6

Figure 6. 1D model simulations for selected parameters. Each panel shows: kymograph of the emerged pattern, where stripes show wave trajectories; trace of the three repressor concentrations over time in the center of the colony; the whole colony showing the final ring pattern. (A) $\bar{\gamma }=0.0$, α = 10,000. (B) $\bar{\gamma }=0.0$, α = 100. (C) $\bar{\gamma }=0.3$, α = 10,000. (D) $\bar{\gamma }=0.3$, α = 100. (E) $\bar{\gamma }=1.0$, α = 10,000. (F) $\bar{\gamma }=1.0$, α = 100.

2.6. 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).

FIGURE 7

Figure 7. Simulations of growing colonies. (A) Colonies with 5,000, 20,000, 35,000, and 50,000 cells, equally spaced 9.3 doubling times apart, with $\bar{\gamma }=0.3$, α=10.000. (B–E) Each panel shows: kymograph of repressor concentrations (51 doublings, approximately 60,000 cells); time dynamics for central cell and peripheral cell in colony; close up of edge of colony at end of experiment. Parameters: (B) $\bar{\gamma }=0$, α = 10, 000. (C) $\bar{\gamma }=0.3$, α = 10, 000. (D) $\bar{\gamma }=0.3$, α = 100. (E) $\bar{\gamma }=1.0$, α = 1, 000.

In order to test the dependency of wavelength and wave speed on protein degradation and maximal expression rate, we simulated a range of $\bar{\gamma }$ and α (kymographs in Figures 7B–E). Our findings matched with the prediction of the 1D model; no waves were formed for $\bar{\gamma }=0$, wave speed increased when $\bar{\gamma }$ was increased, and wavelength increased when α was increased. The spatio-temporal dynamics of each repressor is regulated by protein dilution ($\bar{\gamma }$), 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 $\bar{\gamma }$ (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₀, hence fitting a linear model between the predicted and simulated speed and wavelengths gives an estimate of r₀ for the growing colony. From the wave speed values we obtained r₀ = 11.6 and from wavelength r₀ = 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₀/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).

FIGURE 8

Figure 8. Comparison of wave speed (A) and wavelength (B) between 1D model and growing cell colonies (IBM). The length scale of wave speed and wavelength is set by r₀ and the slope of these plots gives estimates for r₀ 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.

2.7. 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).

FIGURE 9

Figure 9. Traveling waves in microfluidic devices. (A) Simulations of a monolayer of cells in a narrow infinite channel 20 cell diameters wide and one cell diameter tall. 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).

3. 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 $\bar{\gamma }=\gamma /{\mu }_0$, is sufficiently high for the waves to form in the time of the experiment. For $\bar{\gamma }=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 $\bar{\gamma }$, 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₁₀(α) = λ−1 and $v_p=\bar{\gamma }/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 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₀, the characteristic length scale of decay in growth rate. However, the radius of the colony also scales with r₀ 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 spatio-temporal 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.

4. Materials and Methods

Computation and analysis in this work were performed in Python (Van Rossum and Drake, 2009) with the use of the packages NumPy (Oliphant, 2006), SciPy (Virtanen et al., 2020), Pandas (McKinney, 2010), Jupyter (Pérez and Granger, 2007), Matplotlib (Hunter, 2007), Seaborn (Waskom et al., 2017), and NetworkX (Hagberg et al., 2018).

4.1. 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₀ = 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.

4.2. 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₀ 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.

4.3. Kymograph Construction for 1D Model

To obtain values of p_i(r, t) we integrate Equation (11). Starting from some initial colony radius R₀ we initialize p_i(r, 0) for r a regularly spaced lattice on (0, R₀). We use p₂(r, 0) = 5 for all r, that is homogeneous expression of only p₂. 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₀ + 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₀) 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^{\prime },t+\Delta 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.

4.4. 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.

4.5. 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).

4.6. 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.

4.7. 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}=\frac{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 $\lambda =(v_p+\frac{1}{2})T=(v_p+\frac{1}{2})\langle \Delta t\rangle$.

4.8. 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 $\bar{\mu }$. For each combination of parameters we simulated oscillations at the colony edge ($\bar{\mu }=1$) and the colony interior ($\bar{\mu }=0$), and measured the periods T_ext and T_int as described above. These values were then substituted into Equations (13)–(14) 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.

Author Contributions

GY, GV, and TR designed the study and analyzed the data. MM, GY, GV, and TR wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

GV was supported by a scholarship from the Institute for Biological and Medical Engineering, Pontificia Universidad Católica de Chile. GY was supported by Beca Ayudante Doctorando scholarship from the Department of Chemical and Bioprocess Engineering, Pontificia Universidad Católica de Chile. TR, GV, GY, and MM acknowledge financial support from the National Agency for Research and Development (ANID)/PIA/ACT192015.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Gustavo Düring, Luca Ciandrini, Pascal Rogalla, and Ignacio Medina for helpful and stimulating discussions. We also thank the members of the Synthetic Biology Lab for their support and encouragement—Anibal Arce, Kevin Simpson, Tamara Matute, Isaac Nuñez, Fernán Federici, among others. A preprint of this work is available on bioRxiv (Yáñez Feliu et al., 2020).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2020.00893/full#supplementary-material

References

Aboobaker, A. A., Tomancak, P., Patel, N., Rubin, G. M., and Lai, E. C. (2005). Drosophila microRNAs exhibit diverse spatial expression patterns during embryonic development. Proc. Natl. Acad. Sci. U.S.A. 102, 18017–18022. doi: 10.1073/pnas.0508823102

PubMed Abstract | CrossRef Full Text | Google Scholar

Alnahhas, R. N., Winkle, J. J., Hirning, A. J., Karamched, B., Ott, W., Josić, K., et al. (2019). Spatiotemporal dynamics of synthetic microbial consortia in microfluidic devices. ACS Synth. Biol. 8, 2051–2058. doi: 10.1021/acssynbio.9b00146

PubMed Abstract | CrossRef Full Text | Google Scholar

Andersen, K. B., and von Meyenburg, K. (1980). Are growth rates of Escherichia coli in batch cultures limited by respiration? J. Bacteriol. 144, 114–123. doi: 10.1128/JB.144.1.114-123.1980

PubMed Abstract | CrossRef Full Text | Google Scholar

Ausländer, D., Ausländer, S., Pierrat, X., Hellmann, L., Rachid, L., and Fussenegger, M. (2018). Programmable full-adder computations in communicating three-dimensional cell cultures. Nat. Methods 15:57. doi: 10.1038/nmeth.4505

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, C. J., Heisenberg, C. P., and Hiiragi, T. (2017). Coordination of morphogenesis and cell-fate specification in development. Curr. Biol. 27, R1024–R1035. doi: 10.1016/j.cub.2017.07.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Cookson, N. A., Mather, W. H., Danino, T., Mondragõn-Palomino, O., Williams, R. J., Tsimring, L. S., et al. (2011). Queueing up for enzymatic processing: correlated signaling through coupled degradation. Mol. Syst. Biol. 7, 1–9. doi: 10.1038/msb.2011.94

PubMed Abstract | CrossRef Full Text | Google Scholar

Danino, T., Mondragón-Palomino, O., Tsimring, L., and Hasty, J. (2010). A synchronized quorum of genetic clocks. Nature 463, 326–330. doi: 10.1038/nature08753

PubMed Abstract | CrossRef Full Text | Google Scholar

De Jong, H. (2002). Modeling and simulation of genetic regulatory systems: a literature review. J. Comput. Biol. 9, 67–103. doi: 10.1089/10665270252833208

PubMed Abstract | CrossRef Full Text | Google Scholar

Deneke, V. E., and Di Talia, S. (2018). Chemical waves in cell and developmental biology. J. Cell Biol. 217, 1193–1204. doi: 10.1083/jcb.201701158

PubMed Abstract | CrossRef Full Text | Google Scholar

Elowitz, M. B., and Leibler, S. (2000). A synthetic oscillatory network of transcriptional regulators. Nature 403, 335–338. doi: 10.1038/35002125

PubMed Abstract | CrossRef Full Text | Google Scholar

Endy, D., and Brent, R. (2001). Modelling cellular behaviour. Nature 409, 391–395. doi: 10.1038/35053181

PubMed Abstract | CrossRef Full Text | Google Scholar

Espeso, D. R., Martínez-García, E., de Lorenzo, V., and Goñi-Moreno, Á. (2016). Physical forces shape group identity of swimming Pseudomonas putida cells. Front. Microbiol. 7:1437. doi: 10.3389/fmicb.2016.01437

PubMed Abstract | CrossRef Full Text | Google Scholar

Farrell, F., Hallatschek, O., Marenduzzo, D., and Waclaw, B. (2013). Mechanically driven growth of quasi-two-dimensional microbial colonies. Phys. Rev. Lett. 111, 1–5. doi: 10.1103/PhysRevLett.111.168101

PubMed Abstract | CrossRef Full Text | Google Scholar

Fraleigh, S. P., and Bungay, H. R. (1986). Modelling of nutrient gradients in a bacterial colony. Microbiology 132, 2057–2060. doi: 10.1099/00221287-132-7-2057

CrossRef Full Text | Google Scholar

Gardner, T. S., Cantor, C. R., and Collins, J. J. (2000). Construction of a genetic toggle switch in Escherichia coli. Nature 403, 339–342. doi: 10.1038/35002131

PubMed Abstract | CrossRef Full Text | Google Scholar

Gjorevski, N., and Nelson, C. M. (2010). Endogenous patterns of mechanical stress are required for branching morphogenesis. Integr. Biol. 2, 424–434. doi: 10.1039/c0ib00040j

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorochowski, T. E. (2016). Agent-based modelling in synthetic biology. Essays Biochem. 60, 325–336. doi: 10.1042/EBC20160037

PubMed Abstract | CrossRef Full Text | Google Scholar

Grant, M. A., Wacław, B., Allen, R. J., and Cicuta, P. (2014). The role of mechanical forces in the planar-to-bulk transition in growing Escherichia coli microcolonies. J. R. Soc. Interface 11:20140400. doi: 10.1098/rsif.2014.0400

PubMed Abstract | CrossRef Full Text | Google Scholar

Green, A. A., Kim, J., Ma, D., Silver, P. A., Collins, J. J., and Yin, P. (2017). Complex cellular logic computation using ribocomputing devices. Nature 548, 117–121. doi: 10.1038/nature23271

PubMed Abstract | CrossRef Full Text | Google Scholar

Guélon, T., Mathias, J.-D., and Deffuant, G. (2012). Influence of spatial structure on effective nutrient diffusion in bacterial biofilms. J. Biol. Phys. 38, 573–588. doi: 10.1007/s10867-012-9272-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hagberg, A. A., Schult, D. A., and Swart, P. J. (2018). “Exploring network structure, dynamics, and function using networkx,” in Proceedings of the 7th Python in Science Conference, eds T. V. Gäel Varoquaux and J. Millman (Pasadena, CA), 11–15.

Google Scholar

Healy, C. P., and Deans, T. L. (2019). Genetic circuits to engineer tissues with alternative functions. J. Biol. Eng. 13, 1–7. doi: 10.1186/s13036-019-0170-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hilborn, R. C., and Erwin, J. D. (2008). Stochastic coherence in an oscillatory gene circuit model. J. Theor. Biol. 253, 349–354. doi: 10.1016/j.jtbi.2008.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Hopfield, J. J. (1974). Kinetic proofreading: a new mechanism for reducing errors in biosynthetic processes requiring high specificity. Proc. Natl. Acad. Sci. U.S.A. 71, 4135–4139. doi: 10.1073/pnas.71.10.4135

PubMed Abstract | CrossRef Full Text | Google Scholar

Hunter, J. D. (2007). Matplotlib: a 2D graphics environment. Comput. Sci. Eng. 9, 90–95. doi: 10.1109/MCSE.2007.55

CrossRef Full Text | Google Scholar

Johnson, M. B., March, A. R., and Morsut, L. (2017). Engineering multicellular systems: using synthetic biology to control tissue self-organization. Curr. Opin. Biomed. Eng. 4, 163–173. doi: 10.1016/j.cobme.2017.10.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Kantsler, V., McDonald, E. O., Kuey, C., Ghanshyam, M. J., and Asally, M. (2020). Pattern engineering of living bacterial colonies using meniscus-driven fluidic channels. 9, 1277–1283. doi: 10.1021/acssynbio.0c00146

PubMed Abstract | CrossRef Full Text | Google Scholar

Karig, D., Martini, K. M., Lu, T., DeLateur, N. A., Goldenfeld, N., and Weiss, R. (2018). Stochastic Turing patterns in a synthetic bacterial population. Proc. Natl. Acad. Sci. U.S.A. 115, 6572–6577. doi: 10.1073/pnas.1720770115

PubMed Abstract | CrossRef Full Text | Google Scholar

Khain, E., and Sander, L. M. (2006). Dynamics and pattern formation in invasive tumor growth. Phys. Rev. Lett. 96:188103. doi: 10.1103/PhysRevLett.96.188103

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, J., Yin, P., and Green, A. A. (2018). Ribocomputing: cellular logic computation using RNA devices. Biochemistry 57, 883–885. doi: 10.1021/acs.biochem.7b01072

PubMed Abstract | CrossRef Full Text | Google Scholar

Kimmel, C. B., Ballard, W. W., Kimmel, S. R., Ullmann, B., and Schilling, T. F. (1995). Stages of embryonic development of the zebrafish. Dev. Dyn. 203, 253–310. doi: 10.1002/aja.1002030302

PubMed Abstract | CrossRef Full Text | Google Scholar

Klumpp, S. (2011). Growth-rate dependence reveals design principles of plasmid copy number control. PLoS ONE 6:e20403. doi: 10.1371/journal.pone.0020403

PubMed Abstract | CrossRef Full Text | Google Scholar

Klumpp, S., Zhang, Z., and Hwa, T. (2009). Growth rate-dependent global effects on gene expression in bacteria. Cell 139, 1366–1375. doi: 10.1016/j.cell.2009.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Lestas, I., Vinnicombe, G., and Paulsson, J. (2010). Fundamental limits on the suppression of molecular fluctuations. Nature 467, 174–178. doi: 10.1038/nature09333

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Prindle, A., Humphries, J., Gabalda-Sagarra, M., Asally, M., Lee, D. Y. D., et al. (2015). Metabolic co-dependence gives rise to collective oscillations within biofilms. Nature 523, 550–554. doi: 10.1038/nature14660

PubMed Abstract | CrossRef Full Text | Google Scholar

Long, Z., Nugent, E., Javer, A., Cicuta, P., Sclavi, B., Lagomarsino, M. C., et al. (2013). Microfluidic chemostat for measuring single cell dynamics in bacteria. Lab Chip 13, 947–954. doi: 10.1039/c2lc41196b

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, N., Wang, S., and You, L. (2019). Synthetic pattern formation. Biochemistry 58, 1478–1483. doi: 10.1021/acs.biochem.8b01242

CrossRef Full Text | Google Scholar

Mangan, S., Zaslaver, A., and Alon, U. (2003). The coherent feedforward loop serves as a sign-sensitive delay element in transcription networks. J. Mol. Biol. 334, 197–204. doi: 10.1016/j.jmb.2003.09.049

PubMed Abstract | CrossRef Full Text | Google Scholar

Matson, J., and Characklis, W. G. (1976). Diffusion into microbial aggregates. Water Res. 10, 877–885. doi: 10.1016/0043-1354(76)90022-1

CrossRef Full Text | Google Scholar

Matsushita, M., and Fujikawa, H. (1990). Diffusion-limited growth in bacterial colony formation. Phys. A 168, 498–506. doi: 10.1016/0378-4371(90)90402-E

CrossRef Full Text | Google Scholar

McKinney, W. (2010). “Data structures for statistical computing in Python,” in Proceedings of the 9th Python in Science Conference, eds S. van der Walt and J. Millman (Austin, TX), 56–61.

Google Scholar

Murray, J. D. (2002). Mathematical Biology I. An Introduction, Vol. 17 of Interdisciplinary Applied Mathematics, 3 Edn. New York, NY: Springer.

PubMed Abstract | Google Scholar

Neubauer, P., Lin, H., and Mathiszik, B. (2003). Metabolic load of recombinant protein production: inhibition of cellular capacities for glucose uptake and respiration after induction of a heterologous gene in Escherichia coli. Biotechnol. Bioeng. 83, 53–64. doi: 10.1002/bit.10645

PubMed Abstract | CrossRef Full Text | Google Scholar

Niederholtmeyer, H., Sun, Z. Z., Hori, Y., Yeung, E., Verpoorte, A., Murray, R. M., et al. (2015). Rapid cell-free forward engineering of novel genetic ring oscillators. eLife 4, 1–18. doi: 10.7554/eLife.09771

PubMed Abstract | CrossRef Full Text | Google Scholar

Nielsen, A. A., Der, B. S., Shin, J., Vaidyanathan, P., Paralanov, V., Strychalski, E. A., et al. (2016). Genetic circuit design automation. Science 352:aac7341. doi: 10.1126/science.aac7341

PubMed Abstract | CrossRef Full Text | Google Scholar

Nuñez, I. N., Matute, T. F., Del Valle, I. D., Kan, A., Choksi, A., Endy, D., et al. (2017). Artificial symmetry-breaking for morphogenetic engineering bacterial colonies. ACS Synth. Biol. 6, 256–265. doi: 10.1021/acssynbio.6b00149

PubMed Abstract | CrossRef Full Text | Google Scholar

Oliphant, T. E. (2006). A Guide to NumPy, Vol. 1. Cambridge, MA: Trelgol Publishing.

Google Scholar

Osella, M., and Lagomarsino, M. C. (2013). Growth-rate-dependent dynamics of a bacterial genetic oscillator. Phys. Rev. E 87:012726. doi: 10.1103/PhysRevE.87.012726

PubMed Abstract | CrossRef Full Text | Google Scholar

Pérez, F., and Granger, B. E. (2007). IPython: a system for interactive scientific computing. Comput. Sci. Eng. 9, 21–29. doi: 10.1109/MCSE.2007.53

CrossRef Full Text | Google Scholar

Perez-Carrasco, R., Barnes, C. P., Schaerli, Y., Isalan, M., Briscoe, J., and Page, K. M. (2018). Combining a toggle switch and a repressilator within the AC-DC circuit generates distinct dynamical behaviors. Cell Syst. 6, 521–530.e3. doi: 10.1016/j.cels.2018.02.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Potvin-Trottier, L., Lord, N. D., Vinnicombe, G., and Paulsson, J. (2016). Synchronous long-term oscillations in a synthetic gene circuit. Nature 538, 514–517. doi: 10.1038/nature19841

PubMed Abstract | CrossRef Full Text | Google Scholar

Purcell, O., Grierson, C. S., Di Bernardo, M., and Savery, N. J. (2012). Temperature dependence of ssrA-tag mediated protein degradation. J. Biol. Eng. 6:10. doi: 10.1186/1754-1611-6-10

PubMed Abstract | CrossRef Full Text | Google Scholar

Riglar, D. T., Richmond, D. L., Potvin-Trottier, L., Verdegaal, A. A., Naydich, A. D., Bakshi, S., et al. (2019). Bacterial variability in the mammalian gut captured by a single-cell synthetic oscillator. Nat. Commun. 10, 1–12. doi: 10.1038/s41467-019-12638-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Rudge, T. J., Federici, F., Steiner, P. J., Kan, A., and Haseloff, J. (2013). Cell polarity-driven instability generates self-organized, fractal patterning of cell layers. ACS Synth. Biol. 2, 705–714. doi: 10.1021/sb400030p

PubMed Abstract | CrossRef Full Text | Google Scholar

Rudge, T. J., Steiner, P. J., Phillips, A., and Haseloff, J. (2012). Computational Modeling of Synthetic Microbial Biofilms. ACS Synth. Biol. 1, 345–352. doi: 10.1021/sb300031n

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruprecht, V., Monzo, P., Ravasio, A., Yue, Z., Makhija, E., Strale, P. O., et al. (2017). How cells respond to environmental cues - insights from bio-functionalized substrates. J. Cell Sci. 130, 51–61. doi: 10.1242/jcs.196162

PubMed Abstract | CrossRef Full Text | Google Scholar

Sagués, F., Sancho, J. M., and García-Ojalvo, J. (2007). Spatiotemporal order out of noise. Rev. Modern Phys. 79, 829–882. doi: 10.1103/RevModPhys.79.829

CrossRef Full Text | Google Scholar

Salis, H. M., Mirsky, E. A., and Voigt, C. A. (2009). Automated design of synthetic ribosome binding sites to control protein expression. Nat. Biotechnol. 27, 946–950. doi: 10.1038/nbt.1568

PubMed Abstract | CrossRef Full Text | Google Scholar

Santos-Moreno, J., and Schaerli, Y. (2019). Using synthetic biology to engineer spatial patterns. Adv. Biosyst. 3, 1–15. doi: 10.1002/adbi.201800280

PubMed Abstract | CrossRef Full Text | Google Scholar

Sickle, J. J., Ni, C., Shen, D., Wang, Z., Jin, M., and Lu, T. (2020). Integrative circuit-host modeling of a genetic switch in varying environments. Sci. Rep. 10:8383. doi: 10.1038/s41598-020-64921-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, W. P., Davit, Y., Osborne, J. M., Kim, W., Foster, K. R., and Pitt-Francis, J. M. (2017). Cell morphology drives spatial patterning in microbial communities. Proc. Natl. Acad. Sci. U.S.A. 114, E280–E286. doi: 10.1073/pnas.1613007114

PubMed Abstract | CrossRef Full Text | Google Scholar

Steiner, P. J., Williams, R. J., Hasty, J., and Tsimring, L. S. (2016). Criticality and adaptivity in enzymatic networks. Biophys. J. 111, 1078–1087. doi: 10.1016/j.bpj.2016.07.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Stricker, J., Cookson, S., Bennett, M. R., Mather, W. H., Tsimring, L. S., and Hasty, J. (2008). A fast, robust and tunable synthetic gene oscillator. Nature 456, 516–519. doi: 10.1038/nature07389

PubMed Abstract | CrossRef Full Text | Google Scholar

Tamsir, A., Tabor, J. J., and Voigt, C. A. (2011). Robust multicellular computing using genetically encoded nor gates and chemical “wires.” Nature 469, 212–215. doi: 10.1038/nature09565

PubMed Abstract | CrossRef Full Text | Google Scholar

Toda, S., Blauch, L. R., Tang, S. K., Morsut, L., and Lim, W. A. (2018). Programming self-organizing multicellular structures with synthetic cell-cell signaling. Science 361, 156–162. doi: 10.1126/science.aat0271

PubMed Abstract | CrossRef Full Text | Google Scholar

Tuson, H. H., Auer, G. K., Renner, L. D., Hasebe, M., Tropini, C., Salick, M., et al. (2012). Measuring the stiffness of bacterial cells from growth rates in hydrogels of tunable elasticity. Mol. Microbiol. 84, 874–891. doi: 10.1111/j.1365-2958.2012.08063.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Rossum, G., and Drake, F. L. (2009). Python 3 Reference Manual. Scotts Valley, CA: CreateSpace.

Google Scholar

Velardo, M. J., Burger, C., Williams, P. R., Baker, H. V., López, M. C., Mareci, T. H., et al. (2004). Patterns of gene expression reveal a temporally orchestrated wound healing response in the injured spinal cord. J. Neurosci. 24, 8562–8576. doi: 10.1523/JNEUROSCI.3316-04.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Vicsek, T., Cserző, M., and Horváth, V. K. (1990). Self-affine growth of bacterial colonies. Phys. A 167, 315–321. doi: 10.1016/0378-4371(90)90116-A

CrossRef Full Text | Google Scholar

Vilar, J. M., Kueh, H. Y., Barkai, N., and Leibler, S. (2002). Mechanisms of noise-resistance in genetic oscillators. Proc. Natl. Acad. Sci. U.S.A. 99, 5988–5992. doi: 10.1073/pnas.092133899

PubMed Abstract | CrossRef Full Text | Google Scholar

Vining, K. H., and Mooney, D. J. (2017). Mechanical forces direct stem cell behaviour in development and regeneration. Nat. Rev. Mol. Cell Biol. 18, 728–742. doi: 10.1038/nrm.2017.108

PubMed Abstract | CrossRef Full Text | Google Scholar

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Walt, S. J. V. D., et al. (2020). Scipy 1.0: fundamental algorithms for scientific computing in python. Nat. Methods. 17, 261–272. doi: 10.1038/s41592-019-0686-2

CrossRef Full Text | Google Scholar

Wang, W., Li, L., Ding, M., Luo, G., and Liang, Q. (2018). A microfluidic hydrogel chip with orthogonal dual gradients of matrix stiffness and oxygen for cytotoxicity test. Biochip J. 12, 93–101. doi: 10.1007/s13206-017-2202-z

CrossRef Full Text | Google Scholar

Waskom, M., Botvinnik, O., O'Kane, D., Hobson, P., Lukauskas, S., Gemperline, D. C., et al. (2017). Seaborn: v0.8.1. Stanford, CA.

Google Scholar

Winkle, J. J., Igoshin, O. A., Bennett, M. R., Josić, K., and Ott, W. (2018). Modeling mechanical interactions in growing populations of rod-shaped bacteria. Phys. Biol. 14:e110742. doi: 10.1101/110742

PubMed Abstract | CrossRef Full Text | Google Scholar

Wong, A., Wang, H., Poh, C. L., and Kitney, R. I. (2015). Layering genetic circuits to build a single cell, bacterial half adder. BMC Biol. 13:40. doi: 10.1186/s12915-015-0146-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Woods, M. L., Leon, M., Perez-Carrasco, R., and Barnes, C. P. (2016). A statistical approach reveals designs for the most robust stochastic gene oscillators. ACS Synth. Biol. 5, 459–470. doi: 10.1021/acssynbio.5b00179

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, M., and Fussenegger, M. (2018). Designing cell function: assembly of synthetic gene circuits for cell biology applications. Nat. Rev. Mol. Cell Biol. 19, 507–525. doi: 10.1038/s41580-018-0024-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Yáñez Feliu, G. A., Vidal Peña, G., Mũñoz Silva, M. A., and Rudge, T. J. (2020). Novel tunable spatio-temporal patterns from a simple genetic oscillator circuit. bioRxiv [Preprint]. doi: 10.1101/2020.07.06.190199

CrossRef Full Text | Google Scholar

Yeung, E., Dy, A. J., Martin, K. B., Ng, A. H., Del Vecchio, D., Beck, J. L., et al. (2017). Biophysical constraints arising from compositional context in synthetic gene networks. Cell Syst. 5, 11–24.e12. doi: 10.1016/j.cels.2017.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, T., Zhang, J., Yuan, Z., and Chen, L. (2008). Synchronization of genetic oscillators. Chaos 18:037126. doi: 10.1063/1.2978183

PubMed Abstract | CrossRef Full Text | Google Scholar