A biased random walk approach for modeling the collective chemotaxis of neural crest cells

Collective cell migration is a multicellular phenomenon that arises in various biological contexts, including cancer and embryo development. ‘Collectiveness’ can be promoted by cell-cell interactions such as co-attraction and contact inhibition of locomotion. These mechanisms act on cell polarity, pivotal for directed cell motility, through influencing the intracellular dynamics of small GTPases such as Rac1. To model these dynamics we introduce a biased random walk model, where the bias depends on the internal state of Rac1, and the Rac1 state is influenced by cell-cell interactions and chemoattractive cues. In an extensive simulation study we demonstrate and explain the scope and applicability of the introduced model in various scenarios. The use of a biased random walk model allows for the derivation of a corresponding partial differential equation for the cell density while still maintaining a certain level of intracellular detail from the individual based setting.


Introduction
Collective behaviour is ubiquitous in the living world and can be observed on various length-scales in many different cells and organisms [12,14,18,37].One example is collective cell migration, which is involved in wound-healing [13], cancer metastasis [13,18,30] and embryonic development [13,31].Notably, both metastatic cancer cells and neural crest cells undergo epithelium-to-mesenchyme transition [13,43].Hence, studying one can also give insights into the other.
During vertebrate embryogenesis, the coordinated migration of neural crest cells (NCCs) plays a crucial role [13,43].NCCs are a multipotent cell population that undergoes epithelial-to-mesenchymal transition during embryonic development [43].In their epithelial state, they reside in the neural tube, 1 arXiv:2310.01294v1[q-bio.CB] 2 Oct 2023 a precursor tissue of the brain and spinal cord.However, in their mesenchymal state, NCCs either migrate individually or collectively through the embryo, differentiate into various cell types, and contribute to the formation of a multitude of developing tissues and organs [44].Successful migration of NCCs involves a diverse array of processes that vary both between species and between different types of NCCs, such as cranial, vagal, trunk, or sacral NCCs [7].
A combination of both in vivo and in vitro experiments have demonstrated that NC migration involves various contributing mechanisms [5,8,11,40].In in vitro studies with zebrafish and Xenopus, the presence of co-attraction, contact inhibition of locomotion (CIL), and spatial confinement were sufficient for the successful collective migration of NCCs [11].Conversely, in vivo investigations indicate that a chemoattractant, such as Sdf1, is required for the guided migration of NCCs in a particular stage of embryonic development in zebrafish [5].Moreover, alternative in vivo studies involving vagal NCC transplantation into the sacral neuroaxis suggest that a chemoattractive cue is dispensable for vagal and sacral NCCs in chick and quail embryos, as these NCCs were capable of colonising the hindgut by following pathways usually taken by sacral cells [8].The reader is referred to [42,43] for reviews on mechanisms involved in NCC migration.
Using mathematical modeling in this context to simulate biological experiments in silico can contribute to the support or explanation of hypotheses and enhance understanding of the mechanisms involved.Previous approaches to model the collective behaviour of NCC comprise, for example, off-latice individual based approaches [22,23,32], cellular Potts models [17], agent-based models [11,19,32,39,46,47] and discrete models with polygonal cell structures [26,27].The proposed models differ not only in their approach but also in the level of biological detail taken into consideration as well as in the type of mechanisms involved in the collective motion.
In the latter mentioned models [26,27], cell-cell contacts mediate the Rac1 and RhoA levels on the cell membrane, which in return influence the direction of movement.This modeling approach captures more details about intracellular processes involved in cell migration, however, the complexity of the model limits to a pure simulation study, as opposed to less involved models that can potentially admit analytical investigation.In the individual-based models, mentioned above, the cell-cell and cellenvironment interactions directly affect the direction of movement without considering the intracellular Rac1 -RhoA bridge.On the other hand, derivation techniques can be applied to acquire a corresponding continuous equation for the cell density.For an extensive review on modelling techniques for collective cell migration the reader is referred to [9,15].The aim of this work is to combine the advantages of a simple random walk model with the consideration of intracellular processes.Hence, we introduce a biased random walk model where the bias depends on the intracellular Rac1 levels surrounding the cell midpoint.By using a random walk approach we maintain the possibility of deriving a corresponding partial differential equation for the cell density following standard methods, e.g.[29,36].In this paper we present a two-dimensional hybrid lattice model describing the continuous evolution of the intracellular Rac1 dynamics and the discrete-in-time migration of the cells biased by the internal differences in Rac1.In Section 2 the mathematical model is introduced, before the simulation setup is presented in Section 3. The numerical simulation results are presented in Section 4 and discussed in Section 6.In certain cases a representation of this model can be derived in form of a system of partial differential equations and the main principle of the derivation is demonstrated in Section 5.

