To investigate the migratory behavior of the TNBC cell lines HCC38 and Hs578T, we plated these lines in triplicate at 4 different cell densities within 24-well plates (20,000, 50,000, 100,000, and 150,000 cells per well). Subsequently, we performed a Random Cell Migration (RCM) assay (van Roosmalen et al., 2011) using DIC and fluorescence microscopy of Hoechst-stained cells, imaged every 11–13 min for approximately 15 h (Supplementary Figure S1 and Supplementary Video S1). To quantify the migratory behavior of cells, we performed automated cell tracking (Figure 1A, see Section 4.3.1 for details) by first segmenting the nuclei using Watershed Masked Clustering (WMC) (Yan and Verbeek, 2012) and then tracking the segmented nuclei in CellProfiler (Carpenter et al., 2006) using overlap tracking. Because of vignetting following stitching of adjacent images (see Figure 1A DIC + Hoechst) and the high densities of cells in some fields of view (Supplementary Figure S1), some segmentation errors still occurred. Since these can affect the quantification of migration characteristics such as cell speed (Beltman et al., 2009), we compared our automated tracking to manually determined tracks in a subset of wells. Analysis of the two methods of tracking revealed that they resulted in similar estimates for cell speed yet that automated tracking led to slightly lower instantaneous cell speeds compared to manual tracking (Figure 1B). This minor difference could be explained by an overestimation of cell speed due to variability in manual center-of-mass determination (Huth et al., 2010). Therefore, and because overall differences between wells were similar for the two tracking approaches, we continued our analysis using automated tracking.

A particularly striking feature that can be appreciated from the time-lapse videos is that Hs578T cells form “streams” (Supplementary Video S1, lower right; clockwise flow in Figure 1C). To quantify this streaming behavior, we analyzed the migration directions of all cells compared to the directions of all other cells. If cells were to migrate randomly, the average angle between two cell migration directions should approach 90° (Beltman et al., 2009). Consistent with the observation of streams within the videos, close-by cells had a lower average angle between their migration directions than remote cells (Figure 1D). This streaming effect was more pronounced for Hs578T cells than for HCC38 cells and occurred at all densities, although visually, it was mainly apparent at high densities. For both cell lines, but especially for Hs578T, the average angle remained below 90° at all densities, which suggests a preferred direction of migration within wells. We confirmed this finding by polar histograms of cell migration direction (Figure 1E), showing that the migration directions of HCC38 cells were approximately uniformly distributed, whereas Hs578T displayed a clear bias in migration direction. However, such large-scale streams could also be due to stage drift (Beltman et al., 2009). Therefore, we took advantage of having two imaged locations per well (technical replicates), for which it would be expected that the polar plots would look highly similar if the streaming effects were due to stage drift. The individual polar plots of two associated well locations frequently exhibited a directional bias for Hs578T cells, yet this bias was typically different between the locations (Supplementary Figure S2A), strongly suggesting that streaming was not due to stage drift. Besides the presence of large-scale streams, both cell lines exhibited strong local streams, as evident from the strong decrease of migration angles for nearby cell movements compared to remote movements (Figure 1D). Moreover, this difference in angle profiles remained present after correction for large-scale streams by subtracting the net overall displacement from each frame (Supplementary Figure S2B). In conclusion, both HCC38 and Hs578T cells formed local streams at all observed cell densities, and especially for Hs578T these streams occurred at a scale beyond the employed image dimensions.

Visual inspection of images of HCC38 and Hs578T cells (Supplementary Video S1) revealed that at low density HCC38 cells formed clusters (Figure 2A top left), whereas Hs578T cells were less closely apposed to each other, although they may still be in contact via extended pseudopods (Figure 2A top left). At high densities (Figure 2A top right), clustering became less dominant for HCC38 yet remained similar for Hs578T. This clustering is surprising since HCC38 is a basal B cell line which are typically considered mesenchymal-like because of their high Vimentin levels (Figure 2B). Consistent with differential clustering between the two cell lines at low density, Hs578T cells traveled further than HCC38 cells, as visible from tracks with starting points normalized to the origin (Figure 2C). Nevertheless, the adhesion presumably driving HCC38 clustering did not completely prevent them from escaping these clusters, a feature that seems to be mediated by pseudopod-driven attachment (Supplementary Video S1).

