Formation of motile cell clusters in heterogeneous model tumors:
the role of cell-cell alignment

Abstract

Circulating tumor cell clusters play an important role in the metastatic cascade. These clusters can acquire a migratory and more invasive phenotype, and coordinate their motion to migrate as a collective. Before such clusters can form by collectively detaching from a primary tumor, however, the cluster must first aggregate in the tumor interior. The mechanism of this cluster formation process is still poorly understood. One of the possible ways for cells to cluster is by aligning their direction of motion with their neighboring cells. This work aims to investigate the role of this cell-cell alignment interaction on the formation of motile cell clusters inside the bulk of a tumor using computer simulations. We employ a Cellular Potts model in which we model a two-dimensional heterogeneous confluent layer containing both motile and non-motile cells. Our results indicate that the degree of clustering is governed by two distinct processes: the formation of clusters due to the presence of cell-cell alignment interactions among motile cells, and the suppression of clustering due to the presence of the dynamic cellular environment (comprised of the non-motile cells). We find that the largest motile clusters are formed for intermediate alignment strengths, contrary to what is observed for motile cells in free space (that is, unimpeded by a dense cellular environment), in which case stronger cell-cell alignment always leads to larger clustering. Our findings suggest that the presence of a densely-packed cellular environment and strong cell-cell alignment inhibits the formation of large migratory clusters within the primary tumor, providing physical insight into potential factors at play during the early stages of metastasis.

I Introduction

Cancer metastasis is the process by which cancer cells spread from a primary tumor site to distant organs or tissues in the body. This process is responsible for the majority of cancer-related deaths [1]. The way in which cancer metastasizes depends on a large number of factors, including cell-intrinsic properties such as the invasive potential, and extrinsic properties such as the extracellular environment [2, 3]. There is growing evidence that aggressive metastasis is governed by groups of tumor cells rather than individual tumor cells, thus implying that collective cluster migration is the more effective and perhaps even dominating mechanism [4, 5, 6, 7]. As cancerous cells detach from the primary tumor and start to migrate through the body, they are considered to be circulating tumor cells (CTCs) or CTC clusters [8, 9]; Typically, most CTC clusters consist of 2 to 20 cells, are relatively small compared to the size of the tumor, and have a relatively high metastatic potential [6, 10, 11]. Understanding how these clusters form and migrate is essential for developing new strategies to prevent metastasis from occurring.

One of the current hypotheses is that cluster formation starts at the primary tumor [6, 12] via a change in cell phenotype, such as via the epithelial-to-mesenchymal transition (EMT). During the EMT, epithelial cells acquire a more mesenchymal phenotype which results in a loss of cell-cell adhesion or basal-apical orientation, and enhanced migratory and invasive properties [13, 14]. This transition plays a role in processes such as embryonic development and wound healing, but, more crucially, a dysregulated EMT can promote cancer metastasis [13]. A full or partial EMT can already manifest itself in the primary tumor to generate a heterogeneous tumor tissue [15, 16], and can lead to hybrid EMT states where cells possess both epithelial and mesenchymal characteristics. Research has shown that such hybrid states can further enable the migration of CTC clusters during various steps of the metastatic cascade [17, 14, 11]. The migratory properties of the mesenchymal cells are of particular relevance for cancer progression, as they drive the active migration of cells and could initiate the formation of migrating CTC clusters in the primary tumor.

If cells start to metastasize as clusters from the primary tumor, there must be a process that causes and stabilizes the collective migration of these cells. From a physical perspective, there are various mechanisms that can lead to collective cell migration [18, 19] and therefore also give rise to cluster formation. Examples of cellular clustering mechanisms include differences in cell-cell adhesion [20, 21, 22], differences in velocities [22, 23], self-alignment (or velocity alignment) [24], and cell-cell neighboring alignment [25, 26, 27]. Among these, the latter is a distinctive two-body (cell-cell) interaction that, in addition to single-cell properties, can affect the collective behavior in non-trivial ways. Experimentally it is difficult to determine how neighboring alignment interactions govern cellular clustering, but computer simulations afford a systematic means to investigate the effect of these cell-cell interactions on the emergent clustering dynamics. The neighboring alignment mechanism has already shown good agreement with experimental observations in other biological phenomena [28, 29, 27]. However, to the best of our knowledge, the effect of cell-cell alignment on the clustering of cells in a confluent layer has not been studied yet.