Mechanisms involved on the molecular scale
One of the main mechanisms regulating cell polarity and hence the direction of movement for a motile cell is the intracellular Rac1-RhoA system [21].At locations with higher Rac1 activity, actin polymerization is enhanced and the extension of lamellipodia and filopodia is promoted; in other words, protrusions are formed at the leading edge of the cell [21].Areas with a higher concentration of active RhoA form at the trailing edge of the cell, where it controls myosin phosphorylation and actomyosin contractions that pull the back end of the cell forward [21].The resulting cell-polarity is pivotal for directed cell motility [21].
Additionally, it has been shown that Sdf1 is capable of stabilizing and maintaining the polarity of NCCs induced by cell-cell contacts [40].
It has been observed that most neural crest (NC) cell types tend to migrate collectively, rather than individually [16,42].Short-range chemotaxis, up the gradient of the auto-produced complement factor C3a, has been suggested as a mechanism for maintaining cohesiveness within the NC cell cluster [11].Specifically, when C3a binds to its receptor C3aR, it activates the Rac1 pathway [11], inducing sufficient cell polarization in escaping cells to enable them to return to the group [11].
A counterbalance to the co-attraction is a short-ranged repulsion mechanism between colliding cells, referred to as contact inhibition of locomotion (CIL) [20].This phenomenon, first observed and described in vitro in fibroblasts [1,2] causes cells to alter their migration direction upon collision and move in opposing directions.In Xenopus and zebrafish, NCCs have been shown to exhibit CIL in both in vivo and in vitro studies [10].Additionally, it has been demonstrated that the induction of cell polarity through cell-cell contact is pivotal for the collective chemotaxis in NCCs [40].Mechanistically, on the site of cell-cell contact RhoA increases [10] while Rac1 activity decreases concurrently [40].
Other important factors influencing NCC migration include surrounding spatial constraints or chemorepellents [34].Extracellular matrix components, such as Versican or signalling factors such as ephrins, semaphorins and DAN act as chemorepellents and prevent mixing between separate NCC stream and keep the NCC from invading certain tissue [43].In Xenopus, the repulsive signal Versican is located at the border of the NC and by acting as spatial constraint it enhances NC migration [39].
The following articles provide detailed information on chemotaxis in NCC [33], chemical and mechanical signals involved in NCC migration [34], contact inhibition of locomotion [20] and molecular mechanisms involved in regulation of Rac1 and RhoA [41], respectively.

Model description
In this section, we present a novel hybrid model to describe the collective movement of neural crest cells (NCCs).Our model adopts a biased random walk in discrete time and space, where the internal state variables that govern the bias evolve continuously between two discrete jumps.The distinctive feature of our approach lies in the unique cell shape that we consider.Previous models of collective migration and chemotaxis have employed random walks of point particles [36], which enable the derivation of a corresponding partial differential equation for the total cell density [29].However, this approach fails to capture any spatial heterogeneity with respect to an internal cell state.In contrast, other models such as those presented in [26,27] employ a polygon shape to represent the cell and track internal cell states at multiple boundary locations as well as within the cell.While these models capture more of the complexity of internal cell processes, deriving an equation for the total cell density from these models using standard techniques is challenging, necessitating the use of large scale simulations to gain insight.
In our approach we aim to combine the advantages of using point particle models and spatially extended models.To achieve this, we track the center of each cell, which acts as a point particle, and incorporate additional spatial detail through the implementation of a surrounding membrane composed of adjacent grid points for every cell (as depicted in Figure 1).By doing so, we are able to retain information on the internal cell states, while still having the potential to derive continuous equations for the cell densities.
The key model assumptions are as follows: • Domain: Our model assumes that the cells migrate on a two-dimensional discrete lattice, characterized by a fixed grid size.In certain simulations we will also impose spatial constraints to mimic the effects of surrounding tissue or chemorepellents.Specifically, the confinement is represented as a corridor with walls on three sides that are "soft", allowing cells to cross them and exit the corridor.However, this action incurs a negative impact on the internal variable that drives the directional bias of the cell.Thus, cells are biased to migrate away from regions outside the corridor, to account for the constraints imposed by the surrounding tissue [39].
• Number of cells: A population of N cells is considered and cells are placed on the grid with the position of their midpoints at one of the lattice points.We do not consider any cell division or cell death, therefore the number of cells remains constant.
• Shape of cells: Each cell is characterised by the midpoint, representing the central cell position, and four adjacent points, mimicking the cell membrane, see Figure 1.We assume that the cells do not undergo any growth or division, thus the distance between the membrane point and the cell position is consistently equal to the grid size, and the shape of the cell remains constant.
• Movement of cells: Cells undergo a biased random walk with a fixed velocity on the lattice, meaning they move a distance h (grid size) per time step.The direction of movement is biased by the differences in the internal state variable along one axis, resulting in a tendency for the cells to move in directions with higher concentrations of the internal state.We allow the cells to overlap either partially or entirely, as our model does not account for all the spatial dimensions present in reality.Hence, the occurrence of cells with the same location should be interpreted as cells that are stacked on top of one another rather than cells occupying the same physical space.We do not allow for cell rotation in our model.
• Internal state variable: We consider the intracellular molecule Rac1 as the internal state variable that regulates cell motility.The support of Rac1 for each cell is limited to the directly adjacent grid points of the central cell position and changes over time.The directional bias is influenced by the difference of the internal state value along one axis, where a higher difference in Rac1 leads to a higher probability in moving in the direction where more Rac1 is present.The internal variable is subject to the following extracellular mechanisms: -Increase by chemoattractant Sdf1 : We assume that the chemoattractant's profile is constant in time and independent of the cell population.Furthermore, the profile should be non-negative and continuous.We consider different profiles for different simulation scenarios.
-Increase by co-attraction molecule C3a: The process of co-attraction between cells is governed by the chemoattractant molecule C3a, which is synthesized and released by each cell.We assume that these processes occur at a faster time scale than those related to the polarization of the cell or cell movement.
Accordingly, we consider each cell to have its own stationary C3a-profile, which is radially symmetric and non-negative, and decreases with distance from the cell, e.g. as proposed by [11].Note that although the individual C3a profile does not change in time due to any production, diffusion or degradation, the total C3a profile changes in time due to the movement of the cells.We assume that an individual cell can sense C3a released by neighbouring cells if the cells centers have a distance less than R, where C3a surrounding different cells are treated in an additive way.In other words, the concentration of C3a around an individual cell is positive only within a radius R around the cell midpoint.
-Decrease due to contact inhibition of motion: Upon contact with neighbouring cells, the local Rac1 concentration decreases at the point of contact.As "contact" we consider the overlapping of one or more center or membrane points occupied by a cell, see Figure 2c.
-Natural inactivation: We assume that Rac1 is naturally inactivated at a fixed rate.
-Decrease upon contact with spatial constraints: The spatial confinement representing surrounding tissue is implemented as a chemorepellent that decreases the local Rac1 activity.Consequently, the cells are discouraged from leaving the corridor due to the negative impact of the confinement on their internal variable.

