Combining Turing and 3D vertex models reproduces autonomous multicellular morphogenesis with undulation, tubulation, and branching

This study demonstrates computational simulations of multicellular deformation coupled with chemical patterning in the three-dimensional (3D) space. To address these aspects, we proposes a novel mathematical model, where a reaction–diffusion system is discretely expressed at a single cell level and combined with a 3D vertex model. To investigate complex phenomena emerging from the coupling of patterning and deformation, as an example, we employed an activator–inhibitor system and converted the activator concentration of individual cells into their growth rate. Despite the simplicity of the model, by growing a monolayer cell vesicle, the coupling system provided rich morphological dynamics such as undulation, tubulation, and branching. Interestingly, the morphological variety depends on the difference in time scales between patterning and deformation, and can be partially understood by the intrinsic hysteresis in the activator-inhibitor system with domain growth. Importantly, the model can be applied to 3D multicellular dynamics that couple the reaction–diffusion patterning with various cell behaviors, such as deformation, rearrangement, division, apoptosis, differentiation, and proliferation. Thus, the results demonstrate the significant advantage of the proposed model as well as the biophysical importance of exploring spatiotemporal dynamics of the coupling phenomena of patterning and deformation in 3D space.

state can trigger further molecular signaling 14 . Local cell dynamics can therefore be coupled with global tissue dynamics, forming a basis of bidirectional interaction between patterning and deformation at a single cell level.
Mathematical models have been well used for understanding multicellular dynamics [15][16][17][18][19][20][21][22] and have been improved to analyze their 3D dynamics [23][24][25][26][27][28] . We have developed a full 3D vertex model that expresses 3D multicellular dynamics compacted in a monolayer sheet as well as a multilayer aggregate, involving cell rearrangements 29 , division 30 , apoptosis 31 , and viscoelastic behaviors 32 . The models have succeeded in reproducing basic epithelial deformations 33,34 as well as reproducing several developmental phenomena, such as blastocyst formation 35 . Notably, although the intercellular transport of signaling molecules has been expressed in a 3D vertex model 36 , it has not yet been applied to complex patterning caused by reaction-diffusion dynamics. Therefore, combining the Turing and 3D vertex models will aid in the exploration of mechanochemical coupling in multicellular morphogenesis.
In this study, we propose a novel mathematical model that combines the Turing and 3D vertex models, and demonstrate computational simulations of complex phenomena emerging from the coupling of patterning and deformation, in 3D space. In embryogenesis, diffusive molecules can be transduced to various cell behaviors such as deformation, rearrangement, division, apoptosis, differentiation, and proliferation. As an example, an activator-inhibitor system is assumed as a regulatory process of cell proliferation, and local activator concentration is converted into the growth rate of individual cells. By examining the physical parameters of molecular transport coefficients, production and degradation rates, and cell growth rate, we discuss bidirectional effects occurring between patterning and deformation.

Model Framework of Combining Turing and 3D Vertex Models
To analyze 3D multicellular dynamics coupling chemical patterning with mechanical deformation, we develop a mathematical model that combines the Turing and 3D vertex models (Fig. 1a). The Turing model is well known to generate various chemical patterns observed in biological phenomena (Fig. 1b), while the 3D vertex model is a general tool to express mechanical behaviors of 3D multicellular dynamics (Fig. 1c). In the combined model, chemical states of individual cells are regulated by chemical interactions among cells; individual cells generate mechanical forces to deform the tissue according to these chemical states. Simultaneously, the chemical pattern can be dynamically rearranged on the deforming tissue, so as to rewrite the chemical states of individual cells. The combined model therefore enables expression of mechanochemical coupling of patterning and deformation. Details of the combined model framework are described below.

