Kana Fuji

Sakurako Tanida

Masaki Sano

Makiko Nonomura

Daniel Riveline

Hisao Honda

Tetsuya Hiraiwa

Computational approaches for simulating luminogenesis

Keywords: Luminogenesis, Theoretical modeling, Computational simulations

Lumens, liquid-filled cavities surrounded by polarized tissue cells, are elementary units involved in the morphogenesis of organs. Theoretical modeling and computations, which can integrate various factors involved in biophysics of morphogenesis of cell assembly and lumens, may play significant roles to elucidate the mechanisms in formation of such complex tissue with lumens. However, up to present, it has not been documented well what computational approaches or frameworks can be applied for this purpose and how we can choose the appropriate approach for each problem. In this review, we report some typical lumen morphologies and basic mechanisms for the development of lumens, focusing on three keywords -mechanics, hydraulics and geometry -while outlining pros and cons of the current main computational strategies. We also describe brief guidance of readouts, i.e. what we should measure in experiments to make the comparison with the model's assumptions and predictions.

Introduction

Multicellular organisms are composed of tissues with complex shapes and specialized functions. They form cysts, branches, buds, tubes, convoluted structures with invagination for example. As a growth process starting with single cells, it is natural to ask how such complex tissue can emerge. An essential unit enabling such complexity is lumen -a liquid-filled cavity -surrounded by the polarized cells [START_REF] Blasky | Polarized Protein Transport and Lumen Formation During Epithelial[END_REF][START_REF] Navis | Pulling together: Tissue-generated forces that drive lumen morphogenesis[END_REF][START_REF] Honda | The world of epithelial sheets[END_REF]. Lumens are ubiquitous structural features in such complex cell assemblies, and those cell assemblies entirely grow through the dynamics of cell proliferation tightly coupled to the evolution of lumen dynamics. Revealing the processes at play in tissues growth and luminogenesis is indispensable to understand processes of organ formations during morphogenesis of a multicellular organism, as well as to predict and control morphologies of organoids, cell assembly culture mimicking an organ [START_REF] Nikolaev | Homeostatic mini-intestines through scaffold-guided organoid morphogenesis[END_REF][START_REF] Hofer | Engineering organoids[END_REF]. This synergy between cells mechanics and liquid circulation associated to hydraulics calls for theoretical approaches.

Complexity of lumen dynamics coupled with biomechanical and biochemical activities of cells makes it difficult to elucidate the mechanisms behind luminogenesis only by experimental observations. Theoretical investigations relying on mathematical models, which can compute dynamics of lumens and cells and reproduce various experimental data in a unified manner, are essential to challenge such complexity by integrating varieties of observations. To explain mechanisms underlying some real luminogenesis theoretically, the first step may be to find a mathematical model which can reproduce a target phenomenon. Once the reliable model has been realized, the mechanism will be clarified by elaborating its mathematical structure and computing its consequences comprehensively. A question then is how to construct such a mathematical model in a calculable form. Writing down all governing equations and making the model with full details which can perfectly reproduce the phenomenon might sound a reasonable strategy. However, in practice, this may be infeasible due to huge computation costs. Even if possible, because we know that theory must reproduce the phenomenon if all the details are implemented, we cannot learn something significant commensurate with the effort. Rather, the practical and effective strategy is to find a model which is simple enough but can reproduce certain aspects of features which one really wants to address. This purpose requires selecting the appropriate computational approach or framework, which can extract such target features sufficiently but simply enough.

Lumen simulation references

One specific example

Strength Weakness Simple geometry model approach (Sect. 4.1) [START_REF] Dasgupta | Physics of lumen growth[END_REF][START_REF] Duclut | Fluid pumping and active flexoelectricity can promote lumen nucleation in cell assemblies[END_REF][START_REF] Dumortier | Hydraulic fracturing and active coarsening position the lumen of the mouse blastocyst[END_REF][START_REF] Verge-Serandour | A hydro-osmotic coarsening theory of biological cavity formation[END_REF] Ref. [START_REF] Dasgupta | Physics of lumen growth[END_REF] studied dynamics of a lumen between two cells 1. High calculability 2. Can capture the most essential mechanics that we focus on.

1. Difficult to apply to simulate tissues with complex shapes 2. Lacks of some possible self-organized shapes of cells/lumens.

Cellular

Vertex model approach (Sect. 4.2) [START_REF] Honda | Computer simulation of emerging asymmetry in the mouse blastocyst[END_REF][START_REF] Rozman | Collective cell mechanics of epithelial shells with organoidlike morphologies[END_REF][START_REF] Gonay | Modelling of Epithelial Growth, Fission and Lumen Formation During Embryonic Thyroid Development: A Combination of Computational and Experimental Approaches[END_REF] Ref. [START_REF] Honda | Computer simulation of emerging asymmetry in the mouse blastocyst[END_REF] Ref. [START_REF] Akiyama | Numerical study on spindle positioning using phase field method[END_REF] studied how the rules determining the orientation of the spindle affect the lumen shape.

Cellular Potts model approach

(Sect. 4.4) [START_REF] Kim | Computational investigation of epithelial cell dynamic phenotype in vitro[END_REF][START_REF] Kim | A computational approach to resolve cell level contributions to early glandular epithelial cancer progression[END_REF][START_REF] Engelberg | MDCK Cystogenesis Driven by Cell Stabilization within Computational Analogues[END_REF][START_REF] Datta | Molecular Regulation of Lumen Morphogenesis[END_REF][START_REF] Cerruti | Polarity, cell division, and out-of-equilibrium dynamics control the growth of epithelial structures[END_REF][START_REF] Sluka | A Liver-Centric Multiscale Modeling Framework for Xenobiotics[END_REF][START_REF] Belmonte | Virtual-tissue computer simulations define the roles of cell adhesion and proliferation in the onset of kidney cystic disease[END_REF][START_REF] Hirashima | Dynamic modeling of branching morphogenesis of ureteric bud in early kidney development[END_REF][START_REF] Hirashima | Cellular Potts modeling of complex multicellular behaviors in tissue morphogenesis[END_REF][START_REF] Mombach | Bursts and cavity formation in Hydra cells aggregates: experiments and simulations[END_REF] Ref. [START_REF] Cerruti | Polarity, cell division, and out-of-equilibrium dynamics control the growth of epithelial structures[END_REF] studied how the timescale of cell cycles affect the lumen shape 1. Extensibility to add new parameters controling cell-cell/lumen interactions etc.

Dynamics of cells/lumens have a lot of flexibility

1. Cell-cell/-lumen interfaces are inevitably subjected to the fluctuation, and the result is easily affected by lattice points at the interfaces.

Not straightforward to implement force balance explicitly

Approaches associated with meshworking.

Brodland's finite element model, Immersed boundary model (IBM), Subcellular element model (SEM), etc. See other review literature [START_REF] Fletcher | Mechanocellular models of epithelial morphogenesis[END_REF] for more details. For luminogenesis, IBM [START_REF] Rejniak | A Computational Study of the Development of Epithelial Acini: I. Sufficient Conditions for the Formation of a Hollow Structure[END_REF][START_REF] Rejniak | A Computational Study of the Development of Epithelial Acini: II. Necessary Conditions for Structure and Lumen Stability[END_REF][START_REF] Rejniak | Computational investigation of intrinsic and extrinsic mechanisms underlying the formation of carcinoma[END_REF][START_REF] Rejniak | Homeostatic Imbalance in Epithelial Ducts and Its Role in Carcinogenesis[END_REF] and SEM [START_REF] Dokmegang | Computational modelling unveils how epiblast remodelling and positioning rely on trophectoderm morphogenesis during mouse implantation[END_REF][START_REF] Van Liedekerke | Quantitative modeling identifies critical cell mechanics driving bile duct lumen formation[END_REF] have been applied. For example, Ref. [START_REF] Rejniak | A Computational Study of the Development of Epithelial Acini: I. Sufficient Conditions for the Formation of a Hollow Structure[END_REF][START_REF] Rejniak | A Computational Study of the Development of Epithelial Acini: II. Necessary Conditions for Structure and Lumen Stability[END_REF] studied the conditions for normal acinar growth and maintenance of the hollow lumen generated by apoptosis, using IBM and introducing various cell state changes. Ref. [START_REF] Dokmegang | Computational modelling unveils how epiblast remodelling and positioning rely on trophectoderm morphogenesis during mouse implantation[END_REF] studied if the repulsive force at the apical surface of the rosette-shaped epiblast drives luminogenesis, using SEM. (Sect. 4.5) In this article, aiming at making a guide to select the computational approach, we review a few approaches or frameworks (listed up in Table 1) for computational modeling relevant for studies of luminogenesis. [footnote: Conventionally, the word "modeling" of some phenomenon can have double meanings; One is the process to wrap up all the biological/physical/chemical assumptions behind the phenomenon whereas the other is the process to implement those assumptions in a calculable form. In this review, the word "modeling" refers to the latter.] We mention strengths and weaknesses of each approach, which could help the readers to decide the best approach for their own specific purpose. After summarizing fundamental facts on lumens and luminogenesis in Section 2, we review biophysics ideas which determine lumen morphology and steers its dynamics in Section 3. Based on the described knowledge, we explain modeling frameworks which can be applied to simulate lumen dynamics leading to the rich variety of lumen's features. This analytical

Brief facts on lumens

Shapes of lumens