Model formulation
Let x = (x1, x2) ∈ hZ 2 be a point on the two-dimensional discrete lattice with grid size h and let A(x) := {x ± he1, x ± he2} and A 0 (x) := {x} ∪ A(x) denote the set of all adjacent grid points to x, excluding and including x, respectively.Here, ej denotes the j-th unit vector.As for the biological interpretation of these two sets, the former can be seen as the discretized membrane of the cell centered at x and the latter is then the area occupied by the whole cell.
Consider a number of N individual cells located on the lattice and let P i (tn) denote the position of the i-th cell at time tn ∈ τ N with initial position P i (0) = x i 0 , for i = 1, ..., N .Furthermore, let P(t) = {P i (t) : i = 1, ..., N } be the set of all cell positions at time t.We assume that each cell leaves its position within a time interval [tn, tn+1) with probability 1.
The position at time tn of the i-th cell can be written as the sum of all jumps that were made up until that time, where Y i n ∈ {±hej}j=1,2 denotes the outcome of the decision the i-th cell makes for the direction of movement in the time interval [tn−1, tn).
Let C i (x, t) denote the internal state variable of the i-th cell at (x, t) for t ≥ 0. Note, that in comparison to the evolution of a cell position discrete in time, the time variable here is continuous.
Biologically, this variable represents the concentration of the intracellular molecule Rac1 at the cell membrane.The temporal change in this variable is twofold: 1. Every discrete time step, when a jump occurs, the cell midpoint and membrane points change.
Hence, the support of C i (x, t) changes every discrete time step.The support of C i (x, t), i.e. the points where C i (x, t) ̸ = 0, is restricted to the adjacent grid points to the central cell position, 2. Between two discrete jumps, i.e. for t k ≤ t < t k+1 , the internal state will be subject to changes described by the following ordinary differential equation where B(x, t k ) = λ3S i 3 (x, t k ) + λ4 + λ5b(x) and λ1, λ2 represent the activation rate of Rac1 by a chemoattractive and a co-attractive cue, respectively, λ3 and λ5 are the deactivation rates of Rac1 by contact inhibition of locomotion and spatial constraints, respectively, λ4 is the natural deactivation rate of Rac1.The rates λj, for j = 1, . . ., 5, are positive constants.The system is completed by initial conditions where x − Y i k is the previous location of this membrane point.For t = 0 we define In (2), S1(x) ≥ 0 is the profile of the outer chemoattractant Sdf1 which we assume to be constant in time, S i 2 (x, t k ) is the cumulative co-attractive profile for the i-th cell, sensing all the neighbours that are located within a radius R of the i-th cell, i.e.P j (t k ) ∈ P(t k ) \ P i (t k ) such that ∥P j (t k ) − x∥ < R for x ∈ A(P i (t k )) and j = 1, . . ., N .The individual co-attraction profiles are given by with parameters M and w, and S i 3 (x, t k ) acts as a counter of overlapping cells Finally, b(x) describes potential spatial obstacles or constraints.We have that b(x) = 0 wherever cells should be able to move freely and b(x) > 0 at locations of spatial constraints, e.g.where the walls of a corridor are located.
Let T i x→y denote the probability of the i-th cell jumping from position x to position y.We assume that the cells move with constant speed h/τ , which in particular means that remaining at the same position and jumps of size larger than h are impossible and that a cell leaves its position with probability 1 within one discrete time step τ .Formally this means T i x→y = 0 for all y / ∈ A(x) and y∈A(x) The probability of the i-th cell jumping from a position x at time t k to an adjacent position y at time t k+1 will be assumed to have the form where g is an increasing (g ′ (•) > 0) and bounded (limx→±∞ g(x) = ±c) function and positive constants c, α, β, such that α > βc to ensure T i x→y > 0. The function g represents the bias coming from the internal state variable of the cell.The boundedness assumption on g ensures positivity of the probabilities and is biologically reasonable.We choose an increasing function to represent the assumption that a higher difference in the internal state variable leads to an increased polarity and with that an increased probability of moving into a beneficial direction.We will compare two bias functions g(•): 1. Symmetric bias function: We assume that g(•) is an odd function, i.e. that g(−x) = −g(x).
Due to the oddness of the function the positive bias in one direction equals the negative bias in the opposite direction.This leads to a normalization constant 4α in the two spatial dimensions, . Note that the probability of moving vertically or horizontally is equal to 1/2 due to that choice.
This also means that the maximal probability of jumping into any direction is 1/2 in this case.

