Viva in Silico: A position-based dynamics model for microcolony morphology simulation

We present a position-based dynamics model for microcolony growth. In addition to achieving fast and stable simulation of thousands of cells, this model allows for the computation of cell interaction with the environment without sacrificing robustness and predictability. We introduce state-ofthe-art principles of synthetic biology into our framework to enable biologically-informed microcolony pattern formation. We give detailed implementation of growth, communication, and external influences within the system and demonstrate our method for several scenarios which are experimentally verified. Finally, we provide a use case for rapid simulations enabled through our method including parameter search for tuning spoke-based pattern formation utilizing predefined and formulated biological primitives.


Introduction
Pattern formation is a prime example of complex synchronized single and multiple cell behavior. The reproducible emergence of patterns has garnered great interest in the field of synthetic biology, which uses engineering principles to design and build biological components, elucidating new insights on life (Liu et al., 2013). Research advances in synthetic pattern formation are increasingly enabling the ideation and implementation of practical applications associated with programmable morphologies. Within the field of engineered living materials, such applications include nanowire formation (Chen et al., 2015), self-healing construction materials (Wiktor and Jonkers, 2011), and multifunctional tissue constructs (Khademhosseini and Langer, 2016).
As cell-to-cell interactions become more complex and challenging to predict across spatial and temporal scales, biologically-informed computer aided design (bioCAD) tools to study and design cell colony behavior become indispensable. Cell colony behavior is frequently modeled by agent-based methods, where complex global colony characteristics arise from local cell behavior dened on the level of an individual cell. Popular implementations of this approach can be found in frameworks such as CellModeller (Rudge et al., 2012), GRO (Jang et al., 2012), and Biocellion (Kang et al., 2014).
We introduce an agent-based modeling approach using position-based dynamics (PBD) (Müller et al., 2007) for the simulation of local cell interactions. Traditional approaches using classical dynamics attempt to approximate variations in acceleration, velocity, and position caused by momentum using numerical integration as a function of applied forces. In contrast, PBD solves systems of position-based constraints in an iterative manner. The PBD approach has been applied to the simulation of both soft and rigid bodies (Deul et al., 2016), fluids (Macklin and Müller, 2013), and crowds (Weiss et al., 2017).
This results in two advantages over rigid-body or massspring systems commonly used in agent-based modeling for colony behavior. First, the robustness and speed of PBD allows for the simulation of large colonies of millions of cells. Second, simulations become controllable as the positions of cells can be changed without disturbing the stability of the system. This enables simple implementation of interactions of cells with an environment acting upon them.
In this paper we provide a detailed description of our method and subsequently show simulation results as well as related in vitro experiments. To conclude, we show a concrete example of pattern formation implemented with our method showcasing some of its advantages.

Method
We describe our model framework for the simulation of controllable pattern formation in microcolonies. We base it on real-world biological primitives used in the synthetic biology context encompassing growth, division, differentiation, communication through quorum signals, state-switching and cell interaction with the environment.

Growth
We model the cell population as a set of vertices V = {v 1 , ..., v n } where each vertex represents a single cell, or a collection of vertices represents a cell with a shape of the cell being determined by edges E = {e 1 , ..., e j } ∈ V × V . We name this collection of vertices, cells C = {c 1 , ..., c m }, where each element c i is either a single vertex or deter-mined by a surjective mapping C(v i ): V → C, where T = C −1 (c i ) are connect by edges in E. The cells C can be associated with additional properties such as expression rate of signals, cell type, division rate (doubling per hour), mass, velocity, size etc. We collect these as a set of state variables S = {s 1 , ..., s m }.
The goal of the simulation framework is to predict behavior of colonies containing millions of cells over long periods of time while still observing cell behavior at the single cell level. Additionally, simulations should produce patterns fast and stable, in order to perform parameter searches with numerous initial and boundary conditions of the system. Such a system will enable users to guide circuit engineering for final in vitro pattern formation. To achieve this, we model cell dynamics, growth, and differentiation in a positionbased dynamics framework (Algorithm 1). PBD simulations solve systems of non-linear constraints updating position directly in addition to working with a force-based position estimate. Accordingly, cell interactions during growth are mostly modeled here by constraints acting on a neighbourhood of cells. We use contact, friction, cohesion, cell division, and inhibition to mimic the dynamics of an evolving colony. Cell contact is resolved through the distance con- where the position and size of v i are denoted by x i and r i respectively. Cell-to-cell kinematic and static friction are adopted from particle-based schemes (Macklin et al., 2014). To achieve cell cohesion, an XSPH viscosity model is adapted to update vertex velocities (Macklin and Müller, 2013). Cell division occurs at a cell-specific rate and is modeled by adding vertices to the system.
In the case of single vertex cells, vertices are inserted at an offset to the cell position in the direction of a position gradient of the local neighborhood. In the case a cell shape is used for c i , vertices are inserted to T at the bisection of edges exceeding a cell-specific length, and a new cell c i+1+m is formed if the cells total edge-lengths exceed a threshold. In both cases, the cell size is increased over time until the threshold is reached. Details on the modified PBD algorithm can be found in Algorithm 1.