Typical shapes of lumens identified in vivo include spherical, tubular, multi-lumenal and network (Figure 1). The simplest structure is a sphere, e.g. thyroid follicles [START_REF] Nitsch | Suspension-Culture of Separated Follicles Consisting of Differentiated Thyroid Epithelial-Cells[END_REF][START_REF] Toda | Plural Cells Organize Thyroid-Follicles through Aggregation and Linkage in Collagen Gel Culture of Porcine Follicle Cells[END_REF][START_REF] Antonica | Generation of functional thyroid from embryonic stem cells[END_REF] form a spherical central lumen surrounded by a monolayer of epithelial cells. An early mouse embryo forms a lumen surrounded by asymmetric cell layers with both monolayers and multilayers at the blastocyst stage [START_REF] Ryan | Lumen Expansion Facilitates Epiblast-Primitive Endoderm Fate Specification during Mouse Blastocyst Formation[END_REF]. Another common lumen shape is a tubular structure, e.g. sweat glands have unbranched tubes [START_REF] Sato | Sweat Secretion by Human Axillary Apoeccrine Sweat Gland Invitro[END_REF] whereas mammary glands and lungs have branched tubes [START_REF] Hogan | Organogenesis: molecular mechanisms of tubulogenesis[END_REF]. Multilumenal structures within a multilayered stratified epithelium are observed in pancreas [START_REF] Flasse | Pancreas morphogenesis: Branching in and then out[END_REF]. Network luminal structures are observed in angiogenesis [START_REF] Folkman | Angiogenesis in vitro[END_REF][START_REF] Herbert | Molecular control of endothelial cell behaviour during blood vessel morphogenesis[END_REF] and in the lung [START_REF] Liu | Novel role for netrins in regulating epithelial behavior during lung branching morphogenesis[END_REF], where the network is formed by branching, and also in the pancreas [START_REF] Flasse | Pancreas morphogenesis: Branching in and then out[END_REF], where the network is formed by connecting multiple lumens.

These lumen structures have been reproduced in in-vitro systems, and it is reported that the lumen structures are altered by genetic manipulations and environment such as media. For example, MDCK cysts are normally monolayer spheres [START_REF] O'brien | Rac1 orientates epithelial apical polarity through effects on basolateral laminin assembly[END_REF][START_REF] Sigurbjornsdottir | Molecular mechanisms of de novo lumen formation[END_REF], but the Rab8a mutants of MDCK have multi-lumen structures [START_REF] Datta | Molecular Regulation of Lumen Morphogenesis[END_REF]. Moreover, spherical and tubular structures are mixed in some cases. For example, intestinal organoids have a spherical lumen in the initial stage. As they grow, several protrusions appear on cell monolayers, and then form tubular lumens. This results in the crypt-villus structure with multiple tubular lumens extending from the spherical lumen [START_REF] Sato | Single Lgr5 stem cells build crypt-villus structures in vitro without a mesenchymal niche[END_REF].

Epithelial cells

Epithelial cells can serve as physical barrier to allow or prevent the passage of molecules.

Composed of cells which are tightly connected and regulated by signaling networks, epithelia are essential to understand morphogenesis.

Epithelial cells are polarized in apico-basal manner. In 3D cysts, the apical side faces lumen, the basal side is in touch with extra-cellular matrix (ECM) made of crosslinked polymer networks, and the lateral sides are in mechanical contact with neighbouring cells. At the molecular levels, each side can be labelled with specific proteins which can be tracked [START_REF] Martin-Belmonte | Cell-polarity dynamics controls the mechanism of lumen formation in epithelial morphogenesis[END_REF]; e.g. podocalyxin is specifically localized at apical membranes and hence often regarded as an apical marker. Each side can be the site of significant modifications through signalling networks regulating the local cell mechanics and environment [START_REF] Bryant | A molecular network for de novo generation of the apical surface and lumen[END_REF]. Different forms of phosphatidylinositol, which is a component of the lipid bilayer and can be phosphorylated to form phosphatidylinositol bisphosphate (PIP2) and phosphatidylinositol trisphosphate (PIP3), are localized in apical and basolateral sides; PIP3 is at apical membrane whereas PIP2 is at basolateral membrane. Signaling network consisting of Cdc42, aPKC, PAR3, PAR6 etc. is relevant to form apico-basal polarity. Formation of apical membrane is mediated by a characteristic organelle within epithelial cells, called vacuolar apical compartment (VAC); VAC undergoes exocytosis on cell surface near areas of cell-cell contact and fuses to form apical membrane. ECM can induce basal membrane through the signaling via activation of Rac1 (RAS-related GTP binding protein) and integrins (see Ref. [START_REF] Honda | The world of epithelial sheets[END_REF] for more details). The actin cytoskeleton cortex spans the complete cell surface.

The interface of each side is related to specific interactions. At the apical side, the cell faces the lumen with pumps. At the basal sides, integrins bind to proteins of the ECM and the associated acto-myosin cortex will regulate the local mechanical state. At the lateral sides, epithelial cells are connected through various cell-cell junctional structures, such as adherens junctions and tight junctions [START_REF] Zihni | Tight junctions: from simple barriers to multifunctional molecular gates[END_REF][START_REF] Harris | Adherens junctions: from molecules to morphogenesis[END_REF][START_REF] Guillot | Mechanics of Epithelial Tissue Homeostasis and Morphogenesis[END_REF]. Here we focus on tight and adherens junctions in vertebrates. Tight junctions are located at the apical side and they play key roles for the barrier function of epithelia and gating the diffusion of molecules through paracellular space [START_REF] Zihni | Tight junctions: from simple barriers to multifunctional molecular gates[END_REF]. The composites are transmembrane proteins including claudin and occludin as well as adaptor proteins including ZO (zonula occludin) 1/2 etc. At the apical-most side near the tight junction, crumbs complex (Pals1/PatJ/Lin7c/Crb3) locates and, via a membrane compartment, controls cortical actin remodeling [START_REF] Tan | The Mammalian Crumbs Complex Defines a Distinct Polarity Domain Apical of Epithelial Tight Junctions[END_REF]. Adherens junctions are located more at the basal side. They do not work only as adhesive junction but also are involved in dynamic features as follows [START_REF] Harris | Adherens junctions: from molecules to morphogenesis[END_REF][START_REF] Guillot | Mechanics of Epithelial Tissue Homeostasis and Morphogenesis[END_REF]. A main component of adherens junction is transmembrane adhesion proteins E-cadherin, which is connected to the scaffolding bundles of actin filaments through alpha-and beta-catenins among other proteins. These actin filaments can also contain myosin and induce motor-dependent tension. Hence, the cell-cell junction behaves as if it has line tension in space-and time-specific manner, significantly contributing to morphogenesis [START_REF] Guillot | Mechanics of Epithelial Tissue Homeostasis and Morphogenesis[END_REF]. [Footnote: In invertebrate such as Drosophila, epithelial cells have septate junctions instead of tight junctions. Septate junctions are located at more basal side than adherens junctions, whereas tight junctions in a vertebrate are located at more apical side than adherens junctions [START_REF] Tanentzapf | Interactions between the crumbs, lethal giant larvae and bazooka pathways in epithelial polarization[END_REF].]

Processes of luminogenesis

Luminogenesis can be classified into two types. In one class, the cell layer deforms and encloses the environmental fluid. The deforming process is achieved by cell rearrangement and/or folding of the cell layer. This pattern can be found in the dorsal trachea branch of the Drosophila organ system [START_REF] Jazwinska | Epithelial tube morphogenesis during Drosophila tracheal development requires Piopio, a luminal ZP protein[END_REF]. In another type "de novo lumen formation", a lumen is formed inside or between cells. This lumen formation involves the establishment of apico-basal polarity and lumen expansion, as seen in organs, including intestine [START_REF] Bagnat | Genetic control of single lumen formation in the zebrafish gut[END_REF] and pancreas [START_REF] Flasse | Pancreas morphogenesis: Branching in and then out[END_REF], and in in vitro systems of blood vessels [START_REF] Folkman | Angiogenesis in vitro[END_REF] and thyroid follicle [START_REF] Toda | Plural Cells Organize Thyroid-Follicles through Aggregation and Linkage in Collagen Gel Culture of Porcine Follicle Cells[END_REF]. The de novo lumen formation can be further classified into three patterns according to the site of lumen formation: between cells, within a single cell, and a region where cells have died due to apoptosis [START_REF] Sigurbjornsdottir | Molecular mechanisms of de novo lumen formation[END_REF].

