Bridging the gap between single-cell migration and collective dynamics

  1. Florian Thüroff
  2. Andriy Goychuk
  3. Matthias Reiter
  4. Erwin Frey  Is a corresponding author
  1. Ludwig-Maximilians-Universität München, Germany

Abstract

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

Introduction

Cell movements range from uncoordinated ruffling of cell boundaries to the migration of single cells (Ridley et al., 2003) to the collective motions of cohesive cell groups (Friedl and Gilmour, 2009). Single-cell migration enables cells to move towards and between tissue compartments – a process that plays an important role in the inflammation-induced migration of leukocytes (Friedl and Weigelin, 2008). One can distinguish between amoeboid and mesenchymal migration, which are characterized by widely different cell morphologies and adhesive interactions with their respective environments (Friedl, 2004; Lämmermann and Sixt, 2009). Cells may also form cohesive clusters and mobilize as a collective (Trepat et al., 2009; Angelini et al., 2011; Doxzen et al., 2013; Deforet et al., 2014; Vedula et al., 2012; Marel et al., 2014). This last mode of cell migration is known to drive tissue remodelling during embryonic morphogenesis (Lecaudey and Gilmour, 2006) and wound repair (Poujade et al., 2007).

Despite this broad diversity of migration modes, there appears to be a general consensus that all require (to varying degrees) the following factors: (i) Cell polarization, cytoskeletal (re)organization, and force generation driven by the interplay between actin polymerization and contraction of acto-myosin networks. (ii) Cell-cell cohesion and coupling mediated by adherens-junction proteins which are coupled to the cytoskeleton. (iii) Guidance by chemical and physical signals. The basic functionalities implemented by these different factors confer on cells the ability to generate forces, adhere (differentially) to each other and to a substrate, and respond to mechanical and chemical signals. However, a fully mechanistic understanding of how these basic functionalities are integrated into single-cell migration and coordinated multicellular movement is still lacking.

Here, we present a computational model which enables us to study cell migration at various scales, and thus provides an integrative perspective on the basic cell functions that enable the emergence of collective cell migration. While a variety of very successful modeling approaches has been used to describe single-cell dynamics (Mogilner, 2009; Marée et al., 2006; Marée et al., 2012; Shao et al., 2010; Ziebert et al., 2012; Ziebert and Aranson, 2013; Camley et al., 2013; Albert and Schwarz, 2014; Dietrich et al., 2018; Goychuk et al., 2018) or the movements of extended tissues (Szabó et al., 2006; Szabó et al., 2010; Kabla, 2012; Sepúlveda et al., 2013; Basan et al., 2013; Banerjee et al., 2015; Alt et al., 2017; Tarle et al., 2017), these models are hard to reconcile with each other. Models that focus on single cells are typically difficult to extend to larger cell numbers, largely due to their computational complexity. On the other hand, approaches which are designed to capture the dynamics at the scale of entire tissues generally adopt a rather coarse-grained point of view, and are therefore difficult to transfer to single cells or small cell cohorts. At present there are two partly competing and partly complementary approaches to bridge the gap between single-cell migration and collective dynamics, namely phase-field models (Shao et al., 2010; Ziebert et al., 2012; Shao et al., 2012; Camley et al., 2014; Camley and Rappel, 2014; Löber et al., 2015), and cellular Potts models (CPMs) (Szabó et al., 2010; Kabla, 2012; Szabó and Merks, 2013; van Oers et al., 2014; Segerer et al., 2015; Niculescu et al., 2015; Albert and Schwarz, 2016; Rens and Merks, 2017) first introduced by Graner and Glazier (1992).

Box 1.

A simple description of complex cells?

Mammalian cells are made up of around 109 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 𝒪(104) cells), thus permitting investigations of collective dynamics at the scale of tissues. Our model reproduces the most pertinent features of cell migration even in the limiting case of solitary cells, and is compatible with a wealth of experimental evidence derived from both small cell groups and larger collectives made up of several thousand cells. Specifically, by studying the characteristics of single-cell trajectories and of small cell groups confined to circular territories, we demonstrate that persistency of movements is significantly affected by cell stiffness and cell polarizability. Moreover, we investigate the dynamics of tissues in the context of a typical wound-healing assay (Poujade et al., 2007; Trepat et al., 2009; Serra-Picamal et al., 2012), and show that the model exhibits the recurring mechanical waves observed experimentally (Serra-Picamal et al., 2012), a feature which we attribute to the coupling between cell-sheet expansion and cell-density-induced growth inhibition.

Computational model

Model geometry

We consider cells that adhere to a two-dimensional surface, spanned by the coordinates (x,y), through some contact area (Figure 1A). Membrane protrusions and retractions, which determine cell motion and shape (Pollard and Borisy, 2003; Lauffenburger and Horwitz, 1996), correspond to size and shape changes of the surface contact area. We assume that processes that take place at the cell boundary drive cell motion, and therefore disregard the cell body, which extends into the 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, , 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.

Illustration of the computational model with the pertinent simulation steps.

(A) Illustration of a small cell cohort that adheres to a surface ((x,y)-plane). The polarization field, ϵ, is defined on the contact surface with the adhesion plane. The magnitude of the polarization field, which is indicated by the colorbar in Figure (C), encodes the local strength of cell-substrate adhesions and emulates the local mass of force-generating (pushing) cytoskeletal structures. Cell-cell adhesions are indicated in red. (B) Cytoskeletal structures respond to external mechanical stimuli through reaction networks involving different feedback loops. We greatly simplify these complex processes into two prototypic feedback loops, which break detailed balance and drive cell migration, as follows. The polarization field induces membrane protrusions and inhibits retractions. In turn, protrusions increase the polarization field (positive feedback) and therefore the likelihood of further protrusive activity, while retractions decrease the polarization field (negative feedback). In the absence of mechanochemical signals, the polarization field approaches its rest state. (C) Zoom-in to a common boundary shared between the substrate contact areas of three cells (bounded by the red lines), each represented by a contiguous set of occupied grid sites (hexagons). Top left: The upper right corner of the lower left cell (source cell) initiates a protrusion event against a neighboring element in the cell to its right (target cell), as indicated by the arrow, in an attempt to displace it. The success of each such attempted elementary event depends on the balance between contractile forces, cytoskeletal forces, and cell adhesion. Top right: If the protrusion event is successful, then the levels of regulatory factors are increased (decreased) in integer steps, at all lattice sites inside the source (target) cell that lie within a radius R of the accepted protrusion event (as indicated by the plus and minus signs). Bottom right: During the course of one MCS, different levels of regulatory factors accumulate locally within each cell, with positive levels of regulatory factors (green plus signs) promoting a build-up of cytoskeletal structures, negative levels of regulatory factors (red minus signs) causing degradation of cytoskeletal structures, and neutral levels of regulatory factors (white zero signs) causing relaxation towards a resting state, as indicated in the lower left image. The color code indicates local levels of cytoskeletal structures, ϵ.

Coarse-grained cellular mechanics

As discussed above, the configuration of a cell at any given time t is associated with a substrate contact area A(t) and perimeter P(t). We assume that the membrane and cortex deformations of each cell are constrained by the elastic energy

(1) cont(t)=κAA2(t)+κPP2(t),

where κA and κP are cell-type-specific stiffness parameters, similar to the original implementation of the CPM (Graner and Glazier, 1992). If the cell does not form adhesions to the substrate, then membrane and cortex contractility will round up the cell body, thereby collapsing the substrate contact area into a contact point.

Gripping the surface through the cell cytoskeleton

Detachment of the cell from the substrate is counteracted by focal adhesions, where the cell cytoskeleton is connected to the underlying substrate by integrins. Cellular protrusions are driven by outward pushing forces generated by the assembly and disassembly of cytoskeletal structures (Pollard and Borisy, 2003; Mogilner, 2009). As a first approximation, we subsume all of these complex dynamic processes, like the formation/degradation of focal adhesions and the assembly/disassembly of cytoskeletal structures, into a single time-dependent and spatially resolved internal field for each cell, ϵ(𝐱,t). This polarization field emulates the mass of force-generating cytoskeletal structures in the associated hexagon, at position 𝐱, 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

(2) cyto(t)=xϵ(x,t).

The polarization field must vanish at positions that are not occupied by a cell. Therefore, a retraction is associated with an energy penalty due to the loss of a substrate adhesion. Consequently, a protrusion, where one source hexagon ‘conquers’ a nearby target hexagon, is associated with an energy gain due to an increase of the substrate contact area. Here, we assume that the newly incorporated hexagon has the same polarization field as its conqueror.

There are several biological factors that constrain the local density of actin filaments, myosin and focal adhesions, whose limited availability corresponds to an upper bound on the polarization field. Furthermore, we assume that there is some minimal attachment energy associated with adhesions that prevents the cells from detaching from the substrate, which implies a lower bound on the polarization field. This motivates to introduce cell-type-specific bounds for the polarization field: ϵ(x,t)[ϵ0Δϵ/2, ϵ0+Δϵ/2], where ϵ0 is the average polarization field and Δϵ is the maximum cell polarity.

Active self-regulation of the cytoskeleton

Assembly and disassembly of cytoskeletal structures are controlled by a myriad of accessory proteins (Lauffenburger and Horwitz, 1996; Ridley et al., 2003). These regulatory proteins form a reaction network involving different feedback mechanisms, which allow cytoskeletal structures to respond to external mechanical stimuli (Marée et al., 2006; Marée et al., 2012). Furthermore, cytoskeletal structures like integrins play a role in the spatiotemporal control of these regulatory proteins (Schwartz and Shattil, 2000). Here, we refrain from formulating a detailed reaction-diffusion model that accounts for the interactions between all of these contributing players. Instead, we assume that the internal chemistry of the cell will generically produce protein patterns, with a typical length scale R, which locally up- or down-regulate cellular cytoskeleton and focal adhesion (dis)assembly. Then, we greatly simplify these complex processes (Lauffenburger and Horwitz, 1996; Schwartz and Shattil, 2000; Ridley et al., 2003) into two prototypic feedback loops (Figure 1B,C):

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

  2. The polarization field also locally inhibits inward motion of the membrane, by emulating the local adhesion strength of the cell to the substrate. If a membrane retraction is successful, then the loss of substrate adhesions locally further increases cell contractility, leading to a negative feedback on the polarization field within a range R.

In the absence of regulatory signals, we assume that the polarization field decays to a fixed value, ϵϵ0, which corresponds to a resting state of the cell cytoskeleton and focal adhesions. For the sake of keeping our model as simple as possible, we assume that all protein patterns have the same range R, and that the regulation of the cell cytoskeleton and focal adhesions follows a single timescale that corresponds to an update rate μ. Because at heart, our model is only based on generic feedback loops with a certain signaling range R, we would argue that any model with similar feedback should, in general, lead to similar cell behavior. Indeed, mutually repressing feedback loops (Marée et al., 2006) and mutually activating feedback loops (Shao et al., 2010; Ziebert et al., 2012; Albert and Schwarz, 2016) are crucial recurring motifs among multiple cell migration studies. Notably, these theoretical approaches all recover comparable cell behavior even when the model setup seems quite different at first glance:

  1. Cell migration couples mechanochemically to a scalar field (Shao et al., 2010), if stresses in the cell are isotropic; this is analogous to the present study.

  2. Cell migration couples mechanochemically to a vector field (Marée et al., 2006; Ziebert et al., 2012), if stresses in the cell are anisotropic.

  3. Cell migration couples to a single polarity vector (Albert and Schwarz, 2016), if propulsive forces are distributed homogeneously throughout the cell. However, this simplification of the former two cases cannot account for the formation of multiple competing lamellopodia/pseudopods.

