Density-dependent migration characteristics of cancer cells driven by pseudopod coordination

The ability of cancer cells to invade neighboring tissue from primary tumors is an important determinant of metastatic behavior. Quantification of cell migration characteristics such as migration speed and persistence helps to understand the requirements for such invasiveness. One factor that may influence invasion is how local tumor cell density shapes cell migration characteristics, which we here investigate with a combined experimental and computational modeling approach. First, we generated and analyzed time-lapse imaging data on two aggressive Triple-Negative Breast Cancer (TNBC) cell lines, HCC38 and Hs578T, during 2D migration assays at various cell densities. HCC38 cells exhibited a counter-intuitive increase in speed and persistence with increasing density, whereas Hs578T did not exhibit such an increase. Moreover, HCC38 cells exhibited strong cluster formation with active pseudopod-driven migration, especially at low densities, whereas Hs578T cells maintained a dispersed positioning. In order to obtain a mechanistic understanding of the density-dependent cell migration characteristics and cluster formation, we developed realistic spatial simulations using a Cellular Potts Model (CPM) with an explicit description of pseudopod dynamics. Model analysis demonstrated that strong coordination between pseudopods within single cells could explain the experimentally observed increase in speed and persistence with increasing density in HCC38 cells. Thus, the density-dependent migratory behavior could be an emergent property of single-cell characteristics without the need for additional mechanisms. This implies that coordination amongst pseudopods may play a role in the aggressive nature of cancers through mediating dispersal.


Introduction
Breast Cancer (BC) is the most common cancer in women, and one of the main contributors to cancer mortality (WCRF, 2021). The primary cause of cancer mortality is metastasis, yet, because of its complexity, metastasis remains poorly understood (Fares et al., 2020;Suhail et al., 2019). Migration of cancer cells plays a crucial role in the metastatic cascade, not only for the long-range translocation of cells from the primary tumor to potential metastatic sites but also for the short-range dispersal of cells within the tumor, thus allowing accelerated tumor growth (Waclaw et al., 2015;Gallaher, Brown, and Anderson, 2019). A detailed understanding of cancer cell migration is essential to obtain insight into cancer progression and metastasis (Stuelten, Parent, and Montell, 2018), especially since expression of genes associated with migration is strongly associated with breast cancer survival (Nair et al., 2019).
A complicating factor in studying cancer cell migration is that BC is a highly heterogeneous disease. Two methods that are used to subdivide BCs into clinically-relevant subtypes are gene expression profiling and hormone receptor status (Viale, 2012). PAM50, a 50-gene classifier, divides BC into five intrinsic subtypes: luminal A, luminal B, basallike, HER2 (human epidermal growth factor receptor 2)-enriched, and normal-like (Parker et al., 2009;Perou et al., 2000;Sørlie et al., 2001). This classification largely corresponds to classification by estrogen receptor (ER), progesterone receptor (PR), and HER2 status. BCs that are negative for these three receptors are called Triple-Negative Breast Cancers (TNBCs). Claudin-low BC, recently redefined as a BC phenotype rather than intrinsic subtype (Fougner et al., 2020), is characterized by its enrichment for Epithelial-Mesenchymal Transition (EMT) markers and stem cell-like features (Prat, Parker, et al., 2010). Different (breast) cancer cells display a stunning variety in migratory strategies, and various methods have been developed to study these strategies in vitro (Kramer et al., 2013;Pijuan et al., 2019), in vivo (Beerling et al., 2016), and in silico (András Szabó and Merks, 2013;Wu et al., 2014;Huang et al., 2015;Te Boekhorst, Preziosi, and Friedl, 2016). Collective cell migration, where cells migrate in loosely or closely associated clusters, has been extensively studied in morphogenesis (Mayor and Etienne-Manneville, 2016), yet is also highly relevant during cancer metastasis (Rørth, 2009;Friedl et al., 2012). For example, in recent years the existence of intermediate EMT phenotypes is increasingly recognized. Such phenotypes are associated with the collective migration of tumor cell clusters (Brabletz et al., 2018), which can have 23-50 fold increased metastatic potential compared to single cells (Aceto et al., 2014). Despite this attention, many open questions remain regarding the mechanisms at play in such 'social contexts' (Angelini et al., 2011;Vedel et al., 2013). Recently, Jayatilaka et al. (2017) presented experimental evidence that paracrine IL-6/8 signaling amplified by cell density causes fast migration of MDA-MB-231 BC cells. Another approach was taken by Vedel et al. (2013), who studied the social context of 3T3 fibroblast cells using computational modeling, thereby demonstrating how complex collective migratory behavior can be an emergent property of single-cell migration properties. Thus, computational modeling is an invaluable tool to understand experimentally observed cell migration behavior, as hypothesized underlying mechanisms can be studied both at the single and collective level (Te Boekhorst, Preziosi, and Friedl, 2016). Various computational model formalisms for cell migration exist (as reviewed by Van Liedekerke et al., 2015;Te Boekhorst, Preziosi, and Friedl, 2016;Buttenschön and Edelstein-Keshet, 2020). The Cellular Potts Model (CPM) (Graner and Glazier, 1992;Glazier and Graner, 1993) is widely used for this purpose owing to its explicit incorporation of cell shape and flexibility to describe various biomechanical properties (András Szabó and Merks, 2013).
Here we study the migratory behavior of HCC38 and Hs578T, two highly migratory and invasive, claudin-low, basal B TNBC cell lines (Herschkowitz et al., 2007;Prat, Karginova, et al., 2013;Neve et al., 2006;Kao et al., 2009), using time-lapse microscopy and computational modeling with the CPM. To investigate how these cells respond to different social contexts, we plated these cells at various densities and performed 2D cell migration assays using Differential Interference Contrast (DIC) and fluorescence microscopy. In this setting, cell density clearly affected cell migration characteristics such as clustering, speed, and persistence for HCC38 cells, yet not for Hs578T cells. Specifically, at low densities HCC38 cells formed tight clusters, which loosened at high densities; this coincided with increased speed and persistence. We could not reproduce these density effects with published CPM models describing persistent cell migration, yet an extension of a CPM model of pseudopod-driven persistence with strong pseudopod coordination was sufficient to achieve the experimentally observed speed and persistence increase for HCC38 cells. Thus, pseudopodial activity can explain speed and persistence increase with density, provided that the pseudopods of a cell have the ability to affect each other's extension.