The process of luminogenesis in between two cells is well known. In a single cell before cell division, the apical components are distributed uniformly across the plasma membrane. When a cell divides, an apical membrane initiation site (AMIS) is formed at the centrosome [START_REF] Fremont | Membrane Traffic in the Late Steps of Cytokinesis[END_REF]. Then, endocytosis takes up the vesicles containing the apical surface molecules and is transported through microtubules to AMIS. The transported vesicles then fuse with AMIS via SNARE proteins to form the lumen. Once the apico-basal polarity is established, the apical surfaces of the two neighboring cells are separated from each other forming a lumen. The lumen is expanded and grows into a shape that fulfils physiological functions. This expansion is known to occur by two mechanisms. The first is the osmotic pressure generated by ions transported from the inside of the cell by ion channels and pumps on the apical surface [START_REF] Dasgupta | Physics of lumen growth[END_REF]. This osmotic pressure causes the inflow of water into the lumen, which drives lumen expansion. In the Zebrafish intestine, for example, ion transport by Na+/K+-ATPase contributes to lumen expansion, and indeed a large decrease in lumen size due to the lack of ion pump has been reported [START_REF] Bagnat | Genetic control of single lumen formation in the zebrafish gut[END_REF]. The second is the anti-adhesive properties between apical surfaces or the resultant conservation of apical area, due to specific molecules localizing on the apical surface [START_REF] Vasquez | Physical basis for the determination of lumen shape in a simple epithelium[END_REF]. Podocalyxin for MDCK and polysaccharides such as chitin for tracheal cells have negatively charged regions, which are antiadhesive factors and considered to contribute to opening a lumen [START_REF] Sigurbjornsdottir | Molecular mechanisms of de novo lumen formation[END_REF]. In addition, because cellcell adhesion molecules such as E-cadherin are excluded from the apical surface, the apical surface is considered to be non-adhesive [START_REF] Maraspini | Optimization of 2D and 3D cell culture to study membrane organization with STED microscopy[END_REF]. Indeed, it is known that loss of Rab35, which is involved in AMIS formation, can lead to the formation of multi-lumen structures [START_REF] Klinkert | Rab35 GTPase couples cell division with initiation of epithelial apico-basal polarity and lumen opening[END_REF], suggesting that perturbing formation of the apical surface can affect the luminal shape.

In a spatial environment sealed by confluent cells monolayers, a liquid can be enclosed. In this context, the appearance of lumen is essential for morphogenesis. For example, from the elongation of lumens (Fig. 1d,e), liquid may be shared between distant cells, and this could convey critical information for morphogenesis. Reflecting this fact, the laws associated with the transport of fluid should be also important.

We have reported the mechanical coupling between cells through adhesion and force generations together with their specificities at each site. We also described the nucleation and growth of lumens by different mechanisms. These phenomena act in concert and altogether they suggest that three biophysical factors enter synergistically into luminogenesis. Mechanics is involved when cells interact with their neighbours and with their specific environment at each side. Hydraulics will occur through the pumping of fluid in and out in particular between the lumen and the cells. Geometry, i.e. spatial confinement and volume occupancy, such as apical area conservation, will play also important roles. We further discuss these factors in the next section.

Biophysics behind Luminogenesis

Three biophysical factors

Here, we review three biophysical factors, mechanics (surface mechanics), hydraulics (hydraulics and osmosis) and geometry (spatial constraint and volume occupancy). These factors are essential to theoretically study luminogenesis.

Mechanics (Surface mechanics)

For formations of lumens surrounded by cell assemblies, confluent cells are important elements providing mechanical force, as described in Section 2.2. Firstly, cells adhere with each other to maintain the tissue integrity. Secondly, epithelial cells have bundles composed of actin filaments and myosin scaffolding the cortex and its adhesive junctions. Such bundles provide the contractile forces associated with myosin motor activities. For example, such junctional contraction is known to play essential roles for various dynamic morphogenetic processes [START_REF] Bertet | Myosin-dependent junction remodelling controls planar cell intercalation and axis elongation[END_REF][START_REF] Rauzi | Planar polarized actomyosin contractile flows control epithelial junction remodelling[END_REF][START_REF] Sato | Left-right asymmetric cell intercalation drives directional collective cell movement in epithelial morphogenesis[END_REF]; the readers interested in this point can refer to earlier review articles [START_REF] Guillot | Mechanics of Epithelial Tissue Homeostasis and Morphogenesis[END_REF][START_REF] Lecuit | Cell surface mechanics and the control of cell shape, tissue patterns and morphogenesis[END_REF]. Thirdly, cortical actomyosin cytoskeleton generates cortical stress which can be encoded into surface tension [START_REF] Graner | The Forms of Tissues, or Cell-aggregates': D'Arcy Thompson's influence and its limits[END_REF]. These adhesive and actomyosin forces will determine the lumen morphology obeying force balance for the integrated system of ECM-epithelial cells-lumen. In some systems, tissues may regulate their mechanical features to achieve the target morphology and impose the shape of lumens.

It is to be noted that, recently, such mechanics of a cytoskeletal system containing molecular motors, like actin filaments with myosin, has been investigated intensively in terms of the so-called active gel model [START_REF] Prost | Active gel physics[END_REF]. It is desirable to explicitly incorporate the active gel model into the frameworks presented in this article, but it is still challenging mathematically and computationally.

Hydraulics (Hydraulics and osmosis)

The most abundant components of cells and lumens are water molecules. Moreover, cells and lumens contain ions, which provide osmotic pressure and are essential to dictate water flows.

As discussed in Section 2.3, by controlling transports of ions and water molecules through channels (for both ions and water molecules) and pumps (for ions) on cell membranes as well as those flows through the pericellular spaces, tissue can eventually dictate the luminal size and shape. Recently, the importance of such hydraulic features of tissue for the expansion, ripening and size regulation of lumens have been recognized [START_REF] Dumortier | Hydraulic fracturing and active coarsening position the lumen of the mouse blastocyst[END_REF][START_REF] Gin | A model for cyst lumen expansion and size regulation via fluid secretion[END_REF][START_REF] Chan | Hydraulic control of mammalian embryo size and cell fate[END_REF]. In addition, hydraulics also plays key roles in single-cellular dynamics [START_REF] Li | The importance of water and hydraulic pressure in cell dynamics[END_REF] as long recognized along the studies of e.g. shape dynamics of a suspension cell [START_REF] Salbreux | Shape oscillations of non-adhering fibroblast cells[END_REF], cell volume regulation [START_REF] Jiang | Cellular Pressure and Volume Regulation and Implications for Cell Mechanics[END_REF][START_REF] Roffay | Passive coupling of membrane tension and cell volume during active response of cells to osmosis[END_REF] and cell volume regulation associated with deformation of the cell on a substrate [START_REF] Xie | Controlling Cellular Volume via Mechanical and Physical Properties of Substrate[END_REF][START_REF] Venkova | A mechano-osmotic feedback couples cell volume to the rate of cell deformation[END_REF]. It is worth noting that, as well recognized in the studies of plant growth, the concept of hydraulics has been essential in morphogenesis study from early days [START_REF] Thompson | On Growth and Form[END_REF] to present [START_REF] Long | Cellular Heterogeneity in Pressure and Growth Emerges from Tissue Topology and Geometry[END_REF]. Interestingly, the study of organs and organoids where lumens are central promoted the re-introduction of hydraulics in morphogenesis.

Geometry (Spatial confinement and volume occupancy)

Confinement and its associated 3D geometry provide varieties of other important regulation mechanisms. Force balance will have to be respected at the cellular side level, at the cell level during proliferation, at the system level where local constraints may differ in space and in time. Spatial rules mediated by volume occupancy of cortex, cells and components of the systems will regulate cell confinement and local geometry.

First, as we discussed in Section 2.3, regulation of apical domains is important and maintains an apical surface area. Indeed, in the initial lumen formation process of a MDCK cyst, such apical surface constraint was reported to be essential [START_REF] Vasquez | Physical basis for the determination of lumen shape in a simple epithelium[END_REF]. Second, the volume of each cell is an essential geometrical factor which affects the morphology of lumen and cell assembly. Each cell has a finite volume of matter composed of water, proteins, nucleus and organelles, which provides a minimum volume constraint. This volume increases during the cell cycle. In addition, there can be mechanisms for each cell to stop the cell cycle progression at the exit of G1 stage when the cell size is smaller than target volumes [START_REF] Ginzberg | On being the right (cell) size[END_REF], which should affect the cell proliferation and eventually the lumen shape. Third, in the developmental processes of some living organisms, such as mouse blastocyst and early Drosophila embryo, the cell assemblies are surrounded by the rigid wall-like substances, such as the zona pellucida and vitelline membrane [START_REF] Honda | Computer simulation of emerging asymmetry in the mouse blastocyst[END_REF][START_REF] Dicko | Geometry can provide long-range mechanical guidance for embryogenesis[END_REF][START_REF] Munster | Attachment of the blastoderm to the vitelline envelope affects gastrulation of insects[END_REF][START_REF] Bailles | Genetic induction and mechanochemical propagation of a morphogenetic wave[END_REF]. Multicellular assemblies have to satisfy all of these confinement and geometrical constraints when they grow and deform.

Luminogenesis as a combination of the biophysical factors

Luminogenesis will be the result of these force balance considerations, integrating mechanics of cells, hydraulics of lumens and cells, spatial constraints and volume occupancy. These factors may be independent or tightly coupled. This aspect can be taken into account in models. Altogether, self-organisation of the multi-cellular system with lumens will appear, and their emerging shapes will provide experimental readouts to test predictions when relevant readouts and functional dependencies are outlined.

Various theoretical modelling approaches.

In this section, we present several modelling approaches, or frameworks, to describe luminogenesis. We detail the four approaches, including one with a simple morphology model and three with free energy-based model, and comment on some meshworking-based computational approaches. Principles for describing the shapes and dynamics of lumens and cells are distinct in each case. We discuss the strengths and weaknesses of each approach associated with those distinct description principles. We also include a part for the concrete design of experiments associated with their readouts and functional dependencies. The Table 1 summarizes the examples of literature, in which each approach is applied, and the inherent pros and cons of each approach.

A central comparison common between models will be often the tracking of cells and lumens during the growth. The quantitative analysis of cell and lumens shapes will be then possible through various approaches. The selection of numbers for the physical parameters and prefactors involved in equations will also be important to test models and their strengths.

Simple morphology model approach

Model description