3D vertex model for describing multicellular deformation.
In the 3D vertex model, an individual cell shape is represented by a polyhedron, a boundary face between neighboring cells expressed by a polygon, and the entire structure of a 3D cell aggregate expressed by a single network (Fig. 1c). Topology of the network is dynamically rearranged to express the changes in the cell configuration and the number of cells 29 . According to the force balance, the movement of the ith vertex, represented by r i , is expressed by the following equation: The left hand side of Eq. (1) indicates a viscous friction force exerted on the ith vertex, where scalar η i and V fi are a viscous friction coefficient of the ith vertex and vector and a local velocity field around the ith vertex, respectively. Here, η i and V fi are defined as functions of velocities of vertices to express the viscous property 32 . The right hand side of Eq. (1) indicates an energetic force acting on the ith vertex, where U is an effective energy. By supposing U as a function of r i and molecular concentration, the mechanical deformation of cells depends on their chemical patterning. Discrete Turing model for describing multicellular patterning. Although the Turing model is often described in a continuum manner, in biological systems, molecules diffuse within a multicellular structure. For simplification, as an example, we assume that (i) the molecules transport within cells through junctional structures but not the outside of cells, (ii) the molecular diffusivity within individual cells is much faster than those between neighboring cells, and thereby the molecular concentrations within individual cells can be regarded to be homogeneous, (iii) the gap width between neighboring cells is homogeneous and the molecular distribution between neighboring cells is linear in the normal direction to the cell-cell boundary. Based on these assumptions, the molecular distribution is discretized by individual cell compartments, and the molecular flux between neighboring cells can be simply expressed to be proportional to the difference in the molecular concentrations. This discrete description of the Turing model in the single cell level can be directly combined with the 3D vertex model.
According to the mass conservation, the dynamics of the number of the th molecules at the jth cell, represented by  m j , is expressed by the following equation: The first term of the right hand side of Eq. (2) is the flux from the cells neighboring the jth cell, where ∑ k j ( ) cell is the summation for all neighboring cells. Constant μ  is the transport rate of the th molecules between neighboring cells as the diffusion coefficients divided by the gap width. Scalar s jk cc is the boundary area between the jth and kth cells, and scalar v j is the jth cell volume. Both s jk cc and v j are functions of the location vector of vertices. The second term of the right hand side of Eq. (2) is the reaction term at the jth cell, where function R  is a reaction rate of the th  molecules as a function of a set of molecular concentrations at the jth cell, represented by Because the flux in the first term depends on the network topology and vertex locat i ons r i , the chemical patterning of cells depends on their mechanical deformation.
Notably, while Eq. (2) is similar to those in our previous study 36 , the second term in Eq. (2) is newly defined as a function of a set of molecular concentrations in a discrete manner to combine the Turing and 3D vertex models. The relationship between discrete and continuum Turing models is described in Sect. 6 in Appendix.
In this model, the mechanical deformation and the chemical patterning of multiple cells are thereby coupled bidirectionally. The total degrees of freedom of a tissue are represented by the network topology, vertex locations, r i , and the number of molecules, m j  .

Modeling Example: Cell Growth Regulation by an Activator-inhibitor System
To investigate what can emerge from the coupling of multicellular patterning and deformation, we simulated the tissue growth dynamics regulated by an activator-inhibitor system (Fig. 2). Previous studies have shown that the activator-inhibitor process operates growth regions in the 3D tissue structure 13 . For this simulation, we assumed that tissue deformations are simply driven by the mechanical forces generated by cell proliferation (volume growth and division), whose rate has a positive dependence on the intracellular concentration of activator molecules, represented by . Hence, the activator molecules act as a mitogen in this system. Moreover, the activator molecules have a negative feedback through inhibitor molecules, represented by  . Details of the mathematical assumption are described below.

Viscoelastic behaviors of cells. The elastic properties of cells are expressed by the effective energy U in
Eq. (1) as follows: For simplification, we assume that the viscous property is homogeneous in the entire tissue by ignoring the inhomogeneity of tissue structures such as apical actin belts and basement membrane. Hence, viscous property of cells is simply expressed using functions of η i and V fi in Eq. (1) as employed in our previous study 32 . In this model, friction η i is defined as the summation for the all surrounding cells: Activator-inhibitor system. As an example of the Turing model, we employed the activator-inhibitor system suggested by Gierer and Meinhardt 37 . To focus on the coupling of patterning and deformation, we simplify the equations of activator and inhibitor dynamics in a discrete manner. Dynamics of the numbers of activator and inhibitor molecules within the jth cell, represented by m j  and m j  , are expressed as follows: The first terms in the right hand side of Eqs (5) and (6) are flux, where constants μ  and  μ are the transport coefficients of activator and inhibitor molecules, respectively. The second and third terms in the right hand side of Eqs (5) and (6) are production and degradation terms, where constant γ is the reaction rate of the activator and inhibitor molecules and constant ρ u is the reference concentration. The discrete formulation of Eqs (5) and (6) corresponds to the continuum formulation described in Sect. 7 in Appendix.

