Continuum Model for Directed Cell Migration

• The continuum model describes cell distributions as continuous fields governed by PIDEs that combine random motility with non-local guidance cues.
• It simulates collective behaviors such as aggregation, pattern formation, and directed migration on both static and dynamically growing domains.
• The framework integrates density-dependent interactions and appropriate boundary conditions to mimic biological processes like tissue self-organization and morphogenesis.

A continuum model for directed cell migration describes the spatiotemporal evolution of cell distributions as continuous fields, capturing the interplay of random motility, directional cues, and cell–cell or cell–environment interactions at scales where discrete cell identities are not individually resolved. Such models constitute a fundamental theoretical framework for analyzing morphogenesis, tissue self-organization, wound repair, cancer invasion, and organogenesis, enabling robust simulation of collective patterning behaviors, aggregation, and migration on expanding or morphologically evolving domains (Mederacke et al., 10 Sep 2025).

1. Fundamental Structure and Governing Equations

Continuum models for directed cell migration are based on conservation laws for a cell density field $c(\mathbf{x}, t)$ over spatial domain $\Omega$ and are governed by partial integro-differential equations (PIDEs) of the form: $$\frac{\partial c}{\partial t} = -\nabla \cdot \mathbf{j}$$ The total cell flux $\mathbf{j}$ typically consists of:

• Random motility (modeled as isotropic or density-dependent diffusion):

$\mathbf{j}_\mathrm{rand} = - d(c) D \nabla c$

where $D$ is a diffusion coefficient and $d(c)$ allows for nonlinear or pressure-driven motility.

• Directed, non-local migration via density-dependent guidance cues:

$$\mathbf{j}_\mathrm{dir} = a h(c) \mathbf{I}(c,\mathbf{x})$$

with migration strength $a$, weighting function $h(c)$, and a non-local interaction integral

$$\mathbf{I}(c, \mathbf{x}) = \int_0^R \int_{S^{n-1}} g(c(\mathbf{x} + r \boldsymbol{\eta}))\, r^{n-1} \, \boldsymbol{\eta}\, d\boldsymbol{\eta}\, dr$$

where $R$ is the sensing radius, $g(c)$ is a function encoding how local density modulates the guidance, and the integration is over all directions $\boldsymbol{\eta}$ on the unit sphere $S^{n-1}$. This non-local term enables the incorporation of behaviors such as aggregation, contact guidance, chemotaxis, and crowding or saturation effects (Mederacke et al., 10 Sep 2025, Chen et al., 2019).

In dimensionless variables (after scaling with a characteristic length $L_c$), the canonical equation is: $$\frac{\partial c}{\partial \tau} = \nabla \cdot \left( d(c) \nabla c \right) - \alpha \nabla \cdot \left( h(c) \mathbf{I}(c, \xi) \right)$$ where $\alpha = a L_c^{n+1} / D$ is a nondimensional parameter quantifying the ratio of directed to random migration, and $\xi = \mathbf{x} / L_c$.

2. Non-Local Sensing and Density-Dependent Guidance

Non-local terms are essential for modeling realistic cell responses to spatially extended micro-environmental cues. The non-local integral

$$\mathbf{I}(c, \xi) = \int_0^\rho \int_{S^{n-1}} g(c(\xi + r \boldsymbol{\eta}))\, r^{n-1} \, \boldsymbol{\eta}\, d\boldsymbol{\eta}\, dr$$

with $\rho = R/L_c$, mathematically encodes the mechanism by which cells integrate directional information from their surrounding microenvironment within a finite radius. The structure allows for:

• Density-dependent aggregation: e.g., using $h(c) = c(1-c)$ for crowding-limited attraction,
• Pattern formation: cells sort or aggregate based on the averaged densities in their locality,
• Multi-cue integration: the kernel $g(\cdot)$ can be generalized to incorporate contact guidance, chemotaxis, or other signal modalities (Mederacke et al., 10 Sep 2025, Chen et al., 2019, Loy et al., 2020).

In contrast to local models, this framework enables robust simulation of patterning and aggregation phenomena, as the integral operators can capture long-range interactions and spatially extended sensing, which are observed in biological systems.

3. Boundary Conditions and Domain Handling

Two major types of boundary conditions are systematized for simulation:

• Zero-Flux/Neumann Condition: On a physical boundary $\partial\Omega$, the net normal flux vanishes,

$$\mathbf{j} \cdot \mathbf{n} \big|_{\partial\Omega} = 0$$

ensuring conservation of cell mass and precluding artificial inflow or outflow at the edge. This is typically managed at the numerical implementation level by introducing an auxiliary buffer region (with $c=0$), so the non-local integrals remain well-defined even near the boundary (Mederacke et al., 10 Sep 2025).

