Modelling chase-and-run migration in heterogeneous populations

Cell migration is crucial for many physiological and pathological processes. During embryogenesis, neural crest cells undergo coordinated epithelial to mesenchymal transformations and migrate towards various forming organs. Here we develop a computational model to understand how mutual interactions between migrating neural crest cells (NCs) and the surrounding population of placode cells (PCs) generate coordinated migration. According to experimental findings, we implement a minimal set of hypotheses, based on a coupling between chemotactic movement of NCs in response to a placode-secreted chemoattractant (Sdf1) and repulsion induced from contact inhibition of locomotion (CIL), triggered by heterotypic NC–PC contacts. This basic set of assumptions is able to semi-quantitatively recapitulate experimental observations of the characteristic multispecies phenomenon of “chase-and-run”, where the colony of NCs chases an evasive PC aggregate. The model further reproduces a number of in vitro manipulations, including full or partial disruption of NC chemotactic migration and selected mechanisms coordinating the CIL phenomenon. Finally, we provide various predictions based on altering other key components of the model mechanisms. Electronic supplementary material The online version of this article (10.1007/s00285-019-01421-9) contains supplementary material, which is available to authorized users.


Introduction
The neural crest is a transient structure, unique in embryonic vertebrates, that provides various neural derivatives which populate distinct tissues. It forms as a ribbon of ectoderm-derived cells and extends head to tail along the dorsal neural tube (Le Douarin and Kalcheim 1999;Theveneau and Mayor 2012). Notably, a controlled spatiotemporal epithelial-mesenchymal transition (EMT) allows detachment of its component cells, followed by their further migration and dispersion throughout the embryo as individuals, chains, clusters or sheets. Failure of neural crest cell (NC) dispersion/proliferation/differentiation can lead to various embryonic pathologies, collectively termed neurocristopathies (Watt and Trainor 2014). Beyond the critical function in early vertebrate life, neural crest dispersal offers a model system for studying the mechanisms underlying other cell invasion processes, such as those occurring during fibrosis or cancer development.
Various factors are implicated in NC migration and guidance, including homotypic and heterotypic cell-cell interactions, a combination of "go" and "no-go" factors in the extracellular matrix (ECM) and surrounding tissues that act to restrict NC migration along specific pathways, and long-range diffusible factors that provide chemotactic guidance (e.g., McLennan et al. 2010;Shellard and Mayor 2016;Theveneau et al. 2010;Mayor and Theveneau 2013;Tosney 2004;Twitty 1949). In particular, the study by Theveneau et al. (2010) has shown that Xenopus laevis NCs exhibit (positive) chemotaxis in the presence of gradients of the extracellular ligand Sdf1. Specifically, Sdf1 binds to the cell membrane receptor Cxcr4 and promotes intracellular Rac1, a key player in the activation and stabilisation of the cell motility structures (e.g., filopodia, pseudopodia) that lead to cell movement.
As remarked above, NC migration is also regulated by cell-to-cell contact interactions, which can be attractive, as in adhesion, or repulsive, as in contact inhibition of locomotion (CIL). CIL was first identified in vitro more than half a century ago (Abercrombie and Heaysman 1953), when the contact between two migrating fibroblasts was shown to lead to a transient arrest in their motion, a repolarisation and a subsequent reversal of migration heading. CIL therefore acts to promote cell repulsion and, intuitively, it could enhance dispersal. Current interest in CIL has been sparked by demonstrations that it also occurs in vivo, during dispersal of NCs in Xenopus laevis and zebrafish (Carmona-Fontaine et al. 2008;Theveneau et al. 2013). Further discoveries of its operation in cancer cell populations (Astin et al. 2010), developmental macrophages (Stramer et al. 2010) and neural cells (Villar-Cervino et al. 2013) have reinforced its relevance for migration and invasion processes. In vivo, neural crest migration takes place in a highly heterogeneous environment where interactions with surrounding tissues or populations are inevitable. Recent studies of NC migration in Xenopus laevis indicate an active interplay between NCs and the epithelial-type "placode cells" (PCs) that initially lie adjacent to the neural crest . Such heterotypic interplay involves both long-range and contact-mediated interactions: -PCs secrete the diffusible ligand Sdf1, which (as described above) acts as a chemoattractant for NCs and draws them towards PCs -the chase phase of the collective movement of the NC-PC system; -direct contact between NCs and PCs then initiates a CIL response, invoking their movement away from each other -the run phase of the collective movement of the NC-PC system.
Cell-cell contacts are mediated through various signalling pathways, typically triggered by linkage of membrane-bound receptors on adjacent surfaces. In the case of NCs and PCs, cadherin family members (classically associated with adhesion) have been shown to play a significant role in their mutual dynamics. The initially attracting (adhesive) interactions that arise through N-cadherin-N-cadherin binding can subsequently give way to a repelling CIL response, mediated via a downstream signalling process. N-cadherin binding leads in fact to Rac1 downregulation, which in turn suppresses local cell membrane protrusions . Thus, protrusions become biased to the opposite end of the cell membrane from where the contact occurred and the individual cell is repolarised accordingly. Overall, N-cadherins therefore appear to generate both attracting (adhesion-type) and repelling (CIL-type) dynamics. While NCs solely express N-cadherins, placode cells also express E-cadherins, which generate stable homotypic E-cadherin bonds therefore promoting stable PC clustering.
In the case of an in vitro aggregate of NCs juxtaposed against a similar aggregate of PCs, this "chase-and-run" process generates a net movement of the overall system, in which the NC cluster continuously chases PCs and is both repelled by and repels the PC population, see Fig. 1. In vivo, it has the potential to draw the NC population away from the neural crest and into the surrounding tissue, clearing a pathway through the PCs on route.
A number of theoretical works have addressed aspects of neural crest migration/collective migration: we refer to Schumacher et al. (2016) and Szabo and Mayor Fig. 1 Schematic of "chase-and-run" dynamics. a NCs follow the chemoattractant gradient formed from Sdf1 producing PCs . PCs are tightly attached through E-cadherin binding , while NCs are more loosely connected via a combination of a C3a/C3aR co-attraction mechanism and NC-NC CIL responses (Carmona-Fontaine et al. 2008. b N-cadherin cell-cell interactions occur as the NC and PC populations collide, resulting in CIL in these two cell types . c The CIL response results in mutual repulsion. Overall, net movement of the two populations results (2016) for specific reviews on the modelling of neural crest dynamics, and to Camley and Rappel (2017) for a review of collective cell migration models in general. With regard to CIL phenomena, one-to-one interactions between migrating and colliding NCs within quasi-one-dimensional geometries have been modelled in Kulawiak et al. (2016) and Merchant et al. (2018). The former developed a computational phase field model, explicitly accounting for cell shape and intracellular biochemistry. Consistent with experimental observations by Scarpa et al. (2013), different collision outcomes were observed (reversals, sticking and walk-past) according to the balance of different factors, such as adhesion. The model by Merchant et al. (2018) also investigated the outcome of collisions along 1D lines, in a model where each cell's membrane is represented by a closed chain of elastic edges. Rac1/RhoA biochemistry was included via kinetic equations, with CIL and a NC-NC co-attraction process (mediated via the complement factor C3a and its receptor C3aR, see Carmona-Fontaine et al. 2011) incorporated through their action on kinetic terms. Sticking/reversal dynamics at a one to one level have also been investigated in Mayett et al. (2017) in an interacting elastic dimer model, where phase boundaries for the transition from clumping to chase-and-run were described for a NC colliding with a PC.
Extending to multicellular clusters, agent-based models have been analysed in Carmona-Fontaine et al. (2011) and Woods et al. (2014) to study the complementary roles of co-attraction and CIL in a population of homogeneous NC cells. Including both mechanisms allowed aggregates to migrate in an efficient, cohesive and directional manner. Along similar lines, a minimalistic stochastic particle model featuring a chemoattractant-regulated CIL process was shown to be sufficient to explain efficient cluster migration (Camley et al. 2016a), while extensions generated more sophisticated dynamics such as particle rotations about a core (Camley et al. 2016b). In , the ECM molecule versican was shown to display a complex spatio-temporal pattern of tissue expression which, when combined with its capacity to restrict NC migration, was supposed to channel NC invasion along specific routes. A Cellular Potts Model approach that incorporated CIL and co-attraction between NCs alongside versican-mediated inhibition of migration helped verify this hypothesis. The above described model by Merchant et al. (2018) has also been applied to explore NC cluster migration, demonstrating how these two mechanisms generate spontaneous and efficient collective migration.
Beyond these studies targeted at CIL/chase-and-run dynamics, further theoretical studies have addressed NC migration in other systems. Cranial neural crest cell migration in chicken embryos has been explored by colleagues (McLennan et al. 2012, 2015a, b), where agent-based models were used to understand the role of "trailblazer" cells in forging a path through the tissue. Specifically, a trailblazer population undergoes chemotaxis in response to self-created gradients of an externally-produced attractant, while follower cells simply chase the leaders. Other neural crest focussed models have been developed to describe invasion of mouse enteric neural crest cells (Cheeseman et al. 2014;Landman et al. 2007Landman et al. , 2011Simpson et al. 2007) and mouse melanocyte cells (Mort et al. 2016), although in those systems cell migration is augmented with significant proliferation that helps drive the dispersal process.
The principal aim of our work is to model "chase-and-run" collective behaviour observed in vitro in multicellular NC-PC systems. Within this scenario, any cellular growth, birth or death processes appear to be minimal, allowing us to solely focus on the interactions that drive coordinated movement. We propose a hybrid multiscale approach, in which cells are individually described as microscopic/discrete interacting particles and PC-produced Sdf1 is represented by a continuous concentration distribution. Moving beyond the studies described above, we specifically consider chase and run within multicellular and heterogeneous clusters, composed from both NCs and PCs. Further, we model the dynamics of the extracellular chemical substance (Sdf1) via an explicit evolution equation that describes its spatiotemporal dynamics. For computationally manageability and limiting the dimensionality of the parameter set, we formulate a minimalistic set of interactions in order to understand the basic requirements necessary for "chase-and-run" dynamics. In silico experiments reveal the model's capacity to replicate features observed in vitro, both in the case of one cell-to-one cell and cluster-to-cluster interactions. Further, our model reproduces a number of experimental perturbations that target key mechanisms involved during "chase-and-run" dynamics and can be used to make a number of testable predictions.

The model
The overall theoretical framework consists of a proper set of ordinary differential equations (ODEs) for the cell dynamics and a reaction-diffusion (RD) law for Sdf1 kinetics. Our general system involves n P placode cells and n N neural crest cells, restricted so that they migrate across a planar domain Ω ⊂ R 2 (e.g., an experimental Petri dish as used for in vitro assays). Regardless of phenotype, each cell is represented by a discrete pointwise agent and can be identified by its position in space, say x α i (t) ∈ Ω with α ∈ {N , P}, i = 1, . . . , n α and t ∈ R + . However, physical cell sizes are taken into account via defining the evolution equations for cell distribution and movement. Specific scenarios, such as heterotypic interactions between single cells or homotypic aggregates, follow through tuning the sizes of n N and n P . In the rest of this section we motivate and define the underlying equations upon which our model is formulated, introducing a set of parameters in the process. For reference, these parameters are subsequently summarised in Table 1 while Sect. 3.1 lays out our parameter estimation. Specifically, we highlight those estimated from source values and those obtained via a data fitting process.

Cell dynamics
Cell dynamics are described by a first-order ODE system, which is derived from a general second-order particle model following a standard set of simplifying biological considerations. First, cells move in environments characterized by very low Reynolds numbers (Van Liedekerke et al. 2015;Odell et al. 1981), where ballistic locomotion is only briefly maintained and straight-line displacements are shorter than the typical cell dimension. Consequently, we can neglect inertial effects and adopt a first-order model in which cell velocity, rather than acceleration, is proportional to the acting forces. This relation, termed as overdamped force-velocity response, lies at the heart of numerous discrete/ individual-based-model (IBM) approaches (see Drasdo 2003;   Scianna and Preziosi 2012 and references therein for details) and demands that cell dynamics are driven by cell-cell friction and cell-substrate friction. For the planar domain considered here, cell-substrate friction is particularly relevant since the cell will have a larger contact area with the substrate compared to neighbouring cells.
We therefore also neglect cell-cell friction terms so that, as a first approximation, cell dynamics can be described by directly postulating cell velocity contributions and including cell-substrate friction coefficients in their characteristic parameters (see Carrillo et al. 2018 for further details).
The velocity of each cell is then determined by a sum of contributions that stem from its long and short range interactions with other agents, which may be directed or mediated by chemical signalling. Motivated by the experimental literature (Steventon et al. 2014;Szabo and Mayor 2015;Theveneau et al. 2013), the hypotheses at the heart of the model are as follows.
hp 1 -Each cell has an intrinsic resistance to compression, mainly due to the stiffness of its nuclear and perinuclear region. This force acts to repel two cells  Table 1 for values of d r , d c , d a ). NCs and PCs do not reveal significant differences in morphology and dimensions and hence, regardless of cell lineage: d r denotes the mean perinuclear region diameter; d c is the mean cell body diameter; d a describes the maximal extension of the motility structures that protrude from the cell surface (e.g., filopodia, lamellipodia, pseudopodia). Right panel: conceptual representation of velocity contributions, given in Eq.
(1), in the simplified representative case of a heterotypic two-cell system, i.e., composed of the generic j-th PC and the generic i-th NC cell. Each bar indicates when a given velocity component is activated according to d N P i j . Specifically, solid bars refer to the dynamics of the i-th NC, while dashed bars represent velocity components affecting the j-th PC. In both instances, bar thickness denotes how their magnitude varies with cell-cell distance d N P i j separated by distances smaller than d r (the mean perinuclear region diameter), see Fig. 2. hp 2 -Two cells can form adhesive bonds through the activity of transmembrane proteins. In particular, N-cadherins regulate both heterotypic (NC-PC) contacts and homotypic (NC-NC or PC-PC) adhesive interactions, while homotypic PC-PC bonds are further augmented by E-cadherins. The force generated through adhesive bonds acts to attract two cells up to a distance d a , which denotes the maximal extension of protruding motility structures (lamellipodia, filopodia etc.) from the cell's centre, see Fig. 2. hp 3 -Heterotypic N-cadherin bonds trigger a Rac1-dependent contact inhibition of locomotion (CIL) response in pairs of interacting cells, whereby adhesive complexes collapse, cytoskeletons repolarise and individuals move away from each other (the run phase of the process). Specifically, we assume the CIL response begins when the distance between two interacting cells falls below a critical value, d c , which is a measure of the cell diameter, see Fig. 2. The run phase duration is characterised by timescales τ αβ for αβ ∈ {N P, P N }. hp 4 -NCs chemotactically migrate in response to PC-produced Sdf1. This directional locomotion (the chase phase of the process) is subsequently downregulated once NC-PC adhesive complexes form, as well as once the NC is compressed by surrounding cells. Due to the extracellular and diffusible nature of Sdf1, chemotaxis is assumed to operate over arbitrary distances (although with a response magnitude varying with the detected attractant gradient). hp 5 -Each cell, regardless of type, shows innate random wandering. This Brownian crawling is taken to be more pronounced for NCs over PCs, as observed by comparing the dynamics of homotypic aggregates (see in Theveneau et al. 2013, Fig. 2(a)-(b)).
The overall system of ODEs regulating cell behaviour can therefore be written as follows: Let us now describe in more detail the cell velocity contributions introduced above.
Cell resistance to compression The terms v N i, res and v P j, res reproduce short-range intercellular repulsion, activated when the body of a cell is compressed by other individuals. In particular, it is natural to assume that for each individual this velocity contribution: (i) results from the superposition of its pairwise interactions; (ii) is strictly dependent on the relative distance between the cells involved (i.e., metric interactions); and (iii) is directed along the unit vector connecting their positions. Hence, we can write where f αβ res : R + −→ R − , with α, β ∈ {N , P}, are interaction kernels determining the intensity of cell repulsive behaviour. Further, we define for i = 1, . . . , n α , j = 1, . . . , n β , and α, β ∈ {N , P}, where |·| denotes the Euclidean norm. Trivially, d αβ i j is invariant under changes to the order of superscripts and subscripts, whiler αβ i j undergoes sign change with the order of subscripts. In principle, there is a wide range of choice for the form of functions f αβ res , however, certain biological observations aid our selection. First, any contact-dependent interaction vanishes when the individuals involved are too distant. Moreover, cell resistance to compression typically increases as the distance between the interacting cells decreases. Hence, taking advantage of existing literature (see in particular Colombi et al. 2015Colombi et al. , 2017Scianna and Colombi 2017), a plausible option is  4), affect the dynamics of a cell of type α when its distance from a cell of type β is lower than its nuclear dimension, denoted by d r . As defined in Eq. (6), adhesive interactions arise when the relative distance between the two cells is lower than the maximal extension of plasma membrane motility structures, measured by d a . Right panel: representation of the function used to describe NC chemotactic migration ( f N P dist ). As defined in Eq. (10), any chemotactic response of the generic NC starts to be downregulated once it undergoes adhesive interactions with at least one PC (when z ≤ d a ), becoming negligible when the two cells are tightly packed (z ≤ d c ) (color figure online) with, as previously stated, d r being the mean diameter of the cell perinuclear region (see Fig. 2) and taken to be equal for both cell lineages. Cell-cell adhesion The terms v N i, adh and v P j, adh describe homotypic and heterotypic adhesive dynamics due to the activation of selected cadherin molecules. Their form is the sum of pairwise interactions analogous to the repulsive counterpart, i.e., where f αβ adh : are the corresponding kernels. The explicit form of f αβ adh can be established by assuming that a pair of cells must be sufficiently close to adhere. Specifically, we assume adhesions start to form if their mutual distance is below the maximal extension of membrane protrusions, d a (see Fig. 2). Further, adhesion increases as the contact surface between two individuals grows, which follow as the mutual distance drops. Taking all of these considerations into account, for any couple αβ with α, β ∈ {N , P}, we opt for the following family of functions: where F αβ a > 0 is a measure of the maximal adhesive stimulus, occurring for highly overlapping cells (dashed blue curve in the left panel of Fig. 3).

