Contact inhibition of locomotion generates collective cell migration without chemoattractants in an open domain

Neural crest cells are embryonic stem cells that migrate throughout embryos and, at different target locations, give rise to the formation of a variety of tissues and organs. The directional migration of the neural crest cells is experimentally described using a process referred to as contact inhibition of locomotion, by which cells redirect their movement upon the cell–cell contacts. However, it is unclear how the migration alignment is affected by the motility properties of the cells. Here, we theoretically model the migration alignment as a function of the motility dynamics and interaction of the cells in an open domain with a channel geometry. The results indicate that by increasing the influx rate of the cells into the domain a transition takes place from random movement to an organized collective migration, where the migration alignment is maximized and the migration time is minimized. This phase transition demonstrates that the cells can migrate efficiently over long distances without any external chemoattractant information about the direction of migration just based on local interactions with each other. The analysis of the dependence of this transition on the characteristic properties of cellular motility shows that the cell density determines the coordination of collective migration whether the migration domain is open or closed. In the open domain, this density is determined by a feedback mechanism between the flux and order parameter, which characterises the alignment of collective migration. The model also demonstrates that the coattraction mechanism proposed earlier is not necessary for collective migration and a constant flux of cells moving into the channel is sufficient to produce directed movement over arbitrary long distances.

Figure 1

Migration of cells simulated with the contour-based model. (a), (b) Snapshots of cells in an open migratory domain with a channel geometry and dimensions $L_x$ = 40 and $L_y$ = 10, where top and bottom boundaries are reflective. Flux rates $f$ are 0.02 (a) and 0.2 (b). Other simulation parameters are: $R_0=0.5$, $A=0.5$, $\omega =2$, $V_0=1$, $\gamma =0.5$, $\xi =0.1$, and $\eta =0.1$. For each cell shape, red arrow: the instantaneous cell velocity; black circle: cell body with radius $R_0$; blue contour: cell shape; see Eq. (1). (c) Simulations are run long enough to observe that the number of cells in the channel remains steady. (d) Order parameter $\mathrm{\Phi }$ versus time. For $f=0.2$, $\mathrm{\Phi }=0.83\pm 0.04$. For $f=0.02$, $\mathrm{\Phi }=0.39\pm 0.17$.

Figure 2

Migratory dynamics simulations of the particle model. (a)–(c) Snapshots of cells, modelled as self-propelled particles, in an open channel with dimensions $L_x$ = 400 and $L_y$ = 60, where the top and bottom boundaries are reflective. Flux rates are 2 (a), 4(b), and 20 (c) at time 500. Other simulation parameters are $v_0$ = 2, $\eta$ = 1, and $R$ = 1. Each red arrow: velocity of a cell. Blue trajectory: displacement trajectory of a cell randomly selected at the start of the simulation. (d) Order parameter $\mathrm{\Phi }$ versus time. For $f$ = 2, $\mathrm{\Phi }$ = $0.34\pm 0.09$. For $f$ = 4, $\mathrm{\Phi }$ = $0.48\pm 0.07$. For $f$ = 20, $\mathrm{\Phi }$ = $0.80\pm 0.03$. Parameter values correspond to mean $\pm$ standard deviation (SD). (e) Cell density $\rho$ versus the binned length of the domain, where each bin is five length-unit wide. $\rho$ are $0.19\pm 0.05$, $0.28\pm 0.03$, and $0.84\pm 0.02$ for flux rates $f$ = 2, 4, and 10, respectively. (f) Distribution of the migration time $t_{\text{m}}$ needed for cells cross the domain. For $f$ = 2, $t_{\text{m}}$ = $587.07\pm 219.34$. For $f$ = 4, $t_{\text{m}}$ = $419.76\pm 111.19$. For $f$ = 20, $t_{\text{m}}$ = $251.14\pm 14.94$.

Figure 3

