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, $\mathscr{H}$, which is associated with the cell configuration. The cell is perpetually driven out of equilibrium by active reorganization of its actomyosin network and focal adhesions.

Figure 1

Download asset Open asset

First, we investigated the impact of varying levels of cell perimeter stiffness ${\kappa }_P$ and maximum cell polarity $\mathrm{\Delta }ϵ$ 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 ${\displaystyle \hat{\mathbf{v}}(t)\equiv \mathbf{v}(t)/\Vert \mathbf{v}(t)\Vert }$ ($\mathbf{𝐯}$: cell velocity) and (geometrical) center of mass position $\mathbf{𝐑}(t)$ during a total simulation time of ${\displaystyle T_{\text{sim}}={10}^4}$ Monte-Carlo steps (MCS). For each set of parameters, we performed ${\displaystyle 100}$ statistically independent simulations, from which we computed the mean squared displacement, $\text{MSD}(\tau )\equiv ⟨{[\mathbf{𝐑}(t+\tau )-\mathbf{𝐑}(t)]}^2⟩$, and the normalized velocity auto-correlation function, ${\displaystyle C(\tau )\equiv ⟨\hat{\mathbf{v}}(t+\tau )\cdot \hat{\mathbf{v}}(t)⟩}$. Here, ${\displaystyle ⟨\dots ⟩}$ denotes an average with respect to simulation time $t$ as well as over all ${\displaystyle 100}$ independent simulations.

These computer simulations show that the statistics of the migratory patterns is well described by a persistent random walk model (Stokes et al., 1991; Wu et al., 2014) with its two hallmarks: a mean square displacement that exhibits a crossover from ballistic to diffusive motion (Figure 2A), and on sufficiently long time scales an exponential decay of the velocity autocorrelation function $C(\tau )\propto e^{-\tau /{\tau }_p}$ (inset of Figure 2A). We determined the persistence time of directed migration, ${\tau }_p$, by fitting the mean squared displacement with a persistent random walk model. In addition, we also measured cell speed, ${\displaystyle v}$, and cell aspect ratio, $l_{+}/l_{-}$, to further characterize cell motility and shape. Surprisingly, for each of these variables we found a master curve that only depends on the ratio between cell polarizability and cell contractility, $\mathrm{\Delta }ϵ/{\kappa }_P$ (Figure 2B,D,E). This data collapse suggests $\mathrm{\Delta }ϵ/{\kappa }_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

Download asset Open asset