The aim of this work is to identify the effect of cell-cell alignment interactions on (CTC) cluster formation within a model heterogeneous primary tumor. Specifically, we focus on the clustering dynamics of motile mesenchymal cells that reside in a densely-packed tumor environment with many epithelial cells. For simplicity, we model the bulk of the tumor as a two-dimensional (2D) confluent layer, in which the epithelial cells act as a dynamic environment for the migrating cells. As a first approximation to the actual continuous distribution between non-motile and highly motile phenotypes, we consider a binary system with cells that are either completely non-motile or motile. We perform simulations using the Cellular Potts Model (CPM), which is an efficient computational framework to study both collective and single cell dynamics in biologically relevant settings [30, 31, 32, 33, 34, 35]. Indeed, the CPM has already been successfully used in the past to study e.g. invasion of migratory cells into a dynamic non-motile cellular background in the absence of neighboring aligning interactions [24, 23]. To account for cell-cell neighboring alignment among the mesenchymal cells, we draw inspiration from the physics of active matter and introduce a modified Vicsek alignment mechanism [36, 37]. This allows us to study clustering for both strong and weakly aligning cells, thereby covering a wide range of possible cell-cell interaction strengths. In addition to cell-cell interactions, we also vary the persistence time to account for the persistent motion of individual mesenchymal cells [38]. We show that clustering is strongly affected by both the neighboring alignment interactions and the cellular environment through which these cells have to move. Moreover, our results indicate that tuning neighboring interactions has a different effect on clustering than changing the cell-intrinsic persistent motion, even though both relate to the global orientation of motion of the cells. These simulation results improve our understanding of how clusters may form at the onset of cancer metastasis, i.e., in the bulk of a primary tumor.

The remainder of the paper is organized as follows. First, we introduce the Cellular Potts Model, describe the cell-cell alignment mechanism, and discuss the characterization of the clustering. Next, we show how the clustering of migrating cells is affected by the cell-cell alignment interactions and the dense environment. For that, we include the analysis of simulations in free space and in a confluent layer both with and without alignment interactions. Finally, we discuss the interplay between the different effects that simultaneously play a role when including both alignment and the cellular environment in the simulations.

II Materials and methods

II.1 Computational model

To study how migrating clusters form in a densely-packed heterogeneous environment, we employ the CPM as implemented in CompuCell3D [39]. The CPM has originally been formulated by Graner and Glazier to study cell sorting based on differences in cell-cell adhesion [40, 41], but is now widely used to simulate cells and tissue for a broad range of biological phenomena [31, 42, 43, 44]. Briefly, the model represents cells as pixels on a discrete lattice, and the cell dynamics evolves according to the Monte Carlo algorithm [45]. For each trial move, a random pixel (with cell number $\sigma_{i}$) is selected and one of its neighboring pixels (with cell number $\sigma_{j}$) is targeted. To determine whether an pixel-copy attempt changes the target site, the change in the Hamiltonian associated with the pixel-copy is calculated. The probability of copying $\sigma_{i}\rightarrow\sigma_{j}$ is determined by the Metropolis algorithm [46]:

$$P(\sigma_{i}\rightarrow\sigma_{j})=\begin{cases}1&\Delta\mathcal{H}\leq 0\\ \exp\left(-\frac{\Delta\mathcal{H}}{T_{m}}\right)&\Delta\mathcal{H}>0\end{cases}$$ (1)

where $\Delta\mathcal{H}$ is the change in the Hamiltonian due to the pixel-copy and $T_{m}$ is an effective temperature that sets the cell-membrane fluctuations [39]. The timescale of the dynamics is set by a Monte Carlo step (mcs).

Our simulations consist of two cell types, namely passive (non-motile epithelial) cells and active (motile mesenchymal) cells. These two cell types allow us to represent the heterogeneity of the primary tumor, and differentiate between the cell’s motility. The non-motile part of the Hamiltonian $\mathcal{H}_{0}$ for both the active and passive cells is given by [40, 41]

$$\mathcal{H}_{0} =\mathcal{H}_{adh}+\mathcal{H}_{vol}$$

$$=\sum_{\begin{subarray}{c}i,j\\ neighbors\end{subarray}}J_{\alpha_{i},\alpha_{j}}(1-\delta(\sigma_{i},\sigma_{j}))+\sum_{\sigma}\lambda_{V}(V_{\sigma}-V_{t})^{2},$$