To quantify clustering, we employed spatial descriptive statistics by means of the Ripley L function (Ripley, 1977). Ripley L allows objective quantification of clustering compared to fully random placement of objects within a region, which is especially useful at high densities where differences in clustering are hard to determine by eye. Specifically, negative values for the quantity r—L(r) (see Section 4.5 for details) imply clustering, whereas positive values imply preferred dispersion. Beyond the diameter of an average nucleus (approximately 25 μm for Hs578T and 30 μm for HCC38), HCC38 cells clearly exhibited clustering, yet this clustering gradually disappeared at increasing densities (Figure 2D, red). Hs578T cells did not cluster but rather exhibited preferred dispersion (Figure 2D, cyan), suggesting that these cells actively stay away from close apposition. Note that the initial increase in Figure 2D shows short-range dispersion for both cell lines, which is caused by volume exclusion (the small dips in this initial bump could be the result of occasional oversegmentation). Over time HCC38 clusters became less compact, as evident from an upward shift in the Ripley-L curves (Supplementary Figure S3). In conclusion, HCC38 cells formed clusters at low densities, whereas this was not the case for Hs578T cells.

Following our analysis of (collective) migration and clustering, we turned our attention to other aspects of cell migration and whether these depended on cell density and cell type. First, we studied instantaneous speed, for which we investigated the relation with the observed cell densities within wells rather than with the plated densities (Figure 3A), because these might differ due to spatial heterogeneity at different well locations. This speed analysis revealed that HCC38 cells move faster with increasing density (Figure 3A left). In contrast, the speed of Hs578T cells was largely unaffected by cell density (Figure 3A right), i.e., despite substantial experiment-to-experiment variability, there was a similar speed at all densities for each separate experiment.

In addition to cell speed, a short-term measure of migration, we also analyzed the directional persistence of cells, a long-term measure of migration. A commonly used measure of persistence is the confinement ratio [also known as “straightness” (Wortel et al., 2021a) or “meandering index” (Svensson et al., 2018)]. However, since this measure is strongly biased by track duration (Beltman et al., 2009; Gorelik and Gautreau, 2014), it is not suitable for our data which has substantial variation in track duration (Supplementary Figure S5). Instead, we analyzed persistence with Directional Auto Correlation (DAC) (Gorelik and Gautreau, 2014), which represents an unbiased measure of how fast cells lose their direction of migration (see Figure 3A in Gorelik and Gautreau, 2014). Whereas for Hs578T there is no notable difference in the decay of the DAC for different densities, for HCC38 there is a fast decay for the lowest density (Figure 3B). The relation between the DAC and the lag time τ _lag can be characterized by the following exponential decay function: ϕ ⁡ exp − τ lag / τ p .

Here, τ _p is the persistence time and ϕ is the persistent fraction, a measure for the fraction of cells that is persistent (Vedel et al., 2013). Calibration of τ _p and ϕ for the two cell lines (Supplementary Figure S6) allowed us to study the correlations between cell density, speed, and persistence time. However, it should be noted that the optimal parameter fit did not always describe the decrease in the DAC well (see Supplementary Figure S7), so the resulting parameters should be interpreted cautiously, especially ϕ. For Hs578T, an increase in cell density was not associated with an effect on persistence time, but this was the case for HCC38 cells: besides the strong positive correlation between cell density and speed, both cell density and speed also correlated with persistence time (Figure 3C). In conclusion, HCC38 cells strongly increased their speed and persistence for increasing cell densities, whereas the Hs578T cell migration characteristics barely changed for increasing cell densities.