Asymmetric bias function:
In case a cell is confined on three of its sides, the symmetric bias function might yield unsuitably weighted probabilities since the probability of moving into two of those confined sides is still 1/2.Hence, we also consider an asymmetric bias function that only adds to the probability if there is a positive benefit of moving into a certain direction, see Figure 3 (left panel).

Simulation setup
The initial setup: A number of N cells = 10 cells are initially placed as in Figure 3 (right panel).
The initial configuration aims to represent a densely organised cell cluster that has just undergone epithelial to mesenchymal transition, where cells are still maintaining cell-cell contacts at t = 0 but will be able to move independently in the next instant.The pseudo-code can be found in Section A.
We will differentiate between several main simulation set-ups: • Two different bias functions g(•): 1. Symmetric bias function: We use an odd function and positive constants α1, β1 such that ).The choice of this function inherently represents the assumption that there can be both, negative and positive, directional bias.Note that we will chose a different value for α1 and α2 in both these cases, see Table 2, as the value of α1 in the symmetric case has to be sufficiently large to ensure positivity of the probability.
• With or without spatial confinement: In cases with spatial confinement, the domain within which cells are free to move is a corridor displayed in Figure 3 2), where it is made increasingly difficult to move further away from the corridor with the presence of a chemorepelling substance such as Versican in the surrounding tissue increases with distance to the corridor.This assumption is not based on any biological observations and purely practical, ensuring that the cells move back into the corridor with higher probability.
• With or without co-attraction, contact inhibition of locomotion and natural inactivation of Rac1: In cases without one or more of these mechanisms, the corresponding rate in the ODE (2) is set to 0. We distinguish between different cases with different combinations for (in-)active mechanisms, see Table 1.As a benchmark scenario we will consider the case without CIL, co-attraction, natural inactivation and spatial confinement, i.e. the chemoattractant is the only active component for these simulations.This case can be regarded as the chemotaxis of individual cells without any cell-cell interactions.
3. With a chemoattractive profile that can only be sensed by a cell if the cell is already sufficiently polarized due to cell-cell interactions, i.e. if sufficient Rac1 is present [40].This switch is modelled with a Hill function that is dependent on the local Rac1 concentration where K is the threshold for sufficient cell-cell interaction to sense the chemoattractant Sdf1 and n ≥ 1.
• A strong and weak chemotaxis regime: This will be done by a modulation of parameter λ1 in equation (2).

Parametrization of the mathematical model
We follow [45] and assume that cells have a fixed diameter of 40µm and choose the grid size accordingly so that h = 1 ≡ 20µm.NCC migrate with an average speed of approximately 3µm/ min, on the faster end in a solitary state and slower when they are part of a cluster [45].For this reason, we chose the time step τ = 1 ≡ 7 min so that we obtain a constant speed of 20µm/7 min ≈ 2.86µm/ min.
In reality cells slow down upon contact with other cells and then take some time to regain their regular speed.In this model we do not consider any variation in speed, for simplicity.
The corridors, through which NCCs migrate, vary in width depending on the exact type of the NCC as well as on the species [39].Generally, larger clusters migrate through wider tubes [39].We follow [45] and chose the width of the corridor to be 200 µm, which corresponds to 5 cell diameters.
The number of cells in a group of interacting NCCs also depends on the specific NCC type and species, for computational reasons we simulated a cluster comprised of 10 cells.
For the co-attractive profile we assume a steady state distribution at every time step since the diffusion of C3a acts on much faster timescales than the cell migration.According to [45] this steady state distribution can be described by a Bessel function, and approximated by a decaying exponential with an approximate half maximum length of 110 µm.Adapting to our model, this leads to the choice of d = 8 in equation ( 4).The maximum co-attractive strength and the rates of Rac1 activation and deactivation were adapted from [26].

