A quantitative principle to understand 3D cellular connectivity in epithelial tubes

1 A quantitative biophysical principle to explain the 3D cellular 1 connectivity in curved epithelia 2 3 Pedro Gómez-Gálvez1,2,3,4,*,Ώ, Pablo Vicente-Munuera1,2,5Ώ, Samira Anbari6,7Ώ, Antonio Tagua1,2,Ώ, Carmen 4 Gordillo-Vázquez1,2, Jesús A. Andrés-San Román1,2, Daniel Franco-Barranco8,9, Ana M. Palacios1,2, 5 Antonio Velasco1, Carlos Capitán-Agudo1, Clara Grima10, Valentina Annese1,2, Ignacio Arganda- 6 Carreras8,9,11, Rafael Robles10, Alberto Márquez10, Javier Buceta12,*, Luis M. Escudero1,2,* 7 8 1Instituto de Biomedicina de Sevilla (IBiS), Hospital Universitario Virgen del Rocío/CSIC/Universidad de 9 Sevilla and Departamento de Biología Celular, Universidad de Sevilla. 41013 Seville, Spain 10 2Biomedical Network Research Centre on Neurodegenerative Diseases (CIBERNED), Madrid, Spain 11 3Current address: MRC Laboratory of Molecular Biology, Cambridge, UK 12 4Current address: Department of Physiology, Development and Neuroscience, University of 13 Cambridge, UK 14 5Current address: Laboratory for Molecular Cell Biology, University College London. London, UK 15 6Chemical and Biomolecular Engineering Department, Lehigh University. 16 Bethlehem, PA 18018, USA 17 7Current address: Biomedical Engineering Department, Johns Hopkins University. Baltimore, MD 18 21205, USA 19 8Department of Computer Science and Artificial Intelligence, University of the Basque Country 20 (UPV/EHU). San Sebastian, Spain 21 9Donostia International Physics Center (DIPC). San Sebastian, Spain 22 10Departamento de Matemática Aplicada I, Universidad de Sevilla. Seville 41012, Spain 23 11Ikerbasque, Basque Foundation for Science. Bilbao, Spain 24 12Institute for Integrative Systems Biology (I2SysBio), CSIC-UV, Paterna (Valencia) 46980, Spain 25 26 ΏThese authors contributed equally to this work. 27 *Corresponding authors 28 Lead contact: Luis M. Escudero (lmescudero[email protected]s) 29 0DQXVFULSW &OLFNKHUHWRYLHZOLQNHG5HIHUHQFHV 2 SUMMARY 1 Epithelial cell organization and the mechanical stability of tissues are closely related. 2 In this context, it has been recently shown that packing optimization in bended/folded 3 epithelia is achieved by an energy minimization mechanism that leads to a complex 4 cellular shape: the scutoid. Here we focus on the relationship between this shape and 5 the connectivity between cells. We use a combination of computational, experimental, 6 and biophysical approaches to examine how energy drivers affect the three- 7 dimensional (3D) packing of tubular epithelia. We propose an energy-based stochastic 8 model that explain the 3D cellular connectivity. Then, we challenge it by 9 experimentally reducing the cell adhesion. As a result, we observed an increment on 10 the appearance of scutoids that correlated with a decrease of the energy barrier 11 necessary to connect with new cells. We conclude that tubular epithelia satisfy a 12 quantitative biophysical principle, that links tissue geometry and energetics with the 13 average cellular connectivity. 14 15  ƌĞĐŽƌĚ ŽĨ ƚŚŝƐ ƉĂƉĞƌƐ dƌĂŶƐƉĂƌĞŶƚ WĞĞƌ ZĞǀŝĞǁ ƉƌŽĐĞƐƐ ŝƐ ŝŶĐůƵĚĞĚ ŝŶ ƚŚĞ 16 Supplemental Information 17 18 KEYWORDS 19 Tissue/Cellular Biophysics, Computational geometry, Developmental Systems Biology, 20 Mathematical/Biophysical modeling, Bioimage Analysis 21 22 3 INTRODUCTION 1 During the last decades much progress has been achieved in the understanding of 2 the emergence of self-organization in tissues. This problem has been addressed from 3 the viewpoint of energetics considerations (Alt et al., 2017; Canela-Xandri et al., 2011; 4 Fletcher et al., 2014; Misra et al., 2017; Nelson et al., 2005; Siedlik et al., 2017; 5 Sugimura et al., 2016; Trepat et al., 2009), material-like properties (Bi et al., 2015; 6 Campàs et al., 2014; Latorre et al., 2018; Mongera et al., 2018; Pérez-González et al., 7 2019; Yang et al., 2017), and the analysis of the cellular packing (Curran et al., 2017; 8 Farhadifar et al., 2007; Gibson et al., 2006; Gibson et al., 2011; Gómez et al., 2021; 9 Honda, 1978; Lewis, 1928; Mao et al., 2013; Sanchez-Gutierrez et al., 2016; Thompson, 10 1945). As for the latter, the analysis of epithelial surfaces as tessellations of convex 11 polygons has revealed mathematical and physical principles with biological 12 consequences. One well-known example are the implications of the celebrated Euler´s 13 formula, ܸെܧ൅ܨൌ߯ (STAR Methods) (Euler, 1767). This formula implies that cells 14 in packed tissues have, on average, six neighbors (i.e., the average cellular connectivity 15 on a surface reads ۃ݊ଶ஽ۄൌ͸) (Reinhardt, 1918; Wetzel, 1926). This principle has 16 biological consequences, for example, the degree of cellular connectivity regulates the 17 strength of the cell-cell juxtracrine signaling (Guignard et al., 2020; Sharma et al., 2019; 18 Tung et al., 2012). 19 For a long time, the validity of this mathematical concept (i.e., each cell, on average, 20 connects with six neighboring cells) has been assumed in three dimensions (3D): 21 ۃ݊ଶ஽ۄൌ͸ฺۃ݊ଷ஽ۄൌ͸. Such an assumption is rooted in the common idealization of 22 epithelial cells as regular prismatic solids in either planar or bended epithelia. 23 However, the recent discovery of more complex cellular shapes in epithelia, i.e., 24 scutoids, that achieve an efficient 3D tissue packing has set a new paradigm that has 25 not been yet fully explored (Box 1A) (Gómez-Gálvez et al., 2018; Mughal et al., 2018; 26 Rupprecht et al., 2017). Scutoidal cellular shapes are the result of intercalations among 27 cells along the apico-basal axis (Box 1A-C and Fig. 1A). This phenomenon is then a 28 spatial version of the so-called T1 transitions that produce rearrangements of 29 neighboring cells in the plane as a function of time in numerous developmental 30 processes (Box 1B) (Bertet et al., 2004; Irvine and Wieschaus, 1994; Spencer et al., 31 4 2017). Scutoids imply necessarily changes in the neighboring relationship between 1 cells in a 3D spatial context and, consequently, modify the connectivity properties of 2 cells (Box 1C). Still, the analysis of tissue organization in 3D and the corresponding 3 biophysical insight have been hindered by the technical difficulties to accurately 4 segment and 3D-reconstruct cells, especially in curved tissues. In addition, very few 5 computational models account for the presence of apico-basal transitions to 6 investigate 3D self-organization in tissues (Gómez-Gálvez et al., 2018; Ioannou et al., 7 2020; Mughal et al., 2018; Okuda et al., 2019; Rupprecht et al., 2017). Moreover, from 8 an energetics viewpoint, while the appearance of scutoids can be explained by a 9 minimal model based on a surface/line tension minimization mechanism (Gómez- 10 Gálvez et al., 2018; Mughal et al., 2018; Okuda et al., 2019; Rupprecht et al., 2017), the 11 role played by additional energetic contributions to modulate the frequency of apico- 12 basal intercalations is unknown. 13 The analysis of 3D packing is in turn utterly relevant in cubic and columnar monolayer 14 tubular epithelia, where scutoids appear more frequently (Gómez-Gálvez et al., 2018; 15 Gómez et al., 2021; Iruela-Arispe and Beitel, 2013; Sanchez-Corrales et al., 2018). 16 Epithelial tubes are in fact the primary developmental structures in all organisms with 17 bilateral symmetry (Gilbert and Barresi, 2013), and tubulogenesis is fundamental in a 18 broad variety of key developmental processes, including gastrulation and neurulation 19 (Colas and Schoenwolf, 2001; Iruela-Arispe and Beitel, 2013; Leptin and Grunewald, 20 1990; Nelson, 2009; Pilot and Lecuit, 2005; Röper, 2018; Swanson and Beitel, 2006). 21 Furthermore, epithelial tubes are the essential functional unit of many mammalian 22 organs, including glands, components of the digestive apparatus, lungs, and kidney 23 (Huebner and Ewald, 2014). 24 Here, we study the packing and the 3D cellular connectivity properties of epithelial 25 tubes. We analyze the effect of different energetic contributions to modulate the 26 frequency of apico-basal intercalations; demonstrate that the presence of scutoids 27 implies a breakdown of the principle ۃ݊ଷ஽ۄൌ͸; and reveal a quantitative biophysical 28 principle that links the 3D cellular connectivity, energetics, and geometrical descriptors 29 (e.g., tissue curvature/thickness). Our findings are supported by i) a computational 30 model that realistically renders the 3D cellular organization of tubular epithelia 31 5 (including the appearance of scutoids); ii) experimental data of wildtype (wt) and 1 mutant epithelial tubes (DrosophilaƐƐĂůŝǀĂƌLJŐůĂŶĚͿǁŚŽƐĞϯĐĞůůƵůĂƌƐƚƌƵĐƚƵƌĞŚĂƐ 2 been accurately characterized by means of a computer-aided image analysis pipeline. 3 And iii), a biophysical model, supported by mathematical calculations, that connects 4 the tissue energetics with the 3D organization of epithelial tubes. 5 6 RESULTS 7 The Voronoi computational tubular model supports a relationship between energy 8 profiles and the intercalation propensity 9 To understand how the geometry of tubular epithelia and different energy 10 contributions affect the 3D cellular packing and connectivity of these tissues, we 11 designed and implemented a computational epithelial model that follows the 12 principles of Voronoi tessellations (Box 1D, E) (Gómez-Gálvez et al., 2018). In brief, we 13 generated 3D models of epithelial tubes by populating with seeds the apical surface, 14 ȭ௔ (light blue points in Fig. 1B) and implementing normal projections of those seeds 15 up to the basal surface, ȭ௕ (dark blue points in Fig. 1B). Each seed and its projection 16 corresponded to an individual cell of a tube. At each surface section ȭ (from apical to 17 basal) a 2D Voronoi diagram was performed, and the collection of those tessellations 18 rendered the 3D cellular geometry of cells (see details in STAR Methods). We point out 19 that we do not implement any temporal dynamics to the seeds. Thus, our 20 computational model is suited to static epithelial configurations as the ones 21 experimentally reported herein (see below). 22 Epithelial tubes appear in nature with very different thicknesses and cellular 23 arrangements. In order to explore how these features influence the 3D packing 24 properties of tubular epithelia we built diverse in silico Voronoi tubes. First, to 25 investigate the effect of tissue thickness we computed Voronoi tubes with different 26 surface ratios ݏൌܴȀܴ௔ (ܴ and ܴ௔ being the radial coordinate of the tube and the 27 apical radius respectively, Fig. 1C, D). We used ݏ-steps of 0.5 up to ݏൌͳͲ, so we were 28 able to explore 19 different values of the basal radius, ܴ௕ (Fig. 1D). Second, we 29 generated 10 different configurations in terms of the disorder level of the spatial 30 positions of the cellular seeds on the apical surface and the corresponding Voronoi 31 6 tessellations (V1 to V10, Fig. 1D). To that end, we used a fully random Voronoi 1 tessellation, i.e., randomized positions of cellular seeds, as the most disordered 2 pattern (V1). That configuration was made progressively more uniform (i.e., spatially 3 ŽƌĚĞƌĞĚͿĂĨƚĞƌŶŝŶĞƐƵĐĐĞƐƐŝǀĞŝƚĞƌĂƚŝŽŶƐŽĨƚŚĞŚŽŵŽŐĞŶŝnjŝŶŐ>ůŽLJĚƐĂůŐŽƌŝƚŚŵ;Box 4 1E and STAR Methods) (Fig. 1D). The resulting set of 10 different cellular arrangements 5 (V1 to V10) with increasing order properties conforms a Centroidal Voronoi 6 Tessellation (CVT) scale that has been proved useful to analyze the effect of the 7 topological organization of tissues and to simulate different tissues and/or 8 pathological conditions (Sanchez-Gutierrez et al., 2016; Vicente-Munuera et al., 2020). 9 We used the CVT scale to investigate how the average number of apico-basal 10 intercalations per cell, ۃ݅ۄ, changes as a function of the apico-basal coordinate, ݏ, and 11 the disorder level (Fig. 1D). As previously reported, we found that the number of apico- 12 basal transitions (Fig. 1D) and scutoids (Fig. S1) increased with ݏ (Gómez-Gálvez et al., 13 2018). As for the effect of the disorder level, we found that only in the case of fully 14 disordered tubes (i.e., V1: random case), and for low values of ݏ, there are more 15 intercalations, whereas for the rest of cases we observed that ۃ݅ۄ is fairly independent 16 of the CVT scale (Fig. 1D). 17 Energy contributions can be linked to geometric features of the shapes of epithelial 18 cells (Alt et al., 2017), see Box 1F-H. We used the set of Voronoi tubes (V1 to V10) to 19 explore surface tension, elasticity, and apical contractility energies, since these energy 20 contributions have been shown to play key roles in the organization of epithelia (Alt et 21 al., 2017; Farhadifar et al., 2007). As a first step, we estimated the average cellular 22 energy profiles as a function of ݏ in the computational tubular model (Fig. 1E-G). The 23 average surface tension energy (Box 1F) is related to the average lateral area of the 24 cells, ۃܣۄ, and therefore increase with the surface ratio, ݏ. Our results revealed that 25 ۃܣۄ is seemingly independent of the CVT scale (Fig. 1E). Consequently, the average cell 26 surface tension energy profile does not depend on the level of the topological disorder. 27 The average contractile energy (Box 1G) is related to the average and the variance of 28 the apical perimeter, ܮ, (Gilbert and Barresi, 2013; Farhadifar et al., 2007) therefore it 29 does not depend on the surface ratio, ݏ. The Voronoi model revealed that ۃܮۄ is CVT 30 independent, but the apical perimeter fluctuations decrease as the CVT scale increases 31 7 (Fig. 1F). Finally, the average cell elastic energy (Gelbart et al., 2012; Odell et al., 1981) 1 depends on the average and the variance of the cellular volume (Box 1H). Since the 2 average cell volume ۃܸۄ is, by construction, independent of the CVT scale (STAR 3 Methods), the average cellular elastic energy increases with the cellular volume 4 fluctuations, that in turn decrease with the CVT scale (Fig. 1G). 5 In order to evaluate how the appearance of scutoids is modulated by these energy 6 contributions for different values of the tissue thickness, we computed the cross- 7 correlation functions, ܥሺݏሻ, between the average cellular energy profiles, ۃܧ௓ۄ (ܼ being 8 ܣ or ܸ, i.e., surface tension or elastic terms), and the average number of apico-basal 9 intercalations, ۃ݅ۄ (Fig. 1H and STAR Methods). The cross-correlation measures the 10 similarity between two signals as a function of the displacement (or lag) of one signal 11 relative to the other. In our case the displacement/lag refers to the apico-basal 12 coordinate, ݏ, and consequently we inquire into the possibility that energetic 13 contributions either precede or follow the appearance of apico-basal intercalations. 14 Our results indicate that maximum correlations are obtained at zero lag independently 15 of the disorder level and that the appearance of scutoids correlates more significantly 16 with the surface tension energy profile than with the elastic energy: 95% vs. 80% 17 respectively. The latter is in agreement with previous studies that have shown that 18 surface tension energy minimization is the main cause underlying the appearance of 19 scutoids (Gómez-Gálvez et al., 2018; Gómez et al., 2021; Mughal et al., 2018). Also, 20 when assessing the extra effect of the energy input due to the apical contractility term 21 to the surface tension energy, ۃܧ஺ۄ൅ۃܧ௅ۄ, we found that it increases the correlation 22 between energy profile and the number of intercalations up to 98%, but it does not 23 lead to any change in the correlation due to elastic terms, ۃܧ௏ۄ൅ۃܧ௅ۄ (Fig. 1H). 24 We further examined the cross-correlation between the gradient of cellular 25 intercalations along the apico-basal axis, ߲௦ۃ݅ۄൌ߲ۃ݅ۄȀ߲ݏ, and the gradient of the 26 energy, ߲௦ۃܧ௓ۄ. In this way, we evaluated the level of correspondence between the 27 variation of the number of intercalations and the changes of the energy as a function 28 of the radial coordinate, ݏ. We found that, independently of the CVT scale, ߲௦ۃ݅ۄ 29 correlates slightly stronger with changes in the surface tension energy, ߲௦ۃܧ஺ۄ, than 30 with changes of the elastic contribution, ߲௦ۃܧ௏ۄ: ~80% versus ~75% respectively at 31 8 optimal lag (Fig. 1I). Notably, ߲௦ۃ݅ۄ lags behind ߲௦ۃܧ௓ۄ, i.e., the optimal lag for which 1 ܥሺݏሻ is the largest is located at ݏ൐ͲǤ Therefore, energy variations along the apico- 2 basal axis seem to precede changes in the number of intercalations that, in turn, 3 suggests an instructive role of the former over the latter. 4 Summing up, the Voronoi tubular model supports the idea that surface tension 5 energy is the more relevant contribution regulating the appearance of apico-basal 6 intercalations, and suggests that elastic terms play a role, yet less important than 7 surface tension, for modulating the intercalation propensity (see Discussion). 8 9 The Voronoi tubular model suggests a link between 3D tissue packing and energy 10 cues 11 To link quantitatively energy traits and 3D packing, we implemented a benchmark 12 able to simultaneously reveal the existence of apico-basal intercalations (scutoids) and 13 the polygonal distribution of cells at the apical and the basal surfaces. To that end, we 14 computed the probability that cells change their polygonal class between the apical 15 and basal surfaces. Thus, the components (i.e., bins) of this distribution along the 16 diagonal account for cells that have the same polygonal class at apical and basal 17 surfaces, whereas the spreading away from the diagonal ensures the existence of 18 scutoids and, consequently, changes in the cellular 3D connectivity (Fig. 2A and STAR 19 Methods). Our data revealed that, regardless of the value of the tissue thickness (and 20 the CVT scale), the dominant apical-basal polygonal class corresponds to cells with six 21 neighbors (Fig. 2A). As the tissue thickness, ݏ௕ൌܴ௕Ȁܴ௔, increases, more scutoidal 22 shapes with a distinct number of neighbors in apical and basal surfaces appeared. This 23 feature was revealed by the increasing value of the spreading away from the diagonal, 24 ߟଶ (STAR Methods and Table S1). In that regard, in agreement with the results shown 25 in Fig. 1D-G, our data indicates that ߟଶ increases with the tissue thickness and 26 decreases with the CVT scale (Fig. 2B). Also, the cross-correlation analysis between the 27 spreading coefficient and energy profiles agrees with Fig. 1H and reveals that 28 independently of the CVT scale the neighbor exchanges correlate more strongly with 29 the surface tension energy profile than with the elastic contribution: 90% vs. 70% at 30 zero (optimal) lag (Fig. 2C). Furthermore, we computed the average number of total 31 9 contacts between cells (i.e., the average 3D cellular connectivity), ۃ݊ଷ஽ۄ, as a function 1 of the surface ratio (i.e., the radial coordinate) and the Voronoi class (i.e., the level of 2 cellular disorder in the tissue) (Fig. S1). Our data indicated, that cells, on average, are 3 connected to more than six cells, i.e., ۃ݊ଷ஽ۄ൐͸, and the results are quantitatively 4 consistent with a mathematical derivation that shows that ۃ݊ଷ஽ۄ is linearly 5 proportional to the number of apico-basal intercalations ۃ݅ۄ: ۃ݊ଷ஽ۄൌ͸൅ۃ݅ۄȀʹ (STAR 6 Methods and Fig. S1). 7 8 The Voronoi tubular model recapitulates the properties of in vivo epithelial tubes 9 In order to compare the results obtained in our Voronoi computational tubular 10 models against the properties found in real tissues, we implemented a methodological 11 pipeline that combines several image analysis techniques to accurately reconstruct the 12 3D shapes of cells of in vivo epithelial tubes (Arganda-Carreras et al., 2017; Franco- 13 Barranco et al., 2021; Machado et al., 2019) (STAR Methods). We used the Drosophila 14 larval salivary gland, a cubic monolayer epithelium, as a model due to its ideal 15 characteristics to study complex tubular developmental structures (Girdler and Roper, 16 2014) (Fig. 3A). Also, cellular growth and division, as well as possible global tissue 17 deformation processes, do not occur in the Drosophila's salivary gland at the 18 developmental stage of our observations (i.e. the tissue is static); a fact that enables 19 the comparison with the Voronoi computational model. 20 We determined the average basal surface ratio (thickness) of the salivary glands, 21 ۃݏ௕ۄൌͺǤͷേͳǤͳ, the average percentage of scutoids, ͹ʹേͳʹΨ, the average 3D 22 connectivity, ۃ݊ଷ஽ሺݏ௕ሻۄൌ͸Ǥ͸േͲǤʹ, and the average number of apico-basal 23 intercalations per cell, ۃ݅ሺݏ௕ሻۄൌͳǤʹേͲǤ͵, thus confirming in vivo the validity of the 24 formula ۃ݊ଷ஽ۄൌ͸൅ۃ݅ۄȀʹ (STAR Methods and Fig. S1). We also calculated the 25 spreading coefficient of the 3D connectivity, ߟଶൌͳǤʹήͳͲିଶ, (Fig. 3B), and the 2D 26 polygonal distributions in the apical and basal surfaces (Fig. S2). We observed a small, 27 but significant, increase of the number of hexagons on the basal surface of the wt 28 glands (see Fig. S2, Table S1). 29 Further, in order to derive how energy contributions change as a function of the 30 apico-basal coordinate, ݏ, we implemented an algorithm that obtains the concentric 31 16 patterning, and cell fate determination (Guignard et al., 2020; Sharma et al., 2019; 1 Tung et al., 2012). Moreover, recent research has shown that cellular connectivity 2 regulates the viscosity of tissues (Petridou et al., 2021). 3 Therefore, our findings open new ways to draw implications about primary 4 developmental processes in which epithelial bending is essential such as 5 tubulogenesis, gastrulation, neurulation, and early embryo development. In addition, 6 we argue that, while our analyses focus on static tissues, our results could also be 7 relevant to understand active 3D tissue remodeling. Dynamic changes of ߚሺݏሻ would 8 modify the apico-basal intercalation propensity and therefore the material-like 9 properties: the larger ߚሺݏሻ the more solid-like the tissue would behave. 10 Recent studies have confirmed that adhesion-dependent active remodeling can be 11 connected to an increased activity of neighbor exchanges. In particular, loss of function 12 mutants of N-cadherin in the presomitic mesoderm of the zebrafish embryo cause an 13 increase in extracellular spaces and a solid-fluid jamming transition (Mongera et al., 14 2018). In addition, it has been recently shown that the stabilization of E-cad at the 15 cellular junctions in the Drosophila eye drives an increase of tension that can be 16 transmitted across the tissue (Founounou et al., 2021). This tension results in a 17 reinforcement of the solid-like tissue behavior. The salivary glands experiments 18 confirm that a reduction of E-cadherin increases the apico-basal intercalation 19 propensity. Our biophysical model predicted that such an increase of the apico-basal 20 intercalation propensity must be correlated with a decrease of the energy barrier ߚሺݏሻ. 21 Notably, this prediction was confirmed through the biophysical analyses of the ȴEcad 22 samples. 23 Finally, with respect to the applicability of our results to other areas, we expect that 24 the emerging field of organoids will benefit from our discoveries. A precise 25 quantification of the 3D connectivity could then help to understand the lack of 26 reproducibility in organoid production, one of the biggest challenges of the field 27 (Clevers, 2016; Huch et al., 2017; Schutgens et al., 2019). Also, from a medical point of 28 view, it has been recently shown that tissue curvature affects tumor progression due 29 to the imbalance of tensions in apical and basal surfaces of epithelial tubes (Messal et 30 al., 2019). Our study explains how cell energetics affect the 3D packing of these cells 31 17 and therefore may shed light on the mechanism of tumorigenic morphogenesis in 1 tubular organs. 2 3 AUTHOR CONTRIBUTIONS 4 P.G.-G. and P.V.-M implemented the Voronoi tubular model and designed the software 5 for salivary gland image processing and feature extraction. S.A. and J.B. developed the 6 biophysical model and performed the biophysical analyses. C. C.-A., A.T., A.V., C. G.-V., 7 P.G.-G. and P.V.-M segmented and curated the salivary glands confocal images. D.F-B 8 and I.A-C. implemented the 3D U-Net network architecture. J.A.A-S trained the 3D U- 9 net using the wildtype salivary glands, and inferred the cell outlines for the mutant 10 ones. A.T., C. G.-V, A.M.P. and V.A. conducted the biological experiments. C.G., R.R., 11 J.B. and A.M. theoretically described the 3D cell connectivity in Voronoi tubular models 12 depending on the surface ratio anisotropy. L.M.E., J.B. and P.G.-G. designed the figures. 13 L.M.E., J.B. and P.G.-G thought up the study, supervised the experiments, and wrote 14 the paper with input from all authors. All authors participated in the interpretation of 15 results, discussions, and the development of the project. 16 17 ACKNOWLEDGMENTS 18 L.M.E., P.V.-M. and P.G.-G. has been supported by the Ramón y Cajal program 19 (PI13/01347); L.M.E., P.V.-M. and P.G.-G. work is funded by the Ministry of Economy, 20 Industry and Competitiveness grant BFU2016-74975-P co-funded by FEDER funds, and 21 by the Spanish Ministry of Science and Innovation Ministry of Science through grant 22 PID2019-103900GB-I00 and from the LifeHUB Research Network (CSIC). A.T. is funded 23 ďLJĂCŽŶƚƌĂƚŽƉƌĞĚŽĐƚŽƌĂůW/&Ĩrom Universidad de Sevilla. L.M.E. and J.A.A.-S. work 24 is funded by the Junta de Andalucía (Consejería de economía, conocimiento, empresas 25 y Universidad) grant PY18-631 co-funded by FEDER funds. I.A-C. would like to 26 acknowledge the support of the 2020 Leonardo Grant for Researchers and Cultural 27 Creators, BBVA Foundation, and also the University of the Basque Country UPV/EHU 28 grant GIU19/027. S.A. and J.B. have been supported by a Faculty Innovation Grant (FIG) 29 by Lehigh University JB-FIG-2019. J.B. also acknowledges financial support from the 30 Spanish Ministry of Science and Innovation through grant PID2019-105566GB-I00 and 31 from the LifeHUB Research Network (CSIC). 32 18 1 DECLARATION OF INTERESTS 2 The authors declare no competing interests. 3 4 5 19 FIGURE TITLE AND LEGENDS 1 Box 1. Schematic representation of geometrical, topological and biophysical 2 concepts related to epithelial tissues. 3 A: Scutoids are prismatic-like geometric solids bounded between two surfaces (top 4 and bottom). Scutoids are characterized by three main properties: i) The shape of their 5 top and bottom bases, and of every parallel section between them, are polygons. ii) 6 The lateral surfaces of scutoids can be concave and/or convex surfaces such that a set 7 of scutoids can be packed together (laterally) without leaving any empty space. iii) 8 Scutoids have at least one vertex along the top-bottom axis such that when packed 9 together there are changes in the nearest-neighbors relationship. The example shows 10 a scheme of a stereotypical scutoid (left) and four scutoids packed together (right). 11 B: A T1-transition is a tissue rearrangement observed in epithelial surfaces where a 4- 12 cells' motif swaps nearest neighbors along time. T1-transitions enable tissue plasticity 13 through cellular reorganizations that lead, for example, to elongation in developing 14 tissues. 15 C: An apico-basal transition, aka an apico-basal intercalation, is a tissue 16 rearrangement along the apico-basal (top-bottom) axis of cells that lead to new 17 cellular contacts (nearest-neighbors exchange). An apico-basal transition is similar to 18 the T1-transition but instead of developing along the time it does along space. 19 D: A 2D tessellation (aka a 2D mosaic) is a partition of a surface with tiles that do not 20 overlap or leave any gaps. In this example, the tiles are octagons and squares. 21 E: A 2D Voronoi diagram is a particular type of tessellation built by convex polygons 22 (Voronoi cells). These polygons emerge from a set of generator seeds (black points), 23 such that each cell contains the region that is closer to its generating seed. The so- 24 called >ůŽLJĚƐ ĂůŐŽƌŝƚŚŵ makes the seeds of a Voronoi diagram to converge 25 progressively to the centroids (blue points): once a Voronoi diagram is obtained for a 26 set of seeds (black Voronoi diagram, ܸேͿĂŶŝƚĞƌĂƚŝŽŶŽĨƚŚĞ>ůŽLJĚƐĂůŐŽƌŝƚŚŵĐŽŶƐŝƐƚƐ 27 in repeating the Voronoi tessellation by replacing the seeds by the centroids of the 28 Voronoi tiles (blue Voronoi diagram, ܸேାଵͿdŚĞ>ůŽLJĚƐĂůŐŽƌŝƚŚŵŵĂŬĞƐƚŚĞsŽƌŽŶŽŝ 29 20 diagram progressively more ordered in terms of the polygonal distribution of the 1 Voronoi cells. 2 F: The surface tension energy is related to the cell-cell adherence through their 3 lateral area contacts. For each cell, the surface tension energy reads, ܧ஺ൌȦܣ, where 4 Ȧ and ܣ are the effective surface tension parameter and the cellular lateral area 5 respectively. Thus, the average surface tension energy of cells reads ۃܧ஺ۄൌȦۃܣۄൌ 6 ܧۃ஺ۄ and is independent of the fluctuations of ܣ. 7 G: The contractile energy is related to the polarized cortex activity of epithelia cells 8 at the apical surface. The contractile energy reads ܧ௅ൌȞܮଶ, where ܮ stands for the 9 cellular perimeter at the apical surface and Ȟ is the cortical tension energy per unit of 10 cell apical area. As a result, the average cell apical contractile energy increases with 11 the fluctuations of the apical perimeter: ۃܧ௅ۄൌȞۃܮଶۄൌܧۃ௅ۄ൅Ȟߪ௅ଶ where ߪ௅ଶൌ 12 ۃܮଶۄെۃܮۄଶ is the cellular apical perimeter variance. 13 H: The cell elastic energy is related to the volume conservation of cells. The cell elastic 14 energy reads, ܧ௏ൌ௒ଶሺܸെܸ଴ሻଶ, where ܻ is proportional to the Young modulus (a 15 quantification of the relationship between the cellular stress and strain) and, ܸ and ܸ଴ 16 represent the actual and target cellular volumes respectively. The average elastic 17 energy per cell increases with the fluctuations of the volume: ۃܧ௏ۄൌ௒ଶۃሺܸെܸ଴ሻଶۄൌ 18 ܧۃ௏ۄ൅௒ଶߪ௏ଶ where ߪ௏ଶൌۃܸଶۄെۃܸۄଶ is the cellular volume variance. 19 20 21 21 Figure 1. In silico Voronoi epithelial model: energetics and analysis of apico-basal 1 cell intercalations in tubes 2 A: (left) Scutoids entail apico-basal intercalations among packing cells that can be 3 envisioned as spatial T1 transitions to exchange neighbors (right): the green and the 4 red cells are neighbors in the basal surface (but not in the apical surface) while the 5 opposite is true for the blue and the yellow cells. B: In silico Voronoi 3D epithelial 6 models are generated by populating with cell seeds (circles) the apical surface, ߑ௔. The 7 location of cell seeds at any other surface/plane, ߑ, is obtained by implementing apical 8 normal projections, ࢔, up to the basal surface, ߑ௕. At each surface, ߑ, a 2D Voronoi 9 tessellation is performed and the 3D cellular shape of the cell is built upon the 10 collection of these tessellations. C: (top) In the particular case of tubular epithelia, 11 normal projections of apical cell seeds correspond to radial projections, and the 12 thickness/curvature of tubes are characterized by the surface ratio, ݏൌܴȀܴ௔ 13 (apical/basal surfaces: light/dark blue. ܴ: dashed yellow). (bottom) Illustrative 14 rendering of a Voronoi tube cell. D: The so-called CVT scale (iterations of Lloyd's 15 algorithm, STAR Methods) measures the topological disorder in in silico tubular 16 epithelia and leads to different cellular morphologies. In a V1 (Voronoi 1) model cell 17 seeds are randomly distributed on the apical surface to generate a planar Voronoi 18 tessellation (STAR Methods). By applying the Lloyd algorithm iteratively, the apical 19 topological disorder diminishes (top to bottom: V1, V5, and V10 examples). The 20 random location of seeds in V1 implies the emergence of a wide range of different 21 polygonal cell types. As Lloyd's algorithm iterates, the larger the tessellation order is 22 (Box 1E). We observe this progressive ordering from a V1 to V10 as polygonal 23 distributions converge to results with a larger proportion of hexagons and a reduction 24 of the other polygonal shapes (inset: polygon distribution insets for V1, V5, and V10). 25 The average number of apico-basal intercalations per cell in in silico tubes (݊ൌʹͲ 26 realizations per CVT scale, each tube composed by 200 cells) increases with the tissue 27 thickness (surface ratio) but does not change significantly with the CVT scale. The 28 black/red/green arrows correspond to the illustrative examples of the planar apical 29 tessellations shown on the left. E: The average lateral cell area as a function of the CVT 30 and the surface ratio (values normalized to the V1 case at ݏ௕ൌͳͲ, ݊ൌʹͲ per CVT 31 scale) indicates that the average surface tension energy does not depend on the level 32 22 of topological disorder. F: (Left) Green solid circles represent the normalized (with 1 respect to the V1 case) average apical cell perimeter for each Voronoi tube sample as 2 a function of the CVT scale. (Right) Purple solid circles account for the standard 3 deviation of apical cell perimeter for each Voronoi tube sample. Long vertical black 4 lines denote the mean and small segments the standard deviation (݊ൌʹͲ per CVT 5 scale). G: The cellular volume variance is a proxy for the average elastic cell energy 6 (STAR Methods). The latter decreases with the CVT scale and increases with the 7 surface ratio (tissue thickness) (݊ൌʹͲ per CVT scale). Volume values were normalized 8 with respect to the V1 case with ۃܸሺݏ௕ൌͳͲሻۄൌͳ. H: Cross-correlation between 9 average energy profiles along the apico-basal axis, ۃܧ஺ۄ (dark grey), ۃܧ௏ۄ (red), and 10 ۃܧ஺ۄ൅ۃܧ௅ۄ (light grey), and the number of apico-basal intercalations, ۃ݅ۄ. Solid lines 11 stand for the averages among disorder configurations (i.e., CVT scale) and the dotted 12 lines delimit the standard deviation band. ۃܧ௏ۄ൅ۃܧ௅ۄ has not been plotted since the 13 extra contribution of the contractile term, ۃܧ௅ۄ, does not modify the correlation 14 function. I: Cross-correlation between energy gradients (߲௦ۃܧ஺ۄ and ߲௦ۃܧ௏ۄ) and the 15 gradient of intercalations along the apico-basal axis (߲௦ۃ݅ۄ). Color code as in H. 16 23 Figure 2. Three-dimensional packing and connectivity properties of the Voronoi 1 tubular model 2 A: (left) Schematic representation of a 3D histogram that accounts for the probability 3 that cells have ݊௔ (number) of neighbors in the apical surface and ݊௕ neighbors in the 4 basal surface. Cells with the same polygonal class in apical and basal surfaces 5 contribute to bins along the diagonal (red squares). The bins spreading away from the 6 diagonal (green squares) ensures the presence of scutoidal cells. E.g.: the red and 7 green cells shown in the plot contribute to the bins indicated in the 3D histogram (red 8 and green stars respectively). (right) 3D histograms of V5 tubes for increasing values 9 of the surface ratio. A larger value of the spreading coefficient, ߟଶ, (STAR Methods) 10 indicates an increasing number of scutoids. B: Density plot showing the value of the 11 spreading coefficient, ߟଶ, of 3D histograms as a function of the surface ratio and the 12 Voronoi class in in silico tubes (݊ൌʹͲ tubes per CVT scale). C: Cross-correlation 13 between average energy profiles along the apico-basal axis, ۃܧ஺ۄ (dark grey), ۃܧ௏ۄ 14 (red), and ۃܧ஺ۄ൅ۃܧ௅ۄ (light grey), and the spreading coefficient, ߟଶ. Solid lines stand 15 for the averages among disorder configurations (i.e., CVT scale) and the dotted lines 16 delimit the standard deviation band. ۃܧ௏ۄ൅ۃܧ௅ۄ has not been plotted since the extra 17 contribution of the contractile term, ۃܧ௅ۄ, does not modify the correlation function. 18 24 Figure 3. DrosophilaƐ salivary gland analysis 1 A: (top) Full projection of a representative wildtype salivary gland (cell contours 2 stained by Cy3-labeled phalloidin, STAR Methods). (middle) Computer rendering of the 3 segmented salivary gland shown on top. Scale bar ͳͲͲߤ݉. (bottom) 3D rendering of a 4 representative segmented salivary gland. B: Density plot of the 3D distribution of 5 neighbor exchanges between apical and basal surfaces as a function of the number of 6 neighbors in apical, ݊௔, and basal, ݊௕, surfaces (as in Fig. 2A) in wildtype salivary glands 7 (݊ൌʹͲ glands, 979 cells). C: Average profiles of the number of apico-basal 8 intercalations, ۃ݅ۄ, average lateral area, ۃܣۄ, and cellular volume fluctuations, ߪ௏ଶ in in 9 vivo tubes as a function of the apico-basal coordinate, ݏ (݊ൌʹͲ glands, similar to Fig. 10 1D-E, G). D-E: Cross-correlation analysis between energy and intercalation profiles 11 (same as Fig. 1H-I). The error band indicates in this case the variability among 12 experimental samples (݊ൌʹͲ glands). 13 25 Figure 4. Biophysical model: in silico and in vivo results 1 A: (top) Cell intercalations along the apico-basal axis can be visualized as spatial T1 2 transitions (non-reversible: once a neighbor is won it cannot be lost). (bottom) The 3 ƉŽŽƌŐĞƚƌŝĐŚĞƌƉƌŝŶĐŝƉůĞ;STAR Methods) suggests an increasing energetic cost (i.e., 4 a larger activation energy) for recruiting new neighbors as a function of the number of 5 neighbors already won. In our model, ߚሺݏሻ accounts for the energetic cost per 3D 6 neighbor (per apico-basal intercalation) to recruit a new neighbor (STAR Methods). 7 B: The energy landscape shown in B (bottom) can be modeled by a stochastic dynamics 8 (a Kolmogorov rate equation) where cells increase their 3D neighbors with a 9 probability per unit of surface ratio, ݎ௡ǡ௠, that depends on the activation energy, ߚሺݏሻ, 10 ĂďĂƌĞƚƌĂŶƐŝƚŝŽŶƌĂƚĞߙ, and the maximum cell connectivity ܰ௠௔௫ (STAR Methods). 11 C: Results of the Kolmogorov model in V8 (ݏ௕ൌͳǤ͹ͷ) in silico tubes (݊ൌʹͲ). The 12 left/center density plots represent the cellular connectivity distribution (i.e., the 13 fraction/probability of cells with a given number of 3D neighbors) as a function of ݏ 14 obtained in the Voronoi simulation (left) and as predicted by the Kolmogorov model 15 (center); the red circles (left/right) indicate the average number of 3D neighbors per 16 cell ۃ݊ଷ஽ۄǢ the red line (center/right) shows ۃ݊ଷ஽ۄ as obtained by the Kolmogorov 17 model. The density plot on the right shows the difference between the predicted and 18 the actual connectivity distributions and the corresponding error, ߝଶ (magenta lines). 19 D: Same as panel C but results obtained in salivary glands (݊ൌʹͲ glands). The 20 maximum value of ݏ in the analyzed radial sections of the glands is ݏൌ͸Ǥͷ. This value 21 being the largest radial section of the smallest gland (STAR Methods). 22 32 ݊ଷ஽ሺݏሻൌ݉ܽݔሺሼܸሺݏሻሽሻൌ݉݅݊ሺሼܸሺݏሻሽሻ൅݅ሺݏሻ (3) 1 where ሼܸሺݏሻሽൌሼܸሺͳሻǡܸሺͳ൅݀ݏሻǡڮǡܸሺݏ௕ሻሽ and ݅ሺݏሻ denotes the number of 2 intercalation points in the interval ݏאሾͳǡݏ௕ሿ. In the case of a 3D tessellation with ܰ 3 cells, where ܯ of them do not show any intercalation, the total number of 4 accumulated neighbors reads, 5 ݊ଷ஽ሺݏሻൌσ݊ଷ஽ ሺ௝ሻሺݏሻ ே ௝ୀଵ ൌσܸሺ௝ሻሺͳሻ ெ ௝ୀଵ ൅σ݉ܽݔ൫൛ܸሺ௝ሻሺݏሻൟ൯ ேିெ ௝ୀଵ ൌ 6 σܸሺ௝ሻሺͳሻ ெ ௝ୀଵ ൅σ൛݉݅݊൫൛ܸሺ௝ሻሺݏሻൟ൯൅݅ሺ௝ሻሺݏሻൟ ேିெ ௝ୀଵ (4) 7 Given that each intercalation point is shared by four cells, two of them necessarily 8 increase their number of vertices in a given ݏ-plane and two of them decrease their 9 number of vertices. Thus, in the case of a decrease ݉ܽݔ൫൛ܸሺ௝ሻሺݏሻൟ൯ൌܸሺ௝ሻሺͳሻ and in 10 the case of an increase ݉݅݊൫൛ܸሺ௝ሻሺݏሻൟ൯൅݅ሺ௝ሻሺݏሻൌܸሺ௝ሻሺͳሻ൅݅ሺ௝ሻሺݏሻ. Consequently, 11 ݊ଷ஽ሺݏሻൌσܸሺ௝ሻሺͳሻ ே ௝ୀଵ ൅σ݅ሺ௝ሻሺݏሻ ሺேିெሻȀଶ ௝ୀଵ ൌσܸሺ௝ሻሺͳሻ ே ௝ୀଵ ൅ଵଶσ݅ሺ௝ሻሺݏሻ ேିெ ௝ୀଵ (5) 12 where we used the fact that for every intercalation event that increases by one the 13 number of neighbors there is one that decreases the number of neighbors in the same 14 amount; consequently, we can add up all intercalation events and divide by two. Hence 15 the average number of accumulated 3D neighbors, ۃ݊ଷ஽ሺݏሻۄൌ݊ଷ஽ሺݏሻȀܰ reads 16 ۃ݊ଷ஽ሺݏሻۄൌۃܸሺͳሻۄ൅ۃ݅ሺݏሻۄȀʹ; ۃ݅ሺݏሻۄ being the average number of apico-basal 17 intercalations per cell. Finally, by considering that any ݏ-surface, and in particular the 18 apical surface ݏൌͳ, corresponds to a 2D tessellation of convex polygons, ۃܸሺͳሻۄൌ͸ 19 we conclude that, 20 ۃ݊ଷ஽ሺݏሻۄൌ͸൅ଵଶۃ݅ሺݏሻۄ (6) 21 22 dŚĞ ϯ ŶĞŝŐŚďŽƌƐ ĂĐĐƵŵƵůĂƚŝŽŶ ŝŶ ƚƵďƵůĂƌ ĞƉŝƚŚĞůŝĂ ĨŽůůŽǁƐ Ă ƉŽŽƌ ŐĞƚ ƌŝĐŚĞƌ 23 principle 24 In order to investigate additional phenomena that could help to understand how 25 the 3D cellular connectivity is regulated, we computed the net gain of cellular 26 neighbors in epithelial tubes as a function of the 2D polygonal cell class at the apical 27 surface. We observed, both in the salivary glands and in the Voronoi model (in 28 particular in the case V8 (ݏ௕ൌͳǤ͹ͷ) that compares the best with in vivo glands), that 29 33 the smaller the number of neighbors of a cell at the apical surface, the larger the net 1 gain of 3D cellular contacts (Fig. S2). This behavior was also obtained with respect to 2 the 2D polygonal cell class at the basal surface (Fig. S2). These results suggest that, in 3 ƚĞƌŵƐŽĨƚŚĞĐĞůůƵůĂƌƉĂĐŬŝŶŐƚƵďƵůĂƌĞƉŝƚŚĞůŝĂĨŽůůŽǁĂƉŽŽƌŐĞƚƌŝĐŚĞƌƉƌŝŶĐŝƉůĞƚŚĞ 4 less neighbors a cell has in a surface (apical or basal), the larger the net increase of 3D 5 cellular contacts. 6 7 In Voronoi tubes the net gain of 3D neighbors is bounded 8 The ƉŽŽƌŐĞƚƌŝĐŚĞƌ behavior can be justified by mathematical arguments that 9 show that the probability to increase the cellular connectivity necessarily decreases 10 with the number of current neighbors (Fig. S3). Assuming a cylindrical geometry (e.g., 11 epithelial tubes), each point at a given radial surface can be represented into the 12 Cartesian plane; where coordinate ݔ accounts for the cylindrical transversal coordinate 13 and coordinate ݕ for the longitudinal one (see Fig. S3). Thus, if the coordinates of a 14 point (e.g., a Voronoi seed) at the apical surface are given by ሺݔǡݕሻ, the coordinates of 15 that point at a surface with a value of the cylindrical radial expansion ݏאሾͳǡλሻ can 16 be found by defining the function ݂௦ǣԹଶืԹଶ ݂௦ሺݔǡݕሻൌሺݏݔǡݕሻ. Under these 17 conditions, we aim to characterize the seeds that generate scutoids (exchanges in the 18 neighboring relations of seeds) as ݏ changes. 19 Lemma 1. Given three nonʹcolinear points ሼܣǡܤǡܥሽ that define a circle (a nearest- 20 neighbors relation), and another exterior point ܦ, if ݏ൐ͳexists such that ݂௦ሺܦሻ is 21 interior to the circle defined by ሼ݂௦ሺܣሻǡ݂௦ሺܤሻǡ݂௦ሺܥሻሽ, then ܦ is inside of the vertical 22 parabola containing ሼܣǡܤǡܥሽ (Fig. S3). 23 Remark. If two of the three points ሼܣǡܤǡܥሽ are on the same vertical line, then the 24 parabola considered in Lemma 1 degenerates as a vertical strip. Even in this case, the 25 thesis of the Lemma is true if we replace the interior of the parabola by the inside of 26 the strip. 27 Proof. Without loss of generality, we can suppose that ሼܣǡܤǡܥሽ are counterclockwise 28 oriented and that they have Cartesian coordinatesሺܽଵǡܽଶሻ, ሺܾଵǡܾଶሻ and ሺܿଵǡܿଶሻ 29 respectively. 30 34 Thus, the point ܦሺݔǡݕሻ is outside the circle defined by ሼܣǡܤǡܥሽ if, and only if, the 1 sign of the following determinant is negative: 2 ተ ተ ܽଵܽଶܽଵଶ൅ܽଶଶͳ ܾଵܾଶܾଵଶ൅ܾଶଶͳ ܿଵܿଶܿଵଶ൅ܿଶଶͳ ݔ ݕ ݔଶ൅ݕଶͳ ተ ተ ൌ ተ ተ ܽଵܽଶܽଵଶͳ ܾଵܾଶܾଵଶͳ ܿଵܿଶܿଵଶͳ ݔ ݕ ݔଶͳ ተ ተ ൅ ተ ተ ܽଵܽଶܽଶଶͳ ܾଵܾଶܾଶଶͳ ܿଵܿଶܿଶଶͳ ݔ ݕ ݕଶͳ ተ ተ ൏Ͳ (7) 3 For the sake of simplicity, we represent the previous equation as: 4 ݀݁ݐሺࣛሻൌ݀݁ݐሺࣜሻ൅݀݁ݐሺࣝሻ൏Ͳ (8) 5 On the other hand, by considering ݔ and ݕ as variables, the equation ݀݁ݐሺࣛሻൌͲ 6 corresponds to the circle defined by ሼܣǡܤǡܥሽ, and ݀݁ݐሺࣜሻൌͲ corresponds to the 7 vertical parabola defined by the same three points. Consequently, the inequality 8 ݀݁ݐሺࣜሻ൐Ͳ defines the locus of interior points to that parabola. 9 Now, assuming that ݏ൐ͳ exists such that ݂௦ሺܦሻ is interior to the circle defined by 10 ሼ݂௦ሺܣሻǡ݂௦ሺܤሻǡ݂௦ሺܥሻሽ. Then, 11 ተ ተ ݏܽଵܽଶݏଶܽଵଶ൅ܽଶଶͳ ݏܾଵܾଶݏଶܾଵଶ൅ܾଶଶͳ ݏܿଵܿଶݏଶܿଵଶ൅ܿଶଶͳ ݏݔ ݕ ݏଶݔଶ൅ݕଶͳ ተ ተ ൌݏଷ݀݁ݐሺࣜሻ൅ݏ݀݁ݐሺࣝሻ൐Ͳ (9) Or, equivalently, ݏଶ݀݁ݐሺࣜሻ൅݀݁ݐሺࣝሻ൐Ͳ, so, ݏଶ݀݁ݐሺࣜሻ൐െ݀݁ݐሺࣝሻ. If ݀݁ݐሺࣜሻ൏ 12 Ͳ, then ͳ൏ݏଶ൏െௗ௘௧ሺࣝሻ ௗ௘௧ሺࣜሻ and therefore ݀݁ݐሺࣜሻ൐െ݀݁ݐሺࣝሻ. The latter is in 13 contradiction with ݀݁ݐሺࣜሻ൅݀݁ݐሺࣝሻ൏Ͳ. As a result, ݀݁ݐሺࣜሻ൐Ͳ, and the following 14 inequality holds, 15 ݏଶ൐െ݀݁ݐሺࣝሻ ݀݁ݐሺࣜሻ൐ͳ (10) Notice that if the circle defined by ሼܣǡܤǡܥሽ is surrounded by a set of points and we 16 change continuously the parameter ݏ in the interval ሾͳǡλሻ, it is possible to detect the 17 first point touching the circle defined by ሼ݂௦ሺܣሻǡ݂௦ሺܤሻǡ݂௦ሺܥሻሽ. That point can be 18 obtained by computing all the points at ݏൌටെௗ௘௧ሺࣝሻ ௗ௘௧ሺࣜሻ. Hence, the first point 19 contacting the circle will be that with the minimum value of ݏ. 20 35 As for proving that the average of the number of neighbors of a cell induced by a seed 1 grows is bounded as a function of the surface ratio, we state the following proposition: 2 Proposition 1. Given a Voronoi seed representing a cell, if ݊ଷ஽ሺݏሻ is the total number 3 of accumulated cell neighbors as ݏ increases from ݏൌͳ (apical surface) to a given 4 value of ݏ, then ۃ݊ଷ஽ሺݏሻۄ is a bounded function for a finite cylinder. 5 6 Proof. We model the apical surface as the cylinder ʹߨݎൈ݄, where ݎ representes the 7 inner radius and ݄ the length of the cylinder. Given a seed ܣ in that surface, in the 8 corresponding Delaunay triangulation it appears as a point surrounded by triangles 9 defining the neighborhood of ܣ. By Lemma 1, each triangle defines a vertical parabola 10 and a circle. So, any other seed touching ܣ in other layer must be inside of one of the 11 parabolas and outside of all circles (see Fig. S3Ϳ>ĞƚƐĚĞŶŽƚĞ࣬௦ǡ஺ the feasible region 12 for a new neighbor of ܣ in the layer represented by ݏ, i.e., all points inside one of the 13 parabolas and outside all the circles. Thus, if ሺ࣬௦ǡ஺ሻ is the number of seeds in that 14 region that are not neighbors of ܣ in the apical surface, obviously, an upper bound to 15 the number of new neighbors to ܣ is given by ሺ࣬௦ǡ஺ሻ൑ሺ࣬ଵǡ஺ሻ. 16 17 On the other hand, that number of seeds is, in average, proportional to the density 18 of seeds times the area of ࣬௦ǡ஺, therefore, the average number of accumulated 19 neighbors of ܣ, denoted as ۃ݊ଷ஽ሺܣሻۄ, will be bounded by the change of the density of 20 points when growing ݏ, this is to say, 21 ݀ۃ݊ଷ஽ሺܣሻۄ൑ܯڄ ࣬ೞǡಲ ଶగ௦௥ڄ௛݀ݏ (11) 22 where ܯ represents the total number of seeds (i.e., the total number of cells that is 23 a constant) and the quotient is the area of ࣬௦ǡ஺ divided by the area of a given radial 24 layer. In general, it is not possible to integrate equation (4), since the area of ࣬௦ǡ஺ is 25 known only in very few, particular, cases. 26 In the case of a finite cylinder, ۃ݊ଷ஽ሺܣሻۄ൑ሺ࣬௦ǡ஺ሻ൑ሺ࣬ଵǡ஺ሻ leads, summing up to 27 all the seeds and dividing by ܯ, to the upper bound 28 ۃ݊ଷ஽ሺݏሻۄ൑ଵ ெڄσ஺ሺ࣬ଵǡ஺ሻ (12) 29 thus, ۃ݊ଷ஽ሺݏሻۄ is necessarily a bounded function. This expression indicates that the 30 number of new neighbors when increasing ݏ exhausts since the number of cells is a 31 36 resource shared by all the layers. It is possible to obtain an upper bound to ܰ௠௔௫ൌ 1  ௦՜ஶۃ݊ଷ஽ሺݏሻۄ since, after a flip in the Delaunay triangulation, the edge disappearing 2 (i.e., a cell contact loss) can never be recovered in a cylindrical geometry. Thus, ܯή 3 ൫ܰ௠௔௫െ݊ଷ஽ሺͳሻ൯ is bounded by the number of edges that complement the original 4 Delaunay triangulation on the apical surface, that is, 5 ܰ௠௔௫െۃ݊ଷ஽ሺͳሻۄ൑ଵ ெڄቀெሺெିଵሻ ଶെܯۃ௡యವሺଵሻۄ ଶቁൌெିଵ ଶെۃ௡యವሺଵሻۄ ଶ (13) 6 leading to 7 ܰ௠௔௫൑ெିଵ ଶ൅ۃ௡యವሺଵሻۄ ଶ൑ெିଵ ଶ൅͵ൌெାହ ଶ (14) 8 Where we have assumed that ۃ݊ଷ஽ሺͳሻۄൌ͸. The simulations of the computational 9 Voronoi model and the data of the salivary gland show that ܰ௠௔௫ is in fact much 10 smaller that the theoretical bound ெାହ ଶ. 11 12 A Kolmogorov rate equation for the 3D cellular connectivity 13 The equation for how the probability, ܲ௠, of having ݉ accumulated 3D neighbors 14 (i.e., ݉ൌ݊ଷ஽) changes as the surface ratio (apico-basal dimensionless radial 15 coordinate) increases from ݏ to ݏ൅݀ݏ can be described by the following Markov 16 equation (Fig. 4A-B), 17 ܲ௠ሺݏ൅݀ݏሻൌܲ௠ሺݏሻܶ௠ǡ௠൅ܲ௠ିଵሺݏሻܶ௠ିଵǡ௠ (15) 18 where ܶ௡ǡ௠ is the probability of changing the number of neighbors from ݊ to ݉ due to 19 an apico-basal intercalation. Since σܶ௡ǡ௠௠ ൌͳ (normalization of the transition 20 probabilities) and ܶ௡ǡ௠ൌ݂ሺ݊ǡ݉ሻ൛ߜ௡ିଵǡ௠൅ߜ௡ǡ௠ାଵൟ (each intercalation can only 21 possibly induce to win one neighbor) then ܶ௠ǡ௠ൌͳെܶ௠ǡ௠ାଵ and the above Markov 22 equation can be written as a Kolmogorov equation (a.k.a. Master equation): 23 ௗ௉೘ሺ௦ሻ ௗ௦ ൌܲ௠ିଵሺݏሻݎ௠ିଵǡ௠െܲ௠ሺݏሻݎ௠ǡ௠ାଵ (16) 24 where ݎ௡ǡ௠ accounts for the probability of apico-basal intercalations per unit of 25 surface ratio, i.e., ܶ௡ǡ௠ൌݎ௡ǡ௠݀ݏ. We point out that our model accounts for apico-basal 26 intercalations that occur due to curvature effects such that ݏൌܴȀܴ௔ changes along 27 the apico-basal coordinate of cells (i.e., bending, folding). Thus, our model does not 28 37 capture the apico-basal intercalations that develop due to active cellular processes 1 (e.g. cellular extrusion, cell divisions) either in bended tissues or in the case of tissue 2 planar geometries. 3 By following the Eyring model (Eyring, 1935), i.e., if we assume an Arrhenius-like 4 kinetics such that to win neighbors there is an energy cost (see (Bi et al., 2014)) then 5 ݎ௠ǡ௠ାଵൌߙො݁ିοா೘, where ߙො is the so-called pre-exponential factor that modulates the 6 ďĂƌĞĨƌĞƋƵĞŶĐLJŽĨŝŶƚĞƌĐĂůĂƚŝŽŶƐ;ƉĞƌƵŶŝƚŽĨƐƵƌĨĂĐĞƌĂƚŝŽĞdžƉĂŶƐŝŽŶݏ) and οܧ௠ is a 7 dimensionless activation energy (in units of the effective thermal energy associated 8 with membrane fluctuations ߦ (Marmottant et al., 2009) dŚĞ ŽďƐĞƌǀĞĚ ƉŽŽƌ ŐĞƚ 9 ƌŝĐŚĞƌďĞŚĂǀŝŽƌƐƵŐŐĞƐƚƐƚŚĂƚƚŚĞĂĐƚŝǀĂƚŝŽŶĞŶĞƌŐLJοܧ௠, increases with ݉. This can 10 be explained as a result of a cumulative process if we assume that each neighbor that 11 is gained implies to overcome an energy barrier, ߚሺݏሻ, through an apico-basal 12 intercalation. Consequently, ݁ିοா೘ൌς݁ିఉሺ௦ሻ ௡ୀ௠ ௡ୀଵ ൌ݁ି௠ήఉሺ௦ሻ. Thus, ߚሺݏሻ represents 13 the dimensionless activation energy of a cell per 3D neighbor or, in the context of the 14 different energetic contributions reviewed in this manuscript, to the energy barrier 15 required to perform a spatial T1-transition following a surface energy minimization 16 process (Gómez-Gálvez et al., 2018; Mughal et al., 2018). As for the dependence of ߚ 17 on ݏ, the simplest mathematical form that recapitulates the fact that the apical and 18 basal surfaces accumulate more cell-cell adherent complexes (either in wt or mutant 19 phenotypes) is quadratic (Fig. S5): ߚሺݏሻൌߚ଴ቀͳ൅ఋଶሺݏെݏ଴ሻଶቁ. The average, along 20 the apico-basal coordinate, of the energy cost then reads ߚҧൌͳȀሺݏ௕െͳሻ׬ߚሺݏሻ ௦್ ଵ݀ݏ. 21 On the other hand, the mathematical calculations indicate that the intercalation rate 22 ݎ௠ǡ௠ାଵ becomes null for a finite value of ݉ or, alternatively, that the activation energy 23 becomes infinite for a finite value of ݉. Otherwise, the net gain of new neighbors is 24 not bounded. This fact can be accounted for by assuming that the bare frequency is a 25 function of the number of neighbors, ߙොൌߙොሺ݉ሻ, such that ௗఈෝ ௗ௠൏Ͳ and becomes null 26 for a finite value of ݉. For the sake of simplicity, we assume that up to first order in ݉: 27 ߙොൌߙሺܰ௠௔௫െ݉ሻ, where ܰ௠௔௫ is the asymptotic, maximum, number of 3D neighbors 28 a cell can possibly have. Summarizing, the apico-basal intercalation rate ݎ௠ǡ௠ାଵ reads, 29 ݎ௠ǡ௠ାଵൌߙሺܰ௠௔௫െ݉ሻ݁ି௠ఉሺ௦ሻ (17) 30 38 Under these conditions, the Kolmogorov equation reads, 1 ௗ௉೘ሺ௦ሻ ௗ௦ ൌߙ൫ܰ௠௔௫െሺ݉െͳሻ൯݁ିఉሺ௦ሻሺ௠ିଵሻܲ௠ିଵሺݏሻെߙሺܰ௠௔௫െ݉ሻ݁ିఉሺ௦ሻ௠ܲ௠ሺݏሻ 2 (18) 3 On the other hand, the equation satisfied by the average number of accumulated 3D 4 neighbors, ۃ݊ଷ஽ۄൌۃ݉ۄ, reads, 5 ௗۃ௠ሺ௦ሻۄ ௗ௦ ൌσ݉ௗ௉೘ሺ௦ሻ ௗ௦ ௠ൌσݎ௠ǡ௠ାଵܲ௠ሺݏሻ ௠ൌۃݎ௠ǡ௠ାଵۄ (19) 6 We notice that this equation implies an important role of the disorder (i.e. the 7 distribution ܲ௠): even in conditions under which the transition rate, ݎ௠ǡ௠ାଵŝƐůĂƌŐĞ 8 the resulting growth of 3D neighbors, ௗۃ௠ሺ௦ሻۄ ௗ௦ , and, consequently, the net accumulation 9 of 3D neighbors, can be more prominent in conditions where the transition rate is 10 ƐŵĂůůdŽŝůůƵƐƚƌĂƚĞƚŚŝƐĞĨĨĞĐƚǁĞĐŽŶƐŝĚĞƌƚŚĞĨŽůůŽǁŝŶŐĞdžĂŵƉůĞ&ŽƌƚŚĞƐĂŬĞŽĨ 11 simplicity, we evaluate the initial growth of 3D neighbors starting from the apical 12 surface, i.e. we particularize Eq. (19) to the case ݏൌͳ (and hence according to Euler's 13 formula ۃܲ௠ሺͳሻۄൌσ݉ܲ௠ሺͳሻൌ͸ ௠) and consider two possible conditions: a fully 14 ordered (݋) distribution with ܲ௠௢ൌߜ௠ǡ଺ (i.e., all hexagons) and a disordered (݀) 15 condition that combines with equal probability cells with 3, 6, and 9 sides in the apical 16 surface, i.e. ܲ௠ௗൌଵଷ൫ߜ௠ǡଷ൅ߜ௠ǡ଺൅ߜ௠ǡଽ൯. We also assume that ߚሺݏሻ؄ߚҧ (i.e., we 17 approximate the energy cost to gain new 3D neighbors by its average) and that ߚҧ௢൏ 18 ߚҧௗ, and also that ߙ௢൏ߙௗ (Fig. S6). Under these conditions, for the same ܰ௠௔௫, the 19 following holds, ݎ௠ǡ௠ାଵ ௢൐ݎ௠ǡ௠ାଵ ௗ, that is, the transition rate to gain new 3D neighbors 20 is larger in the ordered case than in the disorder case. This is in fact the situation that 21 we observed in the Voronoi tubular model when we estimated the value of the energy 22 barrier to gain new neighbors: ߙ and ߚሺݏሻ decrease as the CVT scale increases even 23 though the surface tension energy is independent of the CVT scale (see Fig. S6, Fig. 1, 24 and Table S1). However, it is possible to find large regions in terms of the values of 25 ߚҧ௢, ߚҧௗ, ߙ௢, and ߙௗ where ௗۃ௠ሺଵሻۄ ௗ௦ ቚௗ൐ௗۃ௠ሺଵሻۄ ௗ௦ ቚ௢ (Fig. S6). That is, the growth of 3D 26 neighbors starting from the apical surface (i.e. ݏൌͳ) in the disorder case can be 27 actually larger than that of the order case despite the fact that the transition rate to 28 gain new 3D neighbors is smaller in the former. 29 39 Also, from Eq. (19), it is possible to infer, approximately, the expected behavior of 1 ۃ݉ሺݏሻۄൌۃ݊ଷ஽ሺݏሻۄ as follows. First, by performing a mean-field-like approximation, 2 i.e., ۃܨሺ݉ሻۄൎܨሺۃ݉ۄሻ, 3 ௗۃ௠ۄ ௗ௦ ൎߙሺܰ௠௔௫െۃ݉ۄሻ݁ିఉሺ௦ሻۃ௠ۄ (20) 4 Second, assuming that ߚሺݏሻ൏ͳ (otherwise it is difficult to justify the observed 5 presence of apico-basal intercalations), 6 ௗۃ௠ۄ ௗ௦ ൎߙሺܰ௠௔௫െۃ݉ۄሻሺͳെߚሺݏሻۃ݉ۄሻ൅ࣩሺߚଶሻ (21) 7 Eq. (21) is formally a logistic-like growth equation that can solved subjected to the 8 condition ۃ݉ሺͳሻۄൌ͸ (the average number of neighbors in the apical surface is ͸). 9 We notice that in this case, the disorder levels of the wt and the mutant glands are 10 similar. Consequently, the accumulation of 3D neighbors only depends on the 11 transition rates, ݎ௠ǡ௠ାଵ, that in turn are larger in the mutant background since the 12 energy barrier decreases. 13 For finding the parameters of the Kolmogorov model, Eq. (19), that better fit in silico 14 tubes and salivary glands we implemented an algorithm that solves, numerically, the 15 set of equations defined by Eq. (19) and the normalization condition σܲ௠ሺݏሻ ஶ ௠ୀଵ ൌͳ 16 to obtain ۃ݉ሺݏሻۄൌσ݉ܲ௠ሺݏሻ ௠. Such algorithm minimizes the error between the 17 curves ۃ݉ሺݏሻۄ obtained in the model and in experiments. 18 The values of the parameters obtained were further used to compare the predicted 19 probability distribution of having ݉ accumulated 3D neighbors for a given value of ݏ: 20 ሼܲ௠ሺݏሻሽ. We evaluated the relative error of this prediction with respect to the actual 21 distribution from data, ܲ௠௔௖௧௨௔௟ሺݏሻ, by computing ߝଶൌଵଶσቀܲ௠௔௖௧௨௔௟ሺݏሻെܲ௠ሺݏሻቁଶ ௠. 22 This quantity is normalized such that in case of the following situation of full 23 disagreement between the distributions, ܲ௠௔௖௧௨௔௟ሺݏሻൌߜ௠ǡ௜ and ܲ௠ሺݏሻൌߜ௠ǡ௝ with ്݅ 24 ݆, provides ߝଶൌͳ (i.e., 100% error). 25 26 Quantitative characterization of spreading in neighbor exchange distributions 27 between apical and basal surfaces 28 40 In order to characterize the spreading away from the diagonal in the neighbor 1 exchange distributions between apical and basal surfaces, e.g., Fig. 2A, we followed 2 the same approach used to quantify intrinsic noise during gene expression processes, 3 see (Elowitz, 2002). Thus, ߟଶൌۃሺ௡ೌି௡್ሻమۄ ଶۃ௡ೌۄۃ௡್ۄ where ۃݖሺ݊௔ǡ݊௕ሻۄൌ 4 σݖሺ݊௔ǡ݊௕ሻ݌ሺ݊௔ǡ݊௕ሻ ௡ೌǡ௡್; ݖ representing any function of ݊௔ and ݊௕ and ݌ሺ݊௔ǡ݊௕ሻ 5 being the probability of neighbor exchange events. We point out that bins in the 6 diagonal do not correspond necessarily to prismatic cells since a fraction of cells can 7 conserved the polygonal class in apical and basal surfaces and yet undergo apico-basal 8 intercalations. 9 10 QUANTIFICATION AND STATISTICAL ANALYSIS 11 Quantification of fluorescence intensity 12 The E-cadherin fluorescence intensity was measured in Fiji by using the Plot Profile 13 tool. We used 3 wt and 3 οEcad representative glands, taking 10 individual 14 measurements for each sample (Fig. S5). We used rectangular ROIs to measure the 15 intensity profiles of lateral cell membrane in the Z-depth where the lumen was visible. 16 In this way, we were able to capture the whole lateral cell membranes from apical to 17 basal. To ensure, a high-quality detection of the cell membrane we developed a 18 maximum Z-projection of those Z-slices where the cell outline of interest and the 19 lumen are clearly visible. Note that the output of the Plot Profile is a 2D plot that 20 ĚŝƐƉůĂLJƐĂĐŽůƵŵŶĂǀĞƌĂŐĞƉůŽƚǁŚĞƌĞƚŚĞܺ axis represents the horizontal distance 21 through the selection (apico-basal cell outline) and the ܻ axis the vertically averaged 22 pixel intensity. 23 Salivary glands measurements 24 We quantified the following geometrical and topological/connectivity descriptors of 25 the segmented salivary glands using a custom-made Matlab code: 26 - Surface ratio (ݏ): Assuming a cylindrical shape for glands, we estimated ݏ by 27 measuring the minimum distance between each cell apical centroid and lumen 28 skeleton (ܴ௔), and measuring ݄, the distance between apical cell centroid and 29 41 cell centroid of an outer cell layer (i.e., basal surface or a concentric layer 1 between apical and basal). Being, the individual surface ratio of a cell, ݏ௖௘௟௟ൌ 2 ோೌ೎೐೗೗ା௛೎೐೗೗ ோೌ೎೐೗೗ , we averaged all the individual cell measurements to calculate the 3 representative ݏ value corresponding to a gland, ݏൌۃݏ௖௘௟௟ۄ. 4 - Cell apical perimeter, lateral surface area, and the cellular volume. 5 - Number of cell contacts: we measured the number of cell neighbors of each 6 cell surface, that is, apical, basal or lateral. In order to remove artefacts, 2 cells 7 must share at least 0.5% of their lateral surfaces area to enable them to be 8 considered as neighbors. ݊௔, ݊௕, and ݊ଷ஽ of each gland were calculated after 9 averaging the number of cells neighbors along the gland. 10 - Percentage of scutoids, average of apico-basal transitions. We quantify the 11 percentage of scutoidal cells that conform the gland and the number of apico- 12 basal transitions in which each cell is involved. 13 14 All the measurements were carried along different concentric radial sections of the 15 glands. We captured the gland thicknesses starting at the apical surface (ݏ௔ൌͳ) and 16 increasing progressively the surface ratio byοݏൌͲǤͷ, until reaching the basal surface 17 (ݏ௕). In this way, the number of the captured radial sections will depend on the ݏ௕ value 18 of each gland. To compare the glands of each phenotype (either wt or mutant) in terms 19 of connectivity related with the surface ratio, we used a maximum radial section 20 common to all the glands. 21 In οEcad mutant glands, due to their phenotype (cells bulge at the basal surface), we 22 removed the bulging tips of cells to quantify the effective surface ratio ݏכ: the 23 maximum value of the surface ratio up to which cells are contacting (that is, 24 ͳ൑ݏ൑ݏכ൑ݏ௕). We noticed that the 3D connectivity of cells is not modified by this 25 approach. To remove the volume of cell tips, we captured all lateral and apical surfaces 26 of cells and we filled each cell volume using the alphaShape Matlab function. 27 28 Voronoi tubular model measurements. 29 We measured the following properties of cells in Voronoi tubular models: area, 30 perimeter and number of sides of cells for a given radial section, and total number 31 48 Petridou, N. I., Corominas-Murtra, B., Heisenberg, C.-P. and Hannezo, E. (2021). 1 Rigidity percolation uncovers a structural basis for embryonic tissue phase 2 transitions. Cell 184, 1914-1928.e19. 3 Pilot, F. and Lecuit, T. (2005). Compartmentalized morphogenesis in epithelia: 4 from cell to tissue shape. Dev Dyn 232, 685ʹ694. 5 Reinhardt, K. (1918). Über die Zerlegung der Ebene in Polygone. 6 Röper, K. (2018). Quantitative Imaging and the Effect of Tissue Topology on 7 Morphogenesis. Dev. Cell 47, 537ʹ538. 8 Rupprecht, J. F., Ong, K. H., Yin, J., Huang, A., Dinh, H. H. Q., Singh, A. P., Zhang, 9 S., Yu, W. and Saunders, T. E. (2017). Geometric constraints alter cell arrangements 10 within curved epithelial tissues. Mol. Biol. Cell 28, 3582ʹ3594. 11 Sanchez-Corrales, Y. E., Blanchard, G. B. and Röper, K. (2018). Radially 12 patterned cell behaviours during tube budding from an epithelium. Elife 7,. 13 Sanchez-Gutierrez, D., Tozluoglu, M., Barry, J. D., Pascual, A., Mao, Y. and 14 Escudero, L. M. (2016). Fundamental physical cellular constraints drive self- 15 organization of tissues. EMBO J 35, 77ʹ88. 16 ^ĄŶĐŚĞnjͲ'ƵƚŝĠƌƌĞnjdŽnjůƵŽŐůƵDĂƌƌLJ:WĂƐĐƵĂůDĂŽzĂŶĚ 17 Escudero, L. M. (2016). &ƵŶĚĂŵĞŶƚĂůƉŚLJƐŝĐĂůĐĞůůƵůĂƌĐŽŶƐƚƌĂŝŶƚƐĚƌŝǀĞƐĞůĨͲ 18 organization of tissues. EMBO J. 35, 77ʹ88. 19 Schindelin, J., Arganda-Carreras, I., Frise, E., Kaynig, V., Longair, M., Pietzsch, T., 20 Preibisch, S., Rueden, C., Saalfeld, S., Schmid, B., et al. (2012). Fiji: an open-source 21 platform for biological-image analysis. Nat. Methods 9, 676ʹ682. 22 Schutgens, F., Rookmaaker, M. B., Margaritis, T., Rios, A., Ammerlaan, C., 23 Jansen, J., Gijzen, L., Vormann, M., Vonk, A., Viveen, M., et al. (2019). Tubuloids 24 derived from human adult kidney and urine for personalized disease modeling. Nat. 25 Biotechnol. 37, 303ʹ313. 26 Shahbazi, M. N., Siggia, E. D. and Zernicka-Goetz, M. (2019). Self-organization of 27 stem cells into embryos: A window on early mammalian development. Science (80-. 28 ). 364, 948ʹ951. 29 Sharma, P., Saraswathy, V. M., Xiang, L. and Furthauer, M. (2019). Delta/Notch 30 49 signaling controls neuroepithelial morphogenesis in the zebrafish spinal cord. 1 bioRxiv 517714. 2 Siedlik, M. J., Manivannan, S., Kevrekidis, I. G. and Nelson, C. M. (2017). Cell 3 Division Induces and Switches Coherent Angular Motion within Bounded Cellular 4 Collectives. Biophys. J. 112, 2419ʹ2427. 5 Spencer, M. A., Jabeen, Z. and Lubensky, D. K. (2017). Vertex stability and 6 topological transitions in vertex models of foams and epithelia. Eur. Phys. J. E 40, 2. 7 Sugimura, K., Lenne, P.-F. F. and Graner, F. (2016). Measuring forces and 8 stresses in situ in living tissues. Development 143, 186ʹ196. 9 Swanson, L. E. and Beitel, G. J. (2006). Tubulogenesis: an inside job. Curr Biol 16, 10 R51-3. 11 Tepass, U., Gruszynski-DeFeo, E., Haag, T. A., Omatyar, L., Török, T. and 12 Hartenstein, V. (1996). shotgun encodes Drosophila E-cadherin and is preferentially 13 required during cell rearrangement in the neurectoderm and other 14 morphogenetically active epithelia. Genes Dev. 10, 672ʹ685. 15 Thompson, D. W. D. (1945). On growth and form. Cambridge university press. 16 Trepat, X., Wasserman, M. R., Angelini, T. E., Millet, E., Weitz, D. A., Butler, J. P. 17 and Fredberg, J. J. (2009). Physical forces during collective cell migration. Nat. Phys. 18 5, 426ʹ430. 19 Tung, J. J., Tattersall, I. W. and Kitajewski, J. (2012). Tips, Stalks, Tubes: Notch- 20 Mediated Cell Fate Determination and Mechanisms of Tubulogenesis during 21 Angiogenesis. Cold Spring Harb. Perspect. Med. 2, a006601ʹa006601. 22 Vicente-Munuera, P., Gómez-Gálvez, P., Tetley, R. J., Forja, C., Tagua, A., 23 Letrán, M., Tozluoglu, M., Mao, Y. and Escudero, L. M. (2020). EpiGraph: an open- 24 source platform to quantify epithelial organization. Bioinformatics 36, 1314ʹ1316. 25 Wetzel, G. (1926). Zur entwicklungsmechanischen Analyse des einfachen 26 prismatischen Epithels. Wilhelm Roux Arch. für Entwicklungsmechanik der Org. 27 Wolny, A., Cerrone, L., Vijayan, A., Tofanelli, R., Barro, A. V., Louveaux, M., 28 Wenzl, C., Strauss, S., Wilson-Sánchez, D., Lymbouridou, R., et al. (2020). Accurate 29 and versatile 3D segmentation of plant tissues at cellular resolution. Elife 9,. 30 50 Yang, X., Bi, D., Czajkowski, M., Merkel, M., Manning, M. L. and Marchetti, M. 1 C. (2017). Correlating cell shape and cellular stress in motile confluent tissues. Proc. 2 Natl. Acad. Sci. 114, 12663ʹ12668. 3 Yang, R., Li, E., Kwon, Y. J., Mani, M. and Beitel, G. J. (2019). QuBiT: a 4 quantitative tool for analyzing epithelial tubes reveals unexpected patterns of 5 organization in the Drosophila trachea. Development 146,. 6 Zallen, J. A. and Zallen, R. (2004). Cell-pattern disordering during convergent 7 extension in Drosophila. J. Phys. Condens. Matter 16, S5073ʹS5080. 8 9 10 11 %R[ A B C D E F G H I <L> σL cell apical perimeter )LJXUH A scutoidal cells 45 6 7 8 polygon class (apical surface) 4 5 6 7 8 polygon class (basal surface) V5 tubular model η=1.2 10-2 . 2η=2.0 10-2 . 2η=2.7 10-2 . 2 na nbnbnb nana s =1.5 bs =2.5 bs =5.5 b probability 0 0.3 Voronoi class (CVT scale) surface ratio (s) spreading coefficient (η ) 2 B na nb Ccross-correlation 0 1 2 3 4 5 6 7 8 9-1-2-3-4-5-6-7-8-9 1 0.8 0.6 0.4 0.2 C(s) surface ratio lag (s) η2 E  A E  V E  AE  L + η2 E  A E  V E  AE  L + )LJXUH )LJXUH )LJXUH )LJXUH 1 Supplemental Information A quantitative biophysical principle to explain the 3D cellular connectivity in curved epithelia Pedro Gómez-Gálvez, Pablo Vicente-Munuera, Samira Anbari, Antonio Tagua, Carmen Gordillo-Vázquez, Jesús A. Andrés-San Román, Daniel Franco-Barranco, Ana M. Palacios, Antonio Velasco, Carlos Capitán-Agudo, Clara Grima, Valentina Annese, Ignacio Arganda-Carreras, Rafael Robles, Alberto Márquez, Javier Buceta, Luis M. Escudero 6XSSOHPHQWDO7H[WDQG)LJXUHV 5 from light to dark blue). In D three Voronoi seeds that correspond to neighboring cells at the apical surface, ݏ= 1, define the triangle ܣܤܥ. Panels E and F track changes in the neighboring relations (accumulated neighbors) of cell ܣ for two increasing values of ݏ: 2 and 4.5 (panels E and F respectively). As shown in D, should a new neighboring cell, ܦ, of cell ܣ appear due to an apico-basal intercalation, then its position must lie inside the vertical parabola defined by the points ܣ, ܤ and ܥ, but outside the circle that these points define (white region), (Lemma 1). Regions accessible to new neighbors are then coded by the green shading in D-F. As ݏ increases, see E, and the cells ܣ, ܤ, ܥ, and ܦ become neighbors, then the parabolas and circles defined by ܣܤܦ and ܣܥܦ restrict the locations of future nearest neighbors. This idea is further reinforced in panel F: winning neighbor ܧ set additional limits to the accessible locations of new neighbors. Thus, the potentiality of a connectivity gain by cell ܣ due to apico-basal intercalations diminishes as the surface ratio increases and eventually becomes null: the number of 3D neighbors of a cell is bounded (STAR Methods). 6 Figure S4. Biophysical model in Voronoi tubes. Related to Figure 4. Applicability of the biophysical model for different cases of the Voronoi model in terms of the CVT scale. Color codes and legends as in Fig. 4D, E. GI H J K A B C D E F 7 Figure S5. Wildtype and οEcad salivary glands features. Related to Figure 5. A: Levels of E-cadh protein expression along the cell apico-basal axis in wildtype (green) and οEcad (orange) glands (݊=30 cells per phenotype) (STAR Methods). The solid lines stand for the mean and the bands for the standard deviation. B-F: Statistical comparison between wildtype (݊=20 glands, green) and οEcad (݊=10 glands, orange) glands properties: cell volume, ݏ௕ כ (effective basal surface ratio), proportion of scutoids, ݊ଷ஽ (number of 3D neighbours) and ݅ (number of apico-basal intercalations), respectively (see Table S1). In each panel, the solid circles represent the average value of the corresponding feature for each gland, long horizontal black lines denote the mean value, and small ones the standard deviation. The protocol for the statistical quantification of data is described in STAR Methods. Significant difference ݌< 0.05 *, ݌< 0.01 **, and ݌< 0.001 ***; No significant difference, ns. G: Density plot of the 3D distribution of neighbor exchanges between apical and basal surfaces as a function of the number of neighbors in apical, ݊௔, and basal, ݊௕, surfaces (as in Fig. 2A) in mutant salivary glands. H: “Poor get richer” principle in οEcad salivary glands: average net gain of neighbors, from apical to basal, as a function of the polygon class in the apical surface. I: Average profiles of the number of apico-basal intercalations, ۃ݅ۄ, average lateral area, ۃܣۄ, and cellular volume fluctuations, ߪ௏ ଶ in in vivo mutant tubes as a function of the apico-basal coordinate, ݏ. J-K: Cross-correlation analysis between energy and intercalation profiles: same as Fig. 3D, E. (݊=10 glands). Apical Basal R Ra βο βd α /α d o A B C 8 Figure S6. Energy cost and bare transition rate in Voronoi tubular models, and scheme of vertex connections of a representative scutoid. Related to STAR Methods. A: Energy cost to gain additional 3D neighbors per cell (top) and the corresponding “bare” transition rate (bottom) in the Voronoi tubular models (V1 to V10, ݏ௕=10) as estimated by the Kolmogorov model. B: A scutoid solid can be simplified as a collection of vertexes (circles) that can be connected by straight edges. The green circle accounts of an apicobasal intercalation point. Alternatively, we can represent the solid as a connected plane graph. In both cases the Euler characteristic is 2. C: In the parameter space ൛ߚҧ௢,ߚҧௗ,ߙௗ/ߙ௢ൟ it is possible to find large regions where the disorder facilitates the accumulation of 3D neighbors (green shaded region) even if the transition rates to gain 3D neighbors in more ordered configurations are larger. The red shaded region indicates in the same parameter space the region where both the transition rates and the accumulation of 3D neighbors are larger in an ordered configuration.