A common approach is to design the simplest theory which predicts and reproduces a specific aspect of the system. This allows us to identify what are the indispensable factors as well as for practical reasons such as analytical solvability and low computational cost. For the luminogenesis study, for example, one may assume the spherical shell of the cell assembly surrounding the spherical lumen(s), or someone else may assume the elliptical shapes or some given shapes represented by only a few variables, etc. For example in Ref. [START_REF] Duclut | Fluid pumping and active flexoelectricity can promote lumen nucleation in cell assemblies[END_REF], the lumen and cellular tissue are assumed as a sphere and the surrounding spherical shell, respectively. In this case, the morphology is identified by only two lengths, radii of lumen and shell's outer boundary.

Another example is found in Ref. [START_REF] Dasgupta | Physics of lumen growth[END_REF], where the cell doublet and a lumen sandwiched in between the two cells are considered. The lumen is modelled by two symmetrical spherical caps. The lumen elongates parallel to the cell-cell contact, and a paracellular space exists out of the lumen region. The lumen shape can be specified by the radius of curvature and a contact angle at the lumen edge.

One key advantage of this approach is its ease for calculations. With this simple morphology assumption, we can often construct the fully mechanically and hydraulically reliable model, which is indeed not always easy for the models with more demanding morphologies. In particular, it is interesting to see how Ref. [START_REF] Dasgupta | Physics of lumen growth[END_REF] modeled the mechanics and hydraulics. They assumed mechanical balances and conservations of ions and waters for both lumen and pericellular cleft region. For lumen, the mechanical balance is implemented by the equations expressing Laplace pressure and the Young-Dupre' law. For pericellular cleft, mechanical balance is given by the equation linking the hydrostatic pressure and cleft height. For conservation laws, ion-and water-fluxes in both lumen and pericellular regions and transmembrane fluxes through channels and pumps are considered. The theory is composed of these set of mechanical and hydraulic equations, and its steady state was analytically derived.

Based on this approach, the authors found the relevance of paracellular lumen leak as a determination factor of luminogenesis.

There are applications of this approach to multiluminal cases [START_REF] Dumortier | Hydraulic fracturing and active coarsening position the lumen of the mouse blastocyst[END_REF][START_REF] Verge-Serandour | A hydro-osmotic coarsening theory of biological cavity formation[END_REF], in which each (micro)lumen is represented by the mechanically-balanced simple morphology, like Ref. [START_REF] Dasgupta | Physics of lumen growth[END_REF] mentioned above, while those (micro)lumens are connected as a graph and interact by hydraulic fluxes with each other.

This approach is sometimes combined with the potential function inspired by the vertex model, which will be explained below, to simulate the shape of epithelial cell monolayer sheets [START_REF] Hannezo | Theory of epithelial sheet morphology in three dimensions[END_REF][START_REF] Yang | Cell fate coordinates mechano-osmotic forces in intestinal crypt formation[END_REF].

The budding (crypt) morphology in an intestine organoid was simulated by this method [START_REF] Yang | Cell fate coordinates mechano-osmotic forces in intestinal crypt formation[END_REF] as well.

Strengths/Weaknesses

This approach provides high calculability like analytical solvability and low computational cost. This approach can clearly demonstrate the essential mechanisms behind the target phenomenon. For example, Ref. [START_REF] Dasgupta | Physics of lumen growth[END_REF] was able to study the influence of paracellular leaks because the fluid transport was fully taken into account. However, since this approach assumes some given shape represented by a few variables ab initio, the possible morphologies of the cell assemblies and lumens are restricted. This prohibits the emergence of unexpected complex morphologies due to self-organization. In other words, the theoretical results can not realize some possible lumen morphology and morphodynamics which can appear in real systems.

Experimental readouts for comparison

For this approach, possible readouts in experiments for comparison highly depend on the specific cases, depending on the desired model variables like lumen radius. As an example, lumen shapes might be tracked in space and in time; e.g. In the context of Ref. [START_REF] Dasgupta | Physics of lumen growth[END_REF], the oscillatory behaviour would be revealed as a function of ion and water transport parameters, ideally measured experimentally as well on the same system. Physical parameters of the model associated to transport and to mechanics can be measured experimentally as well.

Cellular vertex model approach

Model description

To simulate dynamics of a multicellular assembly, the model called 'vertex model' is often used. This is the theoretical framework to describe dynamics of tightly packed objects with flat interfaces in general, or various "cellular" patterns in nature, i.e., not only cells in a biological tissue but also foams, grain aggregates and magnetic domains. The essential assumption of this model is that the shapes of some objects (like cells) are described as polygons (polyhedral in three dimensions). The equilibrium state and dynamics are, in the modern form, described by using equations of locations and motion of vertices, respectively, associated with a potential function and the edge reconnection rules. Historically, these concepts of the vertex model were introduced several tens of years ago [START_REF] Fullman | Metal Interfaces[END_REF] and developed to simulate the dynamics of growing domains such as grain growth and the evolution of foam [START_REF] Nagai | Vertex Dynamics of Two-Dimensional Cellular-Patterns[END_REF][START_REF] Enomoto | TWO-DIMENSIONAL VERTEX MODEL WITH LOCAL FRICTION COEFFICIENT[END_REF][START_REF] Kawasaki | Vertex Models for Two-Dimensional Grain-Growth[END_REF]. Later, the vertex model is applied to the simulation of the biological tissue dynamics as well [START_REF] Nagai | A dynamic cell model for the formation of epithelial tissues[END_REF][START_REF] Honda | A three-dimensional vertex dynamics cell model of spacefilling polyhedra simulating cell behavior in a cell aggregate[END_REF]. Note that the prototypical idea of using vertex motion to simulate the epithelial cells had been already provided in much earlier time [START_REF] Honda | How Much Does the Cell Boundary Contract in a Monolayered Cell Sheet[END_REF][START_REF] Honda | Geometrical Models for Cells in Tissues[END_REF][START_REF] Weliky | The Mechanical Basis of Cell Rearrangement .1. Epithelial Morphogenesis during Fundulus Epiboly[END_REF]; the details for such historical aspects are presented in Refs. [START_REF] Fletcher | Vertex Models of Epithelial Morphogenesis[END_REF][START_REF] Honda | Cell models lead to understanding of multi-cellular morphogenesis consisting of successive self-construction of cells[END_REF]. Further later, cellular vertex models are investigated in more detail from the point of view of statistical physics [START_REF] Farhadifar | The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing[END_REF][START_REF] Staple | Mechanics and remodelling of cell packings in epithelia[END_REF][START_REF] Bi | A density-independent rigidity transition in biological tissues[END_REF].

This model is mathematically formulated in the following way. The conformation of a polygon is represented by specifying the locations of the vertices 𝑟 𝑖 ⃗ ⃗ (for 𝑖 -th vertex) in two or three dimensions, and the junctions < 𝑘𝑙 >, which means that there is an edge between the vertices 𝑘 and 𝑙. The time evolution equation of the vertex model is given by considering the force balance between the force derived from a potential 𝐸({𝑟 𝑖 ⃗ ⃗ }) and the friction force. For example, the simple linear friction case is given by 𝜂 𝑖 𝑑𝑟 𝑖 ⃗ ⃗ /𝑑𝑡 = -𝜕𝐸({𝑟 𝑖 ⃗ ⃗ })/𝜕𝑟 𝑖 with the friction coefficient 𝜂 𝑖 . We can simply interpret this system that the vertices move so that the potential 𝐸({𝑟 𝑖 ⃗ ⃗ }) becomes small, because 𝜕𝐸({𝑟 𝑖 ⃗ ⃗ })/𝜕𝑡 = 𝜕𝐸({𝑟 𝑖 ⃗ ⃗ })/𝜕𝑟 𝑖 • 𝑑𝑟 𝑖 ⃗ ⃗ /𝑑𝑡 = -𝜂 𝑖 |𝑑𝑟 𝑖 ⃗ ⃗ /𝑑𝑡| 2 ≤ 0 according to the above equation of motion. By giving the specific form of 𝐸({𝑟 𝑖 ⃗ ⃗ }), we can simulate various processes. Furthermore, in most cases, we also allow junctional remodelling, which leads to cell rearrangement by cell intercalations. For example in the two-dimensional case, the junctional remodeling is usually implemented in the following manner. When the length of a cell-cell junction becomes smaller than a given threshold value, the edge is rotated by 90 degrees around its midpoint, and simultaneously the five edges connected to either of two vertices at the JR site are reconnected such that the T1 transition is achieved (see e.g. Refs. [START_REF] Nagai | A dynamic cell model for the formation of epithelial tissues[END_REF][START_REF] Weliky | The Mechanical Basis of Cell Rearrangement .1. Epithelial Morphogenesis during Fundulus Epiboly[END_REF]). Such junctional remodelling is essential to describe the cell rearrangement processes.