CIL velocity component The velocity contributions v N P
i, CIL and v P N j, CIL describe heterotypic CIL responses, resulting in mutual cell repulsion. To translate this mechanism into mathematical terms, we first note that the CIL response activates when a neural crest, say i, and a placodal cell, say j, come into contact. In other words, it occurs when their relative distance d N P i j = d P N ji falls below the critical value d c which, as previously explained, gives a measure of the cell diameter (see again Fig. 2), assumed equal for both cell lineages . The durations of the subsequent phases of reflected directional migration are characterised by decay times, say τ N P for NCs and τ P N for PCs. The CIL velocity components v N P i, CIL and v P N j, CIL are therefore assumed to evolve according to: In particular, ω N P and ω P N quantify the size of the instantaneous push, which could (amongst other factors) be dictated by both the number of heterotypic N-cadherin adhesive bonds and the level of activation of the downstream intracellular cascade. u N P i and u P N j define the subsequent direction of the CIL-induced migration that follows re-polarisation. For a two cell scenario, these are aligned with the unit vectorr N P 11 : in particular, u N P 1 = −r N P 11 and u P N 1 = −r P N 11 . In a multicellular scenario, we have instead to account for the possibility of multiple simultaneous heterotypic contacts. Specifically, for h = 1, . . . , n α with α, β ∈ {N , P} and α = β, u αβ h is taken to be the weighted sum of the unit vectorsr αβ hk that connect the h-th cell of type α with the surrounding k-th particle of type β: where ( · ) + = ( · )+| · | 2 , and d αβ hk andr αβ hk have been defined in Eq.
(3). The positive part function ( · ) + is necessary in Eq. (8) to ensure that the CIL response is triggered only when d αβ hk ≤ d c . The above equation implements the assumption that the closest agents are those that generate the largest contribution when establishing the CIL velocity direction component: this is a reasonable hypothesis, since nearest individuals are likely to form a larger number of more extended cadherin complexes, thereby triggering a more significant CIL response.
The adhesive coefficients F N P a and F P N a included in Eq. (7) account for the fact that the intracellular cascade underlying CIL is mediated through N-cadherin bonds. This allows us to then implement selective experimental manipulations, such as the inhibition of N-cadherin adhesive complex bond formation (by setting F N P a = F P N a = 0 µm s −1 ) or the disruption of the downstream intracellular cascade regulating cell re-polarisation (by setting ω N P = ω P N = 0).
Finally, we assume that for the duration of a CIL contribution between the i-th NC and the j-th PC, i.e., v N P i,CIL and v P N j,CIL , adhesive interactions between the two individuals are neglected (i.e., v N i,adh = v P j,adh = 0), even given low intercellular . From a biological point of view, this is reasonable given that CIL triggers the collapse of the adhesive complexes in addition to repolarization of the cell cytoskeleton, see Steventon et al. (2014), Szabo and Mayor (2015), Theveneau et al. (2013).

