In this case study, we discuss mathematical modelling in the context of collective cell migration in the chick cranial neural crest. This work formed part of my doctorate, which was supervised by Professor Ruth Baker and Professor Philip Maini. The discussion that follows is an abridged and modified version of our forthcoming book chapter, Integrating mathematical modelling and experiments in the study of collective neural crest cell migration, which will appear in the edited volume Patterns and Forms in Nature: From Organisms to Collective Behaviours in Biological Systems (Springer Singapore).
Cranial neural crest cells on the move
The cranial neural crest is a highly migratory population of cells that contributes to the formation of facial structure during embryonic development [1]. During development, cranial neural crest cells (CNCCs) delaminate from the hindbrain (rhombomeres r1–7; Figure 1(a)) and migrate collectively in streams towards targets in the branchial arches, of which, in the chick, there are five (intriguingly for a mathematician, these are labelled branchial arches (ba) 1, 2, 3, 4, and 6; ba1–4 are shown in Figure 1(a)). Once CNCCs reach the branchial arches, they begin to specialise into cells that underpin the formation of the craniofacial structure. If the migration of these streams is disrupted, then a range of birth defects, most notably cleft palate, can result [2]. As such, it is highly desirable to better understand the biological mechanisms that underpin this process, so that it may become possible to develop novel therapeutics that can either prevent or counteract the effects of aberrant collective cell migration.
A synergy between mathematical modelling and experiments
It is, relatively speaking, easy to formulate biological hypotheses. To test them, however, is another matter entirely. Biological experiments are costly in terms of both time and money. This means that testing a biological hypothesis is a highly non-trivial endeavour that can require months (potentially years, even decades [5]) of work from a team of experimentalists.
Mathematical modelling, therefore, offers a powerful complement to experiments, as it can be used to rapidly test and refine biological hypotheses, and thus to drastically narrow the search space of experiments in a time- and cost-effective manner. This so-called model-test-refine paradigm has underpinned the study of chick cranial neural crest migration for over two decades, most notably in a collaboration between mathematicians in Oxford, led by Professor Ruth Baker and Professor Philip Maini, and the experimental group of Professor Paul Kulesa (1962–2025), who spent the majority of this collaboration at the Stowers Institute for Medical Research in Kansas City, Missouri, USA.
The collaboration between Oxford and the Stowers Institute began with an attempt to uncover the cellular behaviours underpinning the collective migration of CNCCs from r4 to ba2 in the chick (Figure 1(a)). To this end, a then graduate student, Dr Louise Dyson (now Reader in Epidemiology at the University of Warwick), formulated an agent-based mathematical model of collective CNCC migration [3]. The model was hybrid in that it represented CNCCs as individually resolved agents that interacted with one another and their environment, and vascular endothelial growth factor (VEGF), a chemical towards which CNCCs migrate, as a continuously varying scalar field that evolved according to a non-local partial differential equation (PDE):
\begin{equation} \label{chemicalEquation} \frac{\partial c}{\partial t} + \underset{\text{(I)}}{ \vphantom{ \displaystyle c\sum\nolimits_{i=1}^{N(t)} \frac{\lambda}{2\pi R^{2}} \exp\left[ -\frac{(x-x_i)^{2}+(y-y_i)^{2}}{2R^{2}} \right] } \frac{\partial(ac)}{\partial x} } = \underset{\text{(II)}}{ \vphantom{ \displaystyle c\sum\nolimits_{i=1}^{N(t)} \frac{\lambda}{2\pi R^{2}} \exp\left[ -\frac{(x-x_i)^{2}+(y-y_i)^{2}}{2R^{2}} \right] } D\left( \frac{\partial^{2}c}{\partial x^{2}} + \frac{\partial^{2}c}{\partial y^{2}} \right) } - \underset{\text{(III)}}{ \vphantom{ \displaystyle c\sum\nolimits_{i=1}^{N(t)} \frac{\lambda}{2\pi R^{2}} \exp\left[ -\frac{(x-x_i)^{2}+(y-y_i)^{2}}{2R^{2}} \right] } c\sum\nolimits_{i=1}^{N(t)} \frac{\lambda}{2\pi R^{2}} \exp\left[ -\frac{(x-x_i)^{2}+(y-y_i)^{2}}{2R^{2}} \right] } + \underset{\text{(IV)}}{ \vphantom{ \displaystyle c\sum\nolimits_{i=1}^{N(t)} \frac{\lambda}{2\pi R^{2}} \exp\left[ -\frac{(x-x_i)^{2}+(y-y_i)^{2}}{2R^{2}} \right] } \kappa c(1-c) }. \end{equation} where $c(x, y, t)$ is the VEGF concentration at a position $(x, y)$ and time $t$. Term ($\text{I}$) of this non-local PDE accounts for domain growth with a local tissue expansion velocity $a$, and term ($\text{II}$) represents the diffusion of VEGF with a diffusion coefficient $D$. Term ($\text{III}$) represents cell-induced degradation. Each of the $N(t)$ cells present at time $t$ degrades VEGF at a rate determined by $\lambda$ over a Gaussian neighbourhood of width $R$ centred on its position $(x_i, y_i)$. Term ($\text{IV}$) represents the logistic production of VEGF in the domain at a rate $\kappa$.
On the basis of preliminary simulations, Dyson, Baker, and Maini's model posited the existence of two cell phenotypes (two sets of cellular characteristics and associated behaviours) within CNCC streams. CNCCs at the front of streams were assigned leader phenotypes and migrated according to local concentration gradients in VEGF (which are generated by term ($\text{III}$) in the non-local PDE for VEGF concentration). The mechanism of VEGF sensing (in most iterations of this model, but not the first) incorporated a concentration-dependent detection threshold based on the seminal theoretical work regarding the physics of chemoreception by Berg and Purcell [6]. At the point of movement, leader cells extended a fixed number of filopodia (for a mathematician, antennae formed at their membrane). A filopodial protrusion then sensed a VEGF concentration gradient when
\begin{equation}
\label{sensing}
\frac{c_f(x,y,t)-c_b(x,y,t)}
{c_b(x,y,t)}
\geq
\xi
\left(
\frac{c_0}{c_b(x,y,t)}
\right)^{1/2}.
\end{equation}
where $c_f(x, y, t)$ is the average VEGF concentration along the filopodium, $c_b(x, y, t)$ is the VEGF concentration at the position of a cell's current centre, $c_0$ is the initial uniform concentration of VEGF within the domain, and $\xi>0$ is a dimensionless parameter that determines the sensitivity of chemotactic gradient detection by cells. If a gradient in VEGF was detected by one or more filopodia, then movement occurred in the direction of the highest detected gradient and was otherwise unbiased. All other cells in the model were assigned follower phenotypes, and migrated according to cell–cell guidance cues, forming long connected chains of cells led by leaders (see Figure 1(b) and the video below, in which leader cells are green and follower cells are red).
Motivated by the consistency between simulations of this model and observations of CNCC migration in the chick, Professor Kulesa's experimental group then sought to validate the leader–follower hypothesis by profiling gene expression levels within CNCC streams. They divided streams migrating from r4 towards ba2 into leading and trailing subpopulations according to their distance from r4 (Figure 1(a)), and, using both laser-capture microdissection and fluorescence-activated cell sorting, screened 84 genes implicated in CNCC migration across the two subpopulations. Nineteen of these genes were found to be differentially expressed between leading and trailing CNCC subpopulations [3]. In these experiments, CNCCs in leading positions were found to express genes associated with movement along VEGF gradients at levels higher than cells in trailing populations. Conversely, trailing cells were found to express genes associated with cell–cell adhesion that were distinct from those expressed by leading subpopulations.
These experiments thus demonstrated concurrence with the leader–follower model hypothesis in the chick by showing that the transcriptome of chick CNCCs (the genes they express and the levels at which they express them) varied within streams in a manner correlated with their relative stream positions. The leader–follower model has since formed the basis of many modelling–experiment collaborations (see [7] for a review of work conducted before 2021), and has been used to uncover many additional mechanisms underlying collective CNCC migration in the chick.
But wait...
Despite the success of the leader–follower model paradigm and the great number of insights it has yielded into chick CNCC migration, it is based on a simplification that categorises cells into only two discrete types. A study by Dr Linus Schumacher (Reader in Computational Tissue Dynamics at the University of Edinburgh and another former graduate student of Professor Baker and Professor Maini) re-examined transcriptome-wide single-cell RNA-sequencing data for CNCCs (rather than the 84 genes used to test the leader–follower hypothesis) and found that the underlying transcriptomic variation within CNCC streams is likely continuous, rather than discrete (Figure 2(a)). These findings suggested that CNCCs may occupy a continuous set of states that smoothly interpolate between leader-like and follower-like extremes [8,9]. These observations were, in turn, used to motivate a continuous-phenotype agent-based model of CNCC migration, which, as in the case of the leader–follower model, reproduced typical CNCC stream migration akin to that observed in the chick [9].
This finding is particularly concerning for anyone who works, or has worked, within the leader–follower paradigm, as the idea of a phenotypic spectrum or continuum is, in many respects, diametrically opposed to the idea of a rigid leader–follower binary. During the later stages of my doctorate, I found this problem very interesting. How can one reconcile these two seemingly contradictory sets of biological evidence? Can we trust any of the insights gained from the leader–follower model up to now?
New model, new insights
To interrogate this question further, we devised a new mathematical model [11] that drew inspiration from previous ones [3,9] and additional biological findings pertaining to chick CNCCs [10]. In vitro (i.e. in experiments performed outside of the chick embryo) exposure to VEGF causes CNCCs to express genes that drive movement along VEGF gradients (a leader-like behaviour), while a withdrawal of VEGF induces the opposite response, and sees CNCCs begin to express genes associated with, amongst other things, cell–cell adhesion (a follower-like behaviour). From this, one can infer that VEGF drives the phenotypic heterogeneity we observe within CNCC streams (this much is agreed upon). However, this, in turn, suggests that during VEGF-induced transitions, CNCCs likely occupy intermediate phenotypes distinct from both leader-like and follower-like states.
Our new mathematical model, the spring–polarisation model (Figure 2(b)), retains the hybrid description of individual cells coupled to a continuously varying VEGF field that characterised the work conducted at the outset of the Oxford–Stowers Institute collaboration. All cells obey the same behavioural rules, respond to VEGF gradients (though here, they do so through a slightly more biologically detailed mechanism than in previous models), and interact mechanically with their neighbours via a Hookean spring-like interaction potential. As before, each cell samples the local VEGF landscape through the extension of filopodia, using the concentration-dependent detection principle described above.
The central model variable is a cell-specific polarity vector, $\mathbf{m}_i(t)$, which evolves according to the ordinary differential equation
\begin{equation} \label{polarityEquation} \begin{split} \frac{\mathrm{d}\mathbf{m}_i}{\mathrm{d}t} &= \frac{1}{\tau_{\mathrm{sense}}} \left[ \chi_i(t)\,\overset{\scriptscriptstyle\wedge}{\mathbf{u}}_i(t)-\mathbf{m}_i(t) \right]. \end{split} \end{equation}
Here, $\overset{\scriptscriptstyle\wedge}{\mathbf{u}}_i(t)$ is a unit vector pointing towards the strongest detected VEGF gradient, and $\chi_i\in[0,1]$ measures the strength of the chemotactic response. The timescale $\tau_{\mathrm{sense}}>0$ determines how quickly cell polarity (i.e. the cell-specific polarity vector) relaxes towards the VEGF gradient currently detected by a cell. For a sustained cue of a fixed magnitude and direction, $\mathbf{m}_i$ relaxes towards the vector $\chi_i\overset{\scriptscriptstyle\wedge}{\mathbf{u}}_i$; when detectable cues are absent, it decays exponentially at a rate $1/\tau_{\mathrm{sense}}$. Thus, polarity is an exponentially weighted history of the VEGF gradients detected by cells, which, unlike in prior models, builds a notion of memory into VEGF-induced movement.
In previous leader–follower models, CNCC speed was fixed, and the direction of movement was determined by cell–cell interactions and environmental VEGF-gradient sampling. Here, cell speed and direction are jointly determined by these two mechanisms. Neglecting inertia (and hence, invoking an overdamped assumption) and absorbing drag into the force parameters, the deterministic component of the force driving CNCC movement is a linear combination of VEGF-induced polarity and Hookean pairwise interactions:
\begin{equation} \label{cellMovementEquation} \begin{split} \dot{\mathbf{x}}_i &= \mathbf{F}^{\mathrm{pol}}_i + \sum_{j\in\mathrm{N}_i} \left( \mathbf{F}^{\mathrm{spring}}_{ij} + \mathbf{F}^{\mathrm{rep}}_{ij} \right), \end{split} \end{equation}
where $\mathbf{x}_i$ is the position of cell $i$ and $\mathrm{N}_i$ is its set of neighbouring cells. The active chemotactic force is
\begin{equation} \label{activeForceEquation} \begin{split} \mathbf{F}^{\mathrm{pol}}_i &= f_0\,\|\mathbf{m}_i\|\, (\cos\theta_i,\sin\theta_i), \end{split} \end{equation}
where $f_0>0$ sets the force scale and $\theta_i$ is the cell's heading. When polarity is non-zero, this heading gradually aligns with the polarity direction, $\arg(\mathbf{m}_i)$. Polarity, therefore, determines both the strength and the preferred direction of active movement for each cell. The spring-like interaction forces, $\mathbf{F}^{\mathrm{spring}}_{ij}$, represent transient adhesive bonds that form between nearby cells and break when stretched too far. Short-range repulsion, $\mathbf{F}^{\mathrm{rep}}_{ij}$, accounts for the finite size of cells by strongly penalising overlap.
Numerical simulations of the spring–polarisation model show that, although CNCCs enter the domain in a uniformly unpolarised state, an approximate leader–follower phenotypic binary emerges during migration. After an initial transient, the continued exposure of CNCCs at the front of streams to steep VEGF gradients leads to a high degree of polarisation, and thus, a large magnitude of their VEGF-induced polarity vector, $\mathbf{m}_i$. This, in turn, means that their movement occurs primarily according to the VEGF cues that they sense (very much reminiscent of the leader-like behaviour in prior models). All CNCCs behind the leading edge consistently migrate into regions of shallow VEGF gradients, and thus remain unpolarised, meaning that their movement is determined primarily by pairwise cell–cell interactions (follower-like behaviour). As migration continues, this effect becomes increasingly pronounced, until a stable subset of CNCCs at the leading edge of the stream forms a small population that guides the collective. The above video shows one such example of this; the degree to which CNCCs are polarised is denoted by the strength of their green tint, and their direction of polarisation is denoted by the arrow at their centre.
Thus, our model lends credence to the notion of a leader–follower binary; we can all breathe a collective sigh of relief. Our findings do, however, suggest that this may emerge during, rather than prior to, collective CNCC migration in the chick. This prediction has implications for many aspects of collective CNCC migration in the chick, including how streams may respond to environmental perturbations and how migration may break down under atypical conditions.
Outlook
Combining mathematical modelling and experiments in the study of chick CNCC migration has yielded insights that neither discipline could have yielded alone. It is a case in point of the far-reaching potential of interdisciplinary collaboration and frequent dialogue between theoreticians and experimentalists. Many open questions remain regarding the mechanisms of collective chick CNCC migration, and it is highly likely that some of these mechanisms generalise to other systems and biological contexts, including tumour cell migration. In the case of my doctoral work, I hope that the predictions generated by the spring–polarisation model will eventually be tested for experimental consistency (the model is consistent with previous experiments [11]), which could, in turn, inform the development of future mathematical models of CNCC migration and move us closer to new therapeutic interventions for cranial neurocristopathies.
References
[1] Kulesa, P.M., Bailey, C.M., Kasemeier-Kulesa, J.C. and McLennan, R., 2010. Cranial neural crest migration: new rules for an old road. Developmental Biology, 344(2), pp.543–554. https://doi.org/10.1016/j.ydbio.2010.04.010
[2] Alvizi, L., Nani, D., Brito, L.A., Kobayashi, G.S., Passos-Bueno, M.R. and Mayor, R., 2023. Neural crest E-cadherin loss drives cleft lip/palate by epigenetic modulation via pro-inflammatory gene–environment interaction. Nature Communications, 14(1), 2868. https://doi.org/10.1038/s41467-023-38526-1
[3] McLennan, R., Dyson, L., Prather, K.W., Morrison, J.A., Baker, R.E., Maini, P.K. and Kulesa, P.M., 2012. Multiscale mechanisms of cell migration during development: theory and experiment. Development, 139(16), pp.2935–2944. https://doi.org/10.1242/dev.081471
[4] Johnson, S.W.S., 2025. Mathematical models of single and collective cell migration in the cranial neural crest (Doctoral thesis, University of Oxford). https://doi.org/10.5287/ora-mp4jboenv
[5] Wiser, M.J., Ribeck, N. and Lenski, R.E., 2013. Long-term dynamics of adaptation in asexual populations. Science, 342(6164), pp.1364–1367. https://doi.org/10.1126/science.1243357
[6] Berg, H.C. and Purcell, E.M., 1977. Physics of chemoreception. Biophysical Journal, 20(2), pp.193–219. https://doi.org/10.1016/S0006-3495(77)85544-6
[7] Giniūnaitė, R., Baker, R.E., Kulesa, P.M. and Maini, P.K., 2020. Modelling collective cell migration: neural crest as a model paradigm. Journal of Mathematical Biology, 80(1–2), pp.481–504. https://doi.org/10.1007/s00285-019-01436-2
[8] Morrison, J.A., McLennan, R., Wolfe, L.A., Gogol, M.M., Meier, S., McKinney, M.C., Teddy, J.M., Holmes, L., Semerad, C.L., Box, A.C., Li, H., Hall, K.E., Perera, A.G. and Kulesa, P.M., 2017. Single-cell transcriptome analysis of avian neural crest migration reveals signatures of invasion and molecular transitions. eLife, 6, e28415. https://doi.org/10.7554/eLife.28415
[9] Schumacher, L.J., 2019. Neural crest migration with continuous cell states. Journal of Theoretical Biology, 481, pp.84–90. https://doi.org/10.1016/j.jtbi.2019.01.029
[10] McLennan, R., Schumacher, L.J., Morrison, J.A., Teddy, J.M., Ridenour, D.A., Box, A.C., Semerad, C.L., Li, H., McDowell, W., Kay, D., Maini, P.K., Baker, R.E. and Kulesa, P.M., 2015. VEGF signals induce trailblazer cell identity that drives neural crest migration. Developmental Biology, 407(1), pp.12–25. https://doi.org/10.1016/j.ydbio.2015.08.011
[11] Johnson, S.W.S., Baker, R.E. and Maini, P.K., 2026. The spontaneous emergence of leaders and followers in a mathematical model of cranial neural crest cell migration. arXiv preprint arXiv:2601.00374. https://doi.org/10.48550/arXiv.2601.00374