Collective Cell Mechanics of Small-Organoid Morphologies

The study of organoids, artificially grown cell aggregates with the functionality and small-scale anatomy of real organs, is one of the most active areas of research in biology and biophysics, yet the basic physical origins of their different morphologies remain poorly understood. Here we propose a mechanistic theory of small-organoid morphologies. Using a 3D surface-tension-based vertex model, we reproduce the characteristic shapes, ranging from branched and budded structures to invaginated shapes. We find that the formation of branched morphologies relies strongly on junctional activity, enabling temporary aggregations of topological defects in cell packing. To elucidate our numerical results, we develop an effective elasticity theory, which allows one to estimate the apico-basal polarity from the organoid-scale modulation of cell height. Our work provides a generic interpretation of the observed small-organoid morphologies, highlighting the role of physical factors such as the differential surface tension, cell rearrangements, and tissue growth.


Introduction
In less than a decade, organoids became one of the most interesting topics in cells and tissue biology, organogenesis, developmental biology, and study of disease [1,2]. Organoids are aggregates of cells grown in vitro so as to form miniature replicas of a given organ such as intestine, kidney, lung, and brain [3,4]. In most cases, they consist of a closed shell of sheet-like tissue enclosing a lumen, and they have the same microscopic morphology as the mimicked organ itself, e.g., the villus-crypt pattern of the mammalian intestine. Organoids are often grown from embryonic stem cells driven so as to develop a given tissue identity and then cultured in a medium such as Matrigel. They can also be grown from adult tissue-or stem cells [5,6].
One of the defining features of organoids is their physical form. Their initially simple shape transforms into a given morphology as cells divide and mature [7]. The most common organoid morphologies include budded and branched shapes with various spherical or finger-like protrusions, respectively; some branched shapes develop bifurcating networks. The shape of organoids depends on many factors and processes including the intrinsic physical features of individual cells, cell growth and division rates, cell differentiation, etc. So far, theoretical studies of the role of these factors were mostly carried out by representing the tissue as a continuous concentration field in a model of Cahn-Hilliard type [8,9] or as an ensemble of spherical entities using discrete models [7,10]. These approaches account for the biochemical regulation of organoid growth, yielding invaluable insight into, e.g., potential strategies of anticancer therapy [9]. Recently, a 3D vertex model was employed to study how chemical patterning controlling the local cell growth rate may feed back to the mechanics to determine organoid morphology [11]. Together with results addressing morphogenesis of optic cup organoids [12] and transformations of epithelial shells [13], these insights demonstrate that many features of the observed organoid shapes can be interpreted using simple physical models based on mechanisms arising, naturally, from the underlying cell-and tissue-level biological processes.
Here we enhance such a perspective by analyzing a surface-tension-based vertex model of small-organoid shapes which evolve from a spherical shell of cells. To emphasize the collective mechanical effects, we study model organoids consisting of cells of identical properties; in this respect, our organoids resemble tumor and other spheroids where cell differentiation is often absent. The diversity of the obtained shapes is striking, since their formation relies exclusively on the most generic cell-scale mechanics such as the apicobasal differential surface tension. Our model reproduces almost all experimentally observed morphologies and is further interpreted in terms of a theory of epithelial elasticity, which contains several interesting elements such as the curvature-thickness coupling. Furthermore, we explore the formation of in-plane cell arrangements and find that topological defects in cell packing induced by active rearrangements can act as seeds for branching morphogenesis, leading to out-of-equilibrium branched shapes. We also study the formation of organoid shape during tissue growth through successive cell divisions. Overall, these results reveal generic mechanisms of the formation of shape in small organoids of identical cells.