Chemotaxis velocity component
where d N P i, min := min k=1,...,n P d N P ik . The function f N P dist : R + → [0, 1] is defined as follows as represented in Fig. 3 (right panel). The above law reduces chemotactic movement of NCs as they start to undergo adhesive interactions with at least one PC (i.e., if d N P i, min ≤ d a ) and ensures it becomes negligible for close contact (i.e., if d N P i, min ≤ d c ). Biologically, this could reflect the intracellular integration of chemotactic/adhesion pathways, where the positive promotion of Rac1 (and hence cell protrusions) by chemotaxis is negated following N-cadherin binding ). In Eq. (9), χ N represents a chemotactic sensitivity parameter, e.g., the functionality of Sdf1 membrane receptors Cxcr4, and is taken to be equal for all NCs. To avoid unrealistic cell speeds, v N i, chem is capped by parameter v max , denoting the measured maximal speed of NCs.
Random velocity component The terms v N i, rand and v P j, rand account for the isotropic fluctuations that impact on cell trajectories, as observed in biological experiments by Theveneau et al. (2013) where θ α h (t), with α ∈ {N , P}, is a random angle uniformly distributed over [0, 2π). v N rand and v P rand are constant speeds, estimated accounting for observations in in vitro assays in Theveneau et al. (2013). Of course, while other more sophisticated models for incorporating random movement could be considered, the lack of any detailed biological data motivates us to consider the above form for its simplicity.

Sdf1 kinetics
In addition to Eq. (1), we require an equation for Sdf1 kinetics. In particular, the chemical substance is assumed to be secreted (at a constant rate) by all placode cells, to diffuse through the surrounding environment and to degrade at a constant rate. Its spatiotemporal evolution can therefore be modelled by the following classical reaction-diffusion (RD) equation: where D is a homogeneous diffusion coefficient and ε and c 0 are respectively the decay and production rates. The functions δ x P j (t) describe point sources, centred on the PC positions.