These different modeling approaches (of varying complexity) surprisingly yield a universal phenomenology. The puzzling similarity between these models suggests generic common features that determine cell shape and motility: mechanical constraints like cell elasticity and mechanochemical feedback mechanisms that break detailed balance, maintain cell polarity and drive cell motion.

Intercellular adhesion and friction

In addition to internal remodeling of the cytoskeleton, adhesion of cells to neighboring cells and to the substrate plays a key role in explaining migratory phenotypes (Mogilner, 2009; Friedl and Gilmour, 2009). From a mechanical point of view, the implications of cell adhesion are two-fold:

  1. Cell adhesion supports growth of cell-cell and cell-matrix contacts and may thus be described in terms of effective surface energies. In our computational model, cell-matrix contacts are readily accounted for by the polarization field, ϵ. In addition, we associate the formation of cell-cell adhesions with an energy benefit B, which we call cell-cell adhesion parameter.

  2. Once formed, adhesive bonds anchor the cell to the substrate and to neighboring cells. During cell migration, these anchoring points must continuously be broken up and reassembled (Webb et al., 2002; Gumbiner, 2005) and, hence, provide a constant source of energy dissipation. Therefore, we assume that the cost for rupturing an existing cell-cell adhesion, B+ΔB>B, exceeds the gain from forming a new cell-cell adhesion. Then, the dissipative nature of cell-cell adhesions is accounted for by the cell-cell friction parameter ΔB. Similarly, cell-matrix contacts can also provide a source of dissipation, which is further discussed in Appendix 2.

Environmental cues

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

Tissue growth by cell division

In the description so far, the cells are arrested in the cell cycle (mitostatic). To investigate the effect of cell proliferation on tissue dynamics, we introduce a simplified three-state model of cell division. Cells start off in a quiescent state, in which their properties remain constant over time. The cell sizes fluctuate around an average value determined by the cell properties and the local tissue pressure. Cell growth typically arrests at large cell densities, in a phenomenon coined contact inhibition of proliferation (Stoker and Rubin, 1967; Puliafito et al., 2012; Pavel et al., 2018). Since large cell densities correspond to a small spread area for each individual cell, this implies that cell growth is arrested below a critical threshold size (AT). 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 Tg, 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 (κA and κ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 Tg, cells switch to a deterministic division state of duration Td. 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 (Δϵ0) in order to utilize its cytoskeleton for the separation of the cellular material, leading to mitotic rounding. At the end of the division state, each dividing cell splits into two identical daughter cells, whose properties and parameters are identical to the mother cell’s initial values in the quiescent state. Finally, the daughter cells re-initialize migration from an unpolarized state. For a detailed and more technical description we refer the interested reader to Appendix 1.

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:

  1. The cytoskeletal update rate endows the cellular Potts model with a reference timescale and determines the temporal discretization. In this study, we have set μ=0.1.

  2. The average polarization field ϵ0 encodes the energy gain for creating new cell-substrate adhesions, while the area stiffness κA represents the energy cost for increasing the substrate contact area. Then, the number of hexagons occupied by the cell is proportional to the ratio ϵ0/κA. If we use a desired cell area as reference value, then the ratio ϵ0/κA controls the spatial discretization of the cell. To study the migration of single cells and small cell cohorts, we have set the average polarization field to ϵ0=225 and the area stiffness to κA=0.18.

  3. In cellular Potts models, which are Monte-Carlo simulations, the reference energy of fluctuations is determined by an effective temperature. In this study, we have set kBT1.

Furthermore, we used a large computational grid with 9104 sites and periodic boundary conditions to study the migration of single cells. This leaves three parameters that control cell motility in the absence of cell-substrate dissipation: cell polarizability Δϵ, cell contractility κ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 κP and maximum cell polarity Δϵ 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 v^(t)v(t)/v(t) (𝐯: cell velocity) and (geometrical) center of mass position 𝐑(t) during a total simulation time of Tsim=104 Monte-Carlo steps (MCS). For each set of parameters, we performed 100 statistically independent simulations, from which we computed the mean squared displacement, MSD(τ)[𝐑(t+τ)-𝐑(t)]2, and the normalized velocity auto-correlation function, C(τ)v^(t+τ)v^(t). Here, 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(τ)e-τ/τp (inset of Figure 2A). We determined the persistence time of directed migration, τ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, Δϵ/κP (Figure 2B,D,E). This data collapse suggests Δϵ/κP as a relevant parameter (while cell polarizability and contractility are degenerate parameters), which we will henceforth refer to as specific polarizability.

Figure 2 with 2 supplements see all
Cell shape and persistence of migration as a function of cell polarizability.

(A) Mean-squared displacement (MSD) for single-cell movements at different maximum cell polarity Δϵ (stiffness parameters κP= 0.060, κA= 0.18; average polarization field ϵ0= 225; signaling radius R= 5; cell-substrate dissipation D= 0; cell-substrate adhesion penalty φ= 0; cytoskeletal update rate μ=0.1; 100 independent simulations for each set of parameters). Single cells perform a persistent random walk, i.e. they move ballistically (MSDτ2) for ττp, and diffusively (MSDτ) for ττp. Inset: Normalized velocity auto-correlation function for the same parameters as in the main figure. (B) Persistence time of directed cell migration plotted as a function of maximum cell polarity Δϵ, and perimeter stiffness κP (area stiffness κA=0.18; average polarization field ϵ0=225; signaling radius R=5; cell-substrate dissipation D= 0; cell-substrate adhesion penalty φ= 0; cytoskeletal update rate μ=0.1100 independent simulations for each set of parameters). The persistence time of the random walk increases with increasing cytoskeletal polarity and decreasing perimeter elasticity. (C) Cytoskeletal polarity also controls cell shapes, with crescent cell shapes (long persistence times) being observed at large cytoskeletal polartities, and round cell shapes (short persistence times) at small cytoskeletal polarities. Color code: cell polarization; cf. color bar in Figure 1C. (D) Single cell speed plotted as a function of maximum cell polarity Δϵ, and perimeter stiffness κP. (E) Single cell aspect ratio plotted as a function of maximum cell polarity Δϵ, and perimeter stiffness κP. (F) Speed and persistence time of single cells are correlated with the cell aspect ratio.

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, Δϵ/κP500, below which cells remain immobile (Figure 2B,D,E; grey regions). Above this threshold, the persistence time of directed migration, speed and aspect ratio increase markedly with the specific polarizability (Figure 2B,D,E). In our model, the area and perimeter stiffnesses refer to global and homogeneous cell contractility, while the cell polarization field drives cell migration. As discussed in ‘Gripping the surface through the cell cytoskeleton’, the cell polarization field does not explicitly distinguish between a local extensibility (e.g. due to actin polymerization), a local contractility (due to myosin-induced contraction) of the cytoskeleton or spatially regulated cell-substrate adhesions. For example, if cell migration is driven by actin polymerization, then blebbistatin treatment will decrease the global cell contractility, which we predict to lead to more elongated cells that move faster and exhibit extended episodes of ballistic motion. Indeed, an increase of cell migration speed after blebbistatin treatment was observed for mouse hepatic stellate cells (Liu et al., 2010). Alternatively, cell migration could also be driven by myosin contractility, for example by pulling the cell forward or by locally detaching adhesions. Then, polarizability and contractility concomitantly depend on the ability of the cell to exert forces, which can be inhibited by blebbistatin treatment. If polarizability, Δϵ, and contractility, κP, are equally reduced by a blebbistatin-dependent prefactor, then the specific polarizability, Δϵ/κ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.

Figure 3 with 1 supplement see all
Migratory behavior of single cells as a function of the cell’s signaling radius R at different values for the maximal cytoskeletal polarity Δϵ.

(Stiffness parameters κP=0.060κA=0.18; average polarization field ϵ0=225; cell-substrate dissipation D=0; cell-substrate adhesion penalty φ=0; cytoskeletal update rate μ=0.1100 independent simulations for each set of parameters). (A) The persistence times of directed migration of single cells exhibit a pronounced maximum at an optimal signaling radius, which depends on cell polarizability. (B) The shapes of single cells exhibit a pronounced maximal elongation at an optimal signaling radius, which depends on cell polarizability. (C) The signaling radius critically determines the synchronicity of internal cytoskeletal remodeling processes. Small signaling radii frequently lead to transient formation of mutually independent lamellipodia at different positions around the cell perimeter, thereby interrupting persistent motion (reducing persistence times of directed migration). Large signaling radii lead to structurally stable front-rear polarization profiles across the entire cell body (long persistence times of directed migration). Color code: cell polarization; cf. color bar in Figure 1C. (D) The speed of single cells does not drop to zero even when their persistence time of directed migration vanishes. This indicates single cell rotations. (E) The inverse curvature of the cell trajectories as a function of the signaling radius. (F) Depending on whether a cell migrates along its long axis (top) or short axis (bottom), it has to move a different projected contour length. If each protrusion takes roughly the same amount of time, then migration along the long axis (top; cell has to move a smaller projected contour length) allows for greater cell speeds than migration along the short axis (bottom; cell has to move a larger projected contour length).

Interestingly, because of this universal dependence of all the mentioned quantities on the specific polarizability, our simulations also show that there is a strong correlation between cell shape (aspect ratio) and cell motility (speed and persistence time of directed migration); see Figure 2F. While highly persistent trajectories are observed for cells with ‘crescent’ shapes, more erratic cell motion is typically found for cells with more rounded outlines (Figure 2C). In other words, our computational model predicts that cells which are able to polarize their cytoskeletal structures more strongly will adopt crescent shapes and show a higher degree of persistent cell motion. It would be interesting to further test these predictions by using phenotypic variations in cell shapes like those reported in experiments with keratocytes (Keren et al., 2008); there, the authors also found a correlation between cell shape and speed.

Feedback range determines whether individual cells move persistently or rotate

Moreover, we investigated the influence of different signaling radii R (typical range in which signalling molecules diffuse and mediate feedback mechanisms during a single Monte-Carlo step) on the persistence of single-cell trajectories. Since R is the relevant parameter that controls the spatial organization of lamellipodium formation, its value should strongly affect the statistics of a cell’s trajectory (Figure 3A). Indeed, at small values of R, we observe that the spatial coherence of cytoskeletal rearrangements is low, which frequently results in the disruption of ballistic motion due to the formation of independent lamellipodia in spatially separate sectors of the cell boundary (Figure 3C, lower snapshot). In contrast, at larger values of R, we find that spatial coherence is restored, and the formation of one extended lamellipodium across the cell’s leading edge maintains a distinct front-rear axis of cell polarity (Figure 3C, upper snapshot). However, when the signaling radius is too large compared to the cell size, we find an inhibition of ballistic motion and rounding of the cells as signals originating from one cell edge begin to reach the opposing edge. This effect may also occur when cells in tissue become smaller due to an increase of cell density through proliferation or compression; in other words, this means that the cells become smaller than the typical length scale of the chemical patterns that control cell migration. Then, one would not expect these chemical patterns to form (Hubatsch et al., 2019). Therefore, depending on the cell polarizability (Δϵ), 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 time-lapse videos of a cell at high polarizability (Δϵ/κP=1000; cf. Figure 3—video 1, top row), which show a qualitative shift in cell behavior:

  • For small signaling radii, R=2, short polarization fronts ‘pull’ the cell behind them, allowing for transient polarization and quick but erratic movement along the long axis of the cell.

  • For intermediate signaling radii, R=6, broad and correlated polarization fronts emerge, and both the cell polarization and movement always orient themselves along the short axis of the cell.

  • For large signaling radii R=15, we observed circular motion of the cell; because of the large signaling radius, signals originating from the trailing edge affect the leading edge of the cell and vice versa. Due to this circular motion, the cell exhibits a non-zero speed and a vanishing persistence time of directed migration.

Therefore, we find that the cell can transiently polarize and migrate along its long axis for small signaling radii and for high polarizability. Furthermore, in a broad parameter regime, we find keratocyte-like motion and polarization along the short axis of the cell. Note that we do not consider the formation of stress fibers, which could lead to cell migration along the long axis in a broad parameter regime (Kassianidou et al., 2019). Such stress fibers could be modeled via a nematic field that represents the anisotropic part of the intracellular stress. Our counter-intuitive observation that cell migration along the long axis is faster than cell migration along the short axis can be explained as follows: If the cell migrates along its short axis, then it has to move a greater projected contour length than if it migrates along its long axis (Figure 3F). Considering that each protrusion takes roughly the same amount of time, migration along the long axis allows for greater cell speeds than migration along the short axis, because the cell has to spend less time to move a smaller projected contour length (Figure 3F).

To further characterize the single cell rotations that occur at large signaling radii, we determined the average curvature of the trajectories c=sv^(s)v^(s), where s is the contour length along the corresponding trajectory. Here, we averaged the tangent vector v^(s) over 10 Monte-Carlo steps to integrate out fluctuations that occur on short timescales (the internal dynamics of the cell has an intrinsic time scale of 10 Monte-Carlo steps due to our choice of the cytoskeletal update rate, μ=0.1). We find that the curvature of the trajectories has a pronounced minimum at large signaling radii (where the persistence time of directed migration vanishes), which indicates a transition from straight to circular trajectories (Figure 3E). Such a transition from persistent migration to single cell rotations was previously observed in experiments (Lou et al., 2015; Raynaud et al., 2016) and in theory (Reeves et al., 2018; Allen et al., 2018).

Cell clusters on circular micropatterns

To assess the transition to collective cell motion, we next studied the dynamics of small cell groups confined to circular micropatterns (Huang et al., 2005; Doxzen et al., 2013; Deforet et al., 2014; Segerer et al., 2015). We implemented these structures in silico by setting φ(𝐱)=0 inside a radius r0 and φ(𝐱)- 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 𝐱0 (Huang et al., 2005; Doxzen et al., 2013; Deforet et al., 2014; Segerer et al., 2015). The statistics of these states of rotational motion provide insight into the influence of cellular properties on the group’s ability to coordinate cell movements. To quantify collective rotations, we recorded the average signed angular velocity of the cell cluster ω(t)=e^zv~(t)×R~(t)/R~(t)2𝒞. Here, e^z is the out-of-plane unit vector, 𝒞 denotes an average with respect to the cell population, and v~(t)=v(t)v(t)𝒞 as well as R~=R(t)R(t)𝒞 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, |ω(t)|, and the average cell perimeter P(t)Pα(t)𝒞 were then used to characterize the statistics of collective cell rotation. For each specific choice of simulation parameters, we monitored |ω(t)| and P(t) for a set of 100 statistically independent systems, each of which was observed over Tsim=104 MCS. From these data, we then computed the mean overall rotation speed |ω|, its standard deviation σω, and the standard deviation of the cell perimeter, σP.

Figure 4 illustrates the characteristic properties of collective cell rotations in systems containing |𝒞|=4 cells endowed with varying maximum cell polarity Δϵ and varying cell contractility κ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 Δϵ/κ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 (Δϵ/κP450 in Figure 4A), the rotation speed |ω| (purple curves in Figure 4A) vanishes and the cells are immobile. In this regime, which we term the stagnation phase, or 𝒮-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 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 (σP) and blue (σω) 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).

