Neural crest streaming as an emergent property of tissue interactions during morphogenesis

A fundamental question in embryo morphogenesis is how a complex pattern is established in seemingly uniform tissues. During vertebrate development, neural crest cells differentiate as a continuous mass of tissue along the neural tube and subsequently split into spatially distinct migratory streams to invade the rest of the embryo. How these streams are established is not well understood. Inhibitory signals surrounding the migratory streams led to the idea that position and size of streams are determined by a pre-pattern of such signals. While clear evidence for a pre-pattern in the cranial region is still lacking, all computational models of neural crest migration published so far have assumed a pre-pattern of negative signals that channel the neural crest into streams. Here we test the hypothesis that instead of following a pre-existing pattern, the cranial neural crest creates their own migratory pathway by interacting with the surrounding tissue. By combining theoretical modeling with experimentation, we show that streams emerge from the interaction of the hindbrain neural crest and the neighboring epibranchial placodal tissues, without the need for a pre-existing guidance cue. Our model suggests that the initial collective neural crest invasion is based on short-range repulsion and asymmetric attraction between neighboring tissues. The model provides a coherent explanation for the formation of cranial neural crest streams in concert with previously reported findings and our new in vivo observations. Our results point to a general mechanism of inducing collective invasion patterns.

Introduction Shape often plays an essential role for organ function. Therefore, understanding the process of shape acquisition, called morphogenesis, is crucial to understanding developmental processes and how to prevent their breakdown in pathologies. Studies over the last century identified a handful of universal modules controlling tissue morphogenesis, such as the spreading and thinning of epithelial sheets (epiboly) or convergent extension [1]. Most studies aim to understand morphogenesis without the need to consider environmental effects [2,3] despite the fact that developing tissues interact dynamically with their embryonic environment.
A striking example for the importance of environmental interactions during morphogenesis is the migration of the neural crest (NC). NC cells, an embryonic cell population whose migratory behavior has been likened to cancer invasion, are formed along the neural tube in the ectoderm and undergo an epithelial-to-mesenchymal transition (EMT) to form a single bulk pre-migratory NC tissue [4]. Subsequently, the NC cells are gathered into characteristic streams (Fig 1A and 1B) in which they migrate large distances collectively to contribute to a number of organs and give rise to a multitude of cell types [4]. Proper NC migration relies on environmental cues such as Semaphorin-3F (mouse) [5], or Versican (frogs) [6], Sdf1 (frog and fish) [7], Robo2 (chick) [8] (Fig 1A), as well as combinations of Eph-Ephrins (frogs) [9].
In the head, inhibitory cues such as Semaphorin-3F, Versican, and others are expressed in the epibranchial placodal cells surrounding the migratory streams originating from the hindbrain in the frog Xenopus laevis (Fig 1A, 1C-1E). In addition, EphA4, EphB1, and ephrin-B2 are expressed in the NC and matching mesoderm in Xenopus [9]. Blocking the function of these inhibitors impairs the correct formation of the NC streams in mice and frogs [5,6,9], suggesting that NC streams are shaped by a pre-existing environmental pattern of inhibitors that funnels the NC into correct migratory pathways as suggested for mice [10]. However, prior to migration the inhibitors are not patterned along the presumptive NC streams in the epibranchial region as they are expressed as a continuous band adjacent to the NC (Fig 1F-1I). Therefore, how the initial formation of the NC streams is established remains unknown.
Although several computational studies have been published on NC migration previously [6,[11][12][13][14][15][16], these mostly focused on stream directionality and cell-cell coordination during collective migration. Therefore, the question of how NC streams are generated has remained unexplored.
Here we investigate the idea that initial establishment of NC migratory streams results from the dynamic interactions of the NC and placodal tissues (Fig 1J). Placodal cells at this stage of development are motile and form a loose epithelial tissue in frogs and fish [7]. NC and placodes in these species are engaged in repulsion called contact inhibition of locomotion (CIL), whereby placodal and NC cells re-polarize and move away from one another after coming into contact as described for frogs and fish [17]. In addition, placodes attract adjacent NC cells through Sdf1 chemotaxis. The combination of these activities leads to a persistent directional migration of NC cells in a process named 'chase and run' which is essential for NC migration in frogs and fish [7]. We tested in silico and in vivo whether this chase and run interaction could lead to the initial emergence of NC streams, challenging the need for a pre-pattern to funnel the NC into streams.

Cell-based simulations of neural crest-Placode interactions predict emergent streams
To test whether cellular interactions can lead to formation of streams, we extended our previous model focusing solely of NC migration [6] by introducing placodal cells as detailed in the Materials and Methods section. Assumptions of the model are based on previously described behavior of Xenopus cranial NC. Briefly, two cell types (NC and placodes) are included in a cellular Potts model. Cells are self-propelled along their direction of polarization with higher (NC) or lower (placodes) motility persistence (Fig 2A) resulting in higher cell speeds for NC than placodal cells. Note that in vivo NC and placodal cells use an underlying fibronectin matrix as a migratory substrate. It is known that NC cells secrete a diffusible molecule, complement component C3a, which works as an attractant for the NC themselves in a process named 'co-attraction' (CoA) [11]. In the model CoA of NC is implemented through secretion of a NC chemoattractant (Fig 2B). It has been shown in frog and fish that when two NC cells come into contact they undergo a process termed contact inhibition of locomotion (CIL) [17,18]. In the model, the polarization vector of contacting cells is biased away from the contact, eventually leading to the separation of the two cells (Fig 2C). CIL and CoA together constitute repulsion upon contact and attraction at a distance between NC cells. Similarly, it has been demonstrated in frog that NC cells are attracted to Sdf1 secreted by the placodal cells [17,18], which is implemented in the model (Fig 2D). It is also known that NC cells and placodes in frogs exhibit a short-range repulsion [7]. This is implemented as CIL between these cells (Fig 2E). Additionally, based on the expression patterns of inhibitor molecules surrounding the NC streams ( Fig  1C and 1D), we assume that placodes secrete molecules that act to confine the migration of NC cells (Fig 2E). For simplicity, we model all inhibitory molecules as a single combined repulsive signal with an extremely short diffusion range (half-cell diameter) describing pericellular deposition of the signal. Attraction via Sdf1 and repulsion at a shorter range together result in a chase-and-run behavior between the NC and the placodes, whereby NC approach the placodes due to Sdf1 chemotaxis ('chase' phase), establish contact with the placodes, and are repelled from the contact site causing both NC and placodes to move away from each other ('run' phase) [7]. Placodal cells neither attract each other, nor do they exhibit CIL behavior [7]. Table 1 summarizes the cell-cell interactions in the NC-placode model.
At the start of simulations, placodal cells are initialized in 13 rows and NC cells in 2 rows in a rectangular area with closed boundary conditions ( Fig 2F). To mimic the progressive production of NC cells by EMT, the top region of the simulation area is monitored for empty space. Whenever a cell-free space sufficient to accommodate half a cell is detected, a new NC cell is introduced into the simulation at that space.
Simulations of the model show that NC cells invade the placodal region in distinct streams (Fig 2G and 2H, S1 Movie) without the need for any pre-pattern. While initially all NC cells facing the placodes experience the same micro-environment, parts of the interface later become specified as regions of invasion. Together with the NC streams, a pattern of the inhibitor also emerges in the simulations (Fig 2I, S2 Movie) matching the NC stream pattern similarly as observed in vivo (Fig 1C and 1D).  Streams in the model are characterized using a NC density probability function, similar to a pair correlation function of cell positions, defined as: denotes average over all possible R ! positions where Zð R ! Þ ¼ 1. Therefore, rð r ! Þ approximates the probability of finding two NC cells at position r ! relative to each other. In configurations with streams ( Fig 2G) this function outlines the local neighborhood of an average NC cell reminiscent of a stream ( Fig 2J). Sections along the axes of rð r ! Þ show average stream profiles (Fig 2K and 2L). A characteristic width, W = 36 ± 3 lattice sites (±SEM), and length, L = 79 ± 5 lattice sites (±SEM), are established as the distance of the two loci where the profiles first fall below 0.5, denoting 50% probability of NC cell presence. Time evolution of these measures and the aspect ratio L/W in the simulations shows that streams progressively reduce their width and increase in length up to a point after which the streams stabilize ( Fig  2M); note that the streams do not break into segregated clusters, even if the simulations are run for a longer time. This is accompanied by the saturation of the chemoattractant Sdf1, while inhibitor levels are still increasing (Fig 2N), leading to stabilization of the emergent geometry. This dynamic formation and evolution of streams in the model is similar to the migration of NC observed in vivo in frogs and fish [7,11].
In summary, our model predicts that NC streams may emerge from dynamic interactions between NC and placodal cells in the absence of an underlying guidance pattern.

Ectopic neural crest forms streams and affects surrounding inhibitory signal expression
We tested the prediction of our model that NC streams emerge independently from an underlying guiding pattern as suggested by the prevalent view; that is: a set of pre-existing lines that guide the NC into streams. To test this experimentally, we surgically removed the NC and its neighboring placodal population from pre-migratory stage Xenopus laevis embryos (stage 16-18 N&F [19]) and grafted them in the same embryo in either their normal position and orientation as control ( Fig 3A) or at an angle of 90˚( Fig 3B) within the pharyngeal region. This was done by opening the epidermis to retrieve the NC and placodes. The epidermis and the mesoderm have a relatively low adhesion to the NC and placode tissues, allowing to easily separate them. After re-positioning the NC and placodal tissues the epidermis is re-closed. Importantly, the mesoderm and overlying epidermis were neither transferred nor removed during the surgery. At these stages the prospective placodes form a continuous, approximately 5 cell thick band adjacent to the NC, called preplacodal region; this region lacks any obvious pattern in vertebrates [20]. The position of the NC was assessed by in situ hybridization (ISH) against the NC marker Twist at 0, 8, 12, and 16 hours after the graft healed. If the migration of NC into streams is dependent on a pre-pattern present in the embryo, the rotated NC should not form normal streams as they will not reach this hypothetical pre-pattering ( Fig 3B).
The control NC graft migrated to form streams in the normal orientation and position ( Fig  3C-3F), whereas the rotated NC grafts did not follow the path of normal NC migration, but still formed distinct streams that migrated according to their rotated orientation perpendicularly with respect to the control streams (Fig 3H-3M). These results show that NC stream formation can occur irrespective of the orientation of the under-or overlying tissues and support our model whereby the initial formation of NC streams is not governed by a pre-pattern.
In the following we test the proposed model by exploring its behavior under perturbed conditions in order to support our main results presented above.

Number of streams is determined by migratory region width
To further explore how the closed boundary condition along the AP axis in our model affects stream formation, a series of simulations were performed using a range of system widths. In these simulations the initial number of cell rows was kept constant therefore keeping the same linear cell densities. We observed stream formation irrespective of the system width ( Fig 4A-4C) and found a direct relationship between the number of streams and system width ( Fig 4D).  The in silico results were compared with in vivo behavior by performing a series of graft experiments in which the cranial region was extended by grafting segments of NC and placodes, containing the equivalent region from another embryo ( Fig 4E). Similar to our in silico results and control embryos, the extended segments were able to give rise to NC streams ( Fig  4F-4H). We found a similar relationship between the extension ratio and the number of streams in vivo as predicted by the simulations (Fig 4D).
Together these results indicate that the number of migratory streams is determined by the space available for the NC along the AP axis. Note that these results do not test the lack of prepattern and in themselves could be compatible with a pre-pattern driven stream formation. In vivo, these boundaries are most likely the forming eye from the anterior and Semaphorin-3A at the posterior boundary of the epibranchial region in frogs [21].

Long-range attraction and short-range repulsion are required for stream emergence
The effect of Sdf1 chemotaxis was explored in the model by abolishing the chemotactic effect of Sdf1 on NC cells. This led to the complete absence of streams compared to control simulations (Fig 5A and 5B) indicating that in silico stream formation highly depends on Sdf1 attraction. To confirm the in silico results, the effect of Sdf1 on NC cells was abolished in vivo by inhibiting the expression of the Sdf1 receptor CxcR4 using validated morpholino oligomers (Cxcr4-Mo) (S1 Fig) [7,22]. Embryos injected with the Mo had a severe reduction of NC migration compared with control-Mo injected embryos (Fig 5C and 5D), consistent with our in silico observations.
Simulations with different chemotaxis strength parameters show that stream width exhibits a minimum, while lengths reach a plateau (Fig 5E). At higher chemotaxis strengths increase in length is limited by the compaction of placodal cells and the limitation of the model to represent out-of-plane intercalation of the placodes. We note that at low chemotaxis parameter values this result also demonstrates that EMT in the model does not produce an artificial pressure on the NC population by adding NC cells at the dorsal side (top) of the simulation area which could otherwise force NC invasion into the placode region.
To test the requirement of short-range repulsion for stream formation, the CIL behavior was removed specifically from the placodal cells in simulations. This led to a uniform invasion of the placodal region; however, it severely compromised segmentation of the NC into distinct streams compared to control simulations (Fig 5F and 5A). This result was tested experimentally by grafting whereby wild-type host embryos received graft placodes from embryos in which CIL was inhibited by injection of a dominant negative Dsh construct (DshDEP+), that specifically inhibits in frogs the PCP Wnt signaling, required for CIL [7,17]. We found that treated embryos failed to develop well-formed NC streams compared to controls where the grafts were implanted from wild-type embryos (Fig 5G). Nevertheless, NC cells still invaded the DshDEP+ placodal region, consistently with our simulation predictions.
Next, we adjusted the probability of a CIL response in placodes after collision in simulations and measured the emergent morphologies. The extent of invasion remained mostly unaffected with increasing probability of CIL reaction as signified by the length of streams (Fig 5H). Stream width, on the other hand, increases significantly as placode CIL probability decreases and reaches the maximal width at around 30% CIL probability ( Fig 5H).
Finally, we investigated the interaction of the chemotaxis and contact-repulsion mechanisms in simulations. We mapped width and length of NC streams as rectangles in a 2D space of chemotaxis and CIL response parameters (Fig 5I), where shading represents variability of the measurements (light: high uncertainty, dark: low uncertainty). The simulation results reveal that chemotaxis and CIL counteract one another, yet they are both required for stream formation.
Together these results show that the chase-and-run cellular interaction between NC and placodes is essential for the formation of streams with the chase (Sdf1 chemotaxis) being responsible for powering the invasion and the run (NC-placode repulsion) being responsible for shaping the streams.

Emergent streams are shaped by cell-cell adhesion
To investigate the effect of cell-to-cell adhesion within the NC population on the emergence of streams, simulations were run with varying degrees of attachment between NC cells. Simulations with low adhesion (λ M (NC−NC) = 0) gave rise to thinner streams compared with control simulations (Fig 6A and 6B). As N-cadherin is the major cell adhesion molecule present in Xenopus cranial NC [22][23][24][25], we manipulated the adhesion in these cells by interfering with the levels of N-cadherin at the cell surface. It has been shown that in frogs the level of N-cadherin at the cell contact is controlled by N-cadherin endocytosis, which is regulated by lysophosphatidic acid receptor 2 (LPAR2) [26]. We have previously demonstrated that by changing the levels of LPAR2, adhesion strength of NC cells can be altered in frogs [26]. We decreased NC cell-cell adhesion by injecting LPAR2 mRNA, which leads to LPAR2 overexpression and gives rise to significantly thinner NC streams when compared with NC injected with a control Mo (Fig 6C, 6D and 6G) or with wild-type. This is consistent with our simulations where NC cell adhesion was weakened.
Next, we simulated an increase in cell adhesion between NC cells (λ M (NC−NC) = 2), which produced wide streams (Fig 6E). We increased NC cell-cell adhesion in vivo by expressing a validated LPAR2 antisense Mo (S2 Fig), which in frogs leads to strengthening of NC cell adhesions (S3 Fig) [26]. Streams in LPAR2-Mo injected embryos were significantly thicker compared with embryos expressing control Mo (Fig 6C, 6F and 6G), consistently with the simulations.
Simulations with intermediate adhesion values revealed a strong and non-linear increase of stream width with increasing adhesion (Fig 6H). Stream length, as measured from the NC density probability function rð r ! Þ, remained relatively unaffected by changes in the adhesion apart from the extreme low parameter regime, where smaller length measures result from the thinness of streams (Fig 6H).
To investigate the interplay between placodal CIL, chemotaxis, and NC-NC adhesion during stream formation, we mapped stream geometries as the function of parameters related to these behaviors (Fig 6I and 6J). Our results reveal that adhesion and placodal CIL counteracts one another: stronger adhesion requires stronger CIL for similar streams. The interaction between chemotaxis and adhesion is more complex (Fig 6J). Both strong adhesion and strong chemotaxis leads to bulk invasion, while even a strong chemotaxis can produce thinner streams when adhesion is weakened.
Taken together, our simulations and experiments show that cell-cell adhesion plays an important role in NC stream morphology.

Discussion
Our experimental results show that migratory streams of the Xenopus cranial NC can form independently of their orientation with respect to the underlying mesoderm or endoderm. This questions the prevailing view that cranial NC follow a pre-pattern of environmental cues that is present prior to migration. Although a pre-pattern has been demonstrated in the path of trunk NC of mice and chick, for example by semaphorin-3F and versican [27], such demonstration is lacking in the epibranchial region prior to NC migration. At later stages of migration the cranial NC streams enter the pharyngeal arches, a series of composite structures segmenting the head [28,29]. The pattern of the arches is established through Bmp/FGF produced by the endoderm independently of the NC as shown in chick [30], which has been suggested to contribute to the NC stream pattern [31]. While the pattern of the arches explains the spacing of cranial NC streams at late migratory stages and refines the position of the NC according to the pharyngeal arches, it is unclear how early NC streams in the ectoderm would be guided by the distant endoderm. Indeed, in our ectopic graft experiments stream formation was apparently not hampered by neither the epidermis, nor the endoderm. Furthermore, cranial streams of the NC have also been observed prior to pharyngeal arch morphogenesis in fish [32,33]. The fact that NC streams and the pharyngeal arches seem to be initiated independently of one another [30] opens the question of how are they coordinated at later stages to form integrated structures.
Our findings are consistent with a model where invading NC streams emerge solely from a small set of known cellular interactions of NC and placodal cells (Fig 7). Rather than being tissue-autonomous, this model posits that patterning the NC into migratory streams depends on the presence of epibranchial placodes. This is consistent with previous findings in fish and frogs where NC streams failed to form when formation of placodes was inhibited [7] and may Emergence of neural crest streams provide insight into why segregation of NC into distinct streams occurred together with the appearance of epibranchial placodes during evolution [34]. While our results show that this model is sufficient for forming streams, we cannot exclude the presence of a molecular pattern (inhibitory or attractive) prior to cranial NC migration. However, the effect of such a pattern on migration would be insufficient to explain our results (Fig 3).
This model can account for the observed inhibitory patterns surrounding the cranial streams (for example Fig 1C and 1D) that led to the pre-pattern hypothesis. In our model these guidance cues are expressed by the placodes and, as the placodes are compacted by the invading NC, the pattern of inhibitors emerges at the interface of the two tissues. Consistently with previous findings in chick and frogs [5,6,35], these inhibitors are required in our model for stream formation to provide the necessary repulsion between NC and placodes, but are not required to be patterned prior to migration. This is also in line with recent findings in lamprey showing that Sema3F does not guide NC segmentation directly [36]. Furthermore, in chick it has been show that the migrating NC not only require the presence of an inhibitor, but it also re-shapes the pattern of the inhibitor [37].
Our results agree with a previous work focusing on the difference between trunk and head NC migration in chicken embryos [13] which showed that the NC has a tendency to form streams. This earlier study presented a model in which invasion is driven by a global directional bias in the movement of NC cells which is hampered by a uniform extracellular matrix (ECM) distribution. When cells are allowed to digest the ECM at a faster rate, invasion in that model proceeds in a wider front [13]. This compares well with our simulations where repulsion between NC and placodes is rendered less effective, giving rise to wider NC streams (Figs 5F, 1H and 1I). Formulation of our model allows testing a wider range of parameters, leading to unstructured or even inhibited invasion. While simulations of Wynn and coworkers [13] suggest that higher NC adhesion leads to thinner invasion, our simulated and experimental results reveal an opposite trend whereby more adhesive NC tends to form wider streams ( Fig  6). The different simulation results might stem from differences in model implementations and suggest experimental measurements and perturbations of cell-cell adhesions in the chick system. In support of our model, increasing cell-cell adhesion in Xenopus NC cells was sufficient to reduce their ability to invade narrow tunnels [26] indicating that the ability to form thin streams is inversely correlated with the strength of cell-cell adhesion. Forming thin finger-like structures from an initially wide mass of cells requires intercalation, a process that is impaired by strong cell-cell adhesion in Xenopus [38][39][40]. While our study focused on the cranial NC of Xenopus, investigation of the model parameters (Fig 6) suggest that this model may be suitable to describe similar pattern formation of more mesenchymal cell populations in which cell-cell adhesion is even lower. In this case, our model would predict that a higher inter-tissue attraction (chemotaxis) would be required and a lower inter-tissue repulsion (CIL) would be sufficient.
Distinction between leader and follower cells in NC streams has been reported and was explained by a difference in the gene expression levels based on single-cell sequencing in chick [15]. Although leader cells were fated to become leaders in trunk regions in fish even before migration [41], cells would be required to switch between leader and follower fates within minutes according to another study done in chick [16]. We observed individual cells at the tip of streams in our simulations which invaded the placodal population on their own. Although these cells have the exact same set of instructions (or 'genotype') in the model, their behavior is different from the follower cells. These leaders experience a slightly different environment than the follower cells with higher proportion of placodal neighbors, leading to a noted difference in their behavior. A similar emergent mechanism has been proposed for the selection of leading enteric NC cells in mice during the colonization of the gut where cells that get at the forefront of invasion by chance are later differentiated as leaders due to this initial advantage [42]. Our model therefore provides an alternative explanation for the emergence of persistent leader/follower positions for cranial NC migration.
Collective cell movements forming invasive finger-like structures similar to the migratory NC streams have been studied at the free edge of expanding epithelia in vitro [43,44] and in theory [45][46][47]. This phenomenon is well described by a model where invasion is enhanced at regions of high curvature [48,49]. This assumption is supported by observations in other biological systems as well [50][51][52] and is related to the Mullins-Sekerka instability observed typically at the solidification interface of forming crystals [53]. This instability gives a mathematical description of finger-like growth driven by the gradient of a diffusible substance. In our model, invasion of NC into the placodal region is driven by Sdf1 chemotaxis, proportional to the gradient of Sdf1, and therefore it is possible that streams are formed through such a surface instability. However, we did not observe higher chemoattractant gradients at the tips of NC streams, as would be required by the physical instability, and therefore it is unlikely that this mechanism drives stream formation in our simulations. Our system is more analogous to the pressure driven Saffman-Taylor instability [54,55], also known as viscous fingering, at the interface of two non-mixing fluids (for example honey and water). When the less viscous fluid (water) is pressed into the more viscous one (honey), the invading fluid forms finger-like protrusions with a characteristic width (w). In the simplest case of Newtonian fluids, w is determined by the surface tension between the fluids (γ), the velocity of invasion (V), the difference in fluid viscosities (μ), and the height of confinement in 2D (b) as w � b ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi g=ðmVÞ p . Measurements of tissue viscosity and surface tension in various tissues of Xenopus embryos estimate that γ/μ~1.8μm/min [56]. Considering that NC streams in Xenopus are relatively flat, b~30 −50μm [57], and V~1.6μm/min [6], the Saffman-Taylor instability would predict stream widths in the order of w~50μm. Although this underestimates the measured 120μm stream width (Fig 6G), a better estimate would require measurements of γ and the difference in μ specifically for the NC and placode tissues at these stages of development. Analogy between this physical phenomenon and the appearance of NC stream pattern has previously been suggested [58] and is consistent with recent observations of changes in NC tissue tension during EMT [59]. However, a thorough investigation of the analogy is lacking and the underlying mechanisms may well be different in these seemingly similar phenomena. Formation of NC streams must be robust to accommodate unavoidable variations in embryos without affecting NC migration, therefore the NC system has to be stabilized by other mechanisms. This is partially achieved in our model by the accumulation of inhibitors (Fig 2N). However, our model focuses on the early stages of migration, on the establishment of initial streams. It is very likely that this model works together with other mechanisms, especially at later stages of migration, to refine the initial stream pattern and ensure each stream reliably reproduces its stereotypic shape in the precise location. One such important mechanism might be an interaction with the confining pharyngeal arches to achieve perfect alignment with the NC. Our model provides a numerical description of the NC system based on well-established molecular interactions between cells at a tissue interface. How parameters of our model relate to the physical parameters of the viscous finger instability is an intriguing question to be investigated in the future.
We believe that our results are applicable to all systems where two cell populations exhibit a chase-and-run behavior and are not otherwise constrained by their environment. In neural crest systems such conditions have been described in the cranial regions of Xenopus and zebrafish species [17]. Note that in the trunk neural crest regions the neural crest are restricted to migrate in pre-determined pathways [27,60] that are potentially present before the migration of the neural crest. Therefore, our model would not apply in this situation. Whether neural crest of other species exhibit a chase-and-run-behavior remains to be investigated. A mechanism similar to the chase-and-run has been proposed to drive the patterning of initially mixed populations of motile agents [61]. While such theoretical approaches allow a more complete understanding of the model behavior, simplicity often confounds their applicability and experimental testability.
In summary, we uncovered a general mechanism of collective invasion based on the stream formation of neural crest during early development. Our model generally predicts such patterns at the interface of tissues that exhibit long-range attraction and short-range repulsion. These cellular interactions may be present at the tumor-stroma interface [62] where they could enhance metastasis by facilitating collective invasion and may play a role in other developmental processes.

Embryology
Embryos were obtained as previously described [6]. Briefly, Xenopus laevis females were injected with hCG (500U) for stimulating ovulation. Eggs were fertilized in vitro using sperm macerated in 1×MMR (100 mM NaCl, 2mM KCl, 1mM MgSO 4 , 2mM CaCl 2 , 5 mM Hepes, 100 μM disodium-EDTA pH 7.6). Embryos were dejellied using 2% L-cysteine solution. When injected or for graft experiments, embryos were left to recover in normal amphibian medium (NAM) 3/8 and then transferred and maintained in NAM 1/10 [63]. NC and placode grafts were performed as previously described [7]; briefly, at mid-neurula stages (16-18 N&F) the epidermis covering the neural folds was peeled off and a band of deep ectoderm cells located in the lateral aspect of the neural fold containing prospective NC and the preplacodal region was dissected and grafted into host embryos. In Xenopus embryos NC and placodes are originated from the deep layer of the ectoderm and can be easily recognized by their location within the embryo and by the aspect of their cells, as it is remarkably different to all other surrounding tissues (neural plate, mesoderm, epidermis). Host NC and placode were removed before grafting and the grafted tissue was maintained in place with coverslips for several hours.

Antisense morpholino injections
Although Cxcr4 and LPAR2 Mo were characterized in previous publications [22,26] we proceed to analyze the specificity of these treatments in the context of the phenotypes described here. Inhibition of NC migration by Cxcr4 or LPAR2 Mos was efficiently rescued by co-injection of Cxcr4 or LPAR2 mRNA, respectively (S1 and S2 Figs). The mRNA used for the rescue do not bind the Mo sequence [22,26]. In addition, although we had previously shown that LPAR2 Mo affect cell-cell adhesion, we proceed to perform an additional control here. As Ncadherin is one of the major adhesion molecules expressed by Xenopus migrating neural crest cells [12,69] we showed (S3 Fig) that the effect of LPAR2 Mo is completely rescued by expressing a dominant negative of N-cadherin [n-cadΔ] [26]. This result is consistent with the conclusion that LPAR2 Mo leads to an increase in N-cadherin at the cell junction [26] and therefore to higher cell-cell adhesion.

Imaging
Images were acquired using a stereomicroscope (MZ FLI II; Leica Biosystems) fitted with a Plan 1.0×/0.125 objective and a camera (DFC420; Leica Biosystems). Data were acquired using IM50 v5 software (Leica Biosystems).

Statistical analysis
Experiments were repeated at least three times independently on different weeks using different females and males for fertilization to generate embryos. Simulations were repeated 20 times using different random seeds. Data are represented as mean, spread is represented by SEM. For testing significance, the non-parametric Mann-Whitney u-test was used.

Quantifications
Stream width in experiments was measured using ImageJ line tool on calibrated ISH images. Widths for each stream were measured at five different positions and then averaged. Number of streams was assessed manually.
Chemotaxis, persistent cell adhesions, and polarized cell movement is described through function W as: Parameters λ s , λ M , and λ P set the relative weight of chemotaxis with respect to substance s, the persistent cell adhesions, and polarized cell movement. The first term describes extensiononly chemotaxis [51] where chemotaxis only has an effect on free edge extensions of cells, not at contact sites or during retraction; summation here is over diffusing molecules s (CoA, Sdf1, and inhibitor). r ði;jÞ �! in the second term is the vectorial distance of cells i and j, d 0 is the equilibrium distance of persistently adhered cells (d 0 ¼ 2 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi V T =p p ), and n(i) are the persistently adhered neighbors of cell i. As described previously [6,72] cells in contact are connected with a 10% probability in each MCS and are disconnected once cells lose contact and otherwise with a probability of ðjr i;j ! j À d 0 Þ=100. In the final term in Eq 2, Dr i �! is the displacement of cell i in the last MCS and p i ! is the polarization vector of cell i.
Each cell in the model possesses a polarization vector p i ! which is updated in every MCS as The value of the decay parameter δ P (i) sets the motion persistence of cells and depends on the cell type (NC or placode) and whether the cell is in contact with other cells or not as in [6].
CIL behavior is modelled as follows. When two cells undergoing CIL get in contact, each of the cells can independently switch to CIL state with a probability p(CIL). When a cell is in CIL state, its polarization vector is updated as p i ðt þ 1Þ

Model parameters and initial conditions
We model three diffusing substances: A (for C3a), S (for Sdf1), and I (for inhibitor) (Fig 2B,  2D and 2E). A is secreted by NC [11], while S and I are secreted by placodes [7]; in addition, S is taken up by NC cells. The parameters used for the main simulations are listed in Table 2. Cells were initialized as 5 × 5 patches on a 180 × 95 lattice with 2 rows of NC cells and 13 rows of placodal cells ( Fig 2F) and all concentrations were set to zero. As described in the main text, cell-free space is monitored at the top edge of the simulation area and whenever a region of a cell width and at least half a cell height is detected, a new NC cell is initiated there. Simulations are run for 3600 MCS. To allow comparison with experimental observations, we relate a distance of one lattice site to 3.5μm and one Monte Carlo time step to 10 seconds real time. With these relations, the simulated cell diameter is approximately 17μm, the average stream width emerging in the main simulations is 120μm, the total simulation time (3600MCS) is 10h, and the average simulated cell speed (0.1 pixels / MCS) converts to 2.1μm/min. These values are in agreement with the experimentally observed ones: the diameter of a cell spread on a surface is in the range of 10 − 25μm; the average stream width is 120 μm (Fig 6H); the approximate time required for NC stream formation and migration is 6-10h; the reported NC cell speeds are in the range of 0.5-2 μm/min [6].   Validation: András Szabó.