Neural crest cells (NCCs) are produced at the neural plate during early stages of vertebrate embryogenesis. Then they migrate in different streams to various parts of the body that are taking shape, the cephalic stream to the head to form various facial features, and the enteric and cardiac streams toward the caudal end, and then from the dorsal to the ventral side to form the various organs (Fig. 1). This migration is rapid and precisely controlled by the various chemical signals at different locations, some attracting the NCCs while others repelling them. Thus, NCC migration is a prime example of chemotaxis [1, 2].

---

---

What is of interest to us in this article, however, is not the chemotaxis but the peculiar behavior of spontaneous collective migration, i.e. collective migration of a NCC cluster in the absence of external chemical cues. This has been demonstrated mostly by in vitro experiments. Mayor and coworkers [4, 5] have shown that an NCC cluster spontaneously migrates down a corridor coated with the adhesive molecule fibronectin, over the course of several hours. In vivo evidence also points to spontaneous collective migration at work [6, 7], although the evidence tends to be more circumstantial as it is difficult to ensure the absolute absence of an external chemoattractant. Spontaneous collective migration has been documented for other cell types in confined geometries, e.g. bovine capillary endothelial cells [8] and Madin-Darby canine kidney cells [9]. The question is: how do cells interact within a cluster to maintain persistent directional migration over long distances and long time durations?

Two critical factors in spontaneous collective migration have been identified. The first is contact inhibition of locomotion (CIL), the tendency of NCCs to move away from each other upon contact, a process mediated by N-cadherins and the non-canonical Wnt signalling pathway [10, 11]. Second, co-attraction (COA) describes how NCCs attract each other through the autocrine production of the short-ranged chemoattractant C3a and its receptor C3aR [4]. Computer models have successfully reproduced spontaneous migration of a group of NCCs under the simultaneous action of CIL and COA [12, 13]. However, we ask how CIL and COA, neither having any inherent directional preference, could produce persistent directional migration. What keeps the cell cluster moving down the corridor, instead of reversing course? We started out by building our own computer model that encodes CIL and COA, and it failed to recapitulate persistent collective migration. Once away from the end of the corridor, the group loses direction; its centroid wanders back and forth in essentially a 1-D Brownian motion.

From these observations we came up with the hypothesis that a symmetry-breaking mechanism has to be at work, in addition to CIL and COA. This we have called the persistent of polarity (POP). In recent work published in Developmental Biology [14], we have tested this hypothesis on the basis of the biochemical signaling pathways, as opposed to postulating rules governing cell interaction as commonly done in agent-based or rule-based modeling [12, 13]. This online article summarizes the main findings of [14]. First, CIL and COA arise naturally from the underlying reaction and diffusion of Rho GTPases. Second and more importantly, we demonstrate that CIL and COA, acting on a randomization scheme for cell polarity, produce persistence of cell polarity, and consequently spontaneous collective migration of NCCs. As it turns out, POP is not an additional rule to be posed alongside CIL and COA. It is in fact a natural outcome of CIL and COA, as well as a conduit through which these two fundamental mechanisms give rise to the observed spontaneous collective migration.

2 Methods

Our general conceptualization of the NCC collective migration is as follows. Polarization and protrusion of individual cells are governed by Rho GTPases on the membrane [15], subject to turnover between the membrane-bound and cytoplasmic forms of the signaling proteins. In particular, Rac1 and RhoA have been identified as the key proteins modulating CIL and COA [4]. Rac1 promotes F-actin assembly in lamellipodia at the protrusive front of the cell, while RhoA enhances myosin-induced cell contraction at the rear [16]. For simplicity, our model only accounts for Rac1 and RhoA and omits other GTPases such as Cdc42 and various downstream regulators. The level of active Rac1 and RhoA on the membrane determines the protrusive and contractile forces on the membrane and in turn the deformation and movement of the cell. Cell-cell interaction occurs through modulating each other’s Rac-Rho biochemistry. For example, cell-cell contact inhibits Rac1 and elevates RhoA at the site of contact in both cells. Thus, the cell protrusions retract and the cells move apart. Finally, the Rac-Rho dynamics is subject to a random noise so as to produce the tortuous trajectory commonly seen for single migrating cells.

