Active morphogenesis of patterned epithelial shells

Abstract
Editor's evaluation
Introduction
Model
Results
Discussion
Appendix 1
Appendix 2
Appendix 3
Appendix 4
Appendix 5
Appendix 6
Appendix 7
Appendix 8
Appendix 9
Data availability
References
Article and author information
Metrics

Abstract

Shape transformations of epithelial tissues in three dimensions, which are crucial for embryonic development or in vitro organoid growth, can result from active forces generated within the cytoskeleton of the epithelial cells. How the interplay of local differential tensions with tissue geometry and with external forces results in tissue-scale morphogenesis remains an open question. Here, we describe epithelial sheets as active viscoelastic surfaces and study their deformation under patterned internal tensions and bending moments. In addition to isotropic effects, we take into account nematic alignment in the plane of the tissue, which gives rise to shape-dependent, anisotropic active tensions and bending moments. We present phase diagrams of the mechanical equilibrium shapes of pre-patterned closed shells and explore their dynamical deformations. Our results show that a combination of nematic alignment and gradients in internal tensions and bending moments is sufficient to reproduce basic building blocks of epithelial morphogenesis, including fold formation, budding, neck formation, flattening, and tubulation.

Editor's evaluation

The article provides a physical description of shape transformations of epithelial tissues in three dimensions, subject to active forces generated within the cytoskeleton of the epithelial cells. The work is motivated by organoids and more generally by morphogenesis during development. Therefore, this study is useful not only for developmental biology but also for a general understanding of cellular properties, including membrane mechanics and cell shapes.

https://doi.org/10.7554/eLife.75878.sa0

Introduction

Morphogenesis of embryos and the establishment of body shape rely on the three-dimensional deformation of epithelial sheets which undergo repeated events of expansion, contraction, convergence-extension, invagination, evagination, tubulation, and branching (Gilbert and Barresi, 2020). Tissue folding, for instance, is involved at different steps of embryogenesis (Kominami and Takata, 2004; Sui et al., 2018), organ (Sumigray et al., 2018), or entire organism development (Livshits et al., 2017; Braun and Keren, 2018). Recently, the growth of in vitro organoids, organ-like structures derived from stem cells capable of self-renewal and self-organisation, has revealed the intrinsic ability of biological systems to self-organise into complex structures from simple building blocks (Huch et al., 2017; Kamm et al., 2018; Rossi et al., 2018). Early steps in organoid self-organisation often start through the formation of a hollow, fluid-filled unpatterned sphere, undergoing spontaneous symmetry breaking (Ishihara and Tanaka, 2018) for example, in neural tube (Meinhardt et al., 2014) or intestinal (Serra et al., 2019; Yang et al., 2021) organoids. How this repertoire of shape changes and complex organisation emerges physically is a fundamental question.

Continuum theories of active materials, treating the epithelium as an active liquid crystal, have proven highly successful to achieve an understanding of the mechanics and flows of cellular collective motion. Epithelia cultured in vitro exhibit patterns of orientational order and spontaneous flows which are consistent with predictions from hydrodynamic theories of active matter (Duclos et al., 2017; Duclos et al., 2018; Blanch-Mercader et al., 2021a). Constitutive equations involving a shear decomposition of tissue area and anisotropic elongation into cell shape changes, cell division, and cellular topological transitions can reproduce basic features of the developing Drosophila pupal wing (Etournay et al., 2015; Popović et al., 2017). Recently, several studies established a link between topological defects in tissue order, provided by cell elongation or internal anisotropic cellular structure, and morphogenetic events (Kawaguchi et al., 2017; Saw et al., 2017; Mueller et al., 2019; Maroudas-Sacks et al., 2021).

Here, we propose a description of three-dimensional deformations of a patterned epithelial spheroid, considered as a shell of active liquid crystal. We consider an active elastic shell theory which takes into account in-plane tensions and internal bending moments (Lomholt, 2006; Maitra et al., 2014; Sahu et al., 2017; Salbreux and Jülicher, 2017). Internal bending moments arise from an inhomogeneous distribution of stress across the tissue. Such inhomogeneities can arise from, for example, changes in cytoskeletal organisation along the epithelium apico-basal axis, or from apposed epithelial tissues with different mechanical properties (Braun and Keren, 2018; Maroudas-Sacks et al., 2021). Apico-basal gradients of contractility, for instance, play a key role in morphogenetic processes (Martin and Goldstein, 2014; Sui et al., 2018) and are effectively taken into account here by active bending moments.

We consider an initially spherically symmetric tissue subjected to spatially modulated internal forces. Our rationale is to consider a situation where chemical and mechanical processes are uncoupled, such that cell–cell communication mechanisms ensure symmetry-breaking of the sphere, which is then converted into a pattern of mechanical forces (Ishihara and Tanaka, 2018). We consider a particularly simple pattern where the spherical tissue is decomposed into two regions, subjected to different active forces, and explore shape changes that result from this pattern (Figure 1a). We compare the situation where internal tensions and bending moments are isotropic to a situation where a nematic field, provided by cellular anisotropic structures, orients the internal tensions and bending moments.

Figure 1 with 8 supplements see all

Download asset Open asset

A two-dimensional surface with nematic order represents an epithelial sheet undergoing active deformations.

(a) Schematic of an epithelial tissue with a cellular state pattern. (b) Parametrisation of the axially symmetric shell and its deformation with the flow $v$ , and components of the tension and torque tensors. We note that $m^{ϕ s} = {\bar{m}}^{ϕ ϕ} x$ and $m^{s ϕ} = - {\bar{m}}^{s s} / x$ . (c) Stresses integrated across the thickness of the sheet result in tensions $t_{i j}$ and bending moments $m_{i j}$ acting on the midsurface. Anisotropic and possibly different tensions (dark-blue arrow crosses) on the apical and basal sides of the epithelium result in anisotropies in $t_{i j}$ and ${\bar{m}}_{i j}$ , which can be captured by a nematic order parameter $Q_{i j}$ (e.g. blue rods on the top surface).

Model

Viscoelastic nematic active surface model for epithelial mechanics

We first discuss our mechanical description of the deforming tissue. We represent an epithelium as an active surface flowing with velocity $v$ (Salbreux and Jülicher, 2017). The surface is taken to be elastic with respect to area changes, and fluid with respect to pure shear in the plane of the surface. Indeed, cellular rearrangements can fluidify in-plane epithelial flows by allowing cell elongation and cellular elastic stresses to relax on long time scales (Popović et al., 2017). Here, we consider such long enough time scales of hours to days which are relevant to organoid and developmental morphogenesis (Gilbert and Barresi, 2020). We also assume here that cell division and apoptosis or delamination are not occurring, such that elastic isotropic stresses do not relax (Ranft et al., 2010). Implicitly, we assume that cells have a preferred cell area.

Epithelia typically have a non-negligible thickness compared to characteristic transverse dimensions, and the apical and basal surfaces have different structures and are regulated differently. Notably, the basal surface is in contact with the basal lamina, a layer of extracellular matrix (Khalilgharibi and Mao, 2021). Therefore, a purely two-dimensional representation of epithelial stresses would miss essential aspects of their mechanics. We therefore introduce here the tension tensor $t^{i j}$ , but also the bending moment tensor $m^{i j}$ which captures internal torques arising from differential stresses acting along the surface cross section (Figure 1b and c). We assume that the surface possesses a bending rigidity, captured by a bending modulus $κ$ . When the curvature deviates from a flat layer, a bending moment results from the surface curvature (Equation 6). In addition, active bending moments can arise in the surface (Salbreux and Jülicher, 2017), for instance, due to actomyosin-generated differential active stresses along the apicobasal axis (Messal et al., 2019; Fouchard et al., 2020).

Cellular force generating elements are not necessarily isotropic; for instance, because cytoskeletal structures exhibit a preferred orientation (Martin, 2020) or inhomogeneous distribution across cellular interfaces (Bertet et al., 2004), or because the epithelial cells themselves exhibit an elongation axis (Duclos et al., 2017). Therefore, we introduce a coarse-grained surface nematic order parameter $Q$ which quantifies the average level of orientational order in the tissue. We assume that the nematic order parameter is tangent to the active surface.

Force balance

On a curved surface we define the rotated bending moment tensor ${\bar{m}}^{i j} = - m^{i k} ϵ_{k}^{j}$ , which we adopt for convenience. The local force balance projected on the tangential and normal directions reads (Salbreux and Jülicher, 2017)

\nabla_{i} t^{i j} + C_{i}^{j} t_{n}^{i} = - f^{e x t, j}

\nabla_{i} t_{n}^{i} - C_{i j} t^{i j} = - f_{n}^{e x t} - P,

where notations of differential geometry are introduced in Appendix 1; briefly $C_{i j}$ is the curvature tensor, $g_{i j}$ denotes the metric tensor, and $ϵ_{i j}$ the antisymmetric Levi-Civita tensor, $n$ the vector normal to the surface, $t^{i j}$ is the tangential contribution of the tension tensor and $t_{n}^{i}$ its normal contribution, and $\nabla_{i}$ denotes the covariant derivative on the surface. The tangential and normal torque balance provide the transverse tension and antisymmetric part of the tangent tension tensor:

t_{n}^{i} = \nabla_{k} {\bar{m}}^{k i},

ϵ_{i j} t^{i j} = C_{i j} m^{i j} .

We assume an external force density $f^{e x t} = f_{n}^{e x t} n + f^{e x t, j} e_{j}$ acting on the surface in addition to a difference of hydrostatic (uniform) pressure $P = P_{i n} - P_{o u t}$ , but no external torques (Figure 1). Here, we consider situations at low Reynolds number, where inertial forces may be neglected, and where additional external forces are negligible, such that the surface as a whole is force-free, $\oint_{S} dS f^{e x t} = 0$ . Dissipative couplings to the external fluid are ignored here as the characteristic viscosity of a biological tissue ( $\sim 10^{5}$ Pa s; Marmottant et al., 2009; Guevorkian et al., 2010) is several orders of magnitude larger than that of water ( $10^{- 3}$ Pa s).

Constitutive equations

In line with our hypothesis describing the material properties of an epithelium, we use the following constitutive equations:

t_{s}^{i j} = (2 K u + ζ + (η_{b} - η) v_{k}^{k}) g^{i j} + 2 η v^{i j} + ζ_{n} Q^{i j},

{\bar{m}}^{i j} = (2 κ C_{k}^{k} + ζ_{c} + η_{c b} \frac{D}{D t} C_{k}^{k}) g^{i j} + ζ_{c n} Q^{i j} .

where $t_{s}^{i j}$ is the symmetric part of the tension tensor and, on a curved surface, the strain rate tensor $v^{i j}$ and the corotational time derivative of the curvature tensor $\frac{D}{Dt} C^{i j}$ are given by (Salbreux and Jülicher, 2017)

v^{i j} = \frac{1}{2} (\nabla^{i} v^{j} + \nabla^{j} v^{i}) + C^{i j} v_{n},

\frac{D}{D t} C^{i j} = - \nabla^{i} (\partial^{j} v_{n}) - v_{n} C^{i}_{k} C^{k j} + v_{k} \nabla^{k} C^{i j} + ω_{n} (ϵ^{i k} C_{k}^{j} + ϵ^{j k} C_{k}^{i}),

with $ω_{n} = \frac{1}{2} ϵ^{i j} \nabla_{i} v_{j}$ the normal component of the vorticity. $u$ is the area strain, measuring local changes of area relative to a reference value; a precise definition is introduced in Equation 14. $Q^{i j}$ is a traceless, symmetric tensor characterising nematic orientational order on the surface.

We now discuss these constitutive equations. The surface elastic response is determined by the area elastic modulus $K$ and the bending modulus $κ$ . The dynamical deformations of the surface are characterised by the two-dimensional shear and bulk viscosities $η$ and $η_{b}$ and the bulk bending viscosity $η_{c b}$ . While the shear and bulk viscosities penalise in-plane isotropic and anisotropic deformation rates, the bending viscosity penalises the rate of change of total surface curvature $C_{k}^{k}$ . The bending viscosity dampens normal deformations and prevents bending modes, which would otherwise have no dissipative cost and could result in numerical instabilities.

The remaining contributions to Equations 5; 6 proportional to $ζ$ , $ζ_{n}$ , $ζ_{c}$ , $ζ_{c n}$ correspond to active tensions and bending moments. $ζ$ is an isotropic active surface tension, $ζ_{n}$ is the in-plane nematic active stress, with $ζ_{n} > 0$ usually referred to as the ‘contractile’ active stress and $ζ_{n} < 0$ as the ‘extensile’ active stress (Marchetti et al., 2013). $ζ_{c}$ is the isotropic bending moment, which locally favours a spontaneous curvature $C_{k}^{k} = - ζ_{c} / (2 κ)$ . If the active surface corresponds simply to two parallel layers under surface tension $γ_{a}$ , $γ_{b}$ (such as an epithelium with apical surface tension $γ_{a}$ and basal surface tension $γ_{b}$ ), and separated by a distance $h$ , an active isotropic bending moment $ζ_{c} \sim h (γ_{a} - γ_{b}) / 2$ emerges in the surface to lowest order in the curvature tensor. The term in $ζ_{c n}$ corresponds to an anisotropic active bending moment. In the bilayer picture, where the active surface corresponds to two layers $a$ and $b$ , it could generally arise from differences between the two layers in the level of order $Q_{i j}^{a}$ and $Q_{i j}^{b}$ or in the level of nematic active stress $ζ_{n}^{a}$ and $ζ_{n}^{b}$ . For example, such differences could stem from two contractile (respectively extensile) layers with perpendicular nematic orientations $+ Q_{i j}$ and $- Q_{i j}$ (Figure 1c), or from two layers with parallel nematic order, but one subjected to contractile active stresses and the other to extensile active stresses.

In the absence of external forces, deformations of the epithelial shell are driven by distributions of active tensions and bending moments, which are prescribed on it through the isotropic profiles $ζ (s)$ and $ζ_{c} (s)$ , the anisotropic components proportional to $ζ_{n} (s)$ and $ζ_{c n} (s)$ , and the shape-dependent nematic order parameter.

We note that Equations 5 and 6 can be seen as generic constitutive equations for a nematic active surface with broken up-down symmetry but no broken chiral or planar-chiral symmetry, arising from an expansion in the curvature tensor and in the nematic order parameter $Q_{i j}$ of the tensor $t_{s}^{i j}$ and ${\bar{m}}_{i j}$ (Salbreux and Jülicher, 2017; Salbreux et al., 2022). For simplicity some allowed additional couplings entering the generic constitutive equations have not been taken into account here, notably active contributions to the tension tensor (Equation 5) and bending moment tensor (Equation 6) proportional to the curvature tensor $C_{i j}$ . Salbreux et al., 2022 provide a more general list of possible couplings for active fluid nematic surfaces.

Nematic order parameter

For simplicity here we assume that the nematic order parameter minimises an effective free energy, thus ignoring potential active effects on the ordering (Salbreux et al., 2022). We consider the following effective free energy of the nematic on a curved surface (De Gennes and Prost, 1995; Jiang et al., 2007; Kralj et al., 2011; Pearce et al., 2019):

F = \int d S (\frac{k}{2} (\nabla_{i} Q^{j k}) (\nabla^{i} Q_{j k}) - \frac{a}{4} Q_{i j} Q^{i j} + \frac{a}{16} {(Q_{i j} Q^{i j})}^{2}),

with the Frank elastic constant $k$ , which is assumed to be equal for all distortions. The Landau–de Gennes contribution is chosen such that for $k = 0$ the aligned state with $Q_{i j} Q^{i j} = 2$ is a minimiser for $a > 0$ . Additional coupling terms between the nematic and curvature tensor are not considered here for simplicity (Napoli and Vergori, 2012).

Deformations of a polarised active sphere

We now turn to describe axisymmetric deformations of a closed nematic active surface.

Geometric setup

The epithelium is represented by a thin spherical shell undergoing axisymmetric deformations (Figure 1b). Its two-dimensional midsurface $X (ϕ, s) \in ℝ^{3}$ is parametrised by the arc length coordinate $s \in [0, L]$ and the angle of rotation $ϕ \in [0, 2 π]$ as

X (ϕ, s) = (x (s) \cos ϕ, x (s) \sin ϕ, z (s)) .

The local tangent basis is given by ${e_{ϕ}, e_{s}}$ , and $n$ is the outward-pointing surface normal. The geometry of axisymmetric surfaces is described further in Appendix 1. We require that the metric component $g_{s s} = 1$ , which implies relations between the tangent angle $ψ (s) \in [0, π]$ and the shape functions $x (s)$ and $z (s)$

\partial_{s} x = \cos ψ,

\partial_{s} z = \sin ψ,

which, together with the meridional principal curvature

C_{s}^{s} = \partial_{s} ψ,

are sufficient to reconstruct the surface shape from the curvature $C_{s}^{s}$ . In this axisymmetric setup, the velocity field reads $v = v^{s} e_{s} + v_{n} n$ , with $v^{s}$ the tangential and $v_{n}$ the normal velocities.

The undeformed initial surface is a sphere $S_{0}$ with radius $R_{0}$ _, and all quantities defined on it are denoted with a subscript ‘0’. We define the area strain on a point of the surface as

u = \frac{d S - d S_{0}}{d S_{0}},

where $d S$ is the surface area element at the point considered on the surface, and $d S_{0}$ is the surface area element of the same material point on the sphere. With this definition, $u = 0$ on the initial sphere. We denote $s_{0} (s)$ the arc length position on the undeformed sphere $S_{0}$ of a material point at arc length position $s$ on the deformed sphere. One then has $u = f_{ϕ} f_{s} - 1$ with $f_{s} = \frac{d s}{d s_{0}}$ the meridional stretch and $f_{ϕ} = \frac{x}{x_{0}}$ the circumferential stretch. Integrating $f_{s}^{- 1} = f_{ϕ} / (u + 1)$ yields the arc length reparametrisation $s_{0} (s)$ between the initial and the deformed surface. The Lagrangian time derivative of the area strain (Equation 14) is related to the flow through

\frac{D}{D t} u = (1 + u) v_{k}^{k} .

Nematic order

Here, with axial symmetry, the nematic tensor $Q_{i j}$ has the non-zero component $q = Q_{ϕ}^{ϕ} = - Q_{s}^{s}$ . On the closed shell, the nematic director (Appendix 3), which represents the alignment, will have two +1 topological defects at the poles (Figure 3a) as a consequence of the Poincaré–Hopf theorem (Hopf, 1927). The order parameter $q$ vanishes there, creating defect cores of size $l_{c} = \sqrt{k / a}$ , which is the characteristic nematic length. In this geometry the Euler–Lagrange equation resulting from the free energy (Equation 9) is

\partial_{s}^{2} q = \frac{1}{2 l_{c}^{2}} q (q^{2} - 1) + \frac{\cos ψ}{x} (4 \frac{\cos ψ}{x} q - \partial_{s} q) .

An example solution of Equation 16 on the sphere is shown in Figure 3b. From the two possible states with $q = \pm 1$ in the bulk, respectively, we choose $q = 1$ for reference. This corresponds to circumferential alignment of the nematic order (Figure 3a, right). The sign of the tensions and bending moments is then only controlled by the $ζ$ -prefactors. For example, a nematic tension with $ζ_{n} > 0$ corresponds to circumferential active contraction, resulting in an elongated shape. For nematic bending moments, if one chooses $Q_{i j}$ to represent the order parameter on the outer side of the shell, the sign convention is such that $ζ_{c n} > 0$ , $q > 0$ results in circumferential contraction on the outer side and contraction along the meridians on the inner side of the shell. We note that the shape is only influenced by the order parameter via the active tension $ζ_{n} Q^{i j}$ and the active moment $ζ_{c n} Q^{i j}$ , but is otherwise insensitive to the nematic elastic energy (Equation 9). Minimisation of the Frank free energy by deformations of passive nematic surfaces has been previously discussed (Jiang et al., 2007).