The vertex model is a powerful tool to simulate the epithelial tissue dynamics. In two dimensions, since the apical face of each epithelial cell is surrounded by contractile acto-myosin bundles as we described in Section 2.2, it is natural to apply the essential assumption of the vertex model that each interface is flat [START_REF] Nagai | A dynamic cell model for the formation of epithelial tissues[END_REF]. The cell-monolayer dynamics driven by junctional contractility has been widely studied with the 2D vertex model, e.g. for convergent extension [START_REF] Honda | Two different mechanisms of planar cell intercalation leading to tissue elongation[END_REF][START_REF] Collinet | Local and tissue-scale forces drive oriented junction growth during tissue extension[END_REF] and unidirectional tissue flow [START_REF] Sato | Left-right asymmetric cell intercalation drives directional collective cell movement in epithelial morphogenesis[END_REF][START_REF] Hiraiwa | Wave Propagation of Junctional Remodeling in Collective Cell Movement of Epithelial Tissue: Numerical Simulation Study[END_REF]. As a variation, tissue invagination was also simulated based on the 2D vertex model, by focusing on the cross-section of an epithelial monolayer along the apical-basal axis [START_REF] Wen | Epithelial Folding Driven by Apical or Basal-Lateral Modulation: Geometric Features, Mechanical Inference, and Boundary Effects[END_REF]. In three dimensions, the vertex model assumes polyhedral shapes of cells with flat interfaces. The 3D vertex model has been applied to simulate the tissue morphology and its dynamics extensively [START_REF] Honda | Computer simulation of emerging asymmetry in the mouse blastocyst[END_REF][START_REF] Rozman | Collective cell mechanics of epithelial shells with organoidlike morphologies[END_REF][START_REF] Honda | A three-dimensional vertex dynamics cell model of spacefilling polyhedra simulating cell behavior in a cell aggregate[END_REF][START_REF] Odell | The Mechanical Basis of Morphogenesis .1. Epithelial Folding and Invagination[END_REF][START_REF] Okuda | Reversible network reconnection model for simulating large deformation in dynamic tissue morphogenesis[END_REF][START_REF] Okuda | Vertex dynamics simulations of viscositydependent deformation during tissue morphogenesis[END_REF][START_REF] Okuda | Strain-triggered mechanical feedback in self-organizing optic-cup morphogenesis[END_REF][START_REF] Okuda | Combining Turing and 3D vertex models reproduces autonomous multicellular morphogenesis with undulation, tubulation, and branching[END_REF][START_REF] Messal | Tissue curvature and apicobasal mechanical tension imbalance instruct cancer morphogenesis[END_REF]. Since increasing numbers of new works arise recently, which use cellular vertex models to study biological processes, it is difficult to be complete. The interested readers may refer to other review articles e.g. [START_REF] Fletcher | Vertex Models of Epithelial Morphogenesis[END_REF][START_REF] Honda | Cell models lead to understanding of multi-cellular morphogenesis consisting of successive self-construction of cells[END_REF][START_REF] Okuda | Three-dimensional vertex model for simulating multicellular morphogenesis[END_REF][START_REF] Alt | Vertex models: from cell mechanics to tissue morphogenesis[END_REF][START_REF] Hiraiwa | Mathematical Modeling of Tissue Folding and Asymmetric Tissue Flow during Epithelial Morphogenesis[END_REF].

Multicellular systems with lumens can be also simulated based on the vertex model by assuming that several polygons mimic lumens, which have distinct properties from the cellular polygons [START_REF] Honda | Computer simulation of emerging asymmetry in the mouse blastocyst[END_REF][START_REF] Rozman | Collective cell mechanics of epithelial shells with organoidlike morphologies[END_REF][START_REF] Gonay | Modelling of Epithelial Growth, Fission and Lumen Formation During Embryonic Thyroid Development: A Combination of Computational and Experimental Approaches[END_REF]. As an example, we focus on the simulation of cavity (or lumen) size and location in mouse blastocyst [START_REF] Honda | Computer simulation of emerging asymmetry in the mouse blastocyst[END_REF]. The question was about how embryonic polarity such as positioning asymmetry of cavity, or asymmetric positioning of the inner cell mass, is established in mammals.

To give an answer, Ref. [START_REF] Honda | Computer simulation of emerging asymmetry in the mouse blastocyst[END_REF] applied the 3D vertex model in which one of the polyhedrons represents the cavity to simulate growth dynamics and equilibrium positioning of the cavity in the mouse blastocyst (Fig. 2a). During mammalian development, a fertilized cell repeats cell division and forms an aggregate of cells (morula). The cell aggregate undergoes a remarkable change from the morula stage to the blastocyst stage, whereby blastocyst formation exhibits the appearance of a large lumen in the cell aggregate of the morula. For simulation of a cell aggregate, we usually used the potential 𝐸 consisting of the terms of the surface energy and the elastic energy of the interfaces between neighbouring cells. However, for simulation of the cell aggregate involving the lumen, Ref. [START_REF] Honda | Computer simulation of emerging asymmetry in the mouse blastocyst[END_REF] included the new term of the elastic energy of the lumen volume. When the lumen reaches the target volume, the cell aggregate relaxes. Ref. [START_REF] Honda | Computer simulation of emerging asymmetry in the mouse blastocyst[END_REF] sets the condition that the target volume would increase with time because of secretion of liquid from the cells, and that the lumen volume would continue to increase during the developmental process. A result of simulation is shown in Fig. 2a,b. We can distinguish between two types of cells, cells of the inner cell mass (dotted polyhedrons in Fig. 2c) and cells of the trophectoderm (white polyhedrons in Fig. 2c). Cells of the inner cell mass among the two types of cells cause asymmetric positioning of the lumen [START_REF] Honda | Cell models lead to understanding of multi-cellular morphogenesis consisting of successive self-construction of cells[END_REF]. Thus, these simulations clarified that geometrical constraint at the outer boundary, i.e. the existence of zona pelucida (ZP), or a rigid elliptical capsule around the cell assembly, with a few particular assumptions can autonomously position the cavity properly.

Fig. 2 Computer simulation of mammalian blastocyst morphology. Drawing (top left) shows cross-sections in xz-and yz-planes of simulated aggregates. (a) An example of calculation process to a stable state viewed from cross-sections of the yz-plane. Numbers indicate the time point of the simulation (t). A vertex (a black circle in the sample at t=0) is replaced by a small tetrahedron and is enlarged until its volume reaches half the initial total volume (at 𝑡 = 500). The blastocyst axis keeps changing during

Strengths/Weaknesses

The first advantage of using vertex models is the computational cost. Since this model calculates motions of only vertices, the computational cost is much lower compared with the phase-field reported below. Even analytical arguments are sometimes available [START_REF] Okuda | A Mechanical Instability in Planar Epithelial Monolayers Leads to Cell Extrusion[END_REF][START_REF] Harmand | 3D Shape of Epithelial Cells on Curved Substrates[END_REF]. Another advantage lies in the model assumptions. It is straightforward for the vertex model to take into account junctional contractile forces of epithelial tissue. One drawback stems from its essential hypothesis. Since the curved interfaces cannot be described by the vertex description, the possible lumen-(and cell-)shapes are limited. For example, it is known that there are situations where microlumens are initiated between two cells, which the vertex model cannot deal with. Lumen morphodynamics at the interface of two cells in a cell doublet, which we saw in the Subsection 4.1 with the simple morphology model, cannot be simulated by the vertex model. A variation called bubbly vertex model [START_REF] Ishimoto | Bubbly vertex dynamics: A dynamical and geometrical model for epithelial tissues with curved cell shapes[END_REF] might become a solution to challenge this point. Another drawback is the difficulty to incorporate the density distribution of signal molecules. It is not straightforward to combine the continuum field description of the chemicals and pointposition description of the vertices. Furthermore, cell rearrangement or T1 transition has to be implemented additionally as a rule independent from the equation of motions, thus can introduce an arbitrariness.

Experimental readouts for comparison

In the vertex model, it is important to follow the 2D or 3D contours of cells and lumens in time and in space. This allows quantifying the cell deformation and the cell rearrangement like neighbour exchange or delamination. We do not have to track the detailed shapes to compare with vertex model descriptions, in which cells and lumens are modelled by polygons/polyhedrons. The terms in the energy could contain the physical model parameters like junctional tension. They can be measured by mechanical methods like laser ablation experiments [START_REF] Kong | Experimental validation of force inference in epithelia from cell to tissue scale[END_REF].

Multicellular phase-field model approach

Model description

Another way to represent the multicellular system with lumens is the phase-field model. The basic idea of this framework is to introduce the continuum field in space to describe the shape of the interface (Fig. 3). For example, by specifying the closed interface between intra-and extracellular spaces using this framework, the cell shape can be fully described. To our knowledge, the phase-field model was originally introduced in 1980s [START_REF] Fix | Free Boundary Problems: Theory and Applications[END_REF][START_REF] Langer | MODELS OF PATTERN FORMATION IN FIRST-ORDER PHASE TRANSITIONS[END_REF], followed by several works studying mechanisms behind various morphodynamics in physical systems, including dendritic crystal growth [START_REF] Kobayashi | Modeling and Numerical Simulations of Dendritic Crystal-Growth[END_REF], motion of air bubbles and falling water drops [START_REF] Sussman | A Level Set Approach for Computing Solutions to Incompressible 2-Phase Flow[END_REF], fracture [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF] and grain growth [START_REF] Krill | Computer simulation of 3-D grain growth using a phase-field model[END_REF]. The phase-field model framework was applied to template morphology of a living cell [START_REF] Kockelkoren | Computational approach for modeling intra-and extracellular dynamics[END_REF], and extended to incorporate various intracellular components, such as actin and signaling molecule distributions, and simulate complex behaviors of cells [START_REF] Shao | Computational Model for Cell Morphodynamics[END_REF][START_REF] Ziebert | Model for self-polarization and motility of keratocyte fragments[END_REF][START_REF] Taniguchi | Phase geometries of two-dimensional excitable waves govern self-organized morphodynamics of amoeboid cells[END_REF][START_REF] Saito | Three-dimensional morphodynamic simulations of macropinocytic cups[END_REF]. Multicellular systems are also simulated using such cellular phase-field models [START_REF] Akiyama | Numerical study on spindle positioning using phase field method[END_REF][START_REF] Nonomura | Study on Multicellular Systems Using a Phase Field Model[END_REF][START_REF] Camley | Polarity mechanisms such as contact inhibition of locomotion regulate persistent rotational motion of mammalian cells on micropatterns[END_REF][START_REF] Lober | Collisions of deformable cells lead to collective migration[END_REF][START_REF] Mueller | Emergence of Active Nematic Behavior in Monolayers of Isotropic Cells[END_REF][START_REF] Kuang | Computable early Caenorhabditis elegans embryo with a phase field model[END_REF]. For more details of historical aspects, several review articles are available [START_REF] Ziebert | Computational approaches to substrate-based cell motility[END_REF][START_REF] Moure | Phase-Field Modeling of Individual and Collective Cell Migration[END_REF].