Figure 4 with 6 supplements see all
Phases of collective motion.

(4-cell systems; confinement radius r0=30.6; area stiffness κA=0.18; average polarization field ϵ0=225; signaling radius R=5; cytoskeletal update rate μ=0.1; cell-cell adhesion B=0; cell-cell dissipation ΔB=12; cell-substrate dissipation D=0; cell-substrate adhesion penalty φ=0 (r<r0), φ (r>r0); 100 independent simulations for each set of parameters). (A) Characteristic observables of collective cell rotation at different values of the cell perimeter stiffness parameter κP: mean (|ω|) and standard deviation (σω) of the magnitude of the cell cluster's angular velocity, and the standard deviation of the cell perimeter (σP). The statistics of collective cell motion depends only on the ratio of maximum cell polarity, Δϵ, to cell contractility, κP (specific polarizability). (B) Representative angular trajectories and (C) cell shapes (color code represents cell polarization; cf. Figure 1C) for the different parameter regimes as described in the main text. The cellular dynamics in the different parameter regimes are shown in Figure 4—video 1Figure 4—video 2 and Figure 4—video 3.

At intermediate values of specific polarizability (2-phase), the cellular systems reach a regime of enduring rotational motion, where |ω| varies linearly with the local specific polarizability, and where σP and σω exhibit a rather broad minimum (Figure 4A). In this regime, a range of ‘optimal ratios’ of cytoskeletal to contractile forces sustains stable cell shapes, and sets the stage for the formation of extended lamellipodia and the establishment of permanent front-rear polarizations of cells. As a result, the cells' persistence times of directed migration become very large, rendering cellular rotations strictly unidirectional within the observed time window (Figure 4B). Finally, at large values of the specific polarizability (3-phase), the system’s dynamics is dominated by cytoskeletal forces and the rotational speed |ω| 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 σ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 σω; blue curve in Figure 4A).

Tissue-level dynamics

As an application of our computational model at the tissue level, we considered a setup in which an epithelial cell sheet expands into free space. As in recent experimental studies (Serra-Picamal et al., 2012; Sepúlveda et al., 2013; Trepat et al., 2009; Poujade et al., 2007), we confined cells laterally between two fixed boundaries, within which they proliferated until they reached confluence; in the 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 ϵ0=35. For each set of parameters, we performed and averaged 100 independent simulations.

Figure 5 with 3 supplements see all
Expansion of a confluent epithelial cell sheet after removal of boundaries positioned at x=±175 for two different parameter settings.

(Stiffness parameters κP=0.12κA=0.18; average polarization field ϵ0=35; signaling radius R=2; cytoskeletal update rate μ=0.1; cell-cell adhesion B=12; cell-cell dissipation ΔB=0; cell-substrate dissipation D=0; cell-substrate adhesion penalty φ=0100 independent simulations for each set of parameters). (A–D) Tissue expansion for a migration-dominated setup without explicit cell growth and mitosis. (3300-cell system; maximum cell polarity Δϵ=30). (E–H) Tissue expansion at low density and cell polarizability for a cell sheet comprised of dividing cells. (Initially a 2500-cell system; maximum cell polarity Δϵ=10; growth time Tg=180; division time Td=20; size threshold for cell growth AT=1Aref, where Aref is the size of a solitary cell in equilibrium). (A, E) Snapshots of the polarization field ϵ; cf. Figure 5—video 1 and Figure 5—video 2. (B, F) Kymographs showing the cell density averaged over the y-direction and (top) final snapshots of the cell density profiles. (C, G) Kymographs showing the component σxx of the stress tensor averaged over the y-direction and (top) final snapshots of the stress profiles. (D, H) Kymographs showing the component vx of the cell velocities averaged over the y-direction and (top) final snapshot of the velocity profiles.

We first investigated how a densely packed pre-grown tissue of mitostatic cells with high polarizability (large Δϵ) expands into cell-free space upon removal of the confining boundaries at the tissue’s lateral edges (Figure 5A). As the cells migrate into the cell-free space, we observe a strongly (spatially) heterogeneous decrease in the initially high and uniform cell density and mechanical pressure in the expanding monolayer (Figure 5B,C). This is quite distinct from the behavior of a homogeneous and ideally elastic thin sheet, which would simply show a homogeneous relaxation in density as it relaxes towards its rest state. Moreover, cell polarization and the ensuing active cell migration lead to inhomogeneously distributed traction stresses in the monolayer. After initial expansion of the monolayer, facilitated by high mechanical pressure, the cells at the monolayer edge begin to polarize outwards, which enhances outward front migration. These actively propagating cells exert traction on the trailing cells, and thereby yield a trailing region with negative stress (Figure 5C). Taken together, this gives rise to a characteristic X-shaped pattern in the kymograph of the total mechanical stresses σxxy (Figure 5C). This profile closely resembles the first period of mechanical waves observed experimentally (Serra-Picamal et al., 2012). It illustrates how stress is transferred towards the center of the monolayer when cells are highly motile and collectively contribute to tissue expansion. At the end of the simulated time window, the cell density exhibited a minimum in the center of the sheet (Figure 5B). This is due to stretching of the central group of cells caused by the equally strong traction forces exerted by their migrating neighbors on both sides. Finally, the simulations also show that outward cell velocities increase approximately linearly with the distance from the center, confirming that in this configuration the entire cell sheet contributes to the monolayer expansion (Figure 5D).

To explore the possible range of tissue dynamics and expansion, we also investigated a qualitatively different parameter regime where cells are less densely packed and can also polarize less due to a narrower range of polarizability (Figure 5E). Here, the expansion of the monolayer is mainly driven by cell division, and cells keep dividing until they reach a homeostatic cell density (Figure 5F). Even though cells should typically exceed the threshold size and hence enter the growth phase at different times, we observe that the cell sheet exhibits periodic ‘bursts’ of growth (Figure 5F) coinciding with the total duration of a complete cell cycle (200 MCS) and alternating with cell migration (Figure 5H). These periodic ‘bursts’ can be explained as follows. Initially, the slightly compressed monolayer expands to relieve mechanical pressure. Due to this initial motion, the cells at the monolayer edge begin to polarize outwards. As in the previous case, where cell proliferation is absent (Figure 5A–D), the polarized cells enhance outward front migration and stretch the cells in the bulk of the cell sheet. For the same reasons as before, we observe a typical X-shaped stress pattern in the kymograph (Figure 5G), albeit less pronounced due to the lower polarizability of the cells (cf. Figure 5C). Because a broad region in the monolayer bulk is stretched by the actively migrating cell fronts, these cells exceed the threshold size and begin growing approximately in phase. Once the mechanical pressure of the cell sheet is relieved, it will stop expanding (Figure 5H). However, cell growth and division once more lead to an increase in mechanical pressure (and cell density) in the monolayer (Figure 5F,G). This cycle of migration-dominated monolayer expansion and cell-density-dependent cell growth and division results in a periodic recurrence of the X-shaped stress pattern (Figure 5G), closely resembling the pattern observed in experiments (Serra-Picamal et al., 2012). On a sidenote, the synchronization of the cell division and cell migration phases by the deterministic portion of the cell cycle can be counteracted by introducing additional stochastic terms in the transition between the different phases of the cell cycle (cf. 'Cell proliferation and mitosis' in Appendix 1).