Mechanochemical coupling.
To couple multicellular deformation with patterning, we introduce cell proliferation depending on the activator concentration. During cell proliferation, the equilibrium volume v eqj increases depending on the concentration of activator molecules,  m v / j j . Moreover, when v eqj increases up to v th , the jth cell divides into two daughter cells with In biological systems, the cell growth rate is regulated by the concentration of mitogen. Hence, by regarding the activator as a mitogen, the growth rate of the jth cell volume, represented by , is defined as Cell division is represented by dividing a single polyhedron at a plane; the direction of the dividing plane is regulated to be normal to the longest axis of the individual cell shape, on the plane of the tissue surface. Details of the cell division manner are similar to those of the local regulation used in our previous studies 30,38 . In the division process, morphogen concentrations in daughter cells are set to be the same values with those of the mother cell.

Physical parameter setting.
To focus on the coupling of patterning and deformation, we simplified Eqs (5) and (6) by introducing two parameters χ and φ. Constants χ and φ express spatial and temporal characteristics, respectively, and were defined as functions of  μ , μ  , and γ as χ = μ μ / A I and φ =  μ γ / . Using χ and φ, Eqs (5) and (6) were rewritten as Under the steady state, the terms in the left sides of Eqs (9) and (10) become zero. Therefore, steady spatial patterns can be determined by χ independent on γ, and the time scale of the patterning can be regulated by γ while maintaining the same steady pattern. Based on the linear approximation, the length scale of steady activator patterns is proportional to χ 1/4 φ 1/2 .
To solve Eqs (1, 9 and 10), parameter values were normalized by unit length (v ref ) 1 Physical parameters are summarized in Table 1. By assuming an incompressibility, the cell volume elasticity k v was set to be much larger than the characteristic surface energy of a cell, represented by κ s (v ref ) 2/3 . By assuming a quasi-static process, the cell cycle τ cycle was set to be much larger than the characteristic time of cell deformation, represented by η c /κ s . Moreover, to focus on the coupling phenomena, the physical parameters of the patterning and coupling were fixed; because the steady-state solutions of Eqs (9) and (10) u ref u ref , the switching concentration for cell growth was simply set to be the medial value of the steady state solutions The Hill coefficient α was set to be much larger than unit to approximate Eq. (8) as a step-wise function. The diffusivity of inhibitor that determines the distance between activator domains was simply set to be the square of the length scale of 3 cells (φ ≈ (3(v ref ) 1/3 ) 2 /(4η c /5κ s )). Hence, the patterning and deformation behaviors in this model are determined by two physical parameters γ and χ. Numerical implementation is described in Sect. 8 in Appendix.