HCC38 and Hs578T cell lines both form streams during in vitro imaging
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 (20000, 50000, 100000, and 150000 cells per well). Subsequently, we performed a Random Cell Migration (RCM) assay (Roosmalen et al., 2011) using DIC and fluorescence microscopy of Hoechst-stained cells, imaged every 11-13 minutes for approximately 15 hours (Fig. S1, and Vid. S1). To quantify the migratory behavior of cells we performed automated cell  (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 Fig. 1A DIC + Hoechst) and the high densities of cells in some fields of view (Fig. S1), some segmentation errors still occurred. Since these can affect the quantification of migration characteristics such as cell speed (J. B. Beltman, Marée, and Boer, 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 (Fig. 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' (Vid. S1, lower right; clockwise flow in Fig. 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 degrees (J. B. Beltman, Marée, and Boer, 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 (Fig. 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 degrees at all densities, which suggests a preferred direction of migration within wells. We confirmed this finding by polar histograms of cell migration direction (Fig. 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 (J. B. Beltman, Marée, and Boer, 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 (Fig. 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 (Fig. 1D). Moreover, this difference in angle profiles remained present after correction for large-scale streams by subtracting the net overall displacement from each frame (Fig. 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.

HCC38 cells form dynamic clusters at low densities
Visual inspection of images of HCC38 and Hs578T cells (Vid. S1) revealed that at low density HCC38 cells formed clusters ( Fig. 2A top left), whereas Hs578T cells were less closely apposed to each other, although they may still be in contact via extended pseudopods ( Fig. 2A top left). At high densities ( Fig. 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 (Fig. 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 (Fig. 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 (Vid. 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 Materials & Methods) 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 ( Fig. 2D, red). Hs578T cells did not cluster, but rather exhibited preferred dispersion (Fig. 2D, cyan), suggesting that these cells actively stay away from close apposition. Note that the initial increase in Fig. 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 (Fig. S3). In conclusion, HCC38 cells formed clusters at low densities, whereas this was not the case for Hs578T cells.

HCC38 cells exhibit increasing instantaneous speed and persistence with increasing density
Following our analysis of coordinated 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 ( Fig. 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 (Fig. 3A left). In contrast, the speed of Hs578T cells was largely unaffected by cell density (Fig. 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, Dannenberg, et al., 2019), or 'meandering index' (Svensson et al., 2018)). However, since this measure is strongly biased by track duration (J. B. Beltman, Marée, and Boer, 2009;Gorelik and Gautreau, 2014) it is not suitable for our data which has substantial variation in track duration (Fig. 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 Fig. 3a in Gorelik and Gautreau, 2014). The relation between the DAC and the lag time τ lag can be characterized by the following exponential decay function: 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 (Fig. 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 Fig. 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, matching earlier published results on 3T3 fibroblasts by Vedel et al. (2013). Interestingly, HCC38 cells reacted quite differently to increasing density: besides the strong positive correlation between cell density and speed, both cell density and speed also correlated with persistence time (Fig. 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.

Previous Cellular Potts migration models do not explain observed HCC38 speed-density behavior
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. 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 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 Materials & Methods) cell motion is driven only by adhesion and cell volume requirements, and therefore cells display Brownian motion, i.e. without any preferred direction and/or persistence (Wortel, Niculescu, et al., 2021). To model persistent cell migration realistically several extensions to the CPM have been proposed: 1) the 'basic persistence model' in which membrane extensions that move a cell in the same direction as its target direction (derived from previous moving  Beltman, Marée, Lynch, et al., 2007;Szabó et al., 2010), and 2) the Act model, which provides cells with persistence by modeling local actin dynamics Wortel, Niculescu, et al., 2021). 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 (Fig. 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 mixing is hindered (Fig. S8 and Vid. S2). This suggests that the increase in speed depends on cells being able to mix. 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 (Fig. 4B and Vid. 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 (Fig. 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.

Coordinated pseudopod dynamics can explain density-dependent migratory behavior of HCC38 cells
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 ( Fig. 5A and Vid. 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 (Fig. 5B). Although this model could achieve the dynamic clustering observed in HCC38 (Vid. S4), it could still not reproduce the experimentally observed dependence of speed and persistence on density (Fig. 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 (Fig. 5D).
In order to test whether increased coordination between cellular pseudopods of a cell   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 coordinated migration between cells, 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 to 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 ( Fig. 5F and Vid. S5). Importantly, our extended pseudopod model could explain the experimentally observed speed and persistence time increase with density for HCC38 cells (Fig. 5G, cf. Fig. S6). The simulated collective migration also goes hand-in-hand with cell alignment over whole clusters, which can be appreciated from the streaming quantification (Fig. 5H). Note though that in the simulations there is a direct dependence of the strength of the streaming on cell density, whereas in the experimental data this is not the case (cf. Fig. 1D). This cell-density dependence is less pronounced at lower surface energies J cell,med , for which there is also less long-range alignment (Fig. S9). 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 (Fig. 5I). Strong coordination of pseudopod behavior within cells is thus an attractive explanation for the altered migration patterns with density.

Discussion & Conclusion
TNBC is an aggressive subtype of breast cancer for which targeted therapies are just recently showing some promise (K. E. McCann, Hurvitz, and McAndrew, 2019). Since migration plays a crucial role in the metastatic cascade, more insight into the mechanisms behind TNBC migratory behavior could help identify potential targets for therapeutic intervention. Here we have used a combination of time-lapse microscopy and computational modeling to unravel the migratory behavior of HCC38 and Hs578T, two highly migratory TNBC cell lines. Both cell lines formed streams in our in vitro set-up, yet this was most clear from visual inspection in Hs578T cells. HCC38 cells formed dynamic clusters at low density which became less cohesive at high densities. Furthermore, HCC38 cells exhibited an increase in both speed and persistence time with increasing density. We could not reproduce this density dependence with CPM simulations implementing previously published persistence models, but a pseudopod-driven persistence model with strong coordination between the pseudopods of a cell could reproduce the key features of the experimentally observed HCC38 migratory behavior. Given that HCC38 is a basal B TNBC cell line with very high Vimentin expression (Fig. 2B), one would expect that HCC38 is a mesenchymal cell-line (mentioned as such by Hollestelle et al., 2010;Kim et al., 2019). Thus, it was surprising that HCC38 cells strongly clustered, which is typically indicative of an intermediate EMT phenotype (Bocci, Kumar Jolly, and Onuchic, 2019). A possible explanation is that HCC38 is composed of epithelial and mesenchymal cells at a fixed ratio (as reported by Yamamoto et al., 2017). However, we could not identify subpopulations in our images, nor was this obvious in our single-cell migration analysis. Two indications that HCC38 is in fact a hybrid epithelial/mesenchymal cell line are that HCC38 has (1) high P-Cadherin expression (Kao et al., 2009), indicative of an intermediate EMT phenotype (Ribeiro and Paredes, 2014), and (2) high EpCAM (Epithelial Cell Adhesion Marker) expression (Klijn et al., 2015;Koedoot et al., 2021). Especially the increased EpCAM seems relevant because it has been reported to trigger "the formation of dynamic actin-rich protrusions" (Guillemot et al., 2001). Moreover, following EpCAM overexpression cell interactions are reduced to "sporadic contacts, mainly involving filopodia-like structures" (Litvinov et al., 1997), a description that matches our HCC38 observations (Vid. S1, cf. Fig. 2 in Winter et al. (2003)). This suggests EpCAM could play an important role in shaping pseudopodial interactions between HCC38 cells, and thereby in their migration characteristics.
Computational modeling of pseudopod-driven motility is a long-standing challenge (Schindler et al., 2021), and the incorporation of appropriate pseudopod mechanics in our CPM simulations was not straightforward. For example, based on the experimental images we aimed for long, finger-like extensions; however, for an approximately constant cell volume, such long pseudopods easily pull a cell apart in the CPM. One solution could be to use a compartmentalized CPM, with a separate nucleus and cytoplasm (Scianna and Preziosi, 2021). However, other model formalisms incorporating physical mechanisms in a spatially implicit (e.g., an Agent-Based Model (ABM) as applied in Vedel et al. (2013)), also represent an appropriate way to model coordinated behaviors of pseudopods.
Our finding that HCC38 cells increase their speed and persistence with increasing cell density is somewhat exceptional. Earlier studies have usually reported cell speed to decrease (Angelini et al., 2011) or stay the same (C. P. McCann et al., 2010;Vedel et al., 2013) with increasing density. However, recently it has been shown that in MDA-MB-231, another claudin-low basal B TNBC cell line, paracrine IL-6/8 signaling amplified by cell density does cause faster migration for high than for low densities (Jayatilaka et al., 2017). Other examples of a cell-density-related speed increase include cell motion in endothelial monolayers (Szabó et al., 2010), or confined cell migration (Liu et al., 2015;András Szabó, Melchionda, et al., 2016). The contexts in which these other experiments were executed are somewhat different compared to our experimental setup, in which the speed increase occurred already at quite a low density (Fig. 3A). Nevertheless, it is possible that also in our experimental setting the observed density effects are (partially) due to a density-dependent nutrient gradient or chemotactic/chemokinetic signal. Still, we here showed that such density-dependent migration can occur as an emergent property of pseudopod coordination in single cells without any additional mechanisms. Based on our simulations it seems that cells at low density can get 'stuck' in their respective clusters (Vid. S5), which is similar to the experimental observations (Vid. S1 top left). At high densities in our simulations, the clusters interact more, thereby avoiding rotating clusters, which causes an increase in persistence and speed. Nevertheless, at high densities the differences between simulations and experiments become more pronounced; whereas in the experiments the clusters became less cohesive, the simulations exhibit no difference with respect to cohesion (compare Fig. 1A HCC38 50000 with Fig. 5F  120). This is also reflected in the streams that form during simulations: Within the large migrating clusters that occur at high cell densities, cells' migration directions become aligned over large distances (Fig. 5H). Lowering the surface energy between the cells and the medium J cell,med results in less cohesive clusters, and a shorter-range alignment (Fig. S9). This suggests that cell adhesion might be decreased at high densities compared to low densities, which might also contribute to the high cell speeds at high density (see for example Fig. 5G).
In conclusion, in this study we shed light on the influence of cell density on the migratory behavior of two TNBC cell lines, HCC38 and Hs578T. We could reproduce the experimentally observed density-dependent speed increase in HCC38 cells using a pseudopod-driven CPM with strong coordination amongst pseudopods. A better understanding of the regulatory processes involved in pseudopod formation is urgently needed since they correlate with poor patient survival in multiple cancer types (Jacquemet et al., 2016). Our finding that strong pseudopod coordination can exacerbate the speed and persistence of cancer cells may be a partial explanation for the aggressive nature of such cancers due to high metastatic potential. Additionally, together with a previous report that showed how cell density affects the expression of cell-adhesion molecules (Stanley et al., 1995), the data presented here emphasize the need to include appropriate density-related controls in cell-migration assays.

Cell culture
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 −1 ), with the layout as shown in Table 1. The seeded densities were 20000, 50000, 100000, and 150000 cells per well, which, assuming uniform distribution in the well, corresponds to approximately 100, 250, 500, and 750 cells/mm 2 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.  Table 1: Plate layout for the random cell migration assays, the numbers denote the imaging order. There are two wells per condition, and two positions (technical replicates) per well. Wells were imaged by column in a zig-zag pattern.

Microscopy
To allow nuclear tracking cells were incubated for one hour with Hoechst 33342. 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 20x objective (Plan Apo, Air, numerical aperture (NA) 0.75, working distance (WD) 1.0), a Perfect Focus System (PFS) and a temperature-and CO 2 -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 minutes (experiment 1) or 13 minutes (experiment 2 and 3) for 15 hours. Images are available at 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 150000 condition cells were dying, therefore, we excluded these wells from further analysis.

Image processing
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).

Automated tracking
Subsequently, the resulting TIFFs were processed in a CellProfiler pipeline (CellProfiler Image Analysis Software, V2.1.1, RRID:SCR_007358) (Carpenter et al., 2006), containing the following steps: • Cropping: Following stitching of the images, they contained zero-intensity patches at the edges as a result of (mis)alignment. To avoid problems with segmentation and edge detection later in the pipeline, we cropped the images by 2 pixels at the edges.
• Segmentation: The cropped images were segmented using the WMC approach (Yan and Verbeek, 2012). See Table 2 for the utilized parameters.
• Object identification: We converted the connected components in the segmented images into objects. The resulting objects were filtered on size; we only retained objects with a diameter between 10 px (8 µm) and 40 px (32 µm) for Hs578T, or 50 px (40 µm) for HCC38. Additionally, we discarded objects touching the image border to prevent inaccurate center-of-mass calculation.
• Tracking: We tracked the remaining objects using the Overlap tracking method with a maximal pixel distance of 30.

Manual tracking
To compare our automated tracking to manual tracking, we used MTrackJ (Meijering, Dzyubachyk, and Smal, 2012) in Fiji (Fiji, RRID:SCR_002285) (Schindelin et al., 2012;Rueden et al., 2017) to manually track a representative subset of the wells by clicking the center of mass of each cell in each frame. Although manual tracking is considered the gold standard for tracking (Cordelières et al., 2013), variability in center-of-mass determination (e.g. due to operator fatigue) can cause an overestimation of the actual cell speed (Huth et al., 2010).

Nucleus diameter calculation
Because the cells show high pseudopodal activity the cell diameters are difficult to estimate. Instead, we estimated the nucleus diameters, which are 30 µm and 25 µm for HCC38 and Hs578T. Using the EBImage R package (Pau et al., 2010), we first applied an adaptive threshold on the Hoechst signal, followed by watershed transformation and object feature analysis. Nucleus diameters were estimated as two times the average nuclei radius reported by EBImage, rounded up. These nucleus diameters serve as an approximation for the nearest possible distance between cells.

Directional Autocorrelation
The Directional Auto Correlation (DAC) of all cells was computed by DAC(τ ) = e i j∆t · e i (j+n)∆t , where e i j∆t denotes the normalized direction of motion of cell i at time j∆t, and the angle brackets denote averaging over all cells i and all times j∆t and (j + n)∆t, where ∆t is the sampling time, of which the lag time τ = n∆t is a multiple. We computed the DAC in R using the overallNormDot function, which we contributed to the celltrackR package (Wortel, Dannenberg, et al., 2019).
After removing DAC(0), which is by definition equal to unity, we fitted the exponential decay function φe −τ /τp , which gives an estimate for the weight factor φ (can be interpreted as the fraction of cells that is persistent) and persistence time τ p (Vedel et al., 2013). Since the estimates for φ resulting from parameter calibration were not always reliable (see Fig. S7) we focussed on τ p in our further analysis.

Correlograms
Averaged correlation values for the correlograms were computed using the Fisher transformation. First, the Pearson correlation r for each experiment and simulation replicate was converted into a Fisher's z: where artanh is the inverse hyperbolic tangent. Then z can be averaged and converted back with where tanh is the hyperbolic tangent (Corey, Dunlap, and Burke, 1998).

Clustering analysis with Ripley K
To analyze spatial clustering we used the common transformation on Ripley K (Ripley, 1977) defined as and provided in the spatstat R package (Baddeley, Rubak, and Turner, 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.

Cellular Potts Modeling
In the Cellular Potts Model (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 τ (σ). The energy within cells (i.e., lattice sites having shared σ) is zero, but 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 (MCSs), 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 their 2D Moore neighborhood (i.e., 8 sites). The probability that such an extension is accepted depends on the change in the Hamiltonian: 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 volume constraint which keeps cells within a range of biologically appropriate sizes. Furthermore, δ is the Kronecker delta, λ Vτ is the elastic constant for the volume of cell type τ , V τ is the target volume of cell type τ , and v(σ) is the actual volume of cell σ. The probability p that an extension is accepted depends on the change in the Hamiltonian ∆H as follows: 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.

Existing persistence models
Basic persistence In the basic persistence model implemented in the PersistentMotion plugin in Morpheus, cells have a target direction t based on previous movements which is updated continuously according to where dr = min(1/dt, 1) is the decay rate, with decay time dt in Monte Carlo Step (MCS), and ∆x = x t − x t−1 is the shift of the cell centroid in the previous MCS. For a proposed copy attempt σ(v) → σ(v n ) in update direction s, the additional change in Hamiltonian H due to persistence is computed as: where λ P is the strength of persistence, and v(σ) is the cell volume. For details on the parameters used for the basic persistence model see Table S1.

Act-CPM
In the Act-CPM model, persistence is achieved by recording each lattice site's "actin activity", which depends on the MCSs elapsed since its most recent protrusive activity. Upon a successful copy attempt, the target lattice site is assigned the maximum activity value (Max act ), which decreases every MCS until it reaches zero. By making a copy attempt from an active site into a less active site more favorable, a local positivefeedback mechanism is created which causes persistent motion. For a proposed copy attempt σ(v) → σ(v n ) the additional change in Hamiltonian H is computed as where GM Act (v) is the geometric mean of all activities of lattice sites in the Moore neighborhood of v which share the same cell identifier σ(v) Wortel, Niculescu, et al., 2021). For this study we used the Act-CPM model provided in Morpheus, for details on the parameters used see Table S2.

Pseudopod dynamics
To obtain realistic pseudopod-driven persistence matching our experimental observations, we adapted a published model by (Ariotti et al., 2012). In this model, pseudopod dynamics are realized by extending explicitly described actin fibers from the center-ofmass of the cell in the direction obtained by drawing a value from a von Mises distribution centered around the current movement direction of the cell. This actin fiber can only grow inside the cell, and cell growth directly next to the actin fiber is promoted (via a DeltaH extension), and cell shrinking is prevented. After a maximum growth time the actin fiber is retracted (see SI Ariotti et al., 2012, for details)). See Fig. 5B for an example of this pseudopod-driven motility with a single pseudopod. We implemented this model in a Morpheus plugin, and we adapted it in several ways to implement various processes involved in the coordination of pseudopod behavior: • In HCC38 cells in vitro we observed a 'stickiness' of pseudopods (Vid. S1). To mimic this effect we added an adhesion bonus tip-bonus to the in silico pseudopod tips. The bonus is applied when a proposed cell extension is within the max-distance-for-tip-bonus of one of its pseudopods. Moreover, when this position is also within the max-distance-for-tip-bonus of a pseudopod from a neighboring cell, the bonus is doubled.
• To increase the persistence of the cells, we implemented a pseudopod pulling effect, similar to an effect simulated in Vedel et al. (2013)). Given a copy attempt σ(v) → σ(v n ) in update direction s, the change in Hamiltonian is computed as where F is the pulling force, f is the vector sum of all pseudopods, and v(σ) is the cell volume.
• To increase the alignment of neighboring cells, we implemented a 'touch' strategy inspired by Contact Inhibition of Locomotion (CIL). Specifically, when a pseudopod touches a neighboring cell it is immediately set to the retraction phase allowing pseudopod repolarization. This implementation is executed when the touch-behavior parameter is set to retract. Another touch strategy we implemented includes an additional rule when a pseudopod touches a neighboring cell laterally (i.e. when cos α < .85, with α the angle between the pseudopod direction and the current movement direction of the cell). In that case, the pseudopod is instantly retracted (i.e. skipping the retraction phase), leading to immediate repolarization of the pseudopod. Note that this implementation is executed when the touch-behavior parameter is set to poof-dir. The results presented here are obtained with the poof-dir touch strategy.
See Table S4 for a full list of parameters for the CPM pseudopod extension, and Table S5 for the values used. The code for the Morpheus plugin is available at 10.5281/zenodo.5484491 (Burger and J. Beltman, 2021).

Choice of simulation parameters
To efficiently explore the parameter space, we used the Python Programming Language (RRID:SCR_008394) in Jupyter Notebook (RRID:SCR_018315) (Kluyver et al., 2016), and FitMultiCell (Alamoudi et al., n.d.) (based on pyABC (Schälte et al., 2021) and Morpheus (Starruß et al., 2014)). Based on this extensive exploration, we selected representative parameter sets that qualitatively matched (parts of) the experimental data. The simulated number of cells were similar as the number of cells observed in the experiments (Fig. S4), and all simulations ran for 20000 MCS.

Simulation measurements
We saved the cell positions from the simulations every MCS, and after discarding of the first 1000 mcs to allow for equilibration, we analyzed them in the same way as the experimental cell tracks, except for instantaneous speed, which was estimated based on 50 mcs subtracks.

Parameter Value Description
A 200 × 200 px Simulation lattice size J τ,τ J cell,cell = 1 Surface energies between cell types V τ V cell = 50 px Target volume of cell type λ Vτ λ V cell = 1 Elastic constant for volume P τ P cell = 10 mcs Persistence decay time for cell type λ Pτ λ P cell = 1 Elastic constant for persistence

Parameter Value Description
A 400 × 400 px Simulation lattice size J τ,τ J cell,cell = 1 Surface energies between cell types V τ V cell = 250 px Target volume of cell type λ Vτ λ V cell = 1 Elastic constant for volume S τ S cell = 0.9 Target asphericity of cell type λ Sτ λ S cell = 0.5 Elastic constant for asphericity   Table S4 for a description of these parameters.

Figures
2×10 4 5×10 4 1×10 5 1.5×10 5 HCC38 HS578T Figure S1: Overview of all experimental plated densities for HCC38 and Hs578T cell lines. See Vid. S1 for the corresponding videos.   Figure S7: Example DAC fits for all densities of HCC38. The long-run persistence time is fitted well, but the persistent fraction (which is the intersect with the y-axis) is fitted quite poorly.  Figure S8: Influence of the connectivity constraint on cellular motility simulated with the basic persistence model using the Morpheus PersistentMotion plugin.