Communication
Communication in and across colonies is mediated by secretion of specific molecules and subsequent response when a threshold concentration is reached. To enable cell-to-cell signaling, the substrate is modeled as a regular grid on which signals A n (chemical inducers) are modeled as fluids. In particular, every grid cell queries the cell population for cells which occupy it. The cells carry associated production rates for signals depending on their type and their respective rate is incorporated in a source term s An (x). Diffusion and sourcing is then modeled through the partial differential equation and discretized by a finite difference scheme. Here, D n is the diffusion coefficient of A n and specific to the signal. Furthermore, to enable quorum sensing, cells sample the grids and associated concentration of the signals and values are interpolated bilinearly. Obtained signaling information can be used to make decisions on the cell state.

States
In a recent progression, synthetic finite state machines have been established using combinations of integrases, enzymes that insert or excise DNA from a bacterial genome, to record state in a DNA register and further create phenotypic expressions for each differentiated cell state (Roquet et al., 2016). The Recombinase-based State Machine (RSM) is thus in-dicative of the ability to program cells to carry out higher level serial operations and cellular computation. Furthermore, a bacterial cell may represent both virtual and physical properties-achieving not only a logical progression, but also a change of chemical and physical signals, in relation to an environment, and spatial positioning relative to other cells. For this reason, the state variable S in our framework carries information about cell behavior in the system and is changed by the sampled quorum signals. direct approach for our system to interact with other simulations. In our case, since the signaling process is already modeled on a grid, we can model the advection of signals by an external velocity field and solving a simplified fluid model. We then model the advection of cells using a fluid velocity field. Similar to fluid implicit particle methods this allows one to combine simulation of cells in parallel to the potential mechanics of fluids inside or over the substrate. Contacts and collisions with external geometries can be resolved in a similar fashion, by updating position estimates.

Experiments and Results
We demonstrate our PBD framework through a variety of scenarios representing growth, cell-to-cell communication, and cell-environment interaction.

Growth
We demonstrate how distinct global colony morphologies can be observed from local cell shapes in Figure 1. Cell shape, as shown above, is modeled as a set of connected edges. Initial cells are colored and color information is passed along during division. Figure 1 (A) shows globular cell shapes with diameters of up to 2 µm. Figure 1 (B) shows rod-like cell shapes with lengths of up to 3 µm. Figure 1 (C) shows rod-like cell shapes up to 10 µm in length. Distinct global colony morphologies can be observed from local cell shapes. While globular cells result in uniform spoke-like propagation of cell types from the center outwards, long strand-like cell shapes result in local curling and a fuzzy boundary of the emerging colony. Unlike other methods, our approach is not limited to flat environments. When gravity is incorporated in the position estimate in the PBD scheme, typical 3D cell colony morphologies start to emerge from interaction of cell division and inhibition. Figure 2 shows the development of cell morphologies in R 3 . Initial cells are colored, and color information is passed along during division. Figure 2 (A) shows cells no longer constrained to R 2 and growing on a flat substrate plane under gravity. Here, cell growth is inhibited by the gravity-induced strain on the cell and through the contact with the substrate plane. In Figure 2 (B), gravity and substrate plane are removed and cells divide freely, yet are constrained by the dampening in the system. In Figure 2 (B) (lower right) cell clusters originating from two individual cells are isolated and show spherical cone morphology. This type of simulation is particularly important in the formation of potential arrangements of cells in 3D for living materials. We give simulation durations in Figure 2 in Table 1 for up 100000 cells. All measurements were taken on an Intel Xeon E5-1650v3 at 3.50 GhZ, with 64GB RAM and NVIDIA Quadro M4000 graphics card.  Figure 3 (B) an example of quorum signaling with a recombinant circuit is given. Here, the initial cell type (white) expresses signals A and B. Upon reaching a threshold concentration, cells switch to another cell type amplifying the production of the received signal causing the state-switch. This demonstrates a feedback loop resulting in a stable spoke pattern.

