Directional Collective Cell Migration Emerges as a Property of Cell Interactions

Collective cell migration is a fundamental process, occurring during embryogenesis and cancer metastasis. Neural crest cells exhibit such coordinated migration, where aberrant motion can lead to fatality or dysfunction of the embryo. Migration involves at least two complementary mechanisms: contact inhibition of locomotion (a repulsive interaction corresponding to a directional change of migration upon contact with a reciprocating cell), and co-attraction (a mutual chemoattraction mechanism). Here, we develop and employ a parameterized discrete element model of neural crest cells, to investigate how these mechanisms contribute to long-range directional migration during development. Motion is characterized using a coherence parameter and the time taken to reach, collectively, a target location. The simulated cell group is shown to switch from a diffusive to a persistent state as the response-rate to co-attraction is increased. Furthermore, the model predicts that when co-attraction is inhibited, neural crest cells can migrate into restrictive regions. Indeed, inhibition of co-attraction in vivo and in vitro leads to cell invasion into restrictive areas, confirming the prediction of the model. This suggests that the interplay between the complementary mechanisms may contribute to guidance of the neural crest. We conclude that directional migration is a system property and does not require action of external chemoattractants.


Introduction
The Neural Crest (NC) is a multi-potent cell population that arises at the dorsal midline during embryo development, migrates ventrally through the embryo and is guided by strict migratory pathways [1]. Collective cell migration is an important biological process that occurs during development [2], wound healing [3], cell renewal [4][5][6] and metastasis [7]. Recent efforts have identified the NC as a suitable model for collective cell migration [8,9] and for metastasis, as similarities between the NC and metastatic cancer cells have been observed [10,11]. The mechanisms that regulate collective cell migration are not fully understood, however data suggests cranial NC cell migration both in vivo and in vitro, is regulated by 1) contact inhibition of locomotion (CIL) [12,13], 2) chemotaxis towards a self secreted chemoattractant [14] and 3) a collection of external negative signalling molecules such as Eph/Ephrin and Robo/Slit, for a review see [8]. CIL was discovered by Abercrombie and Heaysman [15,16] and has been extensively studied in a range of experimental systems [17][18][19][20][21]. Mechanically, CIL can be described as a change in motion of individual cells due to contact, and occurs to differing extents in migratory cell types, such as fibroblasts [15][16][17], keratinocytes [22], Drosophila macrophages [23], NC [12,[18][19][20][21] and the PC-3 cancer cell line [16,[22][23][24]. This process has been characterized in Xenopus, chick and zebrafish NC and has been demonstrated as a key mechanism that confers cell polarity, by regulating the activity of small GTPases, and controlling directional migration of the whole NC population [12,25] Attraction between NC cells has been observed to take place concurrently with CIL [14]. The complement factor C3a and its receptor C3aR were found expressed in the migrating NC and previous work has demonstrated that C3a plays a direct role in collective migration, functioning as a homogenous NC secreted chemo-attractant. This chemokine acts to maintain a high cell density through homotypic attraction, a phenomenon called coattraction [14]. As CIL and co-attraction have been described as two microscopic processes with opposite effects on NC cells (repulsion and attraction, respectively), it is not evident how the relative contribution of these two contrary forces could affect directional migration of a group of cells. To better understand whether CIL and co-attraction could control directional migration of NC cells we developed a mathematical model of the process.
In previous studies, agent based models have addressed the transition from disordered to ordered motion in swarming insects [26] and the dynamics of wound healing assays [27]. These models assign rules to individual agents, from which, changes in local interactions lead to phase transitions, such as a parallel to a torus state [28]. Collective migration in real cells has been compared to flocking behavior modeled in the coordinated movement of animal groups and experiments have confirmed shared properties such as local correlation and responses to the local environment [29,30]. Further analysis in agent based methods has led to analytic approximations on the force required to maintain a particular state [31] and quantification of adaptable interactions to the local environment [32]. Force-based models have provided an alternative descriptive modeling approach allowing parameter prediction based on a macroscopic feature such as group alignment. In a study on keratinocytes, long-range order was shown to depend on repulsive and adhesive forces [33] and in a mechanical model, wave propagation has been described in epithelial monolayers [34].
Several models have been proposed [35][36][37] in which NC cells have been assumed to migrate in response to external chemoattractant gradients. Although there is evidence to support the presence of NC chemoattractants in vivo [38][39][40], it has been very well documented that NC cultured in vitro in the absence of any external chemoattractant exhibit directional collective migration [40,41]. The effect of random perturbations in collective migration has been analysed [35] and the stability of NC chains characterized [37]. In the study of Wynn et al., agents were simulated on a grid and parameter analysis was performed on an initial pattern to test chain persistence with leader and follower cells. Further investigation suggested that cell interactions with the ECM, directional bias and cell contact could play a mutual role in the promotion of chain migration [42]. A different study that compared theory and experiment used an off-lattice individual based model combined with a continuous model of vascular endothelial growth factor to predict behaviour of cranial NC migration in chick. The results of this study suggested that a combination of leading and trailing cells are required to ensure cohesive movement and collective response to external signals [36]. Alternative models have demonstrated network formation in the absence of external gradients, for example where simulated cells are cued by strains [43] and in the NC, where rules of movement include agent path reinforcement and repulsion or preference to follow existing axons [44].
In this work, a model of NC collective migration is presented. In a similar fashion to a study on fibroblast migration [45], microscopic parameters are estimated from biological data analysis and through simulation, macroscopic features of migration are predicted and compared with experimental data. In contrast to Vedel et al. [45], where the effect of local parameters on relative simulated cell dynamics was assessed in the form of an autocorrelation function, we focus on both correlations between velocities and collective properties of the group in the form of a time taken to reach, collectively, a target location, allowing us to assess long-range dynamics. When a group of NC cells is plated on fibronectin, they are able to migrate collectively, and with directionality, in the absence of any external signal. Additionally when the leading edge is removed, previously trailing cells continue migrating assuming a leading phenotype [12,14,40]. Previous models have investigated leading and trailing populations [46] but some have employed different mechanistic rules amongst simulated cells [37]. Models that include predefined differences between leading and trailing cells do not take into account the emergent heterogeneity that can arise as a consequence of the dynamical system. In addition, it has previously been shown that a combination of repulsive (CIL) and attractive (co-attraction) forces could generate directional migration [14]; however, this model did not consider the migration parameters of real cells. Although the model can generate an efficient migration when CIL and coattraction are combined, it does not reproduce the real behaviour of cell clusters when only CIL is present or the real behaviour of single cells. Hence, there is a need to construct a model that reflects the biological observations in the cranial NC that will be better suited for comparison with functional experiments.

The model
We present a microscopic model based on a periodic change in polarity, resulting in a change of direction, which we call rotational turning ( Figure 1a) and the processes CIL and co-attraction. Measured properties of these interactions (see Text S1, Figure S1), taken from Xenopus NC cells migrating in vivo and in vitro, are incorporated to the model, which follows the discrete element method [47].
We abstract NC cells to elastic spheres that we refer to as simulated cells. For a population of size N, each simulated cell is equipped with a natural radius r i~r V i [ 1, . . . ,N f gand a ray R CoA in the direction of polarity corresponding to the sensing range of the simulated cell. In addition, each simulated cell is assigned a mass m and intrinsic speed v j j. In the event that contact occurs between simulated cells, normal contact forces are modelled with Hertz contact theory (see Text S1). Data analysis of CIL in vitro confirms that the mechanism of contact inhibition is significantly different from the dynamics of an equal mass normal force rigid body collision. To account for this, the model is modified through the addition of a repolarisation force that acts in a randomly distributed direction at the free edge, see (Figure 1bc). This implementation is different to previous models of swarming that have assumed inelastic collisions [48] and is consistent with experimental data, as the generation of protrusions at the free edge has not only been observed in vitro in Xenopus but also in vivo in Zebrafish, see (Figure 1d-f). Single NC cells observed in vivo periodically change their direction of migration [12,49]. This change in direction of migration is dependent on the direction of their protrusions and can be observed by plotting individual cell tracks or recording cell persistence. To account for this behavior in our model, each simulated cell is assigned two internal clocks that periodically switch on a force due to coattraction and an impulsive force due to rotational turning. Currently these rates are unidentified experimentally. In the event that a simulated cell responds to co-attraction, the simulated cell is subjected to a force proportional to the gradient of the coattraction profile, as the steepness of external gradients have been shown previously to affect cell motility in eukaryotic cells [50] ( Figure 1g-i). Simulations were performed in a 2D continuous geometry, to represent the permissive extra cellular matrix, with a rigid wall at the dorsal border and a repulsive cue at the lateral borders to represent negative signals that are known to be present in the embryo at the border of each NC stream. It is known that some of these molecules are secreted, like semaphorin [8], which would generate a gradient consistent with the model. In the event that a simulated cell responds to the lateral repulsive gradient, it is subjected to a force proportional to the gradient, which is localised at the border. The domain is equipped with a target at the opposite end of the domain from the initial location of the simulated cells see (Figure 1j). When a simulated cell reaches the target, it remains stationary for the remainder of the simulation, which represents the real cells ceasing migrating once they reach the target tissue in the branchial arches. The extraction of the simulated cells facilitates the analysis of efficiency in directional migration by quantifying the number of cells that reach the target (see Text S1).

Time integration
The dynamics of cellular motion are driven by the sum of the applied forces. Simulated cells maintain a 'preferred' self-propelled velocity v v i (t) in the absence of the forces tested. This is an abstraction of the real biological scenario, where the velocity of migration is generated mainly by actin polymerization at the cell front [51] and single cells exhibit differences in their velocity over time. However this simplification allows us to more clearly explore how interaction forces influence the group level dynamics. A simulated cell always moves in the direction of its polarity. The force that governs the migration of a simulated cell is presented below and simulations were performed with the iterative central difference model [47].
Here, m is the mass, € u u u u i is the acceleration of the simulated cell with position vector u u i .
F F T i is the total interaction traction force that will influence a change in the velocity.
are the co-attraction, rotational turning, self-propulsion, contact, contact damping and contact repolarisation forces (see Text S1). C i is the set of indices of simulated cells in contact with simulated cell i (fk : d ik j j: u u i { u u k j j{r i {r k v0g, see Figure 1k). The coefficient Q sets the co-attraction, self-propelled and rotational turning forces to zero when a simulated cell is in contact and the coefficients a i (t),v i (t) and m i (t) are functions of the internal clocks for co-attraction I iCoA (t) and rotational turning I iRT (t) (see Text S1).

Model Calibration
Where possible we have attempted to match model parameters to the control real cell biological data. Following Wynn et al. and Carmona-Fontaine et al. [14,37] baseline parameters were chosen that correspond to physiological conditions and are presented in Table S1. The computational domain was defined with a height of 850mm and width of 217mm. The simulated cell diameter was uniformly defined as 40mm to approximate the cell width observed in vitro and in vivo, with the simulated cell speed estimated from biological data as 3mm every minute (5.0e-8m=s) [12].
To construct a model that represents the microscopic interaction of the real cells during contact inhibition, we analysed three quantitative values to parameterise force equations based on the theory of contact mechanics. These values were the angle before and after CIL, the contact time and the acceleration after contact. During contact and CIL, cellular motion is modelled as a function of the normal contact force and a repolarisation force. These two forces represent the material properties of a cell and the fact that CIL activates a molecular signalling pathway, which affects molecular activity at the free edge, promoting protrusion formation. It is known that protrusions are inhibited at the site of contact, via a mechanism involving cadherins and Rho-GTPases [52][53][54] (for a review see [55]). In addition, the PCP pathway regulates repolarisation at the free edge [12]. The repolarisation force is not present in standard discrete element models. The force acts in the direction of the unit vector connecting the two cells centre's of mass plus a random angle sampled from a uniform distribution with the range {p=2,p=2 ½ (see Text S1). The exact distribution does not have a significant effect on the contact model, see Figure S2a-b. To test the influence of the normal contact force on the collective dynamics of the cells and to understand whether the relative velocity of two cells during contact was a significant factor in the model, the angles between the paired velocities of two biological cells prior to and following contact were analysed by assessing if they were correlated. Where possible microscopic parameters were approximated with real cell data from the literature [12,14,56]. The simplest form of contact is that of a normal force rigid body elastic collision, which we expected would give rise to highly correlated pre-and post-contact angles. First we tested whether these angles were independent in the experimental data. By assuming that the pre-and post-contact angles have a bivariate normal distribution, testing for independence becomes equivalent to testing whether the correlation coefficient r is zero. The hypothesis H 0 that the pre-and post-contact angles exhibited a correlation coefficient r~0 was tested for the whole data set using a two tailed t test with where r is the sample correlation coefficient and n is the sample size. The biological data yielded a sample correlation coefficient of r~0:036 and a statistic t~0:4 for n~128, which suggests that the pre-and post-collision angles are not correlated (Figure 2c). This test was repeated for cells that remained together for at least the mean time yielding r~0:252 and t~0:97, which again supports the hypothesis that pre-collision angles are not correlated with post collision angles. A low correlation indicates that the predominant forces involved in contact inhibition are not due to normal contact forces but some other mechanism, which we define loosely as repolarisation. To compare these results to the model, a normal force elastic body collision scenario was tested. This model neglects the terms F F Cd ik and F F RCIL ik , which constitute the energy dissipation and repolarisation forces. The normal force model alone in the absence of these terms exhibited an r value of 0:97 and a test statistic t~22:62, indicating that pre-collision and post-collision angles are correlated (Figure 2a). When repolarisation and contact damping were included, the model exhibited an r value of 0:25 and t~0:98, indicating that the pre-and post-collision angles are not correlated as in the case of the biological data, see (Figure 2b). We then looked at the time that cells remain in contact. Experimental movies had a frame rate of 5 minutes so that the minimum contact time recorded would be less than 10 minutes. The experimental distribution in time was compared to both the normal force elastic body collisions and the repolarisation model, (Figure 2d-f). During CIL and contact separation, it takes some time for the real cells to regain their default migratory speed. There is variation in this acceleration. To obtain quantitative data on this process, the speed of a real cell upon contact separation was recorded over time. The results confirm that to regain the default speed, cells must accelerate after contact, which suggests that CIL cannot be fully described by a normal force rigid body elastic collision. Speed after contact was recorded for the repolarisation model. In contrast to the experimental data, the simulated cell was unable to accelerate to 7|10 {8 m=s, however this is due to the default speed being set to 5|10 {8 m=s. When the default speed was increased to 7|10 {8 m=s, the simulated cell's speed increased to a value greater than 6|10 {8 m=s, see Figure S2c. Compared with the normal force elastic model, the repolarisation model can better explain the change in speed after contact, (Figure 2g-i). Together these results suggest that rigid body collisions are not sufficient to model contacts between cells, but that our model, which incorporates a novel repolarization force can better do so.
The frequency of the clocks cannot simply be measured. Therefore conclusions on the effect of these parameters were drawn after parameter analysis. We take a baseline rate for chemotactic response to be one response every two seconds and for reorientation, one random change in direction every five minutes.
We have shown that C3a (the co-attractant) forms a stable gradient by binding to fibronectin [14]; this gradient was measured and a 2D mathematical radial diffusion model was assumed. To model the chemoattractant, we assume a steady state distribution at every iteration, as the timescale for diffusion is smaller than the time it takes for a cell to move a significant distance. This steady state distribution can be described as a Bessel function [57], which for simplicity we approximated with a decaying exponential with a half maximum length l ð Þ of 110mm.

Results
Co-attraction and CIL are sufficient and necessary (in silico) for directional collective migration Migration is found to be a qualitative fit to the behaviour observed for real cells [14], for example, in the absence of an external bias the cells migrate in a coordinated fashion leading to the displacement of the group as a whole. This suggests that the model can reproduce directional migration with the functional processes CIL, co-attraction and rotational turning. In the presence of CIL and co-attraction, directional migration occurred as a travelling wave of density, which reproduces the directional migration observed in real cells [12,14]. To test the relationship between directional collective migration, co-attraction and CIL, four cases were considered: (1) 2CIL,2CoA corresponding to an elimination of all processes except rotational turning, (2) +CIL,2CoA representing a complete knockdown of co-attraction, (3) 2CIL,+CoA which tests the model under the assumption that CIL is inhibited and (4) +CIL,+CoA corresponding to the baseline case (Figure 3a-d, Video S1). Out of all four cases +CIL,+CoA produced the most efficient migration through the domain, in which the centre of mass of the group was the most distal at a simulation time of approximately 2 hours. To quantify this efficiency, we define directional migration as the combination of a high coherence and low target time (Text S1), where high coherence corresponds to a value greater than 0.5. Case (4) was unique in displaying these properties, suggesting that both CIL and co-attraction are necessary for directional migration in the model (Figure 3e-h).  In cases where co-attraction is absent, the macroscopic behaviour is similar to diffusion in a bounded domain (Figure 3a-b). To assess if directional migration requires CIL, the model was simulated in the absence of F F RCIL ik and with continual F F p i during contact. In contrast to the dynamics of a co-attraction knockdown, elimination of CIL resulted in minimal displacement of the bulk population (Figure 3c). These results were upheld under further analysis of velocity, where coherence was high under baseline conditions and low when either CIL or co-attraction was impaired (Figure 3e-h). To compare the model predictions with biological data, automated tracking software was used on both the model and experimental data at the same frame rate frequency [58]. Individual simulated cell and cell velocity was tracked and the coherence computed for a control experiment, the baseline parameters and a knockdown of co-attraction in the model. Under baseline conditions, the software calculated a coherence of 0.6252. To analyse the predictive quality of the model, control experimental data was processed and exhibited a coherence of 0.5568. To test the goodness of fit of these values, a simulation lacking coattraction was processed. For this case, the coherence was 0.008, suggesting that baseline parameters are a better fit to control migration and they reproduce the cell behaviour observed in vivo and in vitro (see Figure 3h, Figure S3a-c).

Migration is diffusive if co-attraction is not sufficiently strong
To test the robustness of directional migration sensitivity analysis was performed by considering the effect of one parameter at a time on directional migration. This was implemented for five physiological parameters, consisting of the C3a diffusion length l, the angle by which the simulated cells can deviate during rotational turning, the rates of the internal clocks I iCoA (t), I iRT (t) and the domain length ( Figure 4). Collective migration occurred in baseline simulations and was maintained under small parameter variation. There was variation in the collective target time for baseline parameters between independent simulations. This variation is negligible when compared to the collective target time for the diffusive state and we refer to the time in which the group remains travelling in one direction as the collective flight time of the group (see Figure S3d-f, Video S2a). We performed a Mann-Whitney U test on the coherence between consecutive parameter values presented in Figure 4a-e to test for a difference in medians between consecutive parameter values. In agreement with previous studies, frequent reorientation resulted in a low coherence (1/(RT rate)~3s, Figure 4a,f) and there is evidence to suggest that there is a difference in the median coherence between parameter values 1/(RT rate)~3s and 1/(RT rate)~450s, (p, 0.0005, n = 10). For baseline parameters, different angles by which the simulated cells can deviate did not disrupt collective migration (Figure 4b,g) and there was no significant difference in the median coherence between consecutive parameter values, (p.0.01, n = 10). The effect of co-attraction on group level dynamics was tested by variation of 1/(CoA rate) and the diffusion length l (Figure 4c,d,h,i). The results show that simulated cells fail to directionally migrate when there is an infrequent response to coattraction, or if the co-attraction gradient is too short-range. For example, there is a significant difference in the median coherence between a response parameter of 1/(CoA rate)~0:1s and 1/(CoA rate)~20s, (p,0.0005, n = 10). Similarly for the gradient, there is a significant difference in the median coherence between l~11mm and l~55mm (p,0.005, n = 10). This co-attraction dependent transition between directional migration and dispersion occurred at a spatial occupancy of A 0 A~0 :34, where A 0 is the area occupied by the cells and A is the total area and this density occupancy was held constant across all simulations. In previous studies of epithelial cell populations, group coherence is exhibited with densities greater than 0.2 [33]. This result suggests that mesenchymal cell populations such as the NC, may naturally disperse in the absence of a co-attractant but the response to coattraction regulates this behaviour and allows cells to acquire motion similar to those of epithelial cell types. To characterise the transition from diffusive to directional collective migration, coherence and target times were recorded under variation of the box height H, for weak and strong co-attraction (Figure 4e,j). The coherence was high and the collective target time increased linearly with H for strong co-attraction. In contrast the coherence was low and the target time increased super linearly for weak coattraction (Figure 4e,j). To obtain an upper bound on the rate of co-attraction, we considered small values of 1/(CoA rate). The coherence was recorded for values between 1/(CoA rate)~0:003s and 1/(CoA rate)~2s, (Figure 5a). Coherence was maintained within this range. In contrast to this, the speed of the simulated cells changed, such that the average speed of a simulated cell during 1/(CoA rate)~0:003s is 1:6|10 {8 m=s. This suggests that although coherence is maintained at high rates of co-attraction, efficiency of bulk displacement is reduced. Furthermore, a reduction in speed coincides with a longer collective target time, suggesting that there is an optimal response rate to co-attraction. The collective coherence, speed and collective target time were recorded for five different rates of I iCoA (t) (Figure 5a-c). From this data, we suggest that although there exists a range of I iCoA (t) rates within which directional migration can emerge (see error boundaries, Table S1), the optimal value for the parameters tested coincides with 1/(CoA rate)~0:1s.

Model predictions
Under baseline conditions, simulated cells remained positioned within the permissive regions during migration to the target (Figure 6a). In contrast, when co-attraction was inhibited, we noticed that single simulated cells appear to cross into the lateral restricted regions (Figure 6b). To quantitatively validate this behaviour, we recorded the average number of simulated cells that reside in the restricted region for a response to the boundary signal every 9 s. This number was recorded as a percentage of the population for cases (2) and (4) in the model. On average, the percentage of cells that crossed the lateral border was close to 0% in the baseline condition, whereas it was close to 14% when coattraction was inhibited (blue in Figure 6e). To test whether this unexpected prediction of the model was also found in real cells, an experiment to reduce co-attraction in vivo was performed. An antisense morpholino was used against C3aR to inhibit coattraction. In this experiment, control cells remain positioned within their migratory streams, (Figure 6c) however, invasion of cells into non-permissive area was observed for cells depleted of C3aR (Figure 6d). To quantitatively compare simulated and real cells, an in vitro experiment was performed. NC explants were cultured on corridors of fibronectin flanked by fibronectin-free regions, and time lapse analysis was performed, (see Figure S4). NC cells need fibronectin for their migration as they inefficiently attach to a fibronectin-free substrate. While control cells rarely invaded the fibronectin-free region ( Figure S4a) an important proportion of the C3aR depleted cells moved into that region ( Figure S4b). The percentage of cells invading the prohibited region was similar between the model and the real cells (Figure 6e). The average percentage of cells to cross into the restricted region throughout the simulation was recorded for different boundary clocks I iB (t) (see Text S1). For values of 1/ (I iB (t) rate)~1s, 1/(I iB (t) rate)~10s and 1/(I iB (t) rate)~20s, the average percentage of cells to cross into the restricted region for control simulated cells remained within 5%. For 1/(I iB (t) rate)~100s, the average percentage of cells to cross into the restricted region for control simulated cells was greater than 50%. The model predicts that co-attraction facilitates guidance of the NC by counteracting the migratory force that is sufficient to overcome negative signals. By assigning a force term to the processes of co-attraction in the model, we were able to compare model predictions with functional experiments. For this behaviour, the results of the model and experiment are in agreement.
Emergent behaviour was not limited to stream guidance. Previously, NC explant confrontation was performed to directly test co-attraction. In this experiment, explants are cultured within a distance that is great enough to ensure that no initial contact occurs between the groups. It is known that the groups consistently move toward each other, however the number of cells in each group determines the distance over which co-attraction can act. This property was consistent with the model, where groups of NC cells respond to co-attraction at greater distances than single cells (see Video S3, S4, Figure S5). Differences in velocity were observed between leading and trailing cells when co-attraction occurred at every 0.1 s or less (see Figure S6, Video S2b). Leading and trailing behaviour has previously been shown to occur in the chick and Xenopus NC [12,37,40]. We suggest that this could emerge in a population of identical cells without requiring differentiation of microscopic parameters, as it has been shown for NC migrating in vitro [12].

Discussion
In this study, two processes (CIL and co-attraction) that occur in migratory cells are analysed with an application of the discrete element method. This method is commonly applied in the field of mechanical engineering and here we apply it to cell migration.   Analysis of experimental data on cell collisions shows that normal contact forces alone cannot well describe these collisions. The inclusion of a repolarisation force, generated by the retraction of protrusions at the contact site and the formation of new protrusions at the free edge, allows a much better fit to data.
Coherence and efficiency of bulk displacement of simulated cells indicate that CIL and co-attraction are both sufficient, and necessary, for effective directional migration. Upon inhibition of either process, collective migration is disrupted and the dynamics consist of high cell density meandering clusters in the absence of CIL or low cell density random movement in the absence of coattraction. Qualitative and quantitative measurements of the global dynamics of the simulation are compared with experimental data and application of automated tracking software [58]. We identify a range of possible values for the co-attraction rate. By designing a model parameterised with experimental data at the microscopic scale, we demonstrate that directional migration is robust to small changes in the processes of CIL and co-attraction, however both infrequent and continual co-attraction can disrupt timely directional migration. This model allows us to explore the effect of variation in the microscopic parameters on collective behaviour, to support existing experiments or to make predictions when real experimental values are unknown.
In contrast to previous studies, our model predicts that coattraction contributes to the guidance of the NC by promoting directional migration and inhibiting single cells from migrating into restricted regions. Here, we provide experimental evidence that confirms the model prediction in vivo and in vitro.
Where feasible, processes are parameterised using experimental data. This uncovers a timescale suitable for modelling contact, a dynamical process with resolution on the millisecond scale. To assess and quantify long-range dynamics, simulations are on the scale of hours. To facilitate these intensive simulations, the code is implemented in CUDA (Compute Unified Device Architecture). This allows us to exploit the highly parallel nature of graphical processing units (GPUs) for quantitative analysis.
Previous models, have shed light on the emergent dynamics arising from individual interactions inspired by biological data, however several have depended on the presence of an external chemo-attractant [35][36][37]. By contrast, this study does not invoke an external chemoattractant but considers instead chemotaxis towards a self-secreted chemoattractant C3a and tests the role of co-attraction in collective migration and stream guidance. Biological evidence suggests that individual interactions between cells work together to allow self-organization in migrating clusters and collective migration [12,14]. This study confirms this from a mechanical perspective and suggests that CIL and co-attraction promote migration similar to epithelial directional migration in cell populations that are mesenchymal with low cell-cell adhesion. In addition to CIL and co-attraction, external signalling, such as chemoattractants and chemorepellents have been shown to play a role in NC migration [1]. The model presented here suggests that long-range directional migration is acquired through local tissue specific interactions and permissive cues that have the greatest effect on migrating collectives. This could potentially promote a flexible system that is ready to adapt to external and internal perturbations. . Rapid response with 1/(CoA rate) = 0.008. Angle made with the vertical axes by leading, centre leading, centre, and centre trailing and trailing. As the data did not appear normally distributed a Wilcoxon signed rank test with continuity correction was applied to the leading and trailing data with a p-value,0.001. (d). Rapid response with 1/(CoA rate) = 0.008. Speed of cells partitioned by leading, centre leading, centre, centre trailing and trailing. As the data did not appear normally distributed a Wilcoxon signed rank test with continuity correction was applied to the leading and trailing data with p-value,0.001. (e). Baseline parameters. Angle made with the vertical axes by leading, centre leading, centre, and centre trailing and trailing. The same statistical test used in (c) and (d) indicated no significant difference between leading and trailing populations at baseline. (f). Baseline parameters. Speed of cells partitioned by leading, centre leading, centre, centre trailing and trailing. Wilcoxon signed rank test showed no significant difference between leading and trailing populations at baseline. (TIFF)

Supporting Information
Table S1 List of parameters used in the discrete element model and their values. Values were approximated from either experimental data, or through comparison of emergent behaviour between model and experiment. Where parameters have been chosen from sensitivity analysis, their error bounds are shown. ND represents a scalar parameter. (DOCX) Text S1 Further details of the equations of motion, the analysis and supplementary results. Video S3 Movie showing that a single cell cannot respond to co-attraction at a distance of 517mm, when placed in proximity to a group of 25 simulated cells.

(AVI)
Video S4 Movie showing that co-attraction acts a distance of 517mm,, when two groups contain 25 simulated cells. (AVI)