Repulsive expansion dynamics in colony growth and patterning

Spatial expansion of a population of cells can arise from growth of microorganisms, plant cells, and mammalian cells. It underlies normal or dysfunctional tissue development and it can be exploited as the foundation for programming spatial patterns. This expansion is often driven by continuous growth and division of cells within a colony, which in turn pushes the peripheral cells outward. This process generates a repulsion velocity field at each location within the colony. Here we show that this process can be approximated as coarse-grained repulsive-expansion kinetics. This framework enables accurate and efficient simulation of growth and gene expression dynamics in radially symmetric colonies with homogenous z-directional distribution. It is robust even if cells are not spherical and vary in size. The simplicity of the resulting mathematical framework also greatly facilitates generation of mechanistic insights.


Introduction
Spatial expansion of a population of cells is ubiquitous in biology. It can arise from growth of bacterial or yeast colonies (1,2), development of plant tissues (3), animal tissues (4), or growth of tumors (5)(6)(7). Understanding the expansion dynamics of a cell population is important, since colony expansion is a reflection of many cellular properties, including the substrate nutrient uptake rate, the metabolism, the division and death rates of individual cells. In addition, since colony expansion is usually coupled with the gene expression activity and the metabolic activity, understanding the expansion dynamics is the foundation of understanding the spatial dynamics of gene expression in growing colonies (8,9).
To date, primarily two types of models have been used to simulate spatial dynamics in growing colonies (including gene expression): agent-based models (ABMs) and continuum models based on partial differential equations (PDEs). An ABM focuses on simulating a single cell or a cluster of cells as an agent; the behavior of each agent is constrained by a set of pre-defined rules. The collective behavior of the population emerges from the interactions between agents, through contact or diffusible chemicals. Using the ABM entails assigning pre-defined rules and parameters to each agent.
It provides detailed information about individual agents, as well as the overall system. However, it has several potential limitations. First, certain ABMs require more extensive assumptions about the definition of each agent, how each behaves, and how different agents interact (10) - (11). Many of these assumptions are made based on self-consistency and require parameters that are difficult or impossible to determine independently. Second, by construction, the computational demand by an ABM scales with the number of agents, as well as the complexity of reaction kinetics in each agent. As a result, ABMs are typically computationally expensive, making it difficult or even impossible to simulate the spatiotemporal dynamics of a cell population containing tens of millions of cells (a visible bacterial colony contains about 10-100 million cells). Finally, because of the sheer complexity of inter-agent 4 interactions and the computational cost, it is difficult to deduce intuitive understanding of the properties emerging from the overall system (12).
These limitations can be partially alleviated by using continuum models consisting of PDEs.
However, numerically solving PDEs can still be computationally demanding (though typically less so than ABMs) (13). Even though these tools can generate simulation results that can recapture certain aspects of the experimental results, it is difficult to develop intuitive understandings from the simulations.
In general, the colony expansion can be driven by a combination of internal and external forces.
In many cases, however, the expansion is primarily driven by continuous growth and division of cells in the colony (25). That is, the cells in the interior of a colony collectively push the peripheral cells outward, leading to colony expansion. The driver is the mechanical force generated by growth. This notion has been well recognized in diverse organisms, including bacteria (26,27), yeast (28), plant cells (29,30), and tumor cells (31,32). Recognizing this feature enables us to model the colony expansion and gene expression using a highly simplified framework. In contrast to typical dimensionreduction of a complex model, the repulsive-expansion framework represents a completely different view in examining the spatiotemporal dynamics of expanding cell populations and associated gene expression. In particular, our model directly uses an ODE framework to track the "flow" of individual cells as they experience the steric forces generated due to cell growth and division. By reducing the structural complexity of the model, this framework facilitates development of mechanistic insights into the underlying dynamics of interest.