Results
We model organoids as single-cell-thick epithelial shells enclosing an incompressible fluid-like interior referred to as lumen (Fig. 1a). The cells too are assumed to be incompressible, and they carry a surface energy with three distinct tensions at the inner apical side, the outer basal side, and the lateral sides where they adhere to their neighbors [14]. By considering the surface energy alone, we disregard several elements of cell mechanics such as the acto-myosin cable at the apical surface. At the same time, this simplification allows us to work with fewer model parameters which makes the analysis more accessible. We further assume that all cells are identical in terms of their surface tensions and cell volumes V cell . This allows us to explore a more generic scenario where instead of relying on cell differentiation, complex organoid morphologies arise directly from an interplay between the preferred shape of individual cells and long-range interactions due to lumen incompressibility.
We use a 3D surface-tension-based vertex model [15], where cells are represented by polyhedra with polygonal apical and basal sides and rectangular lateral sides (Fig. 1a); the topologies of the apical and the basal cell networks are identical. Cell shapes are parametrized by the vertex positions r j = (x j , y j , z j ), where x j , y j , and z j are dimensionless coordinates expressed in units of V 1/3 cell . The dimensionless energy of the organoid, expressed in units of Γ l V 2/3 cell (where Γ l is the surface tension at the lateral sides) is given by a sum over all N c cells: Here α and β are the dimensionless tensions of the apical and basal sides, respectively, expressed in units of Γ l (Fig. 1a), whereas a cell . In our model, tissue dynamics includes two concurrent processes: (i) deterministic relaxation of the system in the direction of the energy gradient and (ii) cell rearrangements driven by active junctional noise, which allow our model tissue to fluidize. During each time step, the vertices are first moved according to the overdamped equation of motion dr j /dt = −∇ j w; here ∇ j is the gradient with respect to the dimensionless r j , w is the dimensionless energy, and dimensionless time t is measured in units of τ 0 = (Γ l µ) −1 with µ being vertex mobility. Additionally, cells are allowed to rearrange through T1 transitions ( Fig. 1b and Methods) following the model from Ref. [15]. In particular, active noise at cell-cell junctions due to the stochastic turnover dynamics of molecular motors provides energy fluctuations which can drive the system from one metastable state to another via T1 transitions. These transitions are implemented by swapping the four cells arranged around a given lateral side, which is the physical cell-cell junction.
It is convenient to view the junction projected onto the midplane between the apical and the basal surface-and to treat it as an edge of length halfway between those of the corresponding edges on the apical and the basal surface. Since the thus-redefined cell-cell junctions that have a longer length are associated with apical (α) basal (β) lateral lumen (vlumen)  a higher energy barrier between the two metastable states, they are less likely to undergo the T1 transition. To describe this, we employ the previously proposed threshold-based scheme where the probability of a T1 transition in junctions shorter than the threshold δl is unity, whereas junctions longer than the threshold undergo the transition with a probability k T1 δt/E [15]. Here k T1 is the rate of active T1 transitions measured in units of τ −1 0 , δt = 10 −4 is the time step, and E is the total number of cell-cell junctions. We choose δl = 0.15 (defined as a dimensionless quantity expressed in units of V 1/3 cell ) which is between a quarter and a third of the average junction length; we stress that the results do not depend significantly on this choice ( Supplementary  Fig. 4b). As shown previously, tissues that are subject to this type of active noise behave like viscous fluids with activity-dependent stress-relaxation time scale and the associated effective viscosity [15]. Therefore, an important aspect of the organoid shapes could be their sensitivity to the changes of the active T1 rate k T1 .
In most cases discussed here, the model parameters include only a handful of dimensionless quantities: The apical and basal surface tensions α and β, the volume of the lumen relative to the cell volume v lumen = V lumen /V cell , the total cell number N c , and at most two parameters related to junctional activity.

Phase diagram
We start by studying the shapes of organoids with fixed cell number N c = 300 and lumen volume v lumen = 100 and we first consider the case with no active T1 transitions; this limit is denoted by k (0) T1 = 0 for consistency of notation, with the meaning of k (0) T1 explained in the next paragraph. The initially spherical shape evolves into three characteristic types of morphologies depending on the values of α and β: stomatocyte (cup-like), budded, and spherical ( Fig. 1c and Supplementary Movie 1). In the phase diagram in the (α, β) plane ( Fig. 1e and Supplementary Fig. 1a), spherical shapes occupy the region where the tissue tension defined by α + β is larger than about 1.9, budded shapes emerge at α larger than about 1 and β smaller than about 0.6, whereas stomatocytes occur in the region where α + β is smaller than 1.9 but larger than around 0.9 and α is less than about 1. In many shapes, the characteristic features are well-developed so that the classification is unambiguous, but this is not always the case. For example, a generally spherical model organoid with small buds may be regarded either as a spherical or as a budded shape. Shapes found at the boundary between the stomatocyte and the budded-organoid domain contain features of both morphologies and are therefore best referred to as hybrid ( Supplementary Fig. 1c). Finally, at tissue tensions α + β that are smaller than about 0.9 we find model organoids containing one or more cells of anomalous, non-physical shape, or else the organoids start to self-overlap.
We next analyze the case where relaxation is accompanied by an initially high rate of active cell rearrangements with k T1 (t = 0) = k here t max = 1000 is the total simulation time. The finite duration of junctional activity mimics the natural course of morphogenetic events, which, guided by biochemical signals, begin and end at a predefined time; the linear temporal profile is just a simple example of a gradual decrease of activity. In Section "Active cell rearrangements" we show that this choice does not significantly affect the results. Like in the case with no activity, we obtain the budded, stomatocyte, and spherical shapes. However, the active T1 transitions also give rise to the branched morphologies appearing at the apical and basal tensions somewhat smaller than unity, that is in part of the region occupied by stomatocytes if the junctional activity is absent (Figs. 1d and f, Supplementary Figs. 1b and 2, and Supplementary Movie 2). Like before, the phase diagram also contains hybrid shapes that feature either multiple invaginations, both invaginations and evaginations, or are devoid of any clear features ( Supplementary Fig. 1d); these shapes appear at α ≈ β and tissue tension α + β between about 1 and 2. The domain of non-physical shapes is somewhat larger than in absence of active T1 transitions.
The results obtained using the N c = 300 organoids are quite representative of our model. This is demonstrated by the four sets of shapes in Supplementary Fig. 3 which show that with an appropriate rescaling of the preferred lumen volume v lumen (Methods), the type of organoid morphology does not change significantly with the cell number N c even if N c is halved or doubled. On the other hand, the fixed lumen volume constraint does play an important role in the model morphologies. This constraint is implemented using an auxiliary harmonic volumetric energy term, which ensures that the actual lumen volume does not depart from the preferred value. If the modulus K lumen of the auxiliary term is small enough compared to the typical values of α and β (for example, 0.001) then the actual and the preferred volumes differ considerably and the model organoids lack clear features and are more or less round ( Supplementary Fig. 4a).