Our observation that Hs578T cells exhibit dispersion rather than clustering or random positioning in space seems consistent with our findings that speed and persistence do not depend on cell density for this cell line: the cells essentially migrate as solitary cells at all densities. This matches results of an earlier computational model designed for 3T3 fibroblast migration by Vedel et al. (2013), which showed constant speed and persistence time with increased density and a decreased persistent fraction due to a high collision rate. However, for HCC38 cells the relation between the observed clustering and the dependence of cell migration properties on density is less obvious. Therefore, we employed computational modeling to find the minimal requirements to achieve the observed HCC38 density effects. A natural framework to model cell migration is the Cellular Potts Model (CPM), which is a grid-based formalism where multiple grid elements constitute a cell, and membrane elements stochastically extend and retract on the basis of a Hamiltonian. In the base CPM (see Section 4.6 for details), cell motion is driven only by adhesion and cell volume requirements, therefore, cells display Brownian motion, i.e., without any preferred direction and/or persistence (Wortel et al., 2021b). To model persistent cell migration realistically, several extensions to the CPM have been proposed: 1) the “basic persistence model” (see Section 4.6.1.1 for details) in which membrane extensions that move a cell in the same direction as its target direction (derived from previous movement directions) are favored (Beltman et al., 2007; Szabó et al., 2010), and 2) the Act model (see Section 4.6.1.2 for details), which provides cells with persistence by modeling local actin dynamics (Niculescu et al., 2015; Wortel et al., 2021b).

We explored a wide range of parameter values for both the basic persistence model and the Act model to investigate whether either of these CPM persistence extensions could reproduce the increases in speed and persistence with density observed in the HCC38 cell line. The basic persistence model does exhibit an increase in cell speed for increasing densities, yet there is no corresponding increase in persistence time (Figure 4A). This finding goes against the universal coupling between speed and persistence that has been observed across many cell types in various in vitro and in vivo settings (Maiuri et al., 2015). Interestingly, the speed increase with density turns into a speed decrease with density when a connectivity constraint (Merks et al., 2006) is added that requires all membrane elements of a cell to remain in touch at all times, i.e., when cell merging is hindered (Supplementary Figure S8 and Supplementary Video S2). This suggests that the increase in speed depends on cells being able to merge and move through each other. The Act-CPM model matches the observed data worse than the basic persistence model, exhibiting a decrease in both speed and persistence with increasing density (Figure 4B and Supplementary Video S3). Analysis of stream formation revealed that both CPM persistence models form local streams where cells align over a maximum of approximately three cell diameters (Figure 4C). In conclusion, although published CPM extensions to describe cellular persistence lead to local streaming, these models are not consistent with the density dependence of HCC38 cell migration.

Pseudopod formation is essential for persistent cell migration (Bergert et al., 2012; Van Haastert, 2011), and we observed high pseudopodial activity in the experimental videos, which seemed instrumental in HCC38 cells being able to move between clusters (Figure 5A and Supplementary Video S1). Cells simulated with the CPM extensions implementing either basic persistence or the Act model do exhibit a non-roundish shape with small extensions, but these are far shorter than the experimentally observed pseudopods, raising the question of whether the pseudopods might play a role in the observed density effects. Therefore, we implemented our previously developed CPM extension in Morpheus that was used to simulate the migration of dendritic-shaped tissue-resident memory T cells (Ariotti et al., 2012). In this model, cells form dendrite-like protrusions in the form of organized actin bundles that extend and retract within the cell, and that move the cell in the direction of the protrusions (Figure 5B; see Section 4.6.1.3 for details). Although this model could achieve the dynamic clustering observed in HCC38 (Supplementary Video S4), it could still not reproduce the experimentally observed dependence of speed and persistence on density (Figure 5C). Similar to the Act model, the average speed decreased for increasing cell densities, presumably because the fixed dendrites obstruct each other’s extensions, thereby hampering (collective) migration, resulting in limited stream formation (Figure 5D).