• Periodic Boundary Conditions: Used for toroidal domains, facilitating studies of pattern formation absent edge effects. Implemented via coordinate mappings that wrap integration points back into the principal simulation domain, ensuring both the diffusion and non-local migration terms are consistent across boundaries (Mederacke et al., 10 Sep 2025).

4. Extension to Dynamic and Growing Domains

Continuum PIDEs can be rigorously reformulated into Lagrangian coordinates to efficiently handle tissue growth and deformation. By introducing a stretch factor $\Theta(\tau) = \frac{L_0}{L(\tau)}$ and material coordinates $\Xi = \Theta(\tau) \xi$, domain growth is accommodated without remeshing, and the time-evolution becomes: $$\frac{\partial c}{\partial \tau} = \Theta^2(\tau) \nabla_L \cdot \left( d(c) \nabla_L c \right) - \Theta(\tau) \alpha \nabla_L \cdot \left( h(c) \mathbf{I}_L(c, \Xi, \tau) \right) - \frac{n v c}{L_0 + v\tau}$$ with a correspondingly rescaled non-local interaction

$$\mathbf{I}_L(c, \Xi, \tau) = \int_0^{\Theta(\tau)\rho} \int_{S^{n-1}} g(c(\Xi + r \boldsymbol{\eta}))\, r^{n-1} \, \boldsymbol{\eta}\, d\boldsymbol{\eta}\, dr$$

where $v$ is the domain growth rate and $n$ the spatial dimension (Mederacke et al., 10 Sep 2025). This construction removes advection terms related to uniform domain growth and allows direct comparison of pattern evolution on dynamically expanding tissues (as in organogenesis or developmental morphogenesis).

5. Biological Applications and Pattern Formation

The continuum PIDE framework models:

• Cell sorting and aggregation: Non-local density-dependent signals drive self-organization of cell populations into clusters, stripes, or labyrinthine patterns (Mederacke et al., 10 Sep 2025).
• Collective flows and tissue-level transitions: Incorporating density-dependent guidance cues (e.g., chemotaxis, contact guidance, adhesion) generates stable patterns and tracks emergent flows, classifying regimes that recapitulate developmental processes and tissue patterning (Chen et al., 2019, Falcó et al., 2022).
• Robustness to geometry and dimensionality: The COMSOL Multiphysics implementation supports one-, two-, and three-dimensional domains, accommodating complex topologies and boundary geometries.
• Interplay with growth: When applied in the Lagrangian frame, the framework predicts how domain expansion and dynamic tissue shaping reinforce or disrupt cellular patterns as organ-scale morphogenetic movements unfold (Mederacke et al., 10 Sep 2025).
• Phase transitions and parameter sensitivity: The model’s dimensional parameters (e.g., $\alpha$, $\rho$, the choice of $h(c)$ and $g(c)$) govern transitions between homogeneous spreading, stationary aggregation, and dynamic patterning, enabling systematic exploration of morphogenetic regimes.

6. Implementation Considerations and Computational Strategies

Simulation platforms such as COMSOL Multiphysics® allow direct implementation of PIDEs with built-in integration operators (e.g., intop, diskint, ballint) to evaluate non-local terms efficiently in arbitrary dimensions. The non-locality is handled in practice by:

• Extending the computational domain with a buffer for zero-flux conditions,
• Using built-in periodic mapping for toroidal geometries,
• Adjusting integration regions dynamically for time-dependent or growing domains in the Lagrangian frame,
• Modifying $d(c)$ and $h(c)$ to incorporate various biological effects, such as pressure-limited or crowding-restricted migration.

Key non-dimensional groups (for instance, $\alpha = a L_c^{n+1}/D$) set the relative strength of the directed, density-dependent (aggregation or repulsion) component to random motility. The choices of functional form for $h(c)$ and $g(c)$ enable modeling a range of collective behavior (e.g., logistic crowding, monotonic attraction, saturation).

7. Significance and Scope

Continuum models articulated as PIDEs provide a generalizable, experimentally relevant theoretical framework for investigating directed collective cell migration in morphogenesis, tissue engineering, cancer, and regenerative biology. Their success rests upon:

• Capturing the non-local, collective, and density-dependent nature of cell guidance.
• Providing predictive insight into conditions that produce aggregation, sharp tissue boundaries, labyrinth instabilities, and the stabilization/destruction of patterns under tissue deformation.
• Allowing the simulation of large cell populations and tissue-scale processes, which would be otherwise intractable for agent-based models at sufficient scale.

By bridging microscopic cellular behaviors with macroscopic tissue dynamics, such models are integral to the quantitative systems biology of cellular organization and pattern formation in complex biological contexts (Mederacke et al., 10 Sep 2025, Chen et al., 2019, Falcó et al., 2022).