where $\sigma_{i}$ is the cell number (at pixel $i$) and $\alpha_{i}$ the corresponding cell type. The cell-cell adhesion term is proportional to $J_{\alpha_{i},\alpha_{j}}$, which depends on the cell types of the adjacent pixels. The target volume of each cell is given by $V_{t}$ and the strength of the volume constraint is $\lambda_{V}$. For the mesenchymal cells, we include an additional active force contribution that accounts for cell motility. The total change in the Hamiltonian associated with a pixel-copy attempt thus becomes [32]

$$\Delta\mathcal{H}=\Delta\mathcal{H}_{0}-\sum_{\begin{subarray}{c}\sigma\end{subarray}}\kappa\hat{p}_{\sigma}\cdot\Delta\vec{R}.$$ (3)

Here $\kappa$ represents the magnitude of the active force, which is set to $0$ for the passive cells. The unit vector $\hat{p}_{\sigma}\equiv(\cos(\theta_{\sigma}),\sin(\theta_{\sigma}))$ represents the direction of the active force for a given cell $\sigma$, and $\Delta\vec{R}$ is the vector associated with the center-of-mass displacement of the cell during a proposed pixel-copy. When the motile cell is moving in the direction of the active force, the active force term decreases the total energy in the system, thereby promoting motion in this direction.

To account for cell-cell alignment among the motile cells, we incorporate an alignment mechanism based on the Vicsek model [36]. Briefly, this model describes the movement of particles that align their direction of motion with nearby particles in the presence of white noise. Although there are various possible mechanisms by which cells can induce such coordinated motion, the Vicsek model offers a computationally efficient means to mimic an effective alignment interaction between cells. Here, we use an adaptation of the original model where the strength of the alignment interaction can be varied, as in Ref. [37]. This adapted model allows us to study a wide range of effective alignment interactions and investigate how it changes the emerging dynamics of cells.

Our adapted Vicsek model is governed by an effective alignment strength and the persistence time characterizing the persistent motion of cells. The evolution of the cell’s preferred direction of motion is determined by

$$\theta_{\sigma}(t+\Delta t)=\arg\left(\hat{p}_{\sigma}(t)+\sum_{\sigma^{\prime}}\gamma\hat{p}_{\sigma^{\prime}}(t)\right)+\sqrt{\frac{2}{\tau}}\Gamma(\Delta t).$$ (4)

The term in large parentheses models the alignment interaction between the active cells, where $\gamma$ is a measure for the alignment strength between cell $\sigma$ and the neighboring cells $\sigma^{\prime}$. We assume that only the active cells align their motion with their active neighbors. The alignment interaction is only included between active cells connecting to each other with at least one pair of pixels. Since the alignment interaction is only included for motile cells, the non-motile cells are not considered to have an explicit role in the (re)alignment of the active cell’s polarity vector. The presence of non-motile cells will affect the cell dynamics indirectly. The last term of Eq. 4 models the rotational diffusion of a cell’s polarization. Here, $\tau$ is the persistence time and $\Gamma(\Delta t)$ is a Gaussian white noise term with zero mean and variance $\Delta t$, i.e., $\langle\Gamma(t)\rangle=0,\,\langle\Gamma(t)\Gamma(t^{\prime})\rangle=\delta(t-t^{\prime})$. Note that $\tau$ is a single-cell property, while $\gamma$ governs the alignment interaction between two neighboring cells.

The alignment strength $\gamma$ can be varied between $0$ and $1$ to study how aligning cell-cell interactions affect the clustering. When setting $\gamma=0$, cells migrate independently of other cells and are essentially acting as Active Brownian Particles (ABPs). On the other hand, when $\gamma=1$, the model reduces to a standard Vicsek formulation where a cell’s direction is the average over itself and its direct neighbors. By tuning the alignment strength from $0$ to $1$ we thus cover a large range of effective cell-cell interactions.

II.2 Simulation details