In order to test whether increased interaction between cellular pseudopods of a cell would matter for the dependence of migration on cell density, we adapted the modeled behavior of pseudopods in three ways (Figure 5E; see Section 4.6.2 for details): First, we added an adhesive bonus to the pseudopod tips, because close observation of the experimental videos suggested that the pseudopods allow cells to attach to and pull on each other. Second, to further stimulate collective migration in which cells promote rather than hamper each other’s migration, we implemented a type of Contact Inhibition of Locomotion (CIL) (Stramer and Mayor, 2017), where protrusions that are not aligned with the current overall movement direction of a cell and are touching other cells are quickly retracted and repolarized. Third, we let the pseudopods exert a pulling force on the cell as a whole in their combined direction (similar to Vedel et al., 2013). These three mechanisms together result in collective migration of clusters (Figure 5F and Supplementary Video S5). Importantly, our extended pseudopod model could explain the experimentally observed speed and persistence time increase with density for HCC38 cells (Figure 5G, cf. Supplementary Figure S6). The simulated collective migration also goes hand-in-hand with cell alignment over whole clusters, which can be appreciated from the streaming quantification (Figure 5H). Note though that there is a direct dependence of the strength of the streaming on cell density in the simulations, whereas this is not the case in the experimental data (cf. Figure 1D). This cell-density dependence is less pronounced at lower surface energies J _cell,med, for which there is also less long-range alignment (Supplementary Figure S9).

Finally, we investigated the relative importance of the three added mechanisms affecting pseudopod dynamics (pseudopod pulling, pseudopod adhesion, and pseudopod touch behavior, including CIL) for the relation between cell density, speed, and persistence (Figure 5E). When varying the pulling strength and the adhesive tip bonus, we found that pulling strength primarily increases density-dependent persistence (Supplementary Figure S10A and Supplementary Video S6), although with a low adhesive tip bonus cells remain stuck in rotating clusters (Supplementary Figure S10A and Supplementary Video S6, top two rows). Additionally, with low pull strength, individual cells cannot overcome the high surface energy between cell and medium, J _cell,med (Supplementary Figure S10A and Supplementary Video S6, second column). Hence, with increasing density, fewer cells are stuck, which results in high correlations. In contrast, the effect of pseudopod adhesion is that it promotes dynamic cell behavior inside clusters, without displaying cluster rotation; however, without a pulling force, these clusters barely move collectively (Supplementary Video S6, bottom left), leading to a negative correlation between cell density and persistence (Supplementary Figure S10A, bottom left). A combination of pseudopod pulling and pseudopod adhesion is needed to obtain the dynamic collective migration with a density-dependent speed and persistence time qualitatively matching HCC38 behavior (Supplementary Figure S10A and Supplementary Video S6, bottom right). Varying the touch behavior of pseudopods reveals that this can fine-tune intercellular pseudopod interaction, but that this is not essential for reproducing the observed density-dependent speed and persistence increase; it merely influences the range of surface energies between cells and medium (J _cell,med) for which the CPM simulations exhibit this behavior (Supplementary Figure S10B).

In conclusion, our extended pseudopod model can explain an increase of speed and persistence time with increasing cell density as we observed for HCC38 cells, where other CPM migration models cannot (Figure 5I). In this CPM extension, the presence of CIL promotes these density effects, but pseudopod pulling and adhesion are essential determinants. Thus, pseudopod interaction between cells is an attractive explanation for the HCC38 migration patterns with density.