The repulsive expansion model captures the colony growth dynamics
To develop our framework, we first assume the steric force between cell-cell growth and division are the driving force of colony expansion. At a particular location, the local steric force will generate a steric repulsion velocity field. We then assume that cells are spheres with the same radius.
Therefore, on average, the overall pushing force on each cell is perpendicular to the boundary of a radial symmetry colony. If the profile of this velocity field is known with given initial conditions, one can use sets of ODEs to simulate the moving trajectories of any objects within the colony ( Figure 1A).
In a closely packed colony, the cell division rate depends on the nutrient concentration and the location within the colony. That is, at the same nutrient concentration, if a cell is closer to the colony center, the slower the division rate is (33,34). The rate of change in the colony volume can either be expressed as a function of the normal velocity of the colony boundary or the total number of the cells. In our framework, the velocity field of any given position (within a colony) is expressed as a function of cell division rate and the distance between the cellular location and the colony center (see Supplementary Material for detailed derivations). Knowing this velocity field, in turn, enables the computation of the movement of each cell.
Specifically, numerically solving the governing equations (Eq. 9, see Supplementary Table 1 for parameter values) can generate the dynamics of colony expansion as a function of time. Figure 1B shows a typical set of results on colony radius expansion (solid black line). If we choose a specific location within the colony and imagine there is a cell, we can plot how this cell moves based on the velocity field. Therefore, the moving trajectories of all cells in the colony can be tracked during the expansion (colored lines). This notion could be counter-intuitive if a cell does not exist initially and 6 thus does not have a location in the initial colony. For example, consider a cell born at location " and at time " . If a velocity field trajectory passes this position ( " , " ), the trajectory would predict how the cell would move at any time point later than " . The section of the trajectory before " would be imaginary. However, to simplify presentation, all the trajectories are plotted from time zero in our simulation.
To evaluate the validity of the fundamental conceptual framework of the repulsive expansion model, we used an ABM (35) to simulate colony growth with similar physical constraints. Briefly, cells are modeled as elastic rods that increase in length as they consume nutrients. Cellular growth is modeled by a Monod equation (36) depending on nutrient availability. The total length of cell division is a truncated normal distribution, with values taken over a range of values (Supplementary Table 2).
Repulsive growth is maintained by using a soft particle technique (37), which models cellular interactions as those between overlapping elastic spheres. ABMs were conducted to simulate the growth profiles of a single cell growing into 10,000 cells. The cell trajectories generated by an ABM simulation ( Figure 1D) are qualitatively the same as those generated from the repulsive expansion model. In the ABM simulation, however, each trajectory only starts when a cell is born and does not have a hypothetic section leading to a position at time 0.
In the repulsive expansion model, we do not have cell shape or cell size as basic parameters.
We used an ABM to validate that the colony expansion size will not be affected by cell shape or size ( Figure 2A). While the repulsive expansion model is derived assuming spherical cells, it remains to be reliable even when this assumption is relaxed: the colony expansion dynamic from repulsive expansion model is comparable to the results solved from PDEs or ABM. For example, using repulsive expansion model (see Supplementary Table 3 for parameter values), regardless of the initial colony size, with the same environmental condition, the final colony sizes are the same ( Figure 2B). 7 An important advantage of the repulsive expansion model is that the accuracy of computation does not depend on the resolution of the spatial discretization in initiating the simulation. This property drastically increases the computation efficiency without sacrificing computational accuracy.
Despite differences in the segmentation of the spatial meshes, given the same initial colony size, the moving trajectories within the colony from different initial configurations are the same ( Figure 2C).
This property does not occur in the ABMs or PDE models, where a sufficiently high resolution in discretizing the space is critical to ensure computational accuracy or reliability. However, in the repulsive expansion model, a cell trajectory is directly determined by cell growth rate (Eq 9). Once the cell growth rate is given, the trajectory of each position could be expressed in an analytical solution, which does not require any spatial discretization. This feature will substantially simplify the simulation of cell movements.