The cells’ persistence times of directed migration, speeds and aspect ratios all show a characteristic dependence on the specific cell polarizability. There is a threshold value for the specific polarizability, $\mathrm{\Delta }ϵ/{\kappa }_P\approx 500$, below which cells remain immobile (Figure 2B,D,E; grey regions). Above this threshold, the persistence time of directed migration, speed and aspect ratio increase markedly with the specific polarizability (Figure 2B,D,E). In our model, the area and perimeter stiffnesses refer to global and homogeneous cell contractility, while the cell polarization field drives cell migration. As discussed in ‘Gripping the surface through the cell cytoskeleton’, the cell polarization field does not explicitly distinguish between a local extensibility (e.g. due to actin polymerization), a local contractility (due to 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, $\mathrm{\Delta }ϵ$, and contractility, ${\kappa }_P$, are equally reduced by a blebbistatin-dependent prefactor, then the specific polarizability, $\mathrm{\Delta }ϵ/{\kappa }_P$, and the resulting cell phenomenology should remain unchanged. Indeed, blebbistatin treatment of keratocytes and keratocyte fragments was reported not to affect cell shape and speed to any significant degree (Wilson et al., 2010; Ofer et al., 2011). Therefore, blebbistatin treatment can either increase or decrease cell motility, depending on the cell type and possibly on the specific mechanism that drives cell migration.

Figure 3 with 1 supplement see all

Download asset Open asset

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.

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 ($\mathrm{\Delta }ϵ$), 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 ($\mathrm{\Delta }ϵ/{\kappa }_P=1000$; cf. Figure 3—video 1, top row), which show a qualitative shift in cell behavior:

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

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

• For large signaling radii $R=15$, we observed circular motion of the cell; because of the large signaling radius, signals originating from the trailing edge affect the leading edge of the cell and vice versa. Due to this circular motion, the cell exhibits a 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 ${\displaystyle ⟨c⟩=⟨\Vert {\mathrm{\partial }}_s\hat{\mathbf{v}}(s)\cdot \hat{\mathbf{v}}(s)\Vert ⟩}$, where $s$ is the contour length along the corresponding trajectory. Here, we averaged the tangent vector ${\displaystyle \hat{\mathbf{v}}(s)}$ over ${\displaystyle 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 ${\displaystyle 10}$ Monte-Carlo steps due to our choice of the cytoskeletal update rate, $\mu =0.1$). We find that the curvature of the trajectories has a pronounced minimum at large signaling radii (where the persistence time of directed migration vanishes), which indicates a transition from straight to circular trajectories (Figure 3E). Such a transition from persistent migration to single cell rotations was previously observed in experiments (Lou et al., 2015; Raynaud et al., 2016) and in theory (Reeves et al., 2018; Allen et al., 2018).

Biological cells are highly dynamic entities which constantly change shape and move around in space. To reflect this dynamic behavior computationally, the domain ${𝒟}^{(\alpha )}$ of cell $\alpha$ changes over time. The evolution of cell shape and position, as represented by ${𝒟}^{(\alpha )}$, proceeds via a succession of elementary events. In our numerical model, elementary events come in one of two basic flavors: protrusion events and retraction events. During a protrusion event, cell $\alpha$ (referred to as source cell) incorporates one grid site ${\mathbf{𝐱}}_t$ (referred to as target grid site) from its neighborhood ${𝒩}^{(\alpha )}$,

thereby increasing its cellular domain by one grid site. Here, the neighborhood of cell $\alpha$, ${𝒩}^{(\alpha )}$, is defined as

During a retraction event, source cell $\alpha$ expels one of its membrane grid sites ${\mathbf{𝐱}}_s\in {ℬ}^{(\alpha )}$,

where the set of membrane grid sites ${ℬ}^{(\alpha )}$ is defined as

Protrusion and retraction events are the numerical analogs of cell protrusions and cell retractions.

In implementing the reassignment rules, Equation S5 and Equation S7, we have to take into account that cellular domains must not overlap. For solitary cells moving in free space this does not imply any restrictions, and Equation S5 and Equation S7 apply directly. In the bulk of a confluent monolayer of adhesive cells, however, any protrusion of source cell $\alpha$ into the domain of cell $\beta$ (referred to as target cell) must be accompanied by a corresponding retraction event ${𝒟}_{\text{old}}^{(\beta )}\to {𝒟}_{\text{new}}^{(\beta )}={𝒟}_{\text{old}}^{(\beta )}\setminus \{{\mathbf{𝐱}}_t\}$, where ${\mathbf{𝐱}}_t$ denotes the target grid site annexed by cell $\alpha$. We emphasize, however, that the reverse is not generally true. If source cell $\alpha$ retracts, i.e. loses one of its boundary grid sites ${\mathbf{𝐱}}_s\in {ℬ}^{(\alpha )}$, the lost grid site ${\mathbf{𝐱}}_s$ faces either one of two conceivable fates: If, on the one hand, cohesion among cells is sufficiently strong (cf. section ‘Rupture of cell contacts’ for a definition of the notion ‘sufficiently strong’), then the retraction of cell $\alpha$ exerts a pulling force on one of its neighboring cells $\beta$ (the target cell) and forces the target cell to fill the emerging void at ${\mathbf{𝐱}}_s$, i.e. ${𝒟}_{\text{old}}^{(\beta )}\to {𝒟}_{\text{new}}^{(\beta )}={𝒟}_{\text{old}}^{(\beta )}\cup \{{\mathbf{𝐱}}_s\}$, where ${\mathbf{𝐱}}_s$ denotes the grid site lost by cell $\alpha$. On the other hand, if adhesion between cells is weak, then retraction of the source cell $\alpha$ can lead to a rupture of pre-existing cell contacts between $\alpha$ and other cells at ${\mathbf{𝐱}}_s$, such that the lost grid site ${\mathbf{𝐱}}_s$ becomes free space [$c({\mathbf{𝐱}}_s)=\alpha \ge 0\to c({\mathbf{𝐱}}_s)<0$]. Details on the actual implementation of cell rupture are discussed in section ‘Rupture of cell contacts’.

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 ${\displaystyle t+\mathrm{\Delta }t}$, where we set the time step to ${\displaystyle \mathrm{\Delta }t\equiv 1}$. One MCS consists in a series of attempts to perform elementary events, originating from randomly chosen membrane grid sites of randomly chosen cells. The duration of one MCS, i.e. the actual number of attempted elementary events, is chosen such that each of the cells’ membrane segments is given the opportunity to attempt, on average, one elementary event per MCS. During each MCS, cell domains ${𝒟}^{(\alpha )}$ as well as the numerical values of cell areas $A_{\alpha }$ and perimeters $P_{\alpha }$ are updated ‘on the fly’, while the cells’ polarization fields are updated only once at the end of each MCS; cf. section ‘Cytoskeletal structures and focal adhesion’ for the details of this update rule. The simulation then proceeds along the following Monte-Carlo scheme:

1. Initialize the cell population and define the duration of the simulation, i.e. the number of MCS, ${\displaystyle T_{\text{sim}}}$, 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 ${\displaystyle t=t+\mathrm{\Delta }t}$, where ${\displaystyle \mathrm{\Delta }t\equiv 1}$.

6. Repeat steps 3–5 while ${\displaystyle t<T_{\text{sim}}}$.

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 ${\mathbf{𝐱}}_s\in {ℬ}^{(\alpha )}$ of some cell $\alpha$, referred to as source cell. This membrane grid site will be referred to as source grid site. In addition, each elementary event involves a second grid site which lies in the neighborhood of the source grid site ${\mathbf{𝐱}}_s$ and which is not currently occupied by cell $\alpha$: ${\mathbf{𝐱}}_t\in {𝒩}_s\setminus {𝒟}^{(\alpha )}$. In what follows, this additional grid site ${\mathbf{𝐱}}_t$ will be referred to as target grid site. This grid site may either be an empty substrate site or a membrane site of another cell $\beta$, in which case the respective cell will be referred to as target cell. While the source grid site determines the location of the attempted elementary event, the target grid site determines the direction along which the elementary event is bound to proceed.

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 ${𝒯}_{\text{pro}}$ and retractions ${𝒯}_{\text{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 ${\displaystyle {\mathbf{x}}_t}$, or retracts from ${\displaystyle {\mathbf{x}}_s}$. Then, a decision is made with equal probabilities as to whether one attempts ${\displaystyle {\mathcal{𝒯}}_{\text{pro}}}$ or ${\displaystyle {\mathcal{𝒯}}_{\text{ret}}}$. In a second step, the goal function ${\displaystyle p}$ is used to compute the occurrence probability of the attempted event ${\displaystyle \mathcal{𝒯}}$. Finally, this elementary event ${\displaystyle \mathcal{𝒯}}$ is being accepted with probability ${\displaystyle p(\mathcal{𝒯})}$.

3. Choice of the goal function $p\mathit{}\mathrm{(}\mathrm{T}\mathrm{)}$. 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), ${\displaystyle K=\sum _{\alpha }{P}_{\alpha }/\ell }$, 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 ${\displaystyle \alpha }$ of the cell population. Because the cell membrane represents the border between lattice sites occupied by cell $\alpha$ and unoccupied by cell $\alpha$, the chosen membrane segment automatically defines the source grid site ${\displaystyle {\mathbf{x}}_s\in {\mathcal{ℬ}}^{(\alpha )}}$ and the corresponding target grid site ${\displaystyle {\mathbf{x}}_t\in {\mathcal{𝒩}}^{(\alpha )}\cap {\mathcal{𝒩}}_s}$.

3. With equal probability, choose whether to attempt a protrusion event (${\displaystyle {\mathcal{𝒯}}_{\text{pro}}}$) or a retraction event (${\displaystyle {\mathcal{𝒯}}_{\text{ret}}}$).

4. Compute the prospective acceptance probability $p({𝒯}_{\text{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 $\alpha$ and opponent cell $\beta$; for details see section ‘Cell domain update routine’.

6. If $n<K$, set $n\to n+1$ and then repeat steps 2 through 5.

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

with positive coupling constants ${\displaystyle {\kappa }_P^{(\alpha )}}$ and ${\displaystyle {\kappa }_A^{(\alpha )}}$ characterizing the contractility of cell $\alpha$; for empty substrate sites ($\alpha <0$) we set ${\kappa }_P^{(\alpha )}={\kappa }_A^{(\alpha )}=0$. According to Equation S10, the cell’s ‘contractile energy’ increases with increasing cell perimeter and increasing cell area. The model Hamiltonian ${\mathscr{H}}_{\text{cont}}$ can then be used to specify the contractile contribution to the goal function $p(𝒯)$. To this end, let $\mathrm{\Delta }{\mathscr{H}}_{\text{cont}}(𝒯)$ denote the contractile contribution to the energy difference entailed by accepting an elementary event $𝒯$. Following a standard Metropolis algorithm, we then define

where we set the effective thermal energy to ${\displaystyle k_{\text{B}}T\equiv 1}$. The contractile ‘energy’, Equation S10, is similar to the corresponding energy model commonly used in cellular Potts models (Ouchi et al., 2003). Unlike the standard cellular Potts model, however, where a target area and target perimeter are used to keep the simulated cells from collapsing, the energetic contribution in Equation S10 strictly contracts the cell’s body. As will be detailed in the next section, to counteract these contractile forces, we explicitly model cytoskeletal structures within each cell, which provide outward pushing forces to balance cell contraction.

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

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

1. The scalar polarization field ${\displaystyle ϵ}$ 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: ${\displaystyle ϵ({\mathbf{x}}_n,t)\in [{ϵ}_0-\mathrm{\Delta }ϵ/2,{ϵ}_0+\mathrm{\Delta }ϵ/2]}$. These bounds are cell-type specific. While the upper bound ${\displaystyle {ϵ}_0+\mathrm{\Delta }ϵ/2}$ mainly reflects the limited availability of protein resources, the lower bound ${\displaystyle {ϵ}_0-\mathrm{\Delta }ϵ/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 ${\displaystyle ϵ({\mathbf{x}}_n)}$, is regulated by a myriad of accessory proteins. In our computational model we simplify these complex processes by resorting to a single ‘bookkeeping variable’ which we will refer to as ‘regulatory factors’. Its local level is stored as an integer variable $F({\mathbf{𝐱}}_n)$ for each grid site ${\mathbf{𝐱}}_n\in {𝒟}^{(\alpha )}$. We use $F({\mathbf{𝐱}}_n)$ to implement the overall action of regulatory cytoskeletal proteins in an effective and collective manner. Specifically, since the formation of lamellipodial structures depends on active nucleation promoting factors (Pollard and Borisy, 2003), we assume that positive levels, $F({\mathbf{𝐱}}_n)>0$, reflect local conditions in support of network-assembly, whereas negative levels, $F({\mathbf{𝐱}}_n)<0$, represent predominantly degrading (or disassembly) conditions. For neutral levels, $F({\mathbf{𝐱}}_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 ${\displaystyle F({\mathbf{x}}_n)}$ to implement this propagation of mechanical information, perceived by cell $\alpha$ at its periphery ${ℬ}^{(\alpha )}$, across a certain spatial distance $R$. The local levels of $F({\mathbf{𝐱}}_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 ${\displaystyle {\mathbf{x}}_s\in {\mathcal{ℬ}}^{(\alpha )}}$ (source cell: ${\displaystyle \alpha }$; target cell: ${\displaystyle \beta }$), then for all sites ${\displaystyle {\mathbf{x}}_n}$ within a range ${\displaystyle R}$ (i.e. ${\displaystyle \Vert {\mathbf{x}}_n-{\mathbf{x}}_s\Vert <R}$) the integer variable signifying regulatory factors is incremented up and down for the protruding and the retracting cell, respectively:

• Similarly, if a retraction event has been accepted at the source site ${\displaystyle {\mathbf{x}}_s\in {\mathcal{ℬ}}^{(\alpha )}}$ , and the (local) cell contact between source cell ${\displaystyle \alpha }$ and target cell ${\displaystyle \beta }$ has remained intact, then within a range ${\displaystyle R}$ one applies the inverse update rule:

• If a retraction event has been accepted at the source site ${\displaystyle {\mathbf{x}}_s\in {\mathcal{ℬ}}^{(\alpha )}}$ , and in addition the (local) cell contact between source cell ${\displaystyle \alpha }$ and target cell ${\displaystyle \beta }$ has ruptured, then the regulatory factors are reduced only within a range ${\displaystyle R}$ in the retracting cell:

Finally, if the target grid site ${\mathbf{𝐱}}_t$ is not occupied by any cell (substrate is indicated by ${\displaystyle \beta <0}$) prior to the elementary event, then only the first two lines in the above update scheme apply.

By virtue of the above update scheme, Equation S12, ‘regulatory factors’ are continuously distributed across each cell’s domain ${𝒟}^{(\alpha )}$ as the current MCS progresses. At the end of each MCS, the accumulated (local) values of $F({\mathbf{𝐱}}_n)$ are used to update the local values of the polarization field $ϵ({\mathbf{𝐱}}_n)$ inside each cell $\alpha \ge 0$ (${\mathbf{𝐱}}_n\in {𝒟}^{(\alpha )}$): We assume that for positive values, ${\displaystyle F({\mathbf{x}}_n)>0}$, there is assembly of cytoskeletal structures and $ϵ$ is increased by an amount proportional to the distance of $ϵ$ from its upper bound ${ϵ}_0+\mathrm{\Delta }ϵ/2$:

where the time step is defined as ${\displaystyle \mathrm{\Delta }t\equiv 1}$. Thereby ${ϵ}_0+\mathrm{\Delta }ϵ/2$ is a fixed point of this map and limits the build-up of cytoskeletal structures. In contrast, for negative values, $F({\mathbf{𝐱}}_n)<0$, disassembly prevails, and we assume that $ϵ$ then tends towards its lower bound ${ϵ}_0-\mathrm{\Delta }ϵ/2$:

where the time step is defined as ${\displaystyle \mathrm{\Delta }t\equiv 1}$. Neutral values, $F({\mathbf{𝐱}}_n,t)=0$, lead to relaxation of $ϵ$ towards a resting state

where the time step is defined as ${\displaystyle \mathrm{\Delta }t\equiv 1}$. The parameter $\mu$ signifies the rate at which cytoskeletal structures respond to the regulatory factors $F$. For the parameters and cell sizes used in this work (${\displaystyle {ϵ}_0=\mathcal{𝒪}(100)}$ and ${\displaystyle \mathrm{\Delta }ϵ=\mathcal{𝒪}(10)}$, and each cell occupying approximately ${\displaystyle 1000}$ grid sites) we set $\mu =0.1$.

Appendix 1—figure 2

Download asset Open asset

After this update procedure for $ϵ({\mathbf{𝐱}}_n,t)$ is completed, all regulatory factors are reset, ${\displaystyle F({\mathbf{x}}_n)\to 0,\mathrm{\forall }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 $ϵ({\mathbf{𝐱}}_n)$ is defined only for grid sites ${\mathbf{𝐱}}_n\in {𝒟}^{(\alpha )}$ occupied by an actual cell ($\alpha \ge 0$). To allow for spatial variations of substrate properties, we therefore introduce a second scalar field variable $\phi ({\mathbf{𝐱}}_n)$, which is defined on the entire computational grid. The scalar field $\phi ({\mathbf{𝐱}}_n)$ is taken to measure the local density of anchoring points that a cell might use to form focal adhesions. Although one might consider to treat $\phi$ as a time-dependent field variable, in this work $\phi$ is used to implement static substrate patterns, only. The field $\phi ({\mathbf{𝐱}}_n)$ is thus initialized once at the beginning of the simulation and kept fixed throughout the entire simulation.

Having introduced the fields $ϵ({\mathbf{𝐱}}_n)$ and $\phi ({\mathbf{𝐱}}_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 $\alpha$ at source grid site ${\mathbf{𝐱}}_s\in {ℬ}^{(\alpha )}$. Further, let ${\mathbf{𝐱}}_t$ denote the target grid site and $\beta$ denote the index of the target cell (where as ususal $\beta \ge 0$ indicates that ${\mathbf{𝐱}}_t$ is occupied by an actual cell, $\beta <0$ means that ${\mathbf{𝐱}}_t$ exposes substrate). We then define the ‘polarization energy’ $\mathrm{\Delta }{\mathscr{H}}_{\text{cyto}}(𝒯)$ as follows:

Here, the definition of $\mathrm{\Delta }{\mathscr{H}}_{\text{cyto}}$ is such that the likelihood of cell protrusions is enhanced if the concentration of cytoskeletal structures at the source grid site, $ϵ({\mathbf{𝐱}}_s)$, is larger than the concentration at the target grid site, $ϵ({\mathbf{𝐱}}_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 $ϵ+\phi$. Their associated ‘anchoring effects’ (which increase with growing strength of focal adhesions) promote the formation of cell protrusions against unoccupied substrate sites (third row of Equation S14a), and, correspondingly, impede cell retractions (fourth row of Equation S14a). Note, in particular, that the first two rows of Equation S14a can be obtained from a combined evaluation of the lower two rows. For example, if source cell $\alpha$ annexes ${\displaystyle {\mathbf{x}}_t}$ starting from ${\displaystyle {\mathbf{x}}_s}$, two things need to happen: First, focal adhesions formed by the target cell $\beta$ must be broken, implying a contribution $\mathrm{\Delta }{\mathscr{H}}_{\text{cyto}}=ϵ({\mathbf{𝐱}}_t)+\phi ({\mathbf{𝐱}}_t)$ (fourth row of Equation S14a). Secondly, new focal adhesions are formed by the source cell $\alpha$, implying a contribution $\mathrm{\Delta }{\mathscr{H}}_{\text{cyto}}=-\left[ϵ({\mathbf{𝐱}}_s)+\phi ({\mathbf{𝐱}}_t)\right]$ (third row of Equation S14a). Taking the sum of both contributions gives the expression in the first row of Equation S14a. An analogous line of arguments leads to the expression in the second row of Equation S14a.

The contribution to the goal function $p(𝒯)$ due to the polarization energy $\mathrm{\Delta }{\mathscr{H}}_{\text{cyto}}(𝒯)$ is then defined by

where we set the effective thermal energy to $k_{\text{B}}T\equiv 1$. The characteristic ‘energy scale’ for $\mathrm{\Delta }{\mathscr{H}}_{\text{cyto}}$ is set by the polarization bounds ${ϵ}_0-\mathrm{\Delta }ϵ/2$ and ${ϵ}_0+\mathrm{\Delta }ϵ/2$, which turns out to have important implications for collective cell dynamics, as discussed in the main text.

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_{\alpha ,\beta }$ and $B_{\alpha ,\beta }^{\prime }$ quantifying the system’s change in ‘energy’ upon forming a new contact between cells $\alpha$ and $\beta$ [$B_{\alpha ,\beta }$] and upon breaking a pre-existing contact between those cells by an ‘intruder cell’ $\gamma \ne \alpha ,\beta$ [$B_{\alpha ,\beta }^{\prime }$]. In our computational model, we assume that formation of new cell-cell contacts is energetically favored, and that breaking of pre-existing contacts by intruder cells is energetically penalized. The matrix entries of $B_{\alpha ,\beta }$ and $B_{\alpha ,\beta }^{\prime }$, therefore, have a definite sign, which we take to be positive. The ordering of the cell index pair of $B_{\alpha ,\beta }$ and $B_{\alpha ,\beta }^{\prime }$ is of no physical significance, i.e. the adhesion matrices are symmetric. We also assume that a given cell $\alpha$ does not interact with itself, such that the diagonal elements of the adhesion matrices vanish. Finally, there is no adhesion between cells and empty substrate sites, such that all matrix elements containing a negative cell index vanish. In summary, the adhesion matrices $B_{\alpha ,\beta }$ and $B_{\alpha ,\beta }^{\prime }$ exhibit the following properties:

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

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

To implement the effects of cell-cell adhesion, we compute the ‘energy difference’ $\mathrm{\Delta }{\mathscr{H}}_{\text{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 $\alpha$ attempts a protrusion event ${𝒯}_{\text{pro}}$, involving the source grid site ${\mathbf{𝐱}}_s\in {ℬ}^{(\alpha )}$ and the target grid site ${\mathbf{𝐱}}_t\in {ℬ}^{(\beta )}$, as illustrated in Appendix 1—figure 3A. In this case, cell $\alpha$ acts as intruder cell, since the depicted protrusion event affects three pre-existing contacts between the target cells $\beta$ and a ‘third party’ cell $\gamma$. Acceptance of the depicted protrusion event would have the following energetically relevant effects: (i) All pre-existing contacts between the target cell $\beta$ and third party cell $\gamma \ne \alpha ,\beta$ at the target grid site ${\mathbf{𝐱}}_t$ are torn apart. (ii) New contacts between the source cell $\alpha$ and third party cell $\gamma \ne \alpha ,\beta$ are established. (iii) The length of the cell contact line between source cell $\alpha$ and target cell $\beta$ is changed. Altogether, these three effects lead to the following cell adhesion energy difference,

where $\mathrm{\ell }$ is the length of a hexagon edge. The first term in this expression accounts for the (energetically favored) formation of new cellular contacts, as well as for the remodeling of the interface between source cell $\alpha$ and target cell $\beta$ [points (ii) and (iii)]. The second term measures the (energetically penalized) breaking of pre-existing cell contacts [point (i)] and, therefore, impedes cell $\alpha$’s ability to intrude. Conversely, if source cell $\alpha$ attempts a retraction event ${𝒯}_{\text{ret}}$, then the same reasoning as the one leading to Equation S16a applies, only this time the elementary event proceeds in reverse, i.e. from the target site ${\mathbf{𝐱}}_t$ to the source site ${\mathbf{𝐱}}_s$; cf. Appendix 1—figure 3A:

where $\mathrm{\ell }$ is the length of a hexagon edge.

Appendix 1—figure 3

Download asset Open asset

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 ${\displaystyle {\mathcal{𝒯}}_{\text{pro}}}$, which is immediately followed by its inverse retraction event ${\displaystyle {\mathcal{𝒯}}_{\text{ret}}}$, such that the system’s final configuration is identical to its initial configuration. Using Equation S16 we find for the total adhesive energy contribution to a cyclic process:

and can therefore conclude that

where ${\displaystyle {\mathcal{𝒩}}_t}$ denotes the neighborhood of the grid site which temporarily changes its cell index, and where ${\displaystyle \alpha }$ and ${\displaystyle \beta }$ are the indices of the source and target cells involved in the cyclic process; cf. Appendix 1—figure 3A. Since ${\displaystyle (\mathrm{\Delta }B{)}_{\alpha ,\beta }\ge 0}$, the above adhesive energy contribution is non-negative, thus leading to an amount of energy equal to ${\displaystyle \mathrm{\Delta }{\mathcal{\mathscr{H}}}_{\text{adh}}^{(\mathrm{c}\mathrm{y}\mathrm{c}\mathrm{l})}}$ being dissipated as the cyclic process completes. This leads us to refer to the parameter matrix ${\displaystyle (\mathrm{\Delta }B{)}_{\alpha ,\beta }}$ as dissipation matrix. Second, consider two (infinitely extended) columns of cells in adhesive contact, sliding past each other. This situation is depicted in Appendix 1—figure 3B, where the left column of cells moves (as a whole) upwards by one grid site, while the right column of cells remains stationary. To assess the adhesive energy contribution along the contact line connecting both cell columns, note that the depicted transformation can be implemented by letting each cell in the left column protrude its leading (i.e. upper) edge by one grid site. For each protruding (source) cell ${\displaystyle \alpha }$, 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 ${\displaystyle \beta }$) and the corresponding cell in the right column (third party cell ${\displaystyle \gamma }$) get torn apart, leading to an energetic contribution ${\displaystyle 2B_{\beta ,\gamma }^{\prime }}$. (ii) In return, two new contacts between the protruding (source) cell ${\displaystyle \alpha }$ and cell ${\displaystyle \gamma }$ are being established, leading to a contribution ${\displaystyle -2B_{\alpha ,\gamma }}$ . (iii) Since the length of the contact line between cells in the left column (i.e. between protruding source cell ${\displaystyle \alpha }$ and retracting target cell ${\displaystyle \beta }$) remains unchanged, there’s no further energetic contribution due to adhesive contacts between cells in the left column. Assuming that all cells in the system are of equal types, we write ${\displaystyle B_{\alpha ,\beta }\equiv B}$ and ${\displaystyle B_{\alpha ,\beta }^{\prime }\equiv B^{\prime }}$ (${\displaystyle \alpha \ne \beta }$), and therefore, find

i.e. a non-negative dissipative contribution per cell. The size of the dissipation parameter ${\displaystyle \mathrm{\Delta }B}$ thus introduces a natural means to tune the system’s shear viscosity.

With the above definitions of the adhesive energy changes, Equation S16, we define the contribution of cell adhesion to the goal function ${\displaystyle p(\mathcal{𝒯})}$ as follows:

where we set the effective thermal energy to ${\displaystyle k_{\text{B}}T\equiv 1}$.

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, ${\mathbf{𝐱}}_s$ and ${\mathbf{𝐱}}_t$, have been selected, acceptance of a protrusion [retraction] event causes (among other things like the distribution of regulatory factors) the cell index to be copied from ${\mathbf{𝐱}}_s$ to ${\mathbf{𝐱}}_t$ [from ${\mathbf{𝐱}}_t$ to ${\mathbf{𝐱}}_s$]. In other words, the domain ${𝒟}^{(\alpha )}$ of source cell $\alpha$ annexes ${\mathbf{𝐱}}_t$ [loses ${\mathbf{𝐱}}_s$], while the domain ${𝒟}^{(\beta )}$ of target cell $\beta$ is forced to let go ${\mathbf{𝐱}}_t$ [accommodate ${\mathbf{𝐱}}_s$]. However, if both source and target cells are actual cells, i.e. $\alpha ,\beta \ge 0$, and if the source cell attempts a retraction event, there is one additional possible outcome: If cell cohesion is weak, then the pulling force exerted by the retracting source cell $\alpha$ on its neighboring cells might also result in rupture of all pre-existing contacts between the retracting source cell and its neighboring cells at ${\mathbf{𝐱}}_s$, rather than forcing one of its neighboring cells (the target cell) to fill the void created at ${\mathbf{𝐱}}_s$ once $\alpha$ retracts; cf. rupture event depicted in Appendix 1—figure 2. To test for the occurrence of cell rupture, the total energetic cost of each attempted retraction event between two actual cells is evaluated under two different assumptions: First, we assume that the pulling force exerted by the retracting source cell $\alpha$ on target cell $\beta$ is strong enough to force $\beta$ to fill the void created at ${\mathbf{𝐱}}_s$ (i.e. to accommodate ${\mathbf{𝐱}}_s$), and call this a regular retraction event ${𝒯}_{\text{ret}}$. Secondly, we assume that the retraction of source cell $\alpha$ causes all pre-existing cell-cell contacts of cell $\alpha$ at ${\mathbf{𝐱}}_s$ to rupture, leaving a free substrate site at ${\mathbf{𝐱}}_s$ after the retraction event has been accepted. This latter event will be referred to as rupture event ${𝒯}_{\text{rup}}$. We then compute the total energy differences

and

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

In our discussion so far, a cyclic process that follows up a protrusion event ${\displaystyle {\mathcal{𝒯}}_{\text{pro}}}$ with its inverse retraction event ${\displaystyle {\mathcal{𝒯}}_{\text{ret}}}$, involving a cell ${\displaystyle \alpha }$ 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:

Therefore, a cell that adheres to the substrate at some grid point has to pay a cost ${\displaystyle 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 ${\displaystyle \mathrm{\Delta }ϵ}$ 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 ${\displaystyle ({ϵ}_0-\mathrm{\Delta }ϵ/2)+D}$ to detach at an arbitrary location, resulting in $\mathrm{\Delta }ϵ-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=\mathrm{\Delta }ϵ$ completely inhibit (directed) cell migration. Therefore, we study substrate dissipation in the range ${\displaystyle D\in [0,\mathrm{\Delta }ϵ]}$.