The phase-field model is a continuum model developed to describe interfacial dynamics by a scalar field. Here, we use the notation 𝑢(𝑟 , 𝑡) for the phase field, where 𝑟 and 𝑡 represent space and time, respectively. The values of the field evolve to make the energy potential lower. For multicellular systems, we prepare the multiple numbers of the scalar fields, in which the inside and outside of 𝑖-th cell are represented as 𝑢 𝑖 = 1 and 𝑢 𝑖 = 0, respectively. Although there are various formalisms for the dynamics of phase fields, for example, a standard formalism assumes the dynamics of a field of each cell 𝑖 to obey 𝜕𝑢 𝑖 /𝜕𝑡 = -𝛿𝐹({𝑢 𝑖 (𝑟 , 𝑡)})/𝛿𝑢 𝑖 (𝑟 ) with 𝐹({𝑢 𝑖 (𝑟 , 𝑡)}) = 𝐹 𝑃𝐹 ({𝑢 𝑖 (𝑟 , 𝑡)}) + 𝐹 𝑝ℎ𝑦𝑠 ({𝑢 𝑖 (𝑟 , 𝑡)}), where 𝐹 𝑃𝐹 is the potential to realize the basic assumption of the phase field mentioned above and in Fig. 3A, and 𝐹 𝑝ℎ𝑦𝑠 is the potential to introduce various physical factors, like steric repulsion and adhesions [START_REF] Nonomura | Study on Multicellular Systems Using a Phase Field Model[END_REF]. We skip mathematical details here, and the interested readers may refer to Ref. [START_REF] Nonomura | Study on Multicellular Systems Using a Phase Field Model[END_REF].

Figure 3 Simulations of luminogenesis with phase field. (a) An isolated cell described with the variable u of the phase-field model. The region where 𝑢 ∼ 1 represents the space occupied by a cell. (b) An example of the multicellular simulation with a lumen in 3D. (c) The function of the cortical force generator (CFG) in multicellular simulation. With CFG localized on the cell-cell adhesion area (top), micro-lumens which appeared at the center of each cell division plane gather and form a large single lumen. On the other hand with uniform CFG (bottom), multiple lumens emerged. Source [13].

This phase field model is applicable to luminogenesis simulation [START_REF] Akiyama | Numerical study on spindle positioning using phase field method[END_REF] (Fig. 3b). To consider growing tissues, cell divisions may be introduced as follows. Some conditions for cell division, like timing or cell size, may be assumed and the phase field of each cell satisfying the conditions are divided into two fields representing daughter cells. The division plane may be determined by considering the spindle position and orientation. Ref. [START_REF] Akiyama | Numerical study on spindle positioning using phase field method[END_REF] introduced the cortical force generator, which pulls the spindle poles, and assumed that spindle position and orientation are determined by the tug-of-war of pulling forces from such force generators. The lumen is also represented by a phase field. Although the dynamics of a lumen phase field is assumed similarly to that of a cell, it is different from cells phases in the following aspects. First, luminal fluid grows without target volume whereas a cell is assumed to have an optimum volume. Second, all isolated regions of lumens are described by a single phase field and hence automatically fuse when they contact each other because all luminal fluid is described with the same parameter. As an advantage of the phase-field model approach, we can introduce the microlumen generated at the middle of the two cells (Fig. 3c: left). With these rules of lumen creation and growth, how cavities grow in the tissue can be simulated. Ref. [START_REF] Akiyama | Numerical study on spindle positioning using phase field method[END_REF] performed computer simulations with different distributions of cortical force generators and investigated how their distributions affect morphology-in particular, single central lumen or multi-lumensafter the cyst grows. They found that, when cortical force generators are activated specifically at the lateral cortices, a single central lumen emerges autonomously after the growth. In contrast, when the cortical force generators are rather uniformly distributed, multiple lumens emerge (see Fig. 2c). Thus, they concluded that cortical force generators activated specifically at the lateral cortices result in the formation of a single central lumen.

Strengths/Weaknesses

The phase-field model is appropriate when we need to consider the arbitrary cell shape. Furthermore, this model can be easily extended to include spatio-temporal distribution of specific molecules/proteins which affects cell shape and adhesion, because this model ab initio has the space spanning everywhere both in and out of the cells/lumens. In the phase field model, potential terms encoding for cell adhesion, luminal pressure for example, can be included, and these terms can be changed in simulations and in experiments to test quantitative comparisons in shapes. On the other hand, the phase-field model requires a high computational cost because the space should be gridded finely to simulate this model, and all grid points need to be calculated separately. The phase field model as well as vertex model includes the minimisation of an energy, and the effects which cannot be written as an energy potential are difficult to be included.

We comment on usage of the phase field model for simulation with a particular emphasis of mechanical features beyond their control in simulations and in experiments. First, the phase field model can satisfy Young-Laplace law and can be solved simultaneously with the fluid dynamics equation, which have indeed been demonstrated in e.g. Ref. [START_REF] Vakili | Controlling bubble coalescence in metallic foams: A simple phase field-based approach[END_REF]. Since fluid dynamics should be solved to take into account e.g. leaks of lumen [START_REF] Dasgupta | Physics of lumen growth[END_REF], this advantage is likely relevant for the morphodynamics study. Second, a particular care is needed on a subtlety related with the magnitude of surface tension. The phase interface implemented by the phase field model has an unavoidable surface tension to prepare the smooth profile at the interface. This surface tension vanishes in the narrow interface-width limit ideally. However, computationally, since we cannot take the exact narrow interface-width limit, such "artificial" surface tension cannot be ignored. After 2000, the new computational method, called resharpening in some literature, was proposed to guarantee the constant profile and thickness of the phase field at the interface [START_REF] Olsson | A conservative level set method for two phase flow[END_REF][START_REF] Olsson | A conservative level set method for two phase flow II[END_REF][START_REF] Badillo | Quantitative phase-field modeling for boiling phenomena[END_REF][START_REF] Badillo | Quantitative phase-field modeling for wetting phenomena[END_REF]. It is known that implementing this method also removes the artificial surface tension mentioned above. The resharpening method may be a useful method in the case when we need particular care about the mechanistic relations.

Experimental readouts for comparison

Cells and lumens can be tracked with their contours over time and space. Cell shapes can be quantified back-to-back in simulations and in experiments to be compared back-to-back with controlled parameters for cells cortex, adhesion and pumping, ideally measured also in experiments. Detailed quantities of lumen and cell shapes like interface curvature may help the comparison. For this purpose, lumen shapes, or apical membrane of the cells, might be tracked in space and in time. The typical markers for cell membrane and for cell junctions can be labelled fluorescently and outlines can be skeletonized and quantified.

Cellular Potts model approach

Model description

Computational framework which was and is playing important roles in the biophysics field is Cellular Potts model (CPM). Potts model has been used to simulate grain growth dynamics and coarsening of foam driven by surface tension [Footnote: This is a generalization of the Ising model, and the spin with two states, or up/down, placed on each lattice point is replaced by some multi-state vector.]. Where different grain (foam) is expressed by different state 𝜎 𝑖 (𝑖 represents the coordinate on the lattice point and all points belonging to the same grain have the same value of 𝜎). To represent many grains, a large number of states is required. Large state Potts model was extended to simulate multi-cellular dynamics driven by differential adhesion between different types of cells [START_REF] Graner | Simulation of biological cell sorting using a two-dimensional extended Potts model[END_REF][START_REF] Glazier | Simulation of the Differential Adhesion Driven Rearrangement of Biological Cells[END_REF]. As assumed in VM and PFM, CPM algorithm evolves to minimize the total energy of the system. The total energy (Hamiltonian in CPM) can include cell-cell (medium) interactions, cell volume constraint, and chemotaxis. CPM formulated as follows. Some "feature" (state) 𝜎 (𝜎 = 𝑐𝑒𝑙𝑙 1, 𝑐𝑒𝑙𝑙 2, 𝑐𝑒𝑙𝑙 3, …) is placed on each site of a lattice and each state interacts with the state on the adjacent sites [START_REF] Potts | Some generalized order-disorder transformations[END_REF][START_REF] Wu | The Potts model[END_REF][START_REF] Glazier | Grain growth in three dimensions depends on grain topology[END_REF][START_REF] Sire | Growth Laws for 3d Soap Bubbles[END_REF]. Hamiltonian of CPM is defined as, 𝐻 = ∑ 𝑖,𝑗∈𝑁 𝐽(𝜏(𝑖)𝜏′(𝑗))𝐶 𝑖𝑗 -∑ 𝑐𝑒𝑙𝑙 𝜎 𝜆(𝜎)(𝑣(𝜎) -𝑉 𝜏 (𝜎)) 2 .