---

---

Each cell is represented by a polygon of $N$ vertices connected by elastic edges (Fig. 2), similar to the vertex models widely used for epithelial morphogenesis [17, 18]. Rac1 and RhoA levels are defined on the vertices as well as in the cytoplasm. Through “mesh refinement”, we have chosen $N=16$ for all simulations in the rest of the paper [14]. Each vertex on the cell membrane is subject to 2D forces acting in the plane: a pressure force from the cytoplasm enclosed by the cell membrane, cortical tension in the membrane segments, and a protrusion or contraction force determined by the Rac1 and RhoA levels on the vertex. As a result, the vertices move, without inertia, at a speed determined by the resultant force and a friction factor. We do not explicitly account for focal adhesions between the cells and a substrate.

The biochemical and mechanical components of the model are intimately coupled. For the biochemical component, the model tracks three forms of the Rac1 and RhoA proteins: the active membrane-bound form, the inactive membrane-bound form, and the inactive cytosolic form. Note that the membrane-bound forms are defined on the cell boundary vertices, and may exhibit spatial distributions. In fact, cell polarization will be marked by spatially segregated distributions of active Rac1 and RhoA. The cytosolic levels are functions of time but not space; we assume the cytosol to be well mixed due to fast bulk diffusion [19, 20]. The total amounts of Rac1 and RhoA are each conserved.

On the membrane, the active and inactive forms interconvert with activation and deactivation rates, denoted by $K^{+}$ and $K^{-}$ for Rac1 and ${\kappa }^{+}$ and ${\kappa }^{-}$ for RhoA. Only the inactive form of protein may dissociate from the membrane to diffuse within the cytosol. The membrane association and dissociation rates are denoted by $M^{+}$ and $M^{-}$ for Rac1, and by ${\mu }^{+}$ and ${\mu }^{-}$ for RhoA. Using Fick’s law to compute the 1D diffusion flux $J$ of the active and inactive Rac on the membrane, we discretize the reaction-diffusion equations for Rac1 as follows:

$$\begin{array}{lll}\hfill & J_i^X=-D\left(\frac{R_{i+1}^X∕L_i-R_i^X∕L_i}{l_i}\right)& \hfill \text{(1)}\\ \hfill & \frac{\mathrm{d}{R}_i^a}{\mathrm{d}t}=K^{+}{R}_i^i-K^{-}{R}_i^a+\left(J_{i-1}^a-J_i^a\right),& \hfill \text{(2)}\\ \hfill & \frac{\mathrm{d}{R}_i^i}{\mathrm{d}t}=-K^{+}{R}_i^i+K^{-}{R}_i^a+\left(J_{i-1}^i-J_i^i\right)+\frac{M^{+}{R}^c}{N}-M^{-}{R}_i^i,& \hfill \text{(3)}\\ \hfill & \frac{\mathrm{d}{R}^c}{\mathrm{d}t}=\sum _{i=1}^N\left(-\frac{M^{+}{R}^c}{N}+M^{-}{R}_i^i\right),& \hfill \text{(4)}\end{array}$$

where the subscript $i$ indicates the $i^{\mathrm{\text{th}}}$ vertex on the cell membrane, $l_i$ is the edge length between vertex $i$ and $i+1$, $L_i$ is the average of $l_i$ and $l_{i-1}$, and $D$ is the diffusivity on the membrane. $J_i$ approximates the diffusive flux from vertex $i$ to vertex $i+1$, its superscript $X$ being $a$ or $i$ for the active and inactive forms of Rac. The amounts of Rac1 $R^a$, $R^i$ and the cytosolic $R^c$ are normalized by the total amount of Rac1 in the cell. Similar equations can be written for the normalized amounts of RhoA ${\rho }^a$, ${\rho }^i$ and ${\rho }^c$, but we omit them for brevity. These equations imply conservation of the total amount of Rac1 and RhoA.