Migratory dynamics of noninteracting particles (i.e., $R=0$) with the maximum noise intensity (i.e., $\eta =2\pi$). (a) A snapshot of noninteracting particles in an open domain with dimensions $L_x$ = 400 and $L_y$ = 60 at time 500. Simulation parameters are $f$ = 20, $v_0$ = 2, and $\eta$ = 1. Each red arrow: velocity of a cell. Blue trajectory: displacement trajectory of a cell randomly selected at the start of the simulations. (b) The order parameter $\mathrm{\Phi }$ is equal to $0.14\pm 0.07$ (for $\eta =1$) and $0.08\pm 0.04$ (for $\eta =2\pi$). (c) Cell density $\rho$ versus the binned length of the domain, where each bin is five length-unit wide. Lines: fit of $\rho \left(x\right)=(L_x-x)m$, where the slope $m$, which is equal to $0.14\pm 0.01$ (for $\eta =1$) and $0.22\pm 0.01$ (for $\eta =2\pi$), corresponds to $f/D$; see Eq. (9). Least-squares fitting was performed using MATLAB (version 2017b, The MathWorks, Inc.). (d) The migration time $t_{\text{m}}$ is equal to $1338.35\pm 847.16$ (for $\eta =1$) and $1702.28\pm 933.48$ (for $\eta =2\pi$).

Figure 4

Migratory dynamics of cells versus the influx rate $f$, obtained using the particle model. Order parameter $\mathrm{\Phi }$ (a), migration time $t_{\text{m}}$ (b), and cell density $\rho$ (c) calculated with model parameters $v_0$ = 1 and $R$ = 1 for $\eta$ = 0.5 and 1.5. Open domain dimensions: $L_x$ = 400 and $L_y$ = 60. Dashed line in (b): minimum migration time $L_x/v_0$.

Figure 5

Persistence length of the migrating cells calculated using the particle model. (a) Four example calculations of the mean square of the end-to-end distance $R$ of the displacement trajectories of particle versus the trajectory length $L$. Symbols: $\langle R^2\rangle$ of trajectories generated from the particle model simulations. Curves: fit of Eq. (10) to the symbols, performed using MATLAB (version 2017b, The MathWorks, Inc.). Inset: the persistence length $P$ values correspond to mean $\pm$ standard error (SE). (b) The persistence length $P$ versus the influx rate $f$ (mean $\pm$ SE). Simulation parameters $v_0$ = 1 and $R$ = 1 for $\eta$ = 0.5 and 1.5 with an open domain of dimensions $L_x$ = 400 and $L_y$ = 60.

Figure 6

Migratory dynamics of cells calculated using the particle model. Order parameter $\mathrm{\Phi }$ (top row), migration time $t_{\text{m}}$ (middle row), and cell density $\rho$ (bottom row) versus cell velocity $v_0$ (a)–(c), noise intensity of the cell motility $\eta$ (d)–(f), and width of the open domain $L_y$ (g)–(i). $L_x$ = 200. In (a)–(f), $L_y$ = 60. In (\textrma)–(c), $\eta$ = 1 and $R$ = 1. Dashed curves in (c) calculated as $f/\left(v_0\mathrm{\Phi }\right)$ normalized to $L_y$. In (d)–(f), $v_0$ = 1 and $R$ = 1. In (g)–(i), $v_0$ = 1, $\eta$ = 1, and $R$ = 1. Dashed curves: mean $\mathrm{\Phi }$, $t_{\text{m}}$, and $\rho$ versus $L_y$ for different flux density values, proportional to $L_y$. Shaded region: SD. For example, with $f=1$ and $L_y=15$, one cell enters the domain every four time steps.

Figure 7

Alignment of collective migration versus the cell density calculated using the particle model in (a) open ($400\times 60$) and (b) closed ($60\times 60$) domains. Simulation parameters are $v_0$ = 1 and $R$ = 1 for $\eta$ = 0.5, and 1.5.

Figure 8

The alignment of collective migration versus the cell density calculated using contour-based model in open ($40\times 10$) and closed $(30\times 30)$ domains. Simulation parameters: $R_0=0.5,A=0.5,\omega =2,V_0=1,\gamma =0.1$, $\xi =0.1$, and $\eta =0.05$.