Origin of shape
Aside from the role of active T1 transitions discussed in Section "Active cell rearrangements", the main mechanisms responsible for the formation of the distinct morphologies are (i) the incompatibility of the preferred surface area of the organoid and the enclosed lumen volume and (ii) the apico-basal differential tension. In particular, the preferred cell width-to-height ratio scales with the tissue tension α + β as ∼ (α + β) −1 so that small values of α + β drive the columnar-to-squamous transition [15]. In turn, this increases the organoid surface area at a fixed lumen volume, causing the formation of complex non-spherical morphologies. A similar mechanism was previously studied within a 2D model of tissue cross section, yielding shapes resembling our 3D small-organoid morphologies [16].
The area-volume incompatibility can be conveniently quantified by the reduced volume defined by v = 6 √ πV Here A is the area of the midplane surface located halfway between the basal and the apical side of the organoid and V is the corresponding enclosed volume; v = 1 for a sphere, and its value decreases with increasing area at constant volume. We compute the reduced volumes of all 1328 shapes used to construct the phase diagrams in Figs. 1e and f and find that they depend strongly on the tissue tension α + β ( Fig. 1g and Supplementary Fig. 5a). The numerically obtained dependence agrees well with the analytical prediction where we estimate the organoid midplane surface area by considering a flat epithelium of identical cells with regular hexagonal basal and apical sides [15]. From the force balance along cell height h, ∂w/∂h = 0, and the relation between (dimensionless) cell volume 1, height h, and area a, which reads 1 = ha, we can calculate the equilibrium cell height h 0 = 2 1/3 /3 1/6 (α + β) 2/3 and midplane area a 0 = 3 1 where the volume enclosed by the midplane v midplane = 223 (expressed in units of V c ) is the average taken from the simulated shapes. The relation between v and α + β given by Eq. (3) is plotted by the solid cyan line in Fig. 1g. While the apico-basal differential tension α − β has a rather limited effect on the reduced volume v as witnessed by the limited spread of points in Fig. 1g, we find that it is often crucial for the formation of budded and stomatocyte morphologies. This is because the buds consist of cells with the apical sides smaller than the basal sides whereas invaginations in the stomatocyte shapes require cells where the apical sides are larger than the basal ones. The former cell type is energetically preferred at α − β > 0 whereas the latter is favored at α − β < 0.
An insightful way of interpreting our small-organoid shapes builds on the analogy with lipid vesicles characterized by spontaneous curvature [17]. In particular, (i) both organoids and vesicles are closed shells with a well-defined surface area and enclosed volume as well as fluid-like in-plane order in the case of organoids with junctional activity and (ii) the apico-basal polarity due to the differential tension α − β in tissues leads to a preferred conical shape of cells [14,18] so that the tissue has a certain spontaneous curvature reminiscent of that in asymmetric lipid membranes, e.g., due to inclusions. Although organoids are thick rather than thin shells like vesicles, this comparison is quite illuminating because of both similarities and differences, especially within the context of the theory of elasticity developed in the penultimate section of this paper. As far as the characteristic morphologies are concered, we note that the phase diagrams of organoids and vesicles feature stomatocytes at negative spontaneous curvatures and differential tensions α − β, respectively, whereas the prolate and pear-like vesicles at positive curvatures resemble budded organoids found at positive differential tensions ( Supplementary Fig. 6).
The shapes shown in Fig. 1d as examples of active model organoids are characteristic of a given set of model parameters but they are not unique due to the stochasticity of active cell rearrangements. The resulting shape variability is illustrated by Supplementary Fig. 5c, which shows three branched organoids at α = 0.7 and β = 0.5. All of these shapes have five unevenly articulated branches of somewhat different lengths and diameters. To systematically explore this variability, we simulate 300 instances at the four pairs of α and β from Fig. 1d and we compute the corresponding distributions of reduced volumes. We find that shape variability differs among the four shapes, the spherical and the budded one being the least and the most variable, respectively (Supplementary Figs. 5d-g). The variability of the budded shape is likely due its location in the phase diagram: It lies close to the boundary between the budded and the branched morphologies, and thus in some instances the shape of the buds may be considerably more pronounced than in others. Nevertheless, throughout the phase diagram the standard deviation of the distribution never exceeds ≈ 0.04, which shows that the reduced volumes of the organoid shapes are well-defined ( Supplementary  Fig. 5b). Note that since the budded and the stomatocyte morphologies appear at similar values of α + β, their reduced volumes v are similar as well ( Supplementary Fig. 5a).