Note that the inhomogeneously distributed traction stresses in the monolayer, and its wave-like behavior, ultimately emerge from cell polarization and the ensuing active cell migration. Therefore, these traction patterns would look much less prominent if one were to inhibit cell motility (compare Figure 5C with Figure 5G).

Finally, we investigated which parameters control the roughness of the tissue fronts. We found that increasing cell motility, or increasing cell-cell dissipation leads to rougher front morphologies (Figure 5—figure supplement 1 and ‘Velocity and roughness of spreading tissue’ in Appendix 2). Therefore, we hypothesized that one could observe fingering of cell monolayers by adjusting the parameters accordingly:

  • Increase of cell motility by decreasing the membrane stiffness and at the same time increasing polarizability and signaling radius of the cells.

  • Increase of cell-cell dissipation and slight decrease of cell-cell adhesion.

  • Slower and less homogeneously distributed cell division by increasing the cell threshold size.

Indeed, we then observe a drastic roughening of the cell fronts and small cohorts of cells that coherently move into cell-free space (Figure 6). This roughening is more pronounced if we further increase the threshold size that a cell has to exceed to initiate growth (cf. Figure 6A,E). Analogously to our previous discussion, we observe that an increasing mechanical pressure in the monolayer due to the division of cells initiates outward cell migration (Figure 6B,F). Then, cells in the tissue begin to polarize outwards and coordinate their motion with their neighboring cells, leading to small coordinated cell cohorts. As before, we also find distinct traction force patterns, as recurring waves of high stress travel backwards relative to the leading edges (Figure 6C,G), and distinct recurring velocity patterns (Figure 6D,H).

Figure 6 with 2 supplements see all
Expansion of a confluent epithelial cell sheet after removal of boundaries positioned at x=±175 for two different parameter settings that produce rough tissue fronts.

(Initially a 2500-cell system; stiffness parameters κP=0.10, κA=0.18; average polarization field ϵ0=35; maximum cell polarity Δϵ=20; signaling radius R=5; cytoskeletal update rate μ=0.1; cell-cell adhesion B=5; cell-cell dissipation ΔB=10; cell-substrate dissipation D=0; cell-substrate adhesion penalty φ=0; growth time Tg=180; division time Td=20; 100 independent simulations for each set of parameters). (A–D) Tissue expansion at low density and cell polarizability for a cell sheet comprised of quickly dividing cells. (Size threshold for cell growth AT=1.05Aref, where Aref is the size of a solitary cell in equilibrium). (E–H) Tissue expansion at low density and cell polarizability for a cell sheet comprised of slowly dividing cells. (Size threshold for cell growth AT=1.10Aref, where Aref is the size of a solitary cell in equilibrium). (A, E) Snapshots of the polarization field ϵ; cf. Figure 6—video 1 and Figure 6—video 2. (B, F) Kymographs showing the cell density averaged over the y-direction and (top) final snapshots of the cell density profiles. (C, G) Kymographs showing the component σxx of the stress tensor averaged over the y-direction and (top) final snapshots of the stress profiles. (D, H) Kymographs showing the component vx of the cell velocities averaged over the y-direction and (top) final snapshot of the velocity profiles.

Discussion

In this work, we have proposed a generalization of the cellular Potts model (Graner and Glazier, 1992). The model implements a coarse-grained routine that captures the salient features of cytoskeletal remodeling processes on subcellular scales, while being computationally tractable enough to allow for the simulation of entire tissues containing up to 𝒪(104) cells. We have used the model to study the transition from single-cell to cohort cell migration in terms of the interplay between the pertinent cellular functions. Specifically, we have demonstrated that our model consistently reproduces the dynamics and morphology of motile cells down to the level of solitary cells. Our studies also reveal that cytoskeletal forces (relative to cell contractility), as well as the spatial organization of the cells’ lamellipodia, significantly affect the statistics of cellular trajectories, both in the context of single-cell motion and in cohesive cell groups restricted to circular micropatterns. On larger scales, our simulation results suggest that the dynamics of expanding tissues strongly depends on the specific properties of the constituent cells. If monolayer expansion is driven by active cell migration throughout the tissue, then the cell sheet exhibits typical traction-force patterns and an X-shape in the corresponding kymograph. Additionally, a cell-density-dependent cell growth leads to a periodic recurrence of these traction-force patterns in a cycle of migration-dominated expansion and ’burst’-like cell proliferation.

Taken together, our results further highlight the intricacies of collective cell migration, which involves a multitude of intra- and inter-cellular signaling mechanisms operating at different scales in length and time. Establishing a comprehensive picture that incorporates and elucidates the mechanistic basis of these phenomena remains a pressing and challenging task. The multiscale modeling approach proposed here provides a direct link between subcellular processes and macroscopic dynamic observables, and might thus offer a viable route towards this goal.

Materials and methods

The computational model is described in section ‘Computational model’. The numerical implementation of the model is discussed in detail in Appendix 1. The parameter files and source files associated with the figures are given in Table 1.

Table 1
Source and parameter files used for each figure.

All source and parameter files are found in Source data 1.

FigureSimulation codeProcessing codeParameters
Figure 2CPM_NoDivisionTrajectoryAnalysisSinglesingle_Q
Figure 2—figure supplement 1 (A-D)CPM_NoDivisionTrajectoryAnalysisSinglesingle_DQ
Figure 2—figure supplement 1 (E-H)CPM_NoDivisionTrajectoryAnalysisSinglesingle_DM
Figure 3CPM_NoDivisionTrajectoryAnalysisSinglesingle_R
Figure 4CPM_NoDivisionTrajectoryAnalysisCircularPatternrotation_Q
Figure 4—figure supplement 1CPM_NoDivisionTrajectoryAnalysisCircularPatternrotation_N_R1
Figure 4—figure supplement 2CPM_NoDivisionTrajectoryAnalysisCircularPatternrotation_N_R2
Figure 4—figure supplement 3CPM_NoDivisionTrajectoryAnalysisCircularPatternrotation_N_R3
Figure 5 (A-D)CPM_Divisionwound_nodiv
Figure 5 (E-H)CPM_Divisionwound_div
Figure 5—figure supplement 1 (A-B)CPM_Division_SupplementFrontAnalysiswound_div_A
Figure 5—figure supplement 1 (C-D)CPM_Division_SupplementFrontAnalysiswound_div_D
Figure 5—figure supplement 1 (E, F)CPM_Division_SupplementFrontAnalysiswound_div_Q
Figure 6 (A-D)CPM_Divisionwound_div_fing_1.0
Figure 6 (E-H)CPM_Divisionwound_div_fing_1.1
Appendix 2—figure 1CPM_NoDivisionTrajectoryAnalysisSinglesingle_A

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 model-specific terminology which will be used throughout this section to illustrate the mechanics of our model. 

Computational grid 

The basic data structure that underlies our computational model is referred to as the grid; see Appendix 1—figure 1. The grid itself is implemented as a regular, space-filling lattice with lattice vectors {xi}i=1,,N. Each lattice vector xi is understood to represent its associated Voronoi cell which will be referred to as grid site. To be specific, we consider triangular tilings {xi}i=1,,N, such that each grid site is a hexagon, which is surrounded by 6 nearest-neighbor sites that define the neighborhood 𝒩k of xk:

(S1) 𝒩k={xj | xjis\ nearest\ neighbor\ of\ xk}

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.

(S2) xkxj=1  j𝒩k.

This then implies for the side length and the two-dimensional volume (area) a of each hexagonal grid site: =1/3 and a=332/2.

Appendix 1—figure 1
Illustration of the various sets defining a cell and its environment.

Grid sites occupied by cell α, i.e. its domain 𝒟(α), are indicated in red colors. The cell’s membrane sites, (α), are indicated by the lighter red color, the cell’s immediate neighborhood, 𝒩(α), is indicated in gray. Elementary events involving cell α always involve one grid site in (α) and one grid site in 𝒩(α). For the hexagonal lattices used in this work, each grid site 𝐱k is surrounded by 6 nearest neighbors which we collectively denote by 𝒩k.

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

(S3) 𝒟(α)={xk | c(xk)=α},

where the indicator function c(𝐱k) gives the index of the cell occupying 𝐱k. Here and in the following, we use latin indices to reference lattice sites, and greek indices to reference cells. The set 𝒟(α) used to represent the spatial extension of cell α, will be referred to as the domain of cell α. In our model, each grid site 𝐱k can be occupied by at most one cell (i.e. we do not allow for overlapping cell domains). The absence of cells at 𝐱k is numerically implemented by negative values of the indicator function, c(𝐱k)<0. Following this terminology, the area and the perimeter of cell α are given by:

(S4a) Aα=ak=1Nδα,c(xk)=3322k=1Nδα,c(xk),
(S4b) Pα=k=1Nxl𝒩kδα,c(xk)(1δα,c(xl)).

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 𝒟(α) of cell α changes over time. The evolution of cell shape and position, as represented by 𝒟(α), 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 α (referred to as source cell) incorporates one grid site 𝐱t (referred to as target grid site) from its neighborhood 𝒩(α),

(S5) 𝒟old(α)𝒟new(α)=𝒟old(α){xt},xt𝒩(α),

thereby increasing its cellular domain by one grid site. Here, the neighborhood of cell α, 𝒩(α), is defined as

(S6) 𝒩(α)={xl | minxk𝒟(α)xlxk=1}.

During a retraction event, source cell α expels one of its membrane grid sites 𝐱s(α),

(S7) 𝒟old(α)𝒟new(α)=𝒟old(α){xs},xsold(α),

where the set of membrane grid sites (α) is defined as

(S8) (α)={xk𝒟(α) | minxl𝒩(α)xkxl=1}.

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 α into the domain of cell β (referred to as target cell) must be accompanied by a corresponding retraction event 𝒟old(β)𝒟new(β)=𝒟old(β){𝐱t}, where 𝐱t denotes the target grid site annexed by cell α. We emphasize, however, that the reverse is not generally true. If source cell α retracts, i.e. loses one of its boundary grid sites 𝐱s(α), the lost grid site 𝐱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 α exerts a pulling force on one of its neighboring cells β (the target cell) and forces the target cell to fill the emerging void at 𝐱s, i.e. 𝒟old(β)𝒟new(β)=𝒟old(β){𝐱s}, where 𝐱s denotes the grid site lost by cell α. On the other hand, if adhesion between cells is weak, then retraction of the source cell α can lead to a rupture of pre-existing cell contacts between α and other cells at 𝐱s, such that the lost grid site 𝐱s becomes free space [c(𝐱s)=α0c(𝐱s)<0]. Details on the actual implementation of cell rupture are discussed in section ‘Rupture of cell contacts’.

