Bridging the gap between single-cell migration and collective dynamics

Motivated by the wealth of experimental data recently available, we present a cellular-automaton-based modeling framework focussing on high-level cell functions and their concerted effect on cellular migration patterns. Specifically, we formulate a coarse-grained description of cell polarity through self-regulated actin organization and its response to mechanical cues. Furthermore, we address the impact of cell adhesion on collective migration in cell cohorts. The model faithfully reproduces typical cell shapes and movements down to the level of single cells, yet allows for the efficient simulation of confluent tissues. In confined circular geometries, we find that specific properties of individual cells (polarizability; contractility) influence the emerging collective motion of small cell cohorts. Finally, we study the properties of expanding cellular monolayers (front morphology; stress and velocity distributions) at the level of extended tissues.


Introduction
Cell movements range from uncoordinated ruffling of cell boundaries to the migration of single cells (Ridley et al., 2003) to the collective motions of cohesive cell groups (Friedl and Gilmour, 2009). Single-cell migration enables cells to move towards and between tissue compartments -a process that plays an important role in the inflammation-induced migration of leukocytes (Friedl and Weigelin, 2008). One can distinguish between amoeboid and mesenchymal migration, which are characterized by widely different cell morphologies and adhesive interactions with their respective environments (Friedl, 2004;Lämmermann and Sixt, 2009). Cells may also form cohesive clusters and mobilize as a collective (Trepat et al., 2009;Angelini et al., 2011;Doxzen et al., 2013;Deforet et al., 2014;Vedula et al., 2012;Marel et al., 2014). This last mode of cell migration is known to drive tissue remodelling during embryonic morphogenesis (Lecaudey and Gilmour, 2006) and wound repair (Poujade et al., 2007).
Despite this broad diversity of migration modes, there appears to be a general consensus that all require (to varying degrees) the following factors: (i) Cell polarization, cytoskeletal (re)organization, and force generation driven by the interplay between actin polymerization and contraction of actomyosin networks. (ii) Cell-cell cohesion and coupling mediated by adherens-junction proteins which are coupled to the cytoskeleton. (iii) Guidance by chemical and physical signals. The basic functionalities implemented by these different factors confer on cells the ability to generate forces, adhere (differentially) to each other and to a substrate, and respond to mechanical and chemical signals. However, a fully mechanistic understanding of how these basic functionalities are integrated into single-cell migration and coordinated multicellular movement is still lacking.
Here, we present a computational model which enables us to study cell migration at various scales, and thus provides an integrative perspective on the basic cell functions that enable the emergence of collective cell migration. While a variety of very successful modeling approaches has been used to describe single-cell dynamics (Mogilner, 2009;Marée et al., 2006;Marée et al., 2012;Shao et al., 2010;Ziebert et al., 2012;Ziebert and Aranson, 2013;Camley et al., 2013;Albert and Schwarz, 2014;Dietrich et al., 2018;Goychuk et al., 2018) or the movements of extended tissues (Szabó et al., 2006;Szabó et al., 2010;Kabla, 2012;Sepúlveda et al., 2013;Basan et al., 2013;Banerjee et al., 2015;Alt et al., 2017;Tarle et al., 2017), these models are hard to reconcile with each other. Models that focus on single cells are typically difficult to extend to larger cell numbers, largely due to their computational complexity. On the other hand, approaches which are designed to capture the dynamics at the scale of entire tissues generally adopt a rather coarse-grained point of view, and are therefore difficult to transfer to single cells or small cell cohorts. At present there are two partly competing and partly complementary approaches to bridge the gap between single-cell migration and collective dynamics, namely phase-field models (Shao et al., 2010;Ziebert et al., 2012;Shao et al., 2012;Camley et al., 2014;Camley and Rappel, 2014;Lö ber et al., 2015), and cellular Potts models (CPMs) (Szabó et al., 2010;Kabla, 2012;Szabó and Merks, 2013;van Oers et al., 2014;Segerer et al., 2015;Niculescu et al., 2015;Albert and Schwarz, 2016;Rens and Merks, 2017) first introduced by Graner and Glazier (1992).

Box 1. A simple description of complex cells?
Mammalian cells are made up of around 10 9 interacting proteins (Milo and Phillips, 2015) in an aqueous compartment enclosed by a lipid bilayer membrane. A substantial fraction of these proteins is devoted to the structural support of the cell. The cytoskeletal systems that perform this function also mediate elastic deformations of the cell through stresses induced by motor proteins. Cell migration is enabled by transient, transmembrane attachment of the cytoskeleton to external structures (extracellular matrix or a substrate) via integrins, and regulated by various signaling pathways. To gain insights into such a complex system, we simplify these networks, each comprised of many interacting components, into coarse building blocks, which might seem arbitrary at first, but serve to qualitatively capture generic features of the underlying machinery. These generic and qualitative building blocks allow us to finally arrive at a quantitative description of cell dynamics.
Building on and generalizing the CPM (Graner and Glazier, 1992), we present a cellular automaton model that is designed to capture essential cellular features even in the context of the migration of single cells and of small sets of cells. At the same time, it is computationally efficient for simulations with very large cell numbers (currently up to Oð10 4 Þ cells), thus permitting investigations of collective dynamics at the scale of tissues. Our model reproduces the most pertinent features of cell migration even in the limiting case of solitary cells, and is compatible with a wealth of experimental evidence derived from both small cell groups and larger collectives made up of several thousand cells. Specifically, by studying the characteristics of single-cell trajectories and of small cell groups confined to circular territories, we demonstrate that persistency of movements is significantly affected by cell stiffness and cell polarizability. Moreover, we investigate the dynamics of tissues in the context of a typical wound-healing assay (Poujade et al., 2007;Trepat et al., 2009;Serra-Picamal et al., 2012), and show that the model exhibits the recurring mechanical waves observed experimentally (Serra-Picamal et al., 2012), a feature which we attribute to the coupling between cell-sheet expansion and cell-density-induced growth inhibition.

Computational model Model geometry
We consider cells that adhere to a two-dimensional surface, spanned by the coordinates ðx; yÞ, through some contact area ( Figure 1A). Membrane protrusions and retractions, which determine cell motion and shape (Pollard and Borisy, 2003;Lauffenburger and Horwitz, 1996), correspond to size and shape changes of the surface contact area. We assume that processes that take place at the cell boundary drive cell motion, and therefore disregard the cell body, which extends into the zdirection. In our computational model, we tesselate the available surface into a honeycomb lattice, where each hexagon corresponds to a discrete adhesion between the cell and the substrate. Then, protrusion and retraction events correspond to the gain and loss of hexagons at the boundary of the substrate contact area, respectively. The occurrence of these events is determined by a Monte Carlo scheme gradually minimizing an effective energy, H, which is associated with the cell configuration.
The cell is perpetually driven out of equilibrium by active reorganization of its actomyosin network and focal adhesions.

Coarse-grained cellular mechanics
As discussed above, the configuration of a cell at any given time t is associated with a substrate contact area AðtÞ and perimeter PðtÞ. We assume that the membrane and cortex deformations of each cell are constrained by the elastic energy H cont ðtÞ ¼ k A A 2 ðtÞ þ k P P 2 ðtÞ ; (1) where k A and k P are cell-type-specific stiffness parameters, similar to the original implementation of the CPM (Graner and Glazier, 1992). If the cell does not form adhesions to the substrate, then membrane and cortex contractility will round up the cell body, thereby collapsing the substrate contact area into a contact point. Cytoskeletal structures respond to external mechanical stimuli through reaction networks involving different feedback loops. We greatly simplify these complex processes into two prototypic feedback loops, which break detailed balance and drive cell migration, as follows. The polarization field induces membrane protrusions and inhibits retractions. In turn, protrusions increase the polarization field (positive feedback) and therefore the likelihood of further protrusive activity, while retractions decrease the polarization field (negative feedback). In the absence of mechanochemical signals, the polarization field approaches its rest state. (C) Zoom-in to a common boundary shared between the substrate contact areas of three cells (bounded by the red lines), each represented by a contiguous set of occupied grid sites (hexagons).
Top left: The upper right corner of the lower left cell (source cell) initiates a protrusion event against a neighboring element in the cell to its right (target cell), as indicated by the arrow, in an attempt to displace it. The success of each such attempted elementary event depends on the balance between contractile forces, cytoskeletal forces, and cell adhesion. Top right: If the protrusion event is successful, then the levels of regulatory factors are increased (decreased) in integer steps, at all lattice sites inside the source (target) cell that lie within a radius R of the accepted protrusion event (as indicated by the plus and minus signs). Bottom right: During the course of one MCS, different levels of regulatory factors accumulate locally within each cell, with positive levels of regulatory factors (green plus signs) promoting a build-up of cytoskeletal structures, negative levels of regulatory factors (red minus signs) causing degradation of cytoskeletal structures, and neutral levels of regulatory factors (white zero signs) causing relaxation towards a resting state, as indicated in the lower left image. The color code indicates local levels of cytoskeletal structures, .
Gripping the surface through the cell cytoskeleton Detachment of the cell from the substrate is counteracted by focal adhesions, where the cell cytoskeleton is connected to the underlying substrate by integrins. Cellular protrusions are driven by outward pushing forces generated by the assembly and disassembly of cytoskeletal structures (Pollard and Borisy, 2003;Mogilner, 2009). As a first approximation, we subsume all of these complex dynamic processes, like the formation/degradation of focal adhesions and the assembly/disassembly of cytoskeletal structures, into a single time-dependent and spatially resolved internal field for each cell, ðx; tÞ. This polarization field emulates the mass of force-generating cytoskeletal structures in the associated hexagon, at position x, which results in an effective, locally regulated, adhesion energy between cell and substrate. Consequently, the total energy associated with this polarization field is given by The polarization field must vanish at positions that are not occupied by a cell. Therefore, a retraction is associated with an energy penalty due to the loss of a substrate adhesion. Consequently, a protrusion, where one source hexagon 'conquers' a nearby target hexagon, is associated with an energy gain due to an increase of the substrate contact area. Here, we assume that the newly incorporated hexagon has the same polarization field as its conqueror.
There are several biological factors that constrain the local density of actin filaments, myosin and focal adhesions, whose limited availability corresponds to an upper bound on the polarization field. Furthermore, we assume that there is some minimal attachment energy associated with adhesions that prevents the cells from detaching from the substrate, which implies a lower bound on the polarization field. This motivates to introduce cell-type-specific bounds for the polarization field: ðx; tÞ 2 ½ 0 À D=2; 0 þ D=2, where 0 is the average polarization field and D is the maximum cell polarity.
Active self-regulation of the cytoskeleton Assembly and disassembly of cytoskeletal structures are controlled by a myriad of accessory proteins (Lauffenburger and Horwitz, 1996;Ridley et al., 2003). These regulatory proteins form a reaction network involving different feedback mechanisms, which allow cytoskeletal structures to respond to external mechanical stimuli (Marée et al., 2006;Marée et al., 2012). Furthermore, cytoskeletal structures like integrins play a role in the spatiotemporal control of these regulatory proteins (Schwartz and Shattil, 2000). Here, we refrain from formulating a detailed reaction-diffusion model that accounts for the interactions between all of these contributing players. Instead, we assume that the internal chemistry of the cell will generically produce protein patterns, with a typical length scale R, which locally up-or down-regulate cellular cytoskeleton and focal adhesion (dis)assembly. Then, we greatly simplify these complex processes (Lauffenburger and Horwitz, 1996;Schwartz and Shattil, 2000;Ridley et al., 2003) into two prototypic feedback loops (Figure 1B,C): A. The polarization field locally promotes outward motion of the membrane, because it contains a contribution from the local amount of actin filaments. Membrane protrusions facilitate the formation of substrate adhesions and further polymerization of actin filaments, leading to a positive feedback on the polarization field within a range R. B. The polarization field also locally inhibits inward motion of the membrane, by emulating the local adhesion strength of the cell to the substrate. If a membrane retraction is successful, then the loss of substrate adhesions locally further increases cell contractility, leading to a negative feedback on the polarization field within a range R.
In the absence of regulatory signals, we assume that the polarization field decays to a fixed value, ! 0 , which corresponds to a resting state of the cell cytoskeleton and focal adhesions. For the sake of keeping our model as simple as possible, we assume that all protein patterns have the same range R, and that the regulation of the cell cytoskeleton and focal adhesions follows a single timescale that corresponds to an update rate . Because at heart, our model is only based on generic feedback loops with a certain signaling range R, we would argue that any model with similar feedback should, in general, lead to similar cell behavior. Indeed, mutually repressing feedback loops (Marée et al., 2006) and mutually activating feedback loops (Shao et al., 2010;Ziebert et al., 2012;Albert and Schwarz, 2016) are crucial recurring motifs among multiple cell migration studies. Notably, these theoretical approaches all recover comparable cell behavior even when the model setup seems quite different at first glance: 1. Cell migration couples mechanochemically to a scalar field (Shao et al., 2010), if stresses in the cell are isotropic; this is analogous to the present study. 2. Cell migration couples mechanochemically to a vector field (Marée et al., 2006;Ziebert et al., 2012), if stresses in the cell are anisotropic. 3. Cell migration couples to a single polarity vector (Albert and Schwarz, 2016), if propulsive forces are distributed homogeneously throughout the cell. However, this simplification of the former two cases cannot account for the formation of multiple competing lamellopodia/ pseudopods.
These different modeling approaches (of varying complexity) surprisingly yield a universal phenomenology. The puzzling similarity between these models suggests generic common features that determine cell shape and motility: mechanical constraints like cell elasticity and mechanochemical feedback mechanisms that break detailed balance, maintain cell polarity and drive cell motion.