The biochemical interactions are encoded in the activation and inactivation rates $K^{\pm }$ and ${\kappa }^{\pm }$ according to biological observations and prior modeling in the literature [14]. They give rise to cell polarization, stochastic repolarization (i.e. random changes in migration direction), CIL and COA:

• Polarization. To capture cell polarity, the activation rates of Rac1 and RhoA each reflect the species’ autocatalytic capacity, while their de-activation rates reflect mutual inhibition on the cell membrane [21]. Following [19, 21], we represent the Rac and Rho auto-activation and mutual inhibition through Hill functions, which allow spontaneous polarization of cells with Rac1 and RhoA peaking on opposite sides of a cell. This polarity is the precursor of cell motility.
• Stochastic repolarization. NCCs produce random Rac1-mediated protrusions [10] that compete with existing protrusive fronts and potentially change the cell’s polarization and direction of motion. To model this, we randomly select a subset of vertices and up-regulate the Rac1 activation rate on them every $20$ minutes. The model also accounts for competition between protrusions by a negative feedback of the cortical tension on Rac1 activation [22]. Thus, hotspots of Rac1 activity compete with each other on the membrane, allowing an upstart to replace an existing protrusion on occasion. This concept has been confirmed by recent biological observations on chick embryos [23].
• Contact inhibition of locomotion. Contact between two NCCs is known to activate the non-canonical Wnt signalling pathway, which results in the down-regulation of Rac1 and the up-regulation of RhoA [24], leading to CIL. This is effected in the kinetic equations by CIL factors that elevate the Rac1 deactivation rate $K^{-}$ and the RhoA activation rate ${\kappa }^{+}$ on any vertex that has come sufficiently close to a neighboring cell.
• Co-attraction. We realize COA by up-regulating the Rac1 activation rate on any vertex of a cell that is sufficiently close to a neighboring cell to “sense” the C3a distribution of the latter [4]. This enhances protrusion toward each other between neighbors and produces COA. Since C3a diffuses fast through the extracellular matrix, we assume a steady state exponential distribution of C3a surrounding each NCC  [12].

As in previous vertex models, the mechanical model represents the cell’s shape and movement via the position ${\overrightarrow{r}}_i$ of each vertex and its speed $\frac{\mathrm{d}\overrightarrow{r_i}}{\mathrm{d}t}$: $$\begin{array}{lll}\hfill \eta \frac{\mathrm{d}{\overrightarrow{r}}_i}{\mathrm{d}t}={\overrightarrow{E}}_{i-1}+{\overrightarrow{E}}_i+\left(p_i+F_i\right){\overrightarrow{n}}_i,& & \hfill \text{(5)}\end{array}$$

where ${\overrightarrow{E}}_i$ is the linearly elastic tension along the edge between vertices $i$ and $i+1$, $p_i$ is a cytoplasmic pressure that resists cell area reduction, $F_i$ is the protrusion or retraction force due to actin filaments or myosin motors on the membrane, and ${\overrightarrow{n}}_i$ is the unit outward normal vector at vertex $i$ (Fig. 2). Details of the mechanical model, along with the evaluation of parameters and a study of the model’s sensitivity to the parameters, can be found in the full article [14].

3 Results

The biochemo-mechanical model described above is able to reproduce single-cell behavior such as polarization, motility and stochastic repolarization. Movie A gives an overall impression. Besides, contact inhibition and co-attraction have been captured, as illustrated in Movies B and C, respectively. In this context, POP emerges naturally from the simultaneous action of CIL and COA, as illustrated by the two-cell migration simulation discussed below.