Statistics for numerical simulation results
To describe the two-dimensional simulation results, each simulation setup was ran Nsim = 100 times.
To quantify the migrated distance we use the mean and standard deviation of the mean x1 position evaluated at a given time t.To quantify the dispersion of clusters we will use the mean root-meansquared-deviation of the clusters where P i j (t k ) is the position of the i-th cell in the j-th simulation at time t k , xj(tk) is the mean position of the cluster in the j-th simulation at time t k and N cells is the number of cells used for one simulation.While other measures can be conceived, we utilise these for their relative simplicity.In the Figure 4 the simulation results are summarized by plotting the measure for migrated distance Mx 1 (t) − Mx 1 (0) against the measure for dispersal MRMSD(t) − MRMSD(0).

Results
We provide a summary figure of the simulation results in Figure 4, where the output of each simulation is represented in terms of (migrated distance, cluster dispersal) after 350 minutes; points positioned further to the right indicate higher migrated distances, while points positioned higher up the axis indicate greater dispersal of the initial cluster.Simulation scenarios A-F are described in Table 1, and g1 and g2 are the two different bias functions used.Figure 5 illustrates typical model output for some representative simulations.Figures 6-8 provide more detailed statistics of certain scenarios.In

Spatial confinement promotes collective cell migration
All simulations were conducted with and without spatial constraints, with the results consistently indicating that confinement facilitates directed migration.The spatial constraints were implemented as triggering a decrease in Rac1, i.e. as an inhibition of cell movement factor.In the absence of spatial constraints or a dominant chemoattractive cue, the initially dense cell clusters dispersed in the majority of cases.Therefore, spatial constraints appear to be not only advantageous but also necessary for maintaining a certain degree of cohesiveness within the migrating cluster, considering the other assumptions we made on the system, especially the assumptions made on the co-attractive strength, see Figure 6.

Asymmetric gain function allows for a strong cell polarity and increases migration distance
We conducted simulations with two different bias functions, g1 and g2, to investigate different abilities to maintain stable cell polarizations.When considering g1, the polarization alternates on average every other time step between pointing in a horizontal or vertical direction, whereas g2 allows for a strong polarization, potentially guiding the cell down the most favorable path.Based on these assumptions, we expected better guided movement with g2 as the probability of moving in the most advantageous direction can exceed 1/2.This expectation was confirmed in all simulation experiments, as shown in Figures 6, 7, and 8, where cell clusters migrated further when g2 was used, sometimes by a considerable margin.The most crucial distinction between using an asymmetric and symmetric gain function is evident in cases 'A' and 'B', where interactive mechanisms are inactive, and the chemotactic cue is the only factor influencing the bias.The asymmetric function allows the bias to dominate the probability of moving in the direction of the chemotactic cue, see Figure 5a, while the usage of g1 leads to less directed movement, see Figures 7 and 8

In the absence of a chemoattractant
Prior studies have presented in silico findings regarding experiments that exclude chemoattractive cues, solely focusing on the impact of cell-cell interactions and potential spatial restrictions [11,39].
Our present study employs analogous experiments to assess the consistency of our model with regards to these circumstances.
In the absence of all cell-cell and cell-environment interaction mechanisms, cells move randomly without any preferential direction, and hence the unnormalized rate of movement in either direction is equal to αi + βigi(0), resulting in a normalized probability of 1/4 for all cells.This scenario serves as the base case (BM), and the results remain the same for both gain functions, as shown in Figure 6.
Mechanisms are added one at a time to demonstrate its influence.The addition of spatial constraints promotes cell migration down the corridor, as demonstrated in case 'A' in Figure 6.The natural inactivation rate of Rac1 reduces the 'memory time' of previous interactions, causing cells to migrate slightly less far when natural inactivation (NI) is active, see case 'B' in Figure 6.The inclusion of co-attraction encourages cells to cluster densely, preventing collective or individual migration, see 'C' in Figure 6.However, the addition of CIL promotes the dispersal of individual cells while still maintaining a certain level of cohesiveness within the cluster, see case 'E' in Figure 6.Omitting the natural inactivation of Rac1 slightly increases the migrated distance, see case 'F' in Figure 6, which is similar to the effect observed in case 'B' compared to case 'A'.When the deactivation rate is set to zero, the Rac1 level retains information regarding previous encounters with neighboring cells or the environment.All of the above observations hold true for both g1 and g2.When the co-attraction mechanism is inactive, the system behavior becomes strongly dependent on the form of gain function used.The gain function g2 allows for strong singular cell polarity, and consequently, cells actively disperse, resulting in an increased migration distance on average but with a large standard   deviation.In contrast, the symmetric gain function g1 restricts cell polarity and hence results in less drastic dispersal.In summary, in the absence of a chemoattractant cells migrate farthest when both co-attraction and contact inhibition of locomotion are active.