Results
Turing patterns in 3D cell aggregates. First, to test whether the model can create the periodic pattern of the activator-inhibitor system in 3D cell aggregates, we simulated patterning processes while without cell volume growth ( =  v 0 j eq in Eq. (7)). The arrested tissue morphologies were simply set to be a monolayer cell sheet and a compacted cell aggregate. These tissues are composed of about 200 cells, and the initial molecular concentrations were randomly set to satisfy the mean values m v / j j  = 0.5 and m v / j j  = 0.1. Physical parameters were simply set to be χ = 0.1 in case of the monolayer cell sheet and χ = 0.05 in case of the compacted cell aggregate. Figure 3 shows the steady patterns generated in the 3D cell aggregates, which were driven by the activatorinhibitor system. In case of the monolayer cell sheet, in-plane spot patterns were generated on the sheet (Fig. 3a). In case of the compacted cell aggregate, 3D spot patterns were generated inside the aggregate (Fig. 3b). While the patterning in 2D space has been well studied, the research of the patterning in 3D space has been just started 7,8 . To focus on the coupling phenomena, thereby, we used a monolayer cell sheet as initial tissue morphologies for the analyses.
Patterning hysteresis on deforming tissues. Second, to analyze directional effects from deformation to patterning, we simulated patterning processes on the tissues that grow independently on the molecular concentration. To express cell growth that is independent on the molecular concentration, cell growth rate was set to be constant as λ j = λ ref j in Eq. (8). Moreover, because the 2D spatial patterns driven by the Turing model has been well studied, constant χ that determines the spatial scale of the patterning was fixed as χ = 0.1 and constant γ that determines the time scale of the patterning was varied. Here, the initial tissue morphology was simply set to be a spherical vesicle of a monolayer cell sheet composed of about 200 cells, and the initial tissue patterns were randomly set to satisfy the mean values m v / j j  = 0.5 and m v / j j  = 0.1. As described in Eqs (9) and (10), steady spatial patterns are independent on γ. Consistently, when γ = 1 and 100 (the rate of patterning is faster than or equal to that of deformation), steady patterns were formed as shown in Fig. 4a. In this case, the activator domain size is almost independent on γ, and is almost maintained during tissue growth. Interestingly, the number of activator spots was maintained in the case of γ = 1, whereas it increased in the case of γ = 100. Moreover, when γ = 0.01 (the rate of patterning is much slower than that of deformation), steady patterns were not formed, and activator concentration was gradually blurred (Fig. 4a). This was most likely caused by the dilution effect resulting from cell growth.
To clarify whether the patterns at t = 4.6 cell cycles in Fig. 4a reached steady states, we simulated patterning processes on the tissues with cell growth upto t = 4.6 (as in Fig. 4a) and without cell volume growth after t = 4.6. Figure 4b shows patterns of the tissues at t = 6.6 cell cycles. In the case where γ = 0.01, activator concentrations increased in the whole of the tissue. In the case where γ = 1 and 100, the numbers of activator spots were similar to those at t = 4.6 cell cycles. These results suggest that the patterning shows the dependence on its history, i.e., hysteresis, which becomes weaker in case with higher γ.
Moreover, we also randomize the molecular distribution on the grown tissue at t = 4.6 cell cycles in Fig. 4a, and simulated a repatterning process while arresting cell volume growth, as shown in Fig. 4c. Here, the tissue pattern was randomly set to satisfy the mean values m v / j j  = 0.5 and m v / j j  = 0.1. As a result, a new pattern emerged on the tissue, in which the number of spots is greater than that generated on the growing tissues in Fig. 4b. These results clarified that the hysteresis becomes weaker in case with higher γ. Interestingly, even when increasing γ greatly, the number of spots on the growing tissue (γ = 0.01 in Fig. 4a) is still smaller than that on the arrested tissue (Fig. 4c). Hence, it seems difficult to completely suppress the hysteresis of patterning by regulating γ.
The dependence of the patterning on γ can be understood by the ratio of time scales of patterning to deformation (Fig. 4d). Because the rate of patterning is limited by the reaction rate, it can be affected by γ. On the other hand, because the deformation rate is limited by the rate of cell growth, it can be affected by 1/τ cycle . Where the patterning rate γ is similar to that of deformation 1/τ cycle , patterning shows high hysteresis in addition to a dilution effect in that initial molecular concentrations are strongly maintained and gradually blurred by the dilution effect of cell volume growth (γ = 0.01 in Fig. 4a,d). Where the patterning rate γ is slightly larger than that of deformation 1/τ cycle , patterning shows middle hysteresis in that it is rapid enough to maintain the existing domain areas and slow enough to maintain their positions on tissue (γ = 1 in Fig. 4a,d). Where the patterning rate γ is much larger than that of deformation 1/τ cycle , patterning shows low hysteresis; it is rapid enough to change the existing domain areas, but is too fast to keep their positions on tissue (γ = 100 in Fig. 4a,d). These results highlight the importance of the difference in time scales between patterning and deformation.
Coupling patterning and deformation drives undulation, tubulation, and branching. Third, we investigate the coupling phenomena of patterning and deformation. The initial tissue morphology was simply set to be a spherical vesicle of a monolayer cell sheet composed of about 2,000 cells, whose patterns reached steady states. The steady patterns depend on χ; the model provided various tissue deformations, and several typical deformations such as undulation, tubulation, and branching were shown. Figure 5 shows tubulation deformation processes under the conditions with χ = 0.01 and χ = 0.1, where several bumps sprouted from steady domains of activator regions on the vesicle, and grew into tubes. During tubulation, the activator molecules stayed around the tip of tubes, from which the tubes continuously grew.
As described in Sect. 3.4, the size of activator regions should be approximately proportional to χ 1/4 φ 1/2 ; hence, the size in the case of χ = 0.1 is expected to be about 1.8 times larger than that in the case of χ = 0.01. According to the expectation, the tube diameter varied between thin (χ = 0.01 in Fig. 5a) and thick (χ = 0.1 in Fig. 5b). Moreover, because the size of the activator spots is maintained, the tube diameters are constant under the individual conditions. Since the bending rigidity of tubes depends on their diameters, increasing χ leads to straightening the tubes; while thin tubes became undulated (χ = 0.01 in Fig. 5a), thick tubes became straight (χ = 0.1 in Fig. 5b). Figure 6 shows a branching deformation process, in which several spherical protrusions sprouted from steady domains of activator regions on the vesicle (Fig. 6a), and separated into multiple bumps (Fig. 6b). Interestingly, the molecular concentration was distortedly randomized in the sprouted region as shown in Fig. 6a. This is because of the discrete description of multicellular structure; the activator concentration tends to relax into equilibrium within individual cell compartments before diffusing to the surrounding cells. This non-cooperative behavior of cells tend to form discontinuous molecular distributions in space, which facilitate the domain branching as a perturbation. Therefore, the mechanism of this branching deformation seems partially different from those expected from the domain branching in the continuum space, where molecular concentrations tend to be continuously distributed in space. Figure 7 shows an undulation deformation process, where activator domains dynamically move on the tissue. The pattern can be dynamically remodeled by the disturbance of cell growth. In Figs 5 and 6, the pattern remodeling is much slower than the tissue growth. On the other hand, in Fig. 7, the pattern remodeling is much faster than the tissue growth, and the activator domains rapidly move before forming some protrusions. These results suggest that the coupling of patterning and deformation creates a greater variation of dynamics than expected.
Difference in time scales between patterning and deformation varies a 3D morphology of growing tissues. Fourthly, to understand the general behaviors of tissue patterning and deformation, we calculated the morphology diagram with respect to two physical parameters: χ and γ (Fig. 8). The initial tissue morphology was simply set to be a spherical vesicle of a monolayer cell sheet composed of about 2,000 cells, whose patterns reached steady states. The steady patterns depend on χ as shown in Fig. 8a; with increasing χ, the size of spots increased and the number of spots decreased. Figure 8b shows a morphology diagram of tissue patterns and deformations with respect to χ and γ. Importantly, by comparing the tissues which have the same value of χ, the tissue patterns in Fig. 8b were partially different from the arrested tissues in Fig. 8a, and the tissue morphologies in Fig. 8b varied by γ. For example, the molecular distribution within activator domains where γ = 0.01 in Fig. 8b, was inhomogeneous and more blurred than those of the arrested tissues in Fig.8a. In the case of χ = 0.1, while the tissues where γ = 0.01 and 1 formed thick tube structures, the tissue in case with γ = 100 was undulated. In the case of χ = 0.01, the tissue where γ = 0.01 formed spherical protrusions, while the tissues where γ = 1 and 100 formed thin tube structures.
The dependence of the tissue morphologies on γ can be understood by the hysteresis in patterning. Where patterning shows high hysteresis (γ = 0.01), thereby there tend to remain a disturbance in molecular distribution by cell growth and division at a single cell level, which repatterns activator domains to be separated to form branch structures (Fig. 6). Where patterning shows middle hysteresis (γ = 1) in that it maintains the existing domain areas and their positions on tissue; thus, the middle hysteresis causes the tube structures (Fig. 5). Where patterning shows low hysteresis (γ = 100); it dynamically changes positions of domain areas on tissue; thus, the low hysteresis causes the undulation structures (Fig. 7). These results suggest that the tissue morphologies emerging from the coupling phenomena depend on both spatial and temporal properties of patterning.