Simulation details and parametrisation
All numerical simulations are performed on a bounded square domain of dimensions 700 µm×700 µm, based on the biological images reported in Theveneau et al. (2013). If a cell i touches the boundary ∂Ω at time t ∈ [0, T max ], it is assumed to stop (i.e., v α i (t) = 0); Sdf1 is taken to be absorbed at the boundary, i.e., c(x, t) = 0, ∀ x ∈ ∂Ω and ∀ t ∈ R + . At the start of each simulation we assume zero Sdf1 (i.e., c(x, 0) = 0, ∀ x ∈ Ω) and neglect CIL-related velocity components (i.e., v N P i,CIL (0) = v P N j,CIL (0) = 0, ∀ i = 1, . . . , n N and j = 1, . . . , n P ). Each in silico experiment runs until a characteristic observation time, T max , set at 300 min = 5h. Initial cell distributions are specified for each simulation suite.
While many parameters can be directly estimated from biological data, some have no clear experimental correspondence. Their estimate is, rather, based on a process of sensitivity analysis supported by empirical measurements and observations. Default parameters used in simulations are listed in Table 1, with those obtained via data fitting indicated. Precise details for parameter estimation are given below.
Cell sizes From the biological literature, NCs and PCs appear to have similar morphology. In particular, an analysis of the experimental images in Theveneau et al. (2013) allows us to fix a common diameter of the perinuclear region d r = 20 µm and a maximal filopodia extension at 30 µm, so that d a = 60 µm. The mean body diameter of fully adherent cells (d c ) is then set equal to the intermediate value, i.e., d c = 30 µm.
Parameters related to Sdf1 kinetics and chemotactic velocity components Values for coefficients relevant for Sdf1 kinetics in Eq. (12) are set according to those used in : D = 13.3 µm 2 s −1 , ε = 0.0004 s −1 and c 0 = 0.027 µM s −1 . The maximal chemotaxis velocity v max is instead set at 120 µmh −1 = 0.0333 µm s −1 , according to the biological measurements reported in Theveneau et al. (2010). Finally, the chemotactic sensitivity coefficient χ N is set at 500 µm 2 µM −1 s −1 , on the basis of Colombi et al. (2015), where chemotaxis is stimulated by a chemical factor (VEGF) with a similar diffusion coefficient. Model responses to independent variations of c 0 and χ N will be described in the following section.
Parameters related to cell resistance to compression and cell-cell adhesion As far as we are aware, there is no evidence that placode and neural crest cells have a different nuclear stiffness: therefore, we set F αβ r = F r for any α, β ∈ {N , P}. Regarding cell-cell adhesion, we assume that N-cadherin-N-cadherin bonds generate symmetric adhesive responses in each of the adjacent cells, i.e., we take F N P a = F P N a = F N N a , since both cell lines express N-cadherin molecules. PCs additionally express Ecadherin, thereby generating further strong and stable homotypic bonds. Consequently, it is reasonable to take Beyond these biologically motivated considerations, the estimates of F r and F αβ a must still be treated with care. Specifically, regions of the parameter space exist that can generate a physically unrealistic system evolution, such as "population collapse" in which the individuals condense into an infeasibly close aggregate. Consequently, the set of interaction coefficients must be chosen to ensure that a relaxed configuration can be attained and maintained, characterized by fixed and finite mutual distances: a so-called crystalline pattern. In this respect, it has been shown that the large-time behaviour of particle systems subjected to non-local pairwise interactions is related to the H-stability of the relative kernels and potentials, see Ruelle (1969). In particular, taking advantage of the characterization of H-stable potentials provided in Ruelle (1969), recent works (such as Cañizo et al. 2015;Cañizo and Patacchini 2018;Carrillo et al. 2018) provide a criterion that determines a subregion of the parameter space of the repulsive-adhesive interacting kernels that results in realistic crystalline cell pattern. Accounting for these analytical results, for any pair αβ ∈ {N N, N P, P N , P P}, the parameters F r and F αβ a must satisfy the relation For the above specified cell sizes, this becomes F r /F αβ a > 21.6. Despite its usefulness in preventing unrealistic situations, this criterion still does not indicate the exact pairing of interaction parameters. In this respect it is then useful to remark that cells are unable Fig. 4 Estimate of the heterotypic cell-cell interaction parameters. Evolution of the relative distance and the speed of two cells, of type N and P respectively, initially located at a distance equal to 59 µm, with dynamics regulated by only the interaction velocity component. It is possible to observe that the system behaviour is robust with respect to variations in F r , but sensitive to variations in F N P a = F P N a , which set the time needed to establish crystalline configurations. Note that the blue and dashed red curves, as well as the yellow and dashed purple curves, are indistinguishable (color figure online) to move too rapidly, i.e., |v α res (t) + v α adh (t)| should not exceed reasonable values (e.g., v max ). Preliminary simulations where we switch off all velocity components apart from the interaction terms are therefore performed to highlight the role of cell-cell repulsion/adhesion on system behaviour, for the purposes of parameter estimation.
Specifically, we first consider a two cell system composed of a NC and a PC, in order to focus on heterotypic cell-cell interactions. The NC individual is placed at (200 µm, 350 µm) and the PC at (259 µm, 350 µm), i.e., such that d N P (0) is fractionally below d a . The time evolution of the system is evaluated through d N P 11 (t) and the speeds |v N 1 (t)| = |v P 1 (t)|, which are equal for any t ∈ [0, T max ] as a consequence of our assumption of symmetric pairwise interactions. The results in Fig. 4 show that the behaviour does not change significantly with variations in F r : see the overlap between the first and second curves, or between the third and fourth. Conversely, the system is highly sensitive to variations of F N P a = F P N a , which define the timescale required to establish a crystalline configuration in the multicellular case. We further remark that cell speeds remain plausible under all considered parameter settings.
We next turn to the role of homotypic cell-cell interaction parameters on model outcomes. We consider the behaviour of a colony of 100 PCs for three different cell-cell interaction coefficient parameter settings, each satisfying the H-stability condition: F r = 1 µm s −1 with F P P a = 0.02 µm s −1 ; F r = 1 µm s −1 with F P P a = 0.008 µm s −1 ; F r = 0.4 µm s −1 with F P P a = 0.008 µm s −1 ; and F r = 1 μm s −1 with F P P a = 0.002 µm s −1 . All simulations start from the same initial condition: an almost round cluster of radius ≈ 120 µm placed in the centre of the domain, see Fig. 5a. To focus on adhesive/repulsive stimuli, we recall that we switch off the random velocity component (i.e., v P j,rand (t) = 0, for any j = 1, . . . , n P with t ∈ [0, T max ]) and note that the CIL velocity components will be zero for each PC (i.e., v P j,CIL = 0 for any j = 1, . . . , n P with t ∈ [0, T max ]), due to the absence of NCs. As shown in Fig. 5b-e, in all cases the PC colony generates and maintains a crystalline configuration (without collapse or dispersion), consistent with our parameter selection within the H-stability region of the system. On fixing F r , decreasing F P P a (see Fig. 5b-c, e) generates a more dispersed colony, since reducing adhesive contributions allows further emergence of repulsive effects. The most dispersed configuration we can obtain under adhesive-repulsive interactions alone arises by setting F P P a = 0 µm s −1 , and is characterised by bounded minimal interparticle distances that do not trespass d r due to the definition of f αβ res in Eq. (4) (see Cañizo et al. 2015;Cañizo and Patacchini 2018;Carrillo et al. 2018 and references therein for further details). On the other hand, on fixing F P P a we observe that decreasing F r results in a slight decrease of the equilibrium interparticle mean distance, specifically towards a highly packed colony (with a small degree of cell membrane overlap, see Fig. 5c, d). However, by considering parameter pairs that satisfy the Hstability condition in Eq. (13), a minimal interparticle distance is preserved and finite. It is still necessary to control whether selected parameters pairs generate interparticle distances not significantly lower than d r . In this respect, in Fig. 5 cells are represented by circles of diameter d r to highlight that a sensible minimal intercellular distance is maintained, and we keep such a representation in subsequent figures to demonstrate this important property. Given these experiments, we set F P P a = 0.008 µm s −1 and F r = 1 µm s −1 (Fig. 5c): for consistency with experimental observations, we require the maintenance of a highly packed PC colony, but without cell overlap.
The same study is then performed for the NC aggregate, starting from an equivalent initial condition (Fig. 5f). Fixing F r = 1 µm s −1 (as for the PCs), from our preliminary simulations we set F N N a = 0.002 µm s −1 to generate a more dispersed NC aggregate with respect to the PC cluster (compare panels (c) and (e) in Fig. 5), in accordance to in vitro experiments in Theveneau et al. (2013). According to the assumption that F N N a = F N P a = F P N a , we finally set F N P a = F P N a = 0.002 µm s −1 .
Parameters related to CIL velocity components For realistic estimates of the parameters τ αβ and ω αβ , with αβ ∈ {N P, P N }, that characterise the CIL velocity, we start by considering the experimental literature (see, among others, Fig. 6 (b) in Theveneau et al. 2013 for heterotypic CIL responses). In particular, regardless of cell phenotype, we assume that the Rac1-dependent intracellular cascade that is responsible for the CIL mechanism is activated for some given time such that τ N P = τ P N = 6 min, since we estimate the overall CIL contribution to last ≈ 12 min. Then, accounting for images in Theveneau et al. (2013) relative, as seen, to heterotypic dynamics, we can assume that the instantaneous push resulting from CIL is dictated by the intrinsic migratory abilities of the re-polarized cell: since NCs are known to have more significant motility than PCs, we therefore state ω N P > ω P N .
To estimate their values, we analyse how a two-cell system composed of a single neural crest and a single placode cell evolves upon variations of ω N P and ω P N . To highlight CIL velocity components we switch off the random term while maintaining other parameters as previously estimated (and listed in Table 1). Further, for an accurate fit between the numeric and experimental measures in Scarpa et al. (2013) and Theveneau et al. (2013), we use the following critical quantities for comparison: (i) the cell mutual distance evaluated 12 min post NC-PC collision (hereafter denoted by d CIL ); and (ii) the maximum velocity reached by the two individuals as a consequence of CIL (hereafter denoted by v CIL ). From the data generated by the biological experiments, we can assume that reasonable values for ω N P and ω P N should result in d CIL ≈ 50 µm and v CIL < 3 µm min −1 = 0.05 µm s −1 . From the plot of Fig. 6 we observe that such empirical data can be reasonably matched by setting ω N P = 5000 and ω P N = 2500. Fig. 6 Estimates of the parameters characterising the CIL velocity components. Cell intercellular distance d N P 11 (t) = d P N 11 (t) and speed of the two interacting individuals in three parameter settings: i.e., ω N P = 5000 with ω P N = 2500 (blue case); ω N P = 7500 with ω P N = 5000 (purple case); and ω N P = 10,000 with ω P N = 7500 (yellow case). Squares in the top panel indicate the value of d CIL for each case (color figure online) Fig. 7 Estimate of parameters related to random velocity components. a-c Final configuration (i.e., at t = 300 min) of the NC colony for different settings of v N rand . d Final configuration (i.e., at t = 300 min) of the PC colony for the selected value of v P rand Parameters related to random velocity components We first recall from Eq. (11) that θ N (t) and θ P (t) are random variables uniformly distributed on [0, 2π). The values of v N rand and v P rand instead denote the intensity of Brownian fluctuations. Analysing the experimental literature (see, among others, Fig. 2(a, b) of Theveneau et al. 2013), we can observe that NCs have significantly greater autonomous intrinsic motility (i.e., isotropic oscillations) than PCs. Hence we set v N rand v P rand and perform a series of simulations, starting from the same initial conditions used to estimate homotypic interaction parameters (see Figs. 5a, f). In all cases, cell dynamics are regulated only by homotypic interactions and random velocity contributions and all other contributions are subsequently neglected.
Simulation results under variations of v N rand (Fig. 7a-c) are compared against their experimental counterpart: the best-fit (measured in terms of cell displacement and intercellular distance) emerges for v N rand = 0.6 µm s −1 . The dynamics of placode cell clusters are highly determined by the strong homotypic E-cadherin interactions that maintain cluster shape, and we observe little change upon alteration of v P rand . We therefore choose a value v P rand = 0.05 µm s −1 and remark that simulations are relatively insensitive against variations in this parameter.

