{'fraction_non_alphanumeric': 0.02764699989982971, 'fraction_numerical': 0.007913452869878794, 'mean_word_length': 4.939321832242713, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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.


Introduction Organogenesis is a highly
dynamic process that is tightly regulated during embryogenesis. Many of the individual regulatory components, e.g. signaling molecules and their receptors, as well as their regulatory interactions have been identified in experiments. However, an integrative mechanistic understanding of the regulatory network is missing [1]. Mathematical modeling has a long history in developmental biology [2,3]. Limb development, in particular, has attracted much attention from modellers [4]. Early models were rather simplistic, and to this date most models are still solved on idealized domains that at most qualitatively resemble the physiological domains. However, the geometry can greatly impact the patterning process [5], and it is therefore important to solve these models on physiological domains. COMSOL Multiphysics is a versatile package that provides finite element method (FEM)-based solvers to solve a wide range of partial differential equation (PDE)-based problems on complex domains. We have used COMSOL to solve models of limb development [5,6], bone development [7], ovarian follicle development [8], and branching morphogenesis [9,10,11]. In a series of papers on simulating organogenesis in COMSOL [12,13,14,15], we have discussed methods to efficiently solve models for organogenesis on complex static and growing domains as well as models, which consider cells explicitly. Initially, these models were formulated on idealized geometries. Recently, we have started to take advantage of advancements in imaging techniques, which now provide us with detailed imaging data of organogenesis [15]. This now allows us to simulate our models on realistically growing embryonic domains in COMSOL Multiphysics [16].
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. In the first step, computer readable geometries must be extracted from the 3D images and must then be imported into COMSOL. Many tissues contain clearly defined subdomains with different properties. These can be identified with suitable staining protocols for marker proteins or marker protein expression. We show how complex domains with subdomains can be imported. In a second step, the displacement fields between two consecutive image frames must be calculated and imported into COMSOL. Finally, the imported displacement fields can be used to simulate the domain shape evolution in COMSOL. Given the large number of stages that we use, we implement our models using Matlab LiveLink. We use the ALE method to solve our PDE-based signaling models even on strongly deforming domains. The simulation results can be compared to experimental data, and parameter values can be optimized to obtain an optimal match of model predictions and experimental results [14]. We conclude that the image-based modeling approach allows us to build realistic models of highly dynamic developmental processes, and allows us to study the combined impacts of patterning and growth.

Model Formulation
Our models are defined as a set of n reactiondiffusion equations in the form of: where C i denotes the concentration of component i (n total components), D i its diffusion constant, and Δ refers to the Laplace operator such that D i Δ C i describes the diffusion flux of C i . In case of growing domains, advection and dilution terms have to be added to the reaction-diffusion equations; represents the velocity field of the domain, for more details refer to [15]. The reaction term, R i (C 1 ,…,C n ), is very often nonlinear and describes all reactions of component i, i.e. its production, degradation, and complex formation. For more details refer to [12]. The presence of species can be restricted to parts of the domain, and in that case also some reactions can become spatially restricted.

Regulatory Network
To illustrate our approach, we consider a concrete example. Consider a domain with 3 subdomains as shown in Figure 1, and a regulatory network that involves three components, A, B and C ( Figure 2). All components are assumed to diffuse in the entire domains and to be degraded everywhere. Moreover, we assume that the components A and C are produced only in domain 1 and domain 3 , respectively, whereas component B is produced in all domains. We formulate the sub domains using the unit function, The spatio-temporal dynamics of the aforesaid network can be described by the following reaction terms:

Boundaries and Displacement Fields
Using standard techniques for image segmentation, external and internal boundaries can be extracted [16]. This process can be repeated at different developmental time points to obtain a developmental sequence of shapes [15]. In this study, we consider the two geometries in Figure 1 as our extracted 2D geometries at two subsequent developmental time steps, t and t+1. As can be seen, the entire domain as well as the subdomains deform from t to t+1.
To describe the growing domains, we need to calculate the displacement fields between the two shapes at t and t+1. A range of algorithms can be employed, which have their advantages and disadvantages dependent on the details of the geometries and their deformations (Schwaninger et al., submitted). Here, we use the uniform displacement field algorithm proposed by (Schwaninger et al., submitted): consider a curve at time t, γ t , that is deformed to γ t+1 within the next time step. This algorithm interpolates N points on both curves: is equal for all i, j and T ∈ {t, t +1} . The displacement filed matrix, D, is defined as  Figure 3 shows the displacement fields between two curves at two subsequent time steps t and t+1. The displacement fields are imported into COMSOL as Interpolation function and are later employed in the Moving Mesh (ale) module to describe the domain deformation due to the growth.

Displacement of intersecting Curves
The introduction of subdomains results in intersecting boundary curves (Figure 1). The function Interpolation Curve, that we used to generate the boundaries, does not discriminate between intersection points and other points on the curve. All points are interpolated in the same way. Given the interpolation, there is no guarantee that the intersection point of two curves at time t will be accurately displaced to their prescribed intersection point at time t+1. This issue can cause distorted meshes, inverted meshes, and numerical problems close to the intersection points. In case of spatially restricted variables, this inaccuracy can result in leakage of variables out of their restricted domains. From here on we will refer to the model with these problems as Model 1 .
To deal with this problem, we propose the

{ }
We implemented this algorithm for all intersecting curves ( Figure 4C) and determined the displacement fields for each segment individually. The coordinates of each segment and their corresponding displacement fields were then imported into COMSOL separately. Using the above algorithm, we obtain an accurate mapping of all domains and intersection points. From here on, we will refer to this model as Model 2 .

Results
Model 1 and Model 2 were implemented in COMSOL with the parameter values as given in  Table 1 in the Appendix. Figure 5A-C shows that in Model 1 the intersection points are not displaced correctly. This problem leads to the inverted meshes close to the intersection points ( Figure 5B,C bottom panels), and the subdomains are distorted ( Figure 5B,C top panels) as compared to Figure 1B. Segmenting the curves at the intersection points before deformation (Model 2 ) leads to a higher quality of mesh elements and therefore accurate numerical solution and correct domain deformation ( Figure  5D-F). The inverted meshes that result from inaccuracies in the mappings also lead to differences in the domain areas. Thus, Figure 6 reports the area of domain 1 and domain 2 at different time steps. As can be seen, the final area in Model 1 differs from the real size of the domain at t+1. The inverted meshes also lead to numerical problems and consequently inaccurate solutions. Thus, Figure 7 shows that the expression patterns, i.e. the effective spatial production rates, of A and C in Model 1 and Model 2 differ at time t+1. The components A and C have lower expression in Model 1 than in Model 2 .

Conclusion
In this paper, a framework is presented to simulate PDE models on growing, embryonic domains with subdomains, using COMSOL Multiphysics. Intersecting boundaries are not mapped accurately using the standard COMSOL Interpolation function. We addressed this problem by introducing segmented boundaries.
Using this algorithm and COMSOL's Matlab LiveLink interface, one can implement also complicated domains and large sets of PDEs. This permits the simulation of large, complex regulatory networks on physiological domains. We expect that this will further increase the predictive value of the models, and will allow to better test, improve, and validate the models with experimental data.