The computational model used here is an extended version of a previously introduced phase field model coupled to a stochastic bistable reaction-diffusion process [40, 44]. This model was originally designed to study the migratory behavior of an individual amoeba, such as a D. discoideum cell, and is now extended to include interactions between several cells that are jointly simulated in the same computational domain. The model is composed of a concentration variable c accounting for the biochemical processes occurring in the interior of the cell. The variable c represents a generic activatory signal that can be associated with upstream regulatory components of actin activity, such as activated Ras, PI3K or PIP₃. The dynamics of c is coupled to the evolution of the cell shape such that the cell boundary is protruding outward at locations of high values of c. The cell shape is encoded by an auxiliary phase field ϕ which defines the area of the cell and varies from ϕ = 1 inside to ϕ = 0 outside of the cell. In our simulations, each cell (index i) is represented by its own concentration variable c _i (x, t) and its own phase field ϕ _i (x, t). Different cells interact via short-range repulsion that prevents overlap or adhesion. The phase field evolves according to the following equation, τ ϕ ∂ ϕ i ∂ t = γ ∇ 2 ϕ i − G ′ ϕ i ϵ 2 − β ∫ ϕ i d A − A 0 ∇ ϕ i + α ϕ i c i ∇ ϕ i − m r e p ∑ j = 1 , i ≠ j N ϕ i ϕ j ∇ ϕ i , (1) where the first term on the right-hand side corresponds to the surface energy of the cell membrane related with a free energy which incorporates the derivative of a double well potential G(ϕ) = 18 ϕ ² (1—ϕ)². The surface tension and the width of the cell boundary are parametrized by γ and ϵ, respectively. The second term enforces the cell area to stay close to a typical value A ₀, and the third term represents the active force that acts on the cell membrane, which is generated by the biochemical field c _i . These three terms can be derived from free energy arguments [34]. The fourth term models cell-cell interactions controlled by parameter m _rep , thereby preventing cells from overlapping. Here, N corresponds to the total number of cells in the system. The parameters that control the surface tension, the strength of the area conservation constraint and the active force, γ, β and α, respectively, are kept constant in simulations, see Table 1 for parameter values, taken to reproduce the dynamics of vegetative or starving D. discoideum cells [40].

The dynamics of the biochemical component c _i that diffuses and reacts inside the cell follows a noisy, bistable reaction kinetics, ∂ ϕ i c i ∂ t = ∇ ϕ i D ∇ c i + ϕ i k a c i 1 − c i c i − δ c i − ρ c i + ϕ i 1 − ϕ i ξ i , (2) where k _a is the reaction rate, ρ the degradation rate, and D the diffusivity of the component c _i . The parameter δ controls the evolution of the concentration pattern and introduces a global feedback to maintain a fixed amount of c _i inside the cell, thus imposing a mass-conservation condition in the model: δ c i = δ 0 + M ∫ c i d A − C 0 . (3)

The last term on the right-hand side of Eq. 2 introduces noise, accounting for the stochastic nature of the reaction-diffusion processes occurring within the cell: the stochastic field ξ(x, t) follows an Ornstein-Uhlenbeck dynamics, ∂ ξ ∂ t = − k η ξ + η , (4) where η(x, t) is a Gaussian white noise with zero mean, ⟨η⟩ = 0, and the variance ⟨η(x, t)η(x′, t′)⟩ = 2σ ² δ(x − x′)δ(t − t′). The relaxation rate k _η along with the reaction rate k _a are key parameters in our model that control transitions between different forms of cell motility.