Modeling programmed pattern formation dynamics using the repulsive expansion model
The repulsive expansion model can be readily extended to describe gene expression dynamics that are coupled with colony expansion. To demonstrate this practice, we apply this framework to the analysis of a synthetic pattern-formation circuit that we recently engineered (38,39). The circuit consists of a T7 RNA polymerase (T7) that activates its expression. Upon activation by T7, synthesis of AHL (A) will be mediated, which can diffuse across the cell membrane. When the global AHL concentration surpasses a threshold, intracellular AHL will trigger the activation of the synthesis of T7 lysozyme (L). Lysozyme then binds to the T7 and forms a T7-lysozyme complex (P), therefore inhibiting the T7 binding to the T7 promoter. This T7-lysozyme complex also inhibits T7 transcription.
In this process, the AHL concentration is affected by its initial concentration and the domain size. cases, the models are assumed to be radial symmetry. One study is in a one-dimensional ABM (39) and one used a PDE model (38). case, two lysozyme rings can form given a high initial AHL concentration, which is consistent with the observations in Payne et al. (Figure 3D).
In the study by Cao et al., with an initial AHL concentration, a core and a ring with narrower width formed near the edge of the colony compared to the base case (without initial AHL). Given an initial larger domain, the repulsive expansion model indeed generates a wider core and ring will form in a larger colony ( Figure 4C). Moreover, the simple model also successfully generates the scale invariance reported in Cao et al: the ring width and the colony radius are both proportional to the domain radius, when the latter is of moderate magnitude ( Figure 4D).
In addition to the superior computational efficiency, the simplicity of the modeling framework provides more intuitive insights of circuit dynamics. For example, the dynamics of T7 is observed to be faster and more transient than lysozyme's ( Figure 4B). In Cao et al., the gene expression dynamics and the cellular movements are coupled in the PDEs. However, in the repulsive expansion model, particularly when the metabolic burden is small, the colony expansion could be modeled as independent of circuit dynamics. By separating the cellular movement from gene expression, we can simplify the model into two layers: the trajectories of all cells that reflect cellular movement, and the gene expression dynamic along each trajectory.
Since the repulsive expansion model is simpler and has fewer parameters, it provides a much clearer framework to reveal the essential properties of the system's dynamics. By choosing an initial , when reaches equilibrium, we can draw the phase diagram of T7 and lysozyme. When T7 and lysozyme all start with small initial values, T7 will increase to a peak then decrease, while the significant accumulation of lysozyme happens after the peak of the T7 ( Figure 4E). By combining these gene expression dynamics with the trajectories of all cells within the colony, we can see this correlation between the T7 patterns and lysozyme patterns. Compared with the ABM or the PDE model, the 10 repulsive expansion model provides more intuition when it is used to analyze the dynamics of spatial patterns. In this case, the level of metabolic burden and gene expression capacity are two important factors to generate different types of patterns, which indicates the temporal cooperation between colony growth and gene expression are crucial for forming different types of patterns.
Since the repulsive expansion model describes cell moving trajectory and gene expression separately, this framework also provides direct insights into the relationship between cellular movement and circuit dynamics. In the previous PDE model, there are in total 22 parameters, which is very high dimensional. To reduce the complexity, the parameters that dictate to colony growth were fitted to experimental data. Cao et al. concluded the T7 lysozyme profile is mainly determined by circuit logic and growth dilution. During the colony expansion, near the colony edge, the lysozyme is insufficient to overcome the dilution of the cell growth. Therefore, a strong lysozyme core occurred before the ring formation. We agree the parameters being tested in the repulsive expansion model are a subset of the PDE and ABM's parameters. Repulsive expansion model is not just a one-dimensional reduction of complicated models, but a model with core physic assumption so that redundant and trivial parameters can be eliminated for the researchers to quickly identify the key factors in affecting circuit dynamics.
The previous study has focused on the patterning process when the cell division and expansion rates are high. Given the simplicity of the repulsive expansion model, we can readily test the effects of the cell division and expansion rate when they are small. Interestingly, our new modeling result predicts the core-ring pattern formation with the ring forming before the core ( Figure 4F). Lysozyme accumulates near the colony edge due to the high gene expression capacity at the location, which leads to ring formation. However, due to the circuit dynamic, the accumulating T7 and AHL lead to lysozyme increasing near the center of the colony. This patterning process of outer ring developing earlier than the inner core is analogous to the typical eyespot pattern-formation process in butterfly wings (40,41). On the background of wings, there are a layer of parafocal elements (PFEs), which serve as the pattern units. The eyespot patterns are the PFEs with different color distributions. Many studies have shown the eyespot signals that form inner and outer rings are released at different timepoints. The outer ring forms before the inner ring ( Figure 4G) (42,43). The network and the parameter combination identified by the repulsive expansion model could provide insights of how eyespot patterns are formed in butterfly wings.