Active cell rearrangements
Our active model organoids explore the energy landscape at a fixed cell number by actively rearranging the cells and dissipating the energy through friction with the environment [15]. Due to the complexity of the energy landscape, there is no guarantee that the system reaches the global energy minimum by the end of simulation. Instead, due to the decreasing active T1 rate k T1 , which gradually drops to zero in each simulation, it is plausible that the final small-organoid morphologies are trapped in local energy minima. If true, this would explain why the branched morphologies are only found in our simulations that include active cell rearrangements (Fig. 1f). To test this, we first vary the initial active T1 rate k (0) T1 at fixed α = 0.7 and β = 0.5, which give a branched morphology at k (0) T1 = 200. Importantly, we find that the branched morphology develops only when k (0) T1 is high enough and we find that it is associated with a significantly higher energy compared to the corresponding energy-minimized shape ( Figs. 2a and b). This result shows  that the branched small organoids are inherently out-of-equilibrium shapes that require a sufficiently high degree of junctional activity to form. Furthermore, being trapped in the respective local energy minima, the shapes of branched organoids may significantly depend on the details of the relaxation protocol.
To challenge this possibility, we next replace the linear temporal profile of the active T1 rate by a step-like profile where the initial period with a constant T1 rate k T1 = k (0) T1 > 0 of duration t = 500 is followed by an equally long period with no T1 transitions (k T1 = 0). We again find that the branched morphologies only occur if k (0) T1 is high enough during the initial period and the final shapes stay qualitatively the same as before (Figs. 2e-g). We note that the branched morphology is formed soon after the beginning of the active phase as illustrated by the sequence of shapes in Supplementary Fig. 7. Thus neither the duration nor the precise temporal profile of the activity seem essential for the development of the branches, provided, of course, that the duration is long enough.
We further test how the formation of the branched morphology depends on the implementation of junctional activity by considering an entirely different numerical scheme of active junctional noise. Here T1 transitions arise directly from fluctuations in tensions along cell-cell junctions instead of the active-T1 scheme employed above. To this end, we extend the energy function [Eq. (1)] by a line-tension term describing fluctuations of tensions at cell-cell junctions, i.e., along the direction perpendicular to the apicobasal axis. Here the sum runs over all cell-cell junctions E; l ai and l bi are lengths of the apical edge and the basal edge, respectively, that correspond to the i-th junction, whereas γ i (t) is the effective time-dependent line tension associated with these edges. We then assume that the fluctuations of the line tension γ i (t) are described by the Ornstein-Uhlenbeck process [19,20] so that Here τ is the time scale associated with the turnover of molecular motors (set to unity, τ = 1) and where σ 2 is the long-time variance of the tension fluctuation. If the magnitude of fluctuations σ is sufficiently large, individual edges occasionally shrink to zero length, initiating a T1 transition (Methods). We analyze two variants of the tension-fluctuation scheme of T1 transitions. In the first one, σ decreases linearly from some initial value σ (0) to zero (Figs. 2h-j) whereas in the second it has a step-like temporal profile where a period of constant tension fluctuations with σ = σ (0) is followed by an equally long period with σ = 0 (Figs. 2k-m).
In agreement with our initial scheme of active noise, branched morphologies are formed in both variants provided that the activity parameter σ (0) is large enough (Figs. 2j and m).
In order to better understand the formation of branches, we examine the in-plane tissue structure represented by the topology of the polygonal network of cell-cell junctions. We find that the branched shapes, obtained in the more active organoids, develop very different cell arrangements compared to the less active non-branched shapes (Fig. 2a). In particular, the branched shapes contain many pentagonal and heptagonal cells, which accumulate at the tips of the branches and at their bases, respectively. The distribution of these cells within the organoid can be quantified by plotting their pair correlations as functions of the topological distance d defined as the integer shortest path between two cells; d = 1 for the nearest neighbors, d = 2 for the next-nearest neighbors, etc. (Fig. 2c and Methods). For pentagons, we observe a bimodal distribution with a small peak at d = 2, which corresponds to the distance between pentagons within the same branch, and a large peak at d = 10, which corresponds to the distance between pentagons in different branches. In contrast, the distribution of heptagons, which are mostly located at the bases of the branches, consists of a single peak at d ≈ 6. Furthermore, since the distance between the bases is shorter than that between the tips, heptagons are more likely to be found closer to other heptagons than pentagons are to other pentagons. Note that these pair correlations are qualitatively similar in the budded shapes, which favor pentagons at the buds and heptagons in the saddle-like parts. In contrast, stomatocytes and spheres lack such features and are thus associated with unimodal distributions (Supplementary Fig. 8).
These results suggest that defects in cell packing, established through active cell rearrangements, may be crucial for the formation of the overall shape. We hypothesize that the observed characteristic distributions of pentagonal and heptagonal cells in the branched morphologies are established hand in hand with the shape itself. In particular, the high-activity conditions enable the formation of a temporarily biased distribution of defects in cell packing, in which pentagons and heptagons act as a seeds for positive and negative Gaussian curvature, respectively ( Fig. 2d and Supplementary Movie 3). In turn, once the distribution of the Gaussian curvature is established across the organoid, pentagons become locked at the positive-Gaussian-curvature tips, whereas heptagons are bound to the saddles. As a result, this interplay between the in-plane cell organization and the Gaussian curvature preserves the overall organoid shape and thus causes it to become trapped away from the global energy minimum. In contrast, organoids that are never exposed to a high activity cannot develop the required distribution of pentagons and heptagons, which prevents the formation of multiple localized regions of positive Gaussian curvature necessary to develop branches.
Surprisingly, even though the relation between the local Gaussian curvature of the surface and locations of the defects in packing of constituents is well known in general [21], its potential relevance for organoid (and tissue) development has not yet been pointed out. Furthermore, it is possible that this mechanism plays an important role in cell turnover. Indeed, in one of the fastest renewing organs in mammals-the intestinecells divide at the bases of the villi and are extruded at the tips [10]. It could be that cell extrusion is initiated by the presence of cells with pentagonal bases. These cells may accumulate at the tips of the villi where their basal side is larger than their apical side simply because pentagonal cells are energetically favorable at positive Gaussian curvatures ( Fig. 2d and Supplementary Fig. 9a). With their truncated-pyramid shape, such cells could conceivably more easily delaminate from the tissue due to mechanical reasons. In contrast, cells located at the bases of the villi where the Gaussian curvature is negative are preferentially heptagonal. In these cells, the local saddle-like shape of the tissue is consistent with somewhat anisometric shapes of the apical and the basal side, with the two long axes roughly perpendicular to each other ( Fig. 2d and Supplementary Figs. 9b-d). As the areas of the two sides should be roughly the same, cells of such a form are stabilized against extrusion and act as anchors. Overall, this hypothesis suggests that the curvaturedependent extrusion rate, observed in intestinal epithelia [22], could be to a large extent of mechanical origin. Furthermore, our findings suggest that the basis of this mechanism in organoids may be established spontaneously during shape formation.