Here, 𝐽(𝜏(𝑖)𝜏′(𝑗)) is the interaction energy between cell type 𝜏 and 𝜏 ′ . 𝜏(𝑖) denotes the type of the cell sitting on the site i. 𝐶 𝑖𝑗 =0 if the site i and j belong to the same cell 𝜎, otherwise 𝐶 𝑖𝑗 = 1 and 𝐽(𝜏(𝑖)𝜏′(𝑗)) enumerates the interfacial energy. 𝑣(𝜎) is the volume of cell 𝜎. 𝑉 𝜏(𝜎) is the target volume of the cell type 𝜏 (in 2D, volume should be change to area). In CPM, update dynamics of the states is stochastic (so-called Monte Carlo algorithm) which is distinct from the deterministic dynamics in VM and PFM. At each step, neighbouring sites (I,j) are selected randomly and change of the energy 𝛥𝐻 by flipping the state of 𝑖 from 𝜎 to 𝜎 ′ is calculated If it reduces the energy (𝛥𝐻<0), flip the state from 𝜎 to 𝜎 ′ with the probability 1, otherwise flip with the probability 𝑝 = 𝑒 -𝛥𝐻/𝑇 . Here T is the effective temperature representing fluctuations driving cell boundary motion. Typically, one Monte Carlo time step (MCS) requires 10-20 times the number of lattice points [START_REF] Graner | Simulation of biological cell sorting using a two-dimensional extended Potts model[END_REF][START_REF] Glazier | Simulation of the Differential Adhesion Driven Rearrangement of Biological Cells[END_REF][START_REF] Rens | From energy to cellular forces in the Cellular Potts Model: An algorithmic approach[END_REF][START_REF] Fortuna | CompuCell3D Simulations Reproduce Mesenchymal Cell Migration on Flat Substrates[END_REF][START_REF] Bernadskaya | Supracellular organization confers directionality and mechanical potency to migrating pairs of cardiopharyngeal progenitor cells[END_REF] .

To give an insight, an example of modelling a lumen surrounded by several cells is depicted in Fig. 3(a). Here, a lattice of 15x12 sites in 2D is defined. Each number or character on the site represents the state of the site 𝜎. 7 cells (𝜎 = 1,2, . . .7) are forming a monolayer surrounded by ECM (𝜎 = 𝑒) and contains a lumen denoted by 𝜎 = 0 at the centre. Size of the cell has an average of 11 and fluctuations around it. For ECM (𝜎 = 𝑒), depending on the stiffness of ECM, the volume constraint coefficient 𝜆 can be adjusted. Figure 3(b) and 3(c) shows results of cell sorting simulation by CPM. Here, a rounded aggregate of two types of cells (d=dark and l=light) are immersed in the medium as an initial condition. After 10000 MCS, two types of cells are separated due to differential adhesion. For the spontaneous phase separation, the interaction energies must satisfy the relation, 𝐽(𝑑, 𝑑) < (𝐽(𝑑, 𝑑) + 𝐽(𝑙, 𝑙))/2 < 𝐽(𝑑, 𝑙) < 𝐽(𝑙, 𝑙) < 𝐽(𝑙, 𝑀), 𝐽(𝑑, 𝑀) [START_REF] Glazier | Simulation of the Differential Adhesion Driven Rearrangement of Biological Cells[END_REF]. CPM has been applied for different situations [START_REF] Kim | Computational investigation of epithelial cell dynamic phenotype in vitro[END_REF][START_REF] Kim | A computational approach to resolve cell level contributions to early glandular epithelial cancer progression[END_REF][START_REF] Engelberg | MDCK Cystogenesis Driven by Cell Stabilization within Computational Analogues[END_REF][START_REF] Datta | Molecular Regulation of Lumen Morphogenesis[END_REF][START_REF] Cerruti | Polarity, cell division, and out-of-equilibrium dynamics control the growth of epithelial structures[END_REF][START_REF] Sluka | A Liver-Centric Multiscale Modeling Framework for Xenobiotics[END_REF][START_REF] Belmonte | Virtual-tissue computer simulations define the roles of cell adhesion and proliferation in the onset of kidney cystic disease[END_REF][START_REF] Mombach | Quantitative Comparison between Differential Adhesion Models and Cell Sorting in the Presence and Absence of Fluctuations[END_REF].

Mechanism for cyst lumen formation and size regulation is investigated using CPM [START_REF] Engelberg | MDCK Cystogenesis Driven by Cell Stabilization within Computational Analogues[END_REF][START_REF] Belmonte | Virtual-tissue computer simulations define the roles of cell adhesion and proliferation in the onset of kidney cystic disease[END_REF]. Lumen growth law is assumed that the growth rate of the target volume of lumen is proportional to (𝑁 𝑐𝑒𝑙𝑙𝑠 -𝜅𝑆 𝑙𝑢𝑚𝑒𝑛 ), where 𝑁 𝑐𝑒𝑙𝑙𝑠 is number of cells in contact with the lumen, 𝜅 is a constant related to leakage, 𝑆 𝑙𝑢𝑚𝑒𝑛 is the surface area of the lumen 𝜅 is a constant that defines the degree of leakage and the extent of stretch. Cystogenesis was simulated with CPM starting from a single cell [START_REF] Belmonte | Virtual-tissue computer simulations define the roles of cell adhesion and proliferation in the onset of kidney cystic disease[END_REF] (Fig. 4(d)A). Contact inhibition for proliferation was introduced with defining the cell contact fraction dependent target volume of the cell. The cell contact fraction was defined as the cell's surface area that is in contact with other cells. Cell cleavage plane was chosen as perpendicular to the plane of the epithelium in random orientation with respect to the tubule axis. Cystic lumen nucleates at any location, contacting three apical compartments of cells not currently in contact with lumen. Once created, growth rate of the target volume of the lumen is determined by 𝑁 𝑐𝑒𝑙𝑙𝑠 and 𝑆 𝑙𝑢𝑚𝑒𝑛 . Depending on the balance between cell proliferation and lumen growth, balance between ECM stiffness and cell-cell adhesion, different morphologies appeared, i.e. spherical lumen and stop growing (Fig. 4(d)B), spherical lumen grow indefinitely (Fig. 4(d)C), or split into multiple cysts (Fig. 4(d)D).

Strengths/Weakness

CPM has a strong extensibility to add new parameters controlling cell-cell/-lumen interactions etc. Since Potts model framework is based on the model quite intensively investigated in the statistical physics field historically, we can take advantage of its abundant knowledge. Furthermore, dynamics of shapes of cells/lumens have large flexibility under the CPM framework. There are also drawbacks. Cell-cell/-lumen interfaces are inevitably subjected to the fluctuation, and the result is easily affected by lattice points at the interface. In addition, as same as the multicellular phase field model case (Section 4.3.2), it is not straightforward to implement force balance in CPM.

Experimental readouts for comparison

This is similar to the multicellular phase field model case (Section 4.3.3), but it is considered to be difficult to make a fair comparison with experimental data because of the shaggy interface shape and the difficulty to introduce elasticity of cells. More models have been proposed to simulate multicellular behavior. Many of them are apparently based on the similar description principle, i.e. the decomposition of a cellular body or surface into a meshwork of small mechanical elements or the continuous description, which in practice also is treated as a descritized meshwork in numerical simulations. [footnote: but not limited to. For example, Ref. [START_REF] Camacho-Gómez | A 3D multi-agent-based model for lumen morphogenesis: the role of the biophysical properties of the extracellular matrix[END_REF] applies the agent-based model approach in which each cell is assumed to be a single non-deformable sphere.] For examples, an earlier review article [START_REF] Fletcher | Mechanocellular models of epithelial morphogenesis[END_REF] explained the Brodland's finite element model, in which a single cell in an epithelial monolayer is divided into a few discrete finite elements [START_REF] Chen | Cell-Level Finite Element Studies of Viscous Cells in Planar Aggregates[END_REF][START_REF] Brodland | The Mechanics of Heterotypic Cell Aggregates: Insights From Computer Simulations[END_REF], the immersed boundary model (IBM), in which the dynamics of cell boundaries is immersed in the surrounding fluid mimicking cytoplasm and extracellular matrix [START_REF] Peskin | The immersed boundary method[END_REF][START_REF] Rejniak | An immersed boundary framework for modelling the growth of individual cells: An application to the early tumour development[END_REF][START_REF] Fogelson | Immersed-boundary-type models of intravascular platelet aggregation[END_REF][START_REF] Liu | Rheology of red blood cell aggregation by computer simulation[END_REF][START_REF] Biteau | Maintaining tissue homeostasis: dynamic control of somatic stem cell activity[END_REF][START_REF] Cooper | Numerical Analysis of the Immersed Boundary Method for Cell-Based Simulation[END_REF], and the subcellular element model (SEM), which represents the cell body by the meshwork of discrete elements each of which have viscoelastic properties [START_REF] Newman | Modeling Multicellular Systems Using Subcellular Elements[END_REF][START_REF] Sandersius | Modeling cell rheology with the Subcellular Element Model[END_REF], etc. See Ref. [START_REF] Fletcher | Mechanocellular models of epithelial morphogenesis[END_REF] for more descriptions of these models and their strengths and limitations. Furthermore, in a more general sense, SEM might be regarded as a class of models which describe the cellular body or intracellular structures of each cell, such as a cell surface (membrane and cortex) and a nucleus, by a certain meshwork of many nodes. This might include the model relying on the interconnected particles representing subcellular structures [START_REF] Dokmegang | Computational modelling unveils how epiblast remodelling and positioning rely on trophectoderm morphogenesis during mouse implantation[END_REF][START_REF] Van Liedekerke | Quantitative modeling identifies critical cell mechanics driving bile duct lumen formation[END_REF] and the continuous descriptions of cell surfaces implemented by triangularmesh discretizations or some finite-element discretization [START_REF] Bächer | A Three-Dimensional Numerical Model of an Active Cell Cortex in the Viscous Limit[END_REF][START_REF] Borja Da Rocha | A viscous active shell theory of the cell cortex[END_REF][START_REF] Torres-Sánchez | Interacting active surfaces: a model for threedimensional cell aggregates[END_REF].