Two-cell system
We begin with a two-cell system, i.e., based on the interactions between a single NC and a single PC. We first perform a reference simulation that highlights "normal" dynamics, before replicating the impact of selected experimental disruptions. In each simulation, we initially space the cells 70 µm apart: in this respect, x N 1 (0) = (230 µm, 350 µm) and x P 1 (0) = (300 µm, 350 µm) define the cell starting locations, see Fig. 8a and Fig. 9a. To assess the evolution of the system under various scenarios, our primary measurement will be the intercellular distance, d N P 11 , which is defined as the distance between the point positions marking cell centres and is reported in the left panel of Fig. 10.
Reference simulation As observed from time-lapse images in Fig. 8 (and from Supplementary Movie M1), random wandering dominates until the Sdf1 field forming around the PC reaches the NC. Subsequently, the NC chemotactically migrates towards the PC, generating the chase phase of the process (see Fig. 8c). As the intercellular distance "Chase-and-run" dynamics in the reference simulation. a-c Blue and red circles and lines represent positions and trajectories over [0, t] for the NC and PC, respectively. The diameter of the solid circles reflects d r , while outer circles represent d c : we highlight how CIL is triggered as the intercellular distance drops below d c . Note that, for clarity of presentation, we select only a central portion of the overall domain. d Representation of angles θ N (solid blue angle) and θ P (dashed red angle) for cell direction of motion before and after an NC-PC contact activates the CIL response. The solid light blue and dashed pink arrows represent, respectively, NC and PC direction of motion before contact, i.e., where t c denotes the instantaneous time of the cell-cell collision and Δt = 9 min as in Theveneau et al. (2013). Analogously, the solid blue and dashed red arrows indicate NC and PC direction of motion after contact, i.e., drops below d a , adhesive interactions occur and the two individuals further approach each other, with their mutual distance narrowing accordingly ( Fig. 8d and yellow curve of Fig. 10). When the intercellular distance drops to d c , CIL is triggered (see Fig. 9b and Fig. 10). CIL then induces the NC and PC to "bounce" away one from each other ( Fig. 8e and Fig. 9c), moving apart until quasi-stabilising at a distance dictated by the magnitude of the CIL response (Fig. 9c).
A return to the initial Brownian motion is then observed and the show plays on repeat, see Fig. 8f-h, with chemotactic-driven chase drawing cells together until CIL induces a further bounce.
The post-contact directionality of the cells, shown in Fig. 9d-f, lies in qualitative agreement with the experimental literature (see, for instance, Fig. 6 (c) of Theveneau et al. 2013). Specifically, in Theveneau et al. (2013) angles were calculated using three positions for each of the NC and PC, corresponding to 9 min before a contact, at a contact, and 9 min after a contact. Changes in the direction of motion of the NC and PC, following CIL activation, are thus quantified here by the angles θ α (with α ∈ {N , P}) represented in Fig. 9d, where the solid light blue and dashed pink arrows denote, respectively, the NC and PC direction of motion before contact, while the solid blue and dashed red arrow indicate NC and PC direction of motion after contact. Notice that the measure θ α and the approach used to calculate the corresponding angle measured in Theveneau et al. (2013), for α ∈ {N , P}, are consistent. To account for the random component to cell dynamics, means and variances of θ N and θ P are inferred following 100 realisations of the reference test, see Fig. 9e-f, and the resulting values are consistent with experimental measurements reported in Theveneau et al. (2013).
Overall, there is a net directional cell migration, as captured by the substantial shift in the cell positions from the beginning to end of the observation time, see Fig. 9g. However, the global displacement is strongly affected by the random velocity component, as indicated by the fluctuations in d N P 11 (t) (see the yellow curve in the left panel of Fig. 10). For the sake of completeness, cell paths along the x-axis are evaluated in the absence of Brownian fluctuations: as demonstrated in the right panel Disrupting NC chemotaxis migration We next target the capacity of the NC to chemotactically migrate, by either (i) inhibiting the Sdf1 production from the PC individuals (setting c 0 = 0 µM s −1 in Eq. (12)), or (ii) switching off the chemotactic sensitivity term (setting χ N = 0 µm 2 µM −1 s −1 in Eq. (9)), corresponding to a knock-out of Cxcr4 receptor activity. In both cases (see Fig. 11a, b), the NC is unable to sense the PC and therefore the two cells exhibit uncorrelated random crawling about their initial positions: the intercellular distance hovers about (and can also exceed) its initial value over the simulation timecourse (light and dark green lines in the left panel of Fig. 10). We note that crawling is more pronounced for the NC individual, as expected given that the choice v P rand v N rand in Eq. (11) leads to |v P 1,rand | |v N 1,rand | for any t ∈ R + . Note that while random wandering can potentially bring cells within contact range, and hence induce CIL, this is not typically observed within the simulation timecourse.
Disrupting N-cadherin bonds We then inhibit N-cadherin adhesive bond formation through setting to null the heterotypic adhesive strengths, i.e., F N P a = F P N a = 0 µm s −1 in Eqs. (6)-(7). This disruption acts to shut down adhesive and CIL velocity contributions, thus v N 1,adh (t) = v P 1,adh (t) = v N 1,CIL (t) = v P 1,CIL (t) = 0 for all in Sect. 2). Typical trajectories are shown in Fig. 11c and we chart the intercellular distance and the directional cell displacement in the left panel of Fig. 10 (light blue curve). We observe that the chemotactic-driven chase occurs as normal, with the two cells stabilising at a distance d N P The subsequent run phenomenology is no longer obtained and disruption of the adhesive bonds allows cells to approach closer than d c = 30 µm without bouncing back. This lies in contrast to the reference case, where normal chase behaviour is maintained only for intercellular distances d N P 11 (t) ≥ d c and below this threshold the run phase is invoked. These observations lie in accordance with biological data (see, for instance, Fig. 6 (b) in Theveneau et al. 2013). Fluctuations in intercellular distance are now mainly dictated by the compression-driven repulsive velocity contributions, which enter when the intercellular distance falls below d r = 20 µm. Note that the small perturbations in cell positions stem again from the random velocity components.
Disrupting the CIL response Our final disruption targets directly the intracellular cascade responsible for CIL, setting ω N P = ω P N = 0 in Eq. (7). Consequently, only CIL velocity contributions are neglected, i.e., v N P 1, In particular, we remark that adhesive interactions still occur through the formation of N-cadherin adhesive bonds. Resulting dynamics are captured by cell trajectories (see Fig. 11d) and by the time evolution of intercellular distances (dark blue curve in left panel of Fig. 10). As for disruption of N-cadherin adhesive bonds, the two cells now stabilise at distances d N P 11 (t) ∈ [d r , d a ] without bouncing back. This highlights how an inhibition acting directly and only on the Rac1-dependent CIL pathways (i.e., without affecting the activity of the N-cadherin molecules) is sufficient to disrupt normal "chase-and-run" behaviour, again in noteworthy agreement with the corresponding experiments . A subtle distinction arises in cell behaviour due to disruption of either the Rac1-dependent intracellular cascade or the N-cadherin adhesive bonds, however, and it can be captured by quantifying the magnitude of intercellular distance fluctuations. Specifically, we have smaller fluctuations of d N P 11 (t) in the former scenario due to the normal activity of adhesive interactions for d N P 11 (t) < d a , which continue to exert a degree of control over random and repulsive velocity contributions.