All cells in the simulation are identical, apart from their motility which is only applied to the active cells. The adhesion energies $J_{\alpha_{i},\alpha_{j}}$ are set to $10$ for all cell types. The clustering dynamics can therefore not be governed by cell sorting as in Ref. [40]. The target volume $V_{t}$ and the strength of the volume constraint $\lambda_{V}$ for the cells are set to $100$ pixels and $2.0$, respectively. The simulation box contains $400$ by $400$ pixels, so there are $1600$ cells included in the simulations. Periodic boundary conditions are applied. We assume that the fraction of motile cells $\phi$ is relatively small, with typically a value of $\phi=0.25$. We also ran simulations for a smaller fraction $\phi=0.10$ and observed similar clustering behavior. These fractions are sufficiently large to allow mesenchymal cells to form clusters but not too large to let these motile clusters dominate the behavior of the primary tumor. The properties that govern the motility of the mesenchymal cells are the active force strength $\kappa$, the persistence time $\tau$, and the alignment strength $\gamma$. We set the active force strength $\kappa$ to $50$. This value is sufficiently high to induce cell migration but not too high as to only induce single cell migration [24]. The persistence time and the alignment strength are varied such that we can studied their effect on cluster formation.

To get an estimate for the typical values of the cells and their dynamics in the simulations, we use information about the dynamics of breast cancer cells in dense tissue [47, 48]. A typical area of such breast cancer cells is $400$ $$\mathrm{\SIUnitSymbolMicro m}$^{2}$ and therefore the size of a pixel in the simulation is $2~{}$\mathrm{\SIUnitSymbolMicro m}$$ by $2~{}$\mathrm{\SIUnitSymbolMicro m}$$. To determine the typical time scale, we can use the characteristic velocity of cells. Based on these velocities (see Supplementary Material for details), the time scale of $1$ mcs represents approximately $7$ seconds.

We initialize the system with non-motile epithelial cells of $10$ by $10$ pixels and randomly change a fraction $\phi$ of those cells to mesenchymal. We start the simulation with the Hamiltonian in Eq. II.1, i.e., without the motility term, to allow the initial configuration to equilibrate. Since the cells are identical apart from the motile term, there is no clustering process happening during this initial phase of the simulation. After $2500$ mcs the cell shapes have equilibrated so the complete Hamiltonian in Eq. 3 is used and we start the analysis. The results are averaged over $200$ independent simulations and the data is stored every $200$ mcs. Since we do not know a priori where the mesenchymal cells are formed and whether they are formed randomly in space or concentrate in a small subdomain of the primary tumor, we also ran simulations with a different initial condition in which all motile cells start as a single strongly aligned cluster. The initial distribution of the mesenchymal cells does not affect the long-time steady state that we observe, see Supplementary Material. The steady state is reached after $60~{}000$ mcs for all parameters studied in this work. We therefore collect the relevant steady state data, indicated with a subscript $\infty$, for all simulations after this point in time.

II.3 Characterizing clustering

We investigate four different scenarios to study the role of the dense environment and the cell-cell alignment interactions on the clustering of migrating cells:

1. Active cells in free space with alignment interactions ($\gamma>0$);

2. Active cells in a confluent layer with alignment interactions ($\gamma>0$);

3. Active cells in free space in absence of alignment ($\gamma=0$);

4. Active cells in a confluent layer in absence of alignment ($\gamma=0$).

For each of these simulations, we characterize the collective motion of cells by measuring the global ordering of the active force and the clustering of motile cells.

We measure the global ordering of the active force with an order parameter $P$ [36],

$$P=\frac{1}{N}\left|\sum_{\sigma=1}^{N}\hat{p}_{\sigma}\right|,$$ (5)

where $N$ is the total number of motile cells and $\hat{p}_{\sigma}$ is the unit vector set by the cell’s orientation angle $\theta_{\sigma}$ (Eq. 4). When the active force of all cells points in random directions the order parameter $P$ is equal to $0$, while its value is equal to $1$ when all polarity vectors point in the same direction. Note that the active force vector is different from the emerging velocity in the simulations, since the net movement of the cells is also affected by shape fluctuations and the presence of other cells. The order parameter serves to provide information about the underlying collective order.

Similarly, we measure the clusters in the simulations for the four distinct situations. Clusters are defined as groups of motile cells, as indicated by the yellow boxes in Fig. 1a, and the number of cells in a cluster is given by $S_{i}$. The mean cluster size $S$ is defined as

$$S=\frac{1}{N_{c}}\sum_{i=1}^{N_{c}}S_{i}=\frac{N}{N_{c}},$$ (6)