3.1 The origin of POP

Movie D depicts the persistent directional migration of two cells. At the start of the simulation, the cells are assigned random initial Rac1 and RhoA distributions, and placed next to the left end of the corridor, whose boundaries confine the cells via CIL in a similar way as cell-cell contact. Because of this confinement and the initial asymmetry in the geometry, the two cells develop protrusion fronts toward the right, and the cells start to move to the right ($t=47$ min). What is somewhat surprising is that this motion is sustained for the entire duration of the simulation (10 hours) across some $40$ cell diameters, the initial portion of the migration being shown in Fig. 3. This contrasts the behavior of a single cell that, after the initial movement away from the left end, essentially executes a 1D random walk in the corridor.

---

---

We note first that COA keeps the two cells within close proximity for the entire duration. This ensures the continual interaction between the two through CIL and COA. The effect of CIL on the leading cell is such that it never develops a viable new protrusion in the rear. Any such hotspot, as appears at 82 min, is quickly extinguished by CIL as the budding protrusion comes into contact with the trailing cell (116 min). If such contact weakens the protrusion front of the trailing cell, this effect is short-lived as a slow-down of the trailer will end the contact. The nascent Rac1 hotspot on the rear of the trailing cell, visible at 133 min, cannot compete with the forward protrusion that is reinforced by COA. As a result, the pair continues its directional migration, and similar cycles of interaction repeat in time (e.g. 200 min). Thus, COA and CIL act together to suppress the bursts of Rac1 up-regulation on the cell membrane, which would have produced random repolarization on a single cell. POP arises naturally from the pull-push action due to COA and CIL, and perpetuates the initial asymmetric motion of the cells endowed by the geometric confinement.

This insight answers the main questions that had motivated our model. If CIL and COA are postulated as ad hoc rules on the supracellular scale, they are insufficient for spontaneous collective migration. From this we have hypothesized that POP is a necessary third ingredient. Having built our model from the underlying GTPase dynamics on the intracellular scale, however, we find that POP emerges from the collaboration between COA and CIL, and need not be added separately on this scale. What has been added to our model, and missing from prior rule-based models [12, 13], is the Rac-Rho biochemistry that allows a description and rational explanation of random repolarization. The fact that Rac1 suppression promotes persistence in cell motility is well established experimentally for several cell types [25, 26, 27, 23], and indeed this was what motivated our modeling on the level of intracellular biochemistry.

It is important to note that POP is stochastic in nature. If, by chance, new Rac1 hotspots appear simultaneously on the rear end of both cells, they may overcome the forward polarity and cause the pair to reverse course. Then POP would fail. Whether this occurs depends on the initial Rac1 and RhoA distributions, and on the number of cells in the cluster. We will discuss the robustness of POP next.

3.2 Spontaneous collection migration of larger clusters

The key features of the persistent migration of two cells carry over to the collective migration of clusters comprising more cells. Figure 4 and Movie E demonstrate the spontaneous collective migration of 49 cells down a corridor. The cluster migrates persistently as a whole for about 20 cell diameters during a 10-hour period. Although the cluster area fluctuates, COA keeps the coherence of the group so that the mean area does not increase over time.

---

---

To quantify the cluster size effect on the spontaneous collective migration, we have to distinguish it from the confinement effect of the corridor. With a larger number of cells inside the same corridor, the cluster will be more crowded and this amounts to an effectively stronger confinement. To separate the cluster size and confinement effects, we vary the width of the corridor $w$, scaled by the cell diameter $d$, with the cell number $n$ as $w=\sqrt{n}$. Thus, we always have a square initial array of cells enclosed on three sides by walls, as in Fig. 4, $t=0$. This gives us, in a sense, “equal confinement” among the different $n$.

---

---