Environment
Examples of interaction of the growth model with external influences are shown in Figure 4. Figure 4 (A) presents the mixing of cells and signals within a substrate hardening over time (e.g. cooling of agar over time). From left to right, we show two cell types (A and B), signal A, signal B, and the velocity field of the system. Signal A induces cell-type A   Figure 2.
and signal B induces cell-type B (red and yellow). Initially, cells are randomly distributed (leftmost image). Cells and signals are mixed by the velocity field, which vanishes over time due to increasing viscosity of the substrate, resulting in a final stabilized pattern. (B) Shows the interaction of a spoke forming colony with an external object. Here, we show that inhibition in our system is temporary and cause by contact inhibition of growth. Cells that regain free space and therefore have reduced internal strains are no longer inhibited. This results in the restoration of the global colony morphology.

In Vitro Experiments
Two validations of our simulations through experimental results are given in Figure 5. Figure 5 (A) confirms the results of Figure 4 (A), illustrating the mixing of two recombinant Escherichia coli (DH5 Pro) populations within an agar suspension. The heterogeneous distribution of the two populations is visualized by the expression of either green or red fluorescent protein (GFP or RFP respectively), and was imaged using a Canon 5D camera and a transilluminator (SafeLight, Invitrogen). Figure 5 (B) compares in silico and in vitro results regarding the E. coli colony height (elevation view). In silico, we simulated the growth of a colony-similar to Figure 2 (A)-of circa 1 million cells. In vitro, cells from an overnight culture were grown on Luria-Bertrani (LB) agar plates for 24 hours at 37 • C, and subsequently imaged with a Nikon D3300 camera. Figure 5 (C) compares in silico and in vitro results of Bacillus subtilis (ATCC 6051) growth on a soft agar substrate. Decreased substrate stiffness is correlated with increased cell length and sliding motility (Fall et al., 2006). Cells were grown for 24 hours at 30 • C on 1.0% agar supplemented with nutrient broth and subsequently stained with propidium iodide and SYTO 9 (Invitrogen). Samples were imaged using a Leica TCS SP8 confocal system. The in silico model used was similar to Figure 1 (C), with lower friction between substrate and cell.

Case Study
The speed and robustness of our system make it valuable for parameter searches on pattern forming engineered systems. We leverage our framework for initial predictions of controlled spoke formation in bacteria. We model low basal levels of recombinase expression to control switching frequency and repression of recombinase expression upon switching. This is done through quorum signaling of already differentiated cells repressing recombinase expression thus limiting switching of nearby cells.
We parameterize the simulation with three variables: 1) Switching frequency -the probability of cells switching from their undifferentiated state to a differentiated state, 2) Repressor expression -the rate at which the repressor is expressed by the cells, 3) A random initial configuration of cells. Results are shown in Figure 6. In Figure 6 from left to right, switching frequency is decreased and from top to bottom, repressor expression increased. Sixteen simulations are given in panels, where the top shows colony morphology and the bottom shows repressor concentration. Generally, switching frequency allows control over the frequency of spokes, while repressor expression rates allow localiza- Figure 6: Examples of a parameter search on microcolony spoke formation. From left to right we vary switching frequency; from top to bottom, we change repressor expression. This results in control of frequency of spokes and localization of spokes respectively. tion of spokes. Clearly, switching frequency allows for approximate control over number of spokes formed, and repressor expression allows for tuning of localization of the spokes. However, the dynamics of the system are more nuanced. As an example, high switching frequency and low repressor expression still result in high repressor concentration in the substrate due to delayed In Figure 6 from left to right, switching frequency is decreased and from top to bottom, repressor expression increased. 16 simulations are given in panels, where the top shows colony morphology and the bottom shows repressor concentration. Generally, switching frequency allows control over the frequency of spokes, while repressor expression rates allow localization of spokes. Clearly, switching frequency allows for approximate control over number of spokes formed, and repressor expression allows for tuning of localization of the spokes. However, the dynamics of the system are more nuanced. As an example, high switching frequency and low repressor expression still result in high repressor concentration in the substrate due to delayed repression in the system and the ability of more cells to switch. Additional parameter tuning can help to elucidate the system's response, and may be addressed in future work.

Conclusions and Future Work
In this paper, we introduced a position-based dynamics approach to microcolony morphology simulation. Compared to existing methods the advantages of PBD are large time steps, guaranteed stability, and ease of control. Drawbacks of our method as with all PBD systems is the lack of relation of the model to real world parameters. While the simulations shown herein can be used to predict general behaviour and interactions of cells, we do not aim to simulate real world cell behaviour in its entirety. Additional constraints may be introduced in the growth model such as bending stiffness to extend the range of potential organisms and behaviour. Smart parameter tuning of the simulation model, enabled by the large number of experimental data obtainable within a short amount of time may help to align in silico experiments with in vitro observations.