Simulating Organogenesis in COMSOL Multiphysics: Tissue Patterning with Directed Cell Migration (2509.08930v1)

Abstract: We present a COMSOL Multiphysics implementation of a continuum model for directed cell migration, a key mechanism underlying tissue self-organization and morphogenesis. The model is formulated as a partial integro-differential equation (PIDE), combining random motility with non-local, density-dependent guidance cues to capture phenomena such as cell sorting and aggregation. Our framework supports simulations in one, two, and three dimensions, with both zero-flux and periodic boundary conditions, and can be reformulated in a Lagrangian setting to efficiently handle tissue growth and domain deformation. We demonstrate that COMSOL Multiphysics enables a flexible and accessible implementation of PIDEs, providing a generalizable platform for studying collective cell behavior and pattern formation in complex biological contexts.

• The paper presents a continuum PIDE model that captures directed cell migration and non-local interactions to simulate tissue patterning.
• Numerical implementation in COMSOL employs built-in integral operators and careful meshing to achieve robust simulations in 1D, 2D, and 3D domains.
• Simulation results demonstrate how variables like migration strength, initial density, and domain growth govern emergent spot and labyrinth patterns.

Simulating Organogenesis in COMSOL Multiphysics: Tissue Patterning with Directed Cell Migration

Model Formulation and Theoretical Framework

The paper presents a continuum model for directed cell migration (DCM) formulated as a partial integro-differential equation (PIDE), capturing the interplay between random motility and non-local, density-dependent guidance cues. The cell density field $c(\vec{x}, t)$ evolves according to

$$\frac{\partial c}{\partial t} = -\vec{\nabla} \cdot \vec{j}$$

where the flux $\vec{j}$ is given by

$$\vec{j} = -d(c) D \vec{\nabla} c + a h(c) \vec{I}(c, \vec{x})$$

with $D$ as the diffusion coefficient, $a$ as the strength of directed motility, and $\vec{I}(c, \vec{x})$ representing non-local sensing within a radius $R$. The model incorporates four essential aspects: random/pressure-driven motility, environmental sensing, cell-cell attraction, and local competition for space (crowding). The non-local term $\vec{I}(c, \vec{x})$ integrates cell densities over a neighborhood, enabling the simulation of phenomena such as cell sorting and aggregation.

Non-dimensionalization is performed for computational efficiency, with characteristic length $L_c$ set to the sensing radius $R$, yielding dimensionless parameters for migration strength ($\alpha$), sensing radius ($\rho$), and time ($\tau$). The model is further extended to account for domain growth via a Lagrangian transformation, allowing efficient simulation of expanding tissues without remeshing.

Numerical Implementation in COMSOL Multiphysics

The PIDE is implemented in COMSOL Multiphysics v6.3 using the Coefficient Form PDE interface. The non-local integral terms are realized using built-in operators: intop in 1D, diskint in 2D, and ballint in 3D, with the number of quadrature points ($m$) controlling integration accuracy. Smoothing parameters ($\delta$, $\varepsilon$) are introduced for numerical stability, particularly to avoid singularities at $r=0$.

Boundary conditions are handled via geometric coordinate mapping for periodic domains and buffer zones for zero-flux (Neumann) boundaries. For periodic boundaries, the domain is extended and mapped using COMSOL's General extrusion operator, ensuring the integral remains well-defined near edges and corners. In 3D, this requires 26 distinct mappings to cover all faces, edges, and corners.

Initial conditions consist of a uniform cell density $c_0$ perturbed by random fluctuations sampled from a uniform distribution, breaking symmetry and enabling pattern formation. The solver employs an implicit BDF scheme with adaptive order and error tolerance, and the mesh is discretized using quadratic Lagrangian elements, with mesh size expressed relative to the sensing radius.

Simulation Results and Scalability

Simulations are demonstrated in 1D, 2D, and 3D domains, with both periodic and zero-flux boundary conditions. The emergent patterns range from spots to labyrinthine structures, modulated by migration strength ($\alpha$) and initial density ($c_0$). The model exhibits sensitivity to parameter choices; for example, $c_0$ near 0.5 and large $\alpha$ or $\rho$ can induce numerical instabilities.

Scalability is validated with large-scale simulations (e.g., $L_0 = 400\rho$), involving millions of mesh elements and degrees of freedom, executed on multi-core CPU architectures. The characteristic pattern length is governed by the sensing radius, allowing biologically relevant tissue sizes to be modeled. Mesh resolution analysis reveals stable spot areas for element sizes up to $0.5\rho$, with negligible differences between triangular and quadrilateral elements. Integration accuracy is contingent on the number of quadrature points; fewer than four points in 2D yield artifacts, while five or more produce robust patterns.

Domain growth is shown to modulate pattern selection: slower growth rates yield fewer, larger spots, while rapid expansion increases spot number and reduces size. This effect is attributed to the rescaling of intrinsic patterning wavelengths during growth, consistent across repeated simulations.

Practical and Theoretical Implications

The presented framework enables flexible, accessible simulation of tissue patterning driven by directed cell migration, supporting arbitrary spatial dimensions and boundary conditions. The use of COMSOL's general-purpose finite element environment circumvents the need for custom solvers, facilitating model extension and integration with other biophysical processes (e.g., reaction-diffusion, tissue mechanics).

From a theoretical perspective, the model provides a tractable approach to studying self-organization and morphogenesis, bridging cellular-scale interactions and tissue-scale patterning. The explicit treatment of non-local interactions and domain growth offers insights into developmental processes, wound healing, and cancer invasion. The framework is amenable to further generalization, including heterotypic cell populations, alternative interaction kernels, and coupling with image-derived geometries.

Future Directions

Potential future developments include:

• Integration of additional biophysical mechanisms (e.g., chemotaxis, durotaxis, topotaxis)
• Extension to multi-species or heterotypic cell interactions
• Systematic parameter estimation and model calibration using experimental data
• Coupling with high-resolution imaging for geometry extraction
• GPU acceleration for large-scale, high-dimensional simulations

The model's generalizability and computational efficiency position it as a valuable tool for quantitative studies of organogenesis and collective cell behavior.

Conclusion

This work establishes a robust, scalable PIDE-based modeling framework for directed cell migration in COMSOL Multiphysics, enabling simulation of tissue patterning and organogenesis across diverse biological scenarios. The approach supports efficient handling of non-local interactions, domain growth, and large-scale computations, with direct applicability to developmental biology, tissue engineering, and pathology. The framework's extensibility and integration capabilities suggest broad utility for future research in computational biology and biophysics.