where $N$ is the number of active cells in the simulation, $N_{c}$ is the number of clusters, and $S_{i}$ is the number of cells in cluster $i$. Apart from the mean cluster size, we also study the steady-state cluster size distribution and the breaking and merging of clusters in time. The cluster size distribution is defined as $n_{m}m/N$, where $n_{m}$ is the number of clusters of size $m$. Essentially, it defines the probability of taking a cell and finding it in a cluster of size $m$. To calculate the breaking and merging of clusters, we determine the distribution of cells over the clusters at consecutive time steps and determine how this distribution changes. All of these quantities characterize the clustering of the motile cells in the simulations.

Figure 1: a) Snapshot of the simulation for $\gamma=0.01,\tau=2500,\kappa=50$ at $5000$ mcs. Three different clusters are encircled in yellow. b) Mean cluster size as function of time, calculated with Eq. 6. After an initial increase in the mean cluster size, the system evolves towards a steady state.

III Results

III.1 Active cells in free space with alignment

First, we consider the dynamics of active cells in free space to benchmark the clustering behavior without the dense cellular environment. In this case the clustering is determined purely by the alignment interactions and the rotational diffusion. A snapshot of the simulation for $\gamma=0.01$ and $\tau=2500$ mcs at timestep $5000$ mcs is shown in Fig. 1a, where the system of cells is still evolving towards a steady state. The time evolution of the mean cluster size shows how the average number of cells per cluster grows as a function of time (Fig. 1b). The motile cells are initially randomly distributed in space, while most cells belong to the largest cluster in the steady state.

Figure 2: The global order parameter $P_{\infty}$ (a) and the mean cluster size $S_{\infty}$ (b) as a function of the alignment strength $\gamma$ for $\tau=500$ mcs (orange) and $\tau=2500$ mcs (black). The simulations include $400$ motile cells, which corresponds to an active fraction of $0.25$. The inset shows the correlation between $P_{\infty}$ and $S_{\infty}$. Decreasing the alignment strength decreases both the order parameter and the mean cluster size, and thereby leads to less collective behavior.

We measure the global order parameter $P$ upon varying both the persistence time $\tau$ and the alignment strength $\gamma$ (see Fig. 2a). As $\tau$ is increased, we observe a monotonic increase in the global ordering, fully consistent with the standard Vicsek model [36]. The order parameter also increases as we increase the alignment strength. For $\tau=500$ mcs, we find a transition from a highly disordered state where $P$ is close to $0$ to an ordered state as $\gamma$ approaches $1$. Even though both parameters determine the underlying ordering of the active forces of the cells, their origins are rather different: the alignment strength suppresses fluctuations in the active force orientation of the cells, while the persistent time enhances these fluctuations in the system.

Next, we turn to the mean cluster size in steady state. Fig. 2b shows how the mean cluster size increases as we increase either $\gamma$ or $\tau$. When the order parameter is close to $1$, one single cluster is formed. As we decrease the persistence time and the alignment strength, the mean cluster size decreases, indicating that there are more clusters present in the simulations. There is a also correlation between the mean cluster size and the global order parameter (see inset Fig. 2b). These curves do not show a universal scaling behavior, however, implying that the global order parameter does not solely dictate the steady-state mean cluster size, and varying $\gamma$ and $\tau$ independently leads to different clustering behavior.

III.2 Active cells in a confluent layer with alignment

We now consider the more realistic case of a heterogeneous confluent layer with both active and passive cells. Fig. 3a shows a snapshot of the simulation for $\gamma~{}=~{}0.01$ and $\tau=2500$ mcs at timestep $5000$ mcs. The clusters that are formed after $5000$ mcs are different from those in free space with the same set of parameters. The clusters have a different morphology and there are more smaller clusters present. Video S1 shows the dynamic evolution of clusters for $t=0$ mcs to $t=22~{}500$ mcs. The clusters are more strand-like and irregularly shaped, similar to active cell structures observed in Refs. [24, 23, 49]. The cluster morphology is shaped by the cells coordinating their motion with their neighbors but the cell’s migration direction is also restricted by the passive cells in the surroundings. Indeed, two strongly aligned clusters in the two different surroundings studied in this paper evolve very differently: while the active cells in free space migrate collectively (video S2), the cluster in the confluent layer scatters and forms motile-cell strands (video S3). This phenomenon leads to a significantly smaller mean cluster size than in free space over the same time span, and the magnitude of the mean cluster size does not grow as fast (see Fig. 3b). Moreover, there is an overshoot in the mean cluster size before eventually reaching a steady state. The occurrence of such overshoots has been demonstrated as a characteristic of the non-equilibrium nature of biological systems [50]. Even though the cell dynamics of active cells is governed by the same set of parameters as in free space, their dynamic behavior changes when migrating in a dense confluent layer.