We found that, among them, indeed IBM and SEM have been applied to luminogenesis simulations up to present. IBM was extensively used for this purpose [START_REF] Rejniak | A Computational Study of the Development of Epithelial Acini: I. Sufficient Conditions for the Formation of a Hollow Structure[END_REF][START_REF] Rejniak | A Computational Study of the Development of Epithelial Acini: II. Necessary Conditions for Structure and Lumen Stability[END_REF][START_REF] Rejniak | Computational investigation of intrinsic and extrinsic mechanisms underlying the formation of carcinoma[END_REF][START_REF] Rejniak | Homeostatic Imbalance in Epithelial Ducts and Its Role in Carcinogenesis[END_REF]; e.g., using IBM in two dimension and introducing various cell state changes, like proliferating to polarized (not proliferating) states and apoptotic event and so forth, the conditions for normal acinar growth and maintenance of the hollow lumen generated by apoptosis were investigated [START_REF] Rejniak | A Computational Study of the Development of Epithelial Acini: I. Sufficient Conditions for the Formation of a Hollow Structure[END_REF][START_REF] Rejniak | A Computational Study of the Development of Epithelial Acini: II. Necessary Conditions for Structure and Lumen Stability[END_REF]. SEM can also simulate luminogenesis, by assuming multiple cells and allowing the empty space without any cells but with a specific nature corresponding to luminal space [START_REF] Dokmegang | Computational modelling unveils how epiblast remodelling and positioning rely on trophectoderm morphogenesis during mouse implantation[END_REF][START_REF] Van Liedekerke | Quantitative modeling identifies critical cell mechanics driving bile duct lumen formation[END_REF]. For example, Ref. [START_REF] Dokmegang | Computational modelling unveils how epiblast remodelling and positioning rely on trophectoderm morphogenesis during mouse implantation[END_REF] studied if repulsive force at the apical surface of the rosette-shaped epiblast drives luminogenesis.

Strengths/Weakness

Compared with the Vertex model case (Section 4.2.2), this framework allows more detailed cell shape and its dynamics. Still, this framework can take into account the forces which each element experiences straightforwardly similarly to the Vertex model. However, this model requires a large computational cost (depending on the fineness of the meshwork). In addition, unlike multicellular phase field model approach (Section 4.3.2), this case requires extra spatial variables to consider the fluid flow and/or molecular density dynamics out of the meshworked objects.

Experimental readouts for comparison

One may follow the 2D or 3D contours of cells and lumens in time and in space and measure the mechanical parameters similarly to the Vertex model case (Section 4.2.3), but the detailed shapes can be also compared.

Conclusion

In this article, we reported experimental features and theoretical model approaches applicable to study luminogenesis. We presented several approaches using the simple morphology model (Section 4.1), the vertex model (Section 4.2, Fig. 2), the phase-field model (Section 4.3, Fig. 3) and the cellular Potts model (Section 4.4, Fig4) and some meshwork models (Section 4.5).

A phenomenon can often be reproduced by various modeling frameworks. The selection of frameworks represents still a challenge. Although it is difficult to make a clear guide of how to select the model, there are several practices which we should follow to choose a better model. Firstly, as we saw in Section 4 and Table 1, each model framework has their own strength and weakness. We should select the model which can take advantage of its strength and avoid the weakness. Secondly, the number of variables is important. While the model with more variables can typically describe a wider variety of morphology, such a model may require more computational effort. It is better to use the model with a minimal number of the variables. For example, here let us compare the four model approaches focused on in this article. If the possible morphologies which the target system may take are expected well and all of them are simple, we may apply the simple morphology model approach. If not, we may consider the other model approaches. The phase-field and cellular Potts models can describe the curved interface and simulate the phenomena for which the curved interface is essential (like the case with a spherical lumen sandwiched by only two cells) whereas the vertex model has the advantage of less redundancy as long as all the cell interfaces are expected to be flat and the gap between the interfaces of two neighboring cells is negligible. If the vertex model is not applicable, e.g., because one wants to address the curved interface shape, the phase-field or cellular Potts models may be considered. With the same set of free energy, both the models are expected to result in the similar predictions, but their applicable range of the problem and computational efficiency are different. For example, since cellular Potts model cannot avoid interface fluctuation, a phase-field model may be chosen if one wants to add other variables on the interface, like density distribution of some molecules. Listing up the strengths/weaknesses and the necessary set of variables and comparing numerical simulation efficiencies among the models, as performed in Ref. [START_REF] Osborne | Comparing individual-based approaches to modelling the self-organization of multicellular tissues[END_REF], may be a systematic strategy for the selection.

Practically, the availability of the computer simulation tools may be another axis for the model selection.

There are several open software tools to simulate the phase-field model [START_REF] Hong | Open-Sourcing Phase-Field Simulations for Accelerating Energy Materials Design and Optimization[END_REF]. Open software tools for CPM and the vertex model are also available (CompuCell3D [START_REF] Swat | Multi-scale modeling of tissues using CompuCell3D[END_REF] and Tyssue [START_REF] Theis | Tyssue: an epithelium simulation library[END_REF], respectively). For the other computational frameworks provided in Sect. 4.5, including IBM and SEM, the open software called Chaste can be applicable [START_REF] Mirams | Chaste: An Open Source C++ Library for Computational Physiology and Biology[END_REF].

We close this article after briefly describing the key future challenges in this field. One important challenge is to establish how to incorporate the cellular autonomous regulations of biophysical parameters. Such autonomous regulation may be implemented by the support of other variables describing the density distribution of some molecules. Among them, actomyosin cytoskeleton, which dominates surface mechanics and is often represented by the active gel model (Section. 3.1.1), is now one of the key ingredients which several theory groups are challenging to incorporate into various frameworks. Another challenge is to make a clear guide on how to bridge theoretical parameters and the experimental methods of perturbation, such as changes in luminal pressure, pumping ability, cell-cell adhesion and actomyosin activity.

Figure 1 .

 1 Figure 1. Shapes of lumens. Schematic of (a) a spherical central lumen surrounded by a monolayer, (b) a lumen surrounded by asymmetric cell layers, (c) multi-lumen, (d) unbranched tube, and (e) branched tube. Red lines represent the apical side.

  Fig. 2 Computer simulation of mammalian blastocyst morphology. Drawing (top left) shows cross-sections in xz-and yz-planes of simulated aggregates. (a) An example of calculation process to a stable state viewed from cross-sections of the yz-plane. Numbers indicate the time point of the simulation (t). A vertex (a black circle in the sample at t=0) is replaced by a small tetrahedron and is enlarged until its volume reaches half the initial total volume (at𝑡 = 500). The blastocyst axis keeps changing during 𝑡 = 500 -2000 a.u. until it is localized and stabilized at one end of the long axis of the ellipsoidal ZP (𝑡 = 2000 𝑎. 𝑢.), when it no longer migrates (𝑡 = 2000 -3500 𝑎. 𝑢.). (b) The simulated blastocyst at 𝑡 = 2000 in cross-sectional views of xz-and yzplanes. (c) A stereoscopic view of the blastocyst at 𝑡 = 2000 𝑎. 𝑢. , in which the ZP and some of the TE cells are removed for internal view. (Reused from excerpt on page 1412 and Fig. 4 in [10]).

Figure 4 . 4 . 5

 445 Figure 4. (a) Schematics of configuration in CPM. (b) (c) Cell sorting simulation with CPM. (b) Initial random aggregate. Source [133]. (c) After 10000 MCS.(c)(A-D) Cyst formation simulations starting from a single isolated cell. (A) Snapshots of reference cyst formation corresponding to 12 (1), 19.1 (2), 47.6 (3), and 142.9 (4) h. (B) small-cyst phenotype, (C) large-cyst phenotype, and (D) complex multiple-cyst phenotype. Lumen shown in blue, basal compartment shown in darkred; all other cell compartments were made transparent to enhance the visibility of basal compartment and lumen. Source[START_REF] Belmonte | Virtual-tissue computer simulations define the roles of cell adhesion and proliferation in the onset of kidney cystic disease[END_REF].

Table

  Table serves as the basis for Section 4, where we describe the simple morphology model approach, free energy-based modeling approaches (those with cellular vertex model, multicellular phase-field model and the cellular Potts model) and other approaches.

Acknowledgements

We thank C. J. Chan for a fruitful discussion and T. Hirashima and J. A. Glazier for helpful comments to our manuscript. This work was supported by Mechanobiology Institute (to TH) and HFSP Program Grant "Coupling of cell polarization and differentiation in organoids" (to DR, MS).

Conflicts of Interest

The authors have declared that no conflict of interest exists