Intercellular adhesion and friction
In addition to internal remodeling of the cytoskeleton, adhesion of cells to neighboring cells and to the substrate plays a key role in explaining migratory phenotypes (Mogilner, 2009;Friedl and Gilmour, 2009). From a mechanical point of view, the implications of cell adhesion are two-fold: 1. Cell adhesion supports growth of cell-cell and cell-matrix contacts and may thus be described in terms of effective surface energies. In our computational model, cell-matrix contacts are readily accounted for by the polarization field, . In addition, we associate the formation of cell-cell adhesions with an energy benefit B, which we call cell-cell adhesion parameter. 2. Once formed, adhesive bonds anchor the cell to the substrate and to neighboring cells. During cell migration, these anchoring points must continuously be broken up and reassembled (Webb et al., 2002;Gumbiner, 2005) and, hence, provide a constant source of energy dissipation. Therefore, we assume that the cost for rupturing an existing cell-cell adhesion, B þ DB>B, exceeds the gain from forming a new cell-cell adhesion. Then, the dissipative nature of cell-cell adhesions is accounted for by the cell-cell friction parameter DB. Similarly, cell-matrix contacts can also provide a source of dissipation, which is further discussed in Appendix 2.

Environmental cues
The polarization field, , readily includes contributions from cell-substrate adhesions, which are locally up-or down-regulated by the cell. These cell-substrate adhesions require the abundance of surface ligands, which serve as substrate tethers that the cell can attach to, and which are not necessarily distributed homogeneously. By substrate micropatterning, one can arrange areas where the cell is likely to adhere to the surface, and no-go-areas, where the cell adheres less (or cannot adhere at all). To replicate such environmental cues, we introduce a second scalar field 'ðxÞ, whose value is taken to reflect the relative availability of substrate sites at which focal adhesions between cell and substrate can be formed. Here, we have chosen to model micropatterns as impenetrable walls; we locally add a large energy penalty, ' ( 0, to the polarization field ( ! þ '), that a cell has to pay for trespassing onto a no-go-area. However, it is equally valid to treat ' as a multiplicative constant modulating the polarization field ( ! ' ), where ' ¼ 0 models a local inability of the cell to attach to the substrate. Analogously to cell-cell contacts, we account for the dissipative nature of cell-substrate adhesions by associating the breaking of such contacts with an additional energy cost D.

Tissue growth by cell division
In the description so far, the cells are arrested in the cell cycle (mitostatic). To investigate the effect of cell proliferation on tissue dynamics, we introduce a simplified three-state model of cell division. Cells start off in a quiescent state, in which their properties remain constant over time. The cell sizes fluctuate around an average value determined by the cell properties and the local tissue pressure. Cell growth typically arrests at large cell densities, in a phenomenon coined contact inhibition of proliferation (Stoker and Rubin, 1967;Puliafito et al., 2012;Pavel et al., 2018). Since large cell densities correspond to a small spread area for each individual cell, this implies that cell growth is arrested below a critical threshold size (A T ). Upon exceeding this threshold size due to size fluctuations, cells leave the quiescent state and enter a growth state. The duration of the quiescent state is thus a random variable, whose average value depends on the tissue pressure, and lower pressure (due to a lower cell density) leads to a shorter quiescent state. During the subsequent deterministic growth state of duration T g , cells double all of their cellular material and thus double in size. We model this growth as a gradual decrease in the effective cell contractility (k A and k P ). As there is no a priori reason to assume that a cell's migratory behavior should depend on its size, we constrain the parameters accordingly; this is described in detail in Appendix 2. After having grown for a duration T g , cells switch to a deterministic division state of duration T d . During division, cells strongly contract, which leads to mitotic rounding and a drastic decrease of their contact area with the substrate Lock et al., 2018). In principle, a decrease of cell contact area could also lead to perturbations of the stress field in the monolayer. Here, however, we neglect the decrease of the cell spreading area, as the division phase is short compared to the growth phase. We expect that a drastic increase of cell contractility also leads to a loss of polarity in the cell's migratory machinery. Therefore, each cell reduces its polarizability to zero (D ! 0) in order to utilize its cytoskeleton for the separation of the cellular material, leading to mitotic rounding. At the end of the division state, each dividing cell splits into two identical daughter cells, whose properties and parameters are identical to the mother cell's initial values in the quiescent state. Finally, the daughter cells re-initialize migration from an unpolarized state. For a detailed and more technical description we refer the interested reader to Appendix 1.

Persistent migration of single cells
The macroscopic properties of cell clusters and tissues emerge from an interplay between many individual cells. Then, what determines the mechanical and migratory features of these individual cells? In our computational model, we have studied this question by screening its multidimensional parameter space. For such a brute force approach to be numerically feasible, one must first distinguish relevant parameters (these determine the resulting dynamics) from irrelevant parameters. Specifically, in our extended cellular Potts model, there are reference parameters whose sole purpose is to control the spatial and temporal discretization of the numerical model: 1. The cytoskeletal update rate endows the cellular Potts model with a reference timescale and determines the temporal discretization. In this study, we have set ¼ 0:1. 2. The average polarization field 0 encodes the energy gain for creating new cell-substrate adhesions, while the area stiffness k A represents the energy cost for increasing the substrate contact area. Then, the number of hexagons occupied by the cell is proportional to the ratio 0 =k A . If we use a desired cell area as reference value, then the ratio 0 =k A controls the spatial discretization of the cell. To study the migration of single cells and small cell cohorts, we have set the average polarization field to 0 ¼ 225 and the area stiffness to k A ¼ 0:18. 3. In cellular Potts models, which are Monte-Carlo simulations, the reference energy of fluctuations is determined by an effective temperature. In this study, we have set k B T 1.
Furthermore, we used a large computational grid with 9 Á 10 4 sites and periodic boundary conditions to study the migration of single cells. This leaves three parameters that control cell motility in the absence of cell-substrate dissipation: cell polarizability D, cell contractility k P and signalling radius R. However, it is not clear yet whether all of these are independent relevant parameters. In fact, in the following sections it will become clear that cell polarizability and contractility are degenerate parameters (in the sense that the phenomenology only depends strongly on their ratio, which is the corresponding relevant parameter).

Cell persistence increases with polarizability
First, we investigated the impact of varying levels of cell perimeter stiffness k P and maximum cell polarity D on the cell's migratory patterns ( Figure 2-video 1), at a fixed signaling radius R ¼ 5. To assess the statistics of the cell trajectories, we recorded the cell's orientationvðtÞ vðtÞ=kvðtÞk (v: cell velocity) and (geometrical) center of mass position RðtÞ during a total simulation time of T sim ¼ 10 4 Monte-Carlo steps (MCS). For each set of parameters, we performed 100 statistically independent simulations, from which we computed the mean squared displacement, MSDðtÞ h½Rðt þ tÞ À RðtÞ 2 i, and the normalized velocity auto-correlation function, CðtÞ hvðt þ tÞ ÁvðtÞi. Here, h. . .i denotes an average with respect to simulation time t as well as over all 100 independent simulations.
These computer simulations show that the statistics of the migratory patterns is well described by a persistent random walk model (Stokes et al., 1991;Wu et al., 2014) with its two hallmarks: a mean square displacement that exhibits a crossover from ballistic to diffusive motion (Figure 2A), and on sufficiently long time scales an exponential decay of the velocity autocorrelation function CðtÞ / e Àt=tp (inset of Figure 2A). We determined the persistence time of directed migration, t p , by fitting the mean squared displacement with a persistent random walk model. In addition, we also measured cell speed, v, and cell aspect ratio, l þ =l À , to further characterize cell motility and shape. for single-cell movements at different maximum cell polarity D (stiffness parameters k P ¼ 0:060, k A ¼ 0:18; average polarization field 0 ¼ 225; signaling radius R ¼ 5; cell-substrate dissipation D ¼ 0; cellsubstrate adhesion penalty ' ¼ 0; cytoskeletal update rate ¼ 0:1; 100 independent simulations for each set of parameters). Single cells perform a persistent random walk, i.e. they move ballistically (MSD / t 2 ) for t ( t p , and diffusively (MSD / t) for t ) t p . Inset: Normalized velocity auto-correlation function for the same parameters as in the main figure. (B) Persistence time of directed cell migration plotted as a function of maximum cell polarity D, and perimeter stiffness k P (area stiffness k A ¼ 0:18; average polarization field 0 ¼ 225; signaling radius R ¼ 5; cellsubstrate dissipation D ¼ 0; cell-substrate adhesion penalty ' ¼ 0; cytoskeletal update rate ¼ 0:1; 100 independent simulations for each set of parameters). The persistence time of the random walk increases with increasing cytoskeletal polarity and decreasing perimeter elasticity. (C) Cytoskeletal polarity also controls cell shapes, with crescent cell shapes (long persistence times) being observed at large cytoskeletal polartities, and round cell shapes (short persistence times) at small cytoskeletal polarities. Color code: cell polarization; cf. color bar in Figure 1C.  Surprisingly, for each of these variables we found a master curve that only depends on the ratio between cell polarizability and cell contractility, D=k P ( Figure 2B,D,E). This data collapse suggests D=k P as a relevant parameter (while cell polarizability and contractility are degenerate parameters), which we will henceforth refer to as specific polarizability. The cells' persistence times of directed migration, speeds and aspect ratios all show a characteristic dependence on the specific cell polarizability. There is a threshold value for the specific polarizability, D=k P » 500, below which cells remain immobile ( Figure 2B,D,E; grey regions). Above this threshold, the persistence time of directed migration, speed and aspect ratio increase markedly with the specific polarizability ( Figure 2B,D,E). In our model, the area and perimeter stiffnesses refer to global and homogeneous cell contractility, while the cell polarization field drives cell migration. As discussed in 'Gripping the surface through the cell cytoskeleton', the cell polarization field does not explicitly distinguish between a local extensibility (e.g. due to actin polymerization), a local contractility (due to myosin-induced contraction) of the cytoskeleton or spatially regulated cell-substrate adhesions. For example, if cell migration is driven by actin polymerization, then blebbistatin treatment will decrease the global cell contractility, which we predict to lead to more elongated cells that move faster and exhibit extended episodes of ballistic motion. Indeed, an increase of cell migration speed after blebbistatin treatment was observed for mouse hepatic stellate cells (Liu et al., 2010). Alternatively, cell migration could also be driven by myosin contractility, for example by pulling the cell forward or by locally detaching adhesions. Then, polarizability and contractility concomitantly depend on the ability of the cell to exert forces, which can be inhibited by blebbistatin treatment. If polarizability, D, and contractility, k P , are equally reduced by a blebbistatin-dependent prefactor, then the specific polarizability, D=k P , and the resulting cell phenomenology should remain unchanged. Indeed, blebbistatin treatment of keratocytes and keratocyte fragments was reported not to affect cell shape and speed to any significant degree (Wilson et al., 2010;Ofer et al., 2011). Therefore, blebbistatin treatment can either increase or decrease cell motility, depending on the cell type and possibly on the specific mechanism that drives cell migration.
Interestingly, because of this universal dependence of all the mentioned quantities on the specific polarizability, our simulations also show that there is a strong correlation between cell shape (aspect ratio) and cell motility (speed and persistence time of directed migration); see Figure 2F. While highly persistent trajectories are observed for cells with 'crescent' shapes, more erratic cell motion is typically found for cells with more rounded outlines ( Figure 2C). In other words, our computational model predicts that cells which are able to polarize their cytoskeletal structures more strongly will adopt crescent shapes and show a higher degree of persistent cell motion. It would be interesting to further test these predictions by using phenotypic variations in cell shapes like those reported in experiments with keratocytes (Keren et al., 2008); there, the authors also found a correlation between cell shape and speed.

Feedback range determines whether individual cells move persistently or rotate
Moreover, we investigated the influence of different signaling radii R (typical range in which signalling molecules diffuse and mediate feedback mechanisms during a single Monte-Carlo step) on the persistence of single-cell trajectories. Since R is the relevant parameter that controls the spatial organization of lamellipodium formation, its value should strongly affect the statistics of a cell's trajectory ( Figure 3A). Indeed, at small values of R, we observe that the spatial coherence of cytoskeletal rearrangements is low, which frequently results in the disruption of ballistic motion due to the formation of independent lamellipodia in spatially separate sectors of the cell boundary ( Figure 3C, lower snapshot). In contrast, at larger values of R, we find that spatial coherence is restored, and the formation of one extended lamellipodium across the cell's leading edge maintains a distinct front-rear axis of cell polarity ( Figure 3C, upper snapshot). However, when the signaling radius is too large compared to the cell size, we find an inhibition of ballistic motion and rounding of the cells as signals originating from one cell edge begin to reach the opposing edge. This effect may also occur when cells in tissue become smaller due to an increase of cell density through proliferation or compression; in other words, this means that the cells become smaller than the typical length scale of the chemical patterns that control cell migration. Then, one would not expect these chemical patterns to form (Hubatsch et al., 2019). Therefore, depending on the cell polarizability (D), there is an optimal signaling radius that shows both maximal cell elongation and maximal cell persistence ( Figure 3A,B).
Cells with low polarizability need a large signaling radius to feed the positive feedback mechanism and to form a single large cell front. In contrast, highly polarizable cells can already sustain the positive feedback mechanism with a short signaling radius and easily form at least one (or even multiple competing) short cell front(s). With increasing signaling radius, these cell fronts become increasingly correlated and finally merge. Surprisingly, at small signaling radii, we observed that highly polarizable cells slow down with increasing signaling radius ( Figure 3D; yellow squares and black circles), in contrast to the behavior of cells with low polarizability. Furthermore, at large signaling radii, highly polarizable cells speed up, although their persistence time of directed migration has dropped to small values (cf. Figure (bottom), it has to move a different projected contour length. If each protrusion takes roughly the same amount of time, then migration along the long axis (top; cell has to move a smaller projected contour length) allows for greater cell speeds than migration along the short axis (bottom; cell has to move a larger projected contour length).
The online version of this article includes the following video for figure 3: Figure 3-video 1. Single cell motility and shape for different signaling radii (D ¼ 60, k P ¼ 0:060). https://elifesciences.org/articles/46842#fig3video1 black circles). To find an intuitive explanation for these observations, we inspected time-lapse videos of a cell at high polarizability (D=k P ¼ 1000; cf. Figure 3-video 1, top row), which show a qualitative shift in cell behavior: . For small signaling radii, R ¼ 2, short polarization fronts 'pull' the cell behind them, allowing for transient polarization and quick but erratic movement along the long axis of the cell. . For intermediate signaling radii, R ¼ 6, broad and correlated polarization fronts emerge, and both the cell polarization and movement always orient themselves along the short axis of the cell. . For large signaling radii R ¼ 15, we observed circular motion of the cell; because of the large signaling radius, signals originating from the trailing edge affect the leading edge of the cell and vice versa. Due to this circular motion, the cell exhibits a non-zero speed and a vanishing persistence time of directed migration.
Therefore, we find that the cell can transiently polarize and migrate along its long axis for small signaling radii and for high polarizability. Furthermore, in a broad parameter regime, we find keratocyte-like motion and polarization along the short axis of the cell. Note that we do not consider the formation of stress fibers, which could lead to cell migration along the long axis in a broad parameter regime (Kassianidou et al., 2019). Such stress fibers could be modeled via a nematic field that represents the anisotropic part of the intracellular stress. Our counter-intuitive observation that cell migration along the long axis is faster than cell migration along the short axis can be explained as follows: If the cell migrates along its short axis, then it has to move a greater projected contour length than if it migrates along its long axis ( Figure 3F). Considering that each protrusion takes roughly the same amount of time, migration along the long axis allows for greater cell speeds than migration along the short axis, because the cell has to spend less time to move a smaller projected contour length ( Figure 3F).
To further characterize the single cell rotations that occur at large signaling radii, we determined the average curvature of the trajectories hci ¼ hkq sv ðsÞ ÁvðsÞki, where s is the contour length along the corresponding trajectory. Here, we averaged the tangent vectorvðsÞ over 10 Monte-Carlo steps to integrate out fluctuations that occur on short timescales (the internal dynamics of the cell has an intrinsic time scale of 10 Monte-Carlo steps due to our choice of the cytoskeletal update rate, ¼ 0:1). We find that the curvature of the trajectories has a pronounced minimum at large signaling radii (where the persistence time of directed migration vanishes), which indicates a transition from straight to circular trajectories ( Figure 3E). Such a transition from persistent migration to single cell rotations was previously observed in experiments (Lou et al., 2015;Raynaud et al., 2016) and in theory (Reeves et al., 2018;Allen et al., 2018).

Cell clusters on circular micropatterns
To assess the transition to collective cell motion, we next studied the dynamics of small cell groups confined to circular micropatterns (Huang et al., 2005;Doxzen et al., 2013;Deforet et al., 2014;Segerer et al., 2015). We implemented these structures in silico by setting 'ðxÞ ¼ 0 inside a radius r 0 and 'ðxÞ ! À¥ outside. During each simulation run, the number of cells was also kept constant by deactivating cell division. We previously employed this setup to compare our numerical results with actual experimental measurements, and found very good agreement (Segerer et al., 2015). Here, we generalize these studies and present a detailed analysis of the statistical properties of the collective dynamics of cell groups in terms of the key parameters of the computational model.
When adhesive groups of two or more motile cells are confined on a circular island, they arrange themselves in a state of spontaneous collective migration, which manifests itself in the form of coordinated and highly persistent cell rotations about the island's midpoint x 0 (Huang et al., 2005;Doxzen et al., 2013;Deforet et al., 2014;Segerer et al., 2015). The statistics of these states of rotational motion provide insight into the influence of cellular properties on the group's ability to coordinate cell movements. To quantify collective rotations, we recorded the average signed angular velocity of the cell cluster !ðtÞ ¼ê z Á hṽðtÞ ÂRðtÞ=kRðtÞk 2 i C . Here,ê z is the out-of-plane unit vector, . . . h i C denotes an average with respect to the cell population, andṽðtÞ ¼ vðtÞ À hvðtÞi C as well as R ¼ RðtÞ À hRðtÞi C measure the velocity and position of each cell relative to the cell cluster (we have omitted the indices that identify individual cells for the sake of convenience and clarity). The resulting random variables for the magnitude of the angular velocity of the cell assembly, j!ðtÞj, and the average cell perimeter PðtÞ hP a ðtÞi C were then used to characterize the statistics of collective cell rotation. For each specific choice of simulation parameters, we monitored j!ðtÞj and PðtÞ for a set of 100 statistically independent systems, each of which was observed over T sim ¼ 10 4 MCS. From these data, we then computed the mean overall rotation speed hj!ji, its standard deviation s ! , and the standard deviation of the cell perimeter, s P . Figure 4 illustrates the characteristic properties of collective cell rotations in systems containing jCj ¼ 4 cells endowed with varying maximum cell polarity D and varying cell contractility k P . Analogously to our observations for single cells, the statistical measures shown in Figure 4A do not separately depend on cell contractility and maximum cell polarity, but depend only on the specific polarizability D=k P . Overall, we find that upon increasing the specific polarizability there is a marked transition from a quiescent state to a state where the cells are collectively moving. Below a threshold value for the specific polarizability (D=k P » 450 in Figure 4A), the rotation speed hj!ji (purple curves in Figure 4A) vanishes and the cells are immobile. In this regime, which we term the stagnation phase, or S-phase, cytoskeletal forces are too weak to initiate coherent cell rotation, and the system's dynamics is dominated by relatively strong contractile forces, which tend to arrest the system in a 'low energy' configuration. Beyond this threshold, we identify three distinct phases of collective cell rotation. In the R 1 -phase, we find a steep increase in the average rotation speed and a local maximum in the fluctuations of both cell shape and rotation speed; cf. green (s P ) and blue (s ! ) curves in Figure 4A. Now, cytoskeletal forces are sufficiently large to establish actual membrane  Figure 4. Phases of collective motion. (4-cell systems; confinement radius r 0 ¼ 30:6; area stiffness k A ¼ 0:18; average polarization field 0 ¼ 225; signaling radius R ¼ 5; cytoskeletal update rate ¼ 0:1; cell-cell adhesion    protrusions against the contractile forces, and cells begin to rotate ( Figure 4B,C). However, the contractile forces still dominate, such that cellular interfaces tend to straighten out and lamellipodium formation is sustained only over finite lifetimes. Thus, due to the dominance of contractile forces, the systems frequently experience transient episodes of stagnation and repeatedly change their direction of rotation (cf. blue trajectory in Figure 4B).
At intermediate values of specific polarizability (R 2 -phase), the cellular systems reach a regime of enduring rotational motion, where hj!ji varies linearly with the local specific polarizability, and where s P and s ! exhibit a rather broad minimum ( Figure 4A). In this regime, a range of 'optimal ratios' of cytoskeletal to contractile forces sustains stable cell shapes, and sets the stage for the formation of extended lamellipodia and the establishment of permanent front-rear polarizations of cells. As a result, the cells' persistence times of directed migration become very large, rendering cellular rotations strictly unidirectional within the observed time window ( Figure 4B). Finally, at large values of the specific polarizability (R 3 -phase), the system's dynamics is dominated by cytoskeletal forces and the rotational speed hj!ji saturates at some maximal value. Due to the relatively small contractile forces, cell shapes tend to become unstable, as reflected in the growing variance of the cell perimeter s P (green curve in Figure 4A). These instabilities in cell shape frequently lead to a loss of persistence in the rotational motion of the cells (growing s ! ; blue curve in Figure 4A).

Tissue-level dynamics
As an application of our computational model at the tissue level, we considered a setup in which an epithelial cell sheet expands into free space. As in recent experimental studies (Serra-Picamal et al., 2012;Sepúlveda et al., 2013;Trepat et al., 2009;Poujade et al., 2007), we confined cells laterally between two fixed boundaries, within which they proliferated until they reached confluence; in the ydirection we imposed periodic boundary conditions. Then we removed the boundaries and studied how the cell sheet expands. In order to quantify tissue expansion, we monitored cell density and velocity, as well as the mechanical stresses driving the expansion process. Figure 5 shows our results for two representative parameter regimes that highlight the difference between a dynamics dominated by cell motility in the absence of cell proliferation, and a contrasting regime where cells with low motility grow and divide depending on the local cell density. To simulate large numbers of cells, we decreased the amount of hexagons that are typically occupied by each cell (the simulation cost scales linearly with the summed area of all cells) by setting the average polarization field to 0 ¼ 35. For each set of parameters, we performed and averaged 100 independent simulations.
We first investigated how a densely packed pre-grown tissue of mitostatic cells with high polarizability (large D) expands into cell-free space upon removal of the confining boundaries at the tissue's lateral edges ( Figure 5A). As the cells migrate into the cell-free space, we observe a strongly (spatially) heterogeneous decrease in the initially high and uniform cell density and mechanical pressure in the expanding monolayer ( Figure 5B,C). This is quite distinct from the behavior of a homogeneous and ideally elastic thin sheet, which would simply show a homogeneous relaxation in density as it relaxes towards its rest state. Moreover, cell polarization and the ensuing active cell migration lead to inhomogeneously distributed traction stresses in the monolayer. After initial expansion of the monolayer, facilitated by high mechanical pressure, the cells at the monolayer edge begin to polarize outwards, which enhances outward front migration. These actively propagating cells exert traction on the trailing cells, and thereby yield a trailing region with negative stress ( Figure 5C). Taken together, this gives rise to a characteristic X-shaped pattern in the kymograph of the total mechanical stresses hs xx i y ( Figure 5C). This profile closely resembles the first period of mechanical waves observed experimentally (Serra-Picamal et al., 2012). It illustrates how stress is transferred towards the center of the monolayer when cells are highly motile and collectively contribute to tissue expansion. At the end of the simulated time window, the cell density exhibited a minimum in the center of the sheet ( Figure 5B). This is due to stretching of the central group of cells caused by the equally strong traction forces exerted by their migrating neighbors on both sides. Finally, the simulations also show that outward cell velocities increase approximately linearly with the distance from the center, confirming that in this configuration the entire cell sheet contributes to the monolayer expansion ( Figure 5D).
To explore the possible range of tissue dynamics and expansion, we also investigated a qualitatively different parameter regime where cells are less densely packed and can also polarize less due to a narrower range of polarizability ( Figure 5E). Here, the expansion of the monolayer is mainly driven by cell division, and cells keep dividing until they reach a homeostatic cell density ( Figure 5F). Even though cells should typically exceed the threshold size and hence enter the growth phase at different times, we observe that the cell sheet exhibits periodic 'bursts' of growth ( Figure 5F) coinciding with the total duration of a complete cell cycle (200 MCS) and alternating with cell migration ( Figure 5H). These periodic 'bursts' can be explained as follows. Initially, the slightly compressed monolayer expands to relieve mechanical pressure. Due to this initial motion, the cells at the monolayer edge begin to polarize outwards. As in the previous case, where cell proliferation is absent ( Figure 5A-D), the polarized cells enhance outward front migration and stretch the cells in the bulk of the cell sheet. For the same reasons as before, we observe a typical X-shaped stress pattern in the kymograph ( Figure 5G), albeit less pronounced due to the lower polarizability of the cells (cf. Figure 5C). Because a broad region in the monolayer bulk is stretched by the actively migrating cell Expansion of a confluent epithelial cell sheet after removal of boundaries positioned at  fronts, these cells exceed the threshold size and begin growing approximately in phase. Once the mechanical pressure of the cell sheet is relieved, it will stop expanding ( Figure 5H). However, cell growth and division once more lead to an increase in mechanical pressure (and cell density) in the monolayer ( Figure 5F,G). This cycle of migration-dominated monolayer expansion and cell-densitydependent cell growth and division results in a periodic recurrence of the X-shaped stress pattern ( Figure 5G), closely resembling the pattern observed in experiments (Serra-Picamal et al., 2012). On a sidenote, the synchronization of the cell division and cell migration phases by the deterministic portion of the cell cycle can be counteracted by introducing additional stochastic terms in the transition between the different phases of the cell cycle (cf. 'Cell proliferation and mitosis' in Appendix 1). Note that the inhomogeneously distributed traction stresses in the monolayer, and its wave-like behavior, ultimately emerge from cell polarization and the ensuing active cell migration. Therefore, The online version of this article includes the following video(s) for figure 6: Figure 6-video 1. Weak monolayer roughening (fingering) in motility-dominated tissue with quick proliferation. https://elifesciences.org/articles/46842#fig6video1 Figure 6-video 2. Strong monolayer roughening in motility-dominated tissue with slow proliferation. https://elifesciences.org/articles/46842#fig6video2 these traction patterns would look much less prominent if one were to inhibit cell motility (compare Figure 5C with Figure 5G). Finally, we investigated which parameters control the roughness of the tissue fronts. We found that increasing cell motility, or increasing cell-cell dissipation leads to rougher front morphologies ( Figure 5-figure supplement 1 and 'Velocity and roughness of spreading tissue' in Appendix 2). Therefore, we hypothesized that one could observe fingering of cell monolayers by adjusting the parameters accordingly: . Increase of cell motility by decreasing the membrane stiffness and at the same time increasing polarizability and signaling radius of the cells. . Increase of cell-cell dissipation and slight decrease of cell-cell adhesion. . Slower and less homogeneously distributed cell division by increasing the cell threshold size.
Indeed, we then observe a drastic roughening of the cell fronts and small cohorts of cells that coherently move into cell-free space ( Figure 6). This roughening is more pronounced if we further increase the threshold size that a cell has to exceed to initiate growth (cf. Figure 6A,E). Analogously to our previous discussion, we observe that an increasing mechanical pressure in the monolayer due to the division of cells initiates outward cell migration ( Figure 6B,F). Then, cells in the tissue begin to polarize outwards and coordinate their motion with their neighboring cells, leading to small coordinated cell cohorts. As before, we also find distinct traction force patterns, as recurring waves of high stress travel backwards relative to the leading edges ( Figure 6C,G), and distinct recurring velocity patterns ( Figure 6D,H).

Discussion
In this work, we have proposed a generalization of the cellular Potts model (Graner and Glazier, 1992). The model implements a coarse-grained routine that captures the salient features of cytoskeletal remodeling processes on subcellular scales, while being computationally tractable enough to allow for the simulation of entire tissues containing up to Oð10 4 Þ cells. We have used the model to study the transition from single-cell to cohort cell migration in terms of the interplay between the pertinent cellular functions. Specifically, we have demonstrated that our model consistently reproduces the dynamics and morphology of motile cells down to the level of solitary cells. Our studies also Table 1. Source and parameter files used for each figure.
All source and parameter files are found in Source data 1.     reveal that cytoskeletal forces (relative to cell contractility), as well as the spatial organization of the cells' lamellipodia, significantly affect the statistics of cellular trajectories, both in the context of single-cell motion and in cohesive cell groups restricted to circular micropatterns. On larger scales, our simulation results suggest that the dynamics of expanding tissues strongly depends on the specific properties of the constituent cells. If monolayer expansion is driven by active cell migration throughout the tissue, then the cell sheet exhibits typical traction-force patterns and an X-shape in the corresponding kymograph. Additionally, a cell-density-dependent cell growth leads to a periodic recurrence of these traction-force patterns in a cycle of migration-dominated expansion and 'burst'like cell proliferation. Taken together, our results further highlight the intricacies of collective cell migration, which involves a multitude of intra-and inter-cellular signaling mechanisms operating at different scales in length and time. Establishing a comprehensive picture that incorporates and elucidates the mechanistic basis of these phenomena remains a pressing and challenging task. The multiscale modeling approach proposed here provides a direct link between subcellular processes and macroscopic dynamic observables, and might thus offer a viable route towards this goal.

Materials and methods
The computational model is described in section 'Computational model'. The numerical implementation of the model is discussed in detail in Appendix 1. The parameter files and source files associated with the figures are given in Table 1. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Computational model
In this section, we describe in detail the implementation of our computational model, which has been outlined briefly in the main text. While the biological rationale behind our modeling approach has been discussed in the main text, our focus here is on the technical aspects and the details of the numerical implementation. To facilitate subsequent discussions on implementation details, we start by introducing some model-specific terminology which will be used throughout this section to illustrate the mechanics of our model.

Computational grid
The basic data structure that underlies our computational model is referred to as the grid; see Appendix 1-figure 1. The grid itself is implemented as a regular, space-filling lattice with lattice vectors x i f g i¼1;...;N . Each lattice vector x i is understood to represent its associated Voronoi cell which will be referred to as grid site. To be specific, we consider triangular tilings x i f g i¼1;...;N , such that each grid site is a hexagon, which is surrounded by 6 nearest-neighbor sites that define the neighborhood N k of x k : Overall, the grid represents our general notion of (discretized) space, and each grid site holds information specific to cells as well as to environmental factors. In what follows, distances on this spatial grid will be measured in units of the distance between the midpoints of neighboring lattice sites, i.e.
This then implies for the side length ' and the two-dimensional volume (area) a of each hexagonal grid site:

Representation of biological cells
In the spirit of the cellular Potts model (Graner and Glazier, 1992;Glazier and Graner, 1993), each cell is represented by a simply connected set of lattice sites where the indicator function cðx k Þ gives the index of the cell occupying x k . Here and in the following, we use latin indices to reference lattice sites, and greek indices to reference cells.
The set D ðaÞ used to represent the spatial extension of cell a, will be referred to as the domain of cell a. In our model, each grid site x k can be occupied by at most one cell (i.e. we do not allow for overlapping cell domains). The absence of cells at x k is numerically implemented by negative values of the indicator function, cðx k Þ<0. Following this terminology, the area and the perimeter of cell a are given by:

Model dynamics
Protrusion and retraction of cells where the set of membrane grid sites B ðaÞ is defined as Protrusion and retraction events are the numerical analogs of cell protrusions and cell retractions.
In implementing the reassignment rules, Equation S5 and Equation S7, we have to take into account that cellular domains must not overlap. For solitary cells moving in free space this does not imply any restrictions, and Equation S5 and Equation S7 apply directly. In the bulk of a confluent monolayer of adhesive cells, however, any protrusion of source cell a into the domain of cell b (referred to as target cell) must be accompanied by a corresponding retraction event D ðbÞ old ! D ðbÞ new ¼ D ðbÞ old n fx t g, where x t denotes the target grid site annexed by cell a. We emphasize, however, that the reverse is not generally true. If source cell a retracts, i.e. loses one of its boundary grid sites x s 2 B ðaÞ , the lost grid site x s faces either one of two conceivable fates: If, on the one hand, cohesion among cells is sufficiently strong (cf. section 'Rupture of cell contacts' for a definition of the notion 'sufficiently strong'), then the retraction of cell a exerts a pulling force on one of its neighboring cells b (the target cell) and forces the target cell to fill the emerging void at x s , i.e. D The implementation of a MCS, i.e. the sequence of elementary events, is based on the following general considerations: i. Choice of source and target grid sites. Each elementary event T originates from a membrane grid site x s 2 B ðaÞ of some cell a, referred to as source cell. This membrane grid site will be referred to as source grid site. In addition, each elementary event involves a second grid site which lies in the neighborhood of the source grid site x s and which is not currently occupied by cell a: x t 2 N s n D ðaÞ . In what follows, this additional grid site x t will be referred to as target grid site. This grid site may either be an empty substrate site or a membrane site of another cell b, in which case the respective cell will be referred to as target cell. While the source grid site determines the location of the attempted elementary event, the target grid site determines the direction along which the elementary event is bound to proceed. ii. Monte-Carlo method to generate the system's dynamics. As mentioned above, the actual dynamics of cells in our computational model is driven by a succession of elementary events, whose cumulative effects over time allow cells to change shapes and to move relative to the substrate as well as relative to each other. Following a standard Monte-Carlo procedure, the probability of occurrence of elementary events T is determined by a goal function pðT Þ [cf. point (iii) below]. However, since elementary events come in two basic flavors, protrusions T pro and retractions T ret , their actual occurrence is controlled by a twostep process, once source and target grid sites have been determined: In a first step, two alternative scenarios are proposed where either the source cell protrudes toward x t , or retracts from x s . Then, a decision is made with equal probabilities as to whether one attempts T pro or T ret . In a second step, the goal function p is used to compute the occurrence probability of the attempted event T . Finally, this elementary event T is being accepted with probability pðT Þ. iii. Choice of the goal function pðT Þ. As has been detailed above, we use a goal function pðT Þ to control the occurrence and acceptance of elementary events T . Following the standard cellular Potts model (Graner and Glazier, 1992;Glazier and Graner, 1993), this goal function takes into account the effects of cell contractility and cell-cell adhesion, using, however, a slightly different implementation; cf. sections 'Cell contractility' and 'Cell adhesion'. In addition, we generalized the definition of the goal function pðT Þ to explicitly take into account a simplified model of cytoskeletal structures and the ensuing polarization of cells.
The actual definition of the goal function will be developed in section 'Implementation of cellular traits', where, moreover, details concerning the implementation of the cell polarization model will be discussed.
The implementation of a single MCS loop is then given by the following simulation scheme: 1. Determine the current number of trials per Monte-Carlo step (MCS), K ¼ P a P a =', and set the trial counter n ¼ 0.
2. With equal probability, choose a cell membrane segment (cf. solid black line in Appendix 1- figure 1) from a random cell a of the cell population. Because the cell membrane represents the border between lattice sites occupied by cell a and unoccupied by cell a, the chosen membrane segment automatically defines the source grid site x s 2 B ðaÞ and the corresponding target grid site x t 2 N ðaÞ \ N s . 3. With equal probability, choose whether to attempt a protrusion event (T pro ) or a retraction event (T ret ). 4. Compute the prospective acceptance probability pðT pro=ret Þ corresponding to the attempted event, and decide whether to accept the attempted event on the basis of this probability. 5. If the attempted elementary event has been accepted, then update the cellular domains of source cell a and opponent cell b; for details see section 'Cell domain update routine'. 6. If n<K, set n ! n þ 1 and then repeat steps 2 through 5.

Implementation of cellular traits
In this section, we discuss the various contributions of cellular traits to the overall acceptance probability pðT Þ of an elementary event T . Specifically, our model takes into account cell contractility, the assembly and disassembly of cytoskeletal structures, cell-cell adhesion, and focal adhesions. We will assume that each of these cellular properties contributes independently to the acceptance probability p, such that Anticipating our discussions in section 'Cytoskeletal structures and focal adhesion', the effects due to focal adhesions have been combined with the effects due to assembly and disassembly of cytoskeletal structures in p cyto ðT Þ. In the following sections, we give detailed discussions for each of these contributions, separately.

Cell contractility
In biological cells, membrane fluctuations are constrained by elastic forces and contractile cytoskeletal structures, which play a vital role in cell migration (Alberts et al., 2015;Raucher and Sheetz, 2000;Friedl, 2004). In our computational approach, we take cell contractility into account by assigning a contractile 'energy' (S10) with positive coupling constants k ðaÞ P and k ðaÞ A characterizing the contractility of cell a; for empty substrate sites (a<0) we set k ðaÞ P ¼ k ðaÞ A ¼ 0. According to Equation S10, the cell's 'contractile energy' increases with increasing cell perimeter and increasing cell area. The model Hamiltonian H cont can then be used to specify the contractile contribution to the goal function pðT Þ. To this end, let DH cont ðT Þ denote the contractile contribution to the energy difference entailed by accepting an elementary event T . Following a standard Metropolis algorithm, we then define p cont ðT Þ :¼ exp½ À DH cont ðT Þ=k B T; (S11) where we set the effective thermal energy to k B T 1. The contractile 'energy', Equation S10, is similar to the corresponding energy model commonly used in cellular Potts models (Ouchi et al., 2003). Unlike the standard cellular Potts model, however, where a target area and target perimeter are used to keep the simulated cells from collapsing, the energetic contribution in Equation S10 strictly contracts the cell's body. As will be detailed in the next section, to counteract these contractile forces, we explicitly model cytoskeletal structures within each cell, which provide outward pushing forces to balance cell contraction.

Cytoskeletal structures and focal adhesion
The cytoskeleton plays key roles both in maintaining the mechanical integrity of the cell and in the process of active cell migration (Alberts et al., 2015;Friedl, 2004;Mogilner, 2009). Our model design aims at achieving high computational efficiency to allow for the simulation of very large cell numbers (currently, cell numbers up to Oð10 4 Þ can be achieved at acceptable computation times) and, at the same time, to capture the essential effects of cytoskeletal dynamics to attain meaningful results down to the level of single cells. Thus, instead of accounting for a detailed biochemical description by means of reaction-diffusion networks (Marée et al., 2006;Marée et al., 2012), we resort to a simplified implementation of the most pertinent features of cytoskeletal dynamics. Specifically, we propose a rule-based algorithm to model cytoskeletal structures and to assess the integrated effects of cell polarity, cell contractility and adhesion on the collective dynamics of cells as parts of larger groups.
To this end, we define a scalar field ðx n Þ, x n 2 D ðaÞ , on the domain of each cell a. The local quantity ðx n Þ will be referred to as polarization field and is taken to be a measure for the density of cytoskeletal structures at position x n within the cell's body. The field variable ðx n Þ is dynamically updated as the simulation progresses, reflecting cytoskeletal remodeling. To set up a system of rules underlying the actual implementation of these cytoskeletal remodeling processes, we resort to the following biologically motivated premises: 1. The scalar polarization field is bounded: The dynamics of cytoskeletal remodeling not only depends on the local number (density) of actin monomers and polymers, but also on a multitude of accessory proteins controlling cytoskeleton assembly and disassembly. Several biological factors-including the action of sequestering proteins like thymosin-b4, which act to suppress actin polymerization, limited amounts of nucleating proteins like the activated Arp2/3 complex, and the action of capping proteins-keep the local density of actin filaments bounded. We, therefore, introduce bounds for the polarization field: ðx n ; tÞ 2 ½ 0 À D=2; 0 þ D=2. These bounds are cell-type specific. While the upper bound 0 þ D=2 mainly reflects the limited availability of protein resources, the lower bound 0 À D=2 serves to prevent cells from collapsing. 2. Regulatory proteins affect assembly and disassembly of cytoskeletal structures: The assembly and disassembly of cytoskeletal structures, numerically encoded by ðx n Þ, is regulated by a myriad of accessory proteins. In our computational model we simplify these complex processes by resorting to a single 'bookkeeping variable' which we will refer to as 'regulatory factors'. Its local level is stored as an integer variable Fðx n Þ for each grid site x n 2 D ðaÞ . We use Fðx n Þ to implement the overall action of regulatory cytoskeletal proteins in an effective and collective manner. Specifically, since the formation of lamellipodial structures depends on active nucleation promoting factors (Pollard and Borisy, 2003), we assume that positive levels, Fðx n Þ>0, reflect local conditions in support of network-assembly, whereas negative levels, Fðx n Þ<0, represent predominantly degrading (or disassembly) conditions. For neutral levels, Fðx n Þ ¼ 0, the network gradually restores its rest state.

Feedback between cytoskeletal structures and regulatory factors:
The activities of accessory cytoskeletal proteins which regulate the local levels of cytoskeletal structures are themselves controlled by a number of mechanical and chemical signals received by the cell. Here and in the following, our focus will be on mechanical signals. For example, important regulatory proteins like the Arp2/3 complex are activated locally at the cell membrane, from where they diffuse into the bulk of the cell until they are bound by actin (Kovacs et al., 2002;Pollard and Borisy, 2003;Leckband et al., 2011). Adopting a coarse level of description, this diffusion-degradation dynamics entails a finite range of regulatory proteins, which are activated at the cell's membrane. In our model, we use the integer variable Fðx n Þ to implement this propagation of mechanical information, perceived by cell a at its periphery B ðaÞ , across a certain spatial distance R. The local levels of Fðx n Þ are continuously updated as the MCS loop progresses. The actual update procedure is given by the following set of rules; cf. Appendix 1-figure 2: . if a protrusion event has been accepted at the source site x s 2 B ðaÞ (source cell: a; target cell: b), then for all sites x n within a range R (i.e. kx n À x s k<R) the integer variable signifying regulatory factors is incremented up and down for the protruding and the retracting cell, respectively: Fðx n Þ ! Fðx n Þ þ 1; x n 2 D ðaÞ ; Fðx n Þ À 1; x n 2 D ðbÞ : ( (S12a) . Similarly, if a retraction event has been accepted at the source site x s 2 B ðaÞ , and the (local) cell contact between source cell a and target cell b has remained intact, then within a range R one applies the inverse update rule: Fðx n Þ ! Fðx n Þ À 1; x n 2 D ðaÞ ; Fðx n Þ þ 1; x n 2 D ðbÞ : ( (S12b) . If a retraction event has been accepted at the source site x s 2 B ðaÞ , and in addition the (local) cell contact between source cell a and target cell b has ruptured, then the regulatory factors are reduced only within a range R in the retracting cell: Fðx n Þ ! Fðx n Þ À 1; x n 2 D ðaÞ ; no update; else: ( (S12c) Finally, if the target grid site x t is not occupied by any cell (substrate is indicated by b<0) prior to the elementary event, then only the first two lines in the above update scheme apply.
By virtue of the above update scheme, Equation S12, 'regulatory factors' are continuously distributed across each cell's domain D ðaÞ as the current MCS progresses. At the end of each MCS, the accumulated (local) values of Fðx n Þ are used to update the local values of the polarization field ðx n Þ inside each cell a ! 0 (x n 2 D ðaÞ ): We assume that for positive values, Fðx n Þ>0, there is assembly of cytoskeletal structures and is increased by an amount proportional to the distance of from its upper bound 0 þ D=2: ðx n ; t þ DtÞ ¼ ðx n ; tÞ þ Dt ½ 0 þ D=2 À ðx n ; tÞ; where the time step is defined as Dt 1. Thereby 0 þ D=2 is a fixed point of this map and limits the build-up of cytoskeletal structures. In contrast, for negative values, Fðx n Þ<0, disassembly prevails, and we assume that then tends towards its lower bound 0 À D=2: where the time step is defined as Dt 1. Neutral values, Fðx n ; tÞ ¼ 0, lead to relaxation of towards a resting state ðx n ; t þ DtÞ ¼ ðx n ; tÞ þ Dt ½ 0 À ðx n ; tÞ ; where the time step is defined as Dt 1. The parameter signifies the rate at which cytoskeletal structures respond to the regulatory factors F. For the parameters and cell sizes used in this work ( 0 ¼ Oð100Þ and D ¼ Oð10Þ, and each cell occupying approximately 1000 grid sites) we set ¼ 0:1. Colors indicate local levels of regulatory factors F (blue: F is positive; white: F is zero; red: F is negative; gray: substrate site). Note, in particular, that a substrate grid site has been inserted where cell rupture occurred (row 6, right grid site). The following cases can be distinguished: (i) Grid site x k lies in the zone of influence of only positive (blue circles) or negative (red circles) chemical feedback, in which case the level of regulatory factors is positive or negative, respectively (e.g. red grid sites in row 2, or blue grid sites in row 5). (ii) Grid site x k lies outside of any zone of influence, in which case the level regulatory factors is zero (e.g. white grid sites in row 2). (iii) Grid site x k lies in the zone of influence of equally many positive and negative feedbacks, in which case the level of regulatory factors remains zero (e.g. fourth grid site in row 4). (iv) Grid site x k lies in a zone of predominantly positive or negative feedback, in which case the level of regulatory factors is positive or negative, respectively (e.g. third grid site in row 4). Recall that only the sign of F is of significance to update the cells' polarization field; cf. Equation S13.
After this update procedure for ðx n ; tÞ is completed, all regulatory factors are reset, Fðx n Þ ! 0; 8 n. This prevents 'spurious memory effects' which may arise once the cell's rear reaches its initial leading edge position as time goes on. In essence, resetting regulatory factors upon completion of one MCS implies that the diffusion-degradation dynamics, underlying the distribution of regulatory factors, is fast on the scale of one MCS.
We emphasize that the polarization field ðx n Þ is defined only for grid sites x n 2 D ðaÞ occupied by an actual cell (a ! 0). To allow for spatial variations of substrate properties, we therefore introduce a second scalar field variable 'ðx n Þ, which is defined on the entire computational grid. The scalar field 'ðx n Þ is taken to measure the local density of anchoring points that a cell might use to form focal adhesions. Although one might consider to treat ' as a time-dependent field variable, in this work ' is used to implement static substrate patterns, only. The field 'ðx n Þ is thus initialized once at the beginning of the simulation and kept fixed throughout the entire simulation.
Having introduced the fields ðx n Þ and 'ðx n Þ, we now discuss their impact on the system's dynamics by giving their contribution to the goal function pðT Þ. Suppose that the elementary event T is attempted by a source cell a at source grid site x s 2 B ðaÞ . Further, let x t denote the target grid site and b denote the index of the target cell (where as ususal b ! 0 indicates that x t is occupied by an actual cell, b<0 means that x t exposes substrate). We then define the 'polarization energy' DH cyto ðT Þ as follows: Here, the definition of DH cyto is such that the likelihood of cell protrusions is enhanced if the concentration of cytoskeletal structures at the source grid site, ðx s Þ, is larger than the concentration at the target grid site, ðx t Þ (first row of Equation S14a), and vice versa for cell retractions (second row of Equation S14a). The strength of focal adhesions is taken to be measured by the sum þ '. Their associated 'anchoring effects' (which increase with growing strength of focal adhesions) promote the formation of cell protrusions against unoccupied substrate sites (third row of Equation S14a), and, correspondingly, impede cell retractions (fourth row of Equation S14a). Note, in particular, that the first two rows of Equation S14a can be obtained from a combined evaluation of the lower two rows. For example, if source cell a annexes x t starting from x s , two things need to happen: First, focal adhesions formed by the target cell b must be broken, implying a contribution DH cyto ¼ ðx t Þ þ 'ðx t Þ (fourth row of Equation S14a). Secondly, new focal adhesions are formed by the source cell a, implying a contribution DH cyto ¼ À ðx s Þ þ 'ðx t Þ ½ (third row of Equation S14a). Taking the sum of both contributions gives the expression in the first row of Equation S14a. An analogous line of arguments leads to the expression in the second row of Equation S14a .
The contribution to the goal function pðT Þ due to the polarization energy DH cyto ðT Þ is then defined by p cyto ðT Þ :¼ exp½ À DH cyto ðT Þ=k B T; (S14b) where we set the effective thermal energy to k B T 1. The characteristic 'energy scale' for DH cyto is set by the polarization bounds 0 À D=2 and 0 þ D=2, which turns out to have important implications for collective cell dynamics, as discussed in the main text.

Cell adhesion
To implement the ability of cells to establish cell adhesions across cell-cell interfaces, we use a special form for the respective contribution to the goal function p, which is designed to distinguish between the formation of new and the breakage of existing cellcell adhesion sites.
To this end, we define adhesion matrices B a;b and B 0 a;b quantifying the system's change in 'energy' upon forming a new contact between cells a and b [B a;b ] and upon breaking a preexisting contact between those cells by an 'intruder cell' g 6 ¼ a; b [B 0 a;b ]. In our computational model, we assume that formation of new cell-cell contacts is energetically favored, and that breaking of pre-existing contacts by intruder cells is energetically penalized. The matrix entries of B a;b and B 0 a;b , therefore, have a definite sign, which we take to be positive. The ordering of the cell index pair of B a;b and B 0 a;b is of no physical significance, i.e. the adhesion matrices are symmetric. We also assume that a given cell a does not interact with itself, such that the diagonal elements of the adhesion matrices vanish. Finally, there is no adhesion between cells and empty substrate sites, such that all matrix elements containing a negative cell index vanish. In summary, the adhesion matrices B a;b and B 0 a;b exhibit the following properties: In addition, we assume that the energy cost associated with the breakage of a given cellcell contact exceeds the energetic benefit of its initial formation, i.e.
where equality of both quantities would reproduce the assumption underlying the standard cellular Potts model (Graner and Glazier, 1992;Glazier and Graner, 1993). We shall refer to this property as the 'dissipative nature of cell-cell adhesion'.
To implement the effects of cell-cell adhesion, we compute the 'energy difference' DH adh ðT Þ for any given elementary event T according to the scheme illustrated in Appendix 1-figure 3. One has to distinguish between protrusion and retraction events. First, say that a cell a attempts a protrusion event T pro , involving the source grid site x s 2 B ðaÞ and the target grid site x t 2 B ðbÞ , as illustrated in Appendix 1- figure 3A. In this case, cell a acts as intruder cell, since the depicted protrusion event affects three pre-existing contacts between the target cells b and a 'third party' cell g. Acceptance of the depicted protrusion event would have the following energetically relevant effects: (i) All pre-existing contacts between the target cell b and third party cell g 6 ¼ a; b at the target grid site x t are torn apart.
(ii) New contacts between the source cell a and third party cell g 6 ¼ a; b are established. (iii) The length of the cell contact line between source cell a and target cell b is changed. Altogether, these three effects lead to the following cell adhesion energy difference, where ' is the length of a hexagon edge. The first term in this expression accounts for the (energetically favored) formation of new cellular contacts, as well as for the remodeling of the interface between source cell a and target cell b [points (ii) and (iii)]. The second term measures the (energetically penalized) breaking of pre-existing cell contacts [point (i)] and, therefore, impedes cell a's ability to intrude. Conversely, if source cell a attempts a retraction event T ret , then the same reasoning as the one leading to Equation S16a applies, only this time the elementary event proceeds in reverse, i.e. from the target site x t to the source site x s ; cf. Appendix 1- figure 3A: where ' is the length of a hexagon edge.
cell α protrudes Appendix 1-figure 3. Cell-cell adhesion. (A) Adhesive energy contribution in a cyclic process, where a protrusion of source cell a against target cell b is followed by the inverse retraction event. Both events involve a third party cell g, leading to net energy dissipation after the cyclic process has been completed. Protrusion: (i) Three pre-existing cell-cell contacts between b and g are torn apart (red dashed contacts); (ii) three new contacts between a and g are formed; (iii) the contact length between source cell a and target cell b increases by one unit of length. This implies DH adh ðT pro Þ ¼ ' ð3B 0 b;g À 3B a;g À B a;b Þ. Retraction: (i) Three pre-existing cellcell contacts between a and g are torn apart (red dashed contacts); (ii) three new contacts between b and g are formed; (iii) the contact length between source cell a and target cell b decreases by one unit of length. This implies DH adh ðT ret Þ ¼ ' ð3B 0 a;g À 3B b;g þ B a;b Þ. Altogether, this leads to DH ðcyclÞ adh ¼ DH adh ðT pro Þ þ DH adh ðT ret Þ ¼ ' ð3ðDBÞ a;g þ 3ðDBÞ b;g Þ ! 0, i.e. a (nonnegative) dissipative contribution, whose magnitude depends on the dissipation matrix The condition B 0 >B, Equation S15e, thus implies positive friction associated with cellular shear flows, whose magnitude is proportional to the number of cells sliding past each other. Note that this shear viscosity vanishes for B 0 ¼ B, i.e. for zero dissipation matrix.
We may now use Equation S16a and Equation S16b to illustrate the 'dissipative nature' of adhesive interactions by means of two archetypical examples. First, consider the adhesive energy contribution to any cyclic process. By a cyclic process we mean a sequence of two mutually inverse elementary events, e.g. a protrusion event T pro , which is immediately followed by its inverse retraction event T ret , such that the system's final configuration is identical to its initial configuration. Using Equation S16 we find for the total adhesive energy contribution to a cyclic process: a;cðxjÞ þ ðDBÞ b;cðxjÞ Ã ; (S17) and can therefore conclude that DH ðcyclÞ adh ! 0; where N t denotes the neighborhood of the grid site which temporarily changes its cell index, and where a and b are the indices of the source and target cells involved in the cyclic process; cf. Appendix 1-figure 3A. Since ðDBÞ a;b ! 0, the above adhesive energy contribution is non-negative, thus leading to an amount of energy equal to DH ðcyclÞ adh being dissipated as the cyclic process completes. This leads us to refer to the parameter matrix ðDBÞ a;b as dissipation matrix. Second, consider two (infinitely extended) columns of cells in adhesive contact, sliding past each other. This situation is depicted in Appendix 1- figure 3B, where the left column of cells moves (as a whole) upwards by one grid site, while the right column of cells remains stationary. To assess the adhesive energy contribution along the contact line connecting both cell columns, note that the depicted transformation can be implemented by letting each cell in the left column protrude its leading (i.e. upper) edge by one grid site. For each protruding (source) cell a, this transformation entails to the following energetic effects (cf. discussion above): (i) Two of the pre-existing cell-cell contacts between the source cell's upper neighbor in the left column (target cell b) and the corresponding cell in the right column (third party cell g) get torn apart, leading to an energetic contribution 2B 0 b;g . (ii) In return, two new contacts between the protruding (source) cell a and cell g are being established, leading to a contribution À2B a;g . (iii) Since the length of the contact line between cells in the left column (i.e. between protruding source cell a and retracting target cell b) remains unchanged, there's no further energetic contribution due to adhesive contacts between cells in the left column. Assuming that all cells in the system are of equal types, we write B a;b B and B 0 a;b B 0 (a 6 ¼ b), and therefore, find i.e. a non-negative dissipative contribution per cell. The size of the dissipation parameter DB thus introduces a natural means to tune the system's shear viscosity.
With the above definitions of the adhesive energy changes, Equation S16, we define the contribution of cell adhesion to the goal function pðT Þ as follows: where we set the effective thermal energy to k B T 1. To test for the occurrence of cell rupture, the total energetic cost of each attempted retraction event between two actual cells is evaluated under two different assumptions: First, we assume that the pulling force exerted by the retracting source cell a on target cell b is strong enough to force b to fill the void created at x s (i.e. to accommodate x s ), and call this a regular retraction event T ret . Secondly, we assume that the retraction of source cell a causes all pre-existing cell-cell contacts of cell a at x s to rupture, leaving a free substrate site at x s after the retraction event has been accepted. This latter event will be referred to as rupture event T rup . We then compute the total energy differences

Rupture of cell contacts
under both assumptions (the energy difference associated with accepting T rup can be computed with Equation S16b by using the substrate b ¼ À1 as new target cell) and compare the respective outcomes. If the rupture event is energetically favored over the regular retraction event, i.e. DHðT rup Þ<DHðT ret Þ, then cohesion between cells is weak. In this case, the rupture event T rup , rather than the regular retraction event T ret will be attempted. Otherwise, cohesion is strong and a regular retraction event T ret will be attempted.

Rupture of substrate contacts
In our discussion so far, a cyclic process that follows up a protrusion event T pro with its inverse retraction event T ret , involving a cell a and no third party cells (in other words: no cell-cell contacts are made or broken), will not yield a net energy cost or gain; cf. Equation S14a.
To account for the dissipative nature of cell-substrate contacts, we proceed similarly as when we have considered the disspative nature of cell-cell contacts. We introduce dissipation in substrate adhesion by leaving the Hamiltonian unaltered for protrusion events but adding a penalty for retraction events: DHðT ret Þ ! DHðT ret Þ þ D: Therefore, a cell that adheres to the substrate at some grid point has to pay a cost D to retract from it. In other words, we assume that a fixed amount of energy D is dissipated once the adhesive bonds between a cell and the substrate break.
To keep its overall size across translations, the cell has to gain and lose equal amounts of hexagons, with D as the maximal energy gained by a single gain-and-loss in the absence of dissipation. In the presence of dissipation however, the cell has to pay at least a cost of ð 0 À D=2Þ þ D to detach at an arbitrary location, resulting in D À D as the maximal energy gained by a single gain-and-loss in the presence of dissipation. Thus, while for D ¼ 0 there is no impact of substrate dissipation on cell motility, it will at the latest for D ¼ D completely inhibit (directed) cell migration. Therefore, we study substrate dissipation in the range D 2 ½0; D.
3. Distribute regulatory factors according to Equation S12c.

Cell proliferation and mitosis
While cell proliferation and mitosis play no role in the experimental setup of rotating cell clusters, cell growth and division are observed experimentally in a setup where a sheet of cells expands into free space after removal of a stencil. Therefore, it is essential to include proliferation of cells in the numerical model. How this is done is described in this section.
We distinguish between two phases in the cell cycle, an interphase during which cells roughly double in volume and mitosis, the process of cell division. Even though a further partitioning of the interphase was considered in previous work (Li and Lowengrub, 2014), we do not expect that such a distinction is relevant for our results. In our computational framework cell growth may be implemented by progressively changing any cellular parameter that affects the cell's equilibrium size. The two possible, largely equivalent choices are a successive reduction of the area coupling constant k ðaÞ A or an increase of the average cell polarization ðaÞ 0 . We here employ the first method. We assume that individual cells grow exponentially (Barber et al., 2017) over a well-defined period T g . Additional variability in cell cycle length can be achieved by introducing an additional refractory phase with exponentially distributed waiting times and the average waiting time T r , which we set to T r ¼ 0 in this work. Moreover, we assume that the migratory behavior of a cell should not change significantly as it grows. However, as the cell grows in size by a factor of 2, it also increases its perimeter and the corresponding energy cost for adding new membrane segments roughly by a factor of ffiffi ffi 2 p . Therefore, as we do not scale the polarization field and the resulting energy gains for protrusions during cell growth, we mitigate the increased cost for ruffling the membrane by reducing the perimeter stiffness by a factor of ffiffi ffi 2 p . The quantitative viability of this approach is further discussed in section 'Single cell size'.
To prevent tissue overgrowth, cell proliferation is generally contact inhibited in healthy cells: When the tissue approaches a state where each cell has formed adhesive contacts with the substrate and is completely surrounded by neighbours, cells stop proliferating. In addition, it has been proposed that the pressure or local density in the tissue has a negative impact on the local growth rate (Shraiman, 2005;Ranft et al., 2010). To account for these phenomena in the model, we complement the two cell cycle periods interphase and mitosis by a quiescent cell state during which cell growth is halted. The parameters k A and k P are, therefore, kept constant for a quiescent cell; we denote the corresponding values as k A=0 and k P=0 . There are many possible ways to implement contact inhibition in our computational model. For example, it could be implemented by allowing a quiescent cell to enter the cell cycle triggered by low local cell density, or when a sufficiently large fraction of its membrane length is not in contact with neighbour cells but exposed to free space. In our model it proves numerically advantageous to make a quiescent cell enter the interphase when its area succeeds a certain reference area. We choose this area threshold as A T ¼ r A ref , where the factor r ¼ Oð1Þ relates the threshold size to the equilibrium cell size A ref reached by a free, solitary cell with constant polarization field ¼ 0 . Cells living in a densely packed environment will not exceed the area threshold due to the pressure exerted on them by neighboring cells and can, therefore, not grow. Conversely, cells exposed to free space are more likely to reach this threshold and proliferate. Finally, a growing cell in interphase becomes mitotic after the growth time T g has passed, at which point cell size has roughly doubled with respect to the size in the quiescent period. We assume that cell migration and mitosis are processes that exclude each other. Hence, the positive feedback leading to persistent cell migration is switched off for mitotic cells and the polarization field relaxes to the neutral state 0 according to Equation S13c.
There appears to be no universal set of rules which determine the orientation of the cleavage plane along which cells divide (Minc and Piel, 2012). Rather, for epithelial tissues there are a variety of factors which include local cell geometry and the direction of stress in the tissue (Gibson and Gibson, 2009). Though it is in principle possible to implement any given rule in our computational model, in its present version the axis along which a cell divides is chosen as a random direction through the geometric center of the cell. In case of irregular cell shapes, a separation of the cellular domain into more than two connected components can occur. To prevent violation of topological constraints, in this case the two largest components are considered as descendant cells and the residual grid sites are filled by substrate.
We explicitly account for the finite duration of the mitotic phase T d by keeping the cells in a mitotic state for the aforementioned time period, until the final instantaneous splitting of the cellular domains. After cell division, persistent cell migration of the daughter cells is enabled again. The descendent cells will subsequently re-enter the growing phase if their area exceeds the defined threshold, as mentioned above.
The following list summarizes the steps motivated and explained in the previous paragraphs. These additional steps are performed in a simulation that includes cell proliferation: 1. Assign a state variable s ðaÞ to each cell which encodes the current phase in the cell cycle: . Switch from refractory state to growing state with probability p ¼ 1 À expðÀ1=T r Þ: where the terms in the brackets denote the respective probability.
. Exponential growth in interphase over a period of T g : . Switch from interphase to mitosis: s ðaÞ ðtÞ ¼ 2 for all t 2 ½t À T g ; t ) s ðaÞ ðt þ 1Þ ¼ 3: . During cell division, cell motility is switched off and the polarization field relaxes to the neutral state according to Equation S13c. . Perform cell division, reset area and perimeter stiffness and exit mitotic phase: s ðaÞ ðtÞ ¼ 3 for all t 2 ½t À T d ; t ) divide cell a into cells ða; bÞ . Cell motility is restored after cell division. . If none of the above rules apply, then do not perform any changes.

Numerical computation of stress in a tissue
In the section describing the numerical results on tissue expansion, the stress distribution in the tissue is shown in the kymographs Figure 5(C,G) and Figure 6(C,G). Hereafter we explain how the stress tensor for each cell in the tissue can be computed from the forces acting on the cell's membrane segments in the Monte Carlo simulation. The mean value of the stress tensor in a deformed body can be calculated numerically from the formula s which is a discretized version of the surface integral in Landau et al. (1986). Here, f k is the force acting on the membrane element x k of cell a, com is the position of the element with respect to the center of mass x ðaÞ com of the cell, and the superscripts i and j are Cartesian indices. The forces f k can be computed from the energy differences of all possible protrusion and retraction events originating from x k , where the number sign indicates a sum over substrate grid sites only, i.e. grid sites with cðx l Þ<0, and where DH H cont þ H adh þ H cyto .

Numerical computation of the cell shape
We use two complementary measures for the cell shape. The first is a simple measure for the deviation of an object from a circle (we refer to this as cell extension): It becomes zero if the object is a circle and becomes 1 if the object is a line. The second measure for the cell shape is obtained from a principle components analysis of the cell shape. Specifically, we compute the covariance matrix of the point cloud representing the cell domain D ðaÞ : com denotes the coordinates of element x k of cell a, relative to the cell's center of mass x ðaÞ com ; the superscipts i and j are Cartesian indices. Then, we compute the two eigenvalues l 2 þ (larger eigenvalue) and l 2 À (smaller eigenvalue) of the covariance matrix, which determine the variance of the point distributions along the two principal axes of the cell. Finally, the aspect ratio of the cell is given by l þ =l À .

Parameter screening in silico
In this section we provide additional analysis of the model parameters beyond what is already shown in the main text.
We explore all three rotational phases R 1 , R 2 and R 3 within confinements of varying size and constant cell density. In the R 1 -phase, the cell clusters rotate slowly and frequently reorient their direction of rotation. With increasing cell count, the cell clusters cease to rotate. In the highly coordinated R 2 and R 3 -phases, we find scale-free behavior such that there is always a macroscopic rotation of the whole cell population regardless of the cell count and corresponding confinement size.
We also explore the parameter space of the tissue simulations. There, we find that an increased cell-cell dissipation DB impairs monolayer growth, while at the same time increasing front roughness. Similarly, an increased cell-substrate dissipation D also impairs monolayer growth. In contrast, increasing the maximum cell polarity D improves monolayer growth and also increases front roughness. We thus find that the speed of monolayer expansion depends on whether it is dominated by cell migration or cell proliferation, with the former improving monolayer growth through a better exploration of the cell-free area.

Single cell size
To rationalize our choice of the cell growth algorithm (see section 'Cell proliferation and mitosis'), we have explored the shape and motility of differently sized solitary cells. To this end, we have varied the area stiffness parameter k A for different values of perimeter stiffness k P , while keeping all other parameters constant. We find that the area occupied by the motile cell increases linearly with 1=k A (Appendix 2- figure 1A). In particular, the cell area can be approximated quite well by the area of an immotile and equilibrated cell with uniform ¼ 0 (fit not shown):
Furthermore, we find that with increasing size, and all other parameters constant, cells become rounder, slower, and less persistent (Appendix 2- figure 1B-D). To intuitively explain this phenomenology, let us compare a cell of size A ref with a cell of size r A ref , where r 2 ½1; 2, with the respective area stiffnesses k A and k A =r. While the smaller cell has a perimeter P ref , neglecting geometric effects the larger cell has a larger perimeter of approximately ffiffi r p P ref . Hence, the larger cell has to pay a larger energy cost (roughly by a factor ffiffi r p ) to ruffle its membrane by adding segments and therefore increasing its perimeter. Meanwhile, both the energy gain from the polarization field and the energy cost for increasing the cell area (as A k A ¼ cst:) are the same for both cells. Due to the increased energy cost for adding membrane segments, larger cells finds find it more difficult to polarize, and are therefore rounder, slower and less persistent.
To offset this increased energy cost for adding membrane segments to the cell, one can scale the perimeter stiffness of the larger cell to k P = ffiffi r p , such that P k P » cst. We would then predict that the ratio k 2 P =k A is constant for differently sized cells of similar shape, speed and persistence time. The same relation can also be obtained by realizing that different amounts of grid sites occupied by two otherwise identical cells (in terms of their corresponding Hamiltonian) simply stem from a different discretization of said cells, which is controlled by the parameter k A . Interestingly, we observe such a data collapse for the aspect ratio l þ =l À and the velocity v of the cells onto two respective master curves depending on the ratio k 2 P =k A (Appendix 2- figure 1B,C). While the proposed data collapse for the persistence time of directed migration of a cell (Appendix 2- figure 1D) is somewhat unsatisfactory, this may be owed to the following effect: by keeping R constant we have actually varied the ratio between the area that the signaling molecules typically explore due to diffusion and the area of the cell, R 2 =A. Finally, we speculate that all observed quantities collapse unto respective master curves f ðD ffiffiffiffiffi k A p =k P Þ Á gðR ffiffiffiffiffi k A p Þ.

Single cell shape and dynamics depend on substrate dissipation
We have also studied the effect of cell-substrate dissipation (see section 'Rupture of substrate contacts') on cell morphology and motility. We have varied the substrate dissipation D for different values of maximum cell polarity D and cell perimeter stiffness k P ; however, we were not able to achieve a data collapse in D (Figure 2-figure  supplement 1). We observe that with increasing cell-substrate dissipation, cells become round and cease migrating. This can be illustrated as follows: Consider a situation where the cell conquers a new hexagon at its prospective leading edge. Because the cell on average tends to constrain its area and perimeter while migrating, it consequently needs to lose a different hexagon at its prospective trailing edge. However, this retraction at the trailing edge is energetically penalized and thus cell displacement from its initial position and the positive feedback leading to cell polarization are effectively inhibited. With increasing cell-substrate dissipation, retraction events are further penalized and the cell 'sticks' to the substrate at its trailing edge, preventing persistent motion of the cell. Additionally, to further illustrate the correlation between cell shape and cell migration, we have replotted the values of

Cells in circular confinement
In this section we report on additional parameter studies of the dynamics of cells enclosed in a circular confinement (Figure 4-figure supplement 1, Figure 4-figure supplement 2 and Figure 4-figure supplement 3). Specifically, we investigate how the radius of circular confinement affects the synchronized rotation of the cell population. We performed simulations with a densely populated circular confinement and varied the confinement radius, while keeping the cell density constant. The parameters are chosen such that a population consisting of 4 cells (cf. main text, Figure 4) would rotate in the R 1 , R 2 or R 3 -phase, respectively. We studied the mean angular velocity !ðtÞ ¼ê z Áv a ðtÞ ÂR a ðtÞ kR a ðtÞk 2 * + C ; (S33) averaged over the set C ¼ fa j is not substrateg of all cells in confinement, wherẽ v a ðtÞ ¼ v a ðtÞ À hv a ðtÞi C andR a ðtÞ ¼ R a ðtÞ À hR a ðtÞi C are the velocity and position of cell a relative to the cell cluster, respectively. In the lowly polarizable R 1 -phase, small cell populations rotate in a highly synchronized way, and rotation is maximized for populations of 7 cells per confinement (Figure 4-figure  supplement 1A). As can be inferred from the time traces and snapshots (Figure 4-figure supplement 1B, C), cells in small populations all synchronously move in the same direction at a given time and randomly switch between clockwise and counter-clockwise rotation; the switching rate decreases with increasing size of the cell population. Upon increasing the cell count and concomitantly the confinement size, global rotation of the cell population gradually vanishes (Figure 4-figure supplement 1A).
Unlike in the R 1 -phase, we observe that in the highly polarizable R 2 and R 3 -phases populations of all sizes rotate in a highly synchronized way (Figure 4-figure supplement 2A and Figure 4-figure supplement 3A). There, the dependence of ! j j h i on the population size N can be fitted by a power law of the form ! j j h i / N À1=2 / r À1 0 . This inverse proportionality between the average angular velocity and the confinement size r 0 implies total rotational order, with every cell moving at a constant velocity v rot j j » 0:008 (v rot ?R). Furthermore, in the R 2 , and R 3 -phases we have only scarcely observed switching of the rotational direction of cell clusters; e.g. for 4-cell clusters in the R 3 -phase.
Interestingly, fluctuations in the angular velocity (s ! ) change in a highly non-monotonic fashion with the cell count and concomitantly the confinement size. Certain cell counts exhibit especially high fluctuations of the mean angular velocity (e.g. 5 cells in the R 1 -phase, see

Velocity and roughness of spreading tissue
We have studied the velocity and roughness of spreading tissue, while varying cell-cell dissipation DB, cell-substrate dissipation D and maximum cell polarity D.
First, let us introduce the observables that we are interested in. Let X >=< be the sets of xcoordinates of the left and right outermost edges of the cell sheet. Our in silico setup is axially symmetric with respect to the y-axis. This initial symmetry persists, as the cell fronts advance towards the cell-free area with the same average speed. Hence, it is not needed to consider the two cell fronts separately, and we can instead consider the set of unsigned front positions X :¼ absðX >=< Þ. Then, we define the average front position as x F :¼ EðXÞ and the front roughness as s 2 F :¼ VarðXÞ. In particular, we study the total growth of the tissue over the course of 500 MCS, which is captured by the maximal position of the front, maxðx F Þ, as well as the maximal roughness of the front maxðs F Þ. We have chosen our parameters such that a cell takes a total amount of 200 MCS to divide, provided that it exceeds the threshold size of a solitary reference cell A ref . Because the first daughter cells may only appear after 200 MCS have passed, we exclude this initial period from the measurements of the maximal front position and roughness, respectively. Additionally, we provide some exemplary time traces of the front evolution ( Figure 5-figure supplement B,D,F).
First, we investigated how the monolayer expansion and front roughness depend on cellsubstrate dissipation, DB ( Figure 5-figure supplement 1A, B). Our simulations show that the cell sheet expands slower with increasing cell-cell dissipation DB ( Figure 5-figure  supplement 1A, B), because the dissipation penalizes cells sliding past each other. At the same time, the cell sheet also becomes slightly rougher with increasing cell-cell dissipation DB (inset of Figure 5-figure supplement 1A).
We also investigated how the monolayer expansion and front roughness depend on cellsubstrate dissipation, D ( Figure 5-figure supplement 1C, D). Before we turn to the monolayer, let us recall the observed single-cell behavior in the previous section (see section 'Single cell shape and dynamics depend on substrate dissipation'): for high enough cellsubstrate dissipation D (typically of the same order of magnitude as the maximum cell polarity D) cell migration is switched off (Figure 2-figure supplement 1). Extrapolating the singlecell results, we expect that the same holds also for collectives of cells and that cell migration does not play a role for high cell-substrate dissipation. Indeed, with increasing cell-substrate dissipation, the monolayer expands slower, until this effect appears to saturate at a threshold value D $ » 5 ( Figure 5-figure supplement 1C, D). Following this line of argument, monolayer growth is slowed down if we suppress cell migration and thus move the cell monolayer towards a proliferation-dominated mode of expansion.
What about the inverse? Is the monolayer growth increased if we enhance cell migration and thus move the cell monolayer towards a migration-dominated mode of expansion? To test this hypothesis, we have analyzed how the monolayer growth and front roughness depend on the maximum cell polarity D. As predicted, monolayer growth increases with the maximum cell polarity D ( Figure 5-figure supplement 1E, F), because an increased amount of cells exceed the threshold size to switch to mitosis (cf. the stretching of bulk cells in the monolayer in Figure 5B). Additionally, we also find that the front roughness increases with increasing maximum cell polarity D.