Active profiles

We consider initially spherical epithelial shells containing an active region that drives the deformation. For the steady-state analysis, this region is a circular patch of size $l_{a} \leq L_{0}$ (Figure 1b), such that the active terms are given on $S_{0}$ by step-like profiles, for example

ζ_{c} (s_{0}) = {\begin{cases} ζ_{c}^{0} + δ ζ_{c}, & if s_{0} \in [0, l_{a}] \\ ζ_{c}^{0}, & otherwise \end{cases}

and similarly for $ζ (s_{0})$ , $ζ_{n} (s_{0})$ , and $ζ_{c n} (s_{0})$ . The circular patch deforms with the material points, which reflects that the active properties are associated with a predefined group of cells. If not stated otherwise, the values outside the active region are $ζ^{0} = ζ_{c}^{0} = ζ_{n}^{0} = ζ_{c n}^{0} = 0$ . This passive part of the surface is governed by the constitutive equations 5 and 6, but with vanishing active terms.

In dynamical simulations, active tension and bending moment profiles are defined on the spherical surface at time $t = 0$ using sigmoid functions $f (x, μ, σ)$ of the form

f (x, μ, σ) = 1 - {(1 + e^{- \frac{x - μ}{σ}})}^{- 1},

for their space and time dependence. For instance, the active bending moment profile is defined on $S_{0}$ as

ζ_{c} (s_{0}, t = 0) = (1 - f (t = 0, μ_{t}, σ_{t})) (ζ_{c}^{0} + δ ζ_{c} f (s_{0}, l_{a}, σ_{s}))

as a smooth version of the step-profile Equation 17, and $ζ$ , $ζ_{n}$ , and $ζ_{c n}$ are defined analogously. The profile is then advected with the material points (Figure 1b), while its intensity increases through the time-dependent sigmoid (e.g. Figure 2d).

Figure 2 with 10 supplements see all

Download asset Open asset

Volume

We consider two possibilities for the volume enclosed by the epithelium. In one limit the tissue is assumed to be impermeable and the enclosed volume is treated as an incompressible fluid exerting hydrostatic pressure on the tissue. The volume is conserved when the shell deforms:

V = V_{0},

with the pressure $P$ acting as the Lagrange multiplier.

In the other limit the tissue is fully permeable. At steady state, in this limit the volume can change freely and no pressure acts on the tissue, $P = 0$ . In dynamical simulations, we introduce a volume viscosity $η_{V}$ such that the pressure is coupled to the volume change via

P = - η_{V} \partial_{t} V

where $η_{V}$ is a parameter chosen to be small enough that the internal pressure is small compared to other forces.

Stationary shapes

For given profiles of active tensions and bending moments, steady-state shapes are obtained as solutions of the mechanical equilibrium equations. Those are a system of non-linear ode’s containing the force and torque balances Equations 1–4, the geometric Equations 11–13, the constitutive relations Equations 5–8 and Equation 14 with vanishing velocities $v^{s} = v_{n} = 0$ , and, if applicable, the nematic equilibrium Equation 16.

Dynamical deformations

In the dynamical version of the model a given active profile generates a velocity $v (ϕ, s, t)$ , whose normal part deforms the surface (Figure 1b). The components ${v^{s}, v_{n}}$ of this instantaneous velocity are obtained by solving the force and torque balance Equations 1–4 (derived for the axisymmetric surface in Equations 63–65), together with the constitutive Equations 5–8, on the shape $X (ϕ, s, t)$ . Since $u (s, t)$ and $Q_{i j} (s, t)$ are also given, these constitute a linear system of ode’s. The shape is evolved in time in a Lagrangian approach, in which material points move according to the full-velocity vector $v$ ,

\partial_{t} X = v .

Surface quantities, such as the active profiles and the area strain, are advected accordingly. The nematic order parameter evolves in time quasi-statically, where we assume that it relaxes instantaneously to the solution of Equation 16 written on the deformed surface at time $t$ .

Dimensionless variables

The equations are made dimensionless (marked by tilde) by rescaling tensions by $κ / R_{0}^{2}$ , bending moment densities by $κ / R_{0}$ , lengths by $R_{0}$ , force densities by $κ / R_{0}^{3}$ , viscosities by the two-dimensional shear viscosity $η$ of the epithelium, times by the characteristic time scale $τ_{a} = η R_{0}^{2} / κ$ , and velocities by $R_{0} / τ_{a}$ . This leaves the dimensionless parameters $\tilde{K} = K R_{0}^{2} / κ$ , ${\tilde{l}}_{c} = l_{c} / R_{0}$ , ${\tilde{η}}_{b} = η_{b} / η$ , ${\tilde{η}}_{c b} = η_{c b} R_{0}^{2} / η$ and ${\tilde{η}}_{V} = η_{V} R_{0}^{4} / η$ to be fixed. We choose to set $\tilde{η} = {\tilde{η}}_{b} = 1$ , ${\tilde{η}}_{V} = 10^{- 4}$ for fast relaxation of the volume, and the nematic length scale is set to ${\tilde{l}}_{c} = 0.1$ . Working under the assumptions of linear shell theory for a homogeneous thin shell (Reddy, 2006), one can relate the elastic moduli to each other via the thickness $h$ of the cell layer, and express $\tilde{K} = 12 {(R_{0} / h)}^{2}$ . In simulations we use $\tilde{K} = 1000$ , corresponding to $h / R_{0} \approx 0.1$ , which covers a range of systems from gastrulating embryos (e.g. sea urchin Davidson et al., 1995) to organoids (Serra et al., 2019). Similarly, for the bulk bending viscosity we have ${\tilde{η}}_{c b} \sim {(h / R_{0})}^{2} = 10^{- 2}$ .

Numerical methods

For both the steady-state computation and the dynamics, the resulting sets of ode’s are integrated numerically with the boundary-value-problem solver bvp4c of MATLAB, which implements a fourth-order collocation method on an adaptive spatial grid (Kierzenka and Shampine, 2001). The equations are solved on the full interval $[0, L]$ , and geometrical singularities at the poles are handled using analytical limits at $s = 0, L$ (Appendix 6). Any integral constraint, such as volume conservation, is rewritten as a boundary value problem and added to the system of ode’s to be solved.

The dynamics simulations start with a sphere at time $\tilde{t} = 0$ . We study each of the four active effects separately. The corresponding active profile is switched on smoothly via a sigmoid function in time, such that it reaches its target intensity at $\tilde{t} \approx 0.02$ . The time integration according to Equation 22 is done with an explicit Euler method with adaptive step size via

X^{'} (ϕ, s, t + δ t) = X (ϕ, s, t) + δ t v (ϕ, s, t) .

In order to keep the force and torque balance equations in the form given by Equations 63–65, the updated surface is reparametrised as $X^{'} (ϕ, s^{'}, t + δ t)$ in a new arc length $s^{'} (s)$ which is calculated from the condition $g_{s^{'} s^{'}} = 1$ . The profiles and surface quantities are passed between time steps as spline interpolants.

To produce the diagrams of steady-state shapes, $l_{a}$ is fixed and the control parameter is the difference of the active profile value between the passive and the active regions of the shell, for example, for the profile given in Equation 17 it is $δ ζ_{c}$ . A solution branch is found by starting from the spherical solution at zero difference of active profile, and calculating a sequence of steady-state shapes, progressively increasing the magnitude of the difference in activity. Two different methods are used to construct the solution branch for a sequence of control parameter values. For small values, starting from zero, the solution branch is obtained by making small increments in the control parameter. For larger values we switch to an implicit stepping method, which we developed based on a parametric representation of the solution branch (see Appendix 6 section ‘Construction of solution branches’). This second method allows us to continue the solution branches into regions where the steady-state shapes become non-unique in the control parameter.

Details of the numerical methods can be found in Appendices 6 and 7 for the steady state and the dynamics simulations, respectively.

Results

Epithelia as active membranes: Isotropic active tensions

We first consider deformations of an epithelial shell due to patterns of isotropic active tensions and bending moments. A spatially varying isotropic tension represents a change in the preferred area of the epithelium due to either changes in sheet thickness or cell number (Popović et al., 2017). However, one can show that a step-profile of positive (contractile) tension $ζ > 0$ does not lead, at steady state, to a three-dimensional deformation of the shell away from a spherical shape, which is a consequence of the absence of shear elasticity in our model (Appendix 8). Instead, the epithelium remains spherical and regions with higher tension contract. This leads to a rescaling of the relative active region size $l_{a} / L_{0}$ and, if the volume is free to change, also to a decrease in shell radius (Appendix 8). If the tension becomes negative, a buckling of the surface may occur (Salbreux and Jülicher, 2017). Here, we focus on positive tensions; therefore, if only isotropic active effects are considered, active internal bending moments are required to drive deformations away from the spherical shape.

Epithelia as active shells: Isotropic active bending moments

We now turn to deformations induced by an increasing active bending moment in a spherical cap. In Figure 2a and e, we plot a phase diagram of steady-state shapes as a function of the increased active bending moment $δ ζ_{c}$ and the size of the active region $l_{a}$ . The steady-state deformed shapes are plotted with the active region shown in red and the ‘passive’ region, where $ζ_{c} = 0$ , shown in blue. We can contrast the situation where fluid is free to exchange across the surface and at steady state the difference of pressure across the surface vanishes, $P = 0$ (Figure 2a–d), to the case where the volume enclosed by the surface is constrained to a fixed value (Figure 2e–i).

An isotropic active bending moment (term in $ζ_{c}$ in Equation 6) induces a preferred curvature ${(C^{0})}_{k}^{k} = - \frac{ζ_{c}}{2 κ}$ , such that regions of a spherical shell with $ζ_{c} > 0$ can be expected to flatten or bend inwards. Specifically, a difference of $δ ζ_{c}$ applied at the boundary of the active cap induces a jump in meridional curvature $C_{s}^{s}$ and a local folding of the sheet. Due to the spherical topology, the shape of the whole shell is affected by this fold, as can be seen from the sequences of stationary shapes obtained by increasing $δ ζ_{c}$ for intermediate values of $l_{a} / L_{0}$ (Figure 2a). In particular, for the same value of $δ ζ_{c}$ the active region may bend inward or keep a positive curvature, depending on its size.

When $l_{a} / L_{0}$ is small or close to 1, the resulting shape is characterised by the formation of a bud which form either inwards ( $l_{a} ≪ L_{0}$ ) or outwards $(L_{0} - l_{a} ≪ L_{0})$ . In these cases, for sufficiently large values of $δ ζ_{c}$ the steady-state solution is lost through the formation of a constricting neck. In our simulations the constricting neck is numerically resolved up to values of $\sim 10^{- 3} R_{0}$ ; extrapolation indicates full constriction at a finite $δ ζ_{c}$ (Figure 2i). As the neck radius decreases the principal curvatures at the neck diverge as $C_{s}^{s}, C_{ϕ}^{ϕ} \to \pm \infty$ , such that $C_{k}^{k}$ remains finite (Figure 2i) and therefore the limiting, budded shape is a true steady-state solution. Such a transition is reminiscent of models of lipid membrane vesicles, which can be induced to form a budded shape consisting of two spheres connected by an infinitesimal region called the ideal neck (Seifert et al., 1991; Jülicher and Lipowsky, 1993; Fourcade et al., 1994; Jülicher and Lipowsky, 1996; Seifert, 1997). For lipid membranes the ideal neck condition gives the difference in spontaneous curvature between the two domains at which a vesicle will form two spheres, $1 / R_{1} + 1 / R_{2} = C_{0}$ with $R_{1}$ and $R_{2}$ the radius of the two spheres and $C_{0}$ the spontaneous curvature (Seifert, 1997). Here the choice of constitutive Equations 5 and 6 does not correspond to the Helfrich model, and we find alternative matching conditions for the two regions connected by the infinitesimal neck: we find that $t_{s}^{s}$ changes sign across the neck, while ${\bar{m}}_{s}^{s}$ is continuous. This result can be derived by a scaling analysis around the neck (Appendix 2). In the free volume case, these conditions are satisfied when the active and passive regions are separated by the neck, and have the shapes of spheres with vanishing strain ( $u = 0$ ) and radii $R_{a}$ , $R_{p}$ , related by the condition:

\begin{aligned} \frac{1}{R_{a}} - \frac{1}{R_{p}} = - \frac{δ ζ_{c}}{4 κ}, \\ - \frac{1}{R_{a}} - \frac{1}{R_{p}} = - \frac{δ ζ_{c}}{4 κ}, \end{aligned}

where the change of sign in the second line arises because the active region deforms inward and form a sphere with a negative mean curvature. The additional condition of vanishing strain $u = 0$ gives an additional relation for $R_{1}$ and $R_{2}$ as a function of $l_{a} / L_{0}$ . Combining these conditions determine a curve in the parameter space $δ ζ_{c} R_{0} / κ$ , $l_{a} / L_{0}$ , which matches with the numerically determined curve of neck constriction (Figure 2b). In the fixed volume case, the matching conditions do not result in such a simple shape solution; however, using the same condition as for the free volume case appears to still provide a good approximation of the constriction point for small ( $l_{a} ≪ L_{0}$ ) and close to $L_{0}$ $(L_{0} - l_{a} ≪ L_{0})$ values of $l_{a}$ (Figure 2f). We conclude that infinitesimal neck formation can arise outside of the Helfrich model and that the ideal neck condition which is satisfied there does not generally extend to other models of surface mechanics.

At sufficiently large increase in the active bending moment difference $δ ζ_{c}$ and for intermediate values of $l_{a} / L_{0}$ , a fold in the solution branch in the $(δ ζ_{c}, V)$ -plane appears (Figure 2c). For most values of $l_{a} / L_{0}$ , this fold is associated to the loss of a continuously attainable solution with increasing $δ ζ_{c}$ , and a shape transition (Figure 2b and c). We expect shapes obtained by following the continuous branch of shapes beyond the fold to be unstable (Appendix 9). The (potentially unstable) physical branch eventually stops either through a self-intersection of the sheet at the poles (Figure 2c, $l_{a} / L_{0} = 0.35$ ) or through the constriction of a small neck that develops near the boundary of the passive and active regions and separates the shell into two smaller, approximately spherical compartments (Figure 2c, $l_{a} / L_{0} = 0.31, 0.7$ ). Alternatively the solution branch continues in a sequence of loops and the active region elongates (Figure 2c, $l_{a} / L_{0} = 0.5$ ), forming an increasing number of bubble-like compartments.

Since we follow continuous trajectories of steady-state shapes in parameter space, we cannot directly obtain alternative steady-state solution branches after the shape transition. Therefore, we turn to dynamic simulations where we explicitly calculate flow fields, starting from the reference spherical shape, and evolve the surface shape (Figure 2d) with parameters chosen to be away from the transition in parameter space (Figure 2a). This also allows to resolve the sequence of shapes and velocity fields leading to a given steady-state deformed shape (Figure 2h, $l_{a} / L_{0} = 0.1, δ ζ_{c} = 40 κ / R_{0}$ ). For parameters beyond the shape transition, we find that a small neck can form, separating roughly the active and passive regions, whose radius decreases to 0 over time (Figure 2d). Alternatively, the surface ends up self-intersecting (Figure 2d, $l_{a} / L_{0} = 0.35, δ ζ_{c} = 12.5 κ / R_{0}$ ). We do not find therefore alternative solution branches beyond the shape instability. Since intersection of the surface with itself is described by different physical interactions than considered here, our framework does not answer what would happen beyond the self-intersection line. However, assuming that self-intersection results in fusion and rupture of the apposed two surfaces, active isotropic bending moment difference could in principle drive a change in tissue topology, from one sphere to two ( $l_{a} / L_{0} = 0.85, δ ζ_{c} = 15 κ / R_{0}$ ), or from a sphere to a torus via self-intersection ( $l_{a} / L_{0} = 0.35, δ ζ_{c} = 12.5 κ / R_{0}$ ).

When volume is conserved, deformations are broadly similar but tend to be more localised to the fold at the active boundary (Figure 2e–i). For intermediate values of $l_{a} / L_{0}$ , the shell deforms into locally folded shapes, which eventually self-intersect at large bending moment difference (Figure 2g, $l_{a} / L_{0} = 0.3, 0.7$ , Figure 2h).

Nematic active tensions

We now introduce the nematic order parameter $Q_{i j}$ and consider shape changes driven by contractile or extensile active stress in the active region (Figure 3). As expected, solving for the nematic order parameter profile on the undeformed sphere results in maximal order at the equator and two defects at the poles where the nematic order parameter vanishes, $q = 0$ (Figure 3). Two solutions with $q < 0$ and $q > 0$ can exist; in the following we take the convention that $Q_{ϕ}^{ϕ} = q > 0$ , $Q_{s}^{s} = - q < 0$ , corresponding to circumferential alignment of the order parameter, such that a contractile active stress $(ζ_{n} > 0)$ results in a positive circumferential tension, $t_{ϕ}^{ϕ} > 0$ . Due to invariance of the constitutive equation by exchange $Q_{i j} \to - Q_{i j}$ , $ζ_{n} \to - ζ_{n}$ , the same shape deformations occur when considering meridional alignment of the order parameter ( $q < 0$ ) and exchanging contractile $(ζ_{n} > 0)$ and extensile $(ζ_{n} < 0)$ active stresses.

Figure 3

Download asset Open asset

Nematic order on a sphere.

(a) Two possible configurations for the nematic order parameter $Q_{i j}$ on a sphere with a + 1 topological defect at each pole: meridional (left) or circumferential (right) alignment. The order parameter minimises an effective energy (Equation 9 with $l_{c} = 0.1 R_{0}$ ). (b) Order parameter $q (s) = Q_{ϕ}^{ϕ} (s)$ as a solution of the Euler–Lagrange Equation 16 on a sphere with $R_{0} = 1$ and $l_{c} = 0.1 R_{0}; q = 1$ at the equator and $q = 0$ at the locations of the defects (poles). For uniform $ζ_{n}$ , $ζ_{n} \nabla_{i} Q^{i s}$ is the active nematic contribution to the tangential force balance (Equation 63) and, close to the equator, results in the elongation of the surface along the axis of symmetry for $ζ_{n} > 0$ , and its contraction for $ζ_{n} < 0$ .

As before, we study the cases of vanishing pressure difference across the shell (Figure 4a–e) and constrained volume inside the shell (Figure 4f–i). With a nematic tension profile on the surface, a deformation away from the spherical shape occurs even for homogeneous active nematic tension, $l_{a} / L_{0} = 1$ (Figure 4a–e).

Figure 4 with 13 supplements see all

Download asset Open asset

Deformations of epithelial shells due to nematic tensions, with free (**a–e**) and conserved (**f–i**) volume.

(**a, e**) Shape diagrams. (**b, g**) Details of shape diagram illustrating the behaviour of solution branches. (d) Curvature at the south pole for extensile stress. (**c, e, h, i**) Dynamic simulations of shell shape changes, for parameter values indicated in the phase diagrams (**a, f**). Other parameters: $\tilde{K} = 10^{3}, {\tilde{η}}_{c b} = 10^{- 2}$ , ${\tilde{η}}_{V} = 10^{- 4}$ , ${\tilde{l}}_{c} = 0.1$ .

In the extensile case $ζ_{n} < 0$ (or in the contractile case $ζ_{n} > 0$ if $q < 0$ ), and no pressure difference across the shell, the surface progressively flattens into a flat, double-layered disc (Figure 4b, $l_{a} / L_{0} = 1$ , $ζ_{n} < 0$ ). There is no shape transition occurring; instead, we find that the shape converges to a limit shape as $| ζ_{n} | \to \infty$ (Appendix 4). The limit shape corresponds to two parallel flat discs of radius $R_{d}$ , separated by a distance $2 h$ , connected by a narrow curved region. An asymptotic analysis (Appendix 4) shows that the radius of the disc and the separating distance obey the scaling relations, in the limit $κ ≪ K l_{c}^{2}$ :

R_{d} \sim l_{c}, h \sim {(\frac{κ l_{c}}{K})}^{\frac{1}{3}}

The first relation shows that the limit shape has the size of the characteristic nematic length $l_{c}$ . Physically, for $l_{c} ≪ L_{0}$ , the nematic active tension results in a contraction of the shape, until the shape is sufficiently close to the defect core for the nematic order to ‘dissolve,’ thus limiting further increase in the active tension.

In the contractile case ( $ζ_{n} > 0$ ), the shape elongates until a shape transition is reached, characterised by a fold in the solution branch (Figure 4b, $l_{a} / L_{0} = 1$ , $ζ_{n} > 0$ ). Following the solution branch after the fold eventually gives rise to a sequence of presumably unstable shapes with the formation of a central constricting neck. Intrigued by this result, we performed dynamical simulations for contractile active tensions above the shape transition (Figure 4c and e; Figure 4—figure supplement 1). Dynamic simulations show separation of the shape into two or more compartments via dynamical neck constrictions, with the neck radius vanishing over time (Figure 4—figure supplement 1a). Within the neck, $q \to 0$ as a result of the diverging principal curvatures (as can be seen from the presence of a term ${(\frac{\cos (ψ)}{x} q)}^{2}$ term in the nematic free energy, Equation 102). In particular, for values close to the branch fold (Figure 4c) the dynamics is reminiscent of cell division; however, in contrast to existing models of cell division (Salbreux et al., 2009; Turlier et al., 2014), the constriction appearing here does not require a narrow peak of active stress around the equator to occur. At larger contractile stress (Figure 4e), a narrow, elongated tube forms around the equator. This tube thins out over time, and two symmetric necks emerge and constrict, suggesting that the shape would eventually separate into three topologically separated surfaces (Figure 4—figure supplement 1b).

For $0 < l_{a} / L_{0} < 1$ and extensile stress in the active region $δ ζ_{n} < 0$ , the active region tends to flatten more and more strongly as $| δ ζ_{n} |$ is increased, and the total curvature vanishes at the south pole ( $C_{k}^{k} \to 0$ , Figure 4d). For $0 < l_{a} / L_{0} < 1$ and contractile stress $δ ζ_{n} > 0$ , a fold in the solution branch appears at large value of $δ ζ_{n}$ (Figure 4b and d). Following the solution branch beyond the fold results in a complex trajectory in parameter space, corresponding to successive additions of new bubbles to a linear chain of bubbles within the active region. This bubble chain is observed both with free or constrained volume (Figure 4b and g). Here, we cannot conclude however whether these shapes are unstable. Instead, we consider the shape dynamics for $δ ζ_{n}$ values larger than the shape transition, here at fixed internal volume (Figure 4h and i). Here, a neck forms within the active region and its constriction leads to the separation of a smaller bubble. For small enough $l_{a}$ the smaller bubble appears nematic-free and spherical (Figure 4h, Figure 4—figure supplement 1b). This is consistent with restoration of isotropic state stability which can occur on a sphere whose size becomes smaller or comparable to $l_{c}$ (Appendix 3 section ‘Stability of the isotropic state on a sphere’).

Active nematic bending moments

We now turn to shape deformations resulting from active bending moments oriented along the nematic order $Q_{i j}$ . As for nematic tension, we adopt the convention of nematic alignment along the circumference, $Q_{ϕ}^{ϕ} = q > 0$ ; alignment along the meridians can be studied simply by changing the sign of the active coefficient $ζ_{c n}$ .

We first discuss the case where the nematic active bending moment is homogeneous ( $l_{a} / L_{0} = 1$ ), where there is no difference of pressure across the surface, and where $ζ_{c n} = δ ζ_{c n} < 0$ (Figure 5a–c and g). We find that the sphere deforms into a shape with a central cylindrical part (Figure 5a and b). The length of the cylindrical part increases with increasing value of $| ζ_{c n} |$ . To characterise this, we note that the corresponding steady-state shape solutions have vanishing tensions $t_{s}^{s} = 0$ and $t_{n}^{s} = 0$ everywhere (Figure 5—figure supplement 1) and the force balances Equations 63 and 64 are trivially satisfied. The torque balance Equation 65 reads

2 κ \partial_{s} C_{k}^{k} - ζ_{c n} \partial_{s} q = 2 ζ_{c n} \frac{\cos ψ}{x} q .

Combining Equations 26 and 48 one obtains that $L [C_{s}^{s} - C_{ϕ}^{ϕ} - ζ_{c n} q / (2 κ)] = 0$ , with the operator $ℒ = \partial_{s} + 2 \frac{\cos ψ}{x}$ . Solutions to $ℒ [f] = 0$ have the form $f = A / x^{2}$ with A a constant. The boundary condition that the function $f$ should be finite at the poles requires $A = 0$ , such that

C_{s}^{s} - C_{ϕ}^{ϕ} = \frac{q ζ_{c n}}{2 κ} .

As a result, if the shape has a cylindrical part, in which $C_{s}^{s} = 0$ and $q = 1$ , then the cylinder radius $R_{c}$ is given by

\frac{1}{R_{c}} = - \frac{ζ_{c n}}{2 κ},

and since such solutions are area-preserving, with $u = 0$ , the length of the cylindrical part scales as $L_{c} \sim 1 / R_{c}$ . These relations are in excellent agreement with simulation results for large enough $| ζ_{c n} |$ (Figure 5g).

Figure 5 with 1 supplement see all

Download asset Open asset

Deformations of epithelial shells due to nematic bending moments, with free (**a–c**) and conserved (**d, e**) volume.

(**a, d**) Shape diagrams. (**b, e**) Details of shape diagram illustrating the behaviour of solution branches. (**c, f**) Dynamic simulations of shell shape changes, for parameter values indicated in the phase diagrams (**a, d**). In both cases in (f) the dynamics results in self-intersection. (g) Comparison of curvature and length of the cylindrical tubes for $l_{a} / L_{0} = 1, 0.7, 0.3$ , $δ ζ_{c n} < 0$ with analytical predictions. The tube length is measured on the steady-state shape as the arc length of the deformed active region, $s_{t u b e} = s (s_{0} = l_{a})$ , and the tube curvature as $C_{ϕ}^{ϕ} (s_{t u b e} / 2)$ . Other parameters: $\tilde{K} = 1000, {\tilde{η}}_{c b} = 10^{- 2}$ , ${\tilde{η}}_{V} = 10^{- 4}$ , ${\tilde{l}}_{c} = 0.1$ . In (c), (f), for $δ ζ_{c n}, ζ_{c n} < 0$ the orientation of the director field drawn on the surface (black lines) is set by $- Q_{i j}$ .

When $l_{a} < L_{0}$ , the active region forms an outward cylindrical protrusion (Figure 5a, b and g) whose radius is still well described by Equation 28, replacing $ζ_{c n}$ by $δ ζ_{c n}$ , the value of the active nematic bending moment in the active region (Figure 5g). Using that within the cylindrical protrusion $u = 0$ so that the cylindrical protrusion has the same area as the original active domain and the relation Equation 85 for the size of the active domain, we find that the length of the active protrusion is now given by

L_{c} ≃ \frac{R_{0}^{2}}{R_{c}} (1 - \cos \frac{l_{a}}{R_{0}}) = - \frac{δ ζ_{c n} R_{0}^{2}}{2 κ} (1 - \cos \frac{l_{a}}{R_{0}}),

which is again in excellent agreement with numerical simulation for large $| δ ζ_{c n} |$ and for different values of $l_{a} / L_{0}$ (Figure 5g).

For $ζ_{c n} > 0$ and $l_{a} / L_{0} = 1$ we find erythrocyte-like shapes, where the indentations at the poles become stronger with $ζ_{c n}$ until the two poles touch (Figure 5b). This behaviour remains for $l_{a} < L_{0}$ , resulting in a self-intersection line in the phase diagram (Figure 5a). Here, the shape can take the form of an inner tube entering the spherical shell (Figure 5b), reminiscent of epithelial shape changes observed during sea urchin gastrulation (Ettensohn, 1984).

Interestingly, when $l_{a} / L_{0} < 1$ and the volume is free to change, both signs of $δ ζ_{c n}$ result in a cylindrical appendage forming from the active region. The sign of $δ ζ_{c n}$ determines whether the cylinder forms outside or inside of the remaining, roughly spherical shape. Dynamics simulations confirm that the shapes described above are stable solutions (Figure 5c). At the tip of the emerging cylinder lies the +1 topological defect. For $δ ζ_{c n} < 0$ , when the protrusion grows towards the outside, such a situation is reminiscent of the observation of nematic defects in Hydra, where a set of topological defects, with +1 defects at the tip, have been observed in growing tentacles (Maroudas-Sacks et al., 2021). There, actin layers are perpendicular to each other, with circumferential alignment in the inner cell layer and longitudinal in the outer layer, which would indeed result in $δ ζ_{c n} < 0$ with our sign convention if the layers are contractile.

We now describe surfaces with fixed volume (Figure 5d–f). Here, we do not observe cylindrical shapes or protrusions as in the case of free volume. When $ζ_{c n} < 0$ and $l_{a} = L_{0}$ the surface becomes spindle-like, narrowing at the poles with increasing $| ζ_{c n} |$ . As in the free volume case, when $ζ_{c n} > 0$ the two opposite poles come in contact with each other (Figure 5e); such that subsequent fusion of the poles would lead to an overall toroidal shape of the shell. The shapes become more complex for $l_{a} < L_{0}$ . Shape transitions occur at large $| δ ζ_{c n} |$ , for both $δ ζ_{c n} < 0$ and $δ ζ_{c n} > 0$ (Figure 5e). In the case $δ ζ_{c n} < 0$ , for increasing magnitude of the active bending moment, the shape becomes increasingly curved at the boundary between the passive and active regions, until the solution is lost. In the case $δ ζ_{c n} > 0$ , the shell indents within the active region and the solution branch has a fold. To the right of the fold line in the shape diagram, the steady-state solutions are eventually lost through the formation of a small neck that separates off a smaller, internalised compartment. In contrast to the case of isotropic bending moments, here the sign of $δ ζ_{cn}$ determines whether the active region folds inwards or outwards, independent of the initial size $l_{a} / L_{0}$ . As before, we use dynamics simulations to study the deformations for large $| δ ζ_{cn} |$ (Figure 5f). For both signs of $δ ζ_{cn}$ , these result in shapes that are self-intersecting either along a circle ( $l_{a} / L_{0} = 0.3, δ ζ_{cn} = - 150 κ / R_{0}$ ) or at the poles ( $l_{a} / L_{0} = 0.5, δ ζ_{cn} = 50 κ / R_{0}$ ).

Discussion

In this study of deformations of patterned nematic active surfaces, we have found a diverse zoology of possible shape changes (Figure 6), characterised by budding and neck constrictions, transition of sphere to cylinder, tubulation, and flattening. We find that introduction of a nematic field on the surface greatly increases the space of possible shapes. Overall our work contributes to the characterisation of the ‘morphospace’ which biological systems can explore.

Figure 6

Download asset Open asset

Summary of shape changes obtained through patterning of isotropic and anisotropic active tensions and bending moments.

Active tensions and bending moments are present only in the red region of the surface. For $ζ_{c n} < 0$ the director field orientation (black lines) is set by $- Q_{i j}$ .

Some of our findings recapitulate epithelial deformations observed in biological systems. The flattening observed for an extensile homogeneous nematic surface (Figure 4b, $l_{a} / L_{0} = 1$ ) could in principle lead to merging of the two apposed surfaces into a double-layer for large $| ζ_{n} |$ . Such a process of tissue planarisation appears to occur as an intermediate step in skin organoid formation, where epithelial cysts fuse and merge to form transient bilaterally symmetric structures (Lei et al., 2017). The formation of tubular appendages from nematic bending moments appears to recapitulate growth/regeneration of elongated bodies and tentacles in Hydra (Maroudas-Sacks et al., 2021) and, with an opposite sign, of epithelial invagination during sea urchin embryo gastrulation (Ettensohn, 1984).

The axisymmetric structure we have considered here naturally gives rise to two +1 nematic defects at the poles (Figure 3a). These defects then structure the nematic field and, as a result, the shape changes driven by nematic active tension or bending moments. Such an interplay between topological defect and shape changes is a recurring theme that may play a key role in morphogenesis (Frank and Kardar, 2008; Metselaar et al., 2019; Hoffmann et al., 2021; Blanch-Mercader et al., 2021a; Blanch-Mercader et al., 2021b). In practice +1 nematic defects are unstable to separation into two +1/2 defects; however, it is conceivable that a polar or additional weakly polar field stabilises the +1 defects (Amiri et al., 2022). Extension of the present work beyond axisymmetric structures will allow to distinguish more clearly the purely nematic and polar cases.

Continuum theories for curved surfaces, such as the Helfrich theory, have been extremely successful to describe shape transformations of passive vesicles, including homogeneous or phase-separated vesicles with coexisting domains (Seifert et al., 1991; MacKintosh and Lubensky, 1991; Jülicher and Lipowsky, 1993; Seifert, 1997; Allain et al., 2004; Sens and Turner, 2004; Bassereau et al., 2014). The effect of broken symmetry variables on passive surfaces, arising, for instance, from molecular tilt giving rise to polar order on a lipid membrane, has been considered theoretically (MacKintosh and Lubensky, 1991; Lubensky and Prost, 1992; Park et al., 1992). Continuum theories of active surfaces can similarly allow to study epithelial deformations (Salbreux and Jülicher, 2017; Morris and Rao, 2019; Messal et al., 2019). We note some important differences between the active surface model described here and passive membranes. (i) Our constitutive equations for tensions and bending moments Equations 5 and 6 do not in general derive from a free energy (Salbreux and Jülicher, 2017) and describe a system out-of-equilibrium; (ii) while lipid membranes are nearly incompressible and are usually treated as surfaces with constant area, cells within epithelial tissues can change their area significantly (Latorre et al., 2018), which prompted us to consider a finite area modulus $K$ : for example, simulations with constant volume have relative area changes of up to 20% (Figure 2—figure supplement 2); (iii) patterns of active tensions and bending moments imposed here also do not derive from an energy and are thought to respond to spatiotemporal chemical cues: in contrast, phase-separated domains in passive lipid vesicles obey equilibrium thermodynamics and their size is controlled, for instance, by line tension at the domain boundary (Jülicher and Lipowsky, 1993). In some cases, however, a similarity appears between shape transformations obtained in the active model we study here and the passive Helfrich model. For instance, budding occurring in lipid membranes due to phase separation of domains with different spontaneous curvature (Jülicher and Lipowsky, 1993) is similar to the budding we observe here for different regions with different active isotropic bending moments.

We find here that nematically oriented active bending moments can give rise to spontaneous cylindrical tubes, without external force application (Figure 5). Spontaneous formation of hollow cylindrical vesicles with polar order due to molecular tilt has been discussed Lubensky and Prost, 1992; there the cylindrical shapes are considered to be open and the gain in defect energy allows the open cylinder to be more stable than the spherical shape. In contrast, we find here active surfaces which spontaneously form tubes, but stay closed and keep their topological charge. It has also been reported that a supported bilayer membrane under compression can spontaneously form tubes under negative tension (Staykova et al., 2013). In this work we have chosen to consider only positive isotropic tension; negative isotropic tension could give rise to further buckling instabilities. Models for chiral lipid bilayers in a tilted fluid phase have also predicted tubular shapes (Helfrich and Prost, 1988; Selinger and Schnur, 1993; Selinger et al., 1996; Tu and Seifert, 2007). Here, we have not considered chiral effects. These effects could be introduced by generalising the constitutive Equations 5 and 6, including terms which appear for surfaces with broken planar-chiral or chiral symmetry (Salbreux and Jülicher, 2017).

In contrast to purely elastic models of morphogenesis (Höhn et al., 2015; Haas et al., 2018), we have considered here morphogenetic events occurring on time scales long enough for shear elastic stresses to be relaxed by cell topological rearrangements, such that the tissue exhibits fluid behaviour (Popović et al., 2017). Whether a tissue behaves as an elastic or fluid material on time scales relevant to morphogenesis can in principle be probed experimentally (Mongera et al., 2018).

While we have focused the interpretation of our results to epithelial mechanics, the constitutive Equations 5 and 6 we have considered here are generic and may also describe the large-scale behaviour of active nematics formed with cytoskeletal filaments and motors on a deformable surface (Keber et al., 2014). We considered here, however, a situation where the two-dimensional fluid has area elasticity, whereas cytoskeletal networks can in principle be fluid with respect to both shear and bulk shear due to the turnover of components.

In this study, we have considered chemical and mechanical processes to be uncoupled, except for the profile of active tension or torque being advected with the surface flow. Introducing additional couplings explicitly in this framework will extend the repertoire of shapes considered here. A natural choice is to consider the effect of a chemical undergoing reaction-diffusion on the surface and advected by the fluid, regulating active forces on the surface (Mietke et al., 2019a; Mietke et al., 2019b). Here, we assumed that orientational order relaxes quickly compared to other dynamical processes; in future work, this assumption could be lifted and one could study in particular how chemical regulation could influence the dynamics of orientational order in the tissue. Cells could also be sensing their own curvature and actively adapt their behaviour accordingly (Chen et al., 2019), which could lead to a dependency of the active coupling coefficients $ζ$ , $ζ_{n}$ , $ζ_{c}$ or $ζ_{c n}$ on the trace or determinant of the curvature tensor $C_{i j}$ . It would be interesting to explore shapes arising from such a feedback. Volume conservation at cellular level could also be included explicitly, for instance, by introducing a tissue height field (Morris and Rao, 2019). Finally, we have considered here a tissue with a fixed preferred area, implicitly assuming that the epithelium is not growing. Tissue growth is a key aspect of biological development (Gokhale and Shingleton, 2015; Eder et al., 2017), and cell division and death can fluidify elastic stresses in an epithelium (Ranft et al., 2010); adding regulated growth in the model will be a step forward in our understanding of active morphogenesis of biological tissues.

Appendix 1

Differential geometry of axisymmetric surfaces

General relations of differential geometry

Fundamental tensors

A general framework for the mechanics of active curved surfaces is given in Salbreux and Jülicher, 2017; Salbreux et al., 2022, and we follow the differential geometry notation introduced there. Let $X = X (s^{1}, s^{2})$ be a curved surface embedded in $ℝ^{3}$ and parametrised by the generalised coordinates $s^{1}, s^{2}$ . A local covariant basis in the tangent plane is given by the vectors $e_{i} = \partial_{i} X = \partial X / \partial s^{i}$ and the unit normal vector is constructed as $n = (e_{1} \times e_{2}) / | e_{1} \times e_{2} |$ and chosen to point outwards for a closed surface in our convention. These define the metric tensor $g_{i j} = e_{i} \cdot e_{j}$ and the curvature tensor $C_{i j} = - (\partial_{i} e_{j}) \cdot n = e_{i} \cdot \partial_{j} n$ . The infinitesimal surface and line elements are given by $dS = \sqrt{g} d s^{1} d s^{2}$ and $d l^{2} = g_{i j} d s^{i} d s^{j}$ , where $g$ is the determinant of the metric tensor. The antisymmetric Levi–Civita tensor is defined as $ϵ_{i j} = n \cdot (e_{i} \times e_{j})$ .

Covariant derivatives

The Christoffel symbols of the second kind are obtained from derivatives of the metric as

Γ_{i j}^{k} = \frac{1}{2} g^{k m} (\partial_{j} g_{i m} + \partial_{i} g_{j m} - \partial_{m} g_{i j}),

and Christoffel symbols of the first kind are defined as $Γ_{k i j} = \frac{1}{2} (\partial_{j} g_{k i} + \partial_{i} g_{k j} - \partial_{k} g_{i j})$ . The covariant derivatives of a tangent vector field $f^{i}$ and a tangent tensor field $T^{i j}$ , respectively, are given by

\nabla_{i} f^{j} = \partial_{i} f^{j} + Γ_{i k}^{j} f^{k},

\nabla_{i} T^{j k} = \partial_{i} T^{j k} + Γ_{i l}^{j} T^{l k} + Γ_{i l}^{k} T^{j l}

In the following, we also use the divergence theorem on curved surfaces for a tangent vector field $f$ (Salbreux and Jülicher, 2017),

\int_{S} d S \nabla_{i} f^{i} = \int_{C} d l ν_{i} f^{i} .

Infinitesimal variation of surface quantities

For a small deformation of the surface

δ X = δ X^{i} e_{i} + δ X_{n} n,

the variations of the basis vectors, the normal vector, the metric, and the mixed curvature tensor components are given by (Salbreux and Jülicher, 2017)

δ e_{i} = (\nabla_{i} δ X^{j} + C_{i}^{j} δ X_{n}) e_{j} + (\partial_{i} δ X_{n} - C_{i j} δ X^{j}) n,

δ n = (- \partial_{i} δ X_{n} + C_{i j} δ X^{j}) e^{i},

δ g_{i j} = \nabla_{i} δ X_{j} + \nabla_{j} δ X_{i} + 2 C_{i j} δ X_{n},

\frac{δ \sqrt{g}}{\sqrt{g}} = \frac{1}{2} g^{i j} δ g_{i j} = \nabla^{k} δ X_{k} + C_{k}^{k} δ X_{n},

δ C_{i}^{j} = - \nabla_{i} (\partial^{j} δ X_{n}) + (\nabla_{i} δ X^{k}) C_{k}^{j} - (\nabla^{k} δ X^{j}) C_{i k} + (\nabla_{k} C_{i}^{j}) δ X^{k} - δ X_{n} C_{i k} C^{k j} .

These relations can be used to obtain time derivatives of surface quantities using $δ X = δ t v$ (Lagrangian surface update) or $δ X = δ t v_{n} n$ (where the surface shape is updated with the normal flow only).

Since according to the definition Equation 14 the area strain $u$ can be written as, using the coordinates $(s_{0}, ϕ)$ on the undeformed and deformed surfaces and denoting $g_{0}$ the determinant of the metric of the undeformed surface,

u = \frac{\sqrt{g} - \sqrt{g_{0}}}{\sqrt{g_{0}}},

we obtain from Equation 38 the variation:

δ u = (1 + u) (\nabla_{k} δ X^{k} + C_{k}^{k} δ X_{n}),

which yields the Lagrangian time derivative (Equation 15) in the main text, using $δ X = δ t v$ .

Axisymmetric surfaces

Fundamental tensors

On a surface with axial symmetry about the $z$ -axis, as defined in the main text, the basis vectors and the outward normal are

e_{ϕ} = (\begin{matrix} - x \sin ϕ \\ x \cos ϕ \\ 0 \end{matrix}),

e_{s} = (\begin{matrix} \cos ϕ \cos ψ \\ \sin ϕ \cos ψ \\ \sin ψ \end{matrix}),

n = (\begin{matrix} \cos ϕ \sin ψ \\ \sin ϕ \sin ψ \\ - \cos ψ \end{matrix}),

where we have used $\partial_{s} x = \cos ψ$ and $\partial_{s} z = \sin ψ$ , which can be defined through the requirement that $s$ is an arc length parameter, such that ${| e_{s} |}^{2} = {(\partial_{s} x)}^{2} + {(\partial_{s} z)}^{2} = 1$ . The metric and curvature tensors and the surface element are given by

g_{i j} = (\begin{matrix} x^{2} & 0 \\ 0 & 1 \end{matrix}), C_{i}^{j} = (\begin{matrix} \frac{\sin ψ}{x} & 0 \\ 0 & \partial_{s} ψ \end{matrix}), d S = x d s d ϕ .

In the following, because the metric is diagonal, we will not distinguish between the order of indices for diagonal elements of second-order tensors in mixed coordinates, that is, for a tensor $T$ we use $T_{s}^{s} = T_{s}^{s} = T_{s}^{s}$ and $T_{ϕ}^{ϕ} = T_{ϕ}^{ϕ} = T_{ϕ}^{ϕ}$ . The circumferential and meridional principal curvatures $C_{ϕ}^{ϕ}$ and $C_{s}^{s}$ , and the mean and Gaussian curvatures $H$ and $K$ are given by:

C_{ϕ}^{ϕ} = \frac{\sin ψ}{x}, C_{s}^{s} = \partial_{s} ψ,

H = \frac{1}{2} C_{k}^{k}, K = det C_{i}^{j} = C_{ϕ}^{ϕ} C_{s}^{s} .

Some useful relationships involving the two principal curvatures follow from Equation 46 and from the definitions Equations 11 and 12,

\partial_{s} C_{ϕ}^{ϕ} = \frac{\cos ψ}{x} (C_{s}^{s} - C_{ϕ}^{ϕ}),

\partial_{s} (\frac{\cos ψ}{x}) = - C_{ϕ}^{ϕ} C_{s}^{s} - {(\frac{\cos ψ}{x})}^{2} .

The partial area and partial volume are given by

a (s) = 2 π \int_{0}^{s} d s^{'} x (s^{'}),

v (s) = π \int_{0}^{s} d s^{'} x (s^{'})^{2} \sin ψ (s^{'}) .

The corotational time derivative of the curvature tensor, as defined in Equation 8, has trace

\frac{D C_{k}^{k}}{D t} = - \partial_{s}^{2} v_{n} - \frac{\cos ψ}{x} \partial_{s} v_{n} - v_{n} ((C_{s}^{s})^{2} + (C_{ϕ}^{ϕ})^{2}) + v_{s} \partial_{s} C_{k}^{k} .

Covariant derivatives

Axial symmetry implies that all functions on the surface should be $ϕ$ -independent, and we also consider here vector and tensor fields $f$ , $T$ such that $f^{ϕ} = 0$ and $T^{s ϕ} = T^{ϕ s} = 0$ . The only non-vanishing component of $\partial_{i} g_{j k}$ is $\partial_{s} g_{ϕ ϕ} = 2 x \cos ψ$ ; therefore, the non-zero Christoffel symbols of the first kind are

Γ_{ϕ s ϕ} = Γ_{ϕ ϕ s} = - Γ_{s ϕ ϕ} = \frac{1}{2} \partial_{s} g_{ϕ ϕ} = x \cos ψ,

and the non-zero Christoffel symbols of the second kind are

Γ_{ϕ ϕ}^{s} = - x \cos ψ, Γ_{s ϕ}^{ϕ} = Γ_{ϕ s}^{ϕ} = \frac{\cos ψ}{x} .

The non-zero components of $\nabla_{i} f^{j}$ are

\nabla_{s} f^{s} = \partial_{s} f^{s},

\nabla_{ϕ} f^{ϕ} = \frac{\cos ψ}{x} f^{s},

resulting in the components of the strain rate tensor defined in Equation 7:

v_{s}^{s} = \partial_{s} v^{s} + C_{s}^{s} v_{n},

v_{ϕ}^{ϕ} = \frac{\cos ψ}{x} v^{s} + C_{ϕ}^{ϕ} v_{n},

with trace

v_{k}^{k} = \partial_{s} v^{s} + \frac{\cos ψ}{x} v^{s} + C_{k}^{k} v_{n} .

The non-zero components of $\nabla_{i} T^{j k}$ are

\nabla_{s} T^{s s} = \partial_{s} T^{s s},

\nabla_{s} T^{ϕ ϕ} = \partial_{s} T^{ϕ ϕ} + 2 \frac{\cos ψ}{x} T^{ϕ ϕ},

\begin{array}{rcl} \nabla_{ϕ} T^{s ϕ} & = & - x \cos ψ T^{ϕ ϕ} + \frac{\cos ψ}{x} T^{s s} \\ = & \frac{\cos ψ}{x} (T_{s}^{s} - T_{ϕ}^{ϕ}) \\ = & \nabla_{ϕ} T^{ϕ s}, \end{array}

where $T^{ϕ ϕ} = g^{ϕ ϕ} T_{ϕ}^{ϕ} = \frac{1}{x^{2}} T_{ϕ}^{ϕ}$ was used in the last expression.

Force and torque balance

Axial symmetry implies that the tangential and the normal force balances, Equations 1; 2, and the torque balance Equation 3 can be rewritten as

\partial_{s} t_{s}^{s} = \frac{\cos ψ}{x} (t_{ϕ}^{ϕ} - t_{s}^{s}) - C_{s}^{s} t_{n}^{s} - f_{s}^{e x t},

\partial_{s} t_{n}^{s} = C_{ϕ}^{ϕ} (t_{ϕ}^{ϕ} - t_{s}^{s}) + C_{k}^{k} t_{s}^{s} - \frac{\cos ψ}{x} t_{n}^{s} - f_{n}^{e x t} - P,

\partial_{s} {\bar{m}}_{s}^{s} = \frac{\cos ψ}{x} ({\bar{m}}_{ϕ}^{ϕ} - {\bar{m}}_{s}^{s}) + t_{n}^{s} .

The geometric singularities appearing in Equations 63–65 are removed by an appropriate choice of boundary conditions for the tensions and moments at the poles of the surface. The normal torque balance Equation 4 gives the antisymmetric part of the tension tensor. With the constitutive Equation 6 for the bending moment tensor, $ϵ_{i j} t^{i j} = C_{i j} m^{i j} = ζ_{c n} Q^{i k} ϵ_{k} {}^{j}C_{i j}$ which vanishes for an axisymmetric surface where $Q^{s ϕ} = Q^{ϕ s} = 0$ . Therefore, here $t^{i j} = t_{s}^{i j}$ which is given by the constitutive Equation 5.

Direct expression for the transverse tension on an axisymmetric surface

In Capovilla and Guven, 2002 and Knoche and Kierfeld, 2011, it is shown that on axially symmetric surfaces the normal force balance in Equation 64 can be integrated in a closed form in the presence of a uniform pressure. Here, we generalise this to an arbitrary axially symmetric external force. In analogy to Capovilla and Guven, 2002, consider a piece of surface $S_{1}$ bounded by the south pole and a circle $C$ perpendicular to the axis of symmetry, given by $s = s_{1}$ (Appendix 1—figure 1). The bounding circle $C$ has the line element $d l = x d ϕ$ and the unit normal $v = e_{s}$ , tangent to the surface and pointing outward with respect to $S_{1}$ . The balance of forces acting on $S_{1}$ reads

\oint_{C} d l ν_{i} t^{i} + \int_{S_{1}} d S (f^{e x t} + P n) = 0 .

The contributions to the integral Equation 66 are

\begin{array}{rcl} \oint_{C} d l ν_{i} t^{i} & = & \int_{0}^{2 π} d ϕ x t^{s} \\ = & \int_{0}^{2 π} d ϕ x (t^{s s} e_{s} + t_{n}^{s} n) \\ = & 2 π x (t^{s s} \sin ψ - t_{n}^{s} \cos ψ) e_{z}, \end{array}

\begin{array}{rcl} P \int_{S_{1}} d S n & = & (- 2 π P \int_{0}^{s_{1}} d s x \cos ψ) e_{z} \\ = & (- 2 π P \int_{0}^{x (s_{1})} d x x) e_{z} \\ = & - π P x^{2} e_{z}, \end{array}

\begin{array}{rcl} \int_{S_{1}} d S f^{e x t} & = & 2 π (\int_{0}^{s_{1}} d s x f_{z}^{e x t}) e_{z} \\ = & 2 π I (s_{1}) e_{z}, \end{array}

where we have introduced the integrated external force

I (s) = \int_{0}^{s} d s^{'} x f_{z}^{e x t},

and used the shape Equation 11 in Equation 68. As Equations 67–69 only contribute to the $z$ -component, Equation 66 can be rewritten as

2 π x (t_{s}^{s} \sin ψ - t_{n}^{s} \cos ψ) - π P x^{2} + 2 π I = 0.

From Equation 71 one obtains an expression for the transverse tension for $ψ \neq \frac{π}{2}$ :

t_{n}^{s} = t_{s}^{s} \tan ψ - \frac{1}{2} \frac{x}{\cos ψ} P + \frac{1}{x \cos ψ} I,

and it is easy to confirm, using Equation 63, that this is indeed a solution of the normal force balance given by Equation 64.

Appendix 1—figure 1

Download asset Open asset

Schematic of the surface $S_{1}$ used to derive the integral of the normal force balance.

Behaviour at the poles

The poles of the axisymmetric surface are at $s = 0$ (south pole) and $s = L$ (north pole) and satisfy $x (0) = x (L) = 0$ . Besides, we assumed that the shape has a finite curvature, requiring $ψ = 0$ and $ψ = π$ at the south and north poles. The asymptotic behaviour of the shape is then:

\begin{array}{rcl} x (s) & = & x (0) + \cos ψ |_{s = 0} s - \frac{1}{2} ((\sin ψ) \partial_{s} ψ) |_{s = 0} s^{2} + O (s^{3}) \\ = & s + O (s^{3}), \end{array}

\begin{array}{rcl} x (L - s) & = & x (0) + \cos ψ |_{s = L} (L - s) - \frac{1}{2} ((\sin ψ) \partial_{s} ψ) |_{s = L} (L - s)^{2} + O ((L - s)^{3}) \\ = & s - L + O ((L - s)^{3}) . \end{array}

The limits of geometric singularities of the form $f (s) / x (s)$ , for some function $f (s)$ which vanishes at the pole, follow from L’Hôpital’s rule

lim_{s \to 0, L} \frac{f (s)}{x (s)} = lim_{s \to 0, L} \frac{\partial_{s} f |_{s = 0, L}}{\cos ψ (s) |_{s = 0, L}} = \pm \partial_{s} f |_{s = 0, L} .

where the + sign applies to $s = 0$ and the - sign to $s = L$ . For example, applying Equation 75 to $C_{ϕ}^{ϕ} = \sin ψ / x$ yields that

lim_{s \to 0, L} C_{ϕ}^{ϕ} = C_{s}^{s} |_{s = 0, L} .

Any smooth tangent vector field on the closed axisymmetric surface has to vanish at the poles. For example, since $C_{k}^{k}$ is a scalar field, for its derivative we have

\partial_{s} C_{k}^{k} |_{s = 0, L} = 0.

From Equation 48 we find that at the poles $2 \partial_{s} C_{ϕ}^{ϕ} = \partial_{s} C_{s}^{s}$ , which, together with Equation 77 yields

\partial_{s} C_{ϕ}^{ϕ} |_{s = 0, L} = \partial_{s} C_{s}^{s} |_{s = 0, L} = 0.

This relation, together with Equation 76, implies that at the poles, the surface is locally spherical.

Spherical surface

We give here some of the geometrical quantities defined above for the undeformed initial surface, a sphere with radius $R_{0}$ and south pole at the origin, $X (ϕ, 0) = 0$ :

L_{0} = π R_{0}, 0 \leq s \leq π R_{0},

x (s_{0}) = R_{0} \sin \frac{s_{0}}{R_{0}},

z (s_{0}) = R_{0} (1 - \cos \frac{s_{0}}{R_{0}}),

ψ (s_{0}) = \frac{s_{0}}{R_{0}},

C_{s}^{s} = C_{ϕ}^{ϕ} = H = \frac{1}{R_{0}},

\frac{\cos ψ}{x} = \frac{1}{R_{0} \tan \frac{s_{0}}{R_{0}}},

a (s_{0}) = \frac{A_{0}}{2} (1 - \cos \frac{s_{0}}{R_{0}}),

v (s_{0}) = V_{0} (2 + \cos \frac{s_{0}}{R_{0}}) \sin^{4} \frac{s_{0}}{2 R_{0}},

where the arc length is denoted $s_{0}$ , and $A_{0} = 4 π R_{0}^{2}$ , $V_{0} = (4 / 3) π R_{0}^{3}$ are the area and volume of the sphere.

Appendix 2

Infinitesimal neck

We discuss here the infinitesimal neck appearing in steady-state shapes subjected to isotropic active bending moments ( $ζ = ζ_{n} = ζ_{c n} = 0$ ). In that case, the tensors $t^{i j}$ and ${\bar{m}}^{i j}$ are isotropic and the force and torque balance Equations 63–65 can be written, in the absence of external force other than the pressure $P$ and for $ψ \neq \frac{π}{2}$ :

\partial_{s} (\frac{t_{s}^{s}}{\cos ψ}) = \frac{x \partial_{s} ψ}{2 \cos^{2} ψ} P,

\partial_{s} {\bar{m}}_{s}^{s} = t_{s}^{s} \tan ψ - \frac{1}{2} \frac{x}{\cos ψ} P,

where we have used the transverse tension solution (Equation 72).

Scaling analysis

To analyse the behaviour of these equations near an infinitesimal neck, we now perform a scaling analysis following Fourcade et al., 1994. We consider a region around a nearly closed neck with minimal radius $a$ . At the point of the surface closest to the axis of symmetry, $x = a$ and $ψ = \frac{π}{2}$ . We then scale the arc length coordinate $s$ , the distance of the surface to the axis of symmetry $x$ , and the curvature tensor with $a$ , and introduce $\bar{s} = s / a$ , $\bar{x} = x / a$ and ${\bar{C}}_{k}^{k} = a C_{k}^{k}$ . The force balance equations then become for $ψ \neq \frac{π}{2}$ :

\partial_{\bar{s}} (\frac{t_{s}^{s}}{\cos ψ}) = \frac{a \bar{x} \partial_{\bar{s}} ψ}{2 \cos^{2} ψ} P,

\partial_{\bar{s}} (2 κ {\bar{C}}_{k}^{k}) = - a \partial_{\bar{s}} ζ_{c} + a^{2} t_{s}^{s} \tan ψ - \frac{a^{3}}{2} \frac{\bar{x}}{\cos ψ} P .

For $a \to 0$ , the leading order solution has $t_{s}^{s} / \cos ψ$ and ${\bar{C}}_{k}^{k}$ both constant. Using the relation

{\bar{C}}_{k}^{k} = \partial_{\bar{s}} ψ + \frac{\sin ψ}{\bar{x}} = \frac{1}{\bar{x}} \partial_{\bar{x}} ({\bar{x}}^{2} {\bar{C}}_{ϕ}^{ϕ}),

with ${\bar{C}}_{ϕ}^{ϕ} = a C_{ϕ}^{ϕ}$ , and the conditions ${\bar{C}}_{ϕ}^{ϕ} (\bar{x} = 1) = 1$ and that $\sin ψ$ does not diverge for $| \bar{x} | \to \infty$ , the curvatures have solution ${\bar{C}}_{k}^{k} = 0$ and

C_{ϕ}^{ϕ} = - C_{s}^{s} = \frac{a}{x^{2}}, \cos ψ = \pm \sqrt{1 - \frac{a^{2}}{x^{2}}},

where the sign of $\cos ψ$ changes in the regions towards and away from the neck. Therefore, $\cos ψ$ converges to +1 or -1 away from the neck for $x \to \pm \infty$ . Since $t_{s}^{s} / \cos ψ$ is constant across the infinitesimal neck, $t_{s}^{s}$ also changes sign asymptotically away from the neck.

The next order in $a$ of Equation 90 gives $\partial_{\bar{s}} (2 κ C_{k}^{k} + ζ_{c}) = 0$ which corresponds to ${\bar{m}}_{s}^{s}$ constant across the neck (Figure 2—figure supplement 1).

Analytical solution for the free volume case

In the free volume case, $P = 0$ and the force balance equations admits the solution $t_{s}^{s} = u = 0$ and constant ${\bar{m}}_{s}^{s}$ (Figure 2—figure supplement 1a). Considering a shape with the active and passive regions forming spheres with radii $R_{a}$ and $R_{p}$ separated by an infinitesimal neck, the condition of constant ${\bar{m}}_{s}^{s}$ results in

2 κ [C_{k}^{k}]_{a} + δ ζ_{c} = 2 κ [C_{k}^{k}]_{p},

with $[C_{k}^{k}]_{a}$ and $[C_{k}^{k}]_{p}$ the trace of the curvature tensor in the active and passive regions, and $δ ζ_{c}$ the difference in isotropic active bending moment between the passive and active regions. This results in Equation 24, taking into account that $[C_{k}^{k}]_{a} = \pm 2 / R_{a}$ depending on whether the active region is curved towards the outside or the inside part of the surface. In addition, the condition $u = 0$ results in conservation of area of the active and passive regions compared to the undeformed sphere:

\begin{array}{ll} 4 π R_{a}^{2} = 2 π R_{0}^{2} (1 - \cos \frac{l_{a}}{R_{0}}), \\ 4 π R_{p}^{2} = 2 π R_{0}^{2} (1 + \cos \frac{l_{a}}{R_{0}}), \end{array}

where we have used Equation 85. When $[C_{k}^{k}]_{a} = 2 / R_{a}$ corresponding to the active region towards the outside, and $δ ζ_{c} > 0$ , Equation 24 implies that $R_{p} < R_{a}$ which further requires $l_{a} / L_{0} > 1 / 2$ to satisfy Equation 94.

Combining Equations 93 and 94 gives a condition defining a curve in the parameter space $l_{a} / L_{0}$ , $δ ζ_{c} R_{0} / κ$ :

\frac{δ ζ_{c} R_{0}}{4 κ} = \sqrt{\frac{2}{1 + \cos (\frac{π l_{a}}{L_{0}})}} \mp \sqrt{\frac{2}{1 - \cos (\frac{π l_{a}}{L_{0}})}},

which agrees well with the line of neck constriction determined numerically (Figure 2b).

Appendix 3

Nematic order parameter on an axisymmetric surface

Equilibrium equation

A nematic director on a curved surface is given by a unit tangent vector, for which $\hat{n} = - \hat{n}$ . In an orthonormal frame ${{\hat{e}}_{ϕ}, {\hat{e}}_{s}}$ , where ${\hat{e}}_{ϕ} = e_{ϕ} / x$ and ${\hat{e}}_{s} = e_{s}$ , it is characterised by an angle $α \in [0, π]$ as

\hat{n} = \cos α {\hat{e}}_{ϕ} + \sin α {\hat{e}}_{s} .

The director components in the basis ${e_{ϕ}, e_{s}}$ are then

{\hat{n}}^{ϕ} = \frac{\cos α}{x}, {\hat{n}}^{s} = \sin α .

The traceless and symmetric nematic order parameter $Q^{i j}$ can be constructed from the director $n^{i}$ and a magnitude $S$ as

Q^{i j} = S ({\hat{n}}^{i} {\hat{n}}^{j} - \frac{1}{2} g^{i j}) .

Its components read

\begin{array}{ll} Q^{ϕ ϕ} = \frac{S}{2} \frac{\cos 2 α}{x^{2}}, \\ Q^{s s} = - \frac{S}{2} \cos 2 α, \\ Q^{s ϕ} = Q^{ϕ s} = \frac{S}{2} \frac{\sin 2 α}{x} . \end{array}

We assume here that there is no azimuthal flow on the surface, $v^{ϕ} = 0$ . The $ϕ$ -component of the tangential force balance then reads for constant $ζ_{n} \neq 0$ , $ζ_{c n} \neq 0$ :

\nabla_{i} t^{i ϕ} = (ζ_{n} + ζ_{c n} \frac{\sin ψ}{x}) (\partial_{s} Q^{s ϕ} + 3 \frac{\cos ψ}{x} Q^{s ϕ}) = 0,

which requires $Q^{s ϕ} = Q^{ϕ s} = 0$ excluding divergence of $Q^{s ϕ}$ at the poles; therefore, the only possible orientations for the director, compatible with our assumption of vanishing azimuthal flows, are $α = 0, π / 2$ . Therefore, in the axisymmetric setup we consider only one non-zero component

q = \frac{S}{2} \cos 2 α = Q_{ϕ}^{ϕ} = - Q_{s}^{s},

which can take the values $q = \pm S / 2$ , corresponding to azimuthal or longitudinal orientation of the director $\hat{n}$ , respectively.

The total free energy of the nematic Equation 9 reads in terms of $q$

F = \int d S f = \int d S (k ({(\partial_{s} q)}^{2} + 4 {(\frac{\cos ψ}{x} q)}^{2}) - \frac{a}{2} q^{2} + \frac{a}{4} q^{4}) .

We minimise this energy with respect to $q$ on a given surface. The resulting Euler–Lagrange Equation 16 is obtained from

0 = \frac{δ F}{δ q} = h - \nabla_{i} Π^{i},

where

h = \frac{\partial f}{\partial q}, Π^{i} = \frac{\partial f}{\partial (\partial_{i} q)} .

Here, we have used the definition from the functional derivative, $d F ≃ \int dS \frac{δ F}{δ q} d q$ . Equation 16 can be written as two first-order equations

\partial_{s} q = w,

\partial_{s} w = \frac{1}{2 l_{c}^{2}} q (q^{2} - 1) + \frac{\cos ψ}{x} (4 \frac{\cos ψ}{x} q - w) .

The requirement that Equation 106 should be regular at the poles of the surface results in the two boundary conditions

q (0) = q (L) = 0,

and also implies

w (0) = w (L) = 0.

The limit of Equation 106 at $s = 0$ is given by

\begin{array}{ll} \partial_{s} w (0) & = \frac{1}{2 l_{c}^{2}} q (0) (q (0)^{2} - 1) + lim_{s \to 0} [\frac{\cos (ψ)}{x} (4 \frac{\cos (ψ)}{x} q - w)] \\ = \frac{1}{2 l_{c}^{2}} q (0) (q (0)^{2} - 1) + lim_{s \to 0} [\frac{1}{s} (4 \frac{1}{s} (q (0) + w (0) s + \frac{1}{2} \partial_{s} w (0) s^{2}) - (w (0) + \partial_{s} w (0) s))] \\ = lim_{s \to 0} [\frac{1}{s} (3 w (0) + \partial_{s} w (0) s)] \\ = \partial_{s} w (0), \end{array}

and equivalently at $s = L$ . Therefore, Equation 106 does not provide a limit value for $\partial_{s} w$ at the poles. When solving Equations 105 and 106 numerically, we use the analytical limits at the poles

{(\partial_{s} q, \partial_{s} w)}_{s = 0} = (0, W_{0}),

{(\partial_{s} q, \partial_{s} w)}_{s = L} = (0, W_{L}),

with two free parameters $W_{0}$ and $W_{L}$ , which are introduced in order to ensure all four boundary conditions 107 and 108, of which the second two have to be imposed explicitly for numerical reasons.

Stability of the isotropic state on a sphere

We discuss here the stability of the isotropic state $q = 0$ on a sphere of radius $R_{0}$ to axisymmetric perturbations; a more general analysis can be found in Napoli and Vergori, 2012. We note that for a spherical shape, with $θ = s / R_{0}$ and at first order in $q$ :

\frac{δ F}{δ q} ≃ - \frac{2 k}{R_{0}^{2}} [\partial_{θ}^{2} + \cot θ \partial_{θ} - 4 \cot^{2} θ + \frac{R_{0}^{2}}{2 l_{c}^{2}}] q,

where we have used Equation 103 and the geometrical relations for a sphere given in Appendix 1 section ‘Spherical surface’. A set of eigenfunctions of the differential operator $ℒ = \partial_{θ}^{2} + \cot θ \partial_{θ} - 4 \cot^{2} θ$ is provided by taking derivatives of axisymmetric spherical harmonics, $q_{n} (θ) = q_{n} f_{n} (θ)$ with $f_{n} (θ) = [\partial_{θ}^{2} - \cot θ \partial_{θ}] P_{n} (\cos θ)$ with $P_{n}$ the Legendre polynomial of degree $n$ , for $n \geq 2$ . One then finds

\frac{δ F}{δ q} [q_{n}] ≃ \frac{2 k}{R_{0}^{2}} [n (n + 1) - 4 - \frac{R_{0}^{2}}{2 l_{c}^{2}}] q_{n} f_{n} (θ)

The isotropic state is stable if $n (n + 1) - 4 - \frac{R_{0}^{2}}{2 l_{c}^{2}} \geq 0$ for $n \geq 2$ , or for

\frac{l_{c}}{R_{o}} > \frac{1}{2} .

Appendix 4

Active nematic tension and bilayered disc: Asymptotic analysis

We discuss here an asymptotic analysis for the flat bilayers disc steady-state shapes found for surfaces subjected to vanishing internal pressure, uniform nematic active tension, and for $ζ_{n} < 0$ (Figure 4b). We consider here the limit where $ζ = ζ_{c} = ζ_{c n} = 0$ .

We postulate that the limit shape reached as $| ζ_{n} | \to \infty$ consists of two parallel flat central discs of radius $R_{d}$ , separated by a distance $2 h$ , and connected by a narrow curved region (Appendix 4—figure 1). One denotes $s_{t}$ the arclength of the extreme position of the shape where $x = x_{t}$ is maximal, and $s_{c} = s - s_{t}$ the arclength from this point (Appendix 4—figure 2). The shape is assumed to be symmetric about a plane going through the equator, which imposes $\partial_{s} q (s = s_{t}) = 0$ . One looks for an asymptotic shape which satisfies $h ≪ R_{d}$ ; we show later that this condition requires $κ ≪ K l_{c}^{2}$ .

The force and torque balance Equations 63–65 can then be rewritten

\partial_{s} [\frac{t_{s}^{s}}{\cos ψ}] = \frac{2 ζ_{n} q}{x},

2 κ \partial_{s} [\partial_{s} ψ + \frac{\sin ψ}{x}] = t_{s}^{s} \tan ψ .

with $t_{s}^{s} = 2 K u - ζ_{n} q$ . The last equation implies that at $ψ = π / 2$ (i.e. at the point of the surface with the extremal value of $x$ ), $t_{s}^{s} = 0$ and therefore at this point $u = ζ_{n} q / (2 K)$ . However, our definition of deformation implies that $u > - 1$ . Therefore, as $| ζ_{n} | \to \infty$ , one must have $q \to 0$ (Appendix 4—figure 2), and the equilibrium equation for the nematic order parameter $q$ , Equations 105 and 106 can be linearised.

Introducing a renormalised order parameter $\tilde{q} = - ζ_{n} / (2 K) q$ and renormalised tension ${\tilde{t}}_{s}^{s} = t_{s}^{s} / (2 K) = u + \tilde{q}$ , the force and torque balance equation and the linearised equilibrium equation for the order parameter read

\partial_{s} [\frac{{\tilde{t}}_{s}^{s}}{\cos ψ}] = - \frac{2 \tilde{q}}{x},

\frac{κ}{K} \partial_{s} [\partial_{s} ψ + \frac{\sin ψ}{x}] = {\tilde{t}}_{s}^{s} \tan ψ,

\partial_{s}^{2} \tilde{q} + \frac{\cos ψ}{x} \partial_{s} \tilde{q} - 4 \frac{\cos^{2} ψ}{x^{2}} \tilde{q} = - \frac{1}{2 l_{c}^{2}} \tilde{q},

to which one can add the boundary conditions $\tilde{q} (0) = \tilde{q} (L) = 0$ , $ψ (0) = 0$ , $ψ (L) = π$ and a condition on $u$ which follows from Equation 143:

2 R_{0}^{2} = \int_{0}^{L} \frac{d s x}{1 + u} .

In the asymptotic regime where $h ≪ R_{d}$ , we consider separately the flat central discs and the narrow curved connecting region.

Flat central disc

In the lower, flat part of the deformed shape, one has $ψ = 0$ , $x = s$ and Equation 119 becomes

\partial_{s}^{2} \tilde{q} + \frac{1}{s} \partial_{s} \tilde{q} - \frac{4}{s^{2}} \tilde{q} = - \frac{1}{2 l_{c}^{2}} \tilde{q},

with solution, using $\tilde{q} (s = 0) = 0$ , $\tilde{q} (x) = C_{q} J_{2} (\frac{x}{\sqrt{2} l_{c}})$ , with $C_{q}$ a constant to determine, and $J_{n} (x)$ are the Bessel functions of the first kind. The condition $0 = \partial_{s} q (s = s_{t}) ≃ \partial_{s} q (s = R_{d})$ then yields the expression for the radius of the disc:

R_{d} = \sqrt{2} l_{c} β_{0},

with $β_{0}$ the smallest positive solution of $J_{2}^{'} (β_{0}) = 0$ ( $β_{0} ≃ 3$ ).

With this solution at hand, solving Equation 117 gives

{\tilde{t}}_{s}^{s} (s) = 2 C_{q} (\frac{\sqrt{2} l_{c}}{s} J_{1} (\frac{s}{\sqrt{2} l_{c}}) - \frac{1}{β_{0}} J_{1} (β_{0})),

u (s) = C_{q} [2 (\frac{\sqrt{2} l_{c}}{s} J_{1} (\frac{s}{\sqrt{2} l_{c}}) - \frac{1}{β_{0}} J_{1} (β_{0})) - J_{2} (\frac{s}{\sqrt{2} l_{c}})],

using that $0 = {\tilde{t}}_{s}^{s} (s = s_{t}) ≃ {\tilde{t}}_{s}^{s} (s = R_{d})$ . The constant $C_{q}$ can be determined from Equation 120:

R_{0}^{2} ≃ \int_{0}^{R_{d}} \frac{d s s}{1 + u (s)},

which can be rewritten:

\frac{R_{0}^{2}}{l_{c}^{2}} ≃ 2 \int_{0}^{β_{0}} \frac{d ℓ ℓ}{1 + C_{q} \tilde{u} (ℓ)}

where one has used $u = C_{q} \tilde{u} (\frac{s}{\sqrt{2} l_{c}})$ in the last line, following Equation 124. $\tilde{u} (ℓ)$ is a decreasing function from $ℓ = 0$ to $ℓ = β_{0}$ and $\tilde{u} (β_{0}) < 0$ . In the limit $l_{c} / R_{0} \to 0$ , Equation 126 is satisfied provided that $C_{q} \tilde{u} (β_{0}) \to - 1$ , which sets the constant $C_{q}$ in that limit. Because $\tilde{u} (β_{0}) < 0$ , this implies $C_{q} > 0$ .

In the following, we denote $u_{t} = u (s_{t})$ the deformation reached at the end of the circular plate. Because of the arguments given above, when $l_{c} / R_{0} \to 0$ , $u_{t} \to - 1$ . In general, $C_{q} > 0$ implies that $u_{t} < 0$ : we assume this is the case in the following.

Narrow curved connecting region

In the narrow curved region, at leading order in small $h / R_{d}$ , $\tilde{q} ≃ - u_{t} = | u_{t} |$ and $x = x_{t} ≃ R_{d}$ are homogeneous, and $| (\sin ψ) / x | ≪ | \partial_{s} ψ |$ . The force and torque balance Equations 117 and 118 now give

\partial_{s} (\frac{{\tilde{t}}_{s}^{s}}{\cos ψ}) = - \frac{2 | u_{t} |}{x_{t}},

\frac{κ}{K} \partial_{s}^{2} ψ = {\tilde{t}}_{s}^{s} \tan ψ,

which can be combined to obtain, using $\partial_{s} {\tilde{t}}_{s}^{s} (s_{t}) = 0$ because of the symmetry of the shape:

\frac{κ x_{t}}{2 | u_{t} | K} \partial_{s}^{2} ψ = - s_{c} \sin ψ .

We look for a solution of this equation $ψ (s_{c})$ , with the boundary conditions $ψ (s_{c} = 0) = \frac{π}{2}$ , $ψ (s_{c} \to \infty) = π$ and, using $\partial_{s} z = \sin ψ$ ,

\int_{0}^{\infty} d s_{c} \sin ψ (s_{c}) = h .

It is then helpful to introduce the following differential equation:

\begin{aligned} \partial_{ℓ}^{2} \tilde{ψ} & = - ℓ \sin \tilde{ψ}, \\ \tilde{ψ} (ℓ = 0) & = \frac{π}{2}, \\ \tilde{ψ} (ℓ \to \infty) & = π, \end{aligned}

which admits a solution which can be found numerically. The solution of Equation 129 can then be written

ψ (s_{c}) = \tilde{ψ} (s_{c} {(\frac{2 | u_{t} | K}{κ x_{t}})}^{1 / 3}),

and the constraint Equation 130 gives

h ≃ β_{2} {[\frac{κ R_{d}}{2 | u_{t} | K}]}^{\frac{1}{3}},

where we have introduced $β_{2} = \int_{0}^{\infty} d ℓ \sin \tilde{ψ} (ℓ)$ , with the approximate numerical value $β_{2} ≃ 1.27$ . Together with Equation 122, this analysis indicates that the distance $h$ converges to a constant value as $| ζ_{n} | \to \infty$ (Appendix 4—figure 2). Overall, this analysis predicts the limit value:

\frac{2 h}{R_{d}} ≃ β_{2} {[\frac{2 κ}{| u_{t} | β_{0}^{2} K l_{c}^{2}}]}^{\frac{1}{3}},

such that the condition $h / R_{d} ≪ 1$ implies that the analysis is self-consistent for $κ ≪ K l_{c}^{2}$ .

Appendix 4—figure 1

Download asset Open asset

Schematic of notations used for the asymptotic analysis of a bilayered disc.

Appendix 4—figure 2

Download asset Open asset

Details of shape and nematic profiles for flattened steady-state shapes resulting from a homogeneous nematic tension.

(a) Profile of nematic order parameter q, which decreases for increasing $| ζ_{n} |$ . (b) Distance between the poles of the steady-state solution for different values of $| ζ_{n} |$ and $\tilde{K}$ , and corresponding prediction of Equation 134 (dotted lines). (c) Profile of $ψ (s)$ for different values of $ζ_{n}$ . The profile is invariant with respect to $ζ_{n}$ , for large values of $| ζ_{n} |$ .

Appendix 5

Numerical method for global force balance

In steady-state and dynamical solving, an overall degree of freedom of solid translation of the surface along its axis of symmetry has to be fixed. In steady-state shape calculations, we impose that the south pole is fixed ( $z (s = 0) = 0$ ). In dynamical simulations, we impose that the centroid of the shape does not move (Equation 189). In both cases, these constraints are imposed by introducing a constant dummy force

f^{e x t} = f^{c} e_{z},

which is adjusted to constrain the position of the centre of mass of the shape. The corresponding integrated external force, as defined in Equation 70, is

I (s) = \int_{0}^{s} d s^{'} x f^{c} .

Because at low Reynolds number the total force acting on the surface must vanish, the value of $f^{c}$ should be set to zero by the numerical solver. Indeed, inspection of Equation 72 shows that the conditions ${t_{n}^{s} |}_{s = 0} = {t_{n}^{s} |}_{s = L} = 0$ also imply $I (0) = I (L) = 0$ , and since $f^{c}$ is constant, $f^{c} = 0$ . In practice, $f^{c}$ deviates slightly from zero due to numerical errors.

Appendix 6

Numerical methods to determine steady-state shapes

Stationary shape equations and boundary conditions

The system of stationary shape equations comprises, using the force and torque balance Equations 63–65, the equilibrium Equations 105 and 106 for the nematic order parameter, the definition of the strain Equation 14 and geometrical relations introduced in Appendix 1:

\partial_{s} t_{s}^{s} = 2 \frac{\cos ψ}{x} ζ_{n} q - C_{s}^{s} t_{n}^{s} - f^{c} \sin ψ

\partial_{s} t_{n}^{s} = 2 C_{ϕ}^{ϕ} ζ_{n} q + C_{s}^{s} t_{s}^{s} - \frac{P}{2} - \frac{I}{x^{2}} + f^{c} \cos ψ

\partial_{s} C_{k}^{k} = \frac{1}{2 κ} (t_{n}^{s} - \partial_{s} ζ_{c} + (\partial_{s} ζ_{c n}) q + ζ_{c n} (\partial_{s} q + 2 \frac{\cos ψ}{x} q))

\partial_{s} ψ = C_{s}^{s},

\partial_{s} x = \cos ψ,

\partial_{s} z = \sin ψ,

\partial_{s} s_{0} = \frac{x}{x_{0} (u + 1)},

\partial_{s} q = w,

\partial_{s} w = \frac{1}{2 l_{c}^{2}} q (q^{2} - 1) + \frac{\cos ψ}{x} (4 \frac{\cos ψ}{x} q - w),

\partial_{s} I = x f^{c},

where $C_{ϕ}^{ϕ} = \frac{\sin ψ}{x}$ , $C_{s}^{s} = C_{k}^{k} - C_{ϕ}^{ϕ}$ , $x_{0} = x_{0} (s_{0} (s)) = R_{0} \sin (s_{0} (s) / R_{0})$ , $u = (t_{s}^{s} + ζ_{n} q - ζ) / 2 K$ . Here, Equation 138 has been obtained by using Equation 72 to rewrite the normal force balance (Equation 64) as

\partial_{s} t_{n}^{s} = C_{ϕ}^{ϕ} (t_{ϕ}^{ϕ} - t_{s}^{s}) + C_{s}^{s} t_{s}^{s} - \frac{P}{2} - \frac{I}{x^{2}} + f^{c} \cos ψ,

with the integral of the force $I$ defined in Equation 136. The term $I / x^{2}$ in Equation 147 is regular at the poles by construction since it has the Taylor expansions $I (s) = \frac{f^{c}}{2} s^{2} + O (s^{3})$ at the south pole and $I (L - s) = - \frac{f^{c}}{2} (L - s)^{2} + O ((L - s)^{3})$ at the north pole. Therefore, the limits are $lim_{s \to 0} I / x^{2} = \frac{f^{c}}{2}$ and $lim_{s \to L} I / x^{2} = - \frac{f^{c}}{2}$ .

The geometric boundary conditions for the shape are $ψ (0) = x (0) = x (L) = z (0) = 0$ , $ψ (L) = π$ . Further, we require $t_{n}^{s} (0) = 0$ , $I (0) = I (L) = 0$ , and $s_{0} (0) = 0$ . The equations for the nematic require $q (0) = q (L) = 0$ and $w (0) = w (L) = 0$ , as discussed in Appendix 3. If $l_{a} = L_{0}$ , the unknown length $L$ is determined from the condition $s_{0} (L) = L_{0}$ . If $0 < l_{a} < L_{0}$ , the domain of integration is split into $s \in [0, L_{1}]$ and $s \in [L_{1}, L]$ . The additional condition $s_{0} (L_{1}) = l_{a}$ sets L₁. All unknown functions are matched at the internal boundary $s = L_{1}$ , except for the curvature and strain. These may acquire a jump, $C_{s}^{s} (l_{a}^{+}) - C_{s}^{s} (l_{a}^{-}) = (δ ζ_{c} - δ ζ_{c n} q (l_{a})) / 2 κ$ and $u (l_{a}^{+}) - u (l_{a}^{-}) = (δ ζ - δ ζ_{n} q (l_{a})) / 2 K$ , to ensure continuity in ${\bar{m}}_{s}^{s}$ and $t_{s}^{s}$ .

If applicable, it is convenient to formulate the volume constraint $V = V_{0}$ as a boundary value problem for the partial volume $v (s)$ defined in Equation 51,

\partial_{s} v = π x^{2} \sin ψ, v (0) = 0, v (L) = V_{0} .

Due to the geometric singularities appearing in several of the equations at the poles of the surface $s = 0$ and $s = L$ , we derive the limits of these equations there. Denoting the solution vector by $x (s) = (t_{s}^{s}, t_{n}^{s}, C_{k}^{k}, ψ, x, z, s_{0}, q, w, I, v),$ the limiting expressions of Equations 137–148 at the south and north poles, respectively, are

lim_{s \to 0} \partial_{s} x = (0, \frac{1}{2} (C_{k}^{k} (0) t_{s}^{s} (0) - P + f^{c}), 0, \frac{C_{k}^{k} (0)}{2}, 1, 0, \frac{1}{\sqrt{1 + u (0)}}, 0, W_{0}, 0, 0),

lim_{s \to L} \partial_{s} x = (0, \frac{1}{2} (C_{k}^{k} (L) t_{s}^{s} (L) - P - f^{c}), 0, \frac{C_{k}^{k} (L)}{2}, - 1, 0, \frac{1}{\sqrt{1 + u (L)}}, 0, W_{L}, 0, 0) .

In summary, when considering the full interval and conserved volume, the boundary conditions are

\begin{aligned} at s = 0 : t_{n}^{s} = 0, x = z = ψ = 0, v = 0, s_{0} = 0, q = 0, w = 0, I = 0, \end{aligned}

\begin{array}{ll} at s = L : x = 0, ψ = π, v = V_{0}, s_{0} = L_{0}, q = 0, w = 0, I = 0, \end{array}

and the free parameters are $L, P$ , $f^{c}$ , $W_{0}$ , and $W_{L}$ . Otherwise, if the volume is free to change, then $P = 0$ and the condition $v (L) = V_{0}$ is removed.

In the dimensionless equations, arc lengths are transformed to the unit interval by $\tilde{s} = s / L$ if $l_{a} / L_{0} = 1$ , and to two unit intervals in the case of step-profiles. With $L_{1}$ and $L_{2} = L - L_{1}$ , the dimensionless variables are $\tilde{s} \in [0, 1]$ with $s = \tilde{s} L_{1}$ in the first interval and $\tilde{s} \in [1, 2]$ with $s = L_{1} + L_{2} (\tilde{s} - 1)$ in the second. The boundary value problem given by Equations 137–148 and conditions Equations 151 and 152 in their dimensionless form is solved with the bvp4c solver of MATLAB. The relative and absolute tolerances used in simulations are $t o l_{r e l} = 10^{- 4}$ and $t o l_{a n s} = 10^{- 6}$ , leading to typical adaptive grid sizes of $n g r i d = 100 - 500$ , depending on the shape of the surface. For efficiency, we provide the solver with the analytical Jacobians for the main equations, for the limits at the poles, and for the boundary conditions with respect to the unknowns and the free parameters.

Construction of solution branches

A solution of the mechanical equilibrium equations can be represented by a vector $p \in ℝ^{N}$ , where $ℝ^{N}$ is the vector space spanned by all $N$ (or a subset of) parameters of the model, for example, $P$ and $L$ , the boundary values of the curvature, tension, etc., and the control parameter. A gradual change of the control parameter corresponds to moving along a solution branch in this parameter space. For small increments in the control parameter, the new solution can be obtained numerically using the previous solution as the initial guess for the solver. However, in many cases the solution branch has a fold (e.g. see Figure 2b) and becomes multivalued as a function of the control parameter so that the above method cannot be used.

To continue a solution branch after a fold, we implement a parametric curve approach instead. We denote by $p^{(i)}$ the current state and by $p^{(i - 1)}$ the previous state of the system and approximate by ${\hat{t}}^{(i)} = (p^{(i)} - p^{(i - 1)}) / | p^{(i)} - p^{(i - 1)} |$ the tangent vector at the current state. To find a new solution $p^{(i + 1)}$ on the curve, a step of length $l_{s}$ in the direction of the tangent is taken and the new solution is constrained to lie in a (hyper-)plane perpendicular to the tangent, that is,

((p^{(i + 1)} - p^{(i)}) - l_{s} {\hat{t}}^{(i)}) \cdot {\hat{t}}^{(i)} = 0,

which allows us to eliminate one (e.g. the first) component of p⁽ⁱ⁺¹⁾ via

p_{1}^{(i + 1)} = p_{1}^{(i)} + \frac{l_{s} - (p_{2, \dots, N}^{(i + 1)} - p_{2, \dots, N}^{(i)}) \cdot {\hat{t}}_{2, \dots, N}^{(i)}}{{\hat{t}}_{1}^{(i)}},

provided that ${\hat{t}}_{1}^{(i)} \neq 0$ and where $p_{2, \dots, N}^{(i + 1)}$ denotes the vector $p^{(i + 1)}$ without the first component, etc. The remaining $N - 1$ parameters are determined by the solver such that the new point $p^{(i + 1)}$ is a continuation of the solution branch, which can be achieved for small enough $l_{s}$ . It suffices to include only a subset of the free parameters in this construction, for example, $p = (δ ζ_{c}, P, L_{1})$ for a step-like profile of active moments with conserved volume. Since the elimination (Equation 154) introduces new dependencies of the differential equations and boundary conditions on the free parameters, the analytical Jacobians for the solver are adjusted accordingly.

Appendix 7

Numerical method for the dynamics of active shells

Force and torque balance equations and boundary conditions

The force and torque balance Equations 63–65, together with the constitutive Equations 5–8 and with $f^{e x t} = f^{c} e_{z}$ , can be written as

\begin{array}{ll} \partial_{s}^{2} v^{s} = & - \frac{\cos ψ}{x} (\partial_{s} v^{s} - \frac{\cos ψ}{x} v^{s}) \\ - v_{n} \partial_{s} C_{k}^{k} - \frac{η - η_{b}}{η + η_{b}} C_{s}^{s} C_{ϕ}^{ϕ} v^{s} \\ - \frac{1}{η + η_{b}} (η_{b} C_{k}^{k} + η (C_{s}^{s} - C_{ϕ}^{ϕ})) \partial_{s} v_{n} \\ + \frac{1}{η + η_{b}} (- 2 K \partial_{s} u - \partial_{s} ζ + (\partial_{s} ζ_{n}) q + ζ_{n} (\partial_{s} q) \\ + 2 ζ_{n} q \frac{\cos ψ}{x} - C_{s}^{s} t_{n}^{s} - f^{c} \sin ψ), \end{array}

\partial_{s} t_{n}^{s} = C_{ϕ}^{ϕ} (t_{ϕ}^{ϕ} - t_{s}^{s}) + C_{k}^{k} t_{s}^{s} - \frac{\cos ψ}{x} t_{n}^{s} - P + f^{c} \cos ψ,

\partial_{s} {\bar{m}}_{s}^{s} = 2 ζ_{c n} q \frac{\cos ψ}{x} + t_{n}^{s},

\begin{array}{ll} \partial_{s}^{2} v_{n} = & - \frac{\cos ψ}{x} \partial_{s} v_{n} - ({(C_{ϕ}^{ϕ})}^{2} + {(C_{s}^{s})}^{2}) v_{n} \\ + v^{s} \partial_{s} C_{k}^{k} - \frac{1}{η_{c b}} ({\bar{m}}_{s}^{s} - 2 κ C_{k}^{k} - ζ_{c} + ζ_{c n} q), \end{array}

where in Equation 156 one has to replace

\begin{array}{ll} t_{s}^{s} = & 2 K u + ζ - ζ_{n} q \\ + (η + η_{b}) \partial_{s} v^{s} + (η_{b} - η) \frac{\cos ψ}{x} v^{s} \\ + (η_{b} C_{k}^{k} + η (C_{s}^{s} - C_{ϕ}^{ϕ})) v_{n}, \end{array}

t_{ϕ}^{ϕ} - t_{s}^{s} = 2 (ζ_{n} q + η (\frac{\cos ψ}{x} v^{s} - \partial_{s} v^{s} + (C_{ϕ}^{ϕ} - C_{s}^{s}) v_{n})) .

The solution vector

x (s) = (v^{s}, \partial_{s} v^{s}, t_{n}^{s}, {\bar{m}}_{s}^{s}, v_{n}, \partial_{s} v_{n})

is determined from

\partial_{s} x = L [x],

where the linear operator $ℒ$ is constructed from Equations 155–158 together with two trivial relations relating $\partial_{s} v^{s}$ to $v^{s}$ and $\partial_{s} v^{n}$ to $v^{n}$ . The system of ode’s (Equation 162) is solved on the full interval $[0, L]$ with the boundary conditions

at s = 0 : v^{s} = 0, t_{n}^{s} = 0, \partial_{s} v_{n} = 0,

at s = L : v^{s} = 0, t_{n}^{s} = 0, \partial_{s} v_{n} = 0.

They follow from the requirement that any tangent vector field on the closed axisymmetric surface has to vanish at the poles. Equivalently, these are the conditions required to remove the geometric singularities that appear at the poles in Equations 155–158. We can then derive the well-defined limits

lim_{s \to 0} \partial_{s} x = (\partial_{s} v^{s} (0), 0, \frac{1}{2} (C_{k}^{k} (0) t_{s}^{s} (0) - P + f^{c}), 0, 0, - \frac{1}{4} C_{k}^{k} (0) v_{n} (0) - \frac{{\bar{m}}_{s}^{s} (0) - 2 κ C_{k}^{k} (0) - ζ_{c} (0)}{2 η_{c b}}),

lim_{s \to L} \partial_{s} x = (\partial_{s} v^{s} (L), 0, \frac{1}{2} (C_{k}^{k} (L) t_{s}^{s} (L) - P + f^{c}), 0, 0, - \frac{1}{4} C_{k}^{k} (L) v_{n} (L) - \frac{{\bar{m}}_{s}^{s} (L) - 2 κ C_{k}^{k} (L) - ζ_{c} (L)}{2 η_{c b}}),

where

t_{s}^{s} (0) = 2 K u (0) + ζ (0) + (η + η_{b}) \partial_{s} v^{s} (0) + η_{b} C_{k}^{k} (0) v_{n} (0),

t_{s}^{s} (L) = 2 K u (L) + ζ (L) + (η + η_{b}) \partial_{s} v^{s} (L) + η_{b} C_{k}^{k} (L) v_{n} (L) .

In Equations 165 and 166, we have used that $ζ$ and $u$ , defined as continuous functions on the closed surface, satisfy

\partial_{s} ζ (0) = \partial_{s} ζ (L) = 0, \partial_{s} u (0) = \partial_{s} u (L) = 0.

Prior to solving system Equation 162, at every time step the nematic profiles $q$ and $\partial_{s} q = w$ are determined on the shape $X (ϕ, s, t)$ as solutions of the Euler–Lagrange Equation 16,

\partial_{s} q = w,

\partial_{s} w = \frac{1}{2 l_{c}^{2}} q (q^{2} - 1) \frac{\cos ψ}{x} (4 \frac{\cos ψ}{x} q - w),

with the boundary conditions $q (0) = q (L) = w (0) = w (L) = 0$ , as discussed in Appendix 3.

Constraints

Volume

The volume of a closed surface reads

V = \frac{1}{3} \oint d S X \cdot n .

According to Equations 34, 36, and 38, for a Lagrangian surface update with $\partial_{t} X = v^{i} e_{i} + v_{n} n$ we have

\partial_{t} \sqrt{g} = \sqrt{g} (\nabla_{i} v^{i} + C_{i}^{i} v_{n}),

\partial_{t} n = (- \partial_{i} v_{n} + C_{i j} v^{j}) e^{i} .

This allows to calculate from Equation 172 the rate of change in volume

\partial_{t} V = \frac{1}{3} \oint d s^{1} d s^{2} \partial_{t} (\sqrt{g} X \cdot n)

= \frac{1}{3} \oint d S (v_{n} + (\nabla_{i} v^{i} + C_{i}^{i} v_{n}) X \cdot n + (- \partial_{i} v_{n} + C_{i j} v^{j}) (X \cdot e^{i}))

= \frac{1}{3} \oint d S (v_{n} + C_{i}^{i} v_{n} X \cdot n + C_{i j} v^{j} (X \cdot e^{i}) - v^{k} C_{k}^{i} (X \cdot e_{i}) + v_{n} (2 - C_{i}^{i} X \cdot n))

= \oint d S v_{n},

where we have used two integrations by part and the relations $\partial_{k} (X \cdot n) = X \cdot C_{k}^{i} e_{i}$ , and $\frac{\partial_{i} (\sqrt{g} X \cdot e^{i})}{\sqrt{g}} = 2 - C_{i}^{i} X \cdot n$ .

If the volume is fixed then for all $t$ the dynamics has to satisfy

\partial_{t} V = 0.

On the axisymmetric surface, this results in the integral constraint $\partial_{t} V = 2 π \int d s x v_{n}$ . We define a partial rate of volume change $r_{v} (s) = 2 π \int_{0}^{s} d s^{'} x v_{n}$ , such that the constraint Equation 179 can be written as

\partial_{s} r_{v} = 2 π x v_{n}, r_{v} (0) = r_{v} (L) = 0 .

This is solved simultaneously with the boundary value problem (Equations 162–168), where the pressure $P$ is a free parameter which is required to satisfy both conditions in Equation 180.

On the other hand, if the volume is free to change then pressure is no longer a free parameter, but instead couples the normal force balance (Equation 156) to the rate of change of volume via Equation 21, which can be written:

P = - η_{V} r_{v} (L) .

In this case $P \to 0$ as the dynamics simulation approaches a steady state, in agreement with the direct steady-state calculations in which $P = 0$ .

Rigid-body translation

We note that in the absence of external force, for any flow profile $v (s)$ solution of the force and torque balance Equations 155–158, the flow

\begin{aligned} v^{'} = v + a e_{z} \end{aligned}

\begin{aligned} = v + a \sin ψ (s) e_{s} - a \cos ψ n \end{aligned}

with an arbitrary constant $a$ , is again a solution. The addition of the uniform flow field $a e_{z}$ corresponds to a rigid-body translation in the $z$ -direction which does not affect force balance and therefore makes the task of numerically determining $v$ ill-posed.

To remove this degree of freedom we introduce in each time step a constraint on the translation speed of the centroid. The centroid is equivalent to the centre of mass for a surface with uniform density and is defined as

X_{c} = \frac{\oint d S X (ϕ, s)}{\oint d S} .

For a Lagrangian surface update one obtains

\begin{aligned} \partial_{t} (\oint d S X) & = \oint d S (X \frac{\partial_{t} \sqrt{g}}{\sqrt{g}} + \partial_{t} X) \\ = \oint d S (v_{n} C_{i}^{i} X + v_{n} n), \end{aligned}

\partial_{t} (\oint d S) = \oint d S \frac{\partial_{t} \sqrt{g}}{\sqrt{g}}

= \oint d S v_{n} C_{i}^{i},

where we have used Equation 174 to obtain $\partial_{t} (\sqrt{g}) / \sqrt{g}$ , and the divergence theorem on curved surfaces (Equation 33). It follows that the velocity of the centroid is given by

\begin{aligned} \partial_{t} X_{c} & = \frac{1}{\oint d S} (\partial_{t} (\oint d S X) - X_{c} \partial_{t} (\oint d S)) \\ = \frac{1}{\oint d S} (\oint d S v_{n} (C_{i}^{i} (X - X_{c}) + n)) \end{aligned}

and we want to fix

\partial_{t} X_{c} = 0

for all $t$ . We note that tangential velocity components do not contribute to the centroid displacement.

On the axisymmetric surface $X_{c} = z_{c} e_{z}$ and the constraint Equation 189 becomes

\partial_{t} z_{c} = \frac{1}{\int d s x} (\int d s x v_{n} (C_{i}^{i} (z - z_{c}) - \cos ψ)) = 0.

Analogously to the volume constraint, we introduce a partial centroid velocity $r_{c} (s) = \int_{0}^{s} d s^{'} x v_{n} (C_{i}^{i} (z - z_{c}) - \cos (ψ))$ and add the constraint Equation 190 to the boundary value problem (Equations 162–168) in the form of

\partial_{s} r_{c} = x v_{n} (C_{i}^{i} (z - z_{c}) - \cos ψ), r_{c} (0) = r_{c} (L) = 0.

Note that due to the free parameter $f^{c}$ the number of boundary conditions in the full problem, comprising Equations 162–168 and 191, is consistent.

Time update of the surface

The surface update in each time step $t \to t + δ t$ consists of two sub-steps: first, the material points are advected using a Lagrangian coordinate $s$ as given in Equation 23, then the surface $X^{'} (s, t + δ t)$ is reparametrised in a new arc length coordinate $s^{'}$ , for which $g_{s^{'} s^{'}} = 1$ .

From Equation 22 we find the time updates for the shape descriptors:

\begin{array}{rcl} x (s, t + δ t) & = & x (s, t) + δ t (v_{n} \sin ψ + v^{s} \cos ψ), \\ z (s, t + δ t) & = & z (s, t) + δ t (v^{s} \sin ψ - v_{n} \cos ψ), \\ ψ (s, t + δ t) & = & ψ (s, t) + δ t (- \partial_{s} v_{n} + v^{s} C_{s}^{s}), \end{array}

where the last equation can be shown by taking the time derivative of the normal vector $n$ and using Equation 36. It is convenient to derive the time update for $C_{k}^{k}$ and its derivative from the constitutive Equation 6 and the torque balance (Equation 65),

C_{k}^{k} (s, t + δ t) = C_{k}^{k} (s, t) + δ t \frac{1}{η_{c b}} ({\bar{m}}_{s}^{s} - 2 κ C_{k}^{k} - ζ_{c} + ζ_{c n} q),

\partial_{s} C_{k}^{k} (s, t + δ t) = \partial_{s} C_{k}^{k} (s, t) + δ t \frac{1}{η_{c b}} (t_{n}^{s} + 2 ζ_{c n} q \frac{\cos ψ}{x} + \partial_{s} (- 2 κ C_{k}^{k} - ζ_{c} + ζ_{c n} q))

because the trace of the corotational derivative (Equation 8) is the Lagrangian time evolution of $C_{k}^{k}$ , as can be seen by using Equation 39 with $δ X = δ t v$ . From the variation (Equation 39) the circumferential curvature is updated as

C_{ϕ}^{ϕ} (s, t + δ t) = C_{ϕ}^{ϕ} (s, t) + δ t (- \partial_{s} v_{n} \frac{\cos ψ}{x} - v_{n} {(C_{ϕ}^{ϕ})}^{2} + v^{s} \partial_{s} C_{ϕ}^{ϕ})

and its derivative is obtained from relation Equation 48. In this way only functions which are part of the solution vector (Equation 161) are required for the time updates and we avoid taking numerical gradients of the shape or the velocity fields. Finally, according to Equation 15 we update the area strain as

u (s, t + δ t) = u (s, t) + δ t (1 + u) v_{k}^{k},

with the trace of the strain rate tensor given in Equation 59.

The displaced surface is reparametrised in a new arc length coordinate $s^{'}$ ,

X^{'} (s) = X^{'} (s^{'}),

such that $g_{s^{'} s^{'}} = 1$ , or equivalently $| \partial_{s^{'}} X^{'} (s^{'}) | = 1$ . From

\begin{array}{lcl} \partial_{s^{'}} X^{'} & = & \frac{\partial s}{\partial s^{'}} \partial_{s} X^{'} \\ = & \frac{\partial s}{\partial s^{'}} [\partial_{s} X + δ t \partial_{s} (v^{s} e_{s} + v_{n} n)] \\ = & \frac{\partial s}{\partial s^{'}} [e_{s} (1 + δ t (v_{n} C_{s}^{s} + \partial_{s} v^{s})) + n δ t (\partial_{s} v_{n} - v^{s} C_{s}^{s})] \end{array}

we obtain

{(\frac{\partial s^{'}}{\partial s})}^{2} = 1 + 2 δ t (v_{n} C_{s}^{s} + \partial_{s} v^{s}),

so the relationship between the two arc length parameters is given by the differential equation

\partial_{s} s^{'} = 1 + δ t (v_{n} C_{s}^{s} + \partial_{s} v^{s}) .

We rewrite this equation in terms of a rate of change of arc length $r_{s} = \partial_{t} (s^{'} - s)$ ,

\partial_{s} r_{s} = v_{n} C_{s}^{s} + \partial_{s} v^{s}, r_{s} (0) = 0,

and it is added to the linear system (Equations 162–164). The new arc length is obtained as $s^{'} = s + r_{s} δ t$ and the perimeter length is updated as

L (t + δ t) = L (t) + r_{s} (L (t)) δ t .

For the active profiles defined via sigmoid functions, as given in the main text in Equation 18, the parameters chosen for the simulations are $μ_{t} = 0.01 τ_{a}$ , $σ_{t} = 0.002 τ_{a}$ , $μ_{s} = l_{a}$ , and $σ_{s} = 0.005 L_{0}$ . In the time step $t \to t + δ t$ the bending moment profile, as defined at time $t = 0$ via Equation 19, is updated as

ζ_{c} (s^{'}, t + δ t) = (1 - f (t + δ t, μ_{t}, σ_{t})) (ζ_{c}^{0} + δ ζ_{c} f (s_{0} (s^{'}), l_{a}, σ_{s})),

and an analogous relation holds for $ζ$ , $ζ_{n}$ , and $ζ_{c n}$ . For the spatial dependence, we keep track of the arc length on the initial sphere $s_{0} (s^{'})$ as a function of the arc length on the deformed surface after reparametrisation.

Finally, all surface quantities are saved as spline interpolants on the new arc length $s^{'}$ . For example,

C_{k}^{k} (s^{'}, t + δ t) = C_{k}^{k} (s (s^{'}), t + δ t),

\partial_{s^{'}} C_{k}^{k} (s^{'}, t + δ t) = \frac{\partial s}{\partial s^{'}} \partial_{s} C_{k}^{k} (s (s^{'}), t + δ t),

where the Jacobian prefactor for the derivative is given by Equation 200.

The size of the adaptive time step $δ t$ is determined using a standard step doubling method (Press et al., 2007), where the shape after a full time step, denoted by $X^{(1)} (t + δ t)$ , is compared to that after two half steps, denoted by $X^{(2)} (t + δ t)$ . The relative error is calculated from the shape components and the curvature derivative as

ε_{t} = max {| \frac{x^{(1)} - x^{(2)}}{x^{(1)}} |, | \frac{z^{(1)} - z^{(2)}}{z^{(1)}} |, | \frac{\partial_{s} (C^{(1)})_{k}^{k} - \partial_{s} (C^{(2)})_{k}^{k}}{\partial_{s} (C^{(1)})_{k}^{k}} |}

on a uniform grid which is defined on $X (t)$ and projected onto $X^{(1)}$ and $X^{(2)}$ , whereby the poles and values of $z$ and $\partial_{s} C_{k}^{k}$ , which are too close to zero ( $< 10^{- 3}$ ), are excluded.

Numerical convergence to steady state

In order to validate our simulation method for the dynamics of active surfaces, we analyse the numerical convergence of the dynamics simulations to the steady states obtained from direct calculation (see Appendix 6) for the different active effects. A dynamics simulation is regarded as relaxed to steady state when $max {| {\tilde{v}}_{n} |} < 10^{- 4}$ , which defines $t_{r e l a x}$ . As an example we show in Appendix 7—figure 1 the convergence results for the active bending profile

ζ_{c} (s_{0}) = δ ζ_{c} f (s_{0}, l_{a}, σ_{s})

Appendix 7—figure 1

Download asset Open asset

Convergence analysis of a dynamics simulation to a steady state obtained from direct calculation, for the example shape shown in the inset of (a).

For different $t o l_{t}$ (time step) the error in the shape in (a) and error in the external force integral in (b) are shown.

with $δ ζ_{c} = 80.73 κ / R_{0}$ , $l_{a} = 0.3 L_{0}$ , and conserved volume. This profile results in the folded shape shown in the inset of Appendix 7—figure 1a. The error in the shape (in units of $R_{0}$ ) is calculated as

Δ_{s h a p e} = \frac{1}{N} \sum_{i} \sqrt{{({\tilde{x}}_{i}^{d y n} - {\tilde{x}}_{i}^{s t e a d y})}^{2} + {({\tilde{z}}_{i}^{d y n} - {\tilde{z}}_{i}^{s t e a d y})}^{2}}

on a dimensionless, uniform grid ${\tilde{s}}_{i} \in [0, 1]$ , $i = 1, \dots, N$ , with $N = 1000$ , which is obtained through division by $L$ or $L (t_{r e l a x})$ , respectively, from the corresponding simulation. As discussed in Appendix 5, the deviation of the parameter $f^{c}$ (see Equation 135) from zero characterises the numerical error in the global force balance. As a relative error we plot the value of the parameter $f^{c} (t_{r e l a x})$ normalised by the pressure $P (t_{r e l a x})$ . We find good convergence of the dynamics to the steady-state result with decreasing tolerance $t o l_{t}$ on the time step. Based on these results, we take $t o l_{t} = 10^{- 4}$ for the dynamics simulations shown in the main text. The relative and absolute tolerances for the bvp solver are chosen to be the same as for the direct steady-state calculations: $t o l_{r e l} = 10^{- 4}$ , $t o l_{a b s} = 10^{- 6}$ .

Appendix 8

Isotropic active tension

Consider first a shell with vanishing internal hydrostatic pressure, $P = 0$ , and a step-like tension profile given on the reference surface by

ζ (s_{0}) = {\begin{matrix} ζ^{0} + δ ζ, & if s_{0} \in [0, l_{a}] \\ ζ^{0}, & otherwise . \end{matrix}

One can verify that the spherical shape, given by $x (s) = R \sin (s / R)$ , with $t_{s}^{s} = t_{n}^{s} = 0$ , is a solution of the Equations 137–142. The strain has a jump $δ u = u (l_{a}^{+}) - u (l_{a}^{-}) = δ ζ / (2 K)$ , such that

u (s) = {\begin{cases} - (ζ^{0} + δ ζ) / (2 K), & if s \in [0, l_{a}] \\ - ζ^{0} / (2 K), & otherwise . \end{cases}

Using this to solve Equation 143 for $s \in [0, l_{a}^{'}]$ with $s_{0} (l_{a}^{'}) = l_{a}$ yields

\cos \frac{l_{a}}{R_{0}} - 1 = {(\frac{R}{R_{0}})}^{2} \frac{1}{1 + u (0)} (\cos \frac{l_{a}^{'}}{R} - 1),

and similarly for $s \in [l_{a}^{'}, L]$

\cos \frac{l_{a}}{R_{0}} + 1 = {(\frac{R}{R_{0}})}^{2} \frac{1}{1 + u (L)} (\cos \frac{l_{a}^{'}}{R} + 1) .

These two equations determine the unknown radius $R$ and the deformed active region size $l_{a}^{'}$ as

R = R_{0} \sqrt{\frac{1}{2} \cos \frac{l_{a}}{R_{0}} (u (L) - u (0)) + 1 + \frac{1}{2} (u (0) + u (L))},

l_{a}^{'} = R \arccos (1 + (1 + u (0)) {(\frac{R_{0}}{R})}^{2} (\cos \frac{l_{a}}{R_{0}} - 1)) .

How the ratio $\frac{l_{a}^{'}}{L}$ changes with the tension jump $δ ζ$ is plotted in Appendix 8—figure 1 for $ζ^{0} = 0$ . The above rescaling holds only for $δ ζ < 2 K$ since the active region contracts to a point for $δ ζ \to 2 K$ .

For conserved volume, the spherical shape with $R = R_{0}$ , $t_{n}^{s} = \partial_{s} t_{s}^{s} = 0$ , $P = 2 t_{s}^{s} / R$ is a solution of Equations 137–142. Here, only the relative size of the active patch changes from $l_{a}$ to $l_{a}^{'}$ . As before the strain is piecewise constant with a jump at $l_{a}$ , $δ u = u (L) - u (0) = δ ζ / (2 K)$ , and integrating Equation 143 with $s_{0} (s = l_{a}^{'}) = l_{a}$ , $s_{0} (s = L_{0}) = L_{0}$ gives the additional conditions

\cos \frac{l_{a}}{R_{0}} - 1 = \frac{1}{1 + u (0)} (\cos \frac{l_{a}^{'}}{R_{0}} - 1),

\cos \frac{l_{a}}{R_{0}} + 1 = \frac{1}{1 + u (L)} (\cos \frac{l_{a}^{'}}{R_{0}} + 1),

from which $u (0)$ , $u (L)$ , $l_{a}^{'}$ can be determined.

Appendix 8—figure 1

Download asset Open asset

Rescaling of the size of the active region with active tension difference.

(a) Plot of $l_{a}^{'} / L$ as given by Equations 213; 214. This illustrates the rescaling of the active region size as a function of the tension difference $δ ζ$ , for initial values $l_{a} / L_{0} = 0.1, 0.3, 0.5, 0.7, 0.9$ (blue to purple). (b) Shapes corresponding to two points on the $l_{a} / L_{0} = 0.5$ curve.

Appendix 9

Change of stability at a fold in a solution branch

We discuss here change of stability at a fold in a solution branch (Maddocks, 1987). We consider a dynamical system of the form $\partial_{t} X = F (X, λ)$ , with one control parameter $λ$ . The line of steady state solutions is given by

F (x^{*} (s), λ (s)) = 0

where solutions are parametrised by $s$ . Taking the derivative one obtains ${(\partial_{X} F) |}_{X^{*}} \partial_{s} X^{*} + {(\partial_{λ} F) |}_{X^{*}} \partial_{s} λ = 0$ . We consider a fold in the solution curve where $\partial_{s} λ = 0$ and $| \partial_{s} X^{*} | \neq 0$ . In that case:

(\partial_{x} F) |_{x^{*}} \partial_{s} x^{*} = 0 at the fold,

and the linear stability matrix at the fold ${(\partial_{X} F) |}_{X^{*}}$ thus has (at least) one eigenvalue that changes sign at the fold. Assuming that the solution branch is stable up to the fold, this indicates generically the appearance of an unstable mode; except in the special case where the 0 eigenvalue is a local maximum.

Data availability

The current manuscript is a computational study, so no data have been generated for this manuscript. Modelling code is available at the GitHub repository https://github.com/DianaKhoromskaia/EpithelialShell, (copy archived at swh:1:rev:0838f09f1b2228d8d7da5183fc68f3b49c6ee734).

References

1. Allain J-M
2. Storm C
3. Roux A
4. Ben Amar M
5. Joanny J-F
(2004) Fission of a multiphase membrane tube
Physical Review Letters 93:158104.

https://doi.org/10.1103/PhysRevLett.93.158104
- PubMed
- Google Scholar
(2022) Unifying polar and nematic active matter: emergence and co-existence of half-integer and full-integer topological defects
Journal of Physics A 55:094002.

https://doi.org/10.1088/1751-8121/ac4abe
- Google Scholar
(2014) Bending lipid membranes: experiments after W. helfrich’s model
Advances in Colloid and Interface Science 208:47–57.

https://doi.org/10.1016/j.cis.2014.02.002
- PubMed
- Google Scholar
(2004) Myosin-Dependent junction remodelling controls planar cell intercalation and axis elongation
Nature 429:667–671.

https://doi.org/10.1038/nature02590
- PubMed
- Google Scholar
(2021a) Integer topological defects of cell monolayers: mechanics and flows
Physical Review. E 103:012405.

https://doi.org/10.1103/PhysRevE.103.012405
- PubMed
- Google Scholar
(2021b) Quantifying material properties of cell monolayers by analyzing integer topological defects
Physical Review Letters 126:028101.

https://doi.org/10.1103/PhysRevLett.126.028101
- PubMed
- Google Scholar
1. Braun E
2. Keren K
(2018) Hydra regeneration: closing the loop with mechanical processes in morphogenesis
BioEssays 40:e1700204.

https://doi.org/10.1002/bies.201700204
- PubMed
- Google Scholar
1. Capovilla R
2. Guven J
(2002) Stresses in lipid membranes
Journal of Physics A 35:6233–6247.

https://doi.org/10.1088/0305-4470/35/30/302
- Google Scholar
1. Chen T
2. Callan-Jones A
3. Fedorov E
4. Ravasio A
5. Brugués A
6. Ong HT
7. Toyama Y
8. Low BC
9. Trepat X
10. Shemesh T
11. Voituriez R
12. Ladoux B
(2019) Large-Scale curvature sensing by directional actin flow drives cellular migration mode switching
Nature Physics 15:393–402.

https://doi.org/10.1038/s41567-018-0383-6
- PubMed
- Google Scholar
(1995) How do sea urchins invaginate? using biomechanics to distinguish between mechanisms of primary invagination
Development 121:2005–2018.

https://doi.org/10.1242/dev.121.7.2005
- PubMed
- Google Scholar
Book
1. De Gennes PG
2. Prost J
(1995)
International Series of Monographs on Physics

Oxford, United Kingdom: Oxford University Press.
- Google Scholar
(2017) Topological defects in confined populations of spindle-shaped cells
Nature Physics 13:58–62.

https://doi.org/10.1038/nphys3876
- Google Scholar
(2018) Spontaneous shear flow in confined cellular nematics
Nature Physics 14:728–732.

https://doi.org/10.1038/s41567-018-0099-7
- PubMed
- Google Scholar
(2017) Forces controlling organ growth and size
Mechanisms of Development 144:53–61.

https://doi.org/10.1016/j.mod.2016.11.005
- PubMed
- Google Scholar
1. Etournay R
2. Popović M
3. Merkel M
4. Nandi A
5. Blasse C
6. Aigouy B
7. Brandl H
8. Myers G
9. Salbreux G
10. Jülicher F
11. Eaton S
(2015) Interplay of cell dynamics and epithelial tension during morphogenesis of the Drosophila pupal wing
eLife 4:e07090.

https://doi.org/10.7554/eLife.07090
- Google Scholar
1. Ettensohn CA
(1984) Primary invagination of the vegetal plate during sea urchin gastrulation
American Zoologist 24:571–588.

https://doi.org/10.1093/icb/24.3.571
- Google Scholar
1. Fouchard J
2. Wyatt TP
3. Proag A
4. Lisica A
5. Khalilgharibi N
6. Recho P
7. Suzanne M
8. Kabla A
9. Charras G
(2020) Curling of epithelial monolayers reveals coupling between active bending and tissue tension
PNAS 117:9377–9383.

https://doi.org/10.1073/pnas.1917838117
- PubMed
- Google Scholar
1. Fourcade B
2. Miao L
3. Rao M
4. Wortis M
5. Zia RKP
(1994) Scaling analysis of narrow necks in curvature models of fluid lipid-bilayer vesicles
Physical Review E 49:5276–5286.

https://doi.org/10.1103/PhysRevE.49.5276
- Google Scholar
1. Frank JR
2. Kardar M
(2008) Defects in nematic membranes can buckle into pseudospheres
Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics 77:041705.

https://doi.org/10.1103/PhysRevE.77.041705
- PubMed
- Google Scholar
Book
1. Gilbert MJF
2. Barresi SF
(2020)
Developmental Biology

Oxford, United Kingdom: Oxford University Press.
- Google Scholar
1. Gokhale RH
2. Shingleton AW
(2015) Size control: the developmental physiology of body and organ size regulation
Wiley Interdisciplinary Reviews. Developmental Biology 4:335–356.

https://doi.org/10.1002/wdev.181
- PubMed
- Google Scholar
(2010) Aspiration of biological viscoelastic drops
Physical Review Letters 104:218101.

https://doi.org/10.1103/PhysRevLett.104.218101
- PubMed
- Google Scholar
(2018) The noisy basis of morphogenesis: mechanisms and mechanics of cell sheet folding inferred from developmental variability
PLOS Biology 16:e2005536.

https://doi.org/10.1371/journal.pbio.2005536
- PubMed
- Google Scholar
1. Helfrich W
2. Prost J
(1988) Intrinsic bending force in anisotropic membranes made of chiral molecules
Physical Review. A, General Physics 38:3065–3068.

https://doi.org/10.1103/physreva.38.3065
- PubMed
- Google Scholar
Preprint
(2021) Defect-Mediated Morphogenesis
arXiv.

https://arxiv.org/abs/2105.15200
- Google Scholar
(2015) Dynamics of a Volvox embryo turning itself inside out
Physical Review Letters 114:178101.

https://doi.org/10.1103/PhysRevLett.114.178101
- PubMed
- Google Scholar
1. Hopf H
(1927) Abbildungsklassenn-dimensionaler mannigfaltigkeiten
Mathematische Annalen 96:209–224.

https://doi.org/10.1007/BF01209163
- Google Scholar
(2017) The hope and the hype of organoid research
Development 144:938–941.

https://doi.org/10.1242/dev.150201
- PubMed
- Google Scholar
1. Ishihara K
2. Tanaka EM
(2018) Spontaneous symmetry breaking and pattern formation of organoids
Current Opinion in Systems Biology 11:123–128.

https://doi.org/10.1016/j.coisb.2018.06.002
- Google Scholar
(2007) Vesicle shape, molecular tilt, and the suppression of necks
Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics 76:031908.

https://doi.org/10.1103/PhysRevE.76.031908
- PubMed
- Google Scholar
1. Jülicher F
2. Lipowsky R
(1993) Domain-Induced budding of vesicles
Physical Review Letters 70:2964–2967.

https://doi.org/10.1103/PhysRevLett.70.2964
- Google Scholar
1. Jülicher F
2. Lipowsky R
(1996) Shape transformations of vesicles with intramembrane domains
Physical Review E 53:2670–2683.

https://doi.org/10.1103/PhysRevE.53.2670
- Google Scholar
1. Kamm RD
2. Bashir R
3. Arora N
4. Dar RD
5. Gillette MU
6. Griffith LG
7. Kemp ML
8. Kinlaw K
9. Levin M
10. Martin AC
11. McDevitt TC
12. Nerem RM
13. Powers MJ
14. Saif TA
15. Sharpe J
16. Takayama S
17. Takeuchi S
18. Weiss R
19. Ye K
20. Yevick HG
21. Zaman MH
(2018) Perspective: the promise of multi-cellular engineered living systems
APL Bioengineering 2:040901.

https://doi.org/10.1063/1.5038337
- PubMed
- Google Scholar
(2017) Topological defects control collective dynamics in neural progenitor cell cultures
Nature 545:327–331.

https://doi.org/10.1038/nature22321
- Google Scholar
1. Keber FC
2. Loiseau E
3. Sanchez T
4. DeCamp SJ
5. Giomi L
6. Bowick MJ
7. Marchetti MC
8. Dogic Z
9. Bausch AR
(2014) Topology and dynamics of active nematic vesicles
Science 345:1135–1139.

https://doi.org/10.1126/science.1254784
- PubMed
- Google Scholar
1. Khalilgharibi N
2. Mao Y
(2021) To form and function: on the role of basement membrane mechanics in tissue development, homeostasis and disease
Open Biology 11:200360.

https://doi.org/10.1098/rsob.200360
- PubMed
- Google Scholar
1. Kierzenka J
2. Shampine LF
(2001) A BVP solver based on residual control and the maltab PSE
ACM Transactions on Mathematical Software 27:299–316.

https://doi.org/10.1145/502800.502801
- Google Scholar
1. Knoche S
2. Kierfeld J
(2011) Buckling of spherical capsules
Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics 84:046608.

https://doi.org/10.1103/PhysRevE.84.046608
- PubMed
- Google Scholar
1. Kominami T
2. Takata H
(2004) Gastrulation in the sea urchin embryo: a model system for analyzing the morphogenesis of a monolayered epithelium
Development, Growth & Differentiation 46:309–326.

https://doi.org/10.1111/j.1440-169x.2004.00755.x
- PubMed
- Google Scholar
1. Kralj S
2. Rosso R
3. Virga EG
(2011) Curvature control of valence on nematic shells
Soft Matter 7:670–683.

https://doi.org/10.1039/C0SM00378F
- Google Scholar
1. Latorre E
2. Kale S
3. Casares L
4. Gómez-González M
5. Uroz M
6. Valon L
7. Nair RV
8. Garreta E
9. Montserrat N
10. Del Campo A
11. Ladoux B
12. Arroyo M
13. Trepat X
(2018) Active superelasticity in three-dimensional epithelia of controlled shape
Nature 563:203–208.

https://doi.org/10.1038/s41586-018-0671-4
- PubMed
- Google Scholar
1. Lei M
2. Schumacher LJ
3. Lai YC
4. Juan WT
5. Yeh CY
6. Wu P
7. Jiang TX
8. Baker RE
9. Widelitz RB
10. Yang L
11. Chuong CM
(2017) Self-organization process in newborn skin organoid formation inspires strategy to restore hair regeneration of adult cells
PNAS 114:E7101–E7110.

https://doi.org/10.1073/pnas.1700475114
- PubMed
- Google Scholar
(2017) Structural inheritance of the actin cytoskeletal organization determines the body axis in regenerating Hydra
Cell Reports 18:1410–1421.

https://doi.org/10.1016/j.celrep.2017.01.036
- PubMed
- Google Scholar
1. Lomholt MA
(2006) Mechanics of nonplanar membranes with force-dipole activity
Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics 73:061913.

https://doi.org/10.1103/PhysRevE.73.061913
- PubMed
- Google Scholar
1. Lubensky TC
2. Prost J
(1992) Orientational order and vesicle shape
Journal de Physique II 2:371–382.

https://doi.org/10.1051/jp2:1992133
- Google Scholar
1. MacKintosh FC
2. Lubensky TC
(1991) Orientational order, topology, and vesicle shapes
Physical Review Letters 67:1169–1172.

https://doi.org/10.1103/PhysRevLett.67.1169
- PubMed
- Google Scholar
1. Maddocks JH
(1987) Stability and folds
Archive for Rational Mechanics and Analysis 99:301–328.

https://doi.org/10.1007/BF00282049
- Google Scholar
(2014) Activating membranes
Physical Review Letters 112:258101.

https://doi.org/10.1103/PhysRevLett.112.258101
- PubMed
- Google Scholar
1. Marchetti MC
2. Joanny J-F
3. Ramaswamy S
4. Liverpool TB
5. Prost J
6. Rao M
7. Simha RA
(2013) Hydrodynamics of soft active matter
Reviews of Modern Physics 85:1143–1189.

https://doi.org/10.1103/RevModPhys.85.1143
- Google Scholar
1. Marmottant P
2. Mgharbel A
3. Käfer J
4. Audren B
5. Rieu JP
6. Vial JC
7. van der Sanden B
8. Marée AFM
9. Graner F
10. Delanoë-Ayari H
(2009) The role of fluctuations and stress on the effective viscosity of cell aggregates
PNAS 106:17271–17275.

https://doi.org/10.1073/pnas.0902085106
- PubMed
- Google Scholar
(2021) Topological defects in the nematic order of actin fibres as organization centres of Hydra morphogenesis
Nature Physics 17:251–259.

https://doi.org/10.1038/s41567-020-01083-1
- Google Scholar
1. Martin AC
2. Goldstein B
(2014) Apical constriction: themes and variations on a cellular mechanism driving morphogenesis
Development 141:1987–1998.

https://doi.org/10.1242/dev.102228
- PubMed
- Google Scholar
1. Martin AC
(2020) Self-organized cytoskeletal alignment during Drosophila mesoderm invagination
Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 375:20190551.

https://doi.org/10.1098/rstb.2019.0551
- PubMed
- Google Scholar
1. Meinhardt A
2. Eberle D
3. Tazaki A
4. Ranga A
5. Niesche M
6. Wilsch-Bräuninger M
7. Stec A
8. Schackert G
9. Lutolf M
10. Tanaka EM
(2014) 3D reconstitution of the patterned neural tube from embryonic stem cells
Stem Cell Reports 3:987–999.

https://doi.org/10.1016/j.stemcr.2014.09.020
- PubMed
- Google Scholar
1. Messal HA
2. Alt S
3. Ferreira RM
4. Gribben C
5. Wang VM-Y
6. Cotoi CG
7. Salbreux G
8. Behrens A
(2019) Tissue curvature and apicobasal mechanical tension imbalance instruct cancer morphogenesis
Nature 566:126–130.

https://doi.org/10.1038/s41586-019-0891-2
- PubMed
- Google Scholar
(2019) Topology and morphology of self-deforming active shells
Physical Review Letters 123:208001.

https://doi.org/10.1103/PhysRevLett.123.208001
- PubMed
- Google Scholar
(2019a) Minimal model of cellular symmetry breaking
Physical Review Letters 123:188101.

https://doi.org/10.1103/PhysRevLett.123.188101
- PubMed
- Google Scholar
(2019b) Self-organized shape dynamics of active surfaces
PNAS 116:29–34.

https://doi.org/10.1073/pnas.1810896115
- PubMed
- Google Scholar
1. Mongera A
2. Rowghanian P
3. Gustafson HJ
4. Shelton E
5. Kealhofer DA
6. Carn EK
7. Serwane F
8. Lucio AA
9. Giammona J
10. Campàs O
(2018) A fluid-to-solid jamming transition underlies vertebrate body axis elongation
Nature 561:401–405.

https://doi.org/10.1038/s41586-018-0479-2
- PubMed
- Google Scholar
1. Morris RG
2. Rao M
(2019) Active morphogenesis of epithelial monolayers
Physical Review 100:022413.

https://doi.org/10.1103/PhysRevE.100.022413
- PubMed
- Google Scholar
(2019) Emergence of active nematic behavior in monolayers of isotropic cells
Physical Review Letters 122:048004.

https://doi.org/10.1103/PhysRevLett.122.048004
- PubMed
- Google Scholar
1. Napoli G
2. Vergori L
(2012) Surface free energies for nematic shells
Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics 85:061701.

https://doi.org/10.1103/PhysRevE.85.061701
- PubMed
- Google Scholar
(1992) N -atic order and continuous shape changes of deformable surfaces of genus zero
Europhysics Letters 20:279–284.

https://doi.org/10.1209/0295-5075/20/3/015
- Google Scholar
(2019) Geometrical control of active turbulence in curved topographies
Physical Review Letters 122:168002.

https://doi.org/10.1103/PhysRevLett.122.168002
- PubMed
- Google Scholar
1. Popović M
2. Nandi A
3. Merkel M
4. Etournay R
5. Eaton S
6. Jülicher F
7. Salbreux G
(2017) Active dynamics of tissue shear flow
New Journal of Physics 19:033006.

https://doi.org/10.1088/1367-2630/aa5756
- Google Scholar
Book
(2007)
Numerical Recipes 3rd Edition: The Art of Scientific Computing

Cambridge University Press.
- Google Scholar
1. Ranft J
2. Basan M
3. Elgeti J
4. Joanny JF
5. Prost J
6. Jülicher F
(2010) Fluidization of tissues by cell division and apoptosis
PNAS 107:20863–20868.

https://doi.org/10.1073/pnas.1011086107
- PubMed
- Google Scholar
Book
1. Reddy JN
(2006) Theory and Analysis of Elastic Plates and Shells
Florida, United States: CRC Press.

https://doi.org/10.1201/9780849384165
- Google Scholar
(2018) Progress and potential in organoid research
Nature Reviews. Genetics 19:671–687.

https://doi.org/10.1038/s41576-018-0051-9
- PubMed
- Google Scholar
(2017) Irreversible thermodynamics of curved lipid membranes
Physical Review. E 96:042409.

https://doi.org/10.1103/PhysRevE.96.042409
- PubMed
- Google Scholar
(2009) Hydrodynamics of cellular cortical flows and the formation of contractile rings
Physical Review Letters 103:058102.

https://doi.org/10.1103/PhysRevLett.103.058102
- PubMed
- Google Scholar
1. Salbreux G
2. Jülicher F
(2017) Mechanics of active surfaces
Physical Review. E 96:032404.

https://doi.org/10.1103/PhysRevE.96.032404
- PubMed
- Google Scholar
(2022) Theory of nematic and polar active fluid surfaces
Physical Review Research 4:033158.

https://doi.org/10.1103/PhysRevResearch.4.033158
- Google Scholar
1. Saw TB
2. Doostmohammadi A
3. Nier V
4. Kocgozlu L
5. Thampi S
6. Toyama Y
7. Marcq P
8. Lim CT
9. Yeomans JM
10. Ladoux B
(2017) Topological defects in epithelia govern cell death and extrusion
Nature 544:212–216.

https://doi.org/10.1038/nature21718
- PubMed
- Google Scholar
(1991) Shape transformations of vesicles: phase diagram for spontaneous- curvature and bilayer-coupling models
Physical Review. A, Atomic, Molecular, and Optical Physics 44:1182–1202.

https://doi.org/10.1103/physreva.44.1182
- PubMed
- Google Scholar
1. Seifert U
(1997) Configurations of fluid membranes and vesicles
Advances in Physics 46:13–137.

https://doi.org/10.1080/00018739700101488
- Google Scholar
1. Selinger JV
2. Schnur JM
(1993) Theory of chiral lipid tubules
Physical Review Letters 71:4091–4094.

https://doi.org/10.1103/PhysRevLett.71.4091
- Google Scholar
(1996) Theory of cylindrical tubules and helical ribbons of chiral lipid membranes
Physical Review E 53:3804–3818.

https://doi.org/10.1103/PhysRevE.53.3804
- Google Scholar
1. Sens P
2. Turner MS
(2004) Theoretical model for the formation of caveolae and similar membrane invaginations
Biophysical Journal 86:2049–2057.

https://doi.org/10.1016/S0006-3495(04)74266-6
- PubMed
- Google Scholar
1. Serra D
2. Mayr U
3. Boni A
4. Lukonin I
5. Rempfler M
6. Challet Meylan L
7. Stadler MB
8. Strnad P
9. Papasaikas P
10. Vischi D
11. Waldt A
12. Roma G
13. Liberali P
(2019) Self-Organization and symmetry breaking in intestinal organoid development
Nature 569:66–72.

https://doi.org/10.1038/s41586-019-1146-y
- PubMed
- Google Scholar
1. Staykova M
2. Arroyo M
3. Rahimi M
4. Stone HA
(2013) Confined bilayers passively regulate shape and stress
Physical Review Letters 110:028101.

https://doi.org/10.1103/PhysRevLett.110.028101
- PubMed
- Google Scholar
1. Sui L
2. Alt S
3. Weigert M
4. Dye N
5. Eaton S
6. Jug F
7. Myers EW
8. Jülicher F
9. Salbreux G
10. Dahmann C
(2018) Differential lateral and basal tension drive folding of Drosophila wing discs through two distinct mechanisms
Nature Communications 9:4620.

https://doi.org/10.1038/s41467-018-06497-3
- PubMed
- Google Scholar
(2018) Morphogenesis and compartmentalization of the intestinal crypt
Developmental Cell 45:183–197.

https://doi.org/10.1016/j.devcel.2018.03.024
- PubMed
- Google Scholar
1. Tu ZC
2. Seifert U
(2007) Concise theory of chiral lipid membranes
Physical Review E 76:603.

https://doi.org/10.1103/PhysRevE.76.031603
- Google Scholar
1. Turlier H
2. Audoly B
3. Prost J
4. Joanny J-F
(2014) Furrow constriction in animal cell cytokinesis
Biophysical Journal 106:114–123.

https://doi.org/10.1016/j.bpj.2013.11.014
- PubMed
- Google Scholar
1. Yang Q
2. Xue S-L
3. Chan CJ
4. Rempfler M
5. Vischi D
6. Maurer-Gutierrez F
7. Hiiragi T
8. Hannezo E
9. Liberali P
(2021) Cell fate coordinates mechano-osmotic forces in intestinal crypt formation
Nature Cell Biology 23:733–744.

https://doi.org/10.1038/s41556-021-00700-2
- PubMed
- Google Scholar

Article and author information

Author details

Diana Khoromskaia

The Francis Crick Institute, London, United Kingdom

Contribution
Conceptualization, Software, Investigation, Visualization, Methodology, Writing – original draft

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-2597-6336
Guillaume Salbreux
1. The Francis Crick Institute, London, United Kingdom
2. University of Geneva, Geneva, Switzerland
Contribution
Conceptualization, Supervision, Funding acquisition, Investigation, Methodology, Writing – original draft

For correspondence
Guillaume.Salbreux@unige.ch

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-7041-1292

Funding

Cancer Research UK (C55977/A23342)

Guillaume Salbreux

Cancer Research UK (FC001317)

Guillaume Salbreux

Medical Research Council (FC001317)

Guillaume Salbreux

Wellcome Trust (FC001317)

Guillaume Salbreux

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. For the purpose of Open Access, the authors have applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission.

Acknowledgements

We thank S Grigolon for useful discussions and N Cuny for comments on the manuscript. DK and GS acknowledge support from the Francis Crick Institute, which receives its core funding from Cancer Research UK (FC001317), the UK Medical Research Council (FC001317), and the Wellcome Trust (FC001317), and from the CRUK multidisciplinary project award C55977/A23342.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

1,829

views
290

downloads
15

citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Diana Khoromskaia
Guillaume Salbreux

(2023)

Active morphogenesis of patterned epithelial shells

eLife 12:e75878.

https://doi.org/10.7554/eLife.75878

Share this article

Cite this article

A two-dimensional surface with nematic order represents an epithelial sheet undergoing active deformations.

Deformations of epithelial shells due to active bending moments, with free (a–d) and conserved (e–i) volume.

Nematic order on a sphere.

Deformations of epithelial shells due to nematic tensions, with free (a–e) and conserved (f–i) volume.

Deformations of epithelial shells due to nematic bending moments, with free (a–c) and conserved (d, e) volume.

Summary of shape changes obtained through patterning of isotropic and anisotropic active tensions and bending moments.

Schematic of the surface S1 used to derive the integral of the normal force balance.

Schematic of notations used for the asymptotic analysis of a bilayered disc.

Details of shape and nematic profiles for flattened steady-state shapes resulting from a homogeneous nematic tension.

Convergence analysis of a dynamics simulation to a steady state obtained from direct calculation, for the example shape shown in the inset of (a).

Rescaling of the size of the active region with active tension difference.

Author details

Diana Khoromskaia

Contribution

Competing interests

Guillaume Salbreux

Contribution

For correspondence

Competing interests

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Categories and tags

Further reading

Schematic of the surface $S_{1}$ used to derive the integral of the normal force balance.