Low rate of Rac1 up-regulation by chemoattractant Sdf1
Unsurprisingly, introduction of a chemoattractive cue leads to enhanced migration along the corridor compared to the absence of a chemoattractant, see Figure 7.Note that the chemoattractive profile is constant in the second spatial dimension and, hence, any movement in the second spatial dimension is due to other mechanisms or randomness; we note further that the vertical axis scale differs between Figures 7 and 6.As before, an absence of CIL while co-attraction is active results in high clustering of cells but with the prevention of efficient migration, see case 'C' in Figure 7, and Figure 5b.This observation does not depend on the gain function or chemoattractive profile used.
As previously mentioned, the gain function g2 generates a strong cell polarity and, hence, the bias leading the cells up the gradient of the chemical cue is comparably strong in the absence of other competing interaction mechanisms, see cases 'A' and 'B' in Figure 7a.This observation holds even when the Rac1 activation rate by the chemoattractive profile is dependent on the active Rac1 concentration itself, and the Rac1 activating co-attraction is switched off, see 'A' and 'B' in Figure 7b.
In these cases almost every movement is up the chemical gradient, similar to the case in Figure 5a.
The gain function g1 does not create such a drastic cell polarization and, hence, the average migration down the corridor in the absence of any interaction mechanisms other than the chemical cue is not as efficient here, see 'A' and 'B' in Figure 7a.For both gain functions g1 and g2, adding the co-attraction adds a competing mechanism to the chemical cue and as a result the clusters migrate less far, see cases 'E' and 'F' in Figure 7a.
The most interesting result in this setting is obtained from using the gain function g1 and the Hill-type response to the chemoattractive profile.Here, the cells perform best when both the coattraction and co-repulsion mechanism are active, see cases 'E' and 'F' in Figure 7b.The both cases 'E' and 'F' perform similarly, although the clusters move marginally further if there is no natural inactivation.An example of one realization of case 'E' is shown in Figure 5c.

Higher rate of Rac1 up-regulation by chemoattractant Sdf1
When the Rac1 activation rate by the chemoattractant is increased further, the migration distance increases significantly in case 'C' but only marginally in cases 'E' and 'F', see Figure 8a.In cases 'A' and 'B' the migratory distance did increase only when g1 was used as it was already at its maximum in the g2 case, see Figure 8a.In case 'C' co-attraction and strong chemotaxis together manage to move the tight cluster further down the corridor, see Figures 8a and 8b.
When the Hill-type response to the chemoattractant, i.e.S1(x) = f2(x), and the symmetric gain function g1 are used, the combined effect of CIL and co-attraction causes again a stronger migration, compared to cases 'A' and 'B', but is now similar to when co-attraction only was used, see Figure 8b.this, we provide a proof of concept analysis in which we derive the governing continuous PDE system for the evolving cluster distribution.The subsequent exploration and analysis of the derived system is left for future study.
While we are interested in finding the position of the i-th cell at time t, for a continuous description of the process it is natural to consider the probability distribution of the i-th cell's location at time t.
In a completely unbiased setting in two dimensions, the probability of jumping in either of the four directions would be 1/4 and in the continuous limit we would obtain a pure diffusion equation for the time evolution of the cells.In our setting we assume that the probability of jumping in a direction is biased by an internal state variable C i (x, t) that is non-zero on the cell membrane for x ∈ A(P i (t)).
Let p i (x, t) denote the probability of the i-th cell being located at position x at time t ≥ 0. Initially, we have p i (x i 0 , 0) = 1 and p i (x, 0) = 0 for all x ̸ = x i 0 .The change in probability of the i-th cell being located at x at time t + τ is then where T i x→y is the probability of the i-th cell jumping from its current position x to an adjacent position y ∈ A(x) within a time interval [t, t+τ ), as defined in Section 2. In [36] different dependencies of the transition rates on the bias inducing variable are discussed and corresponding PDEs are derived.
When following the individual trajectory, the position of the centre of a cell and the associated membrane points is clearly defined by P i (t) and A(P i (t)).In the continuous description, while the support of the probability p i (x, t) will initially be just a point, the support of both the centre of the cell and the membrane points will disperse, see Figure 9.In Figure 9 the probability distributions of a single cell (red dots) and its membrane (green dots) are displayed at times t = 0 and t = τ .The green point in the middle of the right hand side of Figure 9 holds membrane information associated with all the four adjacent potential cell locations.Therefore, in order for the master equation formulation to make sense we have to specify the cell midpoint and membrane points associated with it.
Let H i y (P i (t), t) := C i (y, t) denote the value of the internal state variable at (y, t) associated with the midpoint P i (t) of the cell.Then where the first argument of H i is the cell centre and the subscript y is the considered membrane point, and the distance ∥x − y∥ = h.Assuming that τ and h are small, using in (5) the Taylor series expansion about (x, t), and taking the limit as h, τ → 0 such that lim h,τ →0 (h 2 /τ ) =: D is positive and finite, yields Note, that assumption lim h,τ →0 (h 2 /τ ) < ∞ implies h/τ → ∞, hence the effective speed is infinitely large and the solution p i (x, t) of ( 6) is non-zero at positions that are very far away from the starting position after a very short amount of time.Hence, the PDE is not necessarily a good approximation of the biased random walk for small t unless one ignores very small densities.
The master equation for internal state variable at the membrane point y corresponding to the cell where F represents the change of Rac1 in a time interval [t, t+τ ) defined by the mechanisms described in Section 2. The Taylor expansion around (x, t), combined with the limit as h, τ → 0, yields where Si 2 (x, t) and Si 3 (x, t) are smooth interpolations of the functions S i 2 (x, t) and S i 3 (x, t).Equation ( 6), together with the equation for H i , provides the continuous description for the dynamics of the probability p i (x, t) for a cell to be at position x at time t, when initially located at x i 0 .The interactions between cells in a domain are included in Si 2 and Si 3 and the impact of the environment is modeled by S i 1 .While this represents a complicated and coupled system of equations, analysis and numerical simulations for continuous equations may still be more efficient than the discrete model.This lies out of scope of the present article, where the primary aim was an analysis of the underlying random walk model and to provide a 'proof of concept' for the continuous derivation.For further comment, we refer to the discussion.

