{'fraction_non_alphanumeric': 0.030875795731792815, 'fraction_numerical': 0.007094983033736224, 'mean_word_length': 5.207136640557006, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Autonomous epithelial folding induced by an intracellular mechano–polarity feedback loop

Epithelial tissues form folded structures during embryonic development and organogenesis. Whereas substantial efforts have been devoted to identifying mechanical and biochemical mechanisms that induce folding, whether and how their interplay synergistically shapes epithelial folds remains poorly understood. Here we propose a mechano–biochemical model for dorsal fold formation in the early Drosophila embryo, an epithelial folding event induced by shifts of cell polarity. Based on experimentally observed apical domain homeostasis, we couple cell mechanics to polarity and find that mechanical changes following the initial polarity shifts alter cell geometry, which in turn influences the reaction-diffusion of polarity proteins, thus forming a feedback loop between cell mechanics and polarity. This model can induce spontaneous fold formation in silico, recapitulate polarity and shape changes observed in vivo, and confer robustness to tissue shape change against small fluctuations in mechanics and polarity. These findings reveal emergent properties of a developing epithelium under control of intracellular mechano–polarity coupling.


I. INTRODUCTION
Epithelial tissues cover body surfaces and line internal organs and cavities in many metazoan animals. During development, epithelial cells undergo shape changes that orchestrate tissue-scale deformation such as convergent extension, branching and folding [11]. Epithelial cell shape changes typically require mechanical forces that are generated via biochemical processes that occur within the cells themselves, while these biochemical processes are governed by reaction, diffusion and transport of their molecular constituents, which may depend on the boundary conditions, including the cell shape [20]. Cell shape changes can thus influence the biochemical processes that induce these shape changes in the first place, thereby integrating mechanical forces and their underlying biochemical processes into a feedback loop [14]. Such a mechano-biochemical feedback loop has been proposed to play a pivotal role in coordinating cell morphogenesis and tissue reorganization during normal development as well as contributing to pathological conditions that involve dysregulation of tissue growth and cell movement [2,10,32].
Epithelial cells exhibit intrinsic polarity that allows them to perform specialized functions.
Epithelial polarity can be established, and needs to be maintained, by a reaction-diffusion type mechanism that governs the non-uniform distribution of polarity regulators across the cell surface [5]. The partitioning-defective (PAR) polarity proteins display polarized distribution along the apical-basal axis of the cells to regulate the activity of spatially localized structural and signaling proteins, including those that produce, modulate and transmit mechanical forces that control epithelial deformation [24,37]. Thus, modulation of the polarity system may affect cell mechanical properties, inducing cell shape changes that alter tissue architecture.
One such polarity-dependent deformation process is the formation of the epithelial folds that are generated on the dorsal side of the early fruit fly (Drosophila melanogaster) embryo, called dorsal fold formation [38,40]. Distinct from the canonical epithelial folding events that require the motor protein myosin that contracts the actin cytoskeleton, Drosophila dorsal fold formation is initiated by the remodeling of the PAR-dependent epithelial polarity without changes of myosin activity [40]. Although experimental evidence suggests that differential modulation of PAR-dependent polarity between cells is essential for dorsal fold formation and that microtubule-based forces underlie polarity-dependent cell shape changes [38,40], the exact mechanism by which the polarity remodeling process instigates the mechanical alteration to induce folding has not been resolved. Here we postulated that during dorsal fold formation cell shape changes impact the remodeling of epithelial polarity, thereby integrating cell deformation and polarity evolution into a feedback loop. We tested this hypothesis using mathematical modeling and verified our theoretical predictions via quantitative in vivo imaging of the Drosophila embryo. We show that such mechano-polarity feedback allows epithelial cells to autonomously change their shapes, leading to self-organized formation of epithelial folds. Our modeling results thus support the previously proposed hypothesis that external forces from the surroundings are not required for dorsal fold formation [40]. This model allows us to uncover self-organizing principles of tissue reorganization emerging from intracellular mechano-polarity coupling.

II. MATHEMATICAL MODELING
To simulate the polarity-dependent dorsal fold formation in Drosophila, we constructed a three-dimensional (3D) model which consists of biochemical and mechanical modules as described below.
A. The biochemical module for the formation of epithelial apical-basal polarity The plasma membrane of the epithelial cells in the Drosophila blastoderm embryo are polarized along the apical-basal axis (Fig. 1A). The key polarity regulators, atypical protein kinase C (aPKC), Par-1, and Bazooka, are each localized to, and thereby respectively define, the apical and basolateral domains as well as the interface domain between these two [30].
To model epithelial polarity formation, we constructed a minimal reaction-diffusion model following a previous non-epithelial model developed for the one-cell stage C. elegans embryo [12]. Our three-element network model ("A-B-P network"; Fig. 1B) contains the apical and basal-lateral proteins, aPKC and Par-1 (denoted by A and P , respectively), as in the C. elegans model. In addition, we added a third and epithelium-specific element, Bazooka (denoted by B), that specifies the subapical region where cell-cell adherens junctions are positioned [16]. Since the reactions responsible for the polarity formation take place when the polarity regulators associate with the plasma membrane, we assumed that the polarity system can be approximated by a surface reaction that is coupled to a cytoplasmic pool of polarity regulators. Given that the membrane distribution of polarity regulators is approximately rotationally symmetric with respect to the apical-basal axis, we simplified the surface reaction by considering the spatial distribution of regulator concentrations along the apical-basal axis in one dimension.
The evolution of concentrations are thus formulated as follows, where A(s, t), P (s, t), and B(s, t) denote the membrane concentration of aPKC, Par-1, and Bazooka, respectively, at time t and membrane location s along the apical-basal axis (0 ≤ s ≤ L, see Fig. 1B). D A , D P , and D B are respective diffusion coefficients that govern the diffusion of these polarity regulators on the membrane surface.
The second terms in Eq. (1) describe the membrane association or dissociation reaction kinetics of the polarity regulators. These kinetics depend on the molecular interactions among the regulators themselves. Based on previously established polarity protein interactions in various epithelial systems [36], we consider phosphorylation-mediated mutual inhibition of the membrane association between aPKC and Par-1 [8,21], the inhibition of membrane association of Bazooka by both aPKC [7,39] and Par-1 [3], and the membrane recruitment of aPKC by Bazooka [25] (Fig. 1B, left). The reaction terms thus read, wherein k on/of f,A , k on/of f,P and k on/of f,B are the membrane association (or dissociation) rates per unit concentrations in the cytosol (or on the membrane), k XY with X, Y ∈ {A, P, B} is the membrane dissociation rates resulted from the mutual antagonism between regulators X and Y , k AB is the membrane association rate of A due to the presence of B, whereas α, β, γ, δ, and are the parameters of these molecular interactions.
We assumed that the regulators diffuse much faster in the cytosol than on the cell membrane. The cytosolic concentration is thus taken to be uniform in space and can be calculated from the principle of mass conservation, as follows: where ρ total is the total concentration (total number of molecules per unit volume of the cell) with ρ ∈ {A, P, B}, ψ is the surface-to-volume conversion factor (see below for details), whileρ is the average concentration across the membrane, calculated as where L = a + b + is the total length of the membrane region in which regulators are distributed (see Fig. 1B, right), and a, b, and are the respective apical, basal, and lateral lengths of the cell (see Fig. 1A).
B. The dependence of polarity formation on cell geometry The epithelial polarity regulators either reside in the cytosol or are associated with the membrane surface. As indicated in Eq. (3), for a given polarity regulator, the probability of residing in the cytosol versus associating with the membrane depends on the surface-tovolume ratio ψ of the cell, which is calculated as follows: where V cell is the cell volume, and A cell = A a + A b + A is the total surface area of a cell with A a the apical area, A b the basal area, and A the lateral area. Moreover, the half cell perimeter L also affects epithelial membrane polarization as it determines the average membrane concentration of the polarity regulators, as revealed by Eq. (4).

C. The mechanical module for cell shape changes in three dimensions
To simulate 3D epithelial shape formation, we developed a vertex-based mechanical model for a monolayered epithelial tissue consisting of N cells. In the vertex model, the full 3D shape of a cell is described by a hexagonal prism with edges and vertices [1] (Fig. 1A).
Forces acting on the vertices are derived from the derivative of the potential function U , which describes the mechanical properties of the entire tissue. Such forces are balanced by the surrounding viscous drag in an environment of low Reynolds numbers. Thus, the equation of motion for cell shape changes reads, where r j is the position coordinate of vertex j, and η is the friction coefficient. A forcebalanced tissue morphology is achieved when U reaches a minimum.
The potential function that describes the mechanics of the entire tissue takes into account energy contribution from all cell surfaces and cell volumes. Specifically, we considered energy resulting from the surface tensions on the basal and lateral sides of the cells, the volume elasticity that maintains the cell volume during cell shape changes, and the elastic energy of the apical surface given that cytoskeletons (e.g., actomyosin and microtubule networks) are enriched on the apical sides of the cells [27]. The potential function is thus formulated as follows, The first term describes the apical elasticity with elastic modulus K a i , apical perimeterP a i , and preferred perimeterP o i for the i-th cell. The second and third terms describe the tension along the basal and lateral surfaces, respectively, with surface tension λ b i and area A b i for the basal surface, while surface tension λ k and area A k for the lateral surface with an index k that runs across all lateral cell surfaces. The fourth term defines the volume elasticity with elastic modulus K V i , cell volume V i , and preferred volume V o i for the i-th cell. The volume elasticity likely arises from the hydrostatic pressure of the cytoplasm and the intracellular cytoskeletons that resist changes of the cell volume. Note that perimeterP a i , area A b i and A k , as well as volume V i are calculated from the vertex positions, whose dynamics follow Eq. (6).
By using qualitatively distinct mechanical terms, we took into account that the apical and basal surfaces of a polarized epithelial cell possess distinct properties. In particular, in addition to using a tension term for the apical perimeter as was used previously [15], we further considered the elastic properties of the cytoskeletons along the apical perimeter.
Such apical elasticity helps to stabilize cell shape under large deformation.  [38]. It was thus proposed that as the Bazooka position continues to shift towards the basal side of the cells due to Par-1 downregulation, the mechanism that ensures a constant aPKC domain size during the basal shifts of Bazooka enables the reduction of cell height [38]. Because the mechanism underlying the aPKC domain homeostasis remains unclear, we modeled the mechanistic link between cell polarity and mechanics using a phenomenological approach. Specifically, we hypothesized a functional term that defines the polarity-dependent modulation of cell mechanics based on the homeostatic relaxation of the aPKC domain, as follows: where  The Drosophila cellular blastoderm forms its primary epithelium via cellularization. During cellularization, the microtubule-dependent, minus end-directed transport of Bazooka, and its eventual accumulation on the apical side of the cell, function as the key symmetry breaking event that polarizes the cells [18]. To simulate the observed apical accumulation of Bazooka, we used the Heaviside step function for the initial distribution of B, while keeping A and P uniform along the apical-basal axis ( Fig. 2A, dashed lines). We then numerically solved the spatio-temporal distribution of polarity regulators along the cell perimeter with periodic boundary conditions according to Eq. (1) and parameter values given by Table I. As shown in Fig. 2A, the distribution of A, B and P evolves toward a steady-state polarized profile, indicating that the A-B-P network recapitulates epithelial polarization in the early Drosophila embryo. The simulated steady-state profile resembles that observed in the early Drosophila embryo, in the sense that while aPKC (A) and Par-1 (P ) form the polarized apical and basal-lateral domains, respectively, Bazooka (B) is localized as a sharp peak between these two domains.
We next asked whether the A-B-P network model recapitulates features of polarity remodeling previously observed in the experiments. On the dorsal side of cellularizing embryo, down-regulation of Par-1 in a specific population of cells, called the initiating cells, results in basal shifts of the positions of Bazooka and adherens junctions [40]. To simulate this, we introduced a 10% reduction in the Par-1 total concentration to the above steady-state polarity profile in Fig. 2A, approximating the observed decrease of Par-1 in the initiating cells prior to the end of cellularization [40]. Figure 2B shows that the peak position of B exhibits a basal shift of approximately 10% cell height, comparable to the experimental observations [40]. Given that Bazooka co-localizes with adherens junctions in the early Drosophila embryo [17,18], the peak position of B corresponds to the junctional position.
This result thus supports the hypothesis that basal shifts of adherens junctions could result from biochemical remodeling of polarity regulators, i.e. the down-regulation of Par-1 [40].
Previous work also showed that overexpression of Bazooka is sufficient to induce basal shifts of junctions that cause the formation of ectopic dorsal folds [40]. Motivated by this, we increased the total concentration of B in the steady-state profile shown in Fig. 2A, and found that the peak position of B shifts basally by 10% when B concentration is increased by 80% (Fig. 2C), consistent with experimental observations reported previously [40]. The result also raises the question of whether the system exhibits quantitatively different response to Par-1 down-regulation and Bazooka overexpression. To address this question, we explored a wide range of P total and B total changes. For the same magnitude of variation, we observed a greater shift in the peak position of B in response to changes in P total than B total ( Fig. 2D; see the crosses for comparison of a 20% variation). This result predicts that the position of adherens junctions is more sensitive to changes in the concentration of Par-1 than Bazooka.
To test this theoretical prediction, we re-analyzed a previously published dataset on the extent of junctional shifts in response to Par-1 and Bazooka changes [40]. For Par-1 down-regulation, we compared wild-type (WT) neighboring cells (i.e., high Par-1) with WT initiating cells (i.e., low Par-1), whereas for Bazooka overexpression (BazOE), we compared WT neighboring cells (i.e., low Bazooka) with BazOE neighboring cells (i.e., high Bazooka) (see Materials and methods for details). Figure   Our theory thus predicts two distinct types of basal polarity shifts: (1) scaling, in which the relative junctional position remains unchanged during polarity redistribution (e.g., see  our experimental data validate the theoretical prediction on the scaling effect. Note also that for the initiating cells, the initial basal shifts of junctions have been shown to occur without cell shortening [40]. Since ψ is expected to be near constant during these initial basal shifts, these shifts must result from biochemically induced downregulation of Par-1, but not the ψ-dependent remodeling effect. Following such biochemically initiated basal shifts, the initiating cells begin to shorten into a frustum shape, during which it is possible that the increased ψ resulting from the shortening (as will be shown later in Fig. 4A) further enhances the basal shifts. To study how this mechano-polarity feedback loop controls epithelial cell shape, we numerically computed both polarity distribution and cell shape based on Eqs. (1) and (6) with their coupling terms given by Eqs. (5) and (8). Modulation was applied to cells that initially have a stable shape and steady-state polarity distribution. To simulate the condition of Par-1 down-regulation (or Bazooka overexpression), we decreased (or increased) total concentration of P (or B). This induces an initial expansion of the A domain, which subsequently shrinks toward its homeostatic size (Fig. 4A, top). During such homeostatic relaxation, surface-to-volume ratio ψ increases, while half cell perimeter L decreases (Fig. 4A, bottom). Ultimately, cells transform from a tall columnar shape to a short frustum shape, resembling the shape changes that the initiating cells exhibit during Drosophila dorsal fold initiation [40]. Since cells in the Drosophila dorsal fold system initially have columnar shapes, during mechano-polarity feedback coupling, homeostatic relaxation of the apical domain can be most effectively achieved by the alteration of cell height.
We next examined the tissue-scale morphological changes following alterations of cell polarity in this mechano-polarity feedback system. Figure 4B shows that both a decrease in the total concentration of Par-1 (P total ) and an increase in that of Bazooka (B total ) within a localized tissue area induce out-of-plane deformation, recapitulating both normal and ectopic dorsal fold formation [40]. Thus, a local trigger of polarity remodeling can induce spontaneous epithelial folding via mechano-polarity feedback without external forces. These results provide a theoretical foundation for the previously proposed hypothesis that dorsal fold formation is initiated cell-autonomously, whereas external forces from the surrounding tissues are dispensable [40]. This result suggests that the mechano-polarity feedback loop endows cells with a regime within which small mechanical fluctuations, presumably arising stochastically, lead to minuscule, non-consequential morphological changes, whereas active, instructive mechanical alterations exceeding a certain level are required to induce substantial morphological changes, implying that the mechano-chemical coupling could generate a threshold-like response, such that only perturbations above the threshold value can induce sufficiently large deformation. Thus, a system with mechano-polarity coupling could be more resistant to the noise-induced ectopic initiation of morphogenesis. Note that a sharp cell/tissue deformation response was also observed following the perturbation of Par-1 total concentration (Fig. S2). This result suggests that low-amplitude fluctuations of Par-1 concentration at the steady-state do not lead to epithelial fold initiation, whereas genetically programmed down-regulation of Par-1 that provides a sufficiently strong perturbation would be required to trigger epithelial folding.
Notably, as revealed in Fig. 4C (left and middle panels), whereas cell shape change is smaller with mechano-polarity coupling (red line) than without (blue line), for sufficiently strong perturbations, the coupled system generates greater tissue deformation than the system that lacks coupling. This result can be understood as the consequence of the subsequent modulation of basal-lateral mechanics following the initial perturbation of apical mechanics due to mechano-polarity coupling, given that apical modulation induces only a shallow tissue indentation, whereas basal-lateral modulation generates a deeper fold (Fig. S3A). Thus, mechano-polarity coupling allows modulation of mechanics at membrane surfaces that are not the primary target of the initial mechanical stimuli and may promote the efficiency of tissue deformation.
Tissue deformation driven by mechano-polarity feedback is robust against stochasticity Lastly, we sought to determine whether mechanical fluctuations render the self-organized epithelial deformation stochastic as biological systems are intrinsically noisy. We added noise to the mechano-polarity coupling term c i to simulate fluctuations in the rate of conversion from the polarity signal to the mechanical conditions at the cell membrane. The noise is given as follows: where c o i is the deterministic portion of the polarity-to-mechanics conversion rate, while ξ i (t) is the Gaussian white noise that has a zero mean and a standard deviation σ. The value of σ defines the extent of stochasticity in how cell mechanics q i decode the polarity signal (see Eq. (8)).
To evaluate robustness at the cellular level, we quantified the geometric features, ψ and L, of the final cell shape, and the time, τ , it takes to achieve apical domain homeostasis.
We observed that both ψ and L are nearly identical across a wide range of noise strengths σ (Fig. 4D, left and middle), although the variability of relaxation time τ increases with stronger noise (Fig. 4D, right). At the tissue level, we also observed low levels of variability in the epithelial fold depths (Fig. 4B, right). These data indicate that mechano-polarity coupling ensures robust cell shape change and tissue folding despite noise.

IV. DISCUSSION
This work proposes a theoretical model for dorsal fold formation in the Drosophila embryo, an epithelial folding event that is induced by changes in epithelial cell polarity. Our model recapitulates key features of polarity remodeling and links changes in cell polarity to mechanical modulation by phenomenologically mimicking the experimentally observed apical domain homeostasis. Such mechanical modulation induces cell shape changes, altering system geometry on which the reaction-diffusion process of the polarity network depends, thereby forming a feedback loop of mechano-polarity coupling. Simulation with this feedback model indicates that genetically programmed Par-1 down-regulation can trigger spontaneous cell shape changes, leading to self-organized tissue fold formation. Such an intracellular mechano-polarity crosstalk prevents epithelial tissues from undergoing shape changes under small fluctuations in cell mechanics and polarity, which could occur owing to the intrinsic noise that may arise during development. The intrinsic property of robustness against stochastic fluctuations can thus ensure that the initiation of autonomous tissue deformation is spatio-temporally precise and tightly controlled.
The current model contains two key hypotheses for which mechanistic insights remain to be explored. First, we hypothesized that surface-to-volume ratio ψ affects the balance between the cytosolic and membrane pools of the polarity regulators. Previous studies have reported that membrane localization of the PAR proteins depends on their association with the specific lipid moieties in the cell membrane. For instance, Par-1 binds to phosphatidylserine [29], whereas Bazooka mainly associates with phosphatidylinositol [13,26]. Assuming that in the membrane region in which the polarity proteins localize, the percentage presence for each specific type of lipid is near constant, the number of these "membrane receptors" for a given polarity protein would likely scale proportionally with the membrane surface area. The constant membrane lipid composition could thus be the potential mechanism that underlies ψ-dependent cytosol-membrane exchange of the polarity proteins. Note that variation of ψ has also been proposed to create biochemically distinct states (e.g., with regard to the phosphorylation status of a protein) within a cell, which influences cell division timing or the scaling of internal organelles to the overall cell sizes [19,28]. We speculate that alterations of surface-to-volume ratio can be a common mechanism that modulates intracellular reactions, thereby regulating cellular behaviors.
The second hypothesis is that aPKC domain homeostasis underlies mechanical modulation in response to polarity change. As mentioned in the Mathematical Modeling section, a recent study on the initiation of dorsal fold formation reveals that the apical domain of the initiating cells maintains a relatively constant size as the cells undergo shortening to induce fold formation [38]. The study also reveals potential connections between the polarity proteins and the mechanical elements involved in the homeostasis. Specifically, the microtubule minus-end binding protein Patronin was shown to become associated with the cell cortex underneath the apical surface in an aPKC-dependent manner. The cortex-associated Patronin organizes an apical microtubule network that appears to undergo rapid turnover, which depends on Patronin-dependent recruitment of the microtubule-severing enzyme Katanin that likely destabilizes the microtubules. Given that Patronin distribution in the apical domain responds to downregulation of Par-1 during dorsal fold initiation, and the forces that the apical microtubule network exerts on the apical cortex may scale with the cortical density of Patronin, it was proposed that this rapid turnover, coupled with the polarity-dependent cortical localization of Patronin and the apical microtubule network that associates with it, underlies the apical size homeostasis [38]. Patronin could thus be at the nexus of this mechano-polarity coupling in the dorsal fold system. It is worth noting that in the Patronin RNAi embryos the shape of the dorsal epithelial cells become extremely irregular [38]. This observation is consistent with the prediction of our model that without mechano-polarity coupling the cell shapes are more sensitive to fluctuations in mechanics and cell polarity.
We propose that the reduced Patronin activity impairs mechano-polarity coupling such that cells in the epithelium become more prone to the ambient noise in the system.
Although we constructed the current model by closely following properties specifically linked to the Drosophila dorsal fold system, it is likely that the systems-level characteristics of the mechano-biochemical crosstalk, including the autonomy of tissue deformation and the robustness against noise, transcend context-specific biochemical mechanisms of cell polarity and mechanical regulation of cell shapes. In support of this, self-organized behaviors orchestrating cell shapes have been observed in previous models that consider bidirectional feedback between cell mechanics and tissue-scale polarity cues [35] or morphogen patterns [4,6,31]. Moreover, it has been proposed that feedback between intracellular mechanics and extracellular morphogen concentrations leads to robust formation of morphogen concentration profile [4,6], thereby ensuring proper tissue architecture. We envision that self-organization and robustness are general properties of mechano-chemical feedback control during tissue morphogenesis and remodeling. Our modeling framework could potentially serve as a generalizable platform that can be modified to model mechanical-biochemical coupling systems whose mechanistic details are distinct from the dorsal fold system.
The vertex-based modeling represents a powerful framework for simulating tissue morphogenesis in various biological contexts [1,9]. The vertex models have recently been developed to further take into account the regulatory mechanisms of cell mechanics. In contrast to previous 3D models that consider regulation of cell shape changes by a global, tissue-scale signaling pattern [6,31], cell deformation in our model is regulated by a local signal resulting from polarity pattern formation within a cell, which is in turn affected by cell deformation. Thus, our model provides insights into regulation of cell mechanics that arises at the cellular scale, while adjustable to the resultant change of cell shape, and reveals novel tissue properties emerging from such intracellular crosstalk.

MATERIALS AND METHODS
A. Quantitation of junctional positions in the dorsal epithelium of fixed wild-type or Bazooka overexpression embryos Previously published datasets [40] were re-analyzed for the comparison between the effects of Par-1 down-regulation and Bazooka overexpression on the position of adherens junctions (represented as the Bazooka peak position). Par-1 and Bazooka concentration was quantified as previously described [40]. In brief, wild-type embryos and embryos overexpressing Bazooka were fixed, stained with Par-1, Bazooka and Neurotactin (for cell membrane visualization) antibodies, and imaged in 3D on a Leica SP5 scanning confocal system with a 63x glycerol immersion objective. The resultant images were processed, segmented and quantified using the EDGE4D software [23]. Bazooka concentration was quantified as the per voxel intensity within the segmented Bazooka junctional localization. Par-1 concentration was the per voxel intensity within a two-voxel distance from the 3D reconstructed membrane Embryos expressing E-Cadherin-3xGFP [33], Myosin-GFP [34], and Gap43-mCherry [22] were imaged on an Olympus Fluoview FVMPE-RS multi-photon scanning system with a 20X water immersion objective. The resultant images were processed, segmented and quantified using a custom-made ImageJ macro. Because of the slender columnar cell shape, junctional position (i.e., Bazooka peak position) was estimated as the mean vertical distance between the pixels with local maxima of E-Cadherin-3xGFP intensity (i.e., the green circles in Fig. 3D), given that Bazooka and E-Cadherin positions are identical [40], and the basal margin of the segmented autofluorescent intensity of the vitelline membrane (i.e., the red circles in Fig. 3D), defined as the cell apices. Cell height, which approximates L in the theory, was estimated as the mean vertical distance between the bases of the cells, defined as the basal margin of the segmented Myosin-GFP intensity (i.e., the blue circles in Fig. 3D Appendix A: Linking cell surface mechanics to cell and tissue shape via a 3D vertex model Here, we sought to understand how epithelial tissue shape is determined by the mechanical properties of its constituent cells using the 3D vertex model described in Mathematical Modeling section. For a given set of cell mechanics, cell and tissue shape are derived from numerically computed vertex positions according to Eq. (6). The mechanically stable tissue shape is then obtained by minimizing the potential function U , as described in Eq. (7). Using this 3D vertex model, we studied the effect of locally induced modulation of cell mechanics on force-balanced tissue shape (Fig. S3A). For the shape stability analysis, see Appendix B.
Apical constriction of cells at the center of the tissue, highlighted in yellow, produces narrow apices, elongated lateral surface and widened basal membrane, and the tissue forms a shallow indentation (Fig. S3A, upper-left on the left panel). In contrast, basal expansion and lateral shrinkage induce laterally shortened and basally widened cells with minimal changes on the apical side, leading to deep, out-of-plane deformation (Fig. S3A, upper-left on the right panel). These results indicate that quantitative features of cell shape change and tissue deformation depend on the surfaces at which the mechanical modulation is applied, confirming our previous findings using a 2D model [41].
We next analyzed cell deformation in three conventionally categorized cell types, the columnar, cuboidal, and squamous cells (Fig. S3B), to better understand how cells respond to applied mechanical modulation. We found that in the columnar ( ) and cuboidal cells ( ), changes in the apical and basal lengths, a and b, show negative correlation when apical mechanicsP o (Fig. S3B, left column) is altered. In contrast, changes in a and b are positively correlated in the squamous cells ( ). This result indicates that different cell types show distinct deformation pattern in response to a given mechanical modulation.
Furthermore, for a given cell type, cell deformation pattern depends on which cell surface is mechanically modulated. For example, for the columnar cells ( ), a-b and a-h are negatively correlated when apical mechanics is altered (Fig. S3B, left column), whereas ah correlation becomes positive when basal cell mechanics λ b (Fig. S3B, middle column) is modulated, or a-b correlation becomes positive due to modulation of lateral mechanics λ (Fig. S3B, right column). Our 3D vertex model shows that different cell types or modes of mechanical modulation give rise to qualitatively distinct patterns of cell shape alteration, confirming our previous 2D analysis [41]. We suggested previously that experimentally derived correlation pattern on cell shape changes can be used to infer the mechanical principles that induce the shape changes [41]. Our current results suggest that this methodology can be extended to a 3D analysis.
Appendix B: Force balance analysis of 3D epithelial cells According to Eq. (7), the potential energy function of an epithelial cell with the shape of a regular hexagonal prism reads, where being the height and radius of the hexagonal prism, respectively (see Fig. S3B, lower right legend). From Eq. (B1), the force balance condition is given by the following relations: where P V = 2(1 − V ) is the pressure resulting from cell volume elasticity. For a compressed cell (i.e., P V > 0), intracellular pressure pushes outward against shrinkage of compressive cell surfaces, and vice versa. To prevent cells from being distorted into frustum shape, we further considered the following condition: Therefore, from Eqs. (B2) and (B3), the condition for maintaining the shape of a hexagonal prism is The stability of the cell shape can be determined from the Hessian matrix M of U single with the following components, Cell shapes are mechanically stable when the trace and the determinant of M are positive, which gives the stable condition, as follows: Note that K a , γ b , and γ are replaced by h, r, andP o according to Eq. (B4).

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.