Twenty-four hours prior to imaging, HCC38 (ATCC Cat# CRL-2314, RRID:CVCL_1267) and Hs578T (ATCC Cat# HTB-126, RRID:CVCL_0332) cells were seeded in complete medium on 24-well glass bottom plates (Sensoplate, Greiner Bio-One, 662892) coated with collagen (rat tail Type I, 10 μg ml⁻¹), with the layout as shown in Supplementary Table S1. The seeded densities were 20,000, 50,000, 100,000, and 150,000 cells per well, which, assuming uniform distribution in the well, corresponds to approximately 100, 250, 500, and 750 cells/mm². One hour before imaging, live Hoechst was added to the medium, and just before imaging the medium was refreshed (without additional Hoechst). The experiment was performed in triplicate.

To allow nuclear tracking, cells were incubated for 1 h with Hoechst 33342 https://www.sigmaaldrich.com/NL/en/product/sigma/14533. After incubation, the medium was refreshed, and the plate was directly placed on an automated stage of a Nikon Eclipse TI equipped with a fluorescent lamp and ×20 objective [Plan Apo, Air, numerical aperture (NA) 0.75, working distance (WD) 1.0], a Perfect Focus System (PFS) and a temperature- and CO₂-controlled imaging chamber (custom design). Two positions per well were imaged using both fluorescence and DIC microscopy. The plates were imaged at 999 × 999 px (experiment 1) or 948 × 948 px (experiment 2 and 3) using a stitch of 2 × 2 positions with a pixel size of 0.79 µm. The imaging was repeated every 11 min (experiment 1) or 13 min (experiment 2 and 3) for 15 h. Images are available at https://dx.doi.org/10.5281/zenodo.5607734 (Le Dévédec, 2021).

Upon visual inspection of the microscopic images, we noted that for the third experiment of the HCC38 150,000 condition, cells were dying; therefore, we excluded these wells from further analysis.

The imaging processing and analysis consisted of multiple steps. Initially, proprietary Nikon ND2 image files were converted to the Tagged Image File Format (TIFF) using NIS-Elements (NIS-Elements, RRID:SCR_014329).

Tracking data from CellProfiler was imported into R using an in-house developed script (Wink and Burger, 2021; Wink, 2015, Ch. 7) and by fixing track identifiers with the CPTrackR package (Burger et al., 2021a). MTrackJ data were imported using the mdftracks package (Burger, 2021).

Analysis in R (R Project for Statistical Computing, RRID:SCR_001905) (R Core Team, 2018) was performed with RStudio (RStudio, RRID:SCR_000432) (RStudio Team, 2016) and with the packages celltrackR (Wortel et al., 2021a), spatstat (Baddeley et al., 2015), and tidyverse (Wickham, 2017) packages.

To analyze spatial clustering, we used the common transformation on Ripley K (Ripley, 1977) defined as L r = K r π , and provided in the spatstat R package (Baddeley et al., 2015). We subsequently visualized r—L(r) as a function of r such that in case of complete spatial randomness L r = r ⇒ r − L r = 0 , which allows to determine whether clustering (r—L(r) < 0) or dispersion (r—L(r) > 0) occurs.

In the CPM, cells are defined as a collection of lattice sites v ⃗ ∈ Z n with the same cell identifier σ. Each cell also has an associated cell type τ(σ). At sites forming the cell boundaries (referred to as “membrane elements” below), there is a cell-type-dependent surface energy J τ 1 , τ 2 . A simulation consists of a sequence of Monte Carlo Steps (MCS), during which cells attempt to extend membrane elements that would modify the cell identifier of a lattice site σ ( v ⃗ ) into the identifier of one of the neighboring lattice sites σ ( v ⃗ n ) in the 2D Moore neighborhood without the central square (i.e., the 8 sites of the first- and second-order neighbors in 2D). The probability that such an extension is accepted depends on the change in the Hamiltonian: H = ∑ v ⃗ , v ⃗ n J τ σ v ⃗ , τ σ v ⃗ n 1 − δ σ v ⃗ σ v ⃗ n + λ V τ σ ∑ σ v σ − V τ σ 2 , where the first term is the sum of the surface energies over all v ⃗ , v ⃗ n neighbor pairs, and the second term is the elastic area constraint which keeps cells within a range of biologically appropriate sizes. Furthermore, δ is the Kronecker delta, λ V τ is the elastic constant for the area of cell type τ, V _τ is the target area of cell type τ, and v(σ) is the actual area of cell σ.

The probability p that an extension is accepted depends on the change in the Hamiltonian ΔH as follows: p = 1 , for  Δ H ≤ 0 e − Δ H T , for  Δ H > 0 , where T is the temperature (Graner and Glazier, 1992; Glazier and Graner, 1993). We used Morpheus (RRID:SCR_014975, Starruß et al., 2014) for the CPM implementation.