Hysteresis in patterning emerges from the multiple unstable mode in activator-inhibitor system.
As shown in Figs 4-8, the computational simulations combining the Turing and 3D vertex models revealed that the patterning hysteresis partially plays an important role in tissue formations emerging from the coupling phenomena. The patterning hysteresis has been implicitly reported in the Turing pattern resulting from a simple one-dimensional model with growing domain 39 ; under the quasi-static process of domain growth, the wavelength of the pattern can change without increasing the number of spots. Lastly, to understand the patterning hysteresis more clearly, we performed the linear stability analysis of the continuum type of the activator-inhibitor system, described in Sect. 7 in Appendix, in one-dimensional space (0 ≤ x ≤ 2π). The analysis revealed that this model has multiple unstable modes in dispersion relation as shown in Fig. 9a; and this system had the most unstable wavenumber, represented by k max (around 5.5 in this case). Hence, we expect that the wavenumber, represented by k, varies around k max .
To clarify the variation of the wavenumber, we performed the numerical simulation of the continuum type of the activator-inhibitor system, described in Sect. 7 in Appendix, in one-dimensional space (0 ≤ x ≤ 2π). As expected above, when the initial condition of  ρ and  ρ was set as a random distribution with the uniform white noise, the resulting k under the steady state varied from 5.0 to 6.5 depending on the initial condition as shown in Fig. 9b. More interestingly, when the initial condition was set as a periodic distribution with k = 4.0, the resulting k became 4.0 as shown in Fig. 9c, which is different from those resulting from the initial condition with a random distribution in Fig. 9b.
Moreover, to understand the effect of domain growth on patterning, we additionally took into account the domain growth in the above simulation. We set L to be 2π under the initial condition with k = 4.0. Then, we calculated the molecular distribution after expanding L to be 4π slowly enough to satisfy the quasi-static process. As a result, even when k became 2.0 that is out of the unstable region of the dispersion relation in Fig. 9a, the initial number of waves was robustly maintained as shown in Fig. 9d. Thus, the patterning can vary depending on the initial condition and show hysteresis in the process of domain growth.

