Comsol Meg Source Localization

Tissue Patterning With Directed Cell Migration

Malte Mederacke1, 2, Chengyou Yu1, Roman Vetter1, 2, and Dagmar Iber1, 2

Switzerland

2Swiss Institute of Bioinformatics, Schanzenstrasse 44, 4056 Basel, Switzerland

Abstract

We present a COMSOL Multiphysics® implementation of a con- tinuum model for directed cell migration, a key mechanism un- derlying 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 generaliz- able platform for studying collective cell behavior and pattern formation in complex biological contexts.

Keywords: Directed cell migration, tissue patterning, self- organization, partial integro-differential equation, COMSOL

Introduction

In this article, we present our implementation of a contin- uum model for Directed Cell Migration (DCM) , a fun- damental biological mechanism by which cells self-organize

In Comsol

Multiphysics®. We have previously made extensive use of COM- SOL Multiphysics® to simulate biological self-organization through reaction–diffusion mechanisms , tissue mechan- ics , and fluid flow . We implemented these models on growing, deforming, or interacting domains using Arbitrary Lagrangian–Eulerian (ALE) or phase-field methods [7, 8], and represented individual cells with compartmentalized diffu- sion and signaling [5, 11]. Moreover, we implemented efficient methods for parameter estimation , and to use image-derived geometries .

DCM enables the formation of spatial patterns from initially homogeneous conditions, relying solely on interactions on the cellular scale, to pattern tissue on the much larger organ scale.

DCM arises from the interplay between random cellular motility and guided movement, ultimately leading to phenomena such as cell sorting and aggregation. In biological systems, guidance cues for directed migration include gradients in chemical signals (chemotaxis), substrate stiffness (durotaxis), surface topography (topotaxis), and differences in cell-cell adhesion (differential adhesion), among others.

Numerical modeling is particularly valuable in studying these mechanisms, given the complexity of cell behavior and tissue geometry, which often defy analytical treatment. Continuum models are well-suited to capture tissue-scale dynamics while maintaining a connection to analytical insights. To this end, single cell properties are locally averaged and the population of cells is represented by a cell density field whose dynamics are driven by cellular flow. In our approach, DCM is described using partial integro-differential equations (PIDEs), which are capable of incorporating both local diffusion and non-local inter- action terms, essential for modeling realistic cell dynamics. The system employed in this work builds on formulations previously developed .

Solving PIDEs poses significant numerical challenges and typically requires custom-built solvers based on finite volume or finite element methods .

These Implementations

often lack flexibility and demand expert-level knowledge to modify or extend. To address this limitation, we implement the model in COMSOL Multiphysics®, which provides the necessary discrete operators to include integral terms in the PDEs. Our approach supports simulations in one, two, and three spatial dimensions and accommodates both zero-flux (Neumann) and periodic boundary conditions, broadening its applicability to various biological scenarios.

Furthermore, tissue growth is an inherent aspect of develop- mental biology and morphogenesis. To handle domain expan- sion efficiently, we demonstrate how the DCM model can be reformulated in a Lagrangian framework, enhancing computa- tional performance while maintaining model fidelity.

In summary, we present a generalizable and accessible numer- ical framework to study a wide range of directed cell migration phenomena, offering new opportunities for exploring pattern formation and organogenesis through simulation.

Theory

In this section, we introduce the governing equations of the PIDE-based continuum model for directed cell migration. The model describes the temporal evolution of a cell density field

(1)

where ®𝑗denotes the cellular flux. This flux can arise from a combination of random or pressure-driven motility and directed, density-dependent migration. For a single-population model, this leads to the following equation in 𝑛spatial dimensions:

∫

𝑆𝑛−1 𝑔𝑐(®𝑥+ 𝑟®𝜂)
𝑟𝑛−1 ®𝜂𝑑®𝜂𝑑𝑟.

(3)

Here, 𝐷denotes the cellular diffusion coefficient and 𝑎the strength of the directed motility forces. The model is designed to capture four key aspects of collective cell behavior. First, cells exhibit motility, which can be random (diffusion-like) or pressure-driven, thereby contributing to tissue fluidity. Second, cells are able to sense their surroundings. This is represented by the integral ®𝐼(𝑐, ®𝑥), which integrates local densities within a sensing radius 𝑅, corresponding to a line in 1D, a disk in 2D, or a ball in 3D. Third, cells interact and attract one another with strength 𝑎within that radius.

For

simplicity, this attraction is assumed to be spatially uniform here, i.e., independent of distance within the sensing radius. Finally, the domain has a finite carrying capacity, such that cells locally compete for space. This crowding effect gives rise to a population pressure that counteracts aggregation. In the model, this population pressure can be incorporated either into the random motility term 𝑑(𝑐), the directed motility term ℎ(𝑐), or the sensing kernel 𝑔(𝑐). In this technical paper, we focus on the special case 𝑑(𝑐) = 1, ℎ(𝑐) = 𝑐(1 −𝑐), and 𝑔(𝑐) = 1. For an exploration of alternative formulations, see .

Non-Dimensionalization

For computational convenience we non-dimensionalize the

(4)

where 𝐿c is a characteristic length scale.

Here, We Choose