Monte-Carlo scheme

In the spirit of a standard Monte-Carlo scheme, the actual simulation proceeds via a succession of Monte-Carlo steps, where each Monte-Carlo step (MCS) propagates the state of the simulated cell population from time t to time t+Δt, where we set the time step to Δt1. 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 𝒟(α) as well as the numerical values of cell areas Aα and perimeters Pα 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 Monte-Carlo scheme:

  1. Initialize the cell population and define the duration of the simulation, i.e. the number of MCS, Tsim, to be performed.

  2. Set the simulation time t=0.

  3. Perform the next MCS; this step is further detailed below.

  4. Update polarization fields (cf. section ‘Cytoskeletal structures and focal adhesion’).

  5. Set t=t+Δt, where Δt1.

  6. Repeat steps 3–5 while t<Tsim.

The implementation of a MCS, i.e. the sequence of elementary events, is based on the following general considerations:

  1. Choice of source and target grid sites. Each elementary event 𝒯 originates from a membrane grid site 𝐱s(α) of some cell α, 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 𝐱s and which is not currently occupied by cell α: 𝐱t𝒩s𝒟(α). In what follows, this additional grid site 𝐱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 β, 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.

  2. Monte-Carlo method to generate the system’s dynamics. As mentioned above, the actual dynamics of cells in our computational model is driven by a succession of elementary events, whose cumulative effects over time allow cells to change shapes and to move relative to the substrate as well as relative to each other. Following a standard Monte-Carlo procedure, the probability of occurrence of elementary events 𝒯 is determined by a goal function p(𝒯) [cf. point (iii) below]. However, since elementary events come in two basic flavors, protrusions 𝒯pro and retractions 𝒯ret, their actual occurrence is controlled by a two-step 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 xt, or retracts from xs. Then, a decision is made with equal probabilities as to whether one attempts 𝒯pro or 𝒯ret. In a second step, the goal function p is used to compute the occurrence probability of the attempted event 𝒯. Finally, this elementary event 𝒯 is being accepted with probability p(𝒯).

  3. Choice of the goal function p(T). As has been detailed above, we use a goal function p(𝒯) to control the occurrence and acceptance of elementary events 𝒯. Following the standard cellular Potts model (Graner and Glazier, 1992; Glazier and Graner, 1993), this goal function takes into account the effects of cell contractility and cell-cell adhesion, using, however, a slightly different implementation; cf. sections ‘Cell contractility’ and ‘Cell adhesion’. In addition, we generalized the definition of the goal function p(𝒯) to explicitly take into account a simplified model of cytoskeletal structures and the ensuing polarization of cells. The actual definition of the goal function will be developed in section ‘Implementation of cellular traits’, where, moreover, details concerning the implementation of the cell polarization model will be discussed.

The implementation of a single MCS loop is then given by the following simulation scheme:

  1. Determine the current number of trials per Monte-Carlo step (MCS), K=αPα/, and set the trial counter n=0.

  2. With equal probability, choose a cell membrane segment (cf. solid black line in Appendix 1—figure 1) from a random cell α of the cell population. Because the cell membrane represents the border between lattice sites occupied by cell α and unoccupied by cell α, the chosen membrane segment automatically defines the source grid site xs(α) and the corresponding target grid site xt𝒩(α)𝒩s.

  3. With equal probability, choose whether to attempt a protrusion event (𝒯pro) or a retraction event (𝒯ret).

  4. Compute the prospective acceptance probability p(𝒯pro/ret) corresponding to the attempted event, and decide whether to accept the attempted event on the basis of this probability.

  5. If the attempted elementary event has been accepted, then update the cellular domains of source cell α and opponent cell β; for details see section ‘Cell domain update routine’.

  6. If n<K, set nn+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(𝒯) of an elementary event 𝒯. Specifically, our model takes into account cell contractility, the assembly and disassembly of cytoskeletal structures, cell-cell adhesion, and focal adhesions. We will assume that each of these cellular properties contributes independently to the acceptance probability p, such that

(S9) p=min{1, pcontpcytopadh}.

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 pcyto(𝒯). In the following sections, we give detailed discussions for each of these contributions, separately.

Cell contractility

In biological cells, membrane fluctuations are constrained by elastic forces and contractile cytoskeletal structures, which play a vital role in cell migration (Alberts et al., 2015; Raucher and Sheetz, 2000; Friedl, 2004). In our computational approach, we take cell contractility into account by assigning a contractile ‘energy’

(S10) cont=α[κP(α)Pα2+κA(α)Aα2],

with positive coupling constants κP(α) and κA(α) characterizing the contractility of cell α; for empty substrate sites (α<0) we set κP(α)=κA(α)=0. According to Equation S10, the cell’s ‘contractile energy’ increases with increasing cell perimeter and increasing cell area. The model Hamiltonian cont can then be used to specify the contractile contribution to the goal function p(𝒯). To this end, let Δcont(𝒯) denote the contractile contribution to the energy difference entailed by accepting an elementary event 𝒯. Following a standard Metropolis algorithm, we then define

(S11) pcont(𝒯):=exp[Δcont(𝒯)/kBT],

where we set the effective thermal energy to kBT1. 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 𝒪(104) can be achieved at acceptable computation times) and, at the same time, to capture the essential effects of cytoskeletal dynamics to attain meaningful results down to the level of single cells. Thus, instead of accounting for a detailed biochemical description by means of reaction-diffusion networks (Marée et al., 2006; Marée et al., 2012), we resort to a simplified implementation of the most pertinent features of cytoskeletal dynamics. Specifically, we propose a rule-based algorithm to model cytoskeletal structures and to assess the integrated effects of cell polarity, cell contractility and adhesion on the collective dynamics of cells as parts of larger groups.