Multicellular populations system
We next explore "chase-and-run" collective behaviour in a multicellular system, formed by a population of neural crest cells that interact with a placode aggregate. As above, we first consider a reference simulation for the complete model to highlight quantitative and qualitative determinants of the coordinated cell movement. We subsequently describe system phenomenology under selected mechanistic disruptions. For all experiments we consider 200 cells, composed of equal numbers of NC and PCs (i.e., n N = n P = 100), with the same initial set-up: a culture formed from two adjacent aggregates, one of NCs and one of PCs, each arranged in a quasi-round cluster of radius ≈ 120 µm, see Fig. 12a. This initial configuration is similar to the experimental counterpart in Theveneau et al. (2013). As critical measures to quantify system behaviour, we will take the minimal homotypic and heterotypic distances between the generic i-th cell belonging to population α ∈ {N , P} and the other individuals, i.e., d αβ i,min = min j=1,...,n β d αβ i j , with β ∈ {N , P}.
Reference simulation As illustrated in Fig. 12 (and in Supplementary Movie M2), neural crest cells initially undergo collective chemotaxis in the direction of the placode colony, maintaining a quite compact configuration (which mainly stems from the common directional velocity rather than from the low homotypic adhesive interactions). As NC cells at the leading edge of the cluster approach the PC cluster, heterotypic CIL responses are initiated and PC individuals subsequently move away from the NCs, see Fig. 12c. Reverse movements of NC cells, however, are blunted by their ongoing chemotactic responses and by the adhesive interactions (albeit low). As a result, we observe a consistent and collective directional motion of the entire system: the PC cluster is gradually shunted towards the right, pushed and pursued by the NC aggregate, see Fig. 12d. The collective cell phenomenology captured by our simulations (summarised by the trajectories in Fig. 12e-f) lie in general accordance with corresponding in vitro experimental dynamics (cf. Fig. 2(c)-(h) in Theveneau et al. 2013), although we remark that in our numerical simulations the chasing NC cluster becomes increasingly crescent-like in shape. As we show below, the conformation of the NC aggregate subtlety changes upon parameter variation and we will further return to this aspect in the discussion.
Disrupting NC chemotaxis migration Our first virtual knock-out focuses on NC chemotactic movement. We first vary Sdf1 production by the PCs, i.e., we vary c 0 in Eq. (12). Given sufficiently high secretion rates of the chemical (≥ 0.01 µM s −1 ), the NC Fig. 13 Effect of variable Sdf1 production by the PCs. a-f Final configurations, i.e., at t = 300 min. g Displacements of the centres of mass for the two aggregates: the reference value c 0 = 0.027 µM s −1 used in Fig. 12 is written in bold type. h-k The boxplots report the minimal intercellular distance between cells, from left to right: h between a NC and other NCs; i between a NC and the PCs; j between a PC and the NCs; and, k between a PC and other PCs aggregate forms a compact aggregate that chases the PCs, as in the reference scenario, see Fig. 13d-f. For lower values of c 0 (i.e., < 0.01 µM s −1 ), however, the trailing region of the NC colony instead scatters, see Fig. 13a-c. In this case a comet-like tail forms from individual cells that have lost contact with the main mass, due to not receiving a sufficiently high stimulus of the diffusing chemoattractant. In such cases, the collective directional movement of the system is slightly downregulated, probably as a consequence of the lower pushing force that the dispersed NC cluster exerts on the PC island. Notably, however, the formation of such tails can be observed in experiments (cf. Fig. 2 (c) in Theveneau et al. (2013) and its accompanying movies). Finally, and consistent with the two-cell scenario, complete inhibition of Sdf1 production (i.e., c 0 = 0 µM s −1 ) leads to abolition of chase-and-run and the NC and PC cells/colonies crawl about their initial position, as shown in Fig. 13a and in the first column of Fig. 16. Here, NCs (as a whole) do not perceive the PCs and subsequently move randomly, independently with respect to the PCs (with the exception of those initially located close to the PC cluster such that random wandering brings them within adhesive range).  . a-f Final configurations, i.e., at t = 300 min. g Displacements of the centres of mass for the two aggregates: the reference value χ N = 500 µm 2 µM −1 s −1 used in Fig. 12 is written in bold type. h-k Boxplots report the minimal intercellular distance between cells, from left to right: h between a NC and other NCs; i between a NC and PCs; j between a PC and NCs; and, k between a PC and other PCs Addressing Fig. 13h-k, it can be observed that the production of c 0 does not influence homotypic placode intercellular distances, since d N N i,min and d P P i,min are predominantly regulated by adhesion. Neither is the heterotypic distance d P N i,min , since (excluding the complete knock out case) there are always NC cells that reach the PC aggregate, see Fig. 13j. The heterotypic distance d N P i,min , however, presents a more scattered distribution, with a significant number of outliers emerging for low values of c 0 , see Fig. 13i. This phenomenology accompanies the above described trail and results from single cells that have detached from the main NC aggregate, subsequently subjected to greater degrees of random migration.
We now keep Sdf1 production constant and vary instead the expression/activity of NC Sdf1 receptors (i.e., varying χ N in Eq. (9)). As shown in Fig. 14f, high values of χ N (i.e., > 100 µm 2 µM −1 s −1 ) result in the compact chemotactic migration of the NC cluster and CIL-driven coordinated net displacement of the two populations, as in the reference case. A slight downregulation of the NC chemical response (i.e., χ N ∈ [50, 100] µm 2 µM −1 s −1 ) instead promotes dispersion of NCs into the trailing region of the cluster and significant inhibition of the run phase of the process, see Fig. 14c-e. Significantly low values of χ N (say 25 µm 2 µM −1 s −1 ) completely disrupt NC Fig. 15 Effect of blocking Sdf1 production in a given fraction of the PC population. a-f, b1-b4 Final configurations, i.e., at t = 300 min. In particular, a corresponds to the complete knock-out of Sdf1 production in all cells, whereas f is the reference case reported in Fig. 12 where all PCs secrete Sdf1. g Displacement of the center of mass of the two aggregates. b1-b4 Final configurations under specific localisation of the (10%) PCs that produce Sdf1: b1 10% positioned within the north edge of the PC colony; b2 5% positioned within each of the north and south edges; b3 10% positioned within the east edge; and b4 10% positioned within the west edge chemotactic movement and there is negligible overall directional locomotion of the system, see Fig. 14a-d. As highlighted by the first two columns of Fig. 16, complete knock-out of cell sensitivity has effectively the same impact on cell trajectories as a complete knock-out of Sdf1 production. Boxplots for these investigations are similar to those observed following Sdf1 production variation, compare Fig. 13h-k and Fig. 14hk. We do remark, however, that there is a clearer dependence of the heterotypic minimal intercellular distance between an NC and PCs (see d P N imin (t F ) in Fig. 14j), with outliers present only in the transition region when the NC trail forms.
We finally predict the effect of completely blocking Sdf1 production, but only in a given fraction of the placode population. As shown by Fig. 15c-f, chase and run behaviour (as measured by the group displacement) is relatively unaffected under moderate disruptions but becomes significantly reduced when only few placode cells secrete the chemical substance (i.e., 10%), see Fig. 15a, b. It is interesting to note that reducing the number of PCs producing Sdf1 also generates a morphological transition in the NC aggregate, which morphs from a "crescent moon" shape encircling the PCs to an ellipsoid configuration adjacent to the PC cluster (see Fig. 15f) and, finally, a more dispersed cluster accompanied by the comet-tail (see Fig. 15a).
For the simulations reported in Fig. 15a-f, Sdf1 expressing PCs are randomly positioned in the cluster. We next investigate how their localisation within the PC colony impacts on the chase-and-run phenomenon. Specifically, in Fig. 15b1-b4 we consider four distinct cases for the scenario under which only 10% of the PC population produces Sdf1. In Fig. 15b1, all Sdf1-expressing cells are placed within the north edge of the PC colony; in Fig. 15b2, half are placed within the north edge of the PC colony and the remainder within the south edge; in Fig. 15b3, all Sdf1-expressing PCs are within the east edge of the group; and, in Fig. 15b4, all are positioned within the west edge of the PC aggregate. Under the first scenario, the induced asymmetrical pattern of Sdf1 results in NCs moving towards the northern edge of the PC colony, disrupting both the symmetric comet-tail distribution shown in Fig. 15b and the displacement of the PC colony. In Fig. 15b2, NCs organise into a crescent moon that envelopes the PC colony, being attracted towards both the north and the south edges of the aggregate. This results in a slight reduction of net displacement (with respect to Fig. 15b, where PCs producing Sdf1 are randomly distributed). The comet-tail configuration observed in Fig. 15b is preserved in the other cases ( Fig. 15b3-b4) but there are distinct changes to the net migration. For Sdf1-expressing PCs concentrated at the eastern edge, the final distribution is similar to that obtained in Fig. 15a, i.e., in the absence of Sdf1 producing cells, and there is little net migration: NCs are distant from the source of Sdf1 production and minimal chemotactic migration/chase-and-run ensues. On the other hand, see Fig. 15b4, localising the secreting PC population to the west edge allows NCs to organize into a compact aggregate that stimulates chase-andrun, with a displacement comparable to the levels observed in Fig. 15b-c. Thus, while precise positioning of the secreting cells may act to significantly disrupt chaseand-run/cell displacement, any benefits (in terms of increasing the net displamcent) from optimal localisation are marginal (in respect to the randomly distributed case of Fig. 15b). Moreover, we remark that the importance of localisation within the PC colony decreases as the fraction of Sdf1 producing PCs increases.
Disrupting N-cadherin bonds We next inhibit the formation of N-cadherin bonds, setting adhesive strengths F N N a = F N P a = F P N a = 0 µm s −1 while keeping F P P a = 0.008 µm s −1 in Eq. (5) and in Eq. (7). The non-zero PC-PC adhesive parameter stems from their additional expression of E-cadherin. Consequently, we have v N N i,adh (t) = v N P i,adh (t) = v P N j,adh (t) = 0, as well as v N P i,CIL (t) = v P N j,CIL (t) = 0 with i = 1, . . . , n N and j = 1, . . . , n P , for any t ∈ [0, T max ]. Other velocity components and parameters remain the same as the reference simulation. As observed in the third column of Fig. 16 (see also Supplementary Movie M3), the NC aggregate continues to move via chemotaxis towards the PCs; the two populations subsequently form a single quasiround cluster, with a central region characterised by mixing of the two cell lineages. The absence of N-cadherin interactions precludes activation of a CIL response and disrupts the "run" phase of the dynamics. This lies in qualitative agreement with the corresponding biological experiment (cf. Fig. 6 (j) in Theveneau et al. (2013)). However, we remark that the in vitro setting has a three-dimensional element where, despite the 2D plating, a degree of cell overlap occurs at the interface between the two aggregates. This aspect is precluded by our numerical realisation, a consequence of the strictly two-dimensional domain used as a simplification.
Disrupting the CIL response We next target the intracellular cascade responsible for the CIL response, i.e., we set ω N P = ω P N = 0 in Eq. (7). Hence, v N P i,CIL (t) = v P N j,CIL (t) = 0 with i = 1, . . . , n N and j = 1, . . . , n P , for any t ∈ [0, T max ], while the other velocity components/parameters remain as the reference simulation. The system behaviour, captured in the last column of Fig. 16, is analogous to the case of disruption Fig. 16 Final configuration, i.e., at t = 300 min, of NC and PC colonies and cell trajectories in the case of complete knock-outs, namely: inhibition of Sdf1 production by the PCs (leftmost column); inhibition of expression/activity of NC Sdf1 receptors (second column); inhibition of N-cadherin adhesive bond formation (third column); and, disruptions of the Rac1 dependent intracellular cascade responsible for the CIL response (rightmost column) of N-cadherin bonds. Chemotaxis allows the NC aggregate to chase the PC aggregate, as above, and the two populations form a single aggregate in a similar manner to the previous setting, and in qualitative agreement with biological experiments (e.g., see Fig. 6 (j) in Theveneau et al. 2013).