Simulations were performed using finite differences with a spatial and temporal resolution of Δx = 0.15 μm and Δt = 0.002 s, respectively. We used the Euler-Maruyama method for the stochastic integration of the partial differential equations. The total area considered for each cell was kept constant at A ₀ = 113 μm ² corresponding to a circular cell with radius r = 6 μm. The area covered by the biochemical component c _i was maintained at C ₀ = 28 μm ², corresponding to a quarter of the cell area. The size of the grid and the number of cells were varied to explore different packing fractions ρ _R = NA ₀/L ² in the range from 0.24 to 0.78. Initial conditions of the intracellular concentration fields are chosen to produce polarized cells and promote collisions in the case of binary interaction scenarios, and isotropic in the case of multi-particle simulations.

Some characteristic cell morphologies as observed in numerical simulations are displayed in Figure 1. Patches of high concentration of the biochemical component c _i that typically result in extensions/protrusions of the cell boundary are shown in green color in the respective images (cf. Figure 1A). Similar to earlier experimental work, where fluorescently marked patches of PI3K were used to quantify the polarization of cells [60], we take advantage of the distribution of the biochemical field c _i for that purpose here as follows. We measure the total area of the cell (A _T ) and obtain the corresponding centroid coordinate (x _T , y _T ), shown as a white cross in Figure 1. We also measure the total area(s) of the patches covered by high values of the field c _i (A _Ci ) and obtain the centroid coordinates ( x C i , y C i ) of these area(s) (black crosses in Figure 1). If multiple patches exist, we take the weighted average of the individual centroid coordinates to obtain the coordinate (x _CT , y _CT ), highlighted by a yellow cross in Figure 1. Eventually, we define a polarization vector P as the distance between the centroid coordinates (x _T , y _T ) and the coordinates (x _CT , y _CT ). Moreover, a velocity vector of cell propagation V is defined as the distance between the centroids of the cell at time t and at subsequent time t + Δt divided by the numerical time step Δt. In addition, we furthermore calculate the angle θ between the two vectors P and V. In Figures 1B–D, the corresponding definitions of the coordinates, vectors, and the angle are shown.

We divide all possible collisions between two cells that we can generate in numerical simulations, into two specific types: frontal and glancing collisions, see Figure 2A. For both types, we observed three main interaction scenarios with different outcomes as illustrated in Figure 2. In the first scenario, which we call alignment, both cells face and migrate into the same direction after the collision. In the second case, the anti-alignment scenario, cells repel each other and move away in opposite directions after the collision. Finally, upon collision, the motion of the two cells may also stall and both cells may remain on the spot, pushing head-on against each other for an extended period of time, before they eventually get released by a random fluctuation. We call this the stuck/push scenario.

Cell polarity is a fundamental feature that plays a key role in many cellular functions, such as cell growth, division, and migration. In particular, depending on their degree of polarity, motile cells can display several different motion patterns, such as random, oscillatory or persistent movement. For a motile cell, polarity is typically defined based on the leading and trailing edge (head and tail) of the moving cell. In the model used here, it is the biochemical component c _i which triggers the formation of membrane protrusions (pseudopods) at the cell boundary, and is thus responsible for setting the sense of orientation of the moving cell. For this reason, we quantify polarity based on a polarity vector P that measures the asymmetry in the subcellular distribution of the biochemical component c _i , see Section 2.4 for details. To what extent the displacement of the cell is governed by the direction of the polarity vector is the subject of study of this section.

A wide range of cell motion patterns was studied in numerical simulations of our model previously [40]. It was shown that the transition from random to persistent motion is controlled by the model parameter k _a , see Eq. 2. Small values of k _a give rise to an erratic trajectory that is caused by the random appearance of protrusions all around the cell. On the other hand, large values of k _a produce a more persistent trajectory due to the accumulation of protrusions in one region of the membrane that sets the overall direction of motion. Those scenarios were associated with the vegetative and the starvation-developed states of D. discoideum cells [40], respectively.