Discussion
In this paper we have formulated a random walk model for collective cell movement dynamics, motivated by instances of NCC migration during embryonic development.NCCs achieve migration in a variety of ways -from individuals to coordinated clusters, according to cell type and species -and a wide variety of guiding factors have been proposed [38]: extracellular chemoattractant gradients, ligand-mediated co-attraction between cells, contact-induced cell-to-cell repulsion (contact inhibition of locomotion), and confinement that results from inhibiting cues in the surrounding extracellular matrix.Motivated by populations of cranial Xenopus NCC cells that migrate as a collective group [40], we have incorporated these different factors into a random walk model and explored the extent to which different combinations of guiding factors induce migration in a targeted direction, while maintaining a relatively compact configuration.
The model includes a simplistic description of intracellular signalling, whereby each cell is associated with a set of evolution equations for the levels of membrane-localised Rac1 activation that determine the movement probability into different directions; each evolution equation is positively or negatively regulated according to local guidance cues.This level of detail allows certain subtleties to be incorporated, such as, for example, the suggestion that the chemoattractant only stabilizes and reinforces the existing Rac1 protrusions created by cell-cell contacts in cranial NCC in Xenopus [40].In the light of such findings, the model setting that includes a Hill-type response to the chemoattractive profile is perhaps more realistic, as it assumes the chemoattractive cue could only be sensed at a membrane location if Rac1 is sufficiently activated at that point.
The inclusion of intracellular signalling also allows persistence of polarity, where the cell maintains a preferred direction of motion over long time intervals, to emerge naturally.First, slowing the rate of Rac1 inactivation will act to maintain a high Rac1 level in a certain direction, potentially long after the initial polarising event.Coupled to an asymmetric bias function that curbs switching the axis of movement direction, persistent straight line movements can arise.Thus, according to the length of memory encoded by Rac1 inactivation, model dynamics can switch between an uncorrelated and correlated random walk.Note that while the asymmetric bias function allows strong bias, it is potentially unsuitable for cases where weak individual migration towards a chemical attractor is observed.The symmetric bias function instead introduces more randomness, and subsequently requires more cell-cell and cell-environment interactions for efficient collective migration.
In scenarios where cell polarity is less persistent, the combination of CIL and co-attraction, along with the Hill-like Rac1 activation rate by the chemoattractant, performs better than the other scenarios.Co-attraction is the only included cell-cell interaction that can lead to clustering, while CIL is the sole cell-cell interaction that promotes stronger dispersion.Despite this, scenarios in which CIL is included perform better than cases without CIL, suggesting that it might play an important role in promoting robust directed migration.We note that when the chemoattractive cue is relatively strong or, alternatively, when the response rate to the chemoattractant is relatively large, cell-cell interactions become less relevant.Effectively, chemotactic behaviour dominates the dynamics.
Despite our inclusion of some detail at the intracellular level, there are many further aspects that can be considered to refine the model if needed.One potential direction would be to, in addition to the Rac1 dynamics, also consider the RhoA evolution, similar to [26,27].Another alternative would be to incorporate the degradation of the chemoattractant by cells, similar to [22,23,24].This mechanism, in combination with the Hill-like switch for sensing the chemoattractant, might lead to emergent leader-follower behavior since only the outer cells of the cluster would be able to sense the chemoattractant, similar to findings in [22,23,24].If the decision is dependent on the absolute Rac1 level, outer cells would become more confident in their direction, and the other cells would follow due to the co-attraction mechanism.
It is noted that a number of other agent-based models have been developed to describe NCC (and other) cell migration processes.Single site lattice models of cellular automaton type (for examples, [4,28]) offer a relatively simple and computationally efficient framework: a cell occupies a lattice site and probabilistically moves (or proliferates) into an adjacent site according to a set of rules based on the neighbourhood state, e.g.occupancy and interactions with the nearest neighbours.A particular advantage of these models lies in an established framework that allows scaling to a continuous partial differential equation for the density distribution, permitting further analysis, e.g.[35].A disadvantage lies in their relatively course scale.Other NCC models include finer structure, such as vertex-based approaches where cells form evolving polygonal objects [26] or Cellular Potts Model approaches where a cell occupies a connected set of lattice sites [39].These models also have their disadvantages, in particular reliance on simulation instead of analytical insight.To walk a middle path, the model here has assumed that each cell spans 5 lattice sites to form a "+" configuration, where the centre represents the cell centre and ODEs are associated to each neighbouring site to describe the internal dynamics.By retaining the underlying framework of a biased random walk on a lattice, the possibility of course-scaling to a continuous model is saved and we have shown how this can be done, at least in principal.Notably, the resulting PDE for the evolving probability distribution for a cell's position echoes the taxis-type equation obtained for position-jump random walk models in the presence of a chemical signal [36], where the directional bias is proportional to the gradient of a continuous approximation to the internal state variable.The equation for this internal state variable is also represented by a diffusion-advection type equation, where the transport terms for the internal state variable now arise due to the translocation of membrane positions each time a cell moves.
Further analysis of the continuous problem is certainly a point for future effort, but the intricacy of the derived system (coupled systems of diffusion-advection) demands a dedicated study.In particular, we note that the equation for the internal state variable contains "self-gradient-following" terms that could potentially result in negative diffusion and a loss of regularity.Other points to consider for this study would be the state space for the internal variable.The simplest option, which was performed here, involved collapsing the internal state space variable to the central point of the cell: effectively, the distinct spatial structure contained within the four membrane locations is lost, merged into a single PDE for the internal variable.An alternative and interesting option would be to derive four separate internal state PDEs, each corresponding to a different membrane location.While the model will inevitably become more complex, such an approach could allow the spatial structure at the cellular scale to be retained at the macroscopic or continuous scale.

