Spontaneous Spatiotemporal Ordering of Shape Oscillations Enhances Cell Migration

TL;DR

This study uses numerical simulations to show that shape oscillations in cells can spontaneously synchronize and enhance collective migration through purely physical interactions, independent of signaling mechanisms.

Contribution

It introduces a novel physical mechanism where internal shape oscillations promote local synchronization and collective migration without signaling or adhesion.

Findings

Shape oscillations are crucial for local cell rearrangements.

Synchronization of internal clocks enhances collective movement.

Physical interactions alone can induce cell order without signaling.

Abstract

The migration of cells is relevant for processes such as morphogenesis, wound healing, and invasion of cancer cells. In order to move, single cells deform cyclically. However, it is not understood how these shape oscillations influence collective properties. Here we demonstrate, using numerical simulations, that the interplay of directed motion, shape oscillations, and excluded volume enables cells to locally "synchronize" their motion and thus enhance collective migration. Our model captures elongation and contraction of crawling ameboid cells controlled by an internal clock with a fixed period, mimicking the internal cycle of biological cells. We show that shape oscillations are crucial for local rearrangements that induce ordering of internal clocks between neighboring cells even in the absence of signaling and regularization. Our findings reveal a novel, purely physical mechanism through which the internal dynamics of cells influences their collective behavior, which is distinct from well known mechanisms like chemotaxis, cell division, and cell-cell adhesion.

Figures

Equations

