Bridging the gap between singlecell migration and collective dynamics
Abstract
Motivated by the wealth of experimental data recently available, we present a cellularautomatonbased modeling framework focussing on highlevel cell functions and their concerted effect on cellular migration patterns. Specifically, we formulate a coarsegrained description of cell polarity through selfregulated 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). Singlecell migration enables cells to move towards and between tissue compartments – a process that plays an important role in the inflammationinduced 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) Cellcell cohesion and coupling mediated by adherensjunction 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 singlecell 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 singlecell 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 coarsegrained 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 singlecell migration and collective dynamics, namely phasefield 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 $\mathcal{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 singlecell 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 woundhealing assay (Poujade et al., 2007; Trepat et al., 2009; SerraPicamal et al., 2012), and show that the model exhibits the recurring mechanical waves observed experimentally (SerraPicamal et al., 2012), a feature which we attribute to the coupling between cellsheet expansion and celldensityinduced growth inhibition.
Computational model
Model geometry
We consider cells that adhere to a twodimensional 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 $z$direction. 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, $\mathscr{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.
Coarsegrained 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
where ${\kappa}_{A}$ and ${\kappa}_{P}$ are celltypespecific 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.
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 timedependent and spatially resolved internal field for each cell, $\u03f5(\mathbf{\mathbf{x}},t)$. This polarization field emulates the mass of forcegenerating cytoskeletal structures in the associated hexagon, at position $\mathbf{\mathbf{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 celltypespecific bounds for the polarization field: $\u03f5(\mathbf{x},t)\in [{\u03f5}_{0}\mathrm{\Delta}\u03f5/2,\text{}{\u03f5}_{0}+\mathrm{\Delta}\u03f5/2]$, where ${\u03f5}_{0}$ is the average polarization field and $\mathrm{\Delta}\u03f5$ is the maximum cell polarity.
Active selfregulation 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 reactiondiffusion 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 downregulate 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):
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$.
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, $\u03f5\to {\u03f5}_{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 $\mu $. 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:
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.
Cell migration couples mechanochemically to a vector field (Marée et al., 2006; Ziebert et al., 2012), if stresses in the cell are anisotropic.
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 twofold:
Cell adhesion supports growth of cellcell and cellmatrix contacts and may thus be described in terms of effective surface energies. In our computational model, cellmatrix contacts are readily accounted for by the polarization field, $\u03f5$. In addition, we associate the formation of cellcell adhesions with an energy benefit $B$, which we call cellcell adhesion parameter.
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 cellcell adhesion, $B+\mathrm{\Delta}B>B$, exceeds the gain from forming a new cellcell adhesion. Then, the dissipative nature of cellcell adhesions is accounted for by the cellcell friction parameter $\mathrm{\Delta}B$. Similarly, cellmatrix contacts can also provide a source of dissipation, which is further discussed in Appendix 2.
Environmental cues
The polarization field, $\u03f5$, readily includes contributions from cellsubstrate adhesions, which are locally up or downregulated by the cell. These cellsubstrate 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 nogoareas, where the cell adheres less (or cannot adhere at all). To replicate such environmental cues, we introduce a second scalar field $\phi (\mathbf{\mathbf{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, $\phi \ll 0$, to the polarization field ($\u03f5\to \u03f5+\phi $), that a cell has to pay for trespassing onto a nogoarea. However, it is equally valid to treat $\phi $ as a multiplicative constant modulating the polarization field ($\u03f5\to \phi \u03f5$), where $\phi =0$ models a local inability of the cell to attach to the substrate. Analogously to cellcell contacts, we account for the dissipative nature of cellsubstrate 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 threestate 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}_{\text{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 (${\kappa}_{A}$ and ${\kappa}_{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 (Jones et al., 2018; 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 ($\mathrm{\Delta}\u03f5\to 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 reinitialize migration from an unpolarized state. For a detailed and more technical description we refer the interested reader to Appendix 1.
Results
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:
The cytoskeletal update rate endows the cellular Potts model with a reference timescale and determines the temporal discretization. In this study, we have set $\mu =0.1$.
The average polarization field ${\u03f5}_{0}$ encodes the energy gain for creating new cellsubstrate adhesions, while the area stiffness ${\kappa}_{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 ${\u03f5}_{0}/{\kappa}_{A}$. If we use a desired cell area as reference value, then the ratio ${\u03f5}_{0}/{\kappa}_{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 ${\u03f5}_{0}=225$ and the area stiffness to ${\kappa}_{A}=0.18$.
In cellular Potts models, which are MonteCarlo simulations, the reference energy of fluctuations is determined by an effective temperature. In this study, we have set ${k}_{\text{B}}T\equiv 1$.
Furthermore, we used a large computational grid with $9\cdot {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 cellsubstrate dissipation: cell polarizability $\mathrm{\Delta}\u03f5$, cell contractility ${\kappa}_{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 ${\kappa}_{P}$ and maximum cell polarity $\mathrm{\Delta}\u03f5$ 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 orientation $\hat{\mathbf{v}}(t)\equiv \mathbf{v}(t)/\Vert \mathbf{v}(t)\Vert$ ($\mathbf{\mathbf{v}}$: cell velocity) and (geometrical) center of mass position $\mathbf{\mathbf{R}}(t)$ during a total simulation time of $T}_{\text{sim}}={10}^{4$ MonteCarlo steps (MCS). For each set of parameters, we performed $100$ statistically independent simulations, from which we computed the mean squared displacement, $\text{MSD}(\tau )\equiv \u27e8{[\mathbf{\mathbf{R}}(t+\tau )\mathbf{\mathbf{R}}(t)]}^{2}\u27e9$, and the normalized velocity autocorrelation function, $C(\tau )\equiv \u27e8\hat{\mathbf{v}}(t+\tau )\cdot \hat{\mathbf{v}}(t)\u27e9$. Here, $\u27e8\dots \u27e9$ 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(\tau )\propto {e}^{\tau /{\tau}_{p}}$ (inset of Figure 2A). We determined the persistence time of directed migration, ${\tau}_{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. Surprisingly, for each of these variables we found a master curve that only depends on the ratio between cell polarizability and cell contractility, $\mathrm{\Delta}\u03f5/{\kappa}_{P}$ (Figure 2B,D,E). This data collapse suggests $\mathrm{\Delta}\u03f5/{\kappa}_{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, $\mathrm{\Delta}\u03f5/{\kappa}_{P}\approx 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 myosininduced contraction) of the cytoskeleton or spatially regulated cellsubstrate 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, $\mathrm{\Delta}\u03f5$, and contractility, ${\kappa}_{P}$, are equally reduced by a blebbistatindependent prefactor, then the specific polarizability, $\mathrm{\Delta}\u03f5/{\kappa}_{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 MonteCarlo step) on the persistence of singlecell 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 frontrear 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 ($\mathrm{\Delta}\u03f5$), 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 3A,D; blue diamonds, green pentagons, yellow squares and black circles). To find an intuitive explanation for these observations, we inspected timelapse videos of a cell at high polarizability ($\mathrm{\Delta}\u03f5/{\kappa}_{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 nonzero 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 keratocytelike 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 counterintuitive 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 $\u27e8c\u27e9=\u27e8\Vert {\mathrm{\partial}}_{s}\hat{\mathbf{v}}(s)\cdot \hat{\mathbf{v}}(s)\Vert \u27e9$, where $s$ is the contour length along the corresponding trajectory. Here, we averaged the tangent vector $\hat{\mathbf{v}}(s)$ over $10$ MonteCarlo steps to integrate out fluctuations that occur on short timescales (the internal dynamics of the cell has an intrinsic time scale of $10$ MonteCarlo steps due to our choice of the cytoskeletal update rate, $\mu =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 $\phi (\mathbf{\mathbf{x}})=0$ inside a radius ${r}_{0}$ and $\phi (\mathbf{\mathbf{x}})\to \mathrm{\infty}$ 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 ${\mathbf{\mathbf{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 $\omega (t)={\hat{\mathbf{e}}}_{z}\cdot \u27e8\stackrel{~}{\mathbf{v}}(t)\times \stackrel{~}{\mathbf{R}}(t)/\Vert \stackrel{~}{\mathbf{R}}(t){\Vert}^{2}{\u27e9}_{\mathcal{\mathcal{C}}}$. Here, $\hat{\mathbf{e}}}_{z$ is the outofplane unit vector, ${\u27e8\mathrm{\dots}\u27e9}_{\mathcal{C}}$ denotes an average with respect to the cell population, and $\stackrel{~}{\mathbf{v}}(t)=\mathbf{v}(t)\u27e8\mathbf{v}(t){\u27e9}_{\mathcal{\mathcal{C}}}$ as well as $\stackrel{~}{\mathbf{R}}=\mathbf{R}(t)\u27e8\mathbf{R}(t){\u27e9}_{\mathcal{\mathcal{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, $\omega (t)$, and the average cell perimeter $P(t)\equiv {\u27e8{P}_{\alpha}(t)\u27e9}_{\mathcal{C}}$ were then used to characterize the statistics of collective cell rotation. For each specific choice of simulation parameters, we monitored $\omega (t)$ and $P(t)$ for a set of $100$ statistically independent systems, each of which was observed over $T}_{\text{sim}}={10}^{4$ MCS. From these data, we then computed the mean overall rotation speed $\u27e8\omega \u27e9$, its standard deviation ${\sigma}_{\omega}$, and the standard deviation of the cell perimeter, ${\sigma}_{P}$.
Figure 4 illustrates the characteristic properties of collective cell rotations in systems containing $\mathcal{C}=4$ cells endowed with varying maximum cell polarity $\mathrm{\Delta}\u03f5$ and varying cell contractility ${\kappa}_{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 $\mathrm{\Delta}\u03f5/{\kappa}_{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 ($\mathrm{\Delta}\u03f5/{\kappa}_{P}\approx 450$ in Figure 4A), the rotation speed $\u27e8\omega \u27e9$ (purple curves in Figure 4A) vanishes and the cells are immobile. In this regime, which we term the stagnation phase, or $\mathcal{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 ${\mathcal{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 (${\sigma}_{P}$) and blue (${\sigma}_{\omega}$) curves in Figure 4A. Now, cytoskeletal forces are sufficiently large to establish actual membrane 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 (${\mathcal{R}}_{2}$phase), the cellular systems reach a regime of enduring rotational motion, where $\u27e8\omega \u27e9$ varies linearly with the local specific polarizability, and where ${\sigma}_{P}$ and ${\sigma}_{\omega}$ 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 frontrear 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 (${\mathcal{R}}_{3}$phase), the system’s dynamics is dominated by cytoskeletal forces and the rotational speed $\u27e8\omega \u27e9$ 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 ${\sigma}_{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 ${\sigma}_{\omega}$; blue curve in Figure 4A).
Tissuelevel 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 (SerraPicamal 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 $y$direction 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 ${\u03f5}_{0}=35$. For each set of parameters, we performed and averaged $100$ independent simulations.
We first investigated how a densely packed pregrown tissue of mitostatic cells with high polarizability (large $\mathrm{\Delta}\u03f5$) expands into cellfree space upon removal of the confining boundaries at the tissue’s lateral edges (Figure 5A). As the cells migrate into the cellfree 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 Xshaped pattern in the kymograph of the total mechanical stresses ${\u27e8{\sigma}_{xx}\u27e9}_{y}$ (Figure 5C). This profile closely resembles the first period of mechanical waves observed experimentally (SerraPicamal 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 Xshaped 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 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 migrationdominated monolayer expansion and celldensitydependent cell growth and division results in a periodic recurrence of the Xshaped stress pattern (Figure 5G), closely resembling the pattern observed in experiments (SerraPicamal 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 wavelike behavior, ultimately emerge from cell polarization and the ensuing active cell migration. Therefore, 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 cellcell 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 cellcell dissipation and slight decrease of cellcell 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 cellfree 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 coarsegrained 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 $\mathcal{O}({10}^{4})$ cells. We have used the model to study the transition from singlecell 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 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 singlecell 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 tractionforce patterns and an Xshape in the corresponding kymograph. Additionally, a celldensitydependent cell growth leads to a periodic recurrence of these tractionforce patterns in a cycle of migrationdominated 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 intercellular 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.
Appendix 1
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 modelspecific 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, spacefilling lattice with lattice vectors $\left\{{\mathbf{x}}_{i}\right\}}_{i=1,\dots ,N$. Each lattice vector $\mathbf{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 $\left\{{\mathbf{x}}_{i}\right\}}_{i=1,\dots ,N$, such that each grid site is a hexagon, which is surrounded by $6$ nearestneighbor sites that define the neighborhood $\mathcal{\mathcal{N}}}_{k$ of $\mathbf{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 $\ell$ and the twodimensional volume (area) $a$ of each hexagonal grid site: $\ell \phantom{\rule{thinmathspace}{0ex}}=\phantom{\rule{thinmathspace}{0ex}}1/\sqrt{3}$ and $a\phantom{\rule{thinmathspace}{0ex}}=\phantom{\rule{thinmathspace}{0ex}}3\sqrt{3}\phantom{\rule{thinmathspace}{0ex}}{\ell}^{2}/2$.
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({\mathbf{\mathbf{x}}}_{k})$ gives the index of the cell occupying ${\mathbf{\mathbf{x}}}_{k}$. Here and in the following, we use latin indices to reference lattice sites, and greek indices to reference cells. The set ${\mathcal{D}}^{(\alpha )}$ used to represent the spatial extension of cell $\alpha $, will be referred to as the domain of cell $\alpha $. In our model, each grid site ${\mathbf{\mathbf{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 ${\mathbf{\mathbf{x}}}_{k}$ is numerically implemented by negative values of the indicator function, $c({\mathbf{\mathbf{x}}}_{k})<0$. Following this terminology, the area and the perimeter of cell $\alpha $ are given by:
Model dynamics
Protrusion and retraction of cells
Biological cells are highly dynamic entities which constantly change shape and move around in space. To reflect this dynamic behavior computationally, the domain ${\mathcal{D}}^{(\alpha )}$ of cell $\alpha $ changes over time. The evolution of cell shape and position, as represented by ${\mathcal{D}}^{(\alpha )}$, proceeds via a succession of elementary events. In our numerical model, elementary events come in one of two basic flavors: protrusion events and retraction events. During a protrusion event, cell $\alpha $ (referred to as source cell) incorporates one grid site ${\mathbf{\mathbf{x}}}_{t}$ (referred to as target grid site) from its neighborhood ${\mathcal{N}}^{(\alpha )}$,
thereby increasing its cellular domain by one grid site. Here, the neighborhood of cell $\alpha $, ${\mathcal{N}}^{(\alpha )}$, is defined as
During a retraction event, source cell $\alpha $ expels one of its membrane grid sites ${\mathbf{\mathbf{x}}}_{s}\in {\mathcal{B}}^{(\alpha )}$,
where the set of membrane grid sites ${\mathcal{B}}^{(\alpha )}$ 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 $\alpha $ into the domain of cell $\beta $ (referred to as target cell) must be accompanied by a corresponding retraction event ${\mathcal{D}}_{\text{old}}^{(\beta )}\to {\mathcal{D}}_{\text{new}}^{(\beta )}={\mathcal{D}}_{\text{old}}^{(\beta )}\setminus \{{\mathbf{\mathbf{x}}}_{t}\}$, where ${\mathbf{\mathbf{x}}}_{t}$ denotes the target grid site annexed by cell $\alpha $. We emphasize, however, that the reverse is not generally true. If source cell $\alpha $ retracts, i.e. loses one of its boundary grid sites ${\mathbf{\mathbf{x}}}_{s}\in {\mathcal{B}}^{(\alpha )}$, the lost grid site ${\mathbf{\mathbf{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 $\alpha $ exerts a pulling force on one of its neighboring cells $\beta $ (the target cell) and forces the target cell to fill the emerging void at ${\mathbf{\mathbf{x}}}_{s}$, i.e. ${\mathcal{D}}_{\text{old}}^{(\beta )}\to {\mathcal{D}}_{\text{new}}^{(\beta )}={\mathcal{D}}_{\text{old}}^{(\beta )}\cup \{{\mathbf{\mathbf{x}}}_{s}\}$, where ${\mathbf{\mathbf{x}}}_{s}$ denotes the grid site lost by cell $\alpha $. On the other hand, if adhesion between cells is weak, then retraction of the source cell $\alpha $ can lead to a rupture of preexisting cell contacts between $\alpha $ and other cells at ${\mathbf{\mathbf{x}}}_{s}$, such that the lost grid site ${\mathbf{\mathbf{x}}}_{s}$ becomes free space [$c({\mathbf{\mathbf{x}}}_{s})=\alpha \ge 0\to c({\mathbf{\mathbf{x}}}_{s})<0$]. Details on the actual implementation of cell rupture are discussed in section ‘Rupture of cell contacts’.
MonteCarlo scheme
In the spirit of a standard MonteCarlo scheme, the actual simulation proceeds via a succession of MonteCarlo steps, where each MonteCarlo step (MCS) propagates the state of the simulated cell population from time $t$ to time $t+\mathrm{\Delta}t$, where we set the time step to $\mathrm{\Delta}t\equiv 1$. One MCS consists in a series of attempts to perform elementary events, originating from randomly chosen membrane grid sites of randomly chosen cells. The duration of one MCS, i.e. the actual number of attempted elementary events, is chosen such that each of the cells’ membrane segments is given the opportunity to attempt, on average, one elementary event per MCS. During each MCS, cell domains ${\mathcal{D}}^{(\alpha )}$ as well as the numerical values of cell areas ${A}_{\alpha}$ and perimeters ${P}_{\alpha}$ are updated ‘on the fly’, while the cells’ polarization fields are updated only once at the end of each MCS; cf. section ‘Cytoskeletal structures and focal adhesion’ for the details of this update rule. The simulation then proceeds along the following MonteCarlo scheme:
Initialize the cell population and define the duration of the simulation, i.e. the number of MCS, $T}_{\text{sim}$, to be performed.
Set the simulation time $t=0$.
Perform the next MCS; this step is further detailed below.
Update polarization fields (cf. section ‘Cytoskeletal structures and focal adhesion’).
Set $t=t+\mathrm{\Delta}t$, where $\mathrm{\Delta}t\equiv 1$.
Repeat steps 3–5 while $t<{T}_{\text{sim}}$.
The implementation of a MCS, i.e. the sequence of elementary events, is based on the following general considerations:
Choice of source and target grid sites. Each elementary event $\mathcal{T}$ originates from a membrane grid site ${\mathbf{\mathbf{x}}}_{s}\in {\mathcal{B}}^{(\alpha )}$ of some cell $\alpha $, 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 ${\mathbf{\mathbf{x}}}_{s}$ and which is not currently occupied by cell $\alpha $: ${\mathbf{\mathbf{x}}}_{t}\in {\mathcal{N}}_{s}\setminus {\mathcal{D}}^{(\alpha )}$. In what follows, this additional grid site ${\mathbf{\mathbf{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 $\beta $, 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.
MonteCarlo 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 MonteCarlo procedure, the probability of occurrence of elementary events $\mathcal{T}$ is determined by a goal function $p(\mathcal{T})$ [cf. point (iii) below]. However, since elementary events come in two basic flavors, protrusions ${\mathcal{T}}_{\text{pro}}$ and retractions ${\mathcal{T}}_{\text{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 $\mathbf{x}}_{t$, or retracts from $\mathbf{x}}_{s$. Then, a decision is made with equal probabilities as to whether one attempts $\mathcal{\mathcal{T}}}_{\text{pro}$ or $\mathcal{\mathcal{T}}}_{\text{ret}$. In a second step, the goal function $p$ is used to compute the occurrence probability of the attempted event $\mathcal{\mathcal{T}}$. Finally, this elementary event $\mathcal{\mathcal{T}}$ is being accepted with probability $p(\mathcal{\mathcal{T}})$.
Choice of the goal function $p\mathit{}\mathrm{(}\mathrm{T}\mathrm{)}$. As has been detailed above, we use a goal function $p(\mathcal{T})$ to control the occurrence and acceptance of elementary events $\mathcal{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 cellcell 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(\mathcal{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:
Determine the current number of trials per MonteCarlo step (MCS), $K\phantom{\rule{thinmathspace}{0ex}}=\phantom{\rule{thinmathspace}{0ex}}\sum _{\alpha}{P}_{\alpha}/\ell$, and set the trial counter $n=0$.
With equal probability, choose a cell membrane segment (cf. solid black line in Appendix 1—figure 1) from a random cell $\alpha$ of the cell population. Because the cell membrane represents the border between lattice sites occupied by cell $\alpha $ and unoccupied by cell $\alpha $, the chosen membrane segment automatically defines the source grid site $\mathbf{x}}_{s}\in {\mathcal{\mathcal{B}}}^{(\alpha )$ and the corresponding target grid site $\mathbf{x}}_{t}\in {\mathcal{\mathcal{N}}}^{(\alpha )}\cap {\mathcal{\mathcal{N}}}_{s$.
With equal probability, choose whether to attempt a protrusion event ($\mathcal{\mathcal{T}}}_{\text{pro}$) or a retraction event ($\mathcal{\mathcal{T}}}_{\text{ret}$).
Compute the prospective acceptance probability $p({\mathcal{T}}_{\text{pro/ret}})$ corresponding to the attempted event, and decide whether to accept the attempted event on the basis of this probability.
If the attempted elementary event has been accepted, then update the cellular domains of source cell $\alpha $ and opponent cell $\beta $; for details see section ‘Cell domain update routine’.
If $n<K$, set $n\to 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(\mathcal{T})$ of an elementary event $\mathcal{T}$. Specifically, our model takes into account cell contractility, the assembly and disassembly of cytoskeletal structures, cellcell 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}_{\text{cyto}}(\mathcal{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’
with positive coupling constants $\kappa}_{P}^{(\alpha )$ and $\kappa}_{A}^{(\alpha )$ characterizing the contractility of cell $\alpha $; for empty substrate sites ($\alpha <0$) we set ${\kappa}_{P}^{(\alpha )}={\kappa}_{A}^{(\alpha )}=0$. According to Equation S10, the cell’s ‘contractile energy’ increases with increasing cell perimeter and increasing cell area. The model Hamiltonian ${\mathscr{H}}_{\text{cont}}$ can then be used to specify the contractile contribution to the goal function $p(\mathcal{T})$. To this end, let $\mathrm{\Delta}{\mathscr{H}}_{\text{cont}}(\mathcal{T})$ denote the contractile contribution to the energy difference entailed by accepting an elementary event $\mathcal{T}$. Following a standard Metropolis algorithm, we then define
where we set the effective thermal energy to ${k}_{\text{B}}T\equiv 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 $\mathcal{\mathcal{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 reactiondiffusion 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 rulebased 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 $\u03f5({\mathbf{x}}_{n})$, $\mathbf{x}}_{n}\in {\mathcal{\mathcal{D}}}^{(\alpha )$, on the domain of each cell $\alpha$. The local quantity $\u03f5({\mathbf{x}}_{n})$ will be referred to as polarization field and is taken to be a measure for the density of cytoskeletal structures at position $\mathbf{x}}_{n$ within the cell’s body. The field variable $\u03f5({\mathbf{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:
The scalar polarization field $\u03f5$ 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β4, 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: $\u03f5({\mathbf{x}}_{n},t)\in [{\u03f5}_{0}\mathrm{\Delta}\u03f5/2,\text{}{\u03f5}_{0}+\mathrm{\Delta}\u03f5/2]$. These bounds are celltype specific. While the upper bound ${\u03f5}_{0}+\mathrm{\Delta}\u03f5/2$ mainly reflects the limited availability of protein resources, the lower bound ${\u03f5}_{0}\mathrm{\Delta}\u03f5/2$ serves to prevent cells from collapsing.
Regulatory proteins affect assembly and disassembly of cytoskeletal structures: The assembly and disassembly of cytoskeletal structures, numerically encoded by $\u03f5({\mathbf{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({\mathbf{\mathbf{x}}}_{n})$ for each grid site ${\mathbf{\mathbf{x}}}_{n}\in {\mathcal{D}}^{(\alpha )}$. We use $F({\mathbf{\mathbf{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({\mathbf{\mathbf{x}}}_{n})>0$, reflect local conditions in support of networkassembly, whereas negative levels, $F({\mathbf{\mathbf{x}}}_{n})<0$, represent predominantly degrading (or disassembly) conditions. For neutral levels, $F({\mathbf{\mathbf{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 diffusiondegradation 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({\mathbf{x}}_{n})$ to implement this propagation of mechanical information, perceived by cell $\alpha $ at its periphery ${\mathcal{B}}^{(\alpha )}$, across a certain spatial distance $R$. The local levels of $F({\mathbf{\mathbf{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 $\mathbf{x}}_{s}\in {\mathcal{\mathcal{B}}}^{(\alpha )$ (source cell: $\alpha$; target cell: $\beta$), then for all sites $\mathbf{x}}_{n$ within a range $R$ (i.e. $\Vert {\mathbf{x}}_{n}{\mathbf{x}}_{s}\Vert <R$) the integer variable signifying regulatory factors is incremented up and down for the protruding and the retracting cell, respectively:
Similarly, if a retraction event has been accepted at the source site $\mathbf{x}}_{s}\in {\mathcal{\mathcal{B}}}^{(\alpha )$ , and the (local) cell contact between source cell $\alpha$ and target cell $\beta$ has remained intact, then within a range $R$ one applies the inverse update rule:
If a retraction event has been accepted at the source site $\mathbf{x}}_{s}\in {\mathcal{\mathcal{B}}}^{(\alpha )$ , and in addition the (local) cell contact between source cell $\alpha$ and target cell $\beta$ has ruptured, then the regulatory factors are reduced only within a range $R$ in the retracting cell:
Finally, if the target grid site ${\mathbf{\mathbf{x}}}_{t}$ is not occupied by any cell (substrate is indicated by $\beta <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 ${\mathcal{D}}^{(\alpha )}$ as the current MCS progresses. At the end of each MCS, the accumulated (local) values of $F({\mathbf{\mathbf{x}}}_{n})$ are used to update the local values of the polarization field $\u03f5({\mathbf{\mathbf{x}}}_{n})$ inside each cell $\alpha \ge 0$ (${\mathbf{\mathbf{x}}}_{n}\in {\mathcal{D}}^{(\alpha )}$): We assume that for positive values, $F({\mathbf{x}}_{n})>0$, there is assembly of cytoskeletal structures and $\u03f5$ is increased by an amount proportional to the distance of $\u03f5$ from its upper bound ${\u03f5}_{0}+\mathrm{\Delta}\u03f5/2$:
where the time step is defined as $\mathrm{\Delta}t\equiv 1$. Thereby ${\u03f5}_{0}+\mathrm{\Delta}\u03f5/2$ is a fixed point of this map and limits the buildup of cytoskeletal structures. In contrast, for negative values, $F({\mathbf{\mathbf{x}}}_{n})<0$, disassembly prevails, and we assume that $\u03f5$ then tends towards its lower bound ${\u03f5}_{0}\mathrm{\Delta}\u03f5/2$:
where the time step is defined as $\mathrm{\Delta}t\equiv 1$. Neutral values, $F({\mathbf{\mathbf{x}}}_{n},t)=0$, lead to relaxation of $\u03f5$ towards a resting state
where the time step is defined as $\mathrm{\Delta}t\equiv 1$. The parameter $\mu $ signifies the rate at which cytoskeletal structures respond to the regulatory factors $F$. For the parameters and cell sizes used in this work (${\u03f5}_{0}=\mathcal{\mathcal{O}}(100)$ and $\mathrm{\Delta}\u03f5=\mathcal{\mathcal{O}}(10)$, and each cell occupying approximately $1000$ grid sites) we set $\mu =0.1$.
After this update procedure for $\u03f5({\mathbf{\mathbf{x}}}_{n},t)$ is completed, all regulatory factors are reset, $F({\mathbf{x}}_{n})\to 0,\text{}\mathrm{\forall}\phantom{\rule{thinmathspace}{0ex}}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 diffusiondegradation dynamics, underlying the distribution of regulatory factors, is fast on the scale of one MCS.
We emphasize that the polarization field $\u03f5({\mathbf{\mathbf{x}}}_{n})$ is defined only for grid sites ${\mathbf{\mathbf{x}}}_{n}\in {\mathcal{D}}^{(\alpha )}$ occupied by an actual cell ($\alpha \ge 0$). To allow for spatial variations of substrate properties, we therefore introduce a second scalar field variable $\phi ({\mathbf{\mathbf{x}}}_{n})$, which is defined on the entire computational grid. The scalar field $\phi ({\mathbf{\mathbf{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 $\phi $ as a timedependent field variable, in this work $\phi $ is used to implement static substrate patterns, only. The field $\phi ({\mathbf{\mathbf{x}}}_{n})$ is thus initialized once at the beginning of the simulation and kept fixed throughout the entire simulation.
Having introduced the fields $\u03f5({\mathbf{\mathbf{x}}}_{n})$ and $\phi ({\mathbf{\mathbf{x}}}_{n})$, we now discuss their impact on the system’s dynamics by giving their contribution to the goal function $p(\mathcal{T})$. Suppose that the elementary event $\mathcal{T}$ is attempted by a source cell $\alpha $ at source grid site ${\mathbf{\mathbf{x}}}_{s}\in {\mathcal{B}}^{(\alpha )}$. Further, let ${\mathbf{\mathbf{x}}}_{t}$ denote the target grid site and $\beta $ denote the index of the target cell (where as ususal $\beta \ge 0$ indicates that ${\mathbf{\mathbf{x}}}_{t}$ is occupied by an actual cell, $\beta <0$ means that ${\mathbf{\mathbf{x}}}_{t}$ exposes substrate). We then define the ‘polarization energy’ $\mathrm{\Delta}{\mathscr{H}}_{\text{cyto}}(\mathcal{T})$ as follows:
Here, the definition of $\mathrm{\Delta}{\mathscr{H}}_{\text{cyto}}$ is such that the likelihood of cell protrusions is enhanced if the concentration of cytoskeletal structures at the source grid site, $\u03f5({\mathbf{\mathbf{x}}}_{s})$, is larger than the concentration at the target grid site, $\u03f5({\mathbf{\mathbf{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 $\u03f5+\phi $. 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 $\alpha $ annexes $\mathbf{x}}_{t$ starting from $\mathbf{x}}_{s$, two things need to happen: First, focal adhesions formed by the target cell $\beta $ must be broken, implying a contribution $\mathrm{\Delta}{\mathscr{H}}_{\text{cyto}}=\u03f5({\mathbf{\mathbf{x}}}_{t})+\phi ({\mathbf{\mathbf{x}}}_{t})$ (fourth row of Equation S14a). Secondly, new focal adhesions are formed by the source cell $\alpha $, implying a contribution $\mathrm{\Delta}{\mathscr{H}}_{\text{cyto}}=\left[\u03f5({\mathbf{\mathbf{x}}}_{s})+\phi ({\mathbf{\mathbf{x}}}_{t})\right]$ (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(\mathcal{T})$ due to the polarization energy $\mathrm{\Delta}{\mathscr{H}}_{\text{cyto}}(\mathcal{T})$ is then defined by
where we set the effective thermal energy to ${k}_{\text{B}}T\equiv 1$. The characteristic ‘energy scale’ for $\mathrm{\Delta}{\mathscr{H}}_{\text{cyto}}$ is set by the polarization bounds ${\u03f5}_{0}\mathrm{\Delta}\u03f5/2$ and ${\u03f5}_{0}+\mathrm{\Delta}\u03f5/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 cellcell 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}_{\alpha ,\beta}$ and ${B}_{\alpha ,\beta}^{\prime}$ quantifying the system’s change in ‘energy’ upon forming a new contact between cells $\alpha $ and $\beta $ [${B}_{\alpha ,\beta}$] and upon breaking a preexisting contact between those cells by an ‘intruder cell’ $\gamma \ne \alpha ,\beta $ [${B}_{\alpha ,\beta}^{\prime}$]. In our computational model, we assume that formation of new cellcell contacts is energetically favored, and that breaking of preexisting contacts by intruder cells is energetically penalized. The matrix entries of ${B}_{\alpha ,\beta}$ and ${B}_{\alpha ,\beta}^{\prime}$, therefore, have a definite sign, which we take to be positive. The ordering of the cell index pair of ${B}_{\alpha ,\beta}$ and ${B}_{\alpha ,\beta}^{\prime}$ is of no physical significance, i.e. the adhesion matrices are symmetric. We also assume that a given cell $\alpha $ 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}_{\alpha ,\beta}$ and ${B}_{\alpha ,\beta}^{\prime}$ 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 cellcell adhesion’.
To implement the effects of cellcell adhesion, we compute the ‘energy difference’ $\mathrm{\Delta}{\mathscr{H}}_{\text{adh}}(\mathcal{T})$ for any given elementary event $\mathcal{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 $\alpha $ attempts a protrusion event ${\mathcal{T}}_{\text{pro}}$, involving the source grid site ${\mathbf{\mathbf{x}}}_{s}\in {\mathcal{B}}^{(\alpha )}$ and the target grid site ${\mathbf{\mathbf{x}}}_{t}\in {\mathcal{B}}^{(\beta )}$, as illustrated in Appendix 1—figure 3A. In this case, cell $\alpha $ acts as intruder cell, since the depicted protrusion event affects three preexisting contacts between the target cells $\beta $ and a ‘third party’ cell $\gamma $. Acceptance of the depicted protrusion event would have the following energetically relevant effects: (i) All preexisting contacts between the target cell $\beta $ and third party cell $\gamma \ne \alpha ,\beta $ at the target grid site ${\mathbf{\mathbf{x}}}_{t}$ are torn apart. (ii) New contacts between the source cell $\alpha $ and third party cell $\gamma \ne \alpha ,\beta $ are established. (iii) The length of the cell contact line between source cell $\alpha $ and target cell $\beta $ is changed. Altogether, these three effects lead to the following cell adhesion energy difference,
where $\mathrm{\ell}$ 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 $\alpha $ and target cell $\beta $ [points (ii) and (iii)]. The second term measures the (energetically penalized) breaking of preexisting cell contacts [point (i)] and, therefore, impedes cell $\alpha $’s ability to intrude. Conversely, if source cell $\alpha $ attempts a retraction event ${\mathcal{T}}_{\text{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 ${\mathbf{\mathbf{x}}}_{t}$ to the source site ${\mathbf{\mathbf{x}}}_{s}$; cf. Appendix 1—figure 3A:
where $\mathrm{\ell}$ is the length of a hexagon edge.
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 $\mathcal{\mathcal{T}}}_{\text{pro}$, which is immediately followed by its inverse retraction event $\mathcal{\mathcal{T}}}_{\text{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:
and can therefore conclude that
where $\mathcal{\mathcal{N}}}_{t$ denotes the neighborhood of the grid site which temporarily changes its cell index, and where $\alpha$ and $\beta$ are the indices of the source and target cells involved in the cyclic process; cf. Appendix 1—figure 3A. Since $(\mathrm{\Delta}B{)}_{\alpha ,\beta}\ge 0$, the above adhesive energy contribution is nonnegative, thus leading to an amount of energy equal to $\mathrm{\Delta}{\mathcal{\mathscr{H}}}_{\text{adh}}^{(\mathrm{c}\mathrm{y}\mathrm{c}\mathrm{l})}$ being dissipated as the cyclic process completes. This leads us to refer to the parameter matrix $(\mathrm{\Delta}B{)}_{\alpha ,\beta}$ 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 $\alpha$, this transformation entails to the following energetic effects (cf. discussion above): (i) Two of the preexisting cellcell contacts between the source cell’s upper neighbor in the left column (target cell $\beta$) and the corresponding cell in the right column (third party cell $\gamma$) get torn apart, leading to an energetic contribution $2{B}_{\beta ,\gamma}^{\prime}$. (ii) In return, two new contacts between the protruding (source) cell $\alpha$ and cell $\gamma$ are being established, leading to a contribution $2{B}_{\alpha ,\gamma}$ . (iii) Since the length of the contact line between cells in the left column (i.e. between protruding source cell $\alpha$ and retracting target cell $\beta$) 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}_{\alpha ,\beta}\equiv B$ and $B}_{\alpha ,\beta}^{\prime}\equiv {B}^{\prime$ ($\alpha \ne \beta$), and therefore, find
i.e. a nonnegative dissipative contribution per cell. The size of the dissipation parameter $\mathrm{\Delta}B$ 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(\mathcal{\mathcal{T}})$ as follows:
where we set the effective thermal energy to ${k}_{\text{B}}T\equiv 1$.
Rupture of cell contacts
By now, we have introduced all components making up the total acceptance probability $p(\mathcal{T})$, Equation S9. To conclude our discussions concerning the implementation of cellular traits, we highlight one additional aspect of elementary events. So far, the notion of an elementary event can be summarized as follows: Once source and target grid sites, ${\mathbf{\mathbf{x}}}_{s}$ and ${\mathbf{\mathbf{x}}}_{t}$, have been selected, acceptance of a protrusion [retraction] event causes (among other things like the distribution of regulatory factors) the cell index to be copied from ${\mathbf{\mathbf{x}}}_{s}$ to ${\mathbf{\mathbf{x}}}_{t}$ [from ${\mathbf{\mathbf{x}}}_{t}$ to ${\mathbf{\mathbf{x}}}_{s}$]. In other words, the domain ${\mathcal{D}}^{(\alpha )}$ of source cell $\alpha $ annexes ${\mathbf{\mathbf{x}}}_{t}$ [loses ${\mathbf{\mathbf{x}}}_{s}$], while the domain ${\mathcal{D}}^{(\beta )}$ of target cell $\beta $ is forced to let go ${\mathbf{\mathbf{x}}}_{t}$ [accommodate ${\mathbf{\mathbf{x}}}_{s}$]. However, if both source and target cells are actual cells, i.e. $\alpha ,\beta \ge 0$, and if the source cell attempts a retraction event, there is one additional possible outcome: If cell cohesion is weak, then the pulling force exerted by the retracting source cell $\alpha $ on its neighboring cells might also result in rupture of all preexisting contacts between the retracting source cell and its neighboring cells at ${\mathbf{\mathbf{x}}}_{s}$, rather than forcing one of its neighboring cells (the target cell) to fill the void created at ${\mathbf{\mathbf{x}}}_{s}$ once $\alpha $ retracts; cf. rupture event depicted in Appendix 1—figure 2. 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 $\alpha $ on target cell $\beta $ is strong enough to force $\beta $ to fill the void created at ${\mathbf{\mathbf{x}}}_{s}$ (i.e. to accommodate ${\mathbf{\mathbf{x}}}_{s}$), and call this a regular retraction event ${\mathcal{T}}_{\text{ret}}$. Secondly, we assume that the retraction of source cell $\alpha $ causes all preexisting cellcell contacts of cell $\alpha $ at ${\mathbf{\mathbf{x}}}_{s}$ to rupture, leaving a free substrate site at ${\mathbf{\mathbf{x}}}_{s}$ after the retraction event has been accepted. This latter event will be referred to as rupture event ${\mathcal{T}}_{\text{rup}}$. We then compute the total energy differences
and
under both assumptions (the energy difference associated with accepting ${\mathcal{T}}_{\text{rup}}$ can be computed with Equation S16b by using the substrate $\beta =1$ as new target cell) and compare the respective outcomes. If the rupture event is energetically favored over the regular retraction event, i.e. $\mathrm{\Delta}\mathscr{H}({\mathcal{T}}_{\text{rup}})<\mathrm{\Delta}\mathscr{H}({\mathcal{T}}_{\text{ret}})$, then cohesion between cells is weak. In this case, the rupture event ${\mathcal{T}}_{\text{rup}}$, rather than the regular retraction event ${\mathcal{T}}_{\text{ret}}$ will be attempted. Otherwise, cohesion is strong and a regular retraction event ${\mathcal{T}}_{\text{ret}}$ will be attempted.
Rupture of substrate contacts
In our discussion so far, a cyclic process that follows up a protrusion event $\mathcal{\mathcal{T}}}_{\text{pro}$ with its inverse retraction event $\mathcal{\mathcal{T}}}_{\text{ret}$, involving a cell $\alpha$ and no third party cells (in other words: no cellcell contacts are made or broken), will not yield a net energy cost or gain; cf. Equation S14a. To account for the dissipative nature of cellsubstrate contacts, we proceed similarly as when we have considered the disspative nature of cellcell contacts. We introduce dissipation in substrate adhesion by leaving the Hamiltonian unaltered for protrusion events but adding a penalty for retraction events:
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 $\mathrm{\Delta}\u03f5$ as the maximal energy gained by a single gainandloss in the absence of dissipation. In the presence of dissipation however, the cell has to pay at least a cost of $({\u03f5}_{0}\mathrm{\Delta}\u03f5/2)+D$ to detach at an arbitrary location, resulting in $\mathrm{\Delta}\u03f5D$ as the maximal energy gained by a single gainandloss 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=\mathrm{\Delta}\u03f5$ completely inhibit (directed) cell migration. Therefore, we study substrate dissipation in the range $D\in [0,\mathrm{\Delta}\u03f5]$.
Cell domain update routine
Having discussed the implementation concerning the two basic types of elementary events, namely protrusion events ${\mathcal{T}}_{\text{pro}}$ and retraction events ${\mathcal{T}}_{\text{ret}}$, as well as the two subtypes of regular retraction events and rupture events ${\mathcal{T}}_{\text{rup}}$, we can now summarize the cell domain update routine, point 3.5 in section ‘MonteCarlo scheme’. To this end, and in accordance with our previous notation, we use the cell indices $\alpha $ and $\beta $ to denote source and target cell, and ${\mathbf{\mathbf{x}}}_{s}$ and ${\mathbf{\mathbf{x}}}_{t}$ to denote source and target grid site. Moreover, equal signs "$=$" in the following listing are to be interpreted as assignment operators, where the value of the variable on the right hand side of the operator is assigned to the variable on the left hand side. With these preliminary remarks in mind, the cell update routine can be summarized as follows:
If the accepted elementary event is a protrusion event:
Set $\u03f5({\mathbf{\mathbf{x}}}_{t})=\u03f5({\mathbf{\mathbf{x}}}_{s})$ and $F({\mathbf{\mathbf{x}}}_{t})=F({\mathbf{\mathbf{x}}}_{s})$.
${\mathcal{D}}^{(\alpha )}\to {\mathcal{D}}^{(\alpha )}\cup \{{\mathbf{\mathbf{x}}}_{t}\}$.
${\mathcal{D}}^{(\beta )}\to {\mathcal{D}}^{(\beta )}\setminus \{{\mathbf{\mathbf{x}}}_{t}\}$.
Distribute regulatory factors according to Equation S12a.
If the accepted elementary event is a regular retraction event:
Set $\u03f5({\mathbf{\mathbf{x}}}_{s})=\u03f5({\mathbf{\mathbf{x}}}_{t})$ and $F({\mathbf{\mathbf{x}}}_{s})=F({\mathbf{\mathbf{x}}}_{t})$.
Set ${\mathcal{D}}^{(\alpha )}\to {\mathcal{D}}^{(\alpha )}\setminus \{{\mathbf{\mathbf{x}}}_{s}\}$
Set ${\mathcal{D}}^{(\beta )}\to {\mathcal{D}}^{(\beta )}\cup \{{\mathbf{\mathbf{x}}}_{s}\}$
Distribute regulatory factors according to Equation S12b.
If the accepted elementary event is a rupture event:
Set $\u03f5({\mathbf{\mathbf{x}}}_{s})=0$ and $F({\mathbf{\mathbf{x}}}_{s})=0$.
Set ${\mathcal{D}}^{(\alpha )}\to {\mathcal{D}}^{(\alpha )}\setminus \{{\mathbf{\mathbf{x}}}_{s}\}$
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 ${\kappa}_{A}^{(\alpha )}$ or an increase of the average cell polarization $\u03f5}_{0}^{(\alpha )$. We here employ the first method. We assume that individual cells grow exponentially (Barber et al., 2017) over a welldefined 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 $\sqrt{2}$. Therefore, as we do not scale the polarization field $\u03f5$ 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 $\sqrt{2}$. 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 ${\kappa}_{A}$ and ${\kappa}_{P}$ are, therefore, kept constant for a quiescent cell; we denote the corresponding values as $\kappa}_{A/0$ and $\kappa}_{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}_{\text{T}}=r\phantom{\rule{thinmathspace}{0ex}}{A}_{\text{ref}$, where the factor $r=\mathcal{\mathcal{O}}(1)$ relates the threshold size to the equilibrium cell size ${A}_{\text{ref}}$ reached by a free, solitary cell with constant polarization field $\u03f5={\u03f5}_{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 ${\u03f5}_{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 reenter 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:
Assign a state variable ${s}^{(\alpha )}$ to each cell which encodes the current phase in the cell cycle:
(S22) $s}^{(\alpha )}=\{\begin{array}{l}0,\text{quiescent phase}\\ 1,\text{refractory phase}\\ 2,\text{interphase}\\ 3,\text{mitotic phase}\end{array$Compute the equilibrium size ${A}_{\text{ref}}=({\u03f5}_{0}2\pi \sqrt{3}{\kappa}_{P})/(\sqrt{3}{\kappa}_{A})$ of a free, solitary cell with fixed polarization field $\u03f5={\u03f5}_{0}$, which spreads on the substrate used in the simulation.
At the beginning of the simulation, $t=0$, all cells are in the quiescent state, ${s}^{(\alpha )}(0)=0$, and have the following area and perimeter coupling constants, respectively: ${\kappa}_{A}^{(\alpha )}(0)={\kappa}_{A/0}$ and ${\kappa}_{P}^{(\alpha )}(0)={\kappa}_{P/0}$.
After the completion of each Monte Carlo time step $t$, perform one of the following changes for each cell:
Switch from quiescent to refractory state:
(S23) $\displaystyle \begin{array}{rl}& {\displaystyle {s}^{(\alpha )}(t)=0\text{}\wedge \text{}{A}^{(\alpha )}(t){A}_{\text{T}}}\\ & {\displaystyle \Rightarrow {s}^{(\alpha )}(t+1)=1.}\end{array}$Switch from refractory state to growing state with probability $p=1\mathrm{exp}(1/{T}_{r})$:
(S24) $\displaystyle \begin{array}{rl}& {\displaystyle {s}^{(\alpha )}(t)=1}\\ & {\displaystyle \Rightarrow {s}^{(\alpha )}(t+1)=\{\begin{array}{ll}2,& (p),\\ 1,& (1p),\end{array}}\end{array}$
where the terms in the brackets denote the respective probability.
Exponential growth in interphase over a period of ${T}_{g}$:
(S25) $\displaystyle \begin{array}{rl}& {\displaystyle {s}^{(\alpha )}(t)=2}\\ & {\displaystyle \Rightarrow {\kappa}_{A}^{(\alpha )}(t+1)={\kappa}_{A}^{(\alpha )}(t)\cdot (1/2{)}^{1/{T}_{g}}}\\ & {\displaystyle \Rightarrow {\kappa}_{P}^{(\alpha )}(t+1)={\kappa}_{P}^{(\alpha )}(t)\cdot (1/2{)}^{1/(2{T}_{g})}.}\end{array}$Switch from interphase to mitosis:
(S26) $\displaystyle \begin{array}{rl}& {\displaystyle {s}^{(\alpha )}(\tau )=2\text{ for all}\tau \in [t{T}_{g},\text{}t]}\\ & {\displaystyle \Rightarrow {s}^{(\alpha )}(t+1)=3.}\end{array}$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:
(S27) $\displaystyle \begin{array}{rl}& {\displaystyle {s}^{(\alpha )}(\tau )=3\text{ for all}\tau \in [t{T}_{d},\text{}t]}\\ & {\displaystyle \Rightarrow \text{ divide cell}\alpha \text{ into cells}(\alpha ,\beta )}\\ & {\displaystyle \Rightarrow {\kappa}_{A}^{(\alpha ,\beta )}(t+1)={\kappa}_{A/0}}\\ & {\displaystyle \Rightarrow {\kappa}_{P}^{(\alpha ,\beta )}(t+1)={\kappa}_{P/0}}\\ & {\displaystyle \Rightarrow {s}^{(\alpha ,\beta )}(t+1)=0.}\end{array}$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
which is a discretized version of the surface integral in Landau et al. (1986). Here, $\mathbf{f}}_{k$ is the force acting on the membrane element ${\mathbf{\mathbf{x}}}_{k}$ of cell $\alpha $, $\stackrel{~}{\mathbf{x}}}_{k}={\mathbf{x}}_{k}{\mathbf{x}}_{\text{com}}^{(\alpha )$ is the position of the element with respect to the center of mass $\mathbf{x}}_{\text{com}}^{(\alpha )$ of the cell, and the superscripts $i$ and $j$ are Cartesian indices. The forces $\mathbf{f}}_{k$ can be computed from the energy differences of all possible protrusion and retraction events originating from ${\mathbf{\mathbf{x}}}_{k}$,
where the number sign indicates a sum over substrate grid sites only, i.e. grid sites with $c({\mathbf{\mathbf{x}}}_{l})<0$, and where $\mathrm{\Delta}\mathcal{\mathscr{H}}\equiv {\mathcal{\mathscr{H}}}_{\text{cont}}+{\mathcal{\mathscr{H}}}_{\text{adh}}+{\mathcal{\mathscr{H}}}_{\text{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 ${\mathcal{D}}^{(\alpha )}$:
where $\stackrel{~}{\mathbf{x}}}_{k}={\mathbf{x}}_{k}{\mathbf{x}}_{\text{com}}^{(\alpha )$ denotes the coordinates of element $\mathbf{x}}_{k$ of cell $\alpha$, relative to the cell's center of mass $\mathbf{x}}_{\text{com}}^{(\alpha )$; 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}_{}$.
Appendix 2
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 ${\mathcal{R}}_{1}$, ${\mathcal{R}}_{2}$ and ${\mathcal{R}}_{3}$ within confinements of varying size and constant cell density. In the ${\mathcal{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 ${\mathcal{R}}_{2}$ and ${\mathcal{R}}_{3}$phases, we find scalefree 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 cellcell dissipation $\mathrm{\Delta}B$ impairs monolayer growth, while at the same time increasing front roughness. Similarly, an increased cellsubstrate dissipation $D$ also impairs monolayer growth. In contrast, increasing the maximum cell polarity $\mathrm{\Delta}\u03f5$ 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 cellfree 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 ${\kappa}_{A}$ for different values of perimeter stiffness ${\kappa}_{P}$, while keeping all other parameters constant. We find that the area occupied by the motile cell increases linearly with $1/{\kappa}_{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 $\u03f5={\u03f5}_{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 1BD). To intuitively explain this phenomenology, let us compare a cell of size ${A}_{\text{ref}}$ with a cell of size $r\phantom{\rule{thinmathspace}{0ex}}{A}_{\text{ref}}$, where $r\in [1,\text{}2]$, with the respective area stiffnesses $\kappa}_{A$ and ${\kappa}_{A}/r$. While the smaller cell has a perimeter ${P}_{\text{ref}}$, neglecting geometric effects the larger cell has a larger perimeter of approximately $\sqrt{r}\phantom{\rule{thinmathspace}{0ex}}{P}_{\text{ref}}$. Hence, the larger cell has to pay a larger energy cost (roughly by a factor $\sqrt{r}$) to ruffle its membrane by adding segments and therefore increasing its perimeter. Meanwhile, both the energy gain from the polarization field $\u03f5$ and the energy cost for increasing the cell area (as $A\phantom{\rule{thinmathspace}{0ex}}{\kappa}_{A}=\text{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 $\kappa}_{P}/\sqrt{r$, such that $P\phantom{\rule{thinmathspace}{0ex}}{\kappa}_{P}\approx \text{cst}$. We would then predict that the ratio $\kappa}_{P}^{2}/{\kappa}_{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 $\kappa}_{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 $\kappa}_{P}^{2}/{\kappa}_{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(\mathrm{\Delta}\u03f5\sqrt{{\kappa}_{A}}/{\kappa}_{P})\cdot g(R\sqrt{{\kappa}_{A}})$.
Single cell shape and dynamics depend on substrate dissipation
We have also studied the effect of cellsubstrate 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 $\mathrm{\Delta}\u03f5$ and cell perimeter stiffness $\kappa}_{P$; however, we were not able to achieve a data collapse in $D$ (Figure 2—figure supplement 1). We observe that with increasing cellsubstrate 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 cellsubstrate 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 Figure 2—figure supplement 1AC and EG (Figure 2—figure supplement 1D and H, respectively). Here and in the main text, we find that only cells with an aspect ratio larger than $2$ are motile (Figure 2F, Figure 2—figure supplement 1D, H; white regions).
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 $\mathcal{\mathcal{R}}}_{1$, $\mathcal{\mathcal{R}}}_{2$ or $\mathcal{\mathcal{R}}}_{3$phase, respectively. We studied the mean angular velocity
averaged over the set $\mathcal{\mathcal{C}}=\{\alpha \text{}\text{is not substrate}\}$ of all cells in confinement, where $\stackrel{~}{\mathbf{v}}}_{\alpha}(t)={\mathbf{v}}_{\alpha}(t)\u27e8{\mathbf{v}}_{\alpha}(t){\u27e9}_{\mathcal{\mathcal{C}}$ and $\stackrel{~}{\mathbf{R}}}_{\alpha}(t)={\mathbf{R}}_{\alpha}(t)\u27e8{\mathbf{R}}_{\alpha}(t){\u27e9}_{\mathcal{\mathcal{C}}$ are the velocity and position of cell $\alpha$ relative to the cell cluster, respectively.
In the lowly polarizable $\mathcal{\mathcal{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 counterclockwise 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 $\mathcal{\mathcal{R}}}_{1$phase, we observe that in the highly polarizable $\mathcal{\mathcal{R}}}_{2$ and $\mathcal{\mathcal{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 $\u27e8\left\omega \right\u27e9$ on the population size $N$ can be fitted by a power law of the form $\u27e8\left\omega \right\u27e9\propto {N}^{1/2}\propto {r}_{0}^{1}$. 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 $\left{\mathbf{v}}_{\text{rot}}\right\approx 0.008$ ($\mathbf{v}}_{\text{rot}}\perp \stackrel{~}{\mathbf{R}$). Furthermore, in the $\mathcal{\mathcal{R}}}_{2$, and $\mathcal{\mathcal{R}}}_{3$phases we have only scarcely observed switching of the rotational direction of cell clusters; e.g. for $4$cell clusters in the $\mathcal{\mathcal{R}}}_{3$phase.
Interestingly, fluctuations in the angular velocity ($\sigma}_{\omega$) change in a highly nonmonotonic 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 $\mathcal{\mathcal{R}}}_{1$phase, see Figure 4—figure supplement 1A; $3$ or $10$ cells in the $\mathcal{\mathcal{R}}}_{2$phase, see Figure 4—figure supplement 2A; $3$ cells in the $\mathcal{\mathcal{R}}}_{3$phase, see Figure 4—figure supplement 3A). This can likely be attributed to frustration of the cells in the population center (Segerer et al., 2015); cf. $10$ cells in Figure 4—figure supplement 2C.
Velocity and roughness of spreading tissue
We have studied the velocity and roughness of spreading tissue, while varying cellcell dissipation $\mathrm{\Delta}B$, cellsubstrate dissipation $D$ and maximum cell polarity $\mathrm{\Delta}\u03f5$.
First, let us introduce the observables that we are interested in. Let $X}_{>/<$ be the sets of $x$coordinates 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 cellfree 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:=\mathrm{a}\mathrm{b}\mathrm{s}({X}_{>/<})$. Then, we define the average front position as ${x}_{\text{F}}:=\mathbb{E}(X)$ and the front roughness as ${\sigma}_{\text{F}}^{2}:=\mathrm{V}\mathrm{a}\mathrm{r}(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, $\mathrm{m}\mathrm{a}\mathrm{x}({x}_{\text{F}})$, as well as the maximal roughness of the front $\mathrm{m}\mathrm{a}\mathrm{x}({\sigma}_{\text{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}_{\text{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 5figure supplement B,D,F).
First, we investigated how the monolayer expansion and front roughness depend on cellsubstrate dissipation, $\mathrm{\Delta}B$ (Figure 5—figure supplement 1A, B). Our simulations show that the cell sheet expands slower with increasing cellcell dissipation $\mathrm{\Delta}B$ (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 cellcell dissipation $\mathrm{\Delta}B$ (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 singlecell 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 $\mathrm{\Delta}\u03f5$) 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 cellsubstrate dissipation. Indeed, with increasing cellsubstrate dissipation, the monolayer expands slower, until this effect appears to saturate at a threshold value ${D}^{\star}\approx 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 proliferationdominated 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 migrationdominated mode of expansion? To test this hypothesis, we have analyzed how the monolayer growth and front roughness depend on the maximum cell polarity $\mathrm{\Delta}\u03f5$. As predicted, monolayer growth increases with the maximum cell polarity $\mathrm{\Delta}\u03f5$ (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 $\mathrm{\Delta}\u03f5$.
Data availability
We have uploaded the source code used in the main part of our study as well as the one used in the appendix. Furthermore, we have provided the full list of parameters in the figure captions, as well as exemplary parameter files for all applicable figures.
References

Vertex models: from cell mechanics to tissue morphogenesisPhilosophical Transactions of the Royal Society B: Biological Sciences 372:20150520.https://doi.org/10.1098/rstb.2015.0520

Propagating stress waves during epithelial expansionPhysical Review Letters 114:228101.https://doi.org/10.1103/PhysRevLett.114.228101

Details matter: noise and model structure set the relationship between cell size and cell cycle timingFrontiers in Cell and Developmental Biology 5:92.https://doi.org/10.3389/fcell.2017.00092

Periodic migration in a physical model of cells on micropatternsPhysical Review Letters 111:158102.https://doi.org/10.1103/PhysRevLett.111.158102

Velocity alignment leads to high persistence in confined cellsPhysical Review E 89:062705.https://doi.org/10.1103/PhysRevE.89.062705

Guidance of collective cell migration by substrate geometryIntegrative Biology 5:1026–1035.https://doi.org/10.1039/c3ib40054a

Prespecification and plasticity: shifting mechanisms of cell migrationCurrent Opinion in Cell Biology 16:14–23.https://doi.org/10.1016/j.ceb.2003.11.001

Collective cell migration in Morphogenesis, regeneration and CancerNature Reviews Molecular Cell Biology 10:445–457.https://doi.org/10.1038/nrm2720

Interstitial leukocyte migration and immune functionNature Immunology 9:960–969.https://doi.org/10.1038/ni.f.212

Cell topology, geometry, and morphogenesis in proliferating epitheliaCurrent Topics in Developmental Biology 89:87–114.https://doi.org/10.1016/S00702153(09)890042

Simulation of the differential adhesion driven rearrangement of biological cellsPhysical Review E 47:2128–2154.https://doi.org/10.1103/PhysRevE.47.2128

Simulation of biological cell sorting using a twodimensional extended potts modelPhysical Review Letters 69:2013–2016.https://doi.org/10.1103/PhysRevLett.69.2013

Regulation of cadherinmediated adhesion in morphogenesisNature Reviews Molecular Cell Biology 6:622–634.https://doi.org/10.1038/nrm1699

Symmetrybreaking in mammalian cell cohort migration during tissue pattern formation: role of randomwalk persistenceCell Motility and the Cytoskeleton 61:201–213.https://doi.org/10.1002/cm.20077

Cell adhesion is regulated by CDK1 during the cell cycleThe Journal of Cell Biology 217:3203–3218.https://doi.org/10.1083/jcb.201802088

Collective cell migration: leadership, invasion and segregationJournal of the Royal Society Interface 9:3268–3278.https://doi.org/10.1098/rsif.2012.0448

Mechanical modes of 'amoeboid' cell migrationCurrent Opinion in Cell Biology 21:636–644.https://doi.org/10.1016/j.ceb.2009.05.003

Organizing moving groups during morphogenesisCurrent Opinion in Cell Biology 18:102–107.https://doi.org/10.1016/j.ceb.2005.12.001

Mechanotransduction at cadherinmediated adhesionsCurrent Opinion in Cell Biology 23:523–530.https://doi.org/10.1016/j.ceb.2011.08.003

Blebbistatin inhibits contraction and accelerates migration in mouse hepatic stellate cellsBritish Journal of Pharmacology 159:304–315.https://doi.org/10.1111/j.14765381.2009.00477.x

Collisions of deformable cells lead to collective migrationScientific Reports 5:9172.https://doi.org/10.1038/srep09172

Myosin light chain kinase regulates cell polarization independently of membrane tension or rho kinaseThe Journal of Cell Biology 209:275–288.https://doi.org/10.1083/jcb.201409001

Polarization and movement of keratocytes: a multiscale modelling approachBulletin of Mathematical Biology 68:1169–1211.https://doi.org/10.1007/s1153800691317

Flow and diffusion in channelguided cell migrationBiophysical Journal 107:1054–1064.https://doi.org/10.1016/j.bpj.2014.07.017

Predicting division plane position and orientationTrends in Cell Biology 22:193–200.https://doi.org/10.1016/j.tcb.2012.01.003

Mathematics of cell motility: have we got its number?Journal of Mathematical Biology 58:105–134.https://doi.org/10.1007/s0028500801822

Crawling and gliding: a computational model for ShapeDriven cell migrationPLOS Computational Biology 11:e1004280.https://doi.org/10.1371/journal.pcbi.1004280

Improving the realism of the cellular potts model in simulations of biological cellsPhysica A: Statistical Mechanics and Its Applications 329:451–458.https://doi.org/10.1016/S03784371(03)005740

Cell spreading and lamellipodial extension rate is regulated by membrane tensionThe Journal of Cell Biology 148:127–136.https://doi.org/10.1083/jcb.148.1.127

Rotating lamellipodium waves in polarizing cellsCommunications Physics 1:73.https://doi.org/10.1038/s4200501800757

Signaling networks linking integrins and rho family GTPasesTrends in Biochemical Sciences 25:388–391.https://doi.org/10.1016/S09680004(00)016054

Emergence and persistence of collective cell migration on small circular micropatternsPhysical Review Letters 114:228102.https://doi.org/10.1103/PhysRevLett.114.228102

Mechanical waves during tissue expansionNature Physics 8:628–634.https://doi.org/10.1038/nphys2355

Computational model for cell morphodynamicsPhysical Review Letters 105:108104.https://doi.org/10.1103/PhysRevLett.105.108104

Migration of individual microvessel EndothelialCells  StochasticModel and parameter measurementJournal of Cell Science 99 Pt 2:419–430.

Collective cell motion in endothelial monolayersPhysical Biology 7:046007.https://doi.org/10.1088/14783975/7/4/046007

Modeling collective cell migration in geometric confinementPhysical Biology 14:035001.https://doi.org/10.1088/14783975/aa6591

Physical forces during collective cell migrationNature Physics 5:426–430.https://doi.org/10.1038/nphys1269

Mechanical cellmatrix feedback explains pairwise and collective endothelial cell behavior in vitroPLOS Computational Biology 10:e1003774.https://doi.org/10.1371/journal.pcbi.1003774

Model for selfpolarization and motility of keratocyte fragmentsJournal of the Royal Society Interface 9:1084–1092.https://doi.org/10.1098/rsif.2011.0433
Article and author information
Author details
Funding
German Excellence Initiative (NanoSystems Initiative Munich (NIM))
 Erwin Frey
Deutsche Forschungsgemeinschaft (Collaborative Research Center (SFB) 1032 (project B02))
 Erwin Frey
Deutsche Forschungsgemeinschaft (Graduate School of Quantitative Biosciences Munich (QBM))
 Andriy Goychuk
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
FT, AG, MR and EF designed research, performed research, and wrote the paper. AG acknowledges support by a DFG fellowship through the Graduate School of Quantitative Biosciences Munich (QBM). EF acknowledges support by the German Excellence Initiative via the program ‘NanoSystems Initiative Munich’ (NIM) and by the Deutsche Forschungsgemeinschaft (DFG) via Collaborative Research Center (SFB) 1032 (project B02). We thank Felix Kempf, Felix Segerer, Sophia Schaffer and Joachim Rädler for stimulating discussions.
Copyright
© 2019, Thüroff et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 6,763
 views

 1,015
 downloads

 51
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Cell Biology
Classical Gproteincoupled receptor (GPCR) signaling takes place in response to extracellular stimuli and involves receptors and heterotrimeric G proteins located at the plasma membrane. It has recently been established that GPCR signaling can also take place from intracellular membrane compartments, including endosomes that contain internalized receptors and ligands. While the mechanisms of GPCR endocytosis are well understood, it is not clear how well internalized receptors are supplied with G proteins. To address this gap, we use gene editing, confocal microscopy, and bioluminescence resonance energy transfer to study the distribution and trafficking of endogenous G proteins. We show here that constitutive endocytosis is sufficient to supply newly internalized endocytic vesicles with 20–30% of the G protein density found at the plasma membrane. We find that G proteins are present on early, late, and recycling endosomes, are abundant on lysosomes, but are virtually undetectable on the endoplasmic reticulum, mitochondria, and the medialtrans Golgi apparatus. Receptor activation does not change heterotrimer abundance on endosomes. Our findings provide a subcellular map of endogenous G protein distribution, suggest that G proteins may be partially excluded from nascent endocytic vesicles, and are likely to have implications for GPCR signaling from endosomes and other intracellular compartments.

 Cell Biology
 Plant Biology
Plants distribute many nutrients to chloroplasts during leaf development and maturation. When leaves senesce or experience sugar starvation, the autophagy machinery degrades chloroplast proteins to facilitate efficient nutrient reuse. Here, we report on the intracellular dynamics of an autophagy pathway responsible for piecemeal degradation of chloroplast components. Through livecell monitoring of chloroplast morphology, we observed the formation of chloroplast budding structures in sugarstarved leaves. These buds were then released and incorporated into the vacuolar lumen as an autophagic cargo termed a Rubiscocontaining body. The budding structures did not accumulate in mutants of core autophagy machinery, suggesting that autophagosome creation is required for forming chloroplast buds. Simultaneous tracking of chloroplast morphology and autophagosome development revealed that the isolation membranes of autophagosomes interact closely with part of the chloroplast surface before forming chloroplast buds. Chloroplasts then protrude at the site associated with the isolation membranes, which divide synchronously with autophagosome maturation. This autophagyrelated division does not require DYNAMINRELATED PROTEIN 5B, which constitutes the division ring for chloroplast proliferation in growing leaves. An unidentified division machinery may thus fragment chloroplasts for degradation in coordination with the development of the chloroplastassociated isolation membrane.