Discussion
In synthetic biology, tremendous progress has been made in programming single cells or cell populations by using gene circuits. The majority of these are done in a homogenous environment, e.g.
in a liquid culture. Examples include generation of logic operations (44,45), bistability (46,47), and autonomous oscillations (48)(49)(50). Very few gene circuits have been successfully engineered to program self-organized spatial patterns (9,38,39,(51)(52)(53). The lack of progress in this direction is in major part due to the challenge associated with the design of gene circuits that could potentially generate desired function. This challenge, in turn, is due to the intrinsic difficulty in modeling spatial temporal dynamics of cell collectives involving cell growth and gene expression.
This difficulty is not unique to the design of synthetic gene circuits. Modeling analysis of colony expansion and patterning has a long history in endeavors to understand growth of bacterial physiology(1, 2), biofilm formation (54,55), and tumor expansion (5). In general, spatiotemporal dynamics of colony expansion and gene expression (of signaling) are modeled by using ABMs or reaction-diffusion models consisting of PDEs. Both modeling frameworks suffer from two major limitations: high computational cost and difficulty in drawing intuition. Extensive ABMs have been developed to simulate tight-packed expanding colonies consisting of no motile bacteria (56)(57)(58). A one-dimensional molecular dynamic model was used to compare with the ABM model in one of the study (56), by directly accounting for the elastic and friction forces experienced by the cells. Our repulsive-expansion model represents another layer of simplification by focusing on the velocity field as a reflection of contact-based forces between cells.
Our work demonstrates a highly simplified yet accurate modeling framework for a particular class of dynamics involving colony growth and patterning. In essence, the repulsive expansion framework models the "flow" of cells in a growing colony, driving by the force generated by collective growth and division of cells in the colony. One can imagine a perfect radial symmetric colony, given the nutrient concentration and distribution, the colony growth rate of each cell at any position is known. Therefore, the trajectory of any cell with a given initial position can be calculated based on the force the cell experiences. This framework is structurally different from the corresponding PDE model and it does not arise from typical dimensional reduction. Due to its structural simplicity, this modeling framework offers superior computational efficiency as it does not require high resolution of space discretization to ensure computational accuracy. The framework is readily amendable to the modeling of spatial temporal dynamics coupled with or arising from colony growth and expansion. We have illustrated this point by applying the modeling framework to a synthetic gene circuit that has been previously analyzed.
In general, the framework could be applied to a class of spatial distribution problems: the system is radially symmetric, with homogenous z-directional distribution, and with gradient-free chemicals or chemicals that have gradient but exceeds the triggering threshold. These criteria are the cornerstone of the repulsive expansion dynamics. Many examples of biological processes indeed satisfy these conditions. Examples include colony growth and gene expression in a microfluidic chip with a confined chamber height (59-61), a single layer of cell divisions in a microscope slide (62-64), 13 and cells in 3D symmetric spherical growth condition (65,66). Even if these systems can secrete diffusive chemicals as global signaling molecules, our modeling framework is applicable if signaling molecules' diffusion rates are sufficiently fast to form uniform spatial distributions (38,39), or if the molecules accumulate fast enough to exceeds beyond the gene expression triggering threshold (67,68).   The cell growth rate is a decaying function of the production of T7 RNAP and T7 lysozyme.
Right: Gene expression capacity. The x-axis represents the distance from the colony center.
B. Top left to bottom right: Simulated spatial-temporal dynamics of colony radius, AHL, T7 RNAP and T7 lysozyme for varying distance (x-axis) over time (y-axis), respectively. The parameters used in the simulation are listed in Table 4.  A. Simulation assumption with a low metabolic burden and a moderate gene expression capacity. The notations are the same as those in Figure 3.
B. Simulated spatial-temporal dynamics of T7 RNAP and T7 lysozyme for varying distance (x-axis) over time (y-axis). The parameters used in the simulation are listed in Table 4 18 E. Phase diagram of T7 lysozyme (y axis) to T7 RNAP (x axis). The red line shows given an initial condition (T7RNAP=0.1, T7 lysozyme =0), T7 RNAP increases first then decreases, while T7 lysozyme keeps increasing.
F. Simulated spatial-temporal dynamics of T7 RNAP and T7 lysozyme for varying distance (x-axis) over time (y-axis). The parameters used in the simulation are listed in the Table 4