Figure 5 depicts the variation of several indicators of collective migration with the cluster size $n$. In panel (a), the averaged cell-cell separation $S$ is maintained roughly constant for different $n$ values by adjusting the strength of co-attraction. As co-attraction from different neighbors is additive, having a larger $n$ and thus more neighbors will elevate the COA effect. This is normalized by reducing the pairwise COA strength so as to produce the same $S$. Panels (b–e) shows the persistence ratio $R_p$, persistence time $T_p$, average centroid migration speed $V_c$ and the migration intensity $m_I$. Their definitions are given in the caption to Fig. 5. Taken together, these data show that the intensity of persistent collective migration increases with increasing $n$, rapidly for small clusters ($n=\mathrm{}$4  to 9) and mildly for the larger ones.

---

---

The explanation lies in the fact that for small clusters, POP has a greater probability of failure, when the cluster reverses course as a whole. For larger clusters, the chance of POP failure is much reduced. Figure 6 illustrates the stochastic nature of POP fallibility. Take $n=9$ for example. Among the 20 trajectories, the 10 that have traveled the farthest have never suffered a reversal, i.e. $V_c>0$ for the entire duration. For the rest, reversal may first set in as early as $t=\mathrm{}\mathrm{}$ 60 min or as late as 520 min. The failure of POP requires a large portion of the group to develop rear protrusions simultaneously. If only one or a few cells acquire rear hotspots among many, they will be quick extinguished by their neighbors through CIL. This explains why the probability of POP failure declines with increasing $n$. After all, as POP arises from cell-cell interactions, it should be more robust in larger clusters with more frequent and reliable interactions.

Finally, we can make some quantitative comparisons between our model prediction and experimental data. For NCCs migrating in vivo in Xenopus laevis embryos, Szabo et al. 13] reported a collective migration speed of $V_c\approx 1.5$ $\text{µm}$/min. In vitro, the NCCs migrate at $V_c\approx 1$ $\text{µm}$/min down a corridor coated with fibronectin and confined by borders that are rich in versican, an extracellular protein that repels NCCs. With the parameters used here [14], our model predicts an average speed $V_c=1.12$ $\text{µm}$/min for the larger clusters ($n\ge 9$) (Fig. 5d), comparable with the experimental data. Note that in our model, we have adjusted the contractile force factor in $F_i$ (Eq. 5) to produce the experimental speed of 3 $\text{µm}$/min for a single cell in the rectilinear phase of motility. The agreement in cluster migration speed $V_c$ is not the result of parameter fitting. For the persistence ratio $R_p$, Szabo et al. [13] documented $R_p\approx 0.85$ over a time period of 4 hours in vivo, and $R_p\approx 0.87$ in vitro. In Fig. 5, we have predicted an average persistence of $R_p=0.77$ for the larger clusters ($n\ge 9$) over 10 hours, comparable to the experimental values. A caveat about comparison with experiment is that the model has a large number of parameters, not all of which can be ascertained directly from experimental measurements.

In the original article [14], we have also explored how the confinement in the corridor affects the migration of the cell cluster, and found that optimal migration efficiency is achieved when the corridor width equals $\sqrt{n}$. Interested readers can consult the full paper.

4 Discussion

Our goal is to provide an explanation of spontaneous collective migration of neural crest cells in terms of key GTPases that regulate cell polarization and protrusion. From a broad perspective, we can summarize this work by the following three points.

First, we have adopted a modeling approach based on the biochemical signalling pathways as an alternative to the existing paradigms of continuum and agent-based modeling. Continuum models seek a “mean-field” description of the average orientation and movement of many cells without resolving the scales of individual cells [29]. However, it is challenging to represent intercellular interaction within the continuum formulation. Agent-based models address this concern by assigning rules on individual agents that recapitulate known cell-level behavior [8, 12, 13]. While helpful when little is known about the underlying biology, these models tend simply to give back the behavior designed for. The biochemically-based model uses a finer level of resolution. Moreover, it seeks to connect cell-level behavior to intracellular signaling pathways that have been established experimentally. Thus, it can be more general and even simpler than models based on postulated rules. For example, our model predicts cell behaviors ranging from polarization to spontaneous collective migration from a few well-establish principles governing GTPases [30] and cell mechanics [31]. We should note that in recent years, biochemically-based modeling of cell mechanics has gained increasing interest and currency  [32, 33]. This is evidently motivated and enabled by growing experimental knowledge of the role of chemical signaling in complex behaviours such as collective migration.