Figure 1 :
Figure 1: Visualization of a single cell on a two-dimensional lattice: Cell center (grey dot); membrane points (green crosses); Area occupied by a cell (blue circle); Potential locations for jumping (light red dots) (a) Chemoattraction (b) Co-attraction (c) Contact inhibition of locomotion (d) Spatial constraints

Figure 2 :
Figure 2: Schematic of the four mechanisms incorporated in the model, each displayed individually in the absence of the others.Area occupied by one cell is displayed with a color gradient to emphasis the Rac1 polarization, Blue arrows indicate most probable direction of movement.Influences of a) Chemoattractant, b) Co-attractive C3a profiles (blue circles), c) Contact inhibition of locomotion and d) spatial constraints on the internal Rac1 states.

Figure 3 :
Figure 3: (Left) Different gain functions, α 1 = 0.6, α 2 = 0.1, β 1 = β 2 = 1/π.(Right) Initial nondimensionalized configuration in two spacial dimensions.The blue dots mark the initial positions of the 10 cells, the mean position is marked with a red dot.The boundary of the corridor is depicted by red lines located at x 1 = 20, x 2 = 0, x 2 = 10.Note, that only some simulations were conducted in the domain with boundaries, in other cases there was no spatial constraint.
(right panel).The area outside the corridor, with boundary marked by red lines, mimics the presence of chemorepellents or spatial constraints.Note that the cells are technically able to move across the boundary since the presence of corridor walls makes it only less likely but not impossible for them to move there.For these simulations we will use the following corridor function b(x): 20, for x1 < 20, x2 − 10, for x2 > 10, −x2, for x2 < 0, in equation (

Figures 10 ,
Figures 10,11 and 12 we display the temporal evolution of the key measures presented in the summary plot Figure 4.More detailed descriptions follow below.
(a) Individual simulation without cell-cell interactions, large λ1 (b) Individual simulation with inactive CIL, active co-attraction, small λ1 (c) Individual simulation with active CIL, active co-attraction, small λ1

Figure 5 :
Figure 5: Individual realizations of simulations with S 1 (x) = f 1 (x), spatial confinement and gain function g 2 .Initial conditions and positions at t = 350 min are displayed by grey and black dots, respectively.The trajectory of an individual cell is visualized by a purple line.Note that the dimensions are stretched in the vertical direction, which assists visualisation of individual cell positions and tracks.

Figure 6 :
Figure 6: Mean x 1 -position at t = 350 min in different scenarios with or without co-attraction (CoA), contact inhibition of locomotion (CIL), natural Rac1 inactivation (NI) and spatial constraints (SC) in the absence of a chemoattractive substance.

Figure 9 :
Figure 9: Left: Initial position of cell midpoint (red) and membrane points (green) with associated probabilities.Right: After one step in a purely random walk, cell midpoint positions (red) and membrane points (green) with associated probabilities.