Variation in PC E-cadherin expression
We again use our approach in a predictive manner. In particular, we test system dynamics upon variations in the PC E-cadherin expression/activity (evaluated by parameter F P P a ), which, as far as we are aware, has not been systematically explored from an experimental perspective. As in previous perturbation-type experiments, all other model components and coefficients are kept according to the reference case setting. As shown in Fig. 17, PC homotypic adhesiveness strictly correlates with the morphological characteristics of the system. For large values of F P P a (i.e., F N N a = F N P a = F P N a = 0.002 µm s −1 ), the PC cluster behaves in the manner of a quasi-rigid disk and becomes surrounded by the NC aggregate, see Fig. 17c-e. As PC homotypic adhesion is lowered to a similar magnitude to other adhesion contact strengths (see Fig. 17b) both cell clusters undergo coordinated movement while maintaining a quasi-elliptical shape. Finally, when placode self adhesion becomes substantially reduced (i.e., to the point that F P P a F N N a = F N P a = F P N a = 0.002 µm s −1 ), we observe a capacity of the NC population to infiltrate and disperse the PC cluster into an encapsulating crescent-like cluster, see Fig. 17a. We further remark that within the simulation timescales, net directional migration of the system is not significantly affected by variations to PC homotypic adhesiveness.
Variation in the extension of the range of repulsive dynamics We finally provide predictions on how cell dynamics are altered following variation of cell repulsive behaviour.
In particular, an increment in the value of d r can be viewed as a numerical counterpart to an increment in the stiffness of the cytosolic region that surrounds the cell nucleus. The cell cytoskeleton is, in fact, formed by a contractile acto-myosin cortex whose rigidity can be modulated through chemical manipulation. As expected, and highlighted in Fig. 18a, b, under larger values of d r the homotypic clusters become more extended/dispersed. As revealed in Fig. 18c Fig. 18d, e, the heterotypic distances d N P i,min and d P N i,min (effectively measuring the distance between the two cell aggregates) are almost insensitive to variations in d r . This arises from the fact that heterotypic repulsion is, rather, predominantly regulated by CIL, which Numerical simulation with homotypic CIL between NCs, highlighting the excessive scattering of the NC colony with respect to experimental observations. In this case, τ N N = τ N P and ω N N = ω N P ) does not allow NC individuals to become too close to PCs. We note that again, in this case, no other changes have been made to model components from the reference case scenario.