Discussion
Chemical patterning plays a role in regionalizing a tissue, in which a higher-order structure is constructed by integrating different functional cells. The patterning has been simply explained by the Turing model 4 , and observed in biological systems 5 . Despite the difficulty in observing molecular diffusions in vivo, the Turing model has had a significant impact on the understanding of embryonic morphogenesis. Numerical studies over the decades have suggested various kinds of patterns in 2D space 37 , and have led to discovery of further complex patterns in 3D space 7,8 . Currently, descriptions using integral kernel equations have been suggested to explain biological patterns more generally 40 . Moreover, recent studies using continuum models have emphasized the importance of the coupling between patterning and deformation in 3D space 9 . Based on this knowledge, this study has demonstrated computational simulations of 3D multicellular deformation coupled with patterning in the single cell level, which may help connect these basic findings to the understanding of biological systems.
In this study, by combining the Turing and 3D vertex models, a reaction-diffusion process was expressed simply by an activator-inhibitor system coupled with cell growth. The simple coupling system caused various large tissue deformations, which highlighted the importance of the difference in time scales between patterning and deformation. Importantly, branch structures were caused by large tissue deformations, and the inhomogeneous molecular distribution within activator domains was caused by the dilution effect on discrete cells. Therefore, the findings in this study were due to the combined model which describes mechanical and chemical behaviors in the single cell level.
Because of the 3D description, the proposed model can be applied not only to a monolayer cell sheet, but also to a compacted cell aggregate. Moreover, the proposed model can be also applied to the tissue dynamics that can be affected by individual cell deformations; in this study, the deformations at the single cell level have little effects on tissue dynamics since we assumed the homogeneous, isotropic property for cell mechanical behaviors as in Eq.
(3). On the other hand, if we assume inhomogeneous and/or anisotropic property, morphogen patterns would be largely distorted in the tissue level via the flux in Eq. (2). Because of the cell-based description of the model, the method can be applied to the complex phenomena of cell mechanical behaviors (deformation, rearrangement, division, and apoptosis) with chemical behaviors (e.g., differentiation and growth) at the single cell level. Such a cell-based 3D description is useful in the analysis of the complex phenomena of multicellular patterning and deformation. Notably, the applicable area of the method is limited on a single cell level; intracellular patterning can be expressed over a length scale of cell-cell boundaries, which is a unit length of vertex models. The method therefore has the potential for remarkable predictions of undiscovered multicellular dynamics at a single cell level, which could guide future experimental approaches.
In summary, we demonstrated the computational simulation of multicellular deformations, coupled with reaction-diffusion processes in 3D space. As an example, we showed that coupling an activator-inhibitor system with cell growth caused various tissue deformations such as undulation, tubulation, and branching; the cell-based 3D analyses led to discovery of the importance of the time scales of patterning and deformation. The model has a significant applicability to 3D multicellular dynamics, compounding mechanical and chemical cell behaviors at a single cell level. Thus, we anticipate the combined model to open up a new avenue of elucidating yet undiscovered phenomena emerging from mechanochemical coupling in multicellular dynamics. Example results of the numerical simulations for the activator-inhibitor system with respect to the different initial conditions. The initial conditions were set as a random distribution with the uniform white noise, the boundary condition was set as  ρ =  ρ = 1 at x = 0, 2π, and physical parameters were set as  γ D / = 0.01 and  γ D / = 0.1. Red and blue lines indicate activator and inhibitor, respectively. (c) Result of the numerical simulation from the initial condition with dominant number of waves 4. (d) Result of the numerical simulation with domain growth from the initial condition with dominant number of waves 4. The domain size grew from 2π to 4π. In (c and d), the initial condition was set as periodic as  ρ =  ρ = 1 + sin4x.

Symbol Value Description
Δt v 0.