Curvature-thickness coupling
To better understand the obtained shapes, we recast the cell-level energy Eq. (1) in the continuum limit following Ref. [18]. We approximate cell shape by a truncated right square pyramid, a body with in-plane isometry, imposing a fluid-like response to local shear stress. Cell shape is parametrized by the local principal curvatures c 1 and c 2 of the midplane as well as by tissue thickness h, i.e., cell height, measured in units of 1/V 1/3 c and V 1/3 c , respectively (Supplementary Note 1). In the first-order approximation, the dimensionless elastic energy per unit area reads (6) The first two terms together represent the surface tension whereas the third and the fourth term are the local and the Gaussian bending energy per unit area, respectively. Note that the thickness h in Eq. (6) is an independent variable: While it is constant in a flat tissue at α = β, in the nontrivial shapes h is coupled to the local curvature. As a result, the surface tension varies along the surface and so do the local bending modulus √ h/8, the spontaneous curvature c 0 = 2 √ h(α − β), and the Gaussian modulus [i.e., the last term in Eq. (6) divided by the Gaussian curvature c 1 c 2 ]. These features render the tissue energy functional quite different from the usual surface and bending energies of solid or fluid membranes.
The coupling between the mean curvature and the thickness (β − α)(c 1 + c 2 )h/2, contained in the local bending term, is the most important of the above effects. In the more elaborate shapes, this coupling is quite pronounced as illustrated by Figs. 3a and b showing the thickness-curvature portraits of a budded and a stomatocyte shape, respectively. In these two diagrams, each cell is represented by a point indicating the mean curvature c = (c 1 + c 2 )/2 and thickness h of the tissue at this cell (Methods). The positive correlation between thickness and curvature is evident in the budded shape (Fig. 3a) as is the negative correlation in the stomatocyte (Fig. 3b). In the budded α = 1.1, β = 0.5 shape with a positive spontaneous curvature c 0 (i.e. a spontaneous curvature with the same sign as the curvature of a sphere), the taller cells are located at the buds where the mean curvature is largest (Fig. 3a), whereas in the α = 0.5, β = 1.1 stomatocyte where c 0 < 0, they are concentrated at the invagination (Fig. 3b).
To quantify the curvature-thickness coupling throughout the phase diagram, we introduce the relative tissue thickness modulation δ = p · ∆h/h, where ∆h is the amplitude of thickness modulation whereas h is the average cell height in the model organoid; the prefactor p takes the value of +1 if curvature and thickness are correlated (Fig. 3a) or −1 if they are anticorrelated (Fig. 3b). We find that the thickness modulation is sizable in all non-spherical shapes, with |δ| reaching about 0.65 in stomatocytes and 0.35 in budded and branched shapes (Fig. 3c). Consistent with the continuum theory [Eq. (6)], the sign of δ agrees with the sign of the modulus of the curvature-thickness coupling α − β.
We can further analyze the coupling between curvature and thickness by fitting the h(c) plots for individual model organoids by linear functions. We then plot the vertical intercept n of these functions against tissue tension α + β and the slopes k against the differential tension α − β for all 185 model organoid shapes discussed here (Supplementary Fig. 10). We find that the vertical intercept n can be described by the above analytical solution for hexagonal cells in a flat epithelium n ≈ 2 1/3 /3 1/6 (α + β) 2/3 , whereas the slope k can be approximated by λ 1 + (α − β)λ 2 where λ 1 = −0.060 ± 0.004 and λ 2 = 0.45 ± 0.01 are fitting parameters. (These two observations only apply to organoids with α + β < 1.9. Beyond this threshold, organoids approach a spherical shape and their midplane area is defined by the enclosed volume rather than by α and β.) We can then rescale h → h − 2 1/3 /3 1/6 (α + β) 2/3 and c → [λ 1 + (α − β)λ 2 ] c, which allows us to collapse all of the h(c) scatter plots for the model organoids at k (0) T1 = 200 with α + β < 1.9 (Figs. 3d-e); the same can also be done for the organoids at k (0) T1 = 0 in Fig. 1e. While the modulation of cell height is often attributed to cell differentiation [23,24], our model shows that it can appear even in tissues of identical cells where it is a telltale sign of mechanical apico-basal polarity. Similar shape features are often observed in the folded morphologies in a variety of tissues in both vertebrates and invertebrates [14].