Discussion
We have explored "chase-and-run" dynamics in heterogeneous cell populations. Specifically, we have encoded a computational model for the evolution of a system composed of plated neural crest and placode cells. Chemotaxis of NCs in response to a placode-secreted attractant triggers the chase phase and a subsequent CIL response generates the run phase of the process. By limiting to a purposefully small set of assumptions and replicating the in vitro set-up of two populations initially placed in juxtaposing aggregates, our computational simulations have been able to reproduce the experimental findings in Theveneau et al. (2013) in semi-quantitative fashion, i.e., to capture the productive net migration of the cell clusters. Yet, some deficiencies of this minimalist approach are raised by our simulations. We use this discussion to further elaborate, and discuss subsequent extensions to be considered in future modelling.
First, we have taken a simplistic approach to the inclusion of NC-NC homotypic interactions. In the context of CIL, since neural crest cells express N-cadherins one should naturally consider that collisions between two migrating NCs will also generate mutual CIL responses (Carmona-Fontaine et al. 2008). Homotypic CIL interactions between two NC cells 1 can be trivially incorporated in a manner equivalent to the corresponding heterotypic ones: the first equation in Eq. (1) can be augmented with an additional velocity component, say v N N i, CIL , whose dynamics are governed by essentially the same law as for the contribution due to heterotypic NC-PC interactions (i.e., v N P i, CIL , defined in Eq. (7)), but now generated through NC-NC interactions. Simulations analogous to those of Sect. 3.3, see the representative case plotted in Fig. 19, though, highlight a subsequently poor agreement with in vitro outcomes: when homotypic CIL interactions are explicitly incorporated, excessive scattering of the NC colony occurs and "chase-and-run" is clearly impaired.
Experimental and theoretical studies (Camley et al. 2016a, b;Carmona-Fontaine et al. 2011;Merchant et al. 2018;Woods et al. 2014) support an additional NC-NC co-attraction mechanism, in which their secretion of an attracting substance (complement fragment C3a) counterbalances the effects of CIL-induced dispersal and maintains cohesion of the NC cluster. In the interests of model simplicity we have excluded explicit incorporation of this mechanism at present -explicit inclusion would demand an extra equation for the additional molecular species and an accompanying increase in the dimensionality of the parameter space. Rather, this co-attraction has been implicitly included, naively by assuming it counteracts the dispersal effects from homotypic CIL. Of course, an important development for future studies would be explicitly extending the model to account for this and determining in turn how it impacts on chase-and-run dynamics.
We further remark on the crescent-like shape of our chasing NC cluster in the reference simulation (see Fig. 12), which is somewhat different from the rounder aggregate observed in experimental controls, e.g. in Theveneau et al. (2013). Notably, NC crescents are observed under certain experiments scenarios, for example under partial inhibition of CIL (cf. second panel of Supplementary Movie 14 in Theveneau et al. 2013) or when clusters are initially separated by a greater distance (cf. Fig. 2 (h) in Theveneau et al. 2013), suggesting that crescent shapes are by no means unnatural features of the in vitro system. This phenomenon arises naturally from the model mechanism, where chemotaxis allows cells at the back of the NC aggregate to gradually work their way around the group and accumulate at the NC-PC interface. This crescent-like shape becomes more pronounced as the simulation proceeds (cf. Fig. 12 and its accompanying Supplementary Movie, M2, where earlier snapshots show a more circular geometry for the NC cluster), possibly suggesting that our simulations artificially accelerate this process. Notably, we have observed a tendency towards a reduction in the degree of crescent formation as certain parameters are shifted: for example, lower values of PC homotypic adhesion, reductions in the percentage of PCs that secrete Sdf1 and increments in the extension of the intercellular repulsive region. However, a crescent or ellipsoidal shape is not entirely eradicated, suggesting that further mechanisms may act to maintain the round shape of the NC cluster. Of course, a possible candidate lies in explicitly including the above mentioned coattraction mechanism, rather than the naive implicit inclusion here, or through some form of group polarisation/coordination within the NC cluster.
Both one-to-one and many-to-many simulations indicate that "chase-and-run" can generate an overall productive net cell migration, with NCs and PCs significantly shifting from their initial locations. However, the directionality is substantially different in the distinct settings. In one-to-one experiments random motility plays a fundamental role as the orientation of the PC, with respect to the NC, deviates significantly after the CIL response (and compare the axial direction connecting NC to PC in their starting and final locations, see Fig. 9g). In contrast, for cell clusters the Brownian crawling effects are diminished and net movement is predominantly along the axial direction connecting their initial centres of mass. In this respect, we can in principle claim that efficient movement along a specific direction demands cluster to cluster interactions rather than one cell to one cell interplay. Of course, directionality is strongly dependent on the initial placement of the cell aggregates and current simulations only cover distances of the order of 100 μm: in vivo, NC cells may need to travel significantly further and it is unclear whether the present mechanism would allow precise persistent movement without further guidance signals, such as long range chemical attractants or specific substrate patterns.
Parametrisation is a key issue for computational/mathematical models in the biosciences. Ideally, parameters could be obtained directly from biological data and, where possible, we have attempted to source such values (see Table 1). However, even if a parameter is estimated directly from its data, its reliability is questionable if its originates from distinct experimental studies, different developmental processes or from other species. Other parameters may have no clear biological analogue and here our approach has been to estimate such parameters through fitting to biological data (e.g. comparing characteristic run lengths). Building further details of molecular binding events and their subsequent impact on movement via inhibiting/promoting focal adhesion formation would potentially allow closer connection to quantitative data, although clearly such an approach would increase both the complexity of the model and the dimensionality of the parameter set. The principal aims of the current study are more exploratory in nature: employing a "top-down" modelling approach to test the fundamental hypothesis that a relatively minimal set of interactions can capture the biological observations of "chase-and-run". Nevertheless, further insights into the criticality and sensitivity of parameters is desirable and in this regard a more comprehensive sensitive analysis and further analytical investigation would be important.
For practical purposes we have targeted a specific cell system, where model hypotheses can be constrained. Nevertheless, the overall framework is generic and can be adapted to address other systems where similar mechanisms are believed to operate, e.g., in the case of cancer cell populations (Astin et al. 2010) or of developmental macrophages (Stramer et al. 2010). Here, contact-based interactions are mediated through cadherins and quasi-identical responses are assumed in the contacting cells. For other contact-based cell dynamics, such as under Eph/Ephrin signalling (Kania and Klein 2016) or in zebrafish pigment cells (Yamanaka and Kondo 2014), the reaction of each individual involved may be markedly different, involving attraction in the former case and repulsion in the latter. Thus, contact-based interactions can lead to a variety of homotypic (i.e. between cells of the same type) and heterotypic (between cells of different type) movement responses and a theoretical study into how various combinations act to control patterned collective movements of heterogeneous tissue populations would be of keen interest. In this respect, we refer to Painter et al. (2015) for a study using a fully continuous model.
While explicit analytical explorations are often difficult in agent-based models, we remark on two areas that could benefit from further mathematical exploration. First, we have built on recent investigations on the H-stability properties of adhesive/repulsive pairwise interaction potentials (Cañizo et al. 2015;Cañizo and Patacchini 2018;Carrillo et al. 2018;Ruelle 1969), that have allowed us to identify parameter regimes under which the model system evolves to a physically realistic configuration. Such analytical studies have in fact offered us welcome constraints for streamlining the tricky process of parametrisation. Second, we note that the chase time required for a neural crest cell to contact a placode cell (and subsequently initiate CIL) can be interpreted in the context of a mean first passage time problem (e.g., see McKenzie et al. 2009;Redner 2001): the average time needed by some object (in this case, the neural crest cell) to hit a target (the placode cell). A detailed exploration would require extension of existing theory, for example to account for the moving nature of the target and the fact that the bias is driven by chemotaxis. While we defer to a future investigation, such studies could shed light on, for example, how the chase phase is determined by the diffusive and reaction kinetics of the attractant.
The approach proposed in this work has extended previous modelling studies on CIL via an explicit study into how such a mechanism, coupled to other processes, drives the collective migration of heterogeneous interacting cell systems. Further, our computational framework has incorporated an explicit dynamical equation that governs the strength and the duration of the CIL response for each cell (see Eq. (7)). While our parametrisation has been guided here through an experimental-matching procedure, our model is open to future refinements: for example, connecting these parameters to an explicit representation of the intracellular signalling that regulates the CIL response. Such model developments would offer a potential pathway to targeted in silico predictions on the impact in cell behaviour of molecular-level perturbations.