Active morphogenesis of patterned epithelial shells

  1. Diana Khoromskaia
  2. Guillaume Salbreux  Is a corresponding author
  1. The Francis Crick Institute, United Kingdom
  2. University of Geneva, Switzerland


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.


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
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=m¯ϕϕx and msϕ=-m¯ss/x. (c) Stresses integrated across the thickness of the sheet result in tensions tij and bending moments mij 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 tij and m¯ij, which can be captured by a nematic order parameter Qij (e.g. blue rods on the top surface).


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 tij, but also the bending moment tensor mij 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 m¯ij=-mikϵkj, which we adopt for convenience. The local force balance projected on the tangential and normal directions reads (Salbreux and Jülicher, 2017)

(1) itij+Cijtni=fext,j
(2) itniCijtij=fnextP,

where notations of differential geometry are introduced in Appendix 1; briefly Cij is the curvature tensor, gij denotes the metric tensor, and ϵij the antisymmetric Levi-Civita tensor, n the vector normal to the surface, tij is the tangential contribution of the tension tensor and tni its normal contribution, and 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:

(3) tni=km¯ki,
(4) ϵijtij=Cijmij.

We assume an external force density fext=fnextn+fext,jej acting on the surface in addition to a difference of hydrostatic (uniform) pressure P=Pin-Pout, 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, SdSfext=0. Dissipative couplings to the external fluid are ignored here as the characteristic viscosity of a biological tissue (105 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:

(5) tsij=(2Ku+ζ+(ηbη)vkk)gij+2ηvij+ζnQij,
(6) m¯ij=(2κCkk+ζc+ηcbDDtCkk)gij+ζcnQij.

where tsij is the symmetric part of the tension tensor and, on a curved surface, the strain rate tensor vij and the corotational time derivative of the curvature tensor DDtCij are given by (Salbreux and Jülicher, 2017)

(7) vij=12(ivj+jvi)+Cijvn,
(8) DDtCij=i(jvn)vnCikCkj+vkkCij+ωn(ϵikCkj+ϵjkCki),

with ωn=12ϵijivj 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. Qij 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 ηcb. 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 Ckk. 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, ζcn 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 Ckk=ζ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 ζch(γa-γb)/2 emerges in the surface to lowest order in the curvature tensor. The term in ζcn 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 Qija and Qijb or in the level of nematic active stress ζna and ζnb. For example, such differences could stem from two contractile (respectively extensile) layers with perpendicular nematic orientations +Qij and -Qij (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 ζcn(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 Qij of the tensor tsij and m¯ij (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 Cij. 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):

(9) F=dS(k2(iQjk)(iQjk)a4QijQij+a16(QijQij)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 QijQij=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)3 is parametrised by the arc length coordinate s[0,L] and the angle of rotation ϕ[0,2π] as

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

The local tangent basis is given by {eϕ,es}, 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 gss=1, which implies relations between the tangent angle ψ(s)[0,π] and the shape functions x(s) and z(s)

(11) sx=cosψ,
(12) sz=sinψ,

which, together with the meridional principal curvature

(13) Css=sψ,

are sufficient to reconstruct the surface shape from the curvature Css. In this axisymmetric setup, the velocity field reads v=vses+vnn, with vs the tangential and vn the normal velocities.

The undeformed initial surface is a sphere S0 with radius R0, and all quantities defined on it are denoted with a subscript ‘0’. We define the area strain on a point of the surface as

(14) u=dSdS0dS0,

where dS is the surface area element at the point considered on the surface, and dS0 is the surface area element of the same material point on the sphere. With this definition, u=0 on the initial sphere. We denote s0(s) the arc length position on the undeformed sphere S0 of a material point at arc length position s on the deformed sphere. One then has u=fϕfs-1 with fs=dsds0 the meridional stretch and fϕ=xx0 the circumferential stretch. Integrating fs1=fϕ/(u+1) yields the arc length reparametrisation s0(s) between the initial and the deformed surface. The Lagrangian time derivative of the area strain (Equation 14) is related to the flow through

(15) DDtu=(1+u)vkk.

Nematic order

Here, with axial symmetry, the nematic tensor Qij has the non-zero component q=Qϕϕ=-Qss. 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 lc=k/a, which is the characteristic nematic length. In this geometry the Euler–Lagrange equation resulting from the free energy (Equation 9) is

(16) s2q=12lc2q(q21)+cosψx(4cosψxqsq).

An example solution of Equation 16 on the sphere is shown in Figure 3b. From the two possible states with q=±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 Qij to represent the order parameter on the outer side of the shell, the sign convention is such that ζcn>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 ζnQij and the active moment ζcnQij, 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 laL0 (Figure 1b), such that the active terms are given on S0 by step-like profiles, for example

(17) ζc(s0)={ζc0+δζc,if s0[0,la]ζc0,otherwise

and similarly for ζ(s0), ζn(s0), and ζcn(s0). 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=ζc0=ζn0=ζcn0=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

(18) f(x,μ,σ)=1(1+exμσ)1,

for their space and time dependence. For instance, the active bending moment profile is defined on S0 as

(19) ζc(s0,t=0)=(1f(t=0,μt,σt))(ζc0+δζcf(s0,la,σs))

as a smooth version of the step-profile Equation 17, and ζ, ζn, and ζcn 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
Deformations of epithelial shells due to active bending moments, with free (a–d) and conserved (e–i) volume.

(a, e) Shape diagram. (b, f) Details of shape diagram illustrating different behaviours of solution branches. The ideal neck line (green) represents the bending moment difference required to create budded shapes consisting of two spheres with u=0, as given by Equation 24. (c) Examples of solution branches in the (δζc,V)-plane corresponding to four different regions in (b). (g) Examples of solution branches in the (δζc,P)-plane chosen from three different regions in (e). (d, h) Dynamic simulations of shape changes, for parameter values indicated in the shape diagrams (a, e). (i) Neck radius and curvatures at the neck as functions of δζc for the example la/L0=0.04 in (g). Other parameters: K~=103,η~cb=10-2, η~V=10-4.


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:

(20) V=V0,

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

(21) P=ηVtV

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 vs=vn=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 {vs,vn} 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 Qij(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,

(22) tX=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 κ/R02, bending moment densities by κ/R0, lengths by R0, force densities by κ/R03, viscosities by the two-dimensional shear viscosity η of the epithelium, times by the characteristic time scale τa=ηR02/κ , and velocities by R0/τa. This leaves the dimensionless parameters K~=KR02/κ, l~c=lc/R0, η~b=ηb/η, η~cb=ηcbR02/η and η~V=ηVR04/η to be fixed. We choose to set η~=η~b=1, η~V=10-4 for fast relaxation of the volume, and the nematic length scale is set to 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 K~=12(R0/h)2. In simulations we use K~=1000, corresponding to h/R00.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 η~cb(h/R0)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 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 t~0.02. The time integration according to Equation 22 is done with an explicit Euler method with adaptive step size via

(23) X(ϕ,s,t+δt)=X(ϕ,s,t)+δtv(ϕ,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 gss=1. The profiles and surface quantities are passed between time steps as spline interpolants.

To produce the diagrams of steady-state shapes, la 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.


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 la/L0 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 la. 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 (C0)kk=-ζc2κ, 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 Css 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 la/L0 (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 la/L0 is small or close to 1, the resulting shape is characterised by the formation of a bud which form either inwards (laL0) or outwards (L0-laL0). 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 10-3R0; extrapolation indicates full constriction at a finite δζc (Figure 2i). As the neck radius decreases the principal curvatures at the neck diverge as Css,Cϕϕ±, such that Ckk 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/R1+1/R2=C0 with R1 and R2 the radius of the two spheres and C0 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 tss changes sign across the neck, while m¯ss 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 Ra, Rp, related by the condition:

(24) 1Ra1Rp=δζc4κ,1Ra1Rp=δζc4κ,

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 R1 and R2 as a function of la/L0. Combining these conditions determine a curve in the parameter space δζcR0/κ, la/L0, 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 (laL0) and close to L0 (L0-laL0) values of la (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 la/L0, a fold in the solution branch in the (δζc,V)-plane appears (Figure 2c). For most values of la/L0, 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, la/L0=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, la/L0=0.31,0.7). Alternatively the solution branch continues in a sequence of loops and the active region elongates (Figure 2c, la/L0=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, la/L0=0.1,δζc=40κ/R0). 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, la/L0=0.35,δζc=12.5κ/R0). 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 (la/L0=0.85,δζc=15κ/R0), or from a sphere to a torus via self-intersection (la/L0=0.35,δζc=12.5κ/R0).

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 la/L0, the shell deforms into locally folded shapes, which eventually self-intersect at large bending moment difference (Figure 2g, la/L0=0.3,0.7, Figure 2h).

Nematic active tensions

We now introduce the nematic order parameter Qij 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, Qss=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 Qij-Qij, ζn-ζ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.

Nematic order on a sphere.

(a) Two possible configurations for the nematic order parameter Qij 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 lc=0.1R0). (b) Order parameter q(s)=Qϕϕ(s) as a solution of the Euler–Lagrange Equation 16 on a sphere with R0=1 and lc=0.1R0;q=1 at the equator and q=0 at the locations of the defects (poles). For uniform ζn, ζniQis 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, la/L0=1 (Figure 4a–e).

Figure 4 with 13 supplements see all
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: K~=103,η~cb=10-2, η~V=10-4, 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, la/L0=1, ζn<0). There is no shape transition occurring; instead, we find that the shape converges to a limit shape as |ζn| (Appendix 4). The limit shape corresponds to two parallel flat discs of radius Rd, separated by a distance 2h, 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 κKlc2:

(25) Rdlc,h(κlcK)13

The first relation shows that the limit shape has the size of the characteristic nematic length lc. Physically, for lcL0, 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, la/L0=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, q0 as a result of the diverging principal curvatures (as can be seen from the presence of a term (cos(ψ)xq)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<la/L0<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 (Ckk0, Figure 4d). For 0<la/L0<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 la 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 lc (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 Qij. 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 ζcn.

We first discuss the case where the nematic active bending moment is homogeneous (la/L0=1), where there is no difference of pressure across the surface, and where ζcn=δζcn<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 |ζcn|. To characterise this, we note that the corresponding steady-state shape solutions have vanishing tensions tss=0 and tns=0 everywhere (Figure 5—figure supplement 1) and the force balances Equations 63 and 64 are trivially satisfied. The torque balance Equation 65 reads

(26) 2κsCkkζcnsq=2ζcncosψxq.

Combining Equations 26 and 48 one obtains that L[CssCϕϕζcnq/(2κ)]=0, with the operator =s+2cosψx. Solutions to [f]=0 have the form f=A/x2 with A a constant. The boundary condition that the function f should be finite at the poles requires A=0, such that

(27) CssCϕϕ=qζcn2κ .

As a result, if the shape has a cylindrical part, in which Css=0 and q=1, then the cylinder radius Rc is given by

(28) 1Rc=ζcn2κ,

and since such solutions are area-preserving, with u=0, the length of the cylindrical part scales as Lc1/Rc. These relations are in excellent agreement with simulation results for large enough |ζcn| (Figure 5g).

Figure 5 with 1 supplement see all
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 la/L0=1,0.7,0.3, δζcn<0 with analytical predictions. The tube length is measured on the steady-state shape as the arc length of the deformed active region, stube=s(s0=la), and the tube curvature as Cϕϕ(stube/2). Other parameters: K~=1000,η~cb=10-2, η~V=10-4, l~c=0.1. In (c), (f), for δζcn,ζcn<0 the orientation of the director field drawn on the surface (black lines) is set by -Qij.

When la<L0, the active region forms an outward cylindrical protrusion (Figure 5a, b and g) whose radius is still well described by Equation 28, replacing ζcn by δζcn, 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

(29) LcR02Rc(1coslaR0)=δζcnR022κ(1coslaR0),

which is again in excellent agreement with numerical simulation for large |δζcn| and for different values of la/L0 (Figure 5g).

For ζcn>0 and la/L0=1 we find erythrocyte-like shapes, where the indentations at the poles become stronger with ζcn until the two poles touch (Figure 5b). This behaviour remains for la<L0, 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 la/L0<1 and the volume is free to change, both signs of δζcn result in a cylindrical appendage forming from the active region. The sign of δζcn 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 δζcn<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 δζcn<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 ζcn<0 and la=L0 the surface becomes spindle-like, narrowing at the poles with increasing |ζcn|. As in the free volume case, when ζcn>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 la<L0. Shape transitions occur at large |δζcn|, for both δζcn<0 and δζcn>0 (Figure 5e). In the case δζcn<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 δζcn>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 la/L0. 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 (la/L0=0.3,δζcn=-150κ/R0) or at the poles (la/L0=0.5,δζcn=50κ/R0).


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.

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 ζcn<0 the director field orientation (black lines) is set by -Qij.

Some of our findings recapitulate epithelial deformations observed in biological systems. The flattening observed for an extensile homogeneous nematic surface (Figure 4b, la/L0=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 ζcn on the trace or determinant of the curvature tensor Cij. 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(s1,s2) be a curved surface embedded in 3 and parametrised by the generalised coordinates s1,s2. A local covariant basis in the tangent plane is given by the vectors ei=iX=X/si and the unit normal vector is constructed as n=(e1×e2)/|e1×e2| and chosen to point outwards for a closed surface in our convention. These define the metric tensor gij=eiej and the curvature tensor Cij=-(iej)n=eijn. The infinitesimal surface and line elements are given by dS=gds1ds2 and dl2=gijdsidsj, where g is the determinant of the metric tensor. The antisymmetric Levi–Civita tensor is defined as ϵij=n(ei×ej).

Covariant derivatives

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

(30) Γijk=12gkm(jgim+igjmmgij),

and Christoffel symbols of the first kind are defined as Γkij=12(jgki+igkj-kgij). The covariant derivatives of a tangent vector field fi and a tangent tensor field Tij, respectively, are given by

(31) ifj=ifj+Γikjfk,
(32) iTjk=iTjk+ΓiljTlk+ΓilkTjl

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

(33) SdSifi=Cdlνifi.
Infinitesimal variation of surface quantities

For a small deformation of the surface

(34) δX=δXiei+δXnn ,

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)

(35) δei=(iδXj+CijδXn)ej+(iδXnCijδXj)n,
(36) δn=(iδXn+CijδXj)ei,
(37) δgij=iδXj+jδXi+2CijδXn,
(38) δgg=12gijδgij=kδXk+CkkδXn,
(39) δCij=i(jδXn)+(iδXk)Ckj(kδXj)Cik+(kCij)δXkδXnCikCkj.

These relations can be used to obtain time derivatives of surface quantities using δX=δtv (Lagrangian surface update) or δX=δtvnn (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 (s0,ϕ) on the undeformed and deformed surfaces and denoting g0 the determinant of the metric of the undeformed surface,

(40) u=gg0g0,

we obtain from Equation 38 the variation:

(41) δu=(1+u)(kδXk+CkkδXn),

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

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

(42) eϕ=(xsinϕxcosϕ0),
(43) es=(cosϕcosψsinϕcosψsinψ),
(44) n=(cosϕsinψsinϕsinψcosψ),

where we have used sx=cosψ and sz=sinψ, which can be defined through the requirement that s is an arc length parameter, such that |es|2=(sx)2+(sz)2=1. The metric and curvature tensors and the surface element are given by

(45) gij=(x2001),Cij=(sinψx00sψ),dS=xdsdϕ.

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 Tss=Tss=Tss and Tϕϕ=Tϕϕ=Tϕϕ. The circumferential and meridional principal curvatures Cϕϕ and Css, and the mean and Gaussian curvatures H and K are given by:

(46) Cϕϕ=sinψx,Css=sψ,
(47) H=12Ckk,K=detCij=CϕϕCss.

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

(48) sCϕϕ=cosψx(CssCϕϕ),
(49) s(cosψx)=CϕϕCss(cosψx)2.

The partial area and partial volume are given by

(50) a(s)=2π0sdsx(s),
(51) v(s)=π0sdsx(s)2sinψ(s) .

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

(52) DCkkDt=s2vncosψxsvnvn((Css)2+(Cϕϕ)2)+vssCkk.

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 Tsϕ=Tϕs=0. The only non-vanishing component of igjk is sgϕϕ=2xcosψ; therefore, the non-zero Christoffel symbols of the first kind are

(53) Γϕsϕ=Γϕϕs=Γsϕϕ=12sgϕϕ=xcosψ,

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

(54) Γϕϕs=xcosψ,Γsϕϕ=Γϕsϕ=cosψx.

The non-zero components of ifj are

(55) sfs=sfs,
(56) ϕfϕ=cosψxfs,

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

(57) vss=svs+Cssvn,
(58) vϕϕ=cosψxvs+Cϕϕvn,

with trace

(59) vkk=svs+cosψxvs+Ckkvn.

The non-zero components of iTjk are

(60) sTss=sTss,
(61) sTϕϕ=sTϕϕ+2cosψxTϕϕ,
(62) ϕTsϕ=xcosψTϕϕ+cosψxTss=cosψx(TssTϕϕ)=ϕTϕs ,

where Tϕϕ=gϕϕTϕϕ=1x2Tϕϕ 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

(63) stss=cosψx(tϕϕtss)Csstnsfsext,
(64) stns=Cϕϕ(tϕϕtss)+CkktsscosψxtnsfnextP,
(65) sm¯ss=cosψx(m¯ϕϕm¯ss)+tns.

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, ϵijtij=Cijmij=ζcnQikϵkCijj which vanishes for an axisymmetric surface where Qsϕ=Qϕs=0. Therefore, here tij=tsij 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 S1 bounded by the south pole and a circle C perpendicular to the axis of symmetry, given by s=s1 (Appendix 1—figure 1). The bounding circle C has the line element dl=xdϕ and the unit normal v=es, tangent to the surface and pointing outward with respect to S1. The balance of forces acting on S1 reads

(66) Cdlνiti+S1dS(fext+Pn)=0.

The contributions to the integral Equation 66 are

(67) Cdlνiti=02πdϕxts=02πdϕx(tsses+tnsn)=2πx(tsssinψtnscosψ)ez,
(68) PS1dSn=(2πP0s1dsxcosψ)ez=(2πP0x(s1)dxx)ez=πPx2ez,
(69) S1dSfext=2π(0s1dsxfzext)ez=2πI(s1)ez,

where we have introduced the integrated external force

(70) I(s)=0sdsxfzext,

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

(71) 2πx(tsssinψtnscosψ)πPx2+2πI=0.

From Equation 71 one obtains an expression for the transverse tension for ψπ2:

(72) tns=tsstanψ12xcosψP+1xcosψ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
Schematic of the surface S1 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:

(73) x(s)=x(0)+cosψ|s=0s12((sinψ)sψ)|s=0s2+O(s3)=s+O(s3),
(74) x(Ls)=x(0)+cosψ|s=L(Ls)12((sinψ)sψ)|s=L(Ls)2+O((Ls)3)=sL+O((Ls)3).

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

(75) lims0,Lf(s)x(s)=lims0,Lsf|s=0,Lcosψ(s)|s=0,L=±sf|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

(76) lims0,LCϕϕ=Css|s=0,L.

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

(77) sCkk|s=0,L=0.

From Equation 48 we find that at the poles 2sCϕϕ=sCss, which, together with Equation 77 yields

(78) sCϕϕ|s=0,L=sCss|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 R0 and south pole at the origin, X(ϕ,0)=0:

(79) L0=πR0,0sπR0,
(80) x(s0)=R0sins0R0,
(81) z(s0)=R0(1coss0R0),
(82) ψ(s0)=s0R0,
(83) Css=Cϕϕ=H=1R0,
(84) cosψx=1R0tans0R0,
(85) a(s0)=A02(1coss0R0),
(86) v(s0)=V0(2+coss0R0)sin4s02R0,

where the arc length is denoted s0, and A0=4πR02, V0=(4/3)πR03 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=ζcn=0). In that case, the tensors tij and m¯ij 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 ψπ2:

(87) s(tsscosψ)=xsψ2cos2ψP,
(88) sm¯ss=tsstanψ12xcosψ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 ψ=π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 s¯=s/a, x¯=x/a and C¯kk=aCkk. The force balance equations then become for ψπ2:

(89) s¯(tsscosψ)=ax¯s¯ψ2cos2ψP ,
(90) s¯(2κC¯kk)=as¯ζc+a2tsstanψa32x¯cosψP .

For a0, the leading order solution has tss/cosψ and C¯kk both constant. Using the relation

(91) C¯kk=s¯ψ+sinψx¯=1x¯x¯(x¯2C¯ϕϕ),

with C¯ϕϕ=aCϕϕ, and the conditions C¯ϕϕ(x¯=1)=1 and that sinψ does not diverge for |x¯|, the curvatures have solution C¯kk=0 and

(92) Cϕϕ=Css=ax2,cosψ=±1a2x2,

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±. Since tss/cosψ is constant across the infinitesimal neck, tss also changes sign asymptotically away from the neck.

The next order in a of Equation 90 gives s¯(2κCkk+ζc)=0 which corresponds to m¯ss 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 tss=u=0 and constant m¯ss (Figure 2—figure supplement 1a). Considering a shape with the active and passive regions forming spheres with radii Ra and Rp separated by an infinitesimal neck, the condition of constant m¯ss results in

(93) 2κ[Ckk]a+δζc=2κ[Ckk]p,

with [Ckk]a and [Ckk]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 [Ckk]a=±2/Ra 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:

(94) 4πRa2=2πR02(1coslaR0),4πRp2=2πR02(1+coslaR0),

where we have used Equation 85. When [Ckk]a=2/Ra corresponding to the active region towards the outside, and δζc>0, Equation 24 implies that Rp<Ra which further requires la/L0>1/2 to satisfy Equation 94.

Combining Equations 93 and 94 gives a condition defining a curve in the parameter space la/L0, δζcR0/κ:

(95) δζcR04κ=21+cos(πlaL0)21cos(πlaL0),

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 n^=n^. In an orthonormal frame {e^ϕ,e^s}, where e^ϕ=eϕ/x and e^s=es, it is characterised by an angle α[0,π] as

(96) n^=cosαe^ϕ+sinαe^s.

The director components in the basis {eϕ,es} are then

(97) n^ϕ=cosαx,n^s=sinα.

The traceless and symmetric nematic order parameter Qij can be constructed from the director ni and a magnitude S as

(98) Qij=S(n^in^j12gij).

Its components read

(99) Qϕϕ=S2cos2αx2,Qss=S2cos2α,Qsϕ=Qϕs=S2sin2αx.

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 ζn0, ζcn0:

(100) itiϕ=(ζn+ζcnsinψx)(sQsϕ+3cosψxQsϕ)=0,

which requires Qsϕ=Qϕs=0 excluding divergence of Qsϕ 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

(101) q=S2cos2α=Qϕϕ=Qss,

which can take the values q=±S/2, corresponding to azimuthal or longitudinal orientation of the director n^, respectively.

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

(102) F=dSf=dS(k((sq)2+4(cosψxq)2)a2q2+a4q4).

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

(103) 0=δFδq=hiΠi,


(104) h=fq,Πi=f(iq).

Here, we have used the definition from the functional derivative, dFdSδFδqdq. Equation 16 can be written as two first-order equations

(105) sq=w,
(106) sw=12lc2q(q21)+cosψx(4cosψxqw).

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

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

and also implies

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

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

(109) sw(0)=12lc2q(0)(q(0)21)+lims0[cos(ψ)x(4cos(ψ)xqw)]=12lc2q(0)(q(0)21)+lims0[1s(41s(q(0)+w(0)s+12sw(0)s2)(w(0)+sw(0)s))]=lims0[1s(3w(0)+sw(0)s)]=sw(0),

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

(110) (sq,sw)s=0=(0,W0),
(111) (sq,sw)s=L=(0,WL),

with two free parameters W0 and WL, 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 R0 to axisymmetric perturbations; a more general analysis can be found in Napoli and Vergori, 2012. We note that for a spherical shape, with θ=s/R0 and at first order in q:

(112) δFδq2kR02[θ2+cotθθ4cot2θ+R022lc2]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 =θ2+cotθθ-4cot2θ is provided by taking derivatives of axisymmetric spherical harmonics, qn(θ)=qnfn(θ) with fn(θ)=[θ2-cotθθ]Pn(cosθ) with Pn the Legendre polynomial of degree n, for n2. One then finds

(113) δFδq[qn]2kR02[n(n+1)4R022lc2]qnfn(θ)

The isotropic state is stable if n(n+1)4R022lc20 for n2, or for

(114) lcRo>12.

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=ζcn=0.

We postulate that the limit shape reached as |ζn| consists of two parallel flat central discs of radius Rd, separated by a distance 2h, and connected by a narrow curved region (Appendix 4—figure 1). One denotes st the arclength of the extreme position of the shape where x=xt is maximal, and sc=s-st 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 sq(s=st)=0. One looks for an asymptotic shape which satisfies hRd; we show later that this condition requires κKlc2.

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

(115) s[tsscosψ]=2ζnqx,
(116) 2κs[sψ+sinψx]=tsstanψ.

with tss=2Ku-ζnq. The last equation implies that at ψ=π/2 (i.e. at the point of the surface with the extremal value of x), tss=0 and therefore at this point u=ζnq/(2K). However, our definition of deformation implies that u>1. Therefore, as |ζn|, one must have q0 (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 q~=-ζn/(2K)q and renormalised tension t~ss=tss/(2K)=u+q~, the force and torque balance equation and the linearised equilibrium equation for the order parameter read

(117) s[t~sscosψ]=2q~x,
(118) κKs[sψ+sinψx]=t~sstanψ,
(119) s2q~+cosψxsq~4cos2ψx2q~=12lc2q~,

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

(120) 2R02=0Ldsx1+u.

In the asymptotic regime where hRd, 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

(121) s2q~+1ssq~4s2q~=12lc2q~,

with solution, using q~(s=0)=0, q~(x)=CqJ2(x2lc), with Cq a constant to determine, and Jn(x) are the Bessel functions of the first kind. The condition 0=sq(s=st)sq(s=Rd) then yields the expression for the radius of the disc:

(122) Rd=2lcβ0,

with β0 the smallest positive solution of J2(β0)=0 (β03).

With this solution at hand, solving Equation 117 gives

(123) t~ss(s)=2Cq(2lcsJ1(s2lc)1β0J1(β0)),
(124) u(s)=Cq[2(2lcsJ1(s2lc)1β0J1(β0))J2(s2lc)],

using that 0=t~ss(s=st)t~ss(s=Rd). The constant Cq can be determined from Equation 120:

(125) R020Rddss1+u(s),

which can be rewritten:

(126) R02lc220β0d1+Cqu~()

where one has used u=Cqu~(s2lc) in the last line, following Equation 124. u~() is a decreasing function from =0 to =β0 and u~(β0)<0. In the limit lc/R00, Equation 126 is satisfied provided that Cqu~(β0)-1, which sets the constant Cq in that limit. Because u~(β0)<0, this implies Cq>0.

In the following, we denote ut=u(st) the deformation reached at the end of the circular plate. Because of the arguments given above, when lc/R00, ut-1. In general, Cq>0 implies that ut<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/Rd, q~-ut=|ut| and x=xtRd are homogeneous, and |(sinψ)/x||sψ|. The force and torque balance Equations 117 and 118 now give

(127) s(t~sscosψ)=2|ut|xt,
(128) κKs2ψ=t~sstanψ,

which can be combined to obtain, using st~ss(st)=0 because of the symmetry of the shape:

(129) κxt2|ut|Ks2ψ=scsinψ.

We look for a solution of this equation ψ(sc), with the boundary conditions ψ(sc=0)=π2, ψ(sc)=π and, using sz=sinψ,

(130) 0dscsinψ(sc)=h.

It is then helpful to introduce the following differential equation:

(131) 2ψ~=sinψ~,ψ~(=0)=π2,ψ~()=π,

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

(132) ψ(sc)=ψ~(sc(2|ut|Kκxt)1/3),

and the constraint Equation 130 gives

(133) hβ2[κRd2|ut|K]13,

where we have introduced β2=0dsinψ~(), with the approximate numerical value β21.27. Together with Equation 122, this analysis indicates that the distance h converges to a constant value as |ζn| (Appendix 4—figure 2). Overall, this analysis predicts the limit value:

(134) 2hRdβ2[2κ|ut|β02Klc2]13,

such that the condition h/Rd1 implies that the analysis is self-consistent for κKlc2.

Appendix 4—figure 1
Schematic of notations used for the asymptotic analysis of a bilayered disc.
Appendix 4—figure 2
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 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

(135) fext=fcez,

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

(136) I(s)=0sdsxfc.

Because at low Reynolds number the total force acting on the surface must vanish, the value of fc should be set to zero by the numerical solver. Indeed, inspection of Equation 72 shows that the conditions tns|s=0=tns|s=L=0 also imply I(0)=I(L)=0, and since fc is constant, fc=0. In practice, fc 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:

(137) stss=2cosψxζnqCsstnsfcsinψ
(138) stns=2Cϕϕζnq+CsstssP2Ix2+fccosψ
(139) sCkk=12κ(tnssζc+(sζcn)q+ζcn(sq+2cosψxq))
(140) sψ=Css,
(141) sx=cosψ,
(142) sz=sinψ,
(143) ss0=xx0(u+1),
(144) sq=w,
(145) sw=12lc2q(q21)+cosψx(4cosψxqw),
(146) sI=xfc,

where Cϕϕ=sinψx, Css=CkkCϕϕ, x0=x0(s0(s))=R0sin(s0(s)/R0), u=(tss+ζnq-ζ)/2K. Here, Equation 138 has been obtained by using Equation 72 to rewrite the normal force balance (Equation 64) as

(147) stns=Cϕϕ(tϕϕtss)+CsstssP2Ix2+fccosψ,

with the integral of the force I defined in Equation 136. The term I/x2 in Equation 147 is regular at the poles by construction since it has the Taylor expansions I(s)=fc2s2+O(s3) at the south pole and I(Ls)=fc2(Ls)2+O((Ls)3) at the north pole. Therefore, the limits are lims0I/x2=fc2 and limsLI/x2=fc2.

The geometric boundary conditions for the shape are ψ(0)=x(0)=x(L)=z(0)=0, ψ(L)=π. Further, we require tns(0)=0, I(0)=I(L)=0, and s0(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 la=L0, the unknown length L is determined from the condition s0(L)=L0. If 0<la<L0, the domain of integration is split into s[0,L1] and s[L1,L]. The additional condition s0(L1)=la sets L1. All unknown functions are matched at the internal boundary s=L1, except for the curvature and strain. These may acquire a jump, Css(la+)Css(la)=(δζcδζcnq(la))/2κ and u(la+)-u(la-)=(δζ-δζnq(la))/2K, to ensure continuity in m¯ss and tss.

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

(148) sv=πx2sinψ,v(0)=0,v(L)=V0.

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)=(tss,tns,Ckk,ψ,x,z,s0,q,w,I,v), the limiting expressions of Equations 137–148 at the south and north poles, respectively, are

(149) lims0sx=(0,12(Ckk(0)tss(0)P+fc),0,Ckk(0)2,1,0,11+u(0),0,W0,0,0),
(150) limsLsx=(0,12(Ckk(L)tss(L)Pfc),0,Ckk(L)2,1,0,11+u(L),0,WL,0,0).

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

(151) at s=0:tns=0,x=z=ψ=0,v=0,s0=0,q=0,w=0,I=0,
(152) at s=L:x=0,ψ=π,v=V0,s0=L0,q=0,w=0,I=0,

and the free parameters are L,P, fc, W0, and WL. Otherwise, if the volume is free to change, then P=0 and the condition v(L)=V0 is removed.

In the dimensionless equations, arc lengths are transformed to the unit interval by s~=s/L if la/L0=1, and to two unit intervals in the case of step-profiles. With L1 and L2=L-L1, the dimensionless variables are s~[0,1] with s=s~L1 in the first interval and s~[1,2] with s=L1+L2(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 tolrel=10-4 and tolans=10-6, leading to typical adaptive grid sizes of ngrid=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 pN, 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 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 ls 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,

(153) ((p(i+1)p(i))lst^(i))t^(i)=0,

which allows us to eliminate one (e.g. the first) component of p(i+1) via

(154) p1(i+1)=p1(i)+ls(p2,,N(i+1)p2,,N(i))t^2,,N(i)t^1(i),

provided that t^1(i)0 and where p2,,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 ls. It suffices to include only a subset of the free parameters in this construction, for example, p=(δζc,P,L1) 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 fext=fcez, can be written as

(155) s2vs=cosψx(svscosψxvs)vnsCkkηηbη+ηbCssCϕϕvs1η+ηb(ηbCkk+η(CssCϕϕ))svn+1η+ηb(2Ksusζ+(sζn)q+ζn(sq)+2ζnqcosψxCsstnsfcsinψ),
(156) stns=Cϕϕ(tϕϕtss)+CkktsscosψxtnsP+fccosψ,
(157) sm¯ss=2ζcnqcosψx+tns ,
(158) s2vn=cosψxsvn((Cϕϕ)2+(Css)2)vn+vssCkk1ηcb(m¯ss2κCkkζc+ζcnq),

where in Equation 156 one has to replace

(159) tss=2Ku+ζζnq+(η+ηb)svs+(ηbη)cosψxvs+(ηbCkk+η(CssCϕϕ))vn,
(160) tϕϕtss=2(ζnq+η(cosψxvssvs+(CϕϕCss)vn)).

The solution vector

(161) x(s)=(vs,svs,tns,m¯ss,vn,svn)

is determined from

(162) sx=L[x],

where the linear operator is constructed from Equations 155–158 together with two trivial relations relating svs to vs and svn to vn. The system of ode’s (Equation 162) is solved on the full interval [0,L] with the boundary conditions

(163) at s=0:vs=0,tns=0,svn=0,
(164) at s=L:vs=0,tns=0,svn=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

(165) lims0sx=(svs(0),0,12(Ckk(0)tss(0)P+fc),0,0,14Ckk(0)vn(0)m¯ss(0)2κCkk(0)ζc(0)2ηcb),
(166) limsLsx=(svs(L),0,12(Ckk(L)tss(L)P+fc),0,0,14Ckk(L)vn(L)m¯ss(L)2κCkk(L)ζc(L)2ηcb),


(167) tss(0)=2Ku(0)+ζ(0)+(η+ηb)svs(0)+ηbCkk(0)vn(0),
(168) tss(L)=2Ku(L)+ζ(L)+(η+ηb)svs(L)+ηbCkk(L)vn(L).

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

(169) sζ(0)=sζ(L)=0,su(0)=su(L)=0.

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

(170) sq=w,
(171) sw=12lc2q(q21)cosψx(4cosψxqw),

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



The volume of a closed surface reads

(172) V=13dSXn.

According to Equations 34, 36, and 38, for a Lagrangian surface update with tX=viei+vnn we have

(173) tg=g(ivi+Ciivn),
(174) tn=(ivn+Cijvj)ei.

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

(175) tV=13ds1ds2t(gXn)
(176) =13dS(vn+(ivi+Ciivn)Xn+(ivn+Cijvj)(Xei))
(177) =13dS(vn+CiivnXn+Cijvj(Xei)vkCki(Xei)+vn(2CiiXn))
(178) =dSvn,

where we have used two integrations by part and the relations k(Xn)=XCkiei, and i(gXei)g=2CiiXn.

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

(179) tV=0.

On the axisymmetric surface, this results in the integral constraint tV=2πdsxvn. We define a partial rate of volume change rv(s)=2π0sdsxvn, such that the constraint Equation 179 can be written as

(180) srv=2πxvn,rv(0)=rv(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:

(181) P=ηVrv(L).

In this case P0 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

(182) v=v+aez
(183) =v+asinψ(s)esacosψn

with an arbitrary constant a, is again a solution. The addition of the uniform flow field aez 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

(184) Xc=dSX(ϕ,s)dS.

For a Lagrangian surface update one obtains

(185) t(dSX)=dS(Xtgg+tX)=dS(vnCiiX+vnn),
(186) t(dS)=dStgg
(187) =dSvnCii,

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

(188) tXc=1dS(t(dSX)Xct(dS))=1dS(dSvn(Cii(XXc)+n))

and we want to fix

(189) tXc=0

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

On the axisymmetric surface Xc=zcez and the constraint Equation 189 becomes

(190) tzc=1dsx(dsxvn(Cii(zzc)cosψ))=0.

Analogously to the volume constraint, we introduce a partial centroid velocity rc(s)=0sdsxvn(Cii(zzc)cos(ψ)) and add the constraint Equation 190 to the boundary value problem (Equations 162–168) in the form of

(191) src=xvn(Cii(zzc)cosψ),rc(0)=rc(L)=0.

Note that due to the free parameter fc 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 tt+δ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 gss=1.

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

(192) x(s,t+δt)=x(s,t)+δt(vnsinψ+vscosψ),z(s,t+δt)=z(s,t)+δt(vssinψvncosψ),ψ(s,t+δt)=ψ(s,t)+δt(svn+vsCss),

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 Ckk and its derivative from the constitutive Equation 6 and the torque balance (Equation 65),

(193) Ckk(s,t+δt)=Ckk(s,t)+δt1ηcb(m¯ss2κCkkζc+ζcnq),
(194) sCkk(s,t+δt)=sCkk(s,t)+δt1ηcb(tns+2ζcnqcosψx+s(2κCkkζc+ζcnq))

because the trace of the corotational derivative (Equation 8) is the Lagrangian time evolution of Ckk, as can be seen by using Equation 39 with δX=δtv. From the variation (Equation 39) the circumferential curvature is updated as

(195) Cϕϕ(s,t+δt)=Cϕϕ(s,t)+δt(svncosψxvn(Cϕϕ)2+vssCϕϕ)

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

(196) u(s,t+δt)=u(s,t)+δt(1+u)vkk,

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

The displaced surface is reparametrised in a new arc length coordinate s,

(197) X(s)=X(s),

such that gss=1, or equivalently |sX(s)|=1. From

(198) sX=sssX=ss[sX+δts(vses+vnn)]=ss[es(1+δt(vnCss+svs))+nδt(svnvsCss)]

we obtain

(199) (ss)2=1+2δt(vnCss+svs),

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

(200) ss=1+δt(vnCss+svs).

We rewrite this equation in terms of a rate of change of arc length rs=t(s-s),

(201) srs=vnCss+svs,rs(0)=0,

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

(202) L(t+δt)=L(t)+rs(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=la, and σs=0.005L0. In the time step tt+δt the bending moment profile, as defined at time t=0 via Equation 19, is updated as

(203) ζc(s,t+δt)=(1f(t+δt,μt,σt))(ζc0+δζcf(s0(s),la,σs)),

and an analogous relation holds for ζ, ζn, and ζcn. For the spatial dependence, we keep track of the arc length on the initial sphere s0(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,

(204) Ckk(s,t+δt)=Ckk(s(s),t+δt),
(205) sCkk(s,t+δt)=sssCkk(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

(206) εt=max{|x(1)x(2)x(1)|,|z(1)z(2)z(1)|,|s(C(1))kks(C(2))kks(C(1))kk|}

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 sCkk, which are too close to zero (<103), 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{|v~n|}<104, which defines trelax. As an example we show in Appendix 7—figure 1 the convergence results for the active bending profile

(207) ζc(s0)=δζcf(s0,la,σs)
Appendix 7—figure 1
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 tolt (time step) the error in the shape in (a) and error in the external force integral in (b) are shown.

with δζc=80.73κ/R0, la=0.3L0, 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 R0) is calculated as

(208) Δshape=1Ni(x~idynx~isteady)2+(z~idynz~isteady)2

on a dimensionless, uniform grid s~i[0,1], i=1,,N, with N=1000, which is obtained through division by L or L(trelax), 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(trelax) normalised by the pressure P(trelax). We find good convergence of the dynamics to the steady-state result with decreasing tolerance tolt on the time step. Based on these results, we take tolt=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: tolrel=10-4, tolabs=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

(209) ζ(s0)={ζ0+δζ,ifs0[0,la]ζ0,otherwise.

One can verify that the spherical shape, given by x(s)=Rsin(s/R), with tss=tns=0, is a solution of the Equations 137–142. The strain has a jump δu=u(la+)u(la)=δζ/(2K), such that

(210) u(s)={(ζ0+δζ)/(2K),if s[0,la]ζ0/(2K),otherwise.

Using this to solve Equation 143 for s[0,la] with s0(la)=la yields

(211) coslaR01=(RR0)211+u(0)(coslaR1),

and similarly for s[la,L]

(212) coslaR0+1=(RR0)211+u(L)(coslaR+1).

These two equations determine the unknown radius R and the deformed active region size la as

(213) R=R012coslaR0(u(L)u(0))+1+12(u(0)+u(L)),
(214) la=Rarccos(1+(1+u(0))(R0R)2(coslaR01)).

How the ratio laL changes with the tension jump δζ is plotted in Appendix 8—figure 1 for ζ0=0. The above rescaling holds only for δζ<2K since the active region contracts to a point for δζ2K.

For conserved volume, the spherical shape with R=R0, tns=stss=0, P=2tss/R is a solution of Equations 137–142. Here, only the relative size of the active patch changes from la to la. As before the strain is piecewise constant with a jump at la, δu=u(L)u(0)=δζ/(2K), and integrating Equation 143 with s0(s=la)=la, s0(s=L0)=L0 gives the additional conditions

(215) coslaR01=11+u(0)(coslaR01),
(216) coslaR0+1=11+u(L)(coslaR0+1),

from which u(0), u(L), la can be determined.

Appendix 8—figure 1
Rescaling of the size of the active region with active tension difference.

(a) Plot of la/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 la/L0=0.1,0.3,0.5,0.7,0.9 (blue to purple). (b) Shapes corresponding to two points on the la/L0=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 tX=F(X,λ), with one control parameter λ. The line of steady state solutions is given by

(217) F(x(s),λ(s))=0

where solutions are parametrised by s. Taking the derivative one obtains (XF)|X*sX*+(λF)|X*sλ=0. We consider a fold in the solution curve where sλ=0 and |sX*|0. In that case:

(218) (xF)|xsx=0at the fold,

and the linear stability matrix at the fold (XF)|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, (copy archived at swh:1:rev:0838f09f1b2228d8d7da5183fc68f3b49c6ee734).


  1. Book
    1. De Gennes PG
    2. Prost J
    International Series of Monographs on Physics
    Oxford, United Kingdom: Oxford University Press.
  2. Book
    1. Gilbert MJF
    2. Barresi SF
    Developmental Biology
    Oxford, United Kingdom: Oxford University Press.
    1. Knoche S
    2. Kierfeld J
    (2011) Buckling of spherical capsules
    Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics 84:046608.
    1. Lomholt MA
    (2006) Mechanics of nonplanar membranes with force-dipole activity
    Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics 73:061913.
    1. Maddocks JH
    (1987) Stability and folds
    Archive for Rational Mechanics and Analysis 99:301–328.
    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.
    1. Napoli G
    2. Vergori L
    (2012) Surface free energies for nematic shells
    Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics 85:061701.
  3. Book
    1. Press WH
    2. Teukolsky SA
    3. Vetterling WT
    4. Flannery BP
    Numerical Recipes 3rd Edition: The Art of Scientific Computing
    Cambridge University Press.

Article and author information

Author details

  1. Diana Khoromskaia

    The Francis Crick Institute, London, United Kingdom
    Conceptualization, Software, Investigation, Visualization, Methodology, Writing – original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2597-6336
  2. Guillaume Salbreux

    1. The Francis Crick Institute, London, United Kingdom
    2. University of Geneva, Geneva, Switzerland
    Conceptualization, Supervision, Funding acquisition, Investigation, Methodology, Writing – original draft
    For correspondence
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7041-1292


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.


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.

Version history

  1. Preprint posted: November 24, 2021 (view preprint)
  2. Received: November 25, 2021
  3. Accepted: November 18, 2022
  4. Version of Record published: January 17, 2023 (version 1)


© 2023, Khoromskaia and Salbreux

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.


  • 1,520
  • 252
  • 11

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)

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

  1. Diana Khoromskaia
  2. Guillaume Salbreux
Active morphogenesis of patterned epithelial shells
eLife 12:e75878.

Share this article

Further reading

    1. Developmental Biology
    2. Physics of Living Systems
    Raphaël Clément

    Geometric criteria can be used to assess whether cell intercalation is active or passive during the convergent extension of tissue.

    1. Computational and Systems Biology
    2. Physics of Living Systems
    Taegon Chung, Iksoo Chang, Sangyeol Kim
    Research Article

    Locomotion is a fundamental behavior of Caenorhabditis elegans (C. elegans). Previous works on kinetic simulations of animals helped researchers understand the physical mechanisms of locomotion and the muscle-controlling principles of neuronal circuits as an actuator part. It has yet to be understood how C. elegans utilizes the frictional forces caused by the tension of its muscles to perform sequenced locomotive behaviors. Here, we present a two-dimensional rigid body chain model for the locomotion of C. elegans by developing Newtonian equations of motion for each body segment of C. elegans. Having accounted for friction-coefficients of the surrounding environment, elastic constants of C. elegans, and its kymogram from experiments, our kinetic model (ElegansBot) reproduced various locomotion of C. elegans such as, but not limited to, forward-backward-(omega turn)-forward locomotion constituting escaping behavior and delta-turn navigation. Additionally, ElegansBot precisely quantified the forces acting on each body segment of C. elegans to allow investigation of the force distribution. This model will facilitate our understanding of the detailed mechanism of various locomotive behaviors at any given friction-coefficients of the surrounding environment. Furthermore, as the model ensures the performance of realistic behavior, it can be used to research actuator-controller interaction between muscles and neuronal circuits.