Second, we have identified the persistence of polarity, or POP, as an essential factor in enabling spontaneous collective migration of cell clusters. Prior to this work, the extensive experimental studies of Mayor and coworkers had identified CIL and COA as key to spontaneous collective migration. In their recent models, Woods et al. [12] treated the cells as particles that move with inertia, according to Newton’s second law of motion. Szabo et al. [13] used lattice-based rules that preferred the direction of motion of the preceding Monte Carlo step in a cellular Potts model. Although POP appears merely as an incidental feature of these models, we have come to suspect that it may in fact be essential for the appearance of spontaneous collective migration. To test this hypothesis, we built a simpler model that recapitulates CIL and COA but with POP intentionally suppressed by erasing the cell polarity periodically and replacing it by a random “initial distribution” of the Rho GTPases. This model does not exhibit spontaneous collective migration in a corridor; the centroid of the cell cluster executes a 1D random walk once the cluster has moved away from the end of the corridor.

Third, we have demonstrated that POP emerges naturally from the simultaneous action of COA and CIL. Thus, we have not only validated the existing proposal that COA and CIL cause spontaneous collective migration [10, 4], but also placed it on the concrete biological basis of Rac-Rho signaling. COA maintains the integrity of a cell cluster and ensures continual proximity and CIL interaction among neighbors. The two cooperate to suppress Rac1 and perpetuate an initial polarity and direction of migration that may have arisen from asymmetric boundary conditions and confinement. This is consistent with experimental observations that Rac1 suppression by drugs or Syndecan-4 enhances persistence in the polarity and motility in fibroblast and NCCs [25, 26, 27]. Thus, the model offers an explanation for the origin of POP; it is an emergent behavior based on COA and CIL rather than a separate mechanism to be postulated in addition to these two. Conceivably, this hypothesis can be tested by experiments that use drugs to partially suppress various proteins (or their activation levels) in the Rho GTPase signaling pathway.

It is interesting to put the current model in the context of directional motion and symmetry-breaking in vastly different systems, e.g. the active swimming of bacteria and flocking of insects or birds [34]. While alignment and directionality in macroscopic flocks rely on visual cues or pressure waves, we have shown that in the NCC context, it emerges from the pull-push actions of COA and CIL. Generally, one could view our POP as a device for alignment similar to the previously proposed directional persistence of random walk [8], persistence due to particle inertia [12], “persistence decay” for cell polarity [13] and collision-based alignment [35].

References

---

Supplemental Movies (mp4 format):

Movie A (283 KB): A single model cell showing random migration due to stochastic fluctuations in Rac1 activation rates. The lengths of line segments pointing in or out are proportional to the local levels of active GTPases (outwards blue lines for Rac1, inwards red lines for RhoA).

Movie B (1,142 KB): Two model cells undergoing contact inhibition of locomotion (CIL) upon contact in a 1-cell-diameter wide corridor. Co-attraction is turned off in this simulation to highlight the CIL repulsion.

Movie C (10,272 KB): Co-attraction (COA) maintains the integrity of a cell cluster, with nearest neighbors keeping a roughly constant average separation over time. Turning off COA between the cells would lead to dispersion of the group over time.

Movie D (552 KB): Two model cells exhibiting spontaneous collective migration in a 1-cell wide corridor.

Movie E (24,638 KB): A group of 49 cells exhibiting spontaneous collective migration in a 7-cell wide corridor.