## Highlights

- •Homogeneous cells appear heterogeneous due to limited sampling and repeatability
- •Heterogeneity bias increases with attraction/repulsion between cells
- •Movement in confined environments decreases apparent heterogeneity
- •Hypothetical applications in neural crest and
*in vitro*cancer systems

## Summary

## Graphical Abstract

## Introduction

*Drosophila*border cells (

*Drosophila*border cell migration (

*in vivo*or in realistic

*in vitro*assays give rise to rich datasets amenable to measuring distributions of statistics of interest. For example, recent efforts by

*in vivo*and

*in vitro*cell-tracking data. We conclude by discussing alternative measures of heterogeneity, further model extensions, and potential avenues for the development of statistical tests.

## Results

### A Generic Model for Collective Cell Migration

*θ*, of cell

*i*at time

*t*+1 is given by

*r*

_{ij}, between cells

*i*and

*j*: (1) alignment with the movement direction of all neighboring cells (within distance

*r*

_{e}, where ${\mathbf{v}}_{j}^{t}$ is the velocity of cell

*j*at time

*t*, for simplicity assumed to be of fixed speed), scaled by parameter

*α*≥ 0; (2) intercellular forces,

**f**

_{ij}, i.e., attraction/repulsion toward/away from neighboring cells (within distance

*r*

_{0}), scaled by parameter

*β*≥ 0; (3) noise, here chosen to be a random vector on the unit sphere, ${\mathbf{u}}_{i}^{t}$, which scales with the number of neighbors,

*N*

_{i}(within distance

*r*

_{0}, including cell

*i*itself), to represent the uncertainty in the forces from neighboring cells (

*η*≥ 0, and without loss of generality we have set

*η*= 1 in the simulations presented here. Interactions between cells are further restricted to nearest Voronoi neighbors only, even if more cells are within the respective distance cut-off. Figure S1 illustrates the various interaction zones.

### Extension to Three Dimensions

*θ*(from −

*π*to

*π*), and the polar angle,

*ϕ*(from 0 to

*π*), which are updated according to

where the interactions in (Equation 2) are again between nearest neighbors only. Cell positions for cell

*i*are updated via ${\mathbf{x}}_{i}^{t+1}={\mathbf{x}}_{i}^{t}+{v}_{0}\left(\text{cos}{\theta}_{i}^{t+1}\text{sin}{\varphi}_{i}^{t+1},\text{sin}{\theta}_{i}^{t+1}\text{sin}{\varphi}_{i}^{t+1},\text{cos}{\varphi}_{i}^{t+1}\right)$, with fixed speed

*v*

_{0}= 0.05, respecting boundary conditions. For details on the form of the intercellular forces, see Method Details in the STAR Methods. In this paper, we have implemented both free (cells are unconfined) and no-flux (reflective) boundary conditions. A similar 3D model has also been developed by

### Intrinsic Heterogeneity

*x*-axis, without loss of generality). This preferred direction could represent, for example, a chemoattractant gradient or other directional cue. The alignment strength,

*α*, is the same as for the other cells, and other cells still align with the informed cell if they are neighboring them. We chose 10% of cells to be in such a leader state, as similarly small fractions have been reported as sufficient to affect the overall population behavior in migrating cell populations (

### Computational Experiments

*L*= 2. To alleviate the influence of transient behavior dependent on initial conditions, we run simulations for 1,000 time-steps and ignore the first 500 time-steps in our analysis, inspecting the order parameter (see Global Model Behavior in the STAR Methods) plotted against time to validate that this is sufficient to reach a (quasi-)steady state. We find that groups of cells move in an ordered (globally aligned) manner for high values of

*α*and

*β*(Figure S3). Example trajectories (Figure 1) reveal moving streams and slowly translating clusters that break apart as the strength of intercellular forces,

*β*, becomes comparable to the noise strength,

*η*= 1, or lower. At low

*α*and high

*β*, groups of cells are in static, positionally ordered, arrangements (Figure 1).

### Delay Correlation Analysis

*i*and

*j*the directional cross-correlation at lag

*τ*is calculated as

where ⋅ denotes the scalar product and ${\langle \dots \rangle}_{t}$ time-averaging. The peak-delay time,

*τ*

_{C}, for a pair of cells is that for which ${C}_{ij}^{\tau}$ is highest in the observed data. By binning the peak-delay times for all cells in the population, we obtain the sample distribution

*P*(

*τ*

_{C}), so that

*P*(

*τ*

_{1}) = 1 would mean that all pairs of cells have their highest directional cross-correlation at lag

*τ*

_{1}, and

*P*(

*τ*

_{0}) = 0 means no pair of cells has peak correlation at lag

*τ*

_{0}(note that

*P*(

*τ*

_{C}) is symmetric about

*τ*= 0 since ${C}_{ij}^{\tau}={C}_{ji}^{-\tau}$). Thus,

*P*(

*τ*

_{C}) constitutes a population-scale measure of heterogeneity.

*r*

_{0}and have directional cross-correlation of

*C*

_{min}= 0.5 or greater at at least one time point. An illustrative example of the resulting distribution of peak-delay times is shown in Figure 2A, indicating that the width of the distribution increases with

*β*. To test whether some cells showed consistently non-zero lag times (either leading ahead or following behind others), we also calculate the distribution of peak-delay times for the directional correlation between individual cells and all other cells in the population. An example of the individual distributions is shown in Figure 2B. For a more comprehensive view, we calculate the SD of the peak-delay distribution,

*σ*(

*τ*

_{C}), as a measure of the width of the distribution, which can be seen to increase consistently with increasing

*β*over a range of values for

*α*(Figure 2C). The SD is also high for low values of

*β*at some values of

*α*, but this trend is not consistent across the range of

*α*values. Increasing the alignment strength,

*α*, for fixed

*β*does not increase the measure of heterogeneity. In summary, stronger attraction/repulsion between cells leads to increased apparent heterogeneity.

**Biological Applications**

### Neural Crest Migration

*Xenopus*cephalic NC a complementary mechanism has been studied, based on balanced contact inhibition of locomotion and co-attraction (

*in vivo*tracking data from trunk and cranial NC cells, we can thus predict the following: If decreased heterogeneity was measured in the trunk versus head, the populations may, in fact, both be homogeneous. This would suggest that leader-like cell states may not exist

*in vivo*, as the observations (differences in heterogeneity) can be explained without them. If the measure of heterogeneity was greater in the trunk, then both locations could host heterogeneous populations. Such an observation would support the notion that cells adopt leader-like states, or otherwise respond differentially to guidance cues, which could be of biological importance.

### Tumor Invasion

## Discussion

*in vitro*, as well as

*in vivo*, experiments and hypothesized their outcomes.

*In vivo*, cell population heterogeneity is often located at the boundary of a migrating collective (

## STAR★Methods

### Contact for Reagent and Resource Sharing

### Method Details

#### Intercellular Forces

**e**

_{ij}between cells

*i*and

*j*, i.e.,

**f**

_{ij}=

*f*(

*r*

_{ij})

**e**

_{ij}. The magnitude of the force,

*f*(

*r*

_{ij}), was chosen to reflect hard core repulsion at short distances (informally described as a force of infinite magnitude), spring-like attraction-repulsion at intermediate distances, and no force at large distances. Specifically, the magnitude of the force between two cells

*i*and

*j*separated by a distance

*r*

_{ij}is given by

Here,

*r*

_{c}denotes the core radius,

*r*

_{e}the equilibrium cell separation,

*r*

_{a}the attraction radius, and

*r*

_{0}the interaction radius, above which cells cease to exert forces on each other. Default parameter values used here are

*r*

_{0}= 0.2,

*r*

_{e}= 0.5,

*r*

_{a}= 0.8, and

*r*

_{0}= 1. In our numerical implementation we approximate the magnitude of infinite force (volume exclusion) by exp(100). The exponential regime was chosen to represent de-adhesion processes for increasing cell-cell separation. This does not accurately reflect forces for cells forming a new contact (decreasing distances), which would require a force-law exhibiting hysteresis. However, we have opted for the simpler form above as our simulations are of cell collectives in contact and we are here not investigating cases when these become cohesive after being initially separated. Thus, the form of intercellular force chosen is sufficient for our purpose. Figure S2 compares our choice of intercellular force with that of

### Quantification and Statistical Analysis

#### Global Model Behaviour

**v**

_{i}for

*N*cells as

where the second equality holds for constant speed

*v*

_{0}. Thus, a value of Φ = 1 means all cells are perfectly aligned, while a value of Φ = 0 means complete disorder (or the special case of perfect anti-alignment, which we do not observe in our simulations). We use this order parameter to globally characterise the model behaviour, and thus restrict our analysis to instances of the model that correspond to the movement of ordered cell collectives only.

#### Comparison with Experimental Data

*N*

_{c1}= 10 and

*N*

_{c2}= 19 cells. We chose to represent these with

*N =*10 and

*N =*20 cells, respectively, in our computational experiments. To keep the cell density consistent, we chose the (random) initial cell positions to lie within a cube of length

*L*= 0.9 and

*L*= 1.1, respectively. Given these dimensions, our cell speed of

*v*

_{0}= 0.05 approximately matches the displacement of cells (2.0–2.4

*μm*per frame, with an initial cluster diameter of ≈50

*μm*, 2.4/50 = 0.048) seen in the datasets analysed. The experimental distributions were reproduced by reading the height of the bars of Figures 4C and D in

_{c1}= 0.92 and 〈Φ〉

_{c2}= 0.84. Hence, we restricted our comparison to simulation outcomes in which the order parameter was within 5% of this range.

### Data and Software Availability

## Author Contributions

## Acknowledgments

## Supplemental Information

- Document S1. Figures S1–S4

## References

- Cellular heterogeneity: do differences make a difference?.
*Cell.*2010; 141: 559-563 Bouchaud, J.P. and Potters, M.. 2009. Financial applications of random matrix theory: a short review. arXiv, 0910.1205v1, [q-fin.ST].

- Complement fragment C3a controls mutual cell attraction during collective cell migration.
*Dev. Cell.*2011; 21: 1026-1037 - Diffusion of individual birds in starling flocks.
*Proc. Biol. Sci.*2013; 280: 20122484 - Tumour-stromal interactions generate emergent persistence in collective cancer cell migration.
*Interf. Focus.*2013; 3: 20130017 - Collective motion of groups of self-propelled particles following interacting leaders.
*Phys. Stat. Mech. Appl.*2017; 479: 467-477 - Moving and staying together without a leader.
*Phys. Nonlinear Phenom.*2003; 181: 157-170 - Effective guidance of collective migration based on differences in cell states.
*Proc. Natl. Acad. Sci. USA.*2012; 109: 2027-2032 - Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells.
*Cell.*2015; 161: 1187-1201 - Visual analytics of delays and interaction in movement data.
*Int. J. Geogr. Inform. Sci.*2016; 8816: 1-26 - Neutrophil trails guide influenza-specific CD8+ T cells in the airways.
*Science.*2015; 349: aaa4352 - Inference of causal information flow in collective animal behavior.
*IEEE Trans. Mol. Biol. Multi-Scale Commun.*2016; 2: 107-116 - Multiscale mechanisms of cell migration during development: theory and experiment.
*Development.*2012; 139: 2935-2944 - Neural crest migration is driven by a few trailblazer cells with a unique molecular signature narrowly confined to the invasive front.
*Development.*2015; 142: 2014-2025 - VEGF signals induce trailblazer cell identity that drives neural crest migration.
*Dev. Biol.*2015; 407: 12-25 - Hierarchical group dynamics in pigeon flocks.
*Nature.*2010; 464: 890-893 Richardson, T.O., Perony, N., Tessone, C.J., Bousquet, C.A.H., Manser, M.B. and Schweitzer, F.. 2013. Dynamical coupling during collective animal motion. arXiv, 1311.1417v1, [q-bio.QM].

- Leader cells define directionality of trunk, but not cranial, neural crest cell migration.
*Cell Rep.*2016; 15: 2076-2088 - Fellow travellers: emergent properties of collective cell migration.
*EMBO Rep.*2012; 13: 984-991 - Multidisciplinary approaches to understanding collective cell migration in developmental biology.
*Open Biol.*2016; 6: 160056 - Collective motion of mammalian cell cohorts in 3D.
*Integr. Biol. (Camb).*2015; 7: 1526-1533 - Coalescent models for developmental biology and the spatio-temporal dynamics of growing tissues.
*J. R. Soc. Interface.*2016; 13: 20160112 - Collective cell migration guided by dynamically maintained gradients.
*Phys. Biol.*2011; 8: 045004 - Phase transition in the collective migration of tissue cells: experiment and model.
*Phys. Rev. E Stat. Nonlin. Soft Matter Phys.*2006; 74: 061908 - Beyond comparisons of means: understanding changes in gene expression at the single-cell level.
*Genome Biol.*2016; 17: 70 - Directional collective cell migration emerges as a property of cell interactions.
*PLoS One.*2014; 9: e104969 Yllanes, D. and Marchetti, M.C.. 2017. How many dissenters does it take to disorder a flock? arXiv, , 1701.05477v2, [cond-mat.stat-mech].

## Article Info

### Publication History

### Identification

### Copyright

### User License

Elsevier user license |## Permitted

### For non-commercial purposes:

- Read, print & download
- Text & data mine
- Translate the article

## Not Permitted

- Reuse portions or extracts from the article in other works
- Redistribute or republish the final article
- Sell or re-use for commercial purposes

Elsevier's open access license policy