Discussion
While our model yields a diverse catalog of possible small-organoid morphologies, their formation in real systems may rely on processes other than surface tension and active cell rearrangements. In particular, the shape usually develops while cells in an organoid grow and divide. To test how shapes are established during growth and how they compare to those predicted by our original approach, we generalize our model by including volumetric cell growth and cell division (Methods). We simulate growing organoids at fixed cell-growth and cell-division rates, τ d = 2000 being the expected time until a cell next divides, whereas the lumen volume follows Eq. (7) (Methods). We start from a sphere of 100 cells and we monitor the growing organoid until the number of cells reaches 300. We first compute the shapes in absence of activity (i.e, k (0) T1 = 0, where only spontaneous T1 transitions are allowed) and find that in terms of their final reduced volumes, the thus obtained shapes agree well with those from our original model ( Fig. 4a and Supplementary Figs. 5d-g). The spherical, budded, and stomatocyte shapes all develop their characteristic features, although they are somewhat less regular than in their non-growing counterparts due to a different relaxation mechanism (Figs. 4b-d and Supplementary Movie 4). In contrast, the α = 0.7, β = 0.5 shape does not develop branches although its reduced volume is similar to that of the branched shapes obtained in active organoids at a fixed cell number (Figs. 4a, e, and f). These results show that cell-scale activity due to cell divisions is not sufficient to form branched morphologies.
To check if branching can be recovered in growing organoids by including junctional activity, we combined cell division with the previously used active T1 transitions with a linearly decreasing T1 rate [k T1 (t) = Our results show that small organoids can form a variety of shapes relying on a simple collective mechanics of non-differentiated cells, coupled through a fluid-like organoid interior. We explore a scenario where a complex shape forms due to a mismatch between the organoid's preferred surface area and the volume of the enclosed lumen. The four characteristic morphologies-spherical, stomatocyte, budded, and branched-are distinguished by the reduced volume. We also develop a continuum theory and use it to explain the fine features of the shapes (specifically, the organoid-scale cell-height modulation) which could be measured in vitro to extract the apico-basal surface-tension polarity directly from the shape. Overall, our model provides clues for the control of organoid shapes in experiments. This can be done by changing the relative apical and basal tensions of cells, which should depend on the biochemical composition of the environment (i.e., the contents of the lumen and the medium). These tensions are best viewed as effective concepts encompassing all membrane-related effects and processes that shape the cells, and may be altered in several ways including, e.g., the activation of the complex cascade leading to contraction of the apical network [25]. Our results also show that the organoid shape may be affected by suppression or promotion of junctional activity and, importantly, its effect on the topology of the tissue represented by the polygonal tiling seen on the basal or apical plane. A systematic assessment of the different experimental options could conceivably expand the toolbox of the organoid technology and tissue engineering [26].
A natural extension of this work would be to include planar cell polarity, which can induce tubulogenesis and thus branching [27], and cell differentiation, which could contribute to self-organization of cells within organoids through differential adhesion as well as to programmed shape formation through spatial modulation of the apico-basal polarity. To better describe the typical experimental setup in growing organoids, the model could consider tissue growth on curved solid substrates and include viscous dissipation at the tissue-matrix interface, which may promote buckling and branching [28]. In such a case, other cell-scale active processes, e.g., cell motility, could contribute to the collective cell dynamics and the formation of the overall organoid morphology.

Implementation
Simulations of organoid shapes are performed within vertX3D, a package of C++ routines that implement the 3D vertex model of epithelial tissues. The initial configuration of cells is a spherical shell enclosing the lumen volume v lumen . These configurations are constructed by finding the dual of the convex hull to the minimal-energy solution of the Thomson problem on a sphere [29]. The dual determines the positions of the apical vertices, whereas the positions of the basal vertices are calculated by displacing their apical counterparts radially by [3(v lumen To correct the initial cell volumes, which do not agree exactly with the preferred volumes, the organoids are first relaxed at k T1 = 0 for a short time before the active T1 rate is set to k (0) T1 . The preferred cell volumes and the lumen volume are enforced by auxiliary harmonic terms in the energy; the modulus of these terms is very large (100 unless stated otherwise).
In the active-T1 scheme, the T1 transitions are initiated as described in the main text; in the fluctuating tension scheme, they are only initiated if the length of a junction expressed in units of V In organoids with no active T1 transitions (k (0) T1 = 0) simulations are run until t = 2000, i.e., twice as long as in organoids with junctional activity; this is because the stomatocyte shapes at α ≈ β require a longer time to fully develop. In our vertex model, steric repulsion between cells is not implemented and self-overlapping shapes are therefore possible. During post-processing, we check each obtained final shape and discard all that selfoverlap; temporary overlaps may occur during the simulation.
When analyzing the obtained organoid shapes, we compute cell height (i.e., tissue thickness) reported in Fig. 3 as the distance between the centroids of the apical and the basal side, and we approximate the mean curvature of each cell by that of a truncated cone with the same apical area, basal area, and height.
When considering growing organoids, cells that are otherwise subject to the fixed-volume constraint enter the growth perido with a constant probability dP/dt = 1/τ d , where τ d = 2000 is the characteristic time scale of division events. During the growth period, the cell volume is increased linearly in time according to dV /dt = 1/τ g , where τ g = 1 is the cell-growth time scale. As soon as the volume is doubled, the cell is divided with the mitotic plane connecting a random pair of opposing lateral sides. Starting with a spherical organoid of 100 cells, we let the cells divide while increasing the preferred lumen volume according to Eq. (7), substituting N c by the total volume of the tissue.

Volume and active T1 rate rescaling
When comparing morphologies with different cell numbers N c in Supplementary Fig. 3, the active T1 rate k T1 and lumen volume v lumen must be appropriately rescaled for the comparison to be relevant. As the probability for an individual edge to undergo a T1 transition is k T1 δt/E, the most relevant comparison is between model organoids where an individual cell-cell junction has the same probability of undergoing a T1 at both N c in question. We therefore rescale k T 1 → k T 1 E Nc /E 300 ; here E Nc is the number of cell-cell junctions in an organoid of N c cells and E 300 is their number in a 300-cell organoid. The same rescaling is also used when considering growing organoids (Fig. 4) so that edges have an equal probability of undergoing an active T1 transition regardless of the number of cells.
When rescaling lumen volumes, we choose them such that the model organoids have the same α, β threshold for the spherical shape, which can be estimated as follows. As before, we approximate cell height by h 0 = 2 1/3 /3 1/6 (α + β) 2/3 . The threshold for a spherical shape is then the value of tissue tension α + β at which a spherical shell of thickness h 0 around the lumen of volume v lumen (N c ) has volume N c . We then equate the α + β values at the given target N c and the reference N c = 300, giving the condition Here v lumen (300) = 100 is the lumen volume at 300 cells, whereas v lumen (N c ) is the lumen volume at the target cell number N c . This equation can trivially be solved numerically for v lumen (N c ).

Topological pair correlations
The topological pair correlation function g γ (d) encodes how many γ-sided cells may on average be expected at a topological distance d from another γ-sided cell. It is defined as Here the sum runs over all γ-sided cells and N γi (d) is the number of γ-sided cells at a topological distance d from the i-th γ-sided cell. The normalizing prefactor n γ is the total number of γ-sided cells in the organoid. The expression in the brackets is averaged over 300 simulation runs.

Shape anisometry
The long axes of the apical and basal sides of cells in Fig. 2d and Supplementary Fig. 9 are obtained by diagonalizing the gyration tensor computed from the vertices of a side. The anisometry factor of a cell side κ is given by κ = (g 1 − g 2 )/(g 1 + g 2 ), where g 1 and g 2 are the largest and the second largest eigenvalue of the gyration tensor, respectively; as the cell sides are slightly non-planar, the third eigenvalue is finite rather than 0 but still much smaller than g 1 and g 2 .

Data availability
The data containing the computed organoid shapes is available from the corresponding author on reasonable request.

Code availability
The codes used to compute the organoid shapes is available from the corresponding author on reasonable request.    While the initial t = 0 shape is spherical, the t = 250 shape shows that branches are formed rather quickly during the active period. The similarity of the t = 250 and the t = 500 shapes suggests that the duration of the active period is not very important as long as it is long enough, and that further extending it would not affect the shapes very much.  Fig. 1f of the main text. Points with α + β > 1.9 where the model organoids are spherical are denoted by hollow markers. The red line shows the function 2 1/3 /3 1/6 (α + β) 2/3 , its dashed section corresponding to α + β > 1.9. b Slopes k for the same organoids as a function of differential tension α − β. Hollow markers are again for points with α + β > 1.9. Blue line shows the linear fit of the type λ 1 + (α − β)λ 2 .

Supplementary Note 1 3D continuum theory
The derivation of the continuum theory is based on generalizing the approach reported in Ref. [2] to three dimensions. We begin by defining a model cell represented by a truncated pyramid with a rectangular base and a volume of V cell . The model cell is defined by the principal curvatures C 1 and C 2 , height H, and the aspect ratio in the midplane located halfway between the apical and the basal side defined by where ∆s 1 and ∆s 2 are the length and the width of the rectangular midplane cross section, respectively (Supplementary Fig. 11). These four parameters can describe four deformation modes: (i) and (ii) bending in the directions of the principal curvatures, (iii) changes in the local tissue thickness, and (iv) pure shear deformation. The geometrical parameters of the model cell can be related to the length and the width of the apical side ∆a 1 and ∆a 2 , respectively, and the length and the width of the basal sides ∆b 1 and ∆b 2 , respectively, by and ∆a 2 = 1 − C 2 H/2 Upon inserting the relations ∆s 1 = (∆a 1 + ∆b 1 )/2 and ∆s 2 = (∆a 2 + ∆b 2 )/2 into Eq. (1) and using Eqs. (2) and (3), we obtain ∆b 2 = 1 ν 1 + C 2 H/2 Lastly, by combining the expression for the midplane area ∆A m = ∆s 1 ∆s 2 with the above relations, we find that Equations (2)-(5) allow us to express the lengths of edges of the apical and basal sides in terms of C 1 , C 2 , ν, H, and ∆A m .
We now wish to derive the elastic energy per unit area. For a section of tissue with midplane area ∆A m at fixed C 1 , C 2 , ν, and L, the energy reads where Γ a , Γ b , and Γ l are the apical, basal and lateral surface tensions, respectively, ∆A a and ∆A b are the apical and basal areas corresponding to ∆A m , whereas dW l /dA m is the energy density associated with the lateral sides. The energy terms associated with the apical and the basal sides can be rewritten by using ∆A a = ∆a 1 ∆a 2 , and ∆A b = ∆b 1 ∆b 2 along with Eqs. (2)-(5) such that The derivation of the term associated with the lateral sides is more involved. We first calculate the energy of the lateral sides of a single cell with a volume of V cell and fixed C 1 , C 2 , ν, and H. This energy reads where A l1 is the area of the lateral side with edges ∆a 1 and ∆b 1 whereas A l2 is the area of the lateral side with edges ∆a 2 and ∆b 2 . By using simple trigonometry and relations derived above, A l1 and A l2 can be recast as and A l2 = 1 respectively. Now we combine these results with the expression for the volume of a single cell with given ∆a 1 , ∆a 2 , ∆b 1 , ∆b 2 , and H, which reads V cell = ∆a 1 ∆a 2 + 1 2 ∆a 1 ∆b 2 + 1 2 ∆b 1 ∆a 2 ∆b 1 ∆b 2 H 3 = 1 + C 1 C 2 H 2 /12 (1 + C 1 H/2) and we obtain we find that dW l dA m = Γ l 1 + 1 12 Upon inserting this result along with Eq. (7) into Eq. (6), we obtain the final energy density of a tissue section of midplane area ∆A m . To recast it in dimensionless form, we divide it by Γ l V and c 0 = 2(α + β) 1/3 (α − β) The curvature-thickness coupling term (β − α)(c 1 + c 2 )h/2 in Eq. (18), which gives rise to the spontaneous curvature in Eq. (19), can also be interpreted in a different manner. We again consider a fluidized tissue with ν = 1 and we assume that the tissue forms a closed shell as in the main text. We then integrate the energy per unit area over the entire surface, which gives If h does not vary across the shell, then the integral A (c 1 + c 2 )h dA is equal to the difference between the apical and the basal area of the shell ∆A, and we have where f is a shape-dependent function of total tension α + β but not of differential tension β − α. From the form of this expression we see that the differential tension β − α is coupled to ∆A and thus (β − α)/2 effectively acts as an area-difference tension. In this respect, the continuum theory of our model tissues derived above departs from the area-difference-elasticity theory of lipid vesicles where the so-called non-local bending energy gives rise to a term ∝ (∆A − ∆A 0 ) 2 [3].
As an estimate of the error arising from rectangular rather than a hexagonal base, our continuum theory gives the height of a flat epithelium as (α + β) 2/3 , whereas the exact result for a flat epithelium of hexagonal cells is 2 1/3 /3 1/6 (α + β) 2/3 , a difference of only about 5%; presumably, a similar relatively small scale of error can be assumed for c 1 and c 2 , and all terms are higher than first order.

Supplementary Movies
Movie 1. Evolution of shapes of the three model organoids in Fig. 1c of the main text without active T1 transitions (k Note that while in the spherical (α = 1.2, β = 1.2), stomatocyte (α = 0.5, β = 1.1), and budded (α = 1.1, β = 0.5) organoid the shapes resulting from growth and from junctional activity (Fig. 1d of the main text) belong to the same category, the (α = 0.7, β = 0.5) growth-induced shape differs from the branched morphology in Fig. 1d as branching requires the presence of active T1 transitions, which are absent in the growing-organoid model.