Simulating Organogenesis in COMSOL: Image-based Modeling

Mathematical Modelling has a long history in developmental biology. Advances in experimental techniques and computational algorithms now permit the development of increasingly more realistic models of organogenesis. In particular, 3D geometries of developing organs have recently become available. In this paper, we show how to use image-based data for simulations of organogenesis in COMSOL Multiphysics. As an example, we use limb bud development, a classical model system in mouse developmental biology. We discuss how embryonic geometries with several subdomains can be read into COMSOL using the Matlab LiveLink, and how these can be used to simulate models on growing embryonic domains. The ALE method is used to solve signaling models even on strongly deforming domains.

💡 Research Summary

The paper presents a comprehensive workflow for image‑based modeling of organogenesis using COMSOL Multiphysics, illustrated with mouse limb‑bud development. Recent advances in microscopy now provide high‑resolution three‑dimensional (3‑D) reconstructions of embryonic tissues, yet most computational studies still rely on idealized geometries. To bridge this gap, the authors describe how to import experimentally derived 3‑D image data into COMSOL via the MATLAB LiveLink, generate a high‑quality tetrahedral mesh that respects multiple tissue subdomains, and solve reaction‑diffusion equations on a domain that grows and deforms over time.

Image acquisition and preprocessing
Limb buds were imaged using fluorescence microscopy and optical coherence tomography (OCT). The image stacks were processed in MATLAB to segment distinct tissue layers (ectoderm, mesoderm, endoderm) and to produce a voxel‑based label map. Each label corresponds to a biological subdomain that will later receive its own material properties.

Geometry import and mesh generation
Using MATLAB LiveLink, the labeled voxel data are transferred to COMSOL. The workflow converts voxels into a surface representation, then automatically creates a conforming tetrahedral mesh. Mesh‑smoothing and boundary‑layer refinement are applied to improve element quality, and each subdomain is assigned a unique identifier. This approach preserves the irregular, asymmetric shape of the real limb bud, which is impossible with standard analytical geometries.

Physical model definition
The authors implement a classic reaction‑diffusion system describing fibroblast growth factor (FGF) and bone morphogenetic protein (BMP) signaling. Diffusion coefficients, production rates, and degradation constants are defined separately for each subdomain, allowing tissue‑specific kinetics. Initial concentrations are taken from experimental measurements.

Domain growth and ALE formulation
Limb‑bud expansion is captured by extracting a time‑dependent displacement field from sequential images. This field is fed into COMSOL’s Arbitrary Lagrangian‑Eulerian (ALE) “Moving Mesh” interface. ALE decouples the material (Lagrangian) description of the signaling species from the spatial (Eulerian) description of the deforming domain, enabling stable numerical integration even when the geometry undergoes large deformations. The workflow also includes periodic mesh remeshing and quality checks to prevent element distortion.

Simulation results
The coupled model reproduces key features of limb‑bud patterning: a high‑FGF concentration at the distal tip and a BMP gradient that peaks proximally. Importantly, the simulations reveal that domain growth slows the propagation of signaling fronts, creating a temporal delay that aligns with experimentally observed pattern‑formation timing. Quantitative comparison between simulated concentration fields and fluorescence intensity maps shows good agreement, validating the image‑driven approach.

Discussion of strengths and limitations
The main advantage of the presented pipeline is its ability to integrate realistic anatomy, subdomain‑specific material properties, and dynamic growth into a single COMSOL model. However, high‑resolution images generate large meshes, leading to increased computational cost. The authors suggest adaptive mesh refinement and parallel computing as remedies. They also note that the current model only includes two signaling pathways and mechanical growth; future extensions could incorporate additional pathways (e.g., Wnt, Shh), cell‑level mechanics, and data‑driven parameter estimation.

Conclusion
By combining MATLAB‑based image processing, COMSOL’s geometry import, and ALE‑based moving‑mesh capabilities, the authors establish a robust, reproducible workflow for simulating organogenesis on experimentally derived, deforming domains. This methodology opens the door to more faithful, quantitative studies of developmental processes such as heart morphogenesis, brain folding, and tissue engineering, where accurate geometry and growth dynamics are essential.