To this end, we define a scalar field ϵ(xn), xn𝒟(α), on the domain of each cell α. The local quantity ϵ(xn) will be referred to as polarization field and is taken to be a measure for the density of cytoskeletal structures at position xn within the cell’s body. The field variable ϵ(xn) is dynamically updated as the simulation progresses, reflecting cytoskeletal remodeling. To set up a system of rules underlying the actual implementation of these cytoskeletal remodeling processes, we resort to the following biologically motivated premises:

  1. The scalar polarization field ϵ is bounded: The dynamics of cytoskeletal remodeling not only depends on the local number (density) of actin monomers and polymers, but also on a multitude of accessory proteins controlling cytoskeleton assembly and disassembly. Several biological factors—including the action of sequestering proteins like thymosin-β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: ϵ(xn,t)[ϵ0Δϵ/2, ϵ0+Δϵ/2]. These bounds are cell-type specific. While the upper bound ϵ0+Δϵ/2 mainly reflects the limited availability of protein resources, the lower bound ϵ0Δϵ/2 serves to prevent cells from collapsing.

  2. Regulatory proteins affect assembly and disassembly of cytoskeletal structures: The assembly and disassembly of cytoskeletal structures, numerically encoded by ϵ(xn), 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(𝐱n) for each grid site 𝐱n𝒟(α). We use F(𝐱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(𝐱n)>0, reflect local conditions in support of network-assembly, whereas negative levels, F(𝐱n)<0, represent predominantly degrading (or disassembly) conditions. For neutral levels, F(𝐱n)=0, the network gradually restores its rest state.

  3. Feedback between cytoskeletal structures and regulatory factors: The activities of accessory cytoskeletal proteins which regulate the local levels of cytoskeletal structures are themselves controlled by a number of mechanical and chemical signals received by the cell. Here and in the following, our focus will be on mechanical signals. For example, important regulatory proteins like the Arp2/3 complex are activated locally at the cell membrane, from where they diffuse into the bulk of the cell until they are bound by actin (Kovacs et al., 2002; Pollard and Borisy, 2003; Leckband et al., 2011). Adopting a coarse level of description, this diffusion-degradation dynamics entails a finite range of regulatory proteins, which are activated at the cell’s membrane. In our model, we use the integer variable F(xn) to implement this propagation of mechanical information, perceived by cell α at its periphery (α), across a certain spatial distance R. The local levels of F(𝐱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 xs(α) (source cell: α; target cell: β), then for all sites xn within a range R (i.e. xnxs<R) the integer variable signifying regulatory factors is incremented up and down for the protruding and the retracting cell, respectively:

(S12a) F(xn){F(xn)+1,xn𝒟(α),F(xn)1,xn𝒟(β).
  • Similarly, if a retraction event has been accepted at the source site xs(α) , and the (local) cell contact between source cell α and target cell β has remained intact, then within a range R one applies the inverse update rule:

(S12b) F(xn){F(xn)1,xn𝒟(α),F(xn)+1,xn𝒟(β).
  • If a retraction event has been accepted at the source site xs(α) , and in addition the (local) cell contact between source cell α and target cell β has ruptured, then the regulatory factors are reduced only within a range R in the retracting cell:

(S12c) F(xn){F(xn)1,xn𝒟(α),no\ update,else.

Finally, if the target grid site 𝐱t is not occupied by any cell (substrate is indicated by β<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 𝒟(α) as the current MCS progresses. At the end of each MCS, the accumulated (local) values of F(𝐱n) are used to update the local values of the polarization field ϵ(𝐱n) inside each cell α0 (𝐱n𝒟(α)): We assume that for positive values, F(xn)>0, there is assembly of cytoskeletal structures and ϵ is increased by an amount proportional to the distance of ϵ from its upper bound ϵ0+Δϵ/2:

(S13a) ϵ(xn,t+Δt)=ϵ(xn,t)+Δtμ[ϵ0+Δϵ/2ϵ(xn,t)],

where the time step is defined as Δt1. Thereby ϵ0+Δϵ/2 is a fixed point of this map and limits the build-up of cytoskeletal structures. In contrast, for negative values, F(𝐱n)<0, disassembly prevails, and we assume that ϵ then tends towards its lower bound ϵ0-Δϵ/2:

(S13b) ϵ(xn,t+Δt)=ϵ(xn,t)+Δtμ[ϵ0Δϵ/2ϵ(xn,t)],

where the time step is defined as Δt1. Neutral values, F(𝐱n,t)=0, lead to relaxation of ϵ towards a resting state

(S13c) ϵ(xn,t+Δt)=ϵ(xn,t)+Δtμ[ϵ0ϵ(xn,t)],

where the time step is defined as Δt1. The parameter μ signifies the rate at which cytoskeletal structures respond to the regulatory factors F. For the parameters and cell sizes used in this work (ϵ0=𝒪(100) and Δϵ=𝒪(10), and each cell occupying approximately 1000 grid sites) we set μ=0.1.

Appendix 1—figure 2
Distribution of regulatory factors on the basis of accepted elementary events.

For ease of reference, grid rows have been numbered from 1 to 10. Left (A): Solid black lines indicate cells’ membrane positions after acceptance of the respective elementary event; colors indicate cellular domains before the respective elementary event has been accepted (gray: substrate; shades of yellow: cells). Blue and red circular arcs (of radius R) delineate areas of local increase or decrease in the level of regulatory factors, respectively. The following elementary events are depicted: (i) lower cell retracts (two grid sites in row 2); (ii) lower cell protrudes (row 5); (iii) upper cell protrudes (row 10). In addition, the following elementary events occur across the cell-cell boundary: (iv) retraction of upper cell leads to rupture of cell-cell contacts (row 6, right event); (v) either the lower cell protrudes and pushes the upper cell or the upper cell retracts and pulls on the lower cell (row 6, left event). Specifically, event (v) entails mechanical signaling between the upper and lower cell and, therefore, affects the distribution of regulatory factors in both cells. Right (B): Identical copy of the left image (A). Colors indicate local levels of regulatory factors F (blue: F is positive; white: F is zero; red: F is negative; gray: substrate site). Note, in particular, that a substrate grid site has been inserted where cell rupture occurred (row 6, right grid site). The following cases can be distinguished: (i) Grid site xk lies in the zone of influence of only positive (blue circles) or negative (red circles) chemical feedback, in which case the level of regulatory factors is positive or negative, respectively (e.g. red grid sites in row 2, or blue grid sites in row 5). (ii) Grid site xk lies outside of any zone of influence, in which case the level regulatory factors is zero (e.g. white grid sites in row 2). (iii) Grid site xk lies in the zone of influence of equally many positive and negative feedbacks, in which case the level of regulatory factors remains zero (e.g. fourth grid site in row 4). (iv) Grid site xk lies in a zone of predominantly positive or negative feedback, in which case the level of regulatory factors is positive or negative, respectively (e.g. third grid site in row 4). Recall that only the sign of F is of significance to update the cells’ polarization field; cf. Equation S13.

After this update procedure for ϵ(𝐱n,t) is completed, all regulatory factors are reset, F(xn)0, n. This prevents ‘spurious memory effects’ which may arise once the cell’s rear reaches its initial leading edge position as time goes on. In essence, resetting regulatory factors upon completion of one MCS implies that the diffusion-degradation dynamics, underlying the distribution of regulatory factors, is fast on the scale of one MCS.

We emphasize that the polarization field ϵ(𝐱n) is defined only for grid sites 𝐱n𝒟(α) occupied by an actual cell (α0). To allow for spatial variations of substrate properties, we therefore introduce a second scalar field variable φ(𝐱n), which is defined on the entire computational grid. The scalar field φ(𝐱n) is taken to measure the local density of anchoring points that a cell might use to form focal adhesions. Although one might consider to treat φ as a time-dependent field variable, in this work φ is used to implement static substrate patterns, only. The field φ(𝐱n) is thus initialized once at the beginning of the simulation and kept fixed throughout the entire simulation.

Having introduced the fields ϵ(𝐱n) and φ(𝐱n), we now discuss their impact on the system’s dynamics by giving their contribution to the goal function p(𝒯). Suppose that the elementary event 𝒯 is attempted by a source cell α at source grid site 𝐱s(α). Further, let 𝐱t denote the target grid site and β denote the index of the target cell (where as ususal β0 indicates that 𝐱t is occupied by an actual cell, β<0 means that 𝐱t exposes substrate). We then define the ‘polarization energy’ Δcyto(𝒯) as follows:

(S14a) Δcyto(𝒯){ϵ(xt)ϵ(xs),𝒯=^𝒯pro  β0,ϵ(xs)ϵ(xt),𝒯=^𝒯ret  β0,[ϵ(xs)+φ(xt)],𝒯=^𝒯pro  β<0,ϵ(xs)+φ(xt),𝒯=^𝒯ret  β<0,

Here, the definition of Δcyto is such that the likelihood of cell protrusions is enhanced if the concentration of cytoskeletal structures at the source grid site, ϵ(𝐱s), is larger than the concentration at the target grid site, ϵ(𝐱t) (first row of Equation S14a), and vice versa for cell retractions (second row of Equation S14a). The strength of focal adhesions is taken to be measured by the sum ϵ+φ. Their associated ‘anchoring effects’ (which increase with growing strength of focal adhesions) promote the formation of cell protrusions against unoccupied substrate sites (third row of Equation S14a), and, correspondingly, impede cell retractions (fourth row of Equation S14a). Note, in particular, that the first two rows of Equation S14a can be obtained from a combined evaluation of the lower two rows. For example, if source cell α annexes xt starting from xs, two things need to happen: First, focal adhesions formed by the target cell β must be broken, implying a contribution Δcyto=ϵ(𝐱t)+φ(𝐱t) (fourth row of Equation S14a). Secondly, new focal adhesions are formed by the source cell α, implying a contribution Δcyto=-[ϵ(𝐱s)+φ(𝐱t)] (third row of Equation S14a). Taking the sum of both contributions gives the expression in the first row of Equation S14a. An analogous line of arguments leads to the expression in the second row of Equation S14a.

The contribution to the goal function p(𝒯) due to the polarization energy Δcyto(𝒯) is then defined by

(S14b) pcyto(𝒯):=exp[Δcyto(𝒯)/kBT],

where we set the effective thermal energy to kBT1. The characteristic ‘energy scale’ for Δcyto is set by the polarization bounds ϵ0-Δϵ/2 and ϵ0+Δϵ/2, which turns out to have important implications for collective cell dynamics, as discussed in the main text.

Cell adhesion

To implement the ability of cells to establish cell adhesions across cell-cell interfaces, we use a special form for the respective contribution to the goal function p, which is designed to distinguish between the formation of new and the breakage of existing cell-cell adhesion sites.

To this end, we define adhesion matrices Bα,β and Bα,β quantifying the system’s change in ‘energy’ upon forming a new contact between cells α and β [Bα,β] and upon breaking a pre-existing contact between those cells by an ‘intruder cell’ γα,β [Bα,β]. In our computational model, we assume that formation of new cell-cell contacts is energetically favored, and that breaking of pre-existing contacts by intruder cells is energetically penalized. The matrix entries of Bα,β and Bα,β, therefore, have a definite sign, which we take to be positive. The ordering of the cell index pair of Bα,β and Bα,β is of no physical significance, i.e. the adhesion matrices are symmetric. We also assume that a given cell α 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α,β and Bα,β exhibit the following properties:

(S15a) Bα,β=Bβ,α0,
(S15b) Bα,β=Bβ,α0,
(S15c) Bα,α=Bα,α=0,
(S15d) Bα,β=Bα,β=0,if\ α<0  β<0.

In addition, we assume that the energy cost associated with the breakage of a given cell-cell contact exceeds the energetic benefit of its initial formation, i.e.

(S15e) Bα,βBα,β,

where equality of both quantities would reproduce the assumption underlying the standard cellular Potts model (Graner and Glazier, 1992; Glazier and Graner, 1993). We shall refer to this property as the ‘dissipative nature of cell-cell adhesion’.

To implement the effects of cell-cell adhesion, we compute the ‘energy difference’ Δadh(𝒯) for any given elementary event 𝒯 according to the scheme illustrated in Appendix 1—figure 3. One has to distinguish between protrusion and retraction events. First, say that a cell α attempts a protrusion event 𝒯pro, involving the source grid site 𝐱s(α) and the target grid site 𝐱t(β), as illustrated in Appendix 1—figure 3A. In this case, cell α acts as intruder cell, since the depicted protrusion event affects three pre-existing contacts between the target cells β and a ‘third party’ cell γ. Acceptance of the depicted protrusion event would have the following energetically relevant effects: (i) All pre-existing contacts between the target cell β and third party cell γα,β at the target grid site 𝐱t are torn apart. (ii) New contacts between the source cell α and third party cell γα,β are established. (iii) The length of the cell contact line between source cell α and target cell β is changed. Altogether, these three effects lead to the following cell adhesion energy difference,

(S16a) Δadh(𝒯pro)xj𝒩t[Bα,c(xj)δα,c(xj)Bβ,c(xj)]+xj𝒩tBβ,c(xj)(1δα,c(xj)),

where is the length of a hexagon edge. The first term in this expression accounts for the (energetically favored) formation of new cellular contacts, as well as for the remodeling of the interface between source cell α and target cell β [points (ii) and (iii)]. The second term measures the (energetically penalized) breaking of pre-existing cell contacts [point (i)] and, therefore, impedes cell α’s ability to intrude. Conversely, if source cell α attempts a retraction event 𝒯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 𝐱t to the source site 𝐱s; cf. Appendix 1—figure 3A:

(S16b) Δadh(𝒯ret)xj𝒩s[Bβ,c(xj)δβ,c(xj)Bα,c(xj)]+xj𝒩sBα,c(xj)(1δβ,c(xj)),

where is the length of a hexagon edge.

Appendix 1—figure 3
Cell-cell adhesion.

(A) Adhesive energy contribution in a cyclic process, where a protrusion of source cell α against target cell β is followed by the inverse retraction event. Both events involve a third party cell γ, leading to net energy dissipation after the cyclic process has been completed. Protrusion: (i) Three pre-existing cell-cell contacts between β and γ are torn apart (red dashed contacts); (ii) three new contacts between α and γ are formed; (iii) the contact length between source cell α and target cell β increases by one unit of length. This implies Δadh(𝒯pro)=(3Bβ,γ3Bα,γBα,β). Retraction: (i) Three pre-existing cell-cell contacts between α and γ are torn apart (red dashed contacts); (ii) three new contacts between β and γ are formed; (iii) the contact length between source cell α and target cell β decreases by one unit of length. This implies Δadh(𝒯ret)=(3Bα,γ3Bβ,γ+Bα,β). Altogether, this leads to Δadh(cycl)=Δadh(𝒯pro)+Δadh(𝒯 ret)=(3(ΔB)α,γ+3(ΔB)β,γ)0, i.e. a (non-negative) dissipative contribution, whose magnitude depends on the dissipation matrix (ΔB)α,β=Bα,βBα,β0. (B) Shear viscosity due to cell-cell adhesion. Consider two rows of adhesive cells sliding past each other as indicated in the figure (left row of cells moves up by one grid site; colors indicate different cells). The associated adhesion energy change (per cell) reads Δadh/nc=2(BB)0, where nc denotes the number of cells sliding past each other, and where we assumed cells of like type, i.e. Bα,βB and Bα,βB (αβ). The condition B>B, Equation S15e, thus implies positive friction associated with cellular shear flows, whose magnitude is proportional to the number of cells sliding past each other. Note that this shear viscosity vanishes for B=B, i.e. for zero dissipation matrix.

We may now use Equation S16a and Equation S16b to illustrate the ‘dissipative nature’ of adhesive interactions by means of two archetypical examples. First, consider the adhesive energy contribution to any cyclic process. By a cyclic process we mean a sequence of two mutually inverse elementary events, e.g. a protrusion event 𝒯pro, which is immediately followed by its inverse retraction event 𝒯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:

(S17) Δadh(cycl)=xj𝒩t(𝒟(α)𝒟(β))[(ΔB)α,c(xj)+(ΔB)β,c(xj)],
(S18) (ΔB)α,β:=Bα,βBα,β0,

and can therefore conclude that

Δadh(cycl)0,

where 𝒩t denotes the neighborhood of the grid site which temporarily changes its cell index, and where α and β are the indices of the source and target cells involved in the cyclic process; cf. Appendix 1—figure 3A. Since (ΔB)α,β0, the above adhesive energy contribution is non-negative, thus leading to an amount of energy equal to Δadh(cycl) being dissipated as the cyclic process completes. This leads us to refer to the parameter matrix (ΔB)α,β as dissipation matrix. Second, consider two (infinitely extended) columns of cells in adhesive contact, sliding past each other. This situation is depicted in Appendix 1—figure 3B, where the left column of cells moves (as a whole) upwards by one grid site, while the right column of cells remains stationary. To assess the adhesive energy contribution along the contact line connecting both cell columns, note that the depicted transformation can be implemented by letting each cell in the left column protrude its leading (i.e. upper) edge by one grid site. For each protruding (source) cell α, this transformation entails to the following energetic effects (cf. discussion above): (i) Two of the pre-existing cell-cell contacts between the source cell’s upper neighbor in the left column (target cell β) and the corresponding cell in the right column (third party cell γ) get torn apart, leading to an energetic contribution 2Bβ,γ. (ii) In return, two new contacts between the protruding (source) cell α and cell γ are being established, leading to a contribution 2Bα,γ . (iii) Since the length of the contact line between cells in the left column (i.e. between protruding source cell α and retracting target cell β) 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α,βB and Bα,βB (αβ), and therefore, find

(S19) Δadh(visc)=2(BB)2ΔB0,

i.e. a non-negative dissipative contribution per cell. The size of the dissipation parameter Δ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(𝒯) as follows:

(S20) padh(𝒯):=exp[Δadh(𝒯)/kBT],

where we set the effective thermal energy to kBT1.

Rupture of cell contacts

By now, we have introduced all components making up the total acceptance probability p(𝒯), 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, 𝐱s and 𝐱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 𝐱s to 𝐱t [from 𝐱t to 𝐱s]. In other words, the domain 𝒟(α) of source cell α annexes 𝐱t [loses 𝐱s], while the domain 𝒟(β) of target cell β is forced to let go 𝐱t [accommodate 𝐱s]. However, if both source and target cells are actual cells, i.e. α,β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 α on its neighboring cells might also result in rupture of all pre-existing contacts between the retracting source cell and its neighboring cells at 𝐱s, rather than forcing one of its neighboring cells (the target cell) to fill the void created at 𝐱s once α 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 α on target cell β is strong enough to force β to fill the void created at 𝐱s (i.e. to accommodate 𝐱s), and call this a regular retraction event 𝒯ret. Secondly, we assume that the retraction of source cell α causes all pre-existing cell-cell contacts of cell α at 𝐱s to rupture, leaving a free substrate site at 𝐱s after the retraction event has been accepted. This latter event will be referred to as rupture event 𝒯rup. We then compute the total energy differences

Δ(𝒯ret)=Δcont(𝒯ret)+Δcyto(𝒯 ret)+Δadh(𝒯ret)

and

Δ(𝒯rup)=Δcont(𝒯rup)+Δcyto(𝒯 rup)+Δadh(𝒯rup)

under both assumptions (the energy difference associated with accepting 𝒯rup can be computed with Equation S16b by using the substrate β=-1 as new target cell) and compare the respective outcomes. If the rupture event is energetically favored over the regular retraction event, i.e. Δ(𝒯rup)<Δ(𝒯ret), then cohesion between cells is weak. In this case, the rupture event 𝒯rup, rather than the regular retraction event 𝒯ret will be attempted. Otherwise, cohesion is strong and a regular retraction event 𝒯ret will be attempted.

Rupture of substrate contacts

In our discussion so far, a cyclic process that follows up a protrusion event 𝒯pro with its inverse retraction event 𝒯ret, involving a cell α and no third party cells (in other words: no cell-cell contacts are made or broken), will not yield a net energy cost or gain; cf. Equation S14a. To account for the dissipative nature of cell-substrate contacts, we proceed similarly as when we have considered the disspative nature of cell-cell contacts. We introduce dissipation in substrate adhesion by leaving the Hamiltonian unaltered for protrusion events but adding a penalty for retraction events:

(S21) Δ(𝒯ret)Δ(𝒯ret)+D.

Therefore, a cell that adheres to the substrate at some grid point has to pay a cost D to retract from it. In other words, we assume that a fixed amount of energy D is dissipated once the adhesive bonds between a cell and the substrate break.

To keep its overall size across translations, the cell has to gain and lose equal amounts of hexagons, with Δϵ as the maximal energy gained by a single gain-and-loss in the absence of dissipation. In the presence of dissipation however, the cell has to pay at least a cost of (ϵ0Δϵ/2)+D to detach at an arbitrary location, resulting in Δϵ-D as the maximal energy gained by a single gain-and-loss in the presence of dissipation. Thus, while for D=0 there is no impact of substrate dissipation on cell motility, it will at the latest for D=Δϵ completely inhibit (directed) cell migration. Therefore, we study substrate dissipation in the range D[0,Δϵ].

Cell domain update routine

Having discussed the implementation concerning the two basic types of elementary events, namely protrusion events 𝒯pro and retraction events 𝒯ret, as well as the two subtypes of regular retraction events and rupture events 𝒯rup, we can now summarize the cell domain update routine, point 3.5 in section ‘Monte-Carlo scheme’. To this end, and in accordance with our previous notation, we use the cell indices α and β to denote source and target cell, and 𝐱s and 𝐱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:

    1. Set ϵ(𝐱t)=ϵ(𝐱s) and F(𝐱t)=F(𝐱s).

    2. 𝒟(α)𝒟(α){𝐱t}.

    3. 𝒟(β)𝒟(β){𝐱t}.

    4. Distribute regulatory factors according to Equation S12a.

  • If the accepted elementary event is a regular retraction event:

    1. Set ϵ(𝐱s)=ϵ(𝐱t) and F(𝐱s)=F(𝐱t).

    2. Set 𝒟(α)𝒟(α){𝐱s}

    3. Set 𝒟(β)𝒟(β){𝐱s}

    4. Distribute regulatory factors according to Equation S12b.

  • If the accepted elementary event is a rupture event:

    1. Set ϵ(𝐱s)=0 and F(𝐱s)=0.

    2. Set 𝒟(α)𝒟(α){𝐱s}

    3. Distribute regulatory factors according to Equation S12c.

Cell proliferation and mitosis

While cell proliferation and mitosis play no role in the experimental setup of rotating cell clusters, cell growth and division are observed experimentally in a setup where a sheet of cells expands into free space after removal of a stencil. Therefore, it is essential to include proliferation of cells in the numerical model. How this is done is described in this section.

We distinguish between two phases in the cell cycle, an interphase during which cells roughly double in volume and mitosis, the process of cell division. Even though a further partitioning of the interphase was considered in previous work (Li and Lowengrub, 2014), we do not expect that such a distinction is relevant for our results. In our computational framework cell growth may be implemented by progressively changing any cellular parameter that affects the cell’s equilibrium size. The two possible, largely equivalent choices are a successive reduction of the area coupling constant κA(α) or an increase of the average cell polarization ϵ0(α). We here employ the first method. We assume that individual cells grow exponentially (Barber et al., 2017) over a well-defined period Tg. 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 Tr, which we set to Tr=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 2. Therefore, as we do not scale the polarization field ϵ and the resulting energy gains for protrusions during cell growth, we mitigate the increased cost for ruffling the membrane by reducing the perimeter stiffness by a factor of 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 κA and κP are, therefore, kept constant for a quiescent cell; we denote the corresponding values as κA/0 and κ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 AT=rAref, where the factor r=𝒪(1) relates the threshold size to the equilibrium cell size Aref reached by a free, solitary cell with constant polarization field ϵ=ϵ0. Cells living in a densely packed environment will not exceed the area threshold due to the pressure exerted on them by neighboring cells and can, therefore, not grow. Conversely, cells exposed to free space are more likely to reach this threshold and proliferate. Finally, a growing cell in interphase becomes mitotic after the growth time Tg has passed, at which point cell size has roughly doubled with respect to the size in the quiescent period. We assume that cell migration and mitosis are processes that exclude each other. Hence, the positive feedback leading to persistent cell migration is switched off for mitotic cells and the polarization field relaxes to the neutral state ϵ0 according to Equation S13c.

There appears to be no universal set of rules which determine the orientation of the cleavage plane along which cells divide (Minc and Piel, 2012). Rather, for epithelial tissues there are a variety of factors which include local cell geometry and the direction of stress in the tissue (Gibson and Gibson, 2009). Though it is in principle possible to implement any given rule in our computational model, in its present version the axis along which a cell divides is chosen as a random direction through the geometric center of the cell. In case of irregular cell shapes, a separation of the cellular domain into more than two connected components can occur. To prevent violation of topological constraints, in this case the two largest components are considered as descendant cells and the residual grid sites are filled by substrate.

We explicitly account for the finite duration of the mitotic phase Td by keeping the cells in a mitotic state for the aforementioned time period, until the final instantaneous splitting of the cellular domains. After cell division, persistent cell migration of the daughter cells is enabled again. The descendent cells will subsequently re-enter the growing phase if their area exceeds the defined threshold, as mentioned above.

The following list summarizes the steps motivated and explained in the previous paragraphs. These additional steps are performed in a simulation that includes cell proliferation:

  1. Assign a state variable s(α) to each cell which encodes the current phase in the cell cycle:

    (S22) s(α)={0, quiescent\ phase1, refractory\ phase2, interphase3, mitotic\ phase
  2. Compute the equilibrium size Aref=(ϵ0-2π3κP)/(3κA) of a free, solitary cell with fixed polarization field ϵ=ϵ0, which spreads on the substrate used in the simulation.

  3. At the beginning of the simulation, t=0, all cells are in the quiescent state, s(α)(0)=0, and have the following area and perimeter coupling constants, respectively: κA(α)(0)=κA/0 and κP(α)(0)=κP/0.

  4. 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) s(α)(t)=0  A(α)(t)>ATs(α)(t+1)=1.
  • Switch from refractory state to growing state with probability p=1-exp(-1/Tr):

    (S24) s(α)(t)=1s(α)(t+1)={2,(p),1,(1p),

where the terms in the brackets denote the respective probability.

  • Exponential growth in interphase over a period of Tg:

    (S25) s(α)(t)=2κA(α)(t+1)=κA(α)(t)(1/2)1/TgκP(α)(t+1)=κP(α)(t)(1/2)1/(2Tg).
  • Switch from interphase to mitosis:

    (S26) s(α)(τ)=2 \ for \ all \ τ[tTg, t]s(α)(t+1)=3.
  • During cell division, cell motility is switched off and the polarization field relaxes to the neutral state according to Equation S13c.

  • Perform cell division, reset area and perimeter stiffness and exit mitotic phase:

    (S27) s(α)(τ)=3 \ for \ all \ τ[tTd, t] \ divide \ cell \ α \ into \ cells \ (α,β)κA(α,β)(t+1)=κA/0κP(α,β)(t+1)=κP/0s(α,β)(t+1)=0.
  • 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

(S28) σ¯ij(α)=2A(α)xk(α)(fkix~kj+fkjx~ki),

which is a discretized version of the surface integral in Landau et al. (1986). Here, fk is the force acting on the membrane element 𝐱k of cell α, x~k=xkxcom(α) is the position of the element with respect to the center of mass xcom(α) of the cell, and the superscripts i and j are Cartesian indices. The forces fk can be computed from the energy differences of all possible protrusion and retraction events originating from 𝐱k,

(S29) fk=xl𝒩kΔ(𝒯pro)xlxkxlxkxlxkxl𝒩k#Δ(𝒯ret)xkxlxkxlxkxl,

where the number sign indicates a sum over substrate grid sites only, i.e. grid sites with c(𝐱l)<0, and where Δcont+adh+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):

(S30) K=14πAP2.

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 𝒟(α):

(S31) (Cov(𝒟(α)))ij=xk𝒟(α)x~kix~kjxk𝒟(α)1,

where x~k=xkxcom(α) denotes the coordinates of element xk of cell α, relative to the cell's center of mass xcom(α); 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 1, 2 and 3 within confinements of varying size and constant cell density. In the 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 2 and 3-phases, we find scale-free behavior such that there is always a macroscopic rotation of the whole cell population regardless of the cell count and corresponding confinement size.

We also explore the parameter space of the tissue simulations. There, we find that an increased cell-cell dissipation ΔB impairs monolayer growth, while at the same time increasing front roughness. Similarly, an increased cell-substrate dissipation D also impairs monolayer growth. In contrast, increasing the maximum cell polarity Δϵ improves monolayer growth and also increases front roughness. We thus find that the speed of monolayer expansion depends on whether it is dominated by cell migration or cell proliferation, with the former improving monolayer growth through a better exploration of the cell-free area.

Single cell size

To rationalize our choice of the cell growth algorithm (see section ‘Cell proliferation and mitosis’), we have explored the shape and motility of differently sized solitary cells. To this end, we have varied the area stiffness parameter κA for different values of perimeter stiffness κP, while keeping all other parameters constant. We find that the area occupied by the motile cell increases linearly with 1/κA (Appendix 2—figure 1A). In particular, the cell area can be approximated quite well by the area of an immotile and equilibrated cell with uniform ϵ=ϵ0 (fit not shown):

(S32) A=ϵ02π3κP3κA1/κA.
Appendix 2—figure 1
Role of area stiffness κA for cell size and motility.

(A) The cell area increases linearly with 1/κA. The aspect ratio (B), speed (C) and persistence (D) of the cell decrease with increasing cell size. In the simulations, the area elasticity was varied in the interval κA[0.09, 0.18], and the membrane elasticity was chosen from κP{0.054,0.057,0.060,0.063,0.066}. Fixed parameters: average cell polarization field ϵ0=225; maximum cell polarity Δϵ=50; signaling radius R=5; cytoskeletal update rate μ=0.1; cell-substrate dissipation D=0; cell-substrate adhesion penalty φ=0.

Furthermore, we find that with increasing size, and all other parameters constant, cells become rounder, slower, and less persistent (Appendix 2—figure 1B-D). To intuitively explain this phenomenology, let us compare a cell of size Aref with a cell of size rAref, where r[1, 2], with the respective area stiffnesses κA and κA/r. While the smaller cell has a perimeter Pref, neglecting geometric effects the larger cell has a larger perimeter of approximately rPref. Hence, the larger cell has to pay a larger energy cost (roughly by a factor r) to ruffle its membrane by adding segments and therefore increasing its perimeter. Meanwhile, both the energy gain from the polarization field ϵ and the energy cost for increasing the cell area (as AκA=cst.) are the same for both cells. Due to the increased energy cost for adding membrane segments, larger cells finds find it more difficult to polarize, and are therefore rounder, slower and less persistent.

To offset this increased energy cost for adding membrane segments to the cell, one can scale the perimeter stiffness of the larger cell to κP/r, such that PκPcst. We would then predict that the ratio κP2/κ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 κ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 κP2/κ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, R2/A. Finally, we speculate that all observed quantities collapse unto respective master curves f(ΔϵκA/κP)g(RκA).

Single cell shape and dynamics depend on substrate dissipation

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

(S33) ω(t)=e^zv^α(t)×R~α(t)R~α(t)2𝒞,

averaged over the set 𝒞={α | is\ not\ substrate} of all cells in confinement, where v~α(t)=vα(t)vα(t)𝒞 and R~α(t)=Rα(t)Rα(t)𝒞 are the velocity and position of cell α relative to the cell cluster, respectively.

In the lowly polarizable 1-phase, small cell populations rotate in a highly synchronized way, and rotation is maximized for populations of 7 cells per confinement (Figure 4—figure supplement 1A). As can be inferred from the time traces and snapshots (Figure 4—figure supplement 1B, C), cells in small populations all synchronously move in the same direction at a given time and randomly switch between clockwise and counter-clockwise rotation; the switching rate decreases with increasing size of the cell population. Upon increasing the cell count and concomitantly the confinement size, global rotation of the cell population gradually vanishes (Figure 4—figure supplement 1A).

Unlike in the 1-phase, we observe that in the highly polarizable 2 and 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 |ω| on the population size N can be fitted by a power law of the form |ω|N1/2r01. This inverse proportionality between the average angular velocity and the confinement size r0 implies total rotational order, with every cell moving at a constant velocity |vrot|0.008 (vrotR~). Furthermore, in the 2, and 3-phases we have only scarcely observed switching of the rotational direction of cell clusters; e.g. for 4-cell clusters in the 3-phase.

Interestingly, fluctuations in the angular velocity (σω) change in a highly non-monotonic fashion with the cell count and concomitantly the confinement size. Certain cell counts exhibit especially high fluctuations of the mean angular velocity (e.g. 5 cells in the 1-phase, see Figure 4—figure supplement 1A; 3 or 10 cells in the 2-phase, see Figure 4—figure supplement 2A; 3 cells in the 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 cell-cell dissipation ΔB, cell-substrate dissipation D and maximum cell polarity Δϵ.

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 cell-free area with the same average speed. Hence, it is not needed to consider the two cell fronts separately, and we can instead consider the set of unsigned front positions X:=abs(X>/<). Then, we define the average front position as xF:=E(X) and the front roughness as σF2:=Var(X). In particular, we study the total growth of the tissue over the course of 500 MCS, which is captured by the maximal position of the front, max(xF), as well as the maximal roughness of the front max(σ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 Aref. Because the first daughter cells may only appear after 200 MCS have passed, we exclude this initial period from the measurements of the maximal front position and roughness, respectively. Additionally, we provide some exemplary time traces of the front evolution (Figure 5-figure supplement B,D,F).

First, we investigated how the monolayer expansion and front roughness depend on cell-substrate dissipation, ΔB (Figure 5—figure supplement 1A, B). Our simulations show that the cell sheet expands slower with increasing cell-cell dissipation Δ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 cell-cell dissipation ΔB (inset of Figure 5—figure supplement 1A).

We also investigated how the monolayer expansion and front roughness depend on cell-substrate dissipation, D (Figure 5—figure supplement 1C, D). Before we turn to the monolayer, let us recall the observed single-cell behavior in the previous section (see section 'Single cell shape and dynamics depend on substrate dissipation'): for high enough cell-substrate dissipation D (typically of the same order of magnitude as the maximum cell polarity Δϵ) cell migration is switched off (Figure 2—figure supplement 1). Extrapolating the single-cell results, we expect that the same holds also for collectives of cells and that cell migration does not play a role for high cell-substrate dissipation. Indeed, with increasing cell-substrate dissipation, the monolayer expands slower, until this effect appears to saturate at a threshold value D5 (Figure 5—figure supplement 1C, D). Following this line of argument, monolayer growth is slowed down if we suppress cell migration and thus move the cell monolayer towards a proliferation-dominated mode of expansion.

What about the inverse? Is the monolayer growth increased if we enhance cell migration and thus move the cell monolayer towards a migration-dominated mode of expansion? To test this hypothesis, we have analyzed how the monolayer growth and front roughness depend on the maximum cell polarity Δϵ. As predicted, monolayer growth increases with the maximum cell polarity Δϵ (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 Δϵ.

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

  1. Book
    1. Alberts B
    2. Johnson A
    3. Lewis J
    4. Raff M
    (2015)
    Molecular Biology of the Cell (Sixth Edition )
    Garland Science.
    1. Alt S
    2. Ganguly P
    3. Salbreux G
    (2017) Vertex models: from cell mechanics to tissue morphogenesis
    Philosophical Transactions of the Royal Society B: Biological Sciences 372:20150520.
    https://doi.org/10.1098/rstb.2015.0520
  2. Book
    1. Landau LD
    2. Pitaevskii LP
    3. Kosevich AM
    4. Lifshitz EM
    (1986)
    Theory of Elasticity (Third Edition)
    Butterworth-Heinemann.
  3. Book
    1. Milo R
    2. Phillips R
    (2015)
    Cell Biology by the Numbers (First Edition)
    Garland Science.
    1. Stokes CL
    2. Lauffenburger DA
    3. Williams SK
    (1991)
    Migration of individual microvessel Endothelial-Cells - Stochastic-Model and parameter measurement
    Journal of Cell Science 99 Pt 2:419–430.

Article and author information

Author details

  1. Florian Thüroff

    Arnold Sommerfeld Center for Theoretical Physics and Center for NanoScience, Department of Physics, Ludwig-Maximilians-Universität München, Munich, Germany
    Contribution
    Conceptualization, Software, Formal analysis, Validation, Investigation, Visualization, Methodology
    Contributed equally with
    Andriy Goychuk
    Competing interests
    No competing interests declared
  2. Andriy Goychuk

    Arnold Sommerfeld Center for Theoretical Physics and Center for NanoScience, Department of Physics, Ludwig-Maximilians-Universität München, Munich, Germany
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology
    Contributed equally with
    Florian Thüroff
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6776-9437
  3. Matthias Reiter

    Arnold Sommerfeld Center for Theoretical Physics and Center for NanoScience, Department of Physics, Ludwig-Maximilians-Universität München, Munich, Germany
    Contribution
    Software, Formal analysis, Validation, Investigation, Visualization, Methodology
    Competing interests
    No competing interests declared
  4. Erwin Frey

    Arnold Sommerfeld Center for Theoretical Physics and Center for NanoScience, Department of Physics, Ludwig-Maximilians-Universität München, Munich, Germany
    Contribution
    Conceptualization, Resources, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Methodology, Project administration
    For correspondence
    frey@lmu.de
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8792-3358

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

A two-part list of links to download the article, or parts of the article, in various formats.

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)

  1. Florian Thüroff
  2. Andriy Goychuk
  3. Matthias Reiter
  4. Erwin Frey
(2019)
Bridging the gap between single-cell migration and collective dynamics
eLife 8:e46842.
https://doi.org/10.7554/eLife.46842

Share this article

https://doi.org/10.7554/eLife.46842

Further reading

    1. Cell Biology
    Wonjo Jang, Kanishka Senarath ... Nevin A Lambert
    Tools and Resources

    Classical G-protein-coupled 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 medial-trans 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.

    1. Cell Biology
    2. Plant Biology
    Masanori Izumi, Sakuya Nakamura ... Shinya Hagihara
    Research Article

    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 live-cell monitoring of chloroplast morphology, we observed the formation of chloroplast budding structures in sugar-starved leaves. These buds were then released and incorporated into the vacuolar lumen as an autophagic cargo termed a Rubisco-containing 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 autophagy-related division does not require DYNAMIN-RELATED 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 chloroplast-associated isolation membrane.