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

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.

Figures (5)

• Figure 1: Illustration of boundary conditions for a 2D domain.a Geometric coordinate mapping for periodic boundary conditions. b Zero-flux boundary conditions only need a buffer zone around the simulation domain $\Omega$.
• Figure 2: Summary of the patterning space.a Emergent cell density patterns in one-, two-, and three-dimensional domains with periodic boundary conditions. b Variation of the initial concentration $c_0$ and migratory strength $\alpha$ generates a wide range of patterns. Over time, these patterns coarsen via Ostwald ripening, yielding fewer, larger, and more stable structures.
• Figure 3: Exemplary simulation on a large domain. The simulation was performed on a square domain of side length $400\rho$, with initial concentration $c_0 = 0.5$, migratory strength $\alpha = 10$, and integrated until $\tau_\mathrm{end} = 20$. The computational mesh was composed of 1,832,722 elements and 917,970 nodes, resulting in 3,630,085 degrees of freedom. The simulation ran for 16 h and 41 min on 8 cores of two Intel Xeon Gold 6244 CPUs (3.60 GHz).
• Figure 4: Impact of mesh element size and number of integration points.a Example of segmented spots for a mesh with maximum element size of $0.1\rho$. b Scaling of spot area with mesh resolution. c Impact of the number of integration points for the integral on the pattern. For all simulations, zero flux boundary conditions, $\alpha = 10$ and $c_0 = 0.3$ were used. Spots were analyzed at $\tau = 20$.
• Figure 5: Impact of growth speed on patterning. Slower growth (a) leads to fewer, more spaced out spots than faster growth (b) on domains of the same size. Slower growth was simulated with a growth speed of $v = 0.1$ for $\tau_\mathrm{end}=100$, the faster growth with $v = 1.6$ for $\tau_\mathrm{end}=6.25$.