𝐿c = 𝑅, i.e., 𝜌= 1. The non-dimensionalized equation then

∫

𝑆𝑛−1 𝑔𝑐( ®𝜉+ 𝑟®𝜂)
𝑟𝑛−1 ®𝜂𝑑®𝜂𝑑𝑟.

Domain Growth

To avoid computationally expensive remeshing during domain growth, we transform the coordinate system from an Eulerian to a Lagrangian frame. This enables the simulation of domain growth on a static geometry. We begin this transformation with the general form of Eqs. 5 and 6, incorporating advection by a velocity field ®𝑢= ®𝜉′(𝜏) induced by the expansion of the tissue:

(7)

To transform Eq. 7 from the Eulerian into the Lagrangian frame,

(8)

where 𝜃(𝜏) = 𝜕Ξ𝑖/𝜕𝜉𝑖= 𝐿0/𝐿(𝜏) is the stretching factor. We consider here uniform linear domain growth of the form 𝐿(𝜏) = 𝐿0 + 𝑣𝜏.

𝐿0

.

(11)

where ®∇L denotes the gradient in Lagrangian coordinates, i.e., with respect to ®Ξ. Since ®𝑢· ®∇𝑐vanishes in the Lagrangian frame, and the divergence of the velocity becomes ®∇· ®𝑢= 𝑛𝐿′(𝜏)/𝐿(𝜏) for uniform growth in 𝑛dimensions, with 𝐿′(𝜏) = 𝑣, Eq. 7

𝑛𝑣𝑐

𝐿0 + 𝑣𝜏.

(12)

The non-local interaction term ®𝐼L in Lagrangian then reads

∫

𝑆𝑛−1 𝑔𝑐(®Ξ + 𝑟®𝜂)
𝑟𝑛−1 ®𝜂𝑑®𝜂𝑑𝑟.

Multiphysics®

All simulations were performed in COMSOL Multiphysics® v6.3. Eq. 5 or 12 was implemented as a Coefficient Form PDE.

Implementation Of The Integral

In COMSOL Multiphysics®, depending on the dimensionality of the problem, different operators are available to implement the integral. In 1D, it is relatively straightforward to use the intop function together with the smoothed Heaviside function flc2hs to perform the directed integration.

The Functions

diskint(r,expr,m) in 2D and ballint(r,expr,m) in 3D compute the spherical integral of an expression expr over a neighborhood of radius r. m specifies the number of quadra- ture points used in the numerical integration. In our case, r corresponds to the sensing radius 𝜌, and expr is the integrand in Eq. 5. It is the availability of these operators that enables COMSOL Multiphysics® to solve PIDEs.

In 1D, the non-local interaction term is implemented as

(14)

with a small smoothing value 𝛿= 0.001 used for numerical

|𝜉−𝜉′| + 𝜀

.

(15)

The small value 𝜀is added to the denominator to avoid division by zero at 𝑟= 0.

𝑟+ 𝜀

.

1)2 + (𝜉2 −𝜉′

2)2.

(18)

Here, (𝜉1, 𝜉2) are the spatial coordinates of the quadrature

2). These Coordinates Can

be accessed using COMSOL’s dest() function. For example,

The Expression 𝜉1 −𝜉′

1 is implemented as x-dest(x). In 3D, the interaction term is similarly defined using the

2)2 + (𝜉3 −𝜉′

3)2.

(21)

For all simulations shown in the following (except where 𝑚 is explicitly varied), we use 𝑚= 10 and 𝜀= 10−5. Figure 1: Illustration of boundary conditions for a 2D do- main. a Geometric coordinate mapping for periodic boundary conditions. b Zero-flux boundary conditions only need a buffer zone around the simulation domain Ω.

Periodic Boundary Conditions

In the above implementation, it is necessary to ensure that the integral remains well-defined near the domain boundary, where the integration region extends beyond the domain. This was achieved with a coordinate mapping. We illustrate the procedure here in 2D (Fig. 1). For a quadratic domain centered at (0, 0), the primary domain (Sq1) of side length 𝐿0 is embedded into a concentric auxiliary quadratic domain (Sq2) of side length 𝐿0 + 2𝜌. For periodic boundary conditions, the positions from the outer square are then mapped to the inner square with the General extrusion operator, using the following rules in 2D:

(23)

(orange domain), one defines a new General extrusion op- erator. In this case, the three upper domains are selected, with

The

x-expression to 𝑥as we only want to shift the values in 𝑦direction, and the y-expression to 𝑦+ 𝐿0, which maps the selected region onto the orange domain labeled 𝑐map, ranging

(+𝐿0/2, −𝐿0/2). (25)

This procedure is then repeated for all four edges and four corners, totaling in eight mappings in 2D. For 3D, the squares are replaced by cubes centered at (0, 0, 0) and the 𝑧-component is added to the mapping in the form of:

(26)

The total number of coordinate mappings is then 26 (6 faces, 12 edges, 8 corners).

Related Journal Articles & DOI Links

Selected peer-reviewed publications relevant to 12 Lead ECG Acquisition. Click the DOI to access the full paper (may require institutional access).

Why Choose Us?

Bangalore guidance for robotics, Spectre and autonomous systems projects.