Morphological instability at topological defects in a three-dimensional vertex model for spherical epithelia

Epithelial monolayers are a central building block of complex organisms. Topological defects have emerged as important elements for single cell behavior in flat epithelia. Here we theoretically study such defects in a three-dimensional vertex model for spherical epithelia like cysts or intestinal organoids. We find that they lead to the same generic morphological instability to an icosahedral shape as it is known from spherical elastic shells like virus capsids, polymerized vesicles or buckyballs. We derive analytical expressions for the effective stretching and bending moduli as a function of the parameters of the vertex model, in excellent agreement with computer simulations. These equations accurately predict both the buckling of a flat epithelial monolayer under uniaxial compression and the faceting transition around the topological defects in spherical epithelia. We further show that localized apico-basal tension asymmetries allow them to reduce the transition threshold to small system sizes.


Introduction.
Epithelial monolayers are a central element of the architecture of complex organisms.They separate different compartments, can form highly convoluted shapes and have exceptional mechanical properties.In particular, they can quickly undergo transitions between fluid-like and elastic properties [1,2], driven e.g. by cell density or the aspect ratios of single cells.In general, the properties of the single cells are essential to understand the physical properties of epithelial monolayers.Topological defects defined by the neighborhood relations of the single cells have emerged as especially important elements for transformations in epithelial monolayers [3].For example, it has been found that single cells tend to be extruded at such topological defects [4].While topological defects are a natural ingredient of hydrodynamic theories [5], it is challenging to include them in elastic descriptions of monolayers [6][7][8].
Here we show that spherical epithelia like cysts or intestinal organoids are a natural starting point to study the global effects of topological defects in epithelial monolayers.They are experimentally very accessible and of large biomedical relevance [9].Due to Euler's polyhedron theorem, they necessarily have to include twelve pentagons in the close-packed tiling of the spherical surface [10].In order to combine these topological defects with the typical mechanical properties of epithelial monolayers, we employ a three-dimensional (3D) vertex model (VM), in which cells are described as polyhedra with a fixed volume and with polygonal faces contributing to the total energy through surface tensions [11].The 3D VM has been used before for modeling spherical epithelia [6,[12][13][14][15], but coarse-graining procedures have not been able to fully address the role of topological defects in such a setting.
By simulating the 3D VM for complete spherical epithelial shells, we discovered an icosahedral faceting instability.While small shells with few cells have a spherical shape (Fig. 1 conical instabilities at the pentagonal defects (Fig. 1(b)) [16].This transition is well known for two-dimensional elastic crystals [10,17], including virus capsids [10,18], polymerized vesicles [19,20], and buckyballs [21], but has not been described before for spherical epithelia.Our numerical results suggest that a continuum limit exists for the 3D VM that like thin elastic shells contains both stretching and bending energies.Here we show how such a coarse-graining procedure can be performed and that it explains the morphological instability.We further show that the threshold for this instability can be actively controlled by epithelia through apico-basal polarity.
Continuum model.We start with the Hamiltonian of a 3D VM with apical, basal and lateral faces, with surface tensions Γ i and areas A i for apical, basal and lateral faces, respectively (i = a, b, l).The factor 1/2 avoids membrane double counting.Assuming volume V being conserved, we non-dimensionalize energy with Γ l and length with V 1/3 .To derive a thin-plate theory from the 3D VM, we consider the nonlinear theory of moderately bent plates, where the total energy is given by stretching and bending energy contributions [22].For the in-plane stretching energy we assume the usual energy density with two-dimensional strain tensor ε and Lamé coefficients µ and λ.We determine the Lamé coefficients in a flat configuration in a mean-field fashion, following ideas from earlier work on 2D VMs [23][24][25].We consider a constant strain tensor for an individual cell, i.e. we assume strain to vary on a larger length scale than cell size.We assume a regular n-gonal lattice structure, which we will take to be a hexagonal lattice with n = 6.We then determine the energy density for a given strain.The equilibrium height can be obtained via minimization of Eq. ( 1) at constant volume.For the lateral faces we employ an angular averaging method (described in detail in the supplement [26]).In addition we consider non-affine deformations, as previously described for 2D models [25,27].Non-affine relaxations correspond to additional relative deformations of the two sublattices that make up the hexagonal lattice and can be included by allowing for an additional degree of freedom in the mid-plane shape, which allows for force balance at triple membrane junctions via angular relaxation.Note that the π/3 rotational symmetry implies that our 2D continuum model has only two elastic constants, exactly like an isotropic 2D material.
Considering the deformed areas and Taylor-expanding in the principal strains, we find 2µ = λ = Γ a + Γ b , or, equivalently, for the two-dimensional Young's modulus and Poisson ratio, respectively.The Young's modulus does not depend on the lateral tension Γ l , because both Γ a/b and Y are in units of Γ l .The reason is that changes in Γ l will affect height and edge length in such a manner that the energy density stays the same.A 2D Poisson ratio of 1/2 means that the sheet is compressible (incompressible materials in 2D have ν = 1), because it can exchange material between the in-plane and out-of-plane dimensions.In addition, we formulated the stretching energy in the fully nonlinear setting (see supplement [26]).The resulting energy density does not match standard hyperelastic models [28].Thus in the following we restrict ourselves to the first (cubic) correction to the linear theory as obtained by a Taylor expansion.This yields for the Young's modulus and bulk modulus in uniaxial and isotropic stretching, respectively.
For the bending energy density we assume the Helfrich form with bending rigidity κ, mean total curvature H = c + c ′ with the principal curvatures c and c ′ , spontaneous curvature c 0 , saddle splay modulus κ G and Gauss curvature K = cc ′ .Rozman et al. [12] have proposed a method to determine these quantities for the 3D VM using quadratic lattices.We generalized this to hexagonal lattices and adapted it such that we can formulate a theory for moderate bending.Consider a cell which is bent with principal curvatures c and c ′ and with unchanged center height.Volume conservation and curvatures determine the apical, basal and lateral face areas after deformation.Normalizing to the undeformed mid-plane area (for consistency and different from Ref. [12]) and Taylor-expanding with respect to ch and c ′ h yields the bending energy density (see supplement for full derivation [26]).Like for the stretching part, we also consider non-affine deformations, which leads to a correction factor k c .For n = 6 we find Bending with c ̸ = c ′ is accompanied by non-isotropic stretching in the apical and basal planes and angular relaxation, similarly as in non-isotropic stretching of the mid-plane.For c = c ′ such non-affine relaxations do not occur since lateral membrane angles do not change.In this case the energy is identical to the case of k c = 1, but otherwise k c is a numerical factor that we obtain from fitting to the simulation results for cylindrical surfaces.
The results from Eq. ( 5) differ in several essential ways from the known formulae for thin elastic plates.For a thin elastic plate the bending rigidity would scale like κ ∝ Y h 2 (with 2D Young's modulus Y = Y 3D h).With our result for Y and the reference height h ∝ (Γ a + Γ b ) 2/3 (compare supplement [26]), this would lead to κ ∝ (Γ a + Γ b ) 7/3 .As seen in Eq. ( 5), the bending rigidity depends more weakly on the tensions, because we do not have to consider the area changes along the entire height.In fact, compression on one of the sides will lower the energy instead of increasing it, as it is the case in 3D elastic plates.The dependence is much stronger, however, for the saddle-splay modulus, as we cannot compensate for area changes of the polygonal faces via shape changes (e.g. from rectangles to trapezoids) when both curvatures do not vanish.We have the same leading order scaling as we would have for elastic plates because here the apical and basal area changes for a given curvature also enter quadratically in the cell height.The spontaneous curvature depends on the apico-basal tension asymmetry, as this introduces a preferred curvature to minimize the energies.For Γ a = Γ b it vanishes as expected.
Like for a thin plate, the full deformation energy can now be calculated as For moderately bent plates there is an additional coupling between the two terms.A mid-plane deflection f (x, y) will contribute to the strain tensor as Stretching and bending of a flat sheet.To test the continuum theory by computer simulations, we have implemented the 3D VM, Eq. (1), as module in the software suite Chaste [29], similarly to Ref. [11] (details in the supplement [26]).For stretching we implemented a finite- size rectangular plate of hexagonal cells.The monolayer consists of n x and n y cells in the x and y-directions, respectively (Fig. 2(a)).For now we assume Γ a = Γ b = Γ, i.e. we do not consider a spontaneous curvature c 0 from apico-basal polarity.
First we considered isotropic stretching with edge stresses σ xx = σ yy and measured the effective bulk modulus as B = λ + µ = ε xx /2σ xx .Then we considered uniaxial stretching with edge stresses σ yy = 0 and measured Young's modulus as Y = σ xx /ε xx and Poisson ratio as ν = −ε yy /ε xx .Fig. 2(b-c) demonstrates excellent agreement between the simulations results (symbols) and the nonlinear continuum theory (dashed lines).Moreover the elastic moduli B and Y from Eq. (3) (solid lines) correspond exactly to the limiting cases of vanishing strain.The Poisson ratio ν is close to 1/2 as predicted by Eq. (3).In the following, we will mainly discuss the case of linear elasticity, but will come back to our nonlinear results when needed.
Next we simulated bending of cylindrical and spherical surfaces.Fig. 3(a) shows the resulting bending rigidity κ as function of curvature c = 1/R of a cylinder segment with mid-plane radius R. For Γ = 1, we determine the non-affinity correction factor by a fit as k c ≈ 1.26 and show that it is caused by non-affine apical and basal deformations (compare supplement [26]).Our simulation results (symbols) agree well with the prediction from continuum theory, Eq. ( 5) (solid lines).Fig. 3(b) shows the results for the saddle splay modulus κ G .For this stretching contributions, arising from the non-developable spherical deformation, were determined via finite element simulations, implemented with FEniCS [26,30].Again we find good agreement with the theoretical pediction from Eq. (5) (solid lines) for small curvature.The deviations at larger curvature are related to finite-size effects, including overestimation of the stretching energy for differently sized plates.
Buckling of a compressed sheet.Our continuum theory effectively describes the epithelial monolayer as a moderately bent plate.A classical application of such a theory is plate buckling upon in-plane compression, which has also been demonstrated experimentally for epithelial monolayers [31], depending on both stretching and bending.
Fig. 4(a) shows the amplitudes of a simply-supported rectangular plate, compressed along the edges parallel to the y-direction with compressive stress σ xx .We assumed straight (but movable) edges, as if the plate was situated in a movable rigid frame, and found a bifurcation toward a bent state with one half-wave along both axes.The critical stress in the 3D VM is slightly smaller than the continuum expectation, which can be explained by the nonlinearities, which we have neglected in the meanfield model, and by finite-size effects of the plate.The post-buckling amplitude for straight edges can be approximated within thin-plate theory [32] and is shown with solid lines.For this the leading-order Fourier modes are considered and the energy is minimized for these modes (compare supplement [26]).We do see good agreement between this approximation and the 3D VM simulation results.In Fig. 4(b) the mid-plane deflection is shown.For large deflections we see a flattening of the profile with stronger deviations from a sinusoidal leading-order approximation, consistent with real thin elastic plates [33].
Topological defects and icosahedral instability.The elastic framework derived above for epithelial monolayers suggests to study the effect of topological defects on spherical shells in the same manner as usually done for 2D elastic crystals.In our context, the disclinations are the pentagonal cells in the hexagonal monolayer.For such a defect with disclinicity s = 2π/6, the in-plane azimuthal stretching energy of a disc scales quadratically with the radius.This deformation will become unstable toward a conical bending deformation with a logarithmic scaling in the energies for large radii [10,17].Thus elastic shells of sufficient size, like large virus capsids or buckyballs, undergo a shape instability, in which each of the 12 pentagonal defects becomes the corner of an icosahedron.
As already shown in Fig. 1, our computer simulations of the 3D VM demonstrate exactly this scenario.To provide more details, Fig. 5(a) shows the asphericity α = ⟨(R − ⟨R⟩) 2 /⟨R 2 ⟩ of the cell centers as a function of the rescaled quadratic radius.It is well known that the transition depends on the ratio of bending rigidity and Young's modulus, which sets the relevant length scale κ/Y [17].In addition, we introduce a nonlinearity correction k ico .For an s = 2π/6 disclinicity the azimuthal strain is as large as 0.2 and we are in the strongly nonlinear regime, cf.Fig. 2(c).We find that a correction factor of k ico ≈ 1/2 is necessary to account for this.With this scaling all curves in Fig. 5(a) collapse onto the continuum limit (solid line taken from the literature [10]) except for small radii.This deviation can be understood to be a finite size effect as small radii correspond to few cells and large lattice constants compared to the radius.Notice that we indeed find a conical deformation at the pentagonal tips of the icosahedron (Fig. 5(b)), where the inner membrane can even collapse for large Γ and thus large heights (Fig. 5(c)).Experiments suggest Y = 200 kPa µm and V 1/3 ≈ 10 µm [34], resulting in κ ≈ 2200 kPa µm 3 for the VM with Γ = 1.The buckling threshold is known to be k ico Y R 2 crit /κ ≈ 154 [10].The critical radius is thus R crit ≈ 60 µm, which is roughly the size at which intestinal organoids undergo budding [15,35].
In passive elastic shells, topological defects tend to form additional structures such as defect scars, which screen the effect of the single defects [36,37].For flat epithelial monolayers it has already been established that active processes modulate their elastic behavior [11], thus also affecting the role of defects.For spherical epithelia active, apico-basally polarized forces become essential for structure formation, as observed experimentally.For example, in cell extrusion cells are pushed outward through contraction [38,39], and in budding organoids luminal (apical) contraction in buds facilitates curvature generation [15].To study such processes in our context, we add polarity in the pentagonal defect cells and their nearest neighbors, by using finite values for ∆Γ = Γ a − Γ b .Fig. 5(d) shows that such concentration of curvature generation around the topological defects facilitates buckling at smaller radii, allowing for active control of the instability threshold.Thus the instability can already occur in a neighbourhood of a few hexagonal cells, with potential implications for organoid formation and cell extrusion in less structured epithelia.Indeed, such hexagonal regions have been observed experimentally for epithelia with and without curvature [40,41].
In summary, here we have shown with computer simulations and analytical calculations that with growing size, spherical epithelia should undergo the same morphological instability at topological defects that is known also for elastic shells such as virus capsids.Our theory applies as long as the system is sufficiently regular and does not become too heterogeneous (e.g. by cell differentiation) before the threshold is reached.Therefore we expect that experimentally it might be observed best for highly regular epithelia, such as the retinal pigmented epithelium.Indeed, our theory might explain the formation of drusen, which are spherical or conical out-of-plane deformations in the retina linked to makular degeneration [42][43][44].
In the future, it has to be seen how the elastic effects described here will be modulated by the dynamics of epithelial monolayers, both on the level of single cells [45] and on the tissue level [46,47].At any rate, however, our theory demonstrates that topological defects are not only important for single cell behaviour, but have a strong effect on the global properties of epithelial monolayers and thus could mediate long-ranged effects.
Looking at the mid-plane, we can parametrize the area as shown in Fig. S1 and as proposed in Ref. [1].The cell has lengths ℓ x and ℓ y and an internal hexagonal angle of 2π/3.We do not consider bending in the derivation of the elastic constants and thus the apical and basal areas are identical to the mid-plane areas Due to membrane tension the internal angle is assumed to be unchanged when the mid-plane is stretched and as such we may parametrize the lateral areas via the total perimeter of the cell The lateral area is given by A l = P tot h.For the reference state we may calculate the height and edge length with minimal energy by minimizing Eq. (S1) with respect to h, i.e.The mid-plane is subject to a deformation described by the diagonalized strain tensor For the deformed state we may now look at deformed quantities ℓ ′ x , ℓ ′ y , h ′ , which are deformed via the principal strains α and β, cf.Fig. S2(a), as such that We have used the volume conservation V = 1 here.The deformed energy for one cell can then be computed via with (undeformed) hexagonal edge length s.Dividing by the mid-plane area ζ geom s 2 , inserting the reference height and edge length, Eq. (S4), and Taylor expanding with respect to α and β to second order, we obtain the energy density which is the linear elasticity result considered in the main manuscript.
To compare to the energy without non-affine angular relaxation, cf.Fig. S2(b), we consider the deformed lateral areas after affine transformation.The apical and basal areas after deformation are identical to the case with angular relaxation: the deformed mid-plane area is For the lateral faces we employ an angular averaging method: we first consider the edge vector e = s(cos(θ), sin(θ)) with edge angle θ.Then the deformed lateral area is determined from the deformed edge length and volume conservation to read Thus, the total lateral area for one cell is given by and we find for the energy density The behavior for α ̸ = β is different with and without angular relaxation.For the case α = β, where isotropic stretching does not change the angles, we find that both energy densities match, as they should.Note that our assumption on the angular orientation of our lattice with respect to the principal strain directions does not influence the resulting elastic constants.The 6-fold symmetry in the stretching energies implies that the material can be described by exactly two elastic constants in two-dimensional linear elasticity and as such the stretching energy should be valid within the applicability range of linear elasticity for arbitrary directions.In the manuscript we have checked numerical results of the (effective) bulk modulus, Young's modulus and Poisson ratio based on the linear relationship between strain and stress.For effects that depend both on bending and stretching we are usually interested in the stretching energy.To show that our continuum model describes the numerical VM well, we have additionally compared the continuum and numerical stretching energy densities in Fig. S3.We show the energy densities (per unit mid-plane area) divided by the quadratic strain as this is constant in linear elasticity with for isotropic stretching and for uniaxial stretching.
We do find similar results as in the moduli in the main manuscript.The nonlienar cubic model with coefficients as derived in the following section describes the energy well.Linear elasticity describes the behavior for vanishing strain.We find this behavior for both the isotropic and uniaxial stretching cases, cf.Fig. S3(a) and (b), respectively.

II. NONLINEAR TREATMENT OF THE STRETCHING ENERGY
In order to describe the VM in the framework of linear elasticity theory we have considered the stretching energy only up to quadratic order in the principal strains.To understand how the (effective) elastic parameters change with increasing strain we now consider higher orders.Going back to Eq. (S7), we find for the full stretching energy density In nonlinear elasticity theory (NLET) we now consider the nonlinear strain tensor [2] which means that the principal strains for a deformation in which the edges are stretched by 1+α and 1+β, cf.Fig. S1, will be α + α 2 /2 and β + β 2 /2.In NLET we consider the new coordinate of a material point χ(x, y) and must consider the deformation gradient tensor The central quantity for the energy in NLET is the Cauchy-Green tensor which has eigenvalues λ 2 x and λ 2 y -the squared principal stretches -which are related to the principal strains as such that in our case We may now rewrite the stretching energy density, Eq. (S14), in these principal stretches to find Equivalently, we may look at the invariants of the Cauchy-Green tensor to find This energy density does not match the standard set of widely used hyperelastic materials such as Neo-Hookean or Ogden.
To simplify our discussion of the nonlinearities of the elastic coefficients we now consider the next-higher order in the Taylor approximation of the energy density, Eq. (S7), For isotropic stretching we have α = β = ε xx and thus for the energy density in both linear elasticity and our cubic Taylor approximation i.e.
which is shown in the main manuscript.
For uniaxial stretching we have α = ε xx and β = −νε xx .Inserting this again into the energy densities of our Taylor approximation, and of linear elasticity yields The orthogonal strain ε yy = −νε xx and thus also the Poisson ratio minimize the stretching energy (density) for a given ε xx .We may differentiate the energy density with respect to ν to find this minimum and get We want to find the linear dependence of ν on the strain, i.e. we use the ansatz ν = 1/2 + νε xx , and find where we neglected orders quadratic and higher in the strain, i.e.
Inserting this into the energy density we also find which are the functions depicted in the main manuscript.

III. FULL DERIVATION OF MEAN-FIELD BENDING CONSTANTS
We consider the non-dimensional Hamiltonian for a single cell, Eq. (S1), and assume principal curvatures c and c ′ with corresponding radii R and R ′ respectively.The center height h is assumed to not change when the cell is bent, similiarly to Fig. S2(c), Assuming that a mid-plane edge is parallel to the axis with curvature c, the lateral face will form a trapezoid and we find for the apical and basal edge lengths, a and b, respectively, Note that s = (a + b)/2.We have assumed the basal edge length to be smaller for positive curvature here.As the apical and basal areas are related to edge lengths in both perpendicular directions we may write for the areas Noticing that for a cut through the cell at a height of z with respect to the mid-plane the areas can be written accordingly, we find for the bent cell volume with mid-plane area A mid = ζ geom s.For the derivation of the lateral area we consider the curvature orthogonal to the edge c and introduce the geometrical constant η geom , which relates the polygonal edge length to the in-radius (for hexagons η geom = √ 3/2)., which we will denote as t s for the mid-plane.On the apical and basal side we have respectively.The trapezoidal height of the lateral face is To determine the total lateral area we employ an angular averaging scheme: we consider the curvature in the direction with angle θ from the direction of principal curvature c and average the curvature over the n-fold symmetry.The total deformed lateral area is thus where the deformed mid-plane edge length s is determined via volume conservation, Eq. (S30), and the trapezoidal height is given by Eq. (S32).The apical and basal energy for a bent cell is The lateral bent energy for the hexagon is Note that these equations are exact for c = c ′ and for rectangular cells, but not necessarily for other polygonal shapes.We now divide by the undeformed mid-plane area, A (0) mid = 1/h, to obtain the energy density and perform a Taylor expansion for small ch and c ′ h up to second order.We find where we rewrote the result to compare it to the bending energy, Eq. ( 4) in the main manuscript.The first line then corresponds to a mid-plane surface tension energy of the reference state, the second to the mean-curvature contribution to the bending energy, and the last line to the Gauss curvature contribution.We identify κ = 9 8 for the bending rigidity, saddle splay modulus and spontaneous curvature, respectively.Similarly as for uniaxial stretching, we have found that additional corrections of the coefficients are necessary due to non-affine deformations.Non-affine apical and basal deformations indeed justify such a correction in the case c ̸ = c ′ .In Fig. S4 we show cylindrically bent cells in the VM with nodes projected onto the surface before and after energy minimization.These apical and basal deformations are similar to the deformations in uniaxial mid-plane stretching, cf.Fig. S2(b).We thus introduce correction factors k c , k G and k 0 for the bending rigidity, saddle splay modulus and spontaneous curvature, respectively, such that Such non-affine deformations will only be present for cases with c ̸ = c ′ and should not occur for c = c ′ , so the total energy density should be the same with and without corrections for a spherical deformation with c = c ′ .
We fix k c based on numerical results for bent plates with one vanishing curvature.From the equal energies in the case of c = c ′ we can derive for the saddle splay modulus correction i.e.
To obtain the results in the main manuscript, we now insert the reference height (ignoring height changes due to stretching) h min , Eq. (S4), and the geometric constants for hexagons.
For the spontaneous curvature computing such a correction is not as simply done.Assuming the bending energies with and without corrections are identical for c = c ′ , leads to a quadratic equation for the k 0 , i.e.
which has the solution of k 0 = 1 for k c = 1 as it should.With the correction in the mean curvature, k c ̸ = 0, we obtain For k c = 1 we see that k 0+ is valid for c 0 < 2c and k 0− for c 0 > 2c, but for k c ̸ = 1 the continuity of k 0 is lost if the curvature matches the spontaneous curvature here.
As such it makes sense to assume that the curvature at which the minimum in the mean-curvature bending energy occurs is not shifted by the correction, which implies k 0 = 1.However, this will lead to energetic changes.These are a consequence of the assumption of small ch, ch ′ in the above derivation, which will not properly work for c 0 ̸ = 0. Future work with Γ a ̸ = Γ b could investigate possible expansions at H ≈ c 0 insted of H ≈ 0 to mitigate this issue.This is, however, beyond the scope of this work.

IV. NUMERICAL VERTEX MODEL IMPLEMENTATION Formulation of the vertex model
The three-dimensional vertex model (VM) is a collection of polyhedral cells C with corresponding faces F(c) for c ∈ C. The faces are a collection of nodes r i ∈ R 3 , see Fig. S5, and can be separated into three different classes: apical, basal and lateral faces.We assume the apical and basal network topologies to be identical, i.e. the lateral faces are always (possibly non-coplanar) rectangles with four nodes and two lateral edges connecting the corresponding apical connecting oppsoing apical and basal nodes, at the plate boundary are constrained to be parallel to the averaged normal vector of the apical and basal faces of the cells to which the edge belongs.This allows for an orthogonal lateral edge at the boundaries.
We used orthogonal boundaries for finite-size plates, i.e. lateral edges are fixed to be parallel to the averaged apicobasal normal vectors.Simply-supported boundary conditions are implemented by constraining the lateral edges' centers of the outer most nodes on the edges to the initial mid-plane surface (z = h/2).If the edges do not fulfill this condition, the corresponding nodes are moved accordingly at every integration step.This allows for a nonvanishing derivative of the deflection function at the boundaries but constrains the mid-plane position.To implement cylindrically or spherically bent plates all lateral edge centers are constrained to the surface of the corresponding deflection function.To implement straight boundaries of the rectangular plate the lateral edge centers are enforced to be in the corresponding plane defined by the boundary normal and the center of mass of all the outer-most lateral edges.

Icosahedral sphere generation
To generate spheres with icosahedral symmetry the center nodes of a corresponding icosahedron, following the Caspar-Klug construction, are projected onto a sphere.We then create the Delaunay triangulation on the sphere by taking the convex hull [6], which yields the face centers as triangulation points.The corresponding Voronoi vertices are then computed as the circumcenters of the Delaunay triangles on the sphere and serve as the nodes on the luminal side of the sphere.The outer side is created by duplicating the nodes, increasing the radius, and constructing the final mesh.All icosahedral spheres are initialized on the spherical surface and then the energy is minimized.

Spontaneous curvature patterning
For the spontaneous curvature patterning the apical and the basal surface tensions of the pentagonal defect cells and their hexagonal nearest neighbors were chosen as such that the sum of the tensions is unchanged, i.e.
where we used Γ = 1.1 in the manuscript.This will induce a spontaneous curvature without changing the bending rigidities or the elastic moduli in the cells.The apical side is assumed to be luminal, i.e. toward the inside of the (icosahedral) sphere, matching the polarity in experimental situations such as intestinal organoids growing in three-dimensional culture.

V. FINITE ELEMENT METHOD FOR SPHERICALLY BENT RECTANGULAR PLATE
Since the deflection function for a spherical deformation with radius R is known, i.e.
we can solve the thin plate theory, with strain as we would solve the in-plane problem of linear elasticity.We want to solve the following system of equations on Ω = with stress-free edges and u x (0, 0) = u y (0, 0) = 0. We have split the strain tensor into the component resulting from the in-plane deformation εij (u) = (∂ i u j + ∂ j u i )/2 and the component from the deflection function fij = (∂ i f ∂ j f )/2 for notational convenience.
To obtain a weak formulation for a solution with the finite element method we introduce a vector valued test function v and consider where we used partial integration and the boundary condition.This can now be rewritten into the variational formulation, i.e. find u in a suitable vector-valued finite element space such that for all v a(u, v) = L(v). (S53) The bilinear form for our problem reads The right hand side contains all the other terms from the deflection function This linear system was solved in FEniCS [7] for given Ω and R where for both u and v linear Lagrange elements were chosen.

Comparison of continuum and VM height distributions in bent monolayers
Due to the volume conservation mid-plane area changes are related to height and we may calculate the deformed height h ′ via h ′ = h/ det(1 + ε).

VI. POST-BUCKLING APPROXIMATION FOR COMPRESSED RECTANGULAR SIMPLY-SUPPORTED PLATE
For a simply-supported rectangular plate, compressed along the edges parallel to the y-direction with compressive stress σ xx , the critical buckling stress is given by [8] where the multiplicity m describes the number of half-waves along the x-direction and it is chosen such that the critical stress is minimal, e.g.m = 1 for quadratic plates.
To approximate the post-buckling behavior, we consider a rectangular plate with straight edges (i.e. it is in a movable but straight frame).We introduce the half-lengths Lx = L x /2 and Ly = L x /2 and choose the center of the plate as the zero in the coordinate system.The rectangular plate is compressed along the x-direction, i.e. it is loaded along the edges x = ± Lx .
Solving the post-buckling behavior exactly is very difficult and can only be achieved numerically [9].However, an approximation is possible based on the dominant mode in the Fourier series description of the deformations.The in-plane deformations u and the deflection function can be assumed to be [8] where e is given by the strain in x-direction from the compression, i.e. e = −ε xx = σ xx /Y with compressive stress σ xx .The constant a is determined from the condition that the integrated stress along the unloaded edge vanishes and is given by The deflection amplitude A and the in-plane deformation amplitudes C x and C y are unknowns and will be determined via energy minimization.
The stresses can now be determined in the framework of moderately bent plates, cf.Eq. (S50).The stretching energy can be computed via the energy density integral For the bending energy one can show that the Gauss curvature contribution vanishes for the given deflection function and one finds [8] The total energy is thus This is a fourth-order polynomial equation in A. The minimum of this energy function is then determined by the zeros of the derivative with respect to the different parameters.Due to the form with respect to A we then obtain the normal form of a pitchfork bifurcation.As this function is quite involved, we determine the minimum numerically by minimizing the energy with respect to the parameters A, C x and C y .
In order to compare the VM results with this post-buckling bifurcation we need to additionally consider that the in-plane deformations u are considered inside the z = 0 plane due to the assumption of moderate plate bending and not inside the real buckled mid-plane.For thin elastic plates this yields a discontinuous transition to a smaller (in-xy-plane) stiffness in the post-buckling regime, which is known to be approximately [9] Y post ≈ Y /2.
(S62)   S1.Fitting-results for pre-buckling and post-buckling slopes and post-buckling correction factor for Young's modulus kY .
To verify this post-buckling behavior in our VM simulations, we show the compressive stress as function of e = −ε xx for the data from Fig. 4 in the main manuscript in Fig. S7.We observe that indeed the Young's modulus, i.e. the slope, of the VM can be described by two linear functions.We used linear regression on the first three points and the last four points for each Γ in Fig. S7 to determine the slopes before and after buckling and calculated the numerical postbuckling stiffness correction k Y as the ratio of the buckled and non-buckled slopes, cf.Table S1.Using Y post = k Y Y we then determined the bifurcation diagram depicted in Fig. 4 in the main manuscript by calculating e from the applied stress and minimizing the energy, Eq. (S61), to determine A(σ xx ).
To show the agreement of the proposed deformation fields with the main characteristic features of the compressed VM plate, we can consider the deformed height of the cells.It is related to the (hydrodynamic) area-changing strain via h ′ = h/ det(1 + ε).We find det ( S63) for the parameters that were found in the minimization, cf.panel (b).The highest order mode in the deformation, i.e. the compression (and larger height) in the plate corners, is captured by the theoretical approximated continuum deformation.However, we do see deviations in the center of the edges and in the center of the plate, which stem from the neglection of higher order modes in the Fourier representation of the deformation and deflection fields.Notice, however, that the order of the height changes is comparable and as such the stretching energy should be approximated reasonably well.
FIG. 1. Spherical epithelia described by a 3D vertex model experience an icosahedral instability that is well known for spherical elastic shells like virus capsids, polymerized vesicles or buckyballs.(a) A small epithelial shell stays spherical.(b) A large epithelial shell becomes faceted with icosahedral symmetry.Cells are arbitrarily colored in grey and pentagonal cells, which are topological defects, are shown in red.

FIG. 4 .
FIG. 4. Buckling instability in a simply supported 3D VM sheet of size (27, 31) under uniaxial compression.(a) Buckling bifurcations in amplitude with the total compressive stress σxx as control parameter for different tensions Γ.The solid lines are the continuum mechanics results for plate buckling, the dashed line is the unstable unbuckled state, the dotted line is the critical buckling stress, and the symbols are 3D VM simulations.(b) 3D VM simulation of deflected monolayer for Γ = 0.7 and σxx = 0.03, indicated in (a) as diamond marker.Below the 3D image a cross-sectional view is shown.Color indicates the mid-plane deflection of cell.

FIG. 5 .
FIG. 5. Icosahedral instability of 3D VM shells.(a) Asphericity α = ⟨(R−⟨R⟩) 2 /⟨R 2 ⟩ of the shells for different apical/basal surface tensions Γ as function of the average cell-radius.The radii are scaled by the ratio of the Young's modulus Y and bending rigidity κ with a nonlinearity correction kico.The solid curve is the continuum prediction [10].(b) A cut-out of the icosahedral tip for Caspar-Klug indices (23, 0) [16] shows a conical deformation.(c) For large tensions the cell height increases and for large radii the inner membrane collapses at the defect.(d) Apico-basal tension asymmetry ∆Γ in the defects and their nearest neighbors lowers the buckling threshold.Cells with ∆Γ ̸ = 0 are shaded in dark (size (5, 0), Γ = 1.1).
Supplemental material for: Morphological instability at topological defects in a three-dimensional vertex model for spherical epithelia Oliver M. Drozdowski and Ulrich S. Schwarz * Institute for Theoretical Physics and BioQuant, Heidelberg University, 69120 Heidelberg, Germany I. FULL DERIVATION OF THE STRETCHING ENERGY DENSITY We consider the dimensionless vertex model (VM) Hamiltonian with Γ l = 1 = V :

FIG. S1 .
FIG. S1.Parametrization of the mid-plane of a cell undergoing stretching.For the angular relaxed state the initial undeformed angles (left) do not change upon stretching (right) due to membrane tension.
FIG. S2.Deformations on a single cell level determine the stretching and bending energies.(a) We consider a strain tensor ε which varies slowly compared to the size of an individual cell.The deformation imposed on the cells changes the energy compared to the reference configuration.(b) Stretching with an affine lattice deformation is not consistent with force balance at mid-plane junctions connecting the membranes.Subsequent angular relaxation yields the minimal energy shape.(c) A cell is bent with two principal radii of curvature R and R ′ , defining curvatures c and c ′ , respectively.For non-rectangular cells the resulting energy is approximated via angular averaging, see text.The lateral face of the cell with height h has trapezoidal height h and experiences curvature c.
FIG. S3.Comparison of the stretching energy densities of the numerical VM and the continuum model for isotropic stretching (a) and uniaxial stretching (b).Energy densities are devided by the quadratic strain to obtain the constant coefficient from linear elasticity.Symbols and colors are consistent to Fig.3in the main manuscript.Dashed and solid lines correspond to nonlinear cubic Taylor approximation and linear elasticity, respectively.

FIG
FIG. S4.Non-affine transformations for cylindrical bending.Projecting a cell onto a cylindrical deformation (left) does not give the equilibrium cell shape.Further relaxation leads to non-affine apical/basal deformations (right).
Fig. S6(a) shows the expected height, while the bent 3D VM monolayer cell heights are shown in (b).The continuum model predicts the height distribution well, but for the strongly bent monolayer shown in Fig. S6 it slightly underestimates the compression at the edge centers.
FIG. S6.Bending of a rectangular 3D VM monolayer.(a) Cell height for a spherically bent rectangular monolayer with Γ = 1.3, size (31, 35) and c = 1/20 from finite element energy minimization of the continuum theory.(b) The same for the 3D VM in the undeformed mid-plane.

3
FIG. S7.Stress-strain curves for the buckling rectangular VM monolayer from Fig.4(a) in the main manuscript.The pre-and post-buckling behavior is linear with different slopes (Young's moduli).The dashed lines are linear fitting results, summarized in TableS1and the horizontal dotted lines are the theoretical buckling thresholds.

2 .
FIG. S8.Cell heights for a simply supported, compressed and buckled plate with straight edges, which is the same plate as in Fig. 4(b) in the main manuscript.The plate has size (27, 31), tension Γ = 0.7 and compressive stress σxx = 0.03.(a) VM simulation result for cell heights, shown in undeformed mid-plane coordinates.. (b) Theoretically expected height from continuum approximation.

Figure
FigureS8depicts the heights of the VM monolayer shown in Fig.5(b) in the main manuscript, cf.panel (a), and the theoretical height that we expect based on Eq. (S63) for the parameters that were found in the minimization, cf.panel (b).The highest order mode in the deformation, i.e. the compression (and larger height) in the plate corners, is captured by the theoretical approximated continuum deformation.However, we do see deviations in the center of the edges and in the center of the plate, which stem from the neglection of higher order modes in the Fourier representation of the deformation and deflection fields.Notice, however, that the order of the height changes is comparable and as such the stretching energy should be approximated reasonably well.
This research was conducted within the Max Planck School Matter to Life supported by the German Federal Ministry of Education and Research (BMBF) in collaboration with the Max Planck Society.USS is a member of the clusters of excellence Structures (EXC 2181/1-390900948) and 3DMM2O (EXC 2082/1-390761711) funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy as well as of the Interdisciplinary Center for Scientific Computing (IWR) at Heidelberg.
The pre-and post-buckling behavior is linear with different slopes (Young's moduli).The dashed lines are linear fitting results, summarized in TableS1and the horizontal dotted lines are the theoretical buckling thresholds.