$$\{\begin{array}{ll}\xi {\mathbf{v}}_i^{\mathrm{b}}& ={\mathbf{F}}_{\mathrm{mot}}({\mathbf{r}}_i;s_i)+{\mathbf{F}}_i^{\mathrm{b}},\\ \xi {\mathbf{v}}_i^{\mathrm{f}}& =-{\mathbf{F}}_{\mathrm{mot}}({\mathbf{r}}_i;s_i)+{\mathbf{F}}_i^{\mathrm{f}},\end{array}$$

$${\mathbf{v}}_i^{\mathrm{b}}=0,{\mathbf{F}}_{\mathrm{mot}}(\mathbf{r};0)=m\mathbf{r}$$

$${\mathbf{v}}_i^{\mathrm{f}}=0,{\mathbf{F}}_{\mathrm{mot}}(\mathbf{r};1)=0$$

$$r_{\mathrm{ss}}^0=R_{\max }\sqrt{1-2\kappa /m},v_{\mathrm{ss}}^0=r_{\mathrm{ss}}^{0}m/(4\xi ).$$

$${\mathbf{F}}_i^{\mathrm{f}}$$

$${\mathbf{F}}_i^{\mathrm{b}}$$

$${\mathbf{F}}_{\mathrm{Fene}}(\mathbf{r})=\frac{-\kappa }{1-(r/R_{\max }{)}^2}\mathbf{r},$$

$${\mathbf{F}}_{\mathrm{WCA}}^{\alpha \beta }(\mathbf{r})=-24ϵ\left[2{\left(\frac{{\sigma }_{\alpha \beta }}{r}\right)}^{12}+{\left(\frac{{\sigma }_{\alpha \beta }}{r}\right)}^6\right]\frac{\mathbf{r}}{r^2},r<r_{\mathrm{cut}}$$

$$\displaystyle C(x^{\prime},y^{\prime})=\Bigl{<}\sum_{i=1}^{N}\sum_{j=i+1}^{N}$$

$$\displaystyle\times\delta(x^{\prime}_{ij}-x^{\prime})\delta(y^{\prime}_{ij}-y^{\prime})\Bigr{>}$$

$$\displaystyle S(x^{\prime},y^{\prime})=\Bigl{<}\sum_{i=1}^{N}\sum_{j=i+1}^{N}$$

$$\displaystyle\times\delta(x^{\prime}_{ij}-x^{\prime})\delta(y^{\prime}_{ij}-y^{\prime})\Bigr{>}$$

$$P=⟨|\mathbf{e}|⟩,\mathbf{e}=\frac{1}{N}\sum _{i=1}^N\frac{{\mathbf{r}}_i^{\mathrm{bf}}}{\left|{\mathbf{r}}_i^{\mathrm{bf}}\right|},$$

$$v_{\parallel }=⟨\frac{1}{N}\sum _{i=1}^N{\mathbf{v}}_i\cdot \frac{\mathbf{e}}{\left|\mathbf{e}\right|}⟩,$$

$$q=v_{\parallel }\varphi /v_{\mathrm{ss}}^0.$$

$$v_{\perp }=⟨\frac{1}{N}\sum _{i=1}^N\left|{\mathbf{v}}_i-{\mathbf{v}}_i\cdot \frac{\mathbf{e}}{\left|\mathbf{e}\right|}\right|⟩$$

Full text

\corrauthor

[*]thomas.speck@uni-mainz.de \papertypeArticle

Introduction

The collective migration of cells plays an essential role for several biological processes including the formation of embryos, the closure of wounds, and metastasis of tumor cells (1, 2, 3). Motile cells convert available free energy into directed motion, and can thus be regarded as a type of active matter (4, 5). Unlike equilibrium passive systems, active matter can exhibit collective dynamic behavior such as swarming and clustering (6, 7), the study of which has gained much attention from physicists over the last decade. The essence of active systems can be captured with rather minimalistic models, such as the Vicsek model (8) in which isotropic particles interact by aligning their velocity with the average velocity of their neighbors.

Biological cells possess features that are not captured by simple active particles. In fact, cells have an internal structure and a variable shape, they can have different internal states, they can interact through long-range chemical signals, and they can change their phenotype in response to external perturbations. In addition, cells in a multicellular organism differentiate into distinct types during development, leading to different properties corresponding to their biological function. Nevertheless, motile cells can be broadly classified into two types according to their migratory phenotype. Mesenchymal cells migrate slowly and protrude their body into multiple competing lamellipodia (9). Amoeboid cells migrate quickly, cyclically extending their body into a well polarized shape, and have a larger persistence of motion (10, 11). This distinction is not intrinsic as some cells are able to transition between these two migration types, such as tumor cells during invasion (12, 13).

This biological diversity makes it challenging to develop models that are able to reproduce certain phenomena of interest, while remaining general and minimalistic. Several physics-based representations of cells have been proposed in the literature (14, 15). Most prominently, deformable active particles (16, 17, 18, 19, 20), vertex (21, 22) or Voronoi models (23, 24), phase field models (25, 26, 27, 28), active gel models (29, 30), and subcellular element models (31, 32, 33, 34, 35) have been investigated. Among these, subcellular element models offer the advantage of naturally taking into account the internal structure and internal forces, as well as shape deformations within a single cell. This opens up the possibility for studying the effects of such features on large aggregates of cells. To date, studies on collective cell migration have focused primarily on the roles of long-range chemical signaling (36, 37), shape (16, 17, 18, 19, 22, 38), cell-cell (39, 40, 41, 38, 22, 42) and cell-substrate adhesions (43), or cell proliferation (44, 45). To our knowledge, no investigation has addressed directly the role of cyclic shape oscillations on the coherent migration of cells.

We start from the subcellular element model devised in (34), and we focus on modeling the class of motile cells that exhibit ameboid movement. The model features a two-particle representation for each cell, with an active force proportional to its length. Such an active force is able to model Contact Inhibition of Locomotion (CIL), a complex molecular process that was found to play a key role in collective migration of cells (46, 47), for which protrusions are retracted upon exchanging chemical signals during cell-cell contact. Here we explicitly model the crawling dynamics, in which the active force and a contracting force alternate in order to mimic the cycle of expansion and contraction of ameboid cells. Thanks to this explicit crawling, we unveil a general mechanism by which conformational shape changes in ameboid cells give rise to a purely physical mechanism that enhances the collective migration of the whole tissue.

Methods

Model

Cells are represented by two disks, referred to as front and back, see Fig. 1. The dynamics of each disk is modeled to be overdamped, which is a good approximation for a large class of cells moving in a viscous environment and dominated by active forces (48). The equations of motion for the front and back disks of cell $i$ read

[TABLE]

where the forces on the right hand side are decomposed into passive forces $\mathbf{F}_{i}^{\alpha}$ (for details see Interactions) and the cell extensional force dipole $\mathbf{F}_{\mathrm{mot}}$. Here, $\mathbf{v}_{i}^{\alpha}$ is the velocity vector, $\xi$ is the friction coefficient, and $\mathbf{r}_{i}=\mathbf{r}^{\mathrm{f}}_{i}-\mathbf{r}^{\mathrm{b}}_{i}$ is the vector connecting the back disk to the front disk. Finally, $s_{i}$ indicates the stage of cell $i$, which can be either extension ($s_{i}=0$) or contraction ($s_{i}=1$). During the former, the extensional dipole force $\mathbf{F}_{\mathrm{mot}}$ is switched on and the position of the back disk is kept fixed

[TABLE]

with $m$ being a constant parameter called motility that we fix to unity. The active extensional dipole force in Eq. (2) drives the front disk forward along the direction of the cell’s axis, so that the cell can effectively extend its length. During contraction, $\mathbf{F}_{\mathrm{mot}}$ is switched off and the front disk is kept fixed

[TABLE]

so that the internal forces within the cell move the back disk towards the front, resulting in a contraction of the cell. The extension and the contraction stages are performed cyclically over a period $\mathcal{T}$. At the beginning of the simulation each cell $i$ is randomly assigned a starting stage $s_{i}\in\{0,1\}$ and a time $t_{i}\in\{0,\mathcal{T}/2-1\}$ tracking the progression of the cell within its current stage. Such cyclic dynamics induces a crawling motion [Fig. 1], aimed at mimicking the more complex motion of ameboid cells which is regulated by internal cytoskeletal cycles of actin polymerization or blebbing (49). Models of expanding-contracting particles have been studied recently, focusing on the collective motion of passive particles (20) or on the single cell behaviour (35). Our model makes a step further by considering fluctuating volume and active particles, and can be derived as a limiting case of the one described in Ref. (35), where the two elements of the cell are assigned periodically two finite friction coefficients with a generic phase shift.

The active force implemented in Eq. (2) is not constant, but depends linearly on the length of the cell. This choice results in a coupling between shape deformations and migration that is found in biological cells (50): circular cells receive a weaker active force, thus move slower, while elongated cells move faster. Since a cell in a dense environment will have a shorter length, the reduction of the active force in Eq. (2) also models the inhibition of locomotion real cells experience when coming into contact with surrounding cells. This mechanism is often referred to as Contact Inhibition of Locomotion (CIL), and results from a complex interaction that involves both mechanical and chemical signals among neighboring cells (47). In our model, the details of CIL are enclosed in the shape of the dipole force $\mathbf{F}_{\mathrm{mot}}$ described above. As discussed in Ref.(34), this is sufficient to correctly reproduce the effects of CIL on the collective migration, which has been investigated in several works (51, 52, 53).

In Ref. (34) the model is studied in the limit of $\mathcal{T}\to 0$ where no cyclic dynamics is performed, the effective friction coefficients are $2\xi$, and the motility force is applied to the front disk exclusively. Given $R_{\mathrm{max}}$ the maximum physical extension of a cell, and $\kappa$ the parameter controlling its stiffness, it is possible to derive in the limit of $\mathcal{T}\to 0$ the steady-state cell’s length and speed

[TABLE]

From these a migration time scale can be defined as the time it takes for a cell to travel a distance equal to its maximum length, $\tau=R_{\mathrm{max}}/v^{0}_{\mathrm{ss}}$. For more information on how these quantities change upon changing $\mathcal{T}$, see Supplementary Figure S1 and Supplementary Movie 1.

Interactions

The passive forces in Eq. (1) are

[TABLE]

where $\mathbf{F}_{\mathrm{Fene}}$ is the intracellular force between the front and the back disk, $\mathbf{F}^{\alpha\beta}_{\mathrm{WCA}}$ is the interaction between disk $\alpha$ of cell $i$ with disk $\beta$ of cell $j$, $\mathbf{r}_{i}$ is the vector from the back to the front disk of cell $i$, and $\mathbf{r}^{\alpha\beta}_{i}$ is the vector from disk $\alpha$ of cell $i$ to disk $\beta$ of cell $j$. The second term on the right hand side of Eq. (1) represents the internal forces within the cell. The intracellular force is modeled via the Fene force (54),

[TABLE]

where $\kappa$ represents the elasticity of the cell and $R_{\mathrm{max}}$ is its maximum physical extension. The interaction between disks pertaining to different cells is modeled via the repulsive Weeks-Chandler-Andersen force (55)

[TABLE]

with $\sigma_{\alpha\beta}=(\sigma_{\alpha}+\sigma_{\beta})/2$, and $r_{\mathrm{cut}}=2^{1/6}\sigma_{\alpha\beta}$. The relative size of the front and back disk determines the shape of the cell, which we fix to $\sigma_{\mathrm{b}}/\sigma_{\mathrm{f}}=0.8$.

Simulations

We integrate Eq. (1) employing the Euler-Maruyama integrator, with timestep $4\cdot 10^{-3}$ in Lennard-Jones units, corresponding to approximately $3.75\cdot 10^{-4}\tau$ (with $\tau$ is the migration time scale of the model at $\mathcal{T}=0$). Cells are initialized by positioning the front disk on a cubic lattice, and assigning random positions to the back disk around the front disk, taking care of not overlapping different cells. For the simulations with finite $\mathcal{T}>0$, the initial internal time $t_{i}$ and stage $s_{i}$ are assigned randomly using the pseudo-random number generator of the GNU Scientific Library (GSL), so that cells’ cycles are not synchronized. We use a square box for each simulation, and implemented periodic boundary conditions.

For each simulation, data is collected only after the steady state is reached. Each point in the figures shown is obtained by averaging over 15 independent simulations, excluding those simulations that displayed an exceptionally long time to reach the steady state, because of the formation of long lived metastable vortices, see Supplementary Movie 5.

Correlation functions

In order to quantify local correlation/anticorrelation of the internal time $t_{i}$ and stage $s_{i}$ of a test cell $i$ with its surrounding cells, we compute the respective spatial correlation functions $C(x^{\prime},y^{\prime})$ and $S(x^{\prime},y^{\prime})$ as follows

[TABLE]

and

[TABLE]

where $z_{ij}=z_{j}-z_{i}$, and $x^{\prime}_{j},y^{\prime}_{j}$ are the spatial coordinates of the front disk of cell $j$ in the reference frame of cell $i$, in which the $\mathbf{\hat{y}^{\prime}}$ axis coincides with the cell’s orientation, centered at the front disk. Each cosine in the summation of $C(x^{\prime},y^{\prime})$ can give values up to $+1$ if the cell pair $(i,j)$ has close to no difference in the internal times or a difference close to the period $\mathcal{T}$ (correlated), while giving values close to $-1$ if the difference is approximately $\mathcal{T}/2$ (anticorrelated). From this calculation alone, it is not known whether two cells with (anti)correlated internal times possess the same or the opposite stage: $S(x^{\prime},y^{\prime})$ is useful to distinguish between these two different cases.

Results

Collective migration without crawling

We now focus on the emergence of aligned migration and the formation of velocity waves. We first consider the model for $\mathcal{T}=0$ (no crawling) and show that the collective migration shows similar behavior to that of traffic flow. In the next section, using the model at $\mathcal{T}=0$ as a reference, we study how the introduction of an explicit crawling dynamics through the use of a finite $\mathcal{T}>0$ affects the collective behavior of the cellular tissue.

In order to characterize the collective alignment, we compute the polar order parameter,

[TABLE]

where $N$ is the total number of cells, and $\langle\dots\rangle$ denotes an ensemble average performed by averaging over time and independent runs, each starting from different and randomized initial conditions, see Methods. $P$ can take values between zero (random) and unity (perfect alignment). For non-negligible values of $P$, we also compute the average speed along the instantaneous migration direction $\mathbf{e}$,

[TABLE]

where $\mathbf{v}_{i}=(\mathbf{v}^{f}_{i}+\mathbf{v}^{b}_{i})/2$ is the center of mass velocity of cell $i$.

We compute $P$ and $v_{\parallel}$ for a range of densities. In order to have a dimensionless quantity expressing the density, we use the packing fraction $\phi=NA_{\mathrm{ss}}/L^{2}$, where $A_{\mathrm{ss}}$ is the area of a single cell in the steady state, and $L$ is the linear length of the box. Using this definition, freezing occurs at $\phi\approx 1.15$ in the model with $\mathcal{T}=0$. At high packing fractions ($\phi\geq 0.65$) the flow of the tissue develops velocity waves which resemble the waves observed during traffic congestion, where the wave vector is anti-parallel to the migration direction, see Supplementary Movies 2 and 3. Velocity waves have been observed in experiments of tissue expansion (56) and in the aggregation of Dictostelium discoideum during starvation (57). In the latter, velocity waves are accompanied by waves of chemoattractants that are emitted by the cells during starvation.

In our model, velocity waves are purely related to cell-cell interaction and shape deformation upon collision. More generally, the formation of traveling bands in systems of active particles has been observed for Vicksek-like particles (58), soft deformable active particles (59), and colloidal rollers (60), among others. In order to visualize the waves, a convenient plot is the kymograph of the velocities [see Fig. 2(a)]. The kymograph is constructed by averaging the cell’s velocities along the $y$ direction, orthogonal to the progagation direction $x$ of the waves, for different times. As $\phi$ is increased a structure of bands emerges, where layers of slow particles alternate with layers of faster particles. The slope of the bands in the $t$-$x$ plot indicates the direction of the wave, which is opposite to the direction of individual cells, as indicated by the white line in Fig. 2. It is instructive to note that aligned migration and velocity waves arise only when using the motility force that models CIL. The qualitative reason is that a strong dependence of the local velocity on the local density is present, cf. Fig. 2(b). This is critical for the emergence of traffic waves, which arise from a positive feedback between density and velocity (61): if the speed decreases locally, the density will increase locally, which in turn will cause the speed to lower further, and so on.

Clearly, the formation of traffic waves indicates congestion, which is not beneficial for migration. In order to characterize the efficiency quantitatively, we compute the adimensional flux $q$ of cells

[TABLE]

In Fig. 2(c) the flux $q$ is plotted against area fraction $\phi$. The observed behavior is analogous to the one obtained from studies of car traffic (62), where two regimes can be distinguished by the sign of $\partial{q}/\partial{\phi}$: the free flow and the moving jam. In the free flow regime, the flux increases as the density increases and no congestion arises. The flux $q$ reaches a maximum value at a critical packing fraction $\phi_{c}\approx 0.65$, after which velocity waves start to emerge. By increasing the density beyond $\phi_{c}$, cells start to cluster into waiting lines and the flux is reduced (moving jam regime). The shape of the flux-versus-density curve, usually referred to as the fundamental diagram in traffic flow studies, acquires a parabolic shape. This behavior is typical for transport processes dominated by excluded volume and can be understood within the Totally Asymmetric Exclusion Process (TASEP) (63), with the difference that here the cells are not confined to one dimension.

Collective migration with crawling

We now turn our attention to the model where the cyclic dynamics of extension-contraction are explicitly included. By switching $\mathcal{T}$ to finite non-zero values, the dependence of the average migration velocity $v_{\parallel}$ on $\phi$ changes from linear to curved [Fig. 3(a)]. This change is a consequence of the system being able to accommodate more cells given the same space, and with a smaller decrease in migration speed, which implies emerging cooperativity. Velocity waves are now absent and the system is able to maintain high alignment for increasingly higher packing fractions (see Fig. 3(b) and Supplementary Movie 4). The resulting flux of cells is enhanced [Fig. 3(c)], with $\phi_{c}$ shifted to a higher value of $\phi_{c}\simeq 0.80$. At high packing fractions it is even possible to transition from the zero-flux solid to a traveling tissue with finite flux, see Supplementary Figure S2. Note that a similar behavior has been reported in lattice gas models of traffic flow, where the introduction of a stepping kinetic cycle to the TASEP has been found to enhance the flux of molecular motors (64, 65). In those systems, the time spent queuing is used for advancing the internal state, so that molecular motors are immediately ready to move once free space is available, thus enhancing the flux. Although the effect is similar, the microscopic origin in our system is quite different, as we are going to show at the end of this section.

What, then, is the origin of this collective migration in the absence of any explicit communication mechanism? We observe that for finite $\mathcal{T}$, cells tend to spontaneously arrange themselves in such a way that their neighboring cells are either synchronized or anti-synchronized in their expansion/contraction stage, depending on their relative spatial location. This is surprising since the internal clocks of cells run independently. Such spontaneous synchronization facilitates a cell’s movement within its local neighborhood. In order to quantify this phenomenon, we compute two spatial correlation functions: of the internal time $t_{i}+s_{i}\mathcal{T}/2$ and of the internal stage $s_{i}$, both as a function of cell separation (denoted $C$ and $S$, respectively, for details see Methods). We stress that these quantities evolve independently in the simulation since there is no explicit interaction between $t_{i}$ or $s_{i}$ of different cells.

The spatial coordinates $x^{\prime}$ and $y^{\prime}$ are relative to the frame of reference of a tagged cell, where the $\mathbf{\hat{y}^{\prime}}$ axis coincides with its orientation, and the center is placed at its front disk. Fig. 4(a-f) presents the resulting correlations $C(x^{\prime},y^{\prime})$ and $S(x^{\prime},y^{\prime})$ for packing fraction $\phi=1.1$ and for three different cycle durations $\mathcal{T}=\{0.07,0.68,0.90\}\tau$. Blue and green indicate anti-correlation, red and orange indicate correlation, and white indicates no correlation. At $\mathcal{T}\ll\tau$ both $C(x^{\prime},y^{\prime})$ and $S(x^{\prime},y^{\prime})$ show no particular structure except for concentric circles due to the packing of cells. As we increase $\mathcal{T}$ to values close to $\tau$, a pattern develops at short distances from the cell. First, we note that neighboring cells migrating side by side (regions I in Fig. 4(a,d) ) are more likely to be in the same stage with similar internal times. More interestingly, cells directly in front (regions II) are found to be mostly anti-correlated with respect to the internal time, but correlated with respect to stage. Hence, if cell $i$ is at the start (end) of stage $s_{i}$, a cell $j$ directly behind will most likely be at the same stage $s_{j}=s_{i}$ but at the end (start) of it. Cells directly behind the tagged cell (regions III) are anti-correlated both in the internal time and in the stage, meaning that if cell $i$ is found expanding (contracting) then cell $j$ directly in front will most likely be contracting (expanding). One could expect the plots of $C(x^{\prime},y^{\prime})$ and $S(x^{\prime},y^{\prime})$ to be symmetric both in the $y^{\prime}$ and $x^{\prime}$ axis, since when implementing the calculation one treats together the cell pairs $(i,j)$ and $(j,i)$ and $\mathbf{r}_{ij}=-\mathbf{r}_{ji}$. However, this is not the case since one has to first rotate the coordinates to the reference frame ($\mathbf{r}_{ij}\to\mathbf{r}^{\prime}_{ij}$) of each cell.

It is worth noting that, although cells in the system with $\mathcal{T}\ll\tau$ keep their neighbors for very long times (since cells spend a long time queuing), the correlation pattern is zero in constrast to the case $\mathcal{T}=0.9\tau$. In the latter, cells are dynamically changing neighbors all the time, but show remarkably robust ability to maintain the pattern of correlation we observe in Fig. 4. Fig. 4 shows also that a weaker, but nonetheless non-zero, correlation/anti-correlation pattern is present beyond the first coordination shell. Although the structure we observe is short-ranged and does not extend beyond about twice the cell diameter, it is sufficient to globally enhance the migration of the tissue.

What is the microscopic explanation for the enhancement of collective migration and for the spontaneous synchronization of cells, linking the two observations together? First, it is instructive to notice that trajectories of single particles differ qualitatively when increasing $\mathcal{T}$ to $0.9\tau$: instead of queuing and then moving in straight lines when the velocity waves arrive, cells tend to move in a zig-zag fashion, as shown in Fig. 4(g,h). The reason for the zig-zag trajectory is sketched in Fig. 4(i) and shown in Supplementary Movie 4: if a cell during extension finds another extending cell in front, volume exclusion with the back disk allows the first cell to slip past the second, provided $\mathcal{T}$ is large enough and enough space is available. Such a mechanism is reminiscent of lane changes in the traffic flow of cars, and results in a higher probability to find synchronized cells crawling side by side, and asynchronized cells in the front or back. We stress that the adhesion to the substrate is crucial for this lane changing mechanism, since during the crawling cycle cells have their back and front disk’s position fixed sequentially.

We can quantify the effect by computing the cell velocity along the direction orthogonal to the collective migration direction

[TABLE]

Fig. 4(j) shows the results for $v_{\perp}$, the flux $q$, and polarity $P$ for a fixed packing fraction $\phi=1.1$ as a function of $\mathcal{T}$. As expected, $v_{\perp}$ shows an increase in agreement with the increase of flux $q$, which becomes maximal when $\mathcal{T}\approx 0.9\tau$ reaching a value $55\%$ higher than in the case $\mathcal{T}=0$. The plot shows also a steep decay of the flux $q$ and polarity $P$ immediately after $\mathcal{T}\approx\tau$, where the tissue ceases to exhibit coherent motion. In the latter case, the internal cycle of the cell is too long and alignment is disrupted.

Discussion

In this work, we have shown how the introduction of an internal cycle for the movement of single cells strongly influences the collective migration properties of the whole tissue. The internal cycle we use models the ameboid migratory phenotype, which is typical for a wide class of cells and consists of a discrete series of conformational changes of the cell that result in its directed motion. We find that the flux of cells is enhanced and optimized by increasing the duration of the cycle to match a specific value of around $0.9\tau$, where $\tau$ is the migration timescale of the cell in the limit of very short cycles. We have shown that this behavior can be clarified microscopically by the increased ability of cells to perform “lane changes” while crawling, which in turns leads to the emergence of a local correlation/anti-correlation pattern of the cycles of expansion and contraction. The result is an emergent global cooperation among cells in the tissue, through which they achieve better migration properties, which manifests itself through a larger flux without the traffic waves indicating congestion.

Our insights can potentially be exploited to influence and steer cell migration. Through tuning the duration of the internal cycle of cells it becomes possible to either enhance the migration of cells, or, to the contrary, diminish the migration of non-desirable (malignant) cells. Influencing the cycle duration might be accomplished through manipulating the substrate, e.g., periodic stretching of the substrate is able to reorient solitary crawling cells (66, 67, 68, 69, 70, 71, 28). A first step will be to verify our results in experiments that track cells individually and discriminate between expanding and contracting cells, which allows to compute the correlations $C(x^{\prime},y^{\prime})$ and $S(x^{\prime},y^{\prime})$ introduced here. Such investigations will yield insights into the importance of our purely physical mechanisms as compared to other well-known mechanisms like chemotaxis, cell division, and cell-cell adhesion.

Supporting Material

An online supplement to this article can be found by visiting BJ online at link.

Author contributions

MC, TS, and RY designed research. MC, SKS, and JJM worked on the model and simulations, MC performed the simulations and analyzed the data. All authors contributed to interpreting the results and writing the paper.

Acknowledgments

MC is funded by the DFG through the Graduate School “Materials Science in Mainz” (GSC 266) and the collaborative research center TRR 146. MC gratefully acknowledges GSC 266 for funding the research stay in Kyoto during which part of this work was carried out. RY acknowledges the Japan Society for the Promotion of Science (JSPS) KAKENHI (17H01083) grant, and the JSPS bilateral joint research projects. JJM acknowledges the Japan Society for the Promotion of Science (JSPS) Wakate B (17K17825) grant. The calculations were performed using the computational facilities of the University of Kyoto and the supercomputer Mogon at the Johannes Gutenberg University Mainz.

Appendix A Single cell’s cyclic motion

By changing the duration of the contraction-extension cycle, the steady state properties of single cell migration change, see Fig 5 and Supplementary Movie 6. As mentioned in the main text, the steady state velocity of the cell $v^{0}_{\mathrm{ss}}$ decreases, while the collective migration counterintuitively improves, see Main Text and Supplementary Movies 3 and 4. For cycle duration $\mathcal{T}\approx 0.9\tau$, which is where the collective migration is enhanced the most at $\phi=1.1$, $v^{0}_{\mathrm{ss}}$ is $20\%$ lower compared to the case $\mathcal{T}=0$. The averaged (over one cycle) length of a cell at steady state, $r^{0}_{\mathrm{ss}}$, is also decreased upong incresing $\mathcal{T}$. While the average length of the cell is shorter, the amplitude of oscillation is higher. Since the maximum length of the cell $R_{\mathrm{max}}$ is fixed and the steady state speed $v^{0}_{\mathrm{ss}}$ increases, the migration time $\tau=R_{\mathrm{max}}/v^{0}_{\mathrm{ss}}$ increases as well.

Appendix B Transition to a moving tissue

As shown in the main text in Fig. 4(j) for $\phi=1.1$, increasing $\mathcal{T}$ leads to better migration with enhanced flux. We show here that not only the flux $q$ can be increased, but for packing fractions at which the system at $\mathcal{T}=0$ is in the solid phase, $\mathcal{T}$ can drive a phase transition from the static solid into a migrating tissue. In Fig. 6, we show the flux $q$ as function of $\mathcal{T}$ for increasing packing fractions. For all packing fraction studied above $\phi=1.1$, we observe a sharp transition from $q=0$ to a finite value $q>0$. The data for $\phi=1.2$ (represented in green) demonstrates how steep the transition is, even when increasing $\mathcal{T}$ of just 10 simulation time steps. This finding shows how conformational changes in cell tissues at high density are sufficient to trigger collective migration.

Appendix C Metastable vortices

In systems with CIL it is possible to observe long lasting metastable vortices before the occurrence of the steady state collective migration. Such states are characterized by the presence of two vortices of opposite charge for the system sizes we have so far considered (up to $N=10^{4}$), see Supplementary Movie 5. The occurrence of vortices depends solely on the initial conditions, since the model does not employ any random noise, and it is observed more frequently (from 10% to 50% of the samples) in the region of packing fractions $0.65\leq\phi\leq 1.0$. For this reason, we have excluded from the presented data those runs where the vortices where particularly stable to last for the majority of the simulation time.