Figure 3: a) Snapshot of the simulation for $\gamma=0.01,\tau=2500,\kappa=50$ at $5000$ mcs. Five different clusters are encircles in yellow. b) Mean cluster size $S(t)$ as a function of time, calculated with Eq. 6. After an initial increase, the mean cluster size evolves towards a steady state. The mean cluster size shows an optimum due to the presence of the non-motile cells.

Returning to the steady-state analysis, we observe that the order parameter $P$ shows the same trends upon changing $\gamma$ and $\tau$ as for the simulations in free space (see Fig. 4a). The absolute value of the order parameters is slightly lower, but decreasing either the persistence time or the alignment strength decreases the overall global ordering in the system as was the case for active cells in free space. The small decline in $P$ can be explained by the fact that clusters are smaller and therefore fewer interactions between cells occur. Indeed, when we consider the confluent layer with a smaller number of active cells ($\phi=0.10$), and as a consequence, fewer interactions between active cells, the order parameter also becomes smaller (Fig. 4c). Thus, the passive cells have a minimal effect on the global ordering of the active forces.

Figure 4: The global order parameter $P_{\infty}$ (a) and the mean cluster size $S_{\infty}$ (b) at a fraction of $0.25$ as a function of the alignment strength $\gamma$ for $\tau=500$ mcs (orange, square), $\tau=2500$ mcs (black, circle) and $\tau=4000$ mcs (blue, triangle). The inset shows the correlation between $P_{\infty}$ and $S_{\infty}$. The order parameter shows the same trends upon changing $\gamma,\tau$ as for active cells in absence of passive cells. For the mean cluster size, we see an optimum for intermediate alignment strength. (c, d) The same data for an active fraction $\phi=0.10$.

While the order parameter hardly changes in the confluent layer in steady-state, the mean cluster size shows that the clustering of active cells in a confluent layer is significantly different from that in free space. Fig. 4b shows that the mean cluster size scales non-monotonically with the alignment strength for different values of $\tau$. For small alignment strengths ($\gamma<10^{-2}$), there is a regime where both the mean cluster size $S$ and the order parameter $P$ are relatively small. On the other side of the spectrum where the alignment is strongest ($\gamma$ close to $1$), we note that the mean cluster size is smaller compared to the intermediate regime ($\gamma\approx 10^{-2}$) as well. The inset in Fig. 4b shows that the correlation between the order parameter and the mean cluster size no longer scales monotonically. We have verified that the steady-state values do not depend on the initial condition, see Supplementary Material. Fig. 4d shows that the same trend is also observed when $\phi=0.10$. The presence of the passive cells in the confluent layer strongly affects the mean cluster size and leads to the most significant decrease in the mean cluster size for large alignment strength $\gamma$.

Based on the trends observed for $P$ and $S$, our results show that the clustering of active cells in the confluent layer is different from that in free space even though the global ordering of the active force remains almost unaltered. A similar phenomenon of different clustering phases with similar ordering has been detected in the original Vicsek model recently [51]. In our confluent layer, three aspects (instead of two in free space) determine the emerging clustering dynamics. Apart from the persistence time $\tau$ and the alignment strength $\gamma$, we also need to consider the non-motile cells comprising the dense cellular environment. As discussed before, the persistence time and the alignment strength have a different origin, which lead to distinct effects on the clustering behavior. While the alignment interactions suppresses the breaking of clusters by inducing collective migration, the persistence time promotes breaking of clusters by inducing additional noise. They both affect the collective orientation of the active force and thereby determine the clustering directly. These two effects are in competition with the presence of the passive cells that introduce an additional noise term by manifesting themselves as obstacles when cells are clustering.

We can discern the interplay between the three effects by investigating the results for $\phi=0.10$ and $0.25$. For $\phi=0.25$, the number of active cells is relatively high and frequent interactions between cells increase the alignment. Larger persistence times lead to larger clusters in the low $\gamma$-regime, similar to free space. In the high $\gamma$-regime, on the other hand, the persistence time hardly affects the mean cluster size. The alignment and the presence of the cellular environment must dominate the clustering behavior in this regime. The rotational noise is strongly suppressed while the non-motile cells introduce an additional noise term that affects the emerging clustering dynamics. At a lower fraction of active cells, $\phi=0.10$, the number of interactions between active cells is lowered and thereby the overall alignment effect decreases. Since the persistence time is a single cell property, it is not affected by the fraction of active cells in the confluent layer. Therefore, cells have fewer interactions when $\phi=0.1$ so the alignment is less dominant and the rotational noise has a more noticable effect. We can observe this effect by the increase in the mean cluster size upon increasing $\tau$ over the entire range of $\gamma$. For $\phi=0.10$, the environment again plays the most prominent role for large alignment strength by decreasing the mean cluster size for large $\gamma$. Overall, we observe a competition between the alignment strength and the persistence time, and we find that the dense cellular environment strongly suppresses the formation of clusters where the most notable effect is visible for large alignment strengths.

To get more information about the clusters that are formed in the different regimes of $\gamma$, we characterize the cluster size distribution in Fig. 5. First of all, the cluster size distribution is significantly different for active cells in a confluent layer compared to the case in free space. In the dense cellular environment, single cells are detected more frequently. Interestingly, around the optimum in the mean cluster size ($\gamma=10^{-2}$), the cluster size distribution of Fig. 5 reveals that there is a relative abundance of large clusters (relative with respect to other values of $\gamma$) and a relatively low number of smaller clusters. At this value of $\gamma$, we also see that the probability of forming large clusters initially decreases but peaks for the large clusters. This trend, in which many small clusters appear in tandem with large cell clusters, has also been found in a theoretical model by Peruani and co-workers [52, 53, 54]. In their model, the breaking of clusters is fully determined by single-particle detachment. Although our simulations also contain multicellular cluster detachment, the similarities between our simulations and their theoretical predictions indeed show that single-cell and small-cluster detachments are an important breaking mechanism for cell clusters in the confluent layer.

Figure 5: Cluster size distribution $C(m)$ for $\tau=2500$ mcs and three different values for $\gamma=1,10^{-2},10^{-3}$ in the confluent layer. The distribution in free space for $\gamma=10^{-3}$ is included in gray. The presence of smaller clusters is apparent in the confluent layer. Moreover, the largest clusters form at intermediate alignment strength.

III.3 Clustering without alignment

The results presented so far have shown that both the alignment strength and persistence time, as well as the dense cellular environment, determine how clusters form in the primary tumor. To elucidate the role of the dense environment in absence of these alignment interactions, we set $\gamma=0$ and focus on the clustering dynamics solely due to the environment.

Figure 6a compares the mean cluster size as a function of time for confluent and free-space conditions, starting from a random distribution of active cells with $\gamma=0$. It can be seen that the mean cluster size in the confluent layer remains very small and does not change in time, while in free-space conditions the active cells spontaneously undergo clustering at short times, resulting in a much larger mean cluster size. This initial clustering in free space can be attributed to cell-cell adhesion effects, whereby two cells get close together to decrease the total contact surface with the medium. In the confluent layer, however, the mean cluster size remains identical to the random initial configuration, implying that active cells remain randomly distributed in space when cell-cell alignment interactions are absent. In this case, small clusters are only formed due to random transient encounters. From comparing the case $\gamma=0$ with $\gamma>0$, we can thus conclude that alignment is necessary to form motile cell clusters in a dense cellular environment.

Figure 6: Mean cluster size in time (a), the number and breaking and merging events in time (b) and the cluster size distribution after breaking in the steady state for motile cells in both free space and the confluent layer for $\tau=2500$ mcs and $\gamma=0.0$, i.e. without any alignment interactions. The confluent layer enhances the breaking and merging of clusters, promoting the formation of smaller clusters and thereby decreasing the mean cluster size.

To understand how the dynamics of the clusters evolves in time (in the absence of cell-cell alignment), we plot the number of breaking and merging events in Fig. 6b. The number of breaking and merging events increases in the dense cellular environment compared to free space. Even though the mean cluster size in the confluent layer does not change in steady state, the clusters are highly dynamic. Clusters are constantly changing their composition and dynamically forming and breaking in time. Fig. 6c shows the size distribution of clusters after breaking in steady state. In the confluent layer, small clusters are detected more frequently. Thus, in absence of any alignment, the confluent layer enhances the breaking and merging of clusters, promoting the formation of smaller clusters and thereby decreasing the mean cluster size as compared to the same motile cells in free space.

Finally, we can compare the results for $\gamma=0$ with the results for $\gamma>0$ to hypothesize how the interplay between the cell-cell alignment interactions and the dense cellular environment affects the breaking and merging of clusters during the clustering process. Our simulations have shown that the active clusters in the confluent layer fragment more significantly than in free space when cells interact with the same alignment strength in both surroundings. This breaking of clusters in the confluent layer is inherently different from breaking of clusters in free space where clusters fragment due to a lack of alignment. The non-motile cells in the dynamic cellular environment introduce an additional type of noise in the system. The fragmentation of clusters shows similarities with the disruption of local ordering due to the presence of obstacles. These obstacles can either manifest themselves as a small number of non-aligning agents [55] or as physical (non-motile) obstacles [56, 57]. The non-motile cells in the confluent layer effectively acts as a destructive obstacle that introduces a different cluster breaking mechanism compared to the lack of coordinated motion. Together with the dynamics induced by introducing cell-cell alignment interactions, the clustering of motile cells in a dense tumor environment is a complex process by which multiple breaking mechanisms, affected by multiple sources of noise, simultaneously play a role.

IV Conclusion

In this work, we have investigated the role of cell-cell alignment in the formation of migratory cell clusters in the bulk of a model primary tumor. We simulated the cell dynamics of a small fraction of motile cells in a dense non-motile cellular environment using the CPM. To identify the role of the cell-cell alignment interactions in such a confluent environment, we studied how the coordinated motion of cells, or the lack thereof, determines how active cells in dense cellular environments form clusters. These active clusters could develop inside dense tumors and initiate the formation of metastasis-seeding CTC clusters.

Our computational model has allowed us to study both the effect of cell-cell alignment interactions, and the role of the dense cellular environment, separately. In absence of any alignment interactions, the cellular environment promotes the break-up of small clusters. The non-motile cells act as an impediment to the clustering process. When cell-cell alignment interactions are included, next to the cell-intrinsic persistence time $\tau$, cell dynamics is also governed by a cell-cell alignment interaction strength $\gamma$. Both affect the net ordering of the coordinated motion, but the parameters $\gamma$ and $\tau$ affect the clustering that is observed in the confluent layer differently: while increasing the alignment strength promotes alignment and thereby the formation of clusters, increasing the persistence time suppresses the formation of larger clusters by hindering effective alignment. One of the key findings of our work is an optimum in the mean cluster size upon changing $\gamma$ in the confluent layer. At the optimum, with an alignment strength around $\gamma=10^{-2}$, the largest cell clusters are formed and smaller clusters are detected less frequently compared to other $\gamma$ values. The change in the clustering behavior upon increasing $\gamma$ is due to a transition from a regime where breaking is dominated by a lack of alignment interactions to a regime where breaking is dominated by the scattering of motile cells on non-motile obstacles. The formation of large clusters is strongly suppressed in the confluent layer for large $\gamma$. These results suggest that the formation of actively migrating cells clusters inside the primary tumor is impeded by many distinct effects at play simultaneously. This observation could be an important factor to rationalize the fact that the CTC clusters that are detected in experiments are relatively small.

In this work, we have considered a binary model heterogeneous 2D tumor consisting of epithelial and mesenchymal cells. One important aspect that has only been taken into account partially with this set-up is the fact that the epithelial-to-mesenchymal transition is usually a hybrid transition and primary tumors usually contain more that two types of cells. Many cells express both epithelial and mesenchymal markers and therefore can have both epithelial and mesenchymal properties [17]. When considering the next step in the metastatic cascade, i.e. clusters of cells detaching from the primary tumor, Mukherjee and Levine [58] have already shown using CPM simulations that a relatively small heterogeneous tumor leads to the formation of metastasizing clusters with multiple cell phenotypes. The clusters that detach from the primary tumor consist mostly of motile cells but these cells can pull along less motile cells. The mesenchymal clusters that are studied in our work could thus be the early structures that end up forming the clusters that actually detach from a heterogeneous tumor. Further research is necessary to see how the clusters that are formed inside a primary tumor will initiate the formation of circulating tumor cell clusters.