Model description
We use a 2D phase-field model17,20,21 to study the chemotactic motion of a multicellular system on a planar solid substrate. Our model describes the chemoattractant dynamics in the extracellular environment, as well as the motion and deformation of each cell caused by chemotaxis and cell-cell interactions. We denote the problem domain as \(\Omega\) and the chemoattractant concentration as \(q({\varvec{x}},t)\). We consider a system comprised of N cells and define the location and geometry of the \(\alpha\)-th cell with the phase field \(\phi _\alpha ({\varvec{x}},t)\), which takes the value \(\sim\)1 inside the \(\alpha\)-th cell and \(\sim\)0 elsewhere (see Fig. 1A). The evolution of the system is controlled22 by the energy functional \({\mathcal {F}}[\phi _1,…,\phi _N]=\sum _{\alpha =1}^N\mathcal {F_\alpha } [\phi _1,…,\phi _N]\), where$$\begin{aligned} \mathcal {F_\alpha }[\phi _1,…,\phi _N]= & \int _\Omega \gamma \left( \frac{G(\phi _\alpha )}{\epsilon ^2} + |{\nabla \phi _\alpha }|^2 \right) \,\textrm{d}\Omega +\frac{\lambda }{6A_\alpha }\left( A_{\alpha } – \int _{\Omega }\phi _\alpha ^2(3-2\phi _\alpha ) \textrm{d}\Omega \right) ^2 \nonumber \\+ & \int _{\Omega } \textit{g}\frac{15}{\epsilon ^2}\sum _{\beta \ne \alpha }^N \phi _\alpha ^2\phi _\beta ^2\,\textrm{d}\Omega . \end{aligned}$$
(1)
The functional \({\mathcal {F}}_\alpha\) accounts for the membrane surface tension, cell volume (area in 2D) conservation, and cell-cell repulsion, whose strengths are controlled, respectively, by the parameters \(\gamma\), \(\lambda\), and g. The surface tension contribution makes use of the double-well potential \(G(\phi _\alpha )=30\phi _\alpha ^2(1-\phi _\alpha )^2\); see14. The parameter \(\epsilon\) is proportional to the thickness of the diffuse interface (see Fig. 1A) and \({A}_\alpha\) is the area of cell \(\alpha\), which may change over time to capture cell growth, deformation and division. For more details on our modeling approach to cell growth, division, and death see the Supplementary Information. If we denote the velocity of the \(\alpha\)-th cell as \({\varvec{u}}_\alpha\), we can express the evolution equation for \(\phi _\alpha\) and the velocity \({\varvec{u}}_\alpha\) as$$\begin{aligned} \frac{\partial {\phi _\alpha }}{\partial {t}} + {\varvec{u}}_\alpha \cdot \nabla \phi _\alpha = -\Gamma \frac{\delta {\mathcal {F}}}{\delta \phi _\alpha }, \quad \quad \quad \quad {\varvec{u}}_\alpha = {\varvec{u}}_\alpha ^A + \frac{60g}{\xi \epsilon ^2}\int _{\Omega }\phi _\alpha \nabla \phi _\alpha \sum _{\beta \ne \alpha }^N\phi _\beta ^2 \,\text {d}\Omega , \end{aligned}$$
(2)
where the parameter \(\Gamma\) controls the overall kinetics of the driving forces in the functional \({\mathcal {F}}\), \(\xi\) is a constant that represents friction between cell and substrate, \({\varvec{u}}_\alpha ^A\) is the active velocity, which we assume uniform in space, and \(\frac{\delta {\mathcal {F}}}{\delta \phi _\alpha }\) is the variational derivative of \({\mathcal {F}}\) with respect to \(\phi _\alpha\); see Methods. Thus, \({\varvec{u}}_\alpha\) is uniform in space and has two components: a passive velocity that emanates from cell-cell interactions (\(\xi\)-term, see22,23 for more details) and an active velocity \(({\varvec{u}}_\alpha ^A)\) that describes a persistent random walk biased by the chemoattractant. In particular, the direction of motion of cell \(\alpha\) is biased depending on the characteristic gradient and average concentration of chemoattractant detected by the cell (more details in Methods). The evolution equation for the extracellular chemoattractant is written as$$\begin{aligned} \frac{\partial (\varphi q)}{\partial t} = \nabla \cdot (D_q\varphi \nabla q) + \delta _\varphi b_q – \delta _\varphi K_q\frac{q}{K_d+q} – \varphi r_q q , \end{aligned}$$
(3)
where the function \(\varphi ({\varvec{x}},t)\), defined as \(\varphi = 1 – \sum _{\alpha =1}^N \phi _\alpha\), localizes the extracellular environment and \(\delta _\varphi ({\varvec{x}},t)\), defined as \(\delta _\varphi = \frac{30}{\epsilon }\varphi ^2(1-\varphi )^2\), localizes the cells membrane (see Fig. 1A). Eq. (3) models chemoattractant diffusion through the extracellular environment, production and consumption on the cells membrane, and degradation, which are controlled by the diffusion coefficient \((D_q)\), and the production \((b_q)\), consumption \((K_q)\), and degradation \((r_q)\) rates, respectively. \(K_d\) is the disassociation constant. Note that the different transport and biochemical processes are localized to the extracellular environment or the cell membrane by including the variable \(\varphi\) or \(\delta _\varphi\), respectively, into the corresponding terms in Eq. (3). In the limit \(\epsilon \rightarrow 0\), the localization approach is equivalent to a system of PDEs on moving domains with variationally consistent interface conditions24.Fig. 1(A) Diffuse-interface approach. We localize cell \(\alpha\) and the extracellular environment with the phase field \(\phi _\alpha\) and the marker \(\varphi = 1 – \sum _{\alpha } \phi _\alpha\), respectively, which take the value 1 inside the corresponding region and 0 outside. We localize the cells membrane with the marker \(\delta _\varphi\) (see the main text), which is maximum at \(\phi _\alpha = 0.5\) and decreases to 0 as \(\phi _\alpha\) approaches 0 and 1. (B) Clustering strategy. We group multiple distant cells by assigning them to a single cluster field \(\mathbb {C}_i=\sum _j\phi _{i,j}\), where index j identifies cells within cluster i. Cells with the same color in the figure belong to the same cluster. The shaded region surrounding each cell (i, j) represents the computational window \(\Omega _{i,j}\), where Eq. (2) is solved. (C) Reallocation strategy. When two cells in the same cluster approach (cells (1, 1) and (1, 3) in the figure), we switch the allocation of one of them (cell (1, 3)) with a cell from a different cluster (cell (2, 1)) such that the computational windows of cells in the same cluster do not overlap after reallocation (note that cell (3, 1) cannot be transferred to cluster 1).
Dynamic cluster field
Conventional multicellular phase-field modeling strategies involve the solution of a separate equation (\(\phi _\alpha\)) for each cell \(\alpha\), which limits the total number of cells N due to high computational costs. Here, we bypass this limitation by proposing a novel algorithm based on two premises:
(i)
We only need to solve Eq. (2) near cell \(\alpha\) because \(\partial \phi _\alpha /\partial t = 0\) far from the cell.
(ii)
We can use a single phase field to represent multiple cells that are distant from each other.
This is true as long as we keep track of the area-conservation term and velocity \({\varvec{u}}_\alpha\) of each cell \(\alpha\) in Eq. (2). However, when two cells are close to each other, they must be captured by different phase fields in order to compute cell-cell interaction forces (g-terms in Eq. (2)).Premise (i) allows us to solve Eq. (2) on a computational window\(\Omega _\alpha \subset \Omega\) defined as the circle of radius \(R_\alpha +1.5\epsilon\) centered at cell \(\alpha\), where \(R_\alpha\) is the equivalent radius of cell \(\alpha\). Although \(\Omega _\alpha\) is usually a very small portion of \(\Omega\), this definition guarantees that \(\Omega _\alpha\) encloses cell \(\alpha\) and \(\phi _\alpha \approx 0\) outside \(\Omega _\alpha\) (see Fig. 1B). Premise (ii) allows us to group distant cells into a single cluster field\(\mathbb {C}_i({\varvec{x}},t)\). Cells within cluster i must maintain a minimum pairwise distance, which we achieve by preventing the overlap of their computational windows. To define our formulation in terms of cluster fields we rename each phase field \(\phi _\alpha\) as \(\phi _{i,j}\), where index i identifies clusters and index j identifies cells within each cluster. If we denote the number of cells in each cluster i as \(N_i\), we can define the i-th cluster field as \(\mathbb {C}_i = \sum _{j=1}^{N_i}\phi _{i,j}\). Note that \(N_i\) and the total number of clusters \((N_c)\) and cells (N) may change due to rearrangements in the spatial distribution of cells, cell division, and cell death, while \(\sum _{i=1}^{N_c} N_i = N\) always holds true. We also rename the computational windows \(\Omega _\alpha\) as \(\Omega _{i,j}\), in consistentency with the \(\phi _{i,j}\) notation, and group these computational windows as \(\Omega ^{\mathbb {C}}_i=\cup _{j=1}^{N_i}\Omega _{i,j}\). Thus, we can replace \(\phi _\alpha\) by \(\mathbb {C}_i\) in Eqs. (1) and (2) and solve Eq. (2) in the computational window \(\Omega ^{\mathbb {C}}_i\), which reduces the number of PDEs from N to \(N_c\). To compute \({\varvec{u}}_\alpha\) and the area-conservation term in Eq. (2) for each cell \(\alpha\), we track the center of each cell (and, hence, \(\Omega _\alpha\)), which permits us to use \(\mathbb {C}_i\) to compute these terms by restricting the integrals to \(\Omega _\alpha\) (more details in Methods). Finally, the marker \(\varphi\) in Eq. (3) is rewritten as \(\varphi =1-\sum _{i=1}^{N_c}\mathbb {C}_i\).This cell grouping strategy fails when cells within the same cluster approach each other because cell-cell repulsive forces cannot be computed, which would violate premise (ii). To tackle this issue, we implement a reallocation strategy where we reallocate one (or more) of the approaching cells to a different cluster. This procedure involves the cluster exchange between one of the approaching cells and a distant cell from a different cluster, such that the computational windows of cells in the same cluster do not overlap after reallocation (see Fig. 1C). When cluster exchange is not possible due to the spatial distribution of cells, we create an additional cluster field to capture all pairwise cell-cell interactions and modify \(N_c\) accordingly. Our reallocation procedure maintains cell position and velocity. We refer to the coupled cell grouping and reallocation strategy as Dynamic Cluster Field (DCF), which effectively reduces the computational cost of high-resolution simulations of multicellular systems with a large number of cells. The DCF method extends and generalizes the bounding box and reassignment algorithms used in the context of grain growth25,26,27 by redefining the dynamics in terms of the cluster field [see Eq. (10)], which is critical for an effective calculation of the individual cell velocities and areas (see Methods). For more details about DCF and the numerical implementation, see Methods.
Cell co-attraction and proliferation of PC12 cells
PC12 is a cell line that originates from a pheochromocytoma of the rat adrenal medulla. These cells are widely used in neurobiological research because they are homogeneous, exhibit rapid growth and differentiation, and secrete dopamine28,29. Understanding chemotaxis in PC12 cells is important in neurodegenerative diseases29. Here, we focus on the experiments in30, which study the aggregation and proliferation of PC12 cells in a Petri dish via self-generated chemoattractant. The aggregation process depends on the complex dynamics of the guiding signal in the extracellular space, including how cells push the chemoatractant toward the center of the aggregate as they assemble. Thus, a quantitative understanding of this experimental system requires a detailed description of the chemoattractant dynamics accounting for the time evolution of the extracellular space. We demonstrate the ability of our model to quantitatively interrogate this non-confluent multicellular system. We perform 2D simulations on a periodic rectangle of size \(720 \times 550 \upmu \!{\hbox {m}}^{2}\). We initially place 39 circular cells grouped in 13 cluster fields. All cells have an initial radius of \({6}\upmu \hbox {m}\) and are placed at random locations, forming a few cell aggregates. The cells are placed such that the number of cells per unit area and the aggregate size distribution match the experiment. We also assume that, at the beginning of the simulation, all cells have the same physical and biochemical properties but different ages. The initial cell ages are sampled from a uniform random distribution in the interval \([0,24\text { hours}]\). We run 4 simulations with different samples of the initial age and initial location of the cells. In our simulations, the extracellular environment does not contain any chemoattractant at the initial time. After the simulation starts, the cells secrete chemoattractant through their membranes. The internalization of chemoattractant by the cells is neglected. Our model considers cell growth, division, and death using information specific to PC12 cells as reported in30. After cell division, the new cells are in an early-growth stage for 6h and cannot divide. After the early-growth stage, cells may undergo division until they reach an age of 18h. If no division has occurred at this point, cells enter the apoptosis stage, in which they neither grow nor divide. Cells die at age 24h. During the early-growth and apoptosis stages, the cells gradually acquire and lose their ability to sense and secrete chemoattractant. For a detailed description of the cell growth, division, death, and chemotaxis modules of the algorithm, see the Supplementary Information.Fig. 2Cell aggregation driven by autologous chemotaxis. (A) Cells and chemoattractant distribution at \(t={12}{\hbox {h}}\) and \(t={24}{\hbox {h}}\). The unit of q in the plot is nM. (B) Zoom-in showing the merging of three aggregates in inset i (\(t={12}{\hbox {h}}\)) to form one big aggregate in inset ii (\(t={24}{\hbox {h}}\)). We display phase-field ranging from 0.3 to 1 to effectively visualize the chemoattractant distribution in the extracellular region, where \(\varphi\) spans from 0.7 to 1. (C) Aggregate size distribution at \(t={24}{\hbox {h}}\). Experiment observations in red and simulation results in green. The aggregate size is divided in 3-cell intervals. The error bar indicates the minimum and maximum values observed in our simulations. We set \(D_q = {70}{\upmu \!{\hbox {m}}^2\,{\hbox {s}}^{-1}}\), \(b_q={18.8}{\hbox {nM}\,{\upmu \!{\hbox {m}}\,{\hbox {s}}^{-1}}}\), \(r_q={0.108}{\hbox {s}^{-1}}\), \({\overline{u}}_A={0.04167}{\upmu \!{\hbox {m}}\,{\hbox {s}}^{-1}}\), and \(\overline{\Delta T}={60}\hbox {s}\).Figure 2A shows the cells distribution and chemoattractant concentration at \(t={12}\hbox {h}\) and \({24}\hbox {h}\) for one simulation. The results indicate that aggregate formation and growth are caused by 3 main mechanisms: (1) Higher chemotactic gradients generated by cell aggregates attract isolated cells. (2) Merging of different aggregates driven by chemotaxis (the 3 small aggregates in Fig. 2B, inset i form the larger aggregate in Fig. 2B, inset ii). (3) Division of cells within an aggregate. As aggregates form, peaks of chemoattractant concentration emerge at their centers (see Fig. 2A, inset ii). The chemoattractant concentration at the center of the aggregate increases as the aggregate size grows because more cells release chemoattractant in a small region and the cells within the aggregate slow down chemoattractant diffusion through the extracellular environment. This phenomenon leads to higher chemoattractant gradients pointing towards the aggregate center, which hinders the detachment of cells from the aggregate. Large aggregates exhibit minimal movement because of the same phenomenon; while small aggregates and isolated cells are highly motile until they merge with a larger aggregate. We suspect that the absence of adhesion combined with the reduced ability of cells to sense chemotactic gradient during quiescent period following division might have contributed to the gaps between the cells in the aggregate as observed in Fig. 2(A). Most cells that are isolated at \(t={24}\hbox {h}\) (see Fig. 2A) have previously detached from a cluster during the early-growth stage. During the early-growth stage, cell motion transitions from a persistent random walk to chemotactic migration, as cells develop the ability to sense chemical gradients. Thus, newly born cells can move away from the aggregates even in the presence of opposing chemotactic signals. The possibility of resolving chemoattractant dynamics in the time-evolving extracellular space is a unique feature of our algorithm that enables understanding of the complex interplay between the cell location and the attractant transport dynamics that occur in closely-packed cell aggregates.Fig. 2C shows the aggregate size distribution at \(t={24}\hbox {h}\) from the experiments (red) and simulations (green), where we also indicated the minimum and maximum values in our simulations with a blue line. The simulation results are in good agreement with the experimental data30 despite our simulations displaying slightly larger aggregates than the experiment. We believe this mismatch can be reduced by improving the stochastic model for cell division (more details in Supplementary Information). We also observe that the total number of living cells N doubles by \(t\sim {12}{\hbox {h}}\) and quadruples by \(t\sim {24}{\hbox {h}}\), suggesting a 12-hour cell division cycle in agreement with experimental observations. Mov. S1 shows the time evolution of the proliferation and co-attraction of cells. The stationary cells in Mov. S1 represent the dead cells. Our results show that our model effectively captures aggregate formation through the combined processes of proliferation and chemotaxis of PC12 cells.
Multicellular migration through self-generated chemoattractant gradients
Another situation in which cells play an important role in reshaping the distribution of chemoattractants is when they create chemical gradients by breaking down the extracellular chemoattractant6,31,32. This mechanism enables cells to collectively migrate over long distances in a robust manner33,34,35. Several experiments33,36,37 analyze the migration of non-confluent multicellular systems driven by self-generated chemotactic gradients. In particular, experiments in33 use Dictyostelium discoideum (DD), a type of cell widely used to study chemotaxis due to their strong sensing and motile abilities —DD’s can detect gradients of 1% difference in chemoattractant concentration across their membrane38.Fig. 3Self-generated chemotaxis of multicellular systems. Evolution of 400 cells in an initial chemoattractant concentration \((q_0)\) of (A)\({0.1}{\upmu \hbox {M}}\) and (B)\({1}{\upmu \hbox {M}}\). Cells initially occupy the left-hand 250\({\upmu \hbox {M}}\) of the domain and we assume a null q concentration in that region. The top half of each panel shows the cell distribution and chemoattractant concentration at \(t={7.166}{\hbox {h}}\). The bottom half of each panel shows the y-averaged chemoattractant concentration \((|q|_y)\) at \(t={0}{\hbox {h}}\) (red dashed line) and \(t={7.166}{\hbox {h}}\) (red solid line) in our simulations, as well as the fraction of cells along the x-axis in our simulations (purple solid line) and the experiment33 (green dotted line) at the final time. We ignore the contribution of cells present between the left boundary and \(x = {100}\upmu \!{\hbox {m}}\) while computing the fraction of cells to incorporate a well-like region. Green bars indicate the minimum and maximum cell fraction observed in the experiment. Mov. S2 and Mov. S3 display the self-generated chemotaxis of cells for \(q_0 = 0.1\) and \({1}\upmu \!{\hbox {m}}\), respectively. We set \(D_q = {12.5}{\upmu \!{\hbox {m}}^2\,{\hbox {s}}^{-1}}\), \(K_q = {16.5}{\hbox {nM}\,\upmu \!{\hbox {m}}\,\hbox {s}^{-1}}\), \(K_d = {50}\hbox {nM}\), \({\overline{u}}_A={0.072}{\upmu \!{\hbox {m}}\,{\hbox {s}}^{-1}}\), and \(\overline{\Delta T}={20}\hbox {s}\).Previous attempts to simulate these experiments have resorted to agent-based models, in which cells have a fixed geometry and are permeable to the chemoattractant. The assumption that cells are permeable to the attractant greatly simplifies the simulation because the computational domain for chemoattractant dynamics has a simple geometry that is fixed in time but misses key elements of the cell-attractant interplay. Here, we use our model to study the experiments in33 and analyze the self-driven chemotaxis of DD’s under different chemoattractant levels. The duration of the experiments in33 is approximately \({7}{\hbox {h}}\). This short time permits us to neglect cell growth, division, and death. We also neglect chemoattractant decay and production (\(r_q=0\) and \(b_q=0\) in Eq. (3)) and consider consumption only. We study a rectangular planar surface of \(2100 \times {400}{\upmu \!{\hbox {m}}^2}\) and we initially place 400 cells (grouped in 10 cluster fields) of radius \({5}{\upmu \!{\hbox {m}}}\) at random locations in the 250-\({\upmu \!{\hbox {m}}}\) strip in contact with the left boundary. We assume that the initial chemoattractant concentration is zero in that 250-\({\upmu \!{\hbox {m}}}\) strip and is constant in the rest of the domain with concentration \(q_0\) (see red-dashed lines in Fig. 3). We do not monitor the activity of cells within the initial 100-\({\upmu \!{\hbox {m}}}\) length from the left boundary to create a well-like area as described in33. We consider periodic boundary conditions in the y-direction and fixed q and \(\mathbb {C}_i\) on the left and right walls. We scale the parameter values to perform the simulation on a smaller domain than the experiment (see Supplementary Information). To analyze the influence of the chemoattractant level, we run two simulations with different \(q_0\). In particular, we consider \(q_0={0.1}{\upmu \hbox {M}}\) and \(q_0={1}{\upmu \hbox {M}}\), and we plot the results in Figs. 3A and 3B, respectively. The two panels in Fig. 3 show the 2D distribution of cells and q in the top half and the y-averaged chemoattractant concentration (red line) and cell fraction (purple line) in the bottom half at \(t={7.166}{\hbox {h}}\). The plots also display the range of experimental cell fraction (green-dashed line) observed along the x-axis in the bottom half at the end of the experiment.For \(q_0={0.1}{\upmu \hbox {M}}\), our simulation shows a few cells moving to the right wall (Fig. 3A), as observed in the experiment. These moving (or leading) cells deplete the chemoattractant, which monotonically increases from left to right, with \(q\approx 0\) near the left wall (see red line in Fig. 3A). Thus, cells in that region do not perceive any chemotactic signal and slowly spread out driven by a persistent random walk. The magnitude of the chemotactic gradient sensed by the leading cells is low, but high enough to sustain the motion of the leading cells along the channel (see Fig. 6 in Methods).When the chemoattractant level increases to \(q_0={1}{\upmu \hbox {M}}\) (see Fig. 3B), we observe a significantly larger number of cells moving to the right (compare Figs. 3A and B). The higher chemoattractant concentration enables the presence of larger chemoattractant gradients along longer distances, which results in a more populated wave of leading cells. In this case, the gradients and average levels of chemoattractant sensed by the leading wave are in an optimal range for chemotaxis (see Fig. 6 in Methods). For this reason, the leading wave exhibits robust migration towards the right wall. The chemoattractant concentration behind the leading wave is very low so that cells near the left wall perform a non-chemotactic persistent random walk and slowly spread out (Fig. 3B). On varying the values of \(q_0\) between \({0.1}{\upmu \hbox {M}}\) and \({1}{\upmu \hbox {M}}\), we found that the cells migrate strongly towards the right wall when \(q_0 \ge {0.75}{\upmu \hbox {M}}\).The cell fraction distributions at the final time for both values of \(q_0\) show good agreement with the experiment (compare solid-purple and dashed-green lines in Fig. 3), which confirms the potential of our model to capture collective cell migration of non-confluent systems through self-generated chemoattractant gradients.
Cell migration in narrow channels
Cell migration in the human body is restricted by topographical features such as surrounding cells or components of the extracellular matrix39. Several experiments have been performed in microfluidic assays featuring narrow channels to measure gradients and morphological changes that a cell undergoes while crawling in the narrow channel39,40,41 In these experiments, the cells crawling in the confined spaces block the chemoattractant from going inside the channel. Hence, assuming that chemoattractant can freely permeate the cell membrane can produce results that contradict the biophysics of chemotaxis and appear counter-intuitive. Here, we present one such case where leukocyte cells crawl chemotactically through narrow spaces between endothelial cells. Ref.42 conducted experiments showing that migrating neutrophils completely occlude chemoattractant diffusion, resulting in zero concentration of chemoattractant behind the cell (see Fig. 3(C) in Ref.42). The experiments also showed the formation of a linear chemoattractant gradient as soon as the cell leaves the channel. In this section, we replicate these findings to signify the importance of cell membrane’s impermeability to chemoattractant.Fig. 4Effect of cell membrane’s permeability to chemoattractant in narrow channels. (A) State of cells and chemoattractant distribution at \(t=0\). (B) Location and shape of cells, and distribution of chemoattractant for case (i) in which chemoattractant can permeate cell membrane and case (ii) in which the cell membrane is impermeable to chemoattractant at \(t=20,40,\) and \({60}{\hbox {min}}\). For these results \(D_q = {60}{\upmu \hbox {ms}^{-1}}\). The corresponding simulation video can be found in Mov. S5. (C) Time evolution of average chemoattractant concentration inside the channel defined a \({\overline{q}}_\mathrm{{ic}} = \frac{\int _{\Omega _\mathrm{{ic}}} q\, \varphi \, \textrm{d}\Omega }{\int _{\Omega _\mathrm{{ic}}} \varphi \, \textrm{d}\Omega }\), where \(\Omega _\mathrm{{ic}} = \{ (x,y) \,\,|\,\, 0 \le x \le 230,\,\, 20 \le y \le 40 \}\). We use three values of diffusion constants \(D_q = 20,60\) and \({120}{\upmu \hbox {m}^{2}\,{\hbox {s}}^{-1}}\), and set \(r_q = 0\), \(b_q = 0\), \(k_q = 0\), \({\overline{u}}_A={0.07}{\upmu \hbox {m}\,\hbox {s}^{-1}}\), \(\overline{\Delta T}={20}{\hbox {s}}\), \(\Delta x = \Delta y = {0.33}\upmu \hbox {m}\) and \(\Delta t = 0.001, 0.0005\) and \({0.00025}{\hbox {s}}\) for \(D_q = 20,\,60\) and \({120}{\upmu \!{\hbox {m}}^2\,{\hbox {s}}^{-1}}\) respectively.In this study, we uniformly place 4 cells of radius \({14}\upmu \hbox {m}\) grouped into 2 cluster fields in a highly deformed configuration inside a narrow channel of height \({20}\upmu \hbox {m}\) as shown in Fig. 4A. Initially, a uniform chemoattractant of concentration \(q_{\text {initial}} = {500}\hbox {nM}\) is present outside the channel, whereas the interior contains no chemoattractant. We consider a computational domain of size \(330 \times {60}{\upmu \hbox {m}^2}\) with Dirichlet boundary condition for cluster fields at all boundaries. Here, we examine the effect of chemoattractant diffusion using diffusion constants ranging from \({20}{\upmu \hbox {m}^2\,{\hbox {s}}^{-1}}\) to \({120}{\upmu \hbox {m}^2\,{\hbox {s}}^{-1}}\). The total flux of chemoattractant is set as zero on every boundary except the rightmost boundary where we impose \(\nabla q \cdot {\hat{n}} = 0.1\times (q – q_{\text {initial}})\) to maintain a fixed value of \(q_{\text {initial}}\) at the right wall. To solely focus on diffusion, we set the production, consumption, and decay rate of chemoattractant to zero. For simplicity, we define the walls using stationary phase fields: the upper wall (\(\phi _w^u\)) corresponds to the region where \(y \ge 40\upmu \hbox {m}\) and \(x \le 230\upmu \hbox {m}\), and the lower wall (\(\phi _w^l\)) as the region of intersection of the planes \(y \le 20\upmu \hbox {m}\) and \(x \le 230\upmu \hbox {m}\); see supplementary information for details. We also assume that the walls are rigid by setting \(\frac{\partial \phi _w^i}{\partial t} = 0, \text { for } i = u,l\), so that they do not deform in reaction to the force exerted by the cells. Since the cell-substrate adhesion term is absent in our model, we represent the cell compression against the wall by reducing the repulsion constant between cell and wall \(\textit{g}_\mathrm{{cw}}\) by 10 times in comparison to cell-cell repulsion constant \(\textit{g}_n\).To study the effect of the cell membrane’s impermeability to chemoattractant, we conduct two sets of simulations. In case (i), chemoattractant can freely permeate the cell but not the walls, which is enforced by taking \(\varphi = \varphi _{i} = 1 – \phi _w^l – \phi _w^u\). In case (ii) \(\varphi = \varphi _{ii}= 1 – \sum _{i=1}^{N_c}\mathbb {C}_i – \phi _w^l – \phi _w^u\), which guarantees that the cell membrane and walls are impermeable to chemicals. Fig. 4B shows the position and shape of the cells, and chemoattractant distribution at \(t=20,\,40\) and \({60}{\hbox {m}}\) for \(D_q = {60}{\upmu \hbox {m}^2{\hbox {s}}^{-1}}\). Panel (i) in Fig. 4B shows that the chemoattractant freely permeates the cells and achieves a nearly constant value forming a shallow linear gradient inside the channel. This results in a steady motion of every cell present in the channel in the direction of the chemotactic gradient. From Fig. 4B(i), we observe that every cell achieves a velocity in the direction towards the channel’s exit at the onset of stimulation, thus leaving the channel at a nearly uniform frequency. We see that by \({60}{\hbox {min}}\) in Fig. 4B(i) the last cell reached close to the exit of the channel.In contrast, when the cell membrane is impermeable to chemoattractant we observe that the leading cell completely blocks the chemoattractant from reaching the trailing cells. As a result, the cells at the left-hand side of the channel experience random persistent motion; see Fig. 4B(ii). However, as soon as the leading cell exits the channel, the chemoattractant diffuses inside the channel and gets occluded by the first cell it encounters inside the channel. Comparing the results at different times in Fig. 4B(ii), we can conclude that the leading cell in the channel displaces the chemoattractant out of the channel as it migrates towards the exit. This displacement of chemoattractant by the cells is quite significant in case (ii) as the chemoattractant concentration becomes higher than \(q_\mathrm{{initial}}\) as a result of cell movement. This effect was absent in case (i). These results are in good agreement with the experiments done in ref.42, where the leading cell completely blocks the chemoattractant from reaching trailing cells. As a result, the frequency of cells exiting the channel is non-uniform and monotonically increasing. Mov. S4, S5, and S6 demonstrate the time evolution of cells and chemoattractant when \(D_q = {20}{\upmu \!{\hbox {m}}^2\,\hbox {s}^{-1}},\, {60}{\upmu \!{\hbox {m}}^2\,{\hbox {s}}^{-1}},\) and \({120}{\upmu \!{\hbox {m}}^2\,{\hbox {s}}^{-1}}\), which illustrates that in case (i), cells start migrating towards right direction from \(t=0\). In contrast, in case (ii), cells travel towards the channel’s exit only when no other cell is in front of them blocking the channel.Fig. 4C displays how the average chemoattractant concentration inside the channel denoted by \({\overline{q}}_\mathrm{{ic}}\) evolves over time for cases (i) and (ii). We note that \({\overline{q}}_\mathrm{{ic}}\) achieves a steady state value for each diffusion constant when the chemoattractant can freely permeate the cell membrane as indicated by dashed lines in Fig. 4C. Higher diffusion of chemoattractant leads to faster establishment of a steady average chemoattractant value inside the channel. Meanwhile, we observe that \({\overline{q}}_\mathrm{{ic}}\) oscillates with time for case (ii) (see Fig. 4C solid lines), achieving a minimum value of zero whenever the rightmost cell inside the channel reaches the exit. This oscillation is a consequence of the chemoattractant diffusing inside the channel reaching the rightmost cell, and then getting pushed outside the channel as the cell migrates. We also note in Fig. 4C that the peak value of \({\overline{q}}_\mathrm{{ic}}\) in one cycle (i.e., the time between the exit of two consecutive cells) increases over time. This occurs because the trailing cells, which are devoid of chemoattractant, exhibit random motion and nearly stay at the initial location until the chemoattractant has a free path to reach the cell.Our results align with the experimental findings in42 and demonstrate that neglecting the impermeability of cell membrane to chemoattractant can lead to inaccurate results when the space for cell motion and chemoattractant diffusion is restricted. This highlights the importance of our model for predicting cell migration.
Cell migration in complex chemoattractant environments
Another common scenario in which chemoattractant dynamics are complex and need to be accurately predicted to understand chemotaxis is when there are multiple extracellular species (e.g., an attractant and an enzyme) that interact with each other, giving rise to complex patterns. These patterns emerge from nonlinear dynamics and are sensitive to small perturbations. A cellular system that exhibits this behavior is a conglomerate of vascular mesenchymal cells (VMC), which produce bone morphogenetic protein 2 (BMP-2) and matrix carboxyglutamic acid protein (MGP)43,44. BMP-2 is a molecule recognized for its chemotactic effect on VMC45, while MGP is an enzyme that inhibits the production of BMP-246,47. The different sizes and, therefore, diffusivity of BMP-2 and MGP lead to the formation of stable Turing patterns48, where BMP-2 exhibits spatial localization. Here, we show that accurately predicting pattern formation and chemotaxis requires solving the chemoattractant dynamics in the actual time-evolving extracellular space. In contrast, most models of non-confluent multicellular systems have assumed that the cells are permeable to extracellular components to reduce model complexity by avoiding solving the chemoattractant dynamics on a moving domain. Here, we quantify the impact of this assumption by studying cell pattern formation in two scenarios: (A) The cells are assumed to be permeable to the extracellular components; and (B) The dynamics of extracellular species occur in the time-evolving extracellular environment. These two scenarios can be studied in our computational model by using the extracellular marker \(\varphi\) (see Eq. (3)). In case A, we take \(\varphi =\varphi _\text {A}=1\), while in case B we use \(\varphi =\varphi _\text {B}=1-\sum _{i=1}^{N_c}\mathbb {C}_i\). The interactions between BMP-2 (the chemoattractant or activator) and MGP (the inhibitor) can be modeled with the reaction-diffusion equations43,44,49:$$\begin{aligned} \frac{\partial (\varphi q_a)}{\partial t} = \nabla \cdot (D_{a}\varphi \nabla q_a) + \varphi \eta \bigg (\frac{q_a^2}{q_i(1+k q_a^2)} – r_a q_a + S \bigg ), \end{aligned}$$
(4)
$$\begin{aligned} \frac{\partial (\varphi q_i)}{\partial t} = \nabla \cdot (D_{i}\varphi \nabla q_i) + \varphi \eta \big ( \chi q_a^2 – r_i q_i \big ), \end{aligned}$$
(5)
where \(q_j\), \(D_j\) and \(r_j\), for \(j=\{a,i\}\), are the concentration, diffusion coefficient, and decay rate, respectively, of the activator (a) and inhibitor (i). The parameter \(\eta\) controls the overall kinetic rate, k governs the negative feedback between \(q_a\) and \(q_i\), \(\chi\) controls the inhibitor production, and S is the activator source. Because \(D_i>D_a\) and \(q_i\) decreases the production of \(q_a\) (see k-term in Eq. (4)) the system of equations leads to the formation of patterns in the activator concentration. Eqs. (4) and (5) assume spatially uniform production of \(q_a\) and \(q_i\) in the region defined by \(\varphi\), rather than localizing the production to the cells membrane (as done in44,50). This simplification allows us to study the influence of cell permeability in a simpler scenario.Fig. 5Impact of cell permeability to extracellular components during Turing pattern formation. Cells and activator \(q_a\) distribution for (A) case A (cells permeable to extracellular compounds) and (B) case B (cells not permeable) at \(t={2}{\hbox {h}}\) (left) and \(t={4.5}{\hbox {h}}\) (right). (C) Zoom-in of regions of interest in case A (i, top) and case B (ii, bottom). White lines represent the cell membrane. The colormap shows the \(q_a\) distribution outside and inside the cells. (D) Time evolution of \(\mathbb {R}_{\text {A},\text {B}}\) (defined in the main text). We set \(D_i = {30}{\upmu \hbox {m}^2\,\hbox {s}^{-1}}\), \(D_a = {0.3}{\upmu \hbox {m}^2\,\hbox {s}^{-1}}\), \(r_a = r_i = 1\), \(k ={0.72}{\hbox {nM}^{-2}}\), \(\chi ={1}{\hbox {nM}^{-1}}\), \(S={0.02}{\hbox {nM}}\), and \(\eta = {0.00469}{\hbox {s}^{-1}}\), \({\overline{u}}_A={0.04167}{\upmu \hbox {m}\,{\hbox {s}}^{-1}}\), and \(\overline{\Delta T}={20}{\hbox {s}}\).In this instantiation of our model, Eqs. (4) and (5) replace Eq. (3). The equations governing cell dynamics [Eqs. (1) and (2)] remain identical, but chemotaxis is directed by \(q_a\) instead of q. We study the evolution of a multicellular system for \({4.5}{\hbox {h}}\), which permits us to neglect cell growth, division, and death. For both cases A and B, we consider a computational domain of size \(950\times \,{950}{\upmu \!{\hbox {m}}^2}\) with periodic boundary conditions. We initially place 1040 VMC grouped in 13 clusters. All cells have an initial radius of \({6}{\upmu \hbox {m}}\) and a random initial location. The variables \(q_a\) and \(q_i\) are zero in the intracellular space at the initial time. The initial conditions for \(q_a\) and \(q_i\) in the extracellular space are given by \(q_a = {\overline{q}}_a + p_a\) and \(q_i = {\overline{q}}_i + p_i\). Here, \({\overline{q}}_a = {8.36}{\hbox {nM}}\), \({\overline{q}}_i = {7.81}{\hbox {nM}}\) and \(p_a\), \(p_i\) represent random perturbations sampled from a uniform distribution such that \(q_a\) and \(q_i\) deviate up to 10% from their average values \({\overline{q}}_a\) and \({\overline{q}}_i\).Figures 5A and 5B show the pattern formation process for cases A and B, respectively. The results at \(t={2}{\hbox {h}}\) for case A, which assumes that the chemoattractant can diffuse freely into the cells interior, do not show any discernable pattern (Fig. 5A-left). However, in case B the model prediction at \(t={2}{\hbox {h}}\) already shows a strong pattern for cells and activator (Fig. 5B-left). The faster pattern formation in case B is caused by the motion of the cells. The combined action of cell motion and activator pattern formation leads to a positive feedback effect that accelerates the aggregation. Cells displace the chemoattractant, which accumulates in the cell’s front, and concentrate the chemoattractant in certain regions, leading to a stronger chemical signal. This phenomenon is not observed when cells are permeable to extracellular components (case A). At the final time \((t={4.5}{\hbox {h}})\), both cases display distinctive patterns (Figs. 5A-right and B-right). However, in case A, the cells have a more uniform spatial distribution within the activator pattern. One manifestation of this is that the pairwise distance of cells within the pattern is more uniform in case A than in case B. We hypothesize that this behavior is caused by the higher complexity of the chemoattractant distribution in case B. The cells trap chemoattractant in small regions, leading to a larger number of local maxima in the chemical signal, which produces a more random motion within the activator pattern; compare insets plotted in Fig. 5C.A more quantitative comparison between the two cases can be made by computing$$\begin{aligned} \mathbb {R}_{\text {A},\text {B}}(t)=\frac{ \sum _{\alpha =1}^N \left[ \big | {\textbf {c}}_\alpha ^B(t) – {\textbf {c}}_\alpha ^B(0) \big |^2 – \big | {\varvec{c}}_\alpha ^A(t) – {\varvec{c}}_\alpha ^A(0) \big |^2 \right] }{\sum _{\alpha =1}^N \big | {\varvec{c}}_\alpha ^A(t) – {\varvec{c}}_\alpha ^A(0) \big |^2 }. \end{aligned}$$
(6)
Here, \({\varvec{c}}_\alpha ^A(t)\) is the center of cell \(\alpha\) at time t for case A and \({\varvec{c}}_\alpha ^B(t)\) is the same quantity for case B. Thus, \(\mathbb {R}_{\text {A},\text {B}}\) represents the relative difference of the cell’s mean square displacement between cases A and B. Positive values of \(\mathbb {R}_{\text {A},\text {B}}\) imply larger cell displacements with respect to the initial position in case B, compared to case A. We plot the time evolution of \(\mathbb {R}_{\text {A},\text {B}}\) in Fig. 5D. The plot shows that \(\mathbb {R}_{\text {A},\text {B}}\) increases rapidly at first, reaching a maximum at \(t\approx {2}{\hbox {h}}\). This increase is compatible with the faster pattern formation observed in case B: While the cells in case A still sense an approximately uniform activator distribution at \(t\approx {2}{\hbox {h}}\) that limits their mobility, the cells in case B have undergone significant displacement and formed a well-defined pattern. This effect is also observed in Mov. S7, where cells in case B rapidly develop a pattern, while cells in case A begin to form a pattern only after the chemoattractant establishes the Turing pattern. In the time interval from \(t\approx {2}{\hbox {h}}\) to \(t\approx {200}{\hbox {min}}\), the cells in case A form a pattern, whereas the cells in case B remain within the previously established activator pattern. This explains why \(\mathbb {R}_{\text {A},\text {B}}\) decreases in this time interval; see Fig. 5D. For both cases, after \(t\approx {200}{\hbox {min}}\) the cells display similarly low levels of motion within the activator pattern and therefore \(\mathbb {R}_{\text {A},\text {B}}\) remains roughly constant.Our results show that assuming that the chemoattractant can diffuse into the intracellular space leads to an important underestimation of the ability of the cells to form patterns, especially at early times. This emphasizes the importance of solving the chemoattractant dynamics in the extracellular space and highlights the advantages of the DCF approach to model the chemotaxis of nonconfluent multicellular systems.