Here, we studied cell migration patterns for different values of the parameter k _a , focusing on the dynamics of the displacement and polarity vectors to quantitatively analyze the role of cell polarization for the different modes of locomotion (see Figure 3). For every row in Figure 3, the panel in the first column shows a snapshot of a simulated cell with a superposition of the polarization vector P (in black) and the cell propagation vector V (in red). The second column presents the cell trajectories in space. Here, a pronounced change from random to persistent motion can be seen from top to bottom with growing k _a , in line with our earlier results [40]. Thus, cells perform larger explorations in space as they become more persistent for high values of the parameter k _a . In the third column, the correlation of the normalized magnitude of the vectors P and V is displayed, revealing a high degree of correlation for large values of k _a , whereas correlations are low for small k _a . In the fourth column, circular histograms of the angle θ between the vectors P and V are displayed. For small values of k _a both vectors are typically misaligned, while alignment is observed for high values of k _a . The corresponding probability distribution exhibits a peak around zero that becomes increasingly pronounced for growing values of k _a . These results show that the vectors P and V are more likely to align for growing values of k _a , resulting in an increased probability for values of θ close to zero.

In short, the dynamics of the cell propagation vector V, the polarization vector P and the angle θ between them reflects the more persistent and less random motion of cells as the parameter value of k _a increases. To summarize the dependencies on the parameter k _a , we show the mean of the absolute values of the vectors P and V as well as the angle θ between them as a function of k _a in Figure 4. With increasing parameter k _a , larger mean values of the polarization vector P are observed in Figure 4A. This is associated with increasing mean values of the velocity of polarized cells, see the magnitudes of the vector V in Figure 4B. Figure 4C finally shows that not only the absolute values of P and V increase, but also their alignment is more pronounced with increasing k _a , so that the mean angle θ decreases. For a better comparison, box plots of the angle θ are presented in Figure 4D. In summary, persistent motion (high values of k _a ) is characterized by larger magnitudes of the vectors P and V and by smaller angles θ between them, i.e. by increased alignment of the polarity and the displacement directions. On the other hand, random motion (small values of k _a ) results from smaller magnitudes of P and V and larger angles θ, indicating the absence of alignment due to the irregular distribution of patches of c _i that trigger random protrusion all around the cell boundary.

Before we address the dynamics of multiple interacting cells, we first focused on the interactions between two cells. As detailed in Section 2.1, we extended the phase field description of a single amoeboid cell by including a repulsive force between them; any kind of adhesive force is neglected. For the following simulations, we have chosen the model parameters in the regime of persistent motion, in particular, the reaction rate was set to k _a = 5s ⁻¹. We modeled a pair of cells close to each other on a square grid for a period of 90 s with the purpose of making them interact and analyzing the collision dynamics between them. From our simulations, we could distinguish three types of collision scenarios that we termed alignment, anti-alignment and stuck/push, illustrated in Figure 2 (see also Section 2.5). In Figure 5, we display representative series of snapshots for each of the three scenarios observed in simulations of glancing collisions (Figures 5A–C) as well as for head-on collisions (Figures 5D–F). We found that the avoidance between cells (anti-alignment) was the most frequently observed case in our simulations, see Figures 5B,E. The second most frequent scenario was the alignment of cells, as can be seen in Figures 5A,D. Finally, the stuck/push scenario was only rarely observed, see Figures 5C,F. The two cells, which get stuck upon collision, push against each other for an extended period of time, before the shape of one of them is strongly distorted, as a consequence of which the heads-on pushing configuration is destabilized and they continue moving in different directions. However, as we will show in the next sections, this case will be more frequently observed when the density of cells is increased.

We furthermore analyzed the collision-induced dynamics of the cells in more detail in Figure 6, following the same order as in Figure 5. In the first column, we display snapshots of the two interacting cells for every studied case, including the polarity vector P and velocity vector V for each cell. In the second column, representations of the cell trajectories are shown, revealing clear differences between the different collision scenarios. In the alignment case (A and D), for example, trajectories tend to be parallel, while a crossing of tracks is seen for the anti-alignment scenario (B and E). Only small displacements are obtained in the stuck/push case (C and F) as cells impede each others motion upon head-on collision. The third column displays again the correlation of the normalized magnitudes of the vectors P and V. The fourth column shows the distributions of the angle θ between the polarity and velocity vectors P and V for the two interacting cells—wider distributions for the angle are seen in the stuck/push scenario.

For all the cases, the magnitude of both vectors tend to correlate, regardless of the orientation of the cells. This behavior is more clearly observed in the alignment and anti-alignment cases, while in the stuck/push case this tendency is weaker. As mentioned before, these distributions are wider for the stuck/push case, indicating stronger fluctuations and less correlations in the orientations of the polarity and velocity vectors P and V. The results discussed above are summarized by box plots of the angle θ between P and V for the different scenarios shown in Figure 7. The alignment and anti-alignment scenarios generate only small differences in the orientations of the vectors P and V during the interactions and, therefore, P and V are typically aligned. In contrast, larger deviations of the relative orientation of the vectors and, consequently, larger values of θ occur in the stuck/push scenario, indicating that vectors P and V are not aligned.

We have furthermore used the cross correlation as a tool to show if a relationship between velocity and polarization vectors among the cells can be observed in order to clarify the results described above. However, the cross-correlation of the vectors P ₁ and P ₂, where the index indicates the different cells, has not revealed any conclusive information because of the different time scales of the repolarization, and the cross-correlation of vectors V ₁ and V ₂ only shows significant differences for the stuck/push conditions. The differences of collision types are most evident from the Supplementary Movies SM1–SM6 of the head-on and glancing binary collisions.

We now move to larger numbers of interacting cells. From the various cell-cell interactions within the ensemble of cells, we were able to identify the previously observed three collision scenarios: alignment, anti-alignment and stuck/push, see Supplementary Figure S1.

We varied the system size L, keeping the number of cells N = 25 fixed, to assess the density dependence of motility characteristics, the relevance of collision scenarios and the resulting collective pattern formation. The packing fraction ρ _R = NA ₀/L ², i.e. the relative fraction of surface area covered by cells also referred to confluency, is a central control parameter which determines the density of cells in the simulation box. At high cell densities, jamming is expected to be more relevant as the probability of cell-cell collisions increases. To study these effects, we simulated five different domain sizes; Figure 8 shows exemplarily the gradual transition from a dilute (Figure 8A and Supplementary Movie SM7) to a dense system (Figure 8E and Supplementary Movie SM8).

In high density scenarios, stuck/push interactions are more commonly observed as cells trying to follow their own trajectory collide with others, thereby competing for free space. This behavior is expected due to the higher packing fraction ρ _R . For the chosen model parameters, cells tend to remain fairly polarized—the mean magnitude of the polarity vector does only weakly depend on ρ _R (cf. Figure 8F). The mean velocity, however, decreases significantly as the cell density is increased (see Figure 8G) due to the lack of free space in a dense system. As a result of frequent cell-cell collisions, the angle between polarity and displacement vectors increases on average as a function of the packing fraction ρ _R as shown in Figure 8H. This is in line with our previous observation that velocity V and cell polarity vector P tend to be misaligned as a result of stuck/push interactions (see Figure 6F and the corresponding discussion in Section 3.2). As cells remain polarized—most of the biochemical component c _i is concentrated in one certain part of cells—but collisions change their velocities, the vectors P and V cease to be aligned and the angle between them, thus, increases.

The numerical results suggest that stuck/push collision scenarios are more relevant at higher cell densities. As cells impede each others motion as a result of stalling, these collisions may lead to the formation of clusters composed of immobile cells. This mechanism is reminiscent of the explanation of motility-induced phase separation (MIPS) as observed in self-propelled discs [23] that show phase separation at high density because particles hamper each others motion upon head-on collision. We want to highlight, however, that the clustering dynamics observed in deformable cells is different from classical MIPS: as cells are deformable, the local stress acting on one cell is anisotropic. This effect is even more pronounced as clusters grow in size and, consequently, the local pressure increases. As a result, clusters frequently break apart, thereby giving rise to a complex and dynamic clustering dynamics. To quantify clustering, we measured the cluster size distribution numerically (Figure 9). We consider two cells to belong to the same cluster if the distance between them is smaller than 12 μm which corresponds roughly to twice the average radius of a cell (in the absence of interactions). At low densities, a few collisions lead to the transient formation of small groups of cells. Group sizes range from a few cells to more than half of the total number of cells in the system (Figure 9A). In contrast, larger groups are more frequently formed in dense systems (Figures 9B–D). Accordingly, the cluster size distribution broadens as the packing fraction is increased and adopts even a bimodal shape in the high density limit (Figure 9E). A direct comparison of the transition from a unimodal to a bimodal shape can be seen in Figure 9F, where panels corresponding to the lowest and highest packing fraction are shown together for a better visualization. This structural change of the cluster size distribution eventually allows to define a critical packing fraction above which clustering sets in, similar to the clustering transition observed in ensembles of self-propelled rods [20, 61].

We furthermore quantify the random transport of cells within an ensemble by the mean square displacement ⟨ Δ x ( t ) 2 ⟩ , i.e. the average displacement of a cell evaluated at different time lags. The mean square displacement decreases monotonically with the packing fraction ρ _R as cells can move more persistently at lower cell densities. In all cases, we observe a transition from a ballistic regime ( ⟨ Δ x ( t ) 2 ⟩ ∼t ²) at short time scales to a diffusive regime ( ⟨ Δ x ( t ) 2 ⟩ ∼t) in the long-time limit (see Figure 10). We measured the density dependence of the effective diffusion coefficient D of cells by fitting Fürth’s formula, x t − x t = 0 2 = 4 D t − τ 1 − e − t / τ , (5) to the numerically obtained curves. The diffusion coefficient decreases as the packing fraction is increased (cf. Figure 10C). Moreover, the fits show that the crossover timescale τ decreases as the density of particle increases, due to a reduction of the mean-free path and a higher collision frequency at high packing fraction.

We double-checked that equivalent results are obtained when the packing fraction is increased by changing the particle number in a system of fixed size of L = 108 μm (data not shown) instead of fixing the particle number and decreasing the side length of the simulation box as discussed above (Figures 8–10).

We conclude the results section by assessing the relevance of the finite system size in numerical simulations. We performed a finite-size scaling analysis, varying the number of cells and system size such that the density or, equivalently, the packing fraction is kept constant. Five representative snapshots of the simulated systems are shown in Figures 11A–E, where the number of cells is 25, 36, 49, 64 and 81, respectively, and the system size was adjusted from L = 60 μm to L = 108 μm accordingly. As shown in Figures 11F,G, the measured value of the mean polarity and velocity does not reveal a significant system size dependence. Furthermore, the angle θ between the vectors P and V in all cases does not change on average as the number of cells is varied (cf. Figure 11H).

The tendency of cells to align or impede each others motion upon collision in stuck/push configurations implies a tendency towards cluster formation in the system (see Supplementary Movies SM9, SM10 for simulations with 64 and 49 cells, respectively). As the system size increases, larger clusters of cells tend to form as shown in Figure 11. We quantitatively investigated the clustering dynamics by measuring the cluster size distribution in the stationary state; a comparison of all the studied cases is displayed in Figure 12. In contrast to the numerical experiment discussed in Section 3.3, where the density was increased, the shape of the cluster size distribution is now qualitatively independent of the system size as the cell number and system size is increased simultaneously such that the packing fraction ρ _R remains constant. The cluster size distribution turns out to be bimodal for the considered packing fraction (Figures 12A–E). Figure 12F shows the probability of the rescaled cluster size for the system with the highest (N = 81) and smallest number of cells (N = 25). Note, however that the results display large fluctuations and conclusions have to be drawn with caution.