Active morphogenesis of patterned epithelial shells
Abstract
Shape transformations of epithelial tissues in three dimensions, which are crucial for embryonic development or in vitro organoid growth, can result from active forces generated within the cytoskeleton of the epithelial cells. How the interplay of local differential tensions with tissue geometry and with external forces results in tissue-scale morphogenesis remains an open question. Here, we describe epithelial sheets as active viscoelastic surfaces and study their deformation under patterned internal tensions and bending moments. In addition to isotropic effects, we take into account nematic alignment in the plane of the tissue, which gives rise to shape-dependent, anisotropic active tensions and bending moments. We present phase diagrams of the mechanical equilibrium shapes of pre-patterned closed shells and explore their dynamical deformations. Our results show that a combination of nematic alignment and gradients in internal tensions and bending moments is sufficient to reproduce basic building blocks of epithelial morphogenesis, including fold formation, budding, neck formation, flattening, and tubulation.
Editor's evaluation
The article provides a physical description of shape transformations of epithelial tissues in three dimensions, subject to active forces generated within the cytoskeleton of the epithelial cells. The work is motivated by organoids and more generally by morphogenesis during development. Therefore, this study is useful not only for developmental biology but also for a general understanding of cellular properties, including membrane mechanics and cell shapes.
https://doi.org/10.7554/eLife.75878.sa0Introduction
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.
Model
Viscoelastic nematic active surface model for epithelial mechanics
We first discuss our mechanical description of the deforming tissue. We represent an epithelium as an active surface flowing with velocity (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 , but also the bending moment tensor 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 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 , which we adopt for convenience. The local force balance projected on the tangential and normal directions reads (Salbreux and Jülicher, 2017)
where notations of differential geometry are introduced in Appendix 1; briefly is the curvature tensor, denotes the metric tensor, and the antisymmetric Levi-Civita tensor, the vector normal to the surface, is the tangential contribution of the tension tensor and its normal contribution, and 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:
We assume an external force density acting on the surface in addition to a difference of hydrostatic (uniform) pressure , 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, . Dissipative couplings to the external fluid are ignored here as the characteristic viscosity of a biological tissue ( Pa s; Marmottant et al., 2009; Guevorkian et al., 2010) is several orders of magnitude larger than that of water ( Pa s).
Constitutive equations
In line with our hypothesis describing the material properties of an epithelium, we use the following constitutive equations:
where is the symmetric part of the tension tensor and, on a curved surface, the strain rate tensor and the corotational time derivative of the curvature tensor are given by (Salbreux and Jülicher, 2017)
with the normal component of the vorticity. is the area strain, measuring local changes of area relative to a reference value; a precise definition is introduced in Equation 14. 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 and the bending modulus . The dynamical deformations of the surface are characterised by the two-dimensional shear and bulk viscosities and and the bulk bending viscosity . 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 . 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 , , , correspond to active tensions and bending moments. is an isotropic active surface tension, is the in-plane nematic active stress, with usually referred to as the ‘contractile’ active stress and as the ‘extensile’ active stress (Marchetti et al., 2013). is the isotropic bending moment, which locally favours a spontaneous curvature . If the active surface corresponds simply to two parallel layers under surface tension , (such as an epithelium with apical surface tension and basal surface tension ), and separated by a distance , an active isotropic bending moment emerges in the surface to lowest order in the curvature tensor. The term in corresponds to an anisotropic active bending moment. In the bilayer picture, where the active surface corresponds to two layers and , it could generally arise from differences between the two layers in the level of order and or in the level of nematic active stress and . For example, such differences could stem from two contractile (respectively extensile) layers with perpendicular nematic orientations and (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 and , the anisotropic components proportional to and , 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 of the tensor and (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 . 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):
with the Frank elastic constant , which is assumed to be equal for all distortions. The Landau–de Gennes contribution is chosen such that for the aligned state with is a minimiser for . 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 is parametrised by the arc length coordinate and the angle of rotation as
The local tangent basis is given by , and is the outward-pointing surface normal. The geometry of axisymmetric surfaces is described further in Appendix 1. We require that the metric component , which implies relations between the tangent angle and the shape functions and
which, together with the meridional principal curvature
are sufficient to reconstruct the surface shape from the curvature . In this axisymmetric setup, the velocity field reads , with the tangential and the normal velocities.
The undeformed initial surface is a sphere with radius , and all quantities defined on it are denoted with a subscript ‘0’. We define the area strain on a point of the surface as
where is the surface area element at the point considered on the surface, and is the surface area element of the same material point on the sphere. With this definition, on the initial sphere. We denote the arc length position on the undeformed sphere of a material point at arc length position on the deformed sphere. One then has with the meridional stretch and the circumferential stretch. Integrating yields the arc length reparametrisation between the initial and the deformed surface. The Lagrangian time derivative of the area strain (Equation 14) is related to the flow through
Nematic order
Here, with axial symmetry, the nematic tensor has the non-zero component . 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 vanishes there, creating defect cores of size , which is the characteristic nematic length. In this geometry the Euler–Lagrange equation resulting from the free energy (Equation 9) is
An example solution of Equation 16 on the sphere is shown in Figure 3b. From the two possible states with in the bulk, respectively, we choose 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 corresponds to circumferential active contraction, resulting in an elongated shape. For nematic bending moments, if one chooses to represent the order parameter on the outer side of the shell, the sign convention is such that , 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 and the active moment , 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 (Figure 1b), such that the active terms are given on by step-like profiles, for example
and similarly for , , and . 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 . 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 using sigmoid functions of the form
for their space and time dependence. For instance, the active bending moment profile is defined on as
as a smooth version of the step-profile Equation 17, and , , and 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).
Volume
We consider two possibilities for the volume enclosed by the epithelium. In one limit the tissue is assumed to be impermeable and the enclosed volume is treated as an incompressible fluid exerting hydrostatic pressure on the tissue. The volume is conserved when the shell deforms:
with the pressure 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, . In dynamical simulations, we introduce a volume viscosity such that the pressure is coupled to the volume change via
where 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 , and, if applicable, the nematic equilibrium Equation 16.
Dynamical deformations
In the dynamical version of the model a given active profile generates a velocity , whose normal part deforms the surface (Figure 1b). The components 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 . Since and 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 ,
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 .
Dimensionless variables
The equations are made dimensionless (marked by tilde) by rescaling tensions by , bending moment densities by , lengths by , force densities by , viscosities by the two-dimensional shear viscosity of the epithelium, times by the characteristic time scale , and velocities by . This leaves the dimensionless parameters , , , and to be fixed. We choose to set , for fast relaxation of the volume, and the nematic length scale is set to . 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 of the cell layer, and express . In simulations we use , corresponding to , 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 .
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 , and geometrical singularities at the poles are handled using analytical limits at (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 . 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 . The time integration according to Equation 22 is done with an explicit Euler method with adaptive step size via
In order to keep the force and torque balance equations in the form given by Equations 63–65, the updated surface is reparametrised as in a new arc length which is calculated from the condition . The profiles and surface quantities are passed between time steps as spline interpolants.
To produce the diagrams of steady-state shapes, 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 . A solution branch is found by starting from the spherical solution at zero difference of active profile, and calculating a sequence of steady-state shapes, progressively increasing the magnitude of the difference in activity. Two different methods are used to construct the solution branch for a sequence of control parameter values. For small values, starting from zero, the solution branch is obtained by making small increments in the control parameter. For larger values we switch to an implicit stepping method, which we developed based on a parametric representation of the solution branch (see Appendix 6 section ‘Construction of solution branches’). This second method allows us to continue the solution branches into regions where the steady-state shapes become non-unique in the control parameter.
Details of the numerical methods can be found in Appendices 6 and 7 for the steady state and the dynamics simulations, respectively.
Results
Epithelia as active membranes: Isotropic active tensions
We first consider deformations of an epithelial shell due to patterns of isotropic active tensions and bending moments. A spatially varying isotropic tension represents a change in the preferred area of the epithelium due to either changes in sheet thickness or cell number (Popović et al., 2017). However, one can show that a step-profile of positive (contractile) tension 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 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 and the size of the active region . The steady-state deformed shapes are plotted with the active region shown in red and the ‘passive’ region, where , 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, (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 in Equation 6) induces a preferred curvature , such that regions of a spherical shell with can be expected to flatten or bend inwards. Specifically, a difference of applied at the boundary of the active cap induces a jump in meridional curvature 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 for intermediate values of (Figure 2a). In particular, for the same value of the active region may bend inward or keep a positive curvature, depending on its size.
When is small or close to 1, the resulting shape is characterised by the formation of a bud which form either inwards () or outwards . In these cases, for sufficiently large values of 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 ; extrapolation indicates full constriction at a finite (Figure 2i). As the neck radius decreases the principal curvatures at the neck diverge as , such that 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, with and the radius of the two spheres and 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 changes sign across the neck, while 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 () and radii , , related by the condition:
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 gives an additional relation for and as a function of . Combining these conditions determine a curve in the parameter space , , 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 () and close to values of (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 and for intermediate values of , a fold in the solution branch in the -plane appears (Figure 2c). For most values of , this fold is associated to the loss of a continuously attainable solution with increasing , 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, ) 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, ). Alternatively the solution branch continues in a sequence of loops and the active region elongates (Figure 2c, ), 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, ). 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, ). 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 (), or from a sphere to a torus via self-intersection ().
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 , the shell deforms into locally folded shapes, which eventually self-intersect at large bending moment difference (Figure 2g, , Figure 2h).
Nematic active tensions
We now introduce the nematic order parameter 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, (Figure 3). Two solutions with and can exist; in the following we take the convention that , , corresponding to circumferential alignment of the order parameter, such that a contractile active stress results in a positive circumferential tension, . Due to invariance of the constitutive equation by exchange , , the same shape deformations occur when considering meridional alignment of the order parameter () and exchanging contractile and extensile active stresses.
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, (Figure 4a–e).
In the extensile case (or in the contractile case if ), and no pressure difference across the shell, the surface progressively flattens into a flat, double-layered disc (Figure 4b, , ). There is no shape transition occurring; instead, we find that the shape converges to a limit shape as (Appendix 4). The limit shape corresponds to two parallel flat discs of radius , separated by a distance , 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 :
The first relation shows that the limit shape has the size of the characteristic nematic length . Physically, for , 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 (), the shape elongates until a shape transition is reached, characterised by a fold in the solution branch (Figure 4b, , ). 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, as a result of the diverging principal curvatures (as can be seen from the presence of a term 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 and extensile stress in the active region , the active region tends to flatten more and more strongly as is increased, and the total curvature vanishes at the south pole (, Figure 4d). For and contractile stress , a fold in the solution branch appears at large value of (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 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 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 (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 . As for nematic tension, we adopt the convention of nematic alignment along the circumference, ; alignment along the meridians can be studied simply by changing the sign of the active coefficient .
We first discuss the case where the nematic active bending moment is homogeneous (), where there is no difference of pressure across the surface, and where (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 . To characterise this, we note that the corresponding steady-state shape solutions have vanishing tensions and everywhere (Figure 5—figure supplement 1) and the force balances Equations 63 and 64 are trivially satisfied. The torque balance Equation 65 reads
Combining Equations 26 and 48 one obtains that , with the operator . Solutions to have the form with A a constant. The boundary condition that the function should be finite at the poles requires , such that
As a result, if the shape has a cylindrical part, in which and , then the cylinder radius is given by
and since such solutions are area-preserving, with , the length of the cylindrical part scales as . These relations are in excellent agreement with simulation results for large enough (Figure 5g).
When , the active region forms an outward cylindrical protrusion (Figure 5a, b and g) whose radius is still well described by Equation 28, replacing by , the value of the active nematic bending moment in the active region (Figure 5g). Using that within the cylindrical protrusion 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
which is again in excellent agreement with numerical simulation for large and for different values of (Figure 5g).
For and we find erythrocyte-like shapes, where the indentations at the poles become stronger with until the two poles touch (Figure 5b). This behaviour remains for , 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 and the volume is free to change, both signs of result in a cylindrical appendage forming from the active region. The sign of 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 , 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 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 and the surface becomes spindle-like, narrowing at the poles with increasing . As in the free volume case, when 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 . Shape transitions occur at large , for both and (Figure 5e). In the case , 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 , 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 determines whether the active region folds inwards or outwards, independent of the initial size . As before, we use dynamics simulations to study the deformations for large (Figure 5f). For both signs of , these result in shapes that are self-intersecting either along a circle () or at the poles ().
Discussion
In this study of deformations of patterned nematic active surfaces, we have found a diverse zoology of possible shape changes (Figure 6), characterised by budding and neck constrictions, transition of sphere to cylinder, tubulation, and flattening. We find that introduction of a nematic field on the surface greatly increases the space of possible shapes. Overall our work contributes to the characterisation of the ‘morphospace’ which biological systems can explore.
Some of our findings recapitulate epithelial deformations observed in biological systems. The flattening observed for an extensile homogeneous nematic surface (Figure 4b, ) could in principle lead to merging of the two apposed surfaces into a double-layer for large . 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 : 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 , , or on the trace or determinant of the curvature tensor . 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 be a curved surface embedded in and parametrised by the generalised coordinates . A local covariant basis in the tangent plane is given by the vectors and the unit normal vector is constructed as and chosen to point outwards for a closed surface in our convention. These define the metric tensor and the curvature tensor . The infinitesimal surface and line elements are given by and , where is the determinant of the metric tensor. The antisymmetric Levi–Civita tensor is defined as .
Covariant derivatives
The Christoffel symbols of the second kind are obtained from derivatives of the metric as
and Christoffel symbols of the first kind are defined as . The covariant derivatives of a tangent vector field and a tangent tensor field , respectively, are given by
In the following, we also use the divergence theorem on curved surfaces for a tangent vector field (Salbreux and Jülicher, 2017),
Infinitesimal variation of surface quantities
For a small deformation of the surface
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)
These relations can be used to obtain time derivatives of surface quantities using (Lagrangian surface update) or (where the surface shape is updated with the normal flow only).
Since according to the definition Equation 14 the area strain can be written as, using the coordinates on the undeformed and deformed surfaces and denoting the determinant of the metric of the undeformed surface,
we obtain from Equation 38 the variation:
which yields the Lagrangian time derivative (Equation 15) in the main text, using .
Axisymmetric surfaces
Fundamental tensors
On a surface with axial symmetry about the -axis, as defined in the main text, the basis vectors and the outward normal are
where we have used and , which can be defined through the requirement that is an arc length parameter, such that . The metric and curvature tensors and the surface element are given by
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 we use and . The circumferential and meridional principal curvatures and , and the mean and Gaussian curvatures and are given by:
Some useful relationships involving the two principal curvatures follow from Equation 46 and from the definitions Equations 11 and 12,
The partial area and partial volume are given by
The corotational time derivative of the curvature tensor, as defined in Equation 8, has trace
Covariant derivatives
Axial symmetry implies that all functions on the surface should be -independent, and we also consider here vector and tensor fields , such that and . The only non-vanishing component of is ; therefore, the non-zero Christoffel symbols of the first kind are
and the non-zero Christoffel symbols of the second kind are
The non-zero components of are
resulting in the components of the strain rate tensor defined in Equation 7:
with trace
The non-zero components of are
where 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
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, which vanishes for an axisymmetric surface where . Therefore, here 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 bounded by the south pole and a circle perpendicular to the axis of symmetry, given by (Appendix 1—figure 1). The bounding circle has the line element and the unit normal , tangent to the surface and pointing outward with respect to . The balance of forces acting on reads
The contributions to the integral Equation 66 are
where we have introduced the integrated external force
and used the shape Equation 11 in Equation 68. As Equations 67–69 only contribute to the -component, Equation 66 can be rewritten as
From Equation 71 one obtains an expression for the transverse tension for :
and it is easy to confirm, using Equation 63, that this is indeed a solution of the normal force balance given by Equation 64.
Behaviour at the poles
The poles of the axisymmetric surface are at (south pole) and (north pole) and satisfy . Besides, we assumed that the shape has a finite curvature, requiring and at the south and north poles. The asymptotic behaviour of the shape is then:
The limits of geometric singularities of the form , for some function which vanishes at the pole, follow from L’Hôpital’s rule
where the + sign applies to and the - sign to . For example, applying Equation 75 to yields that
Any smooth tangent vector field on the closed axisymmetric surface has to vanish at the poles. For example, since is a scalar field, for its derivative we have
From Equation 48 we find that at the poles , which, together with Equation 77 yields
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 and south pole at the origin, :
where the arc length is denoted , and , 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 (). In that case, the tensors and are isotropic and the force and torque balance Equations 63–65 can be written, in the absence of external force other than the pressure and for :
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 . At the point of the surface closest to the axis of symmetry, and . We then scale the arc length coordinate , the distance of the surface to the axis of symmetry , and the curvature tensor with , and introduce , and . The force balance equations then become for :
For , the leading order solution has and both constant. Using the relation
with , and the conditions and that does not diverge for , the curvatures have solution and
where the sign of changes in the regions towards and away from the neck. Therefore, converges to +1 or -1 away from the neck for . Since is constant across the infinitesimal neck, also changes sign asymptotically away from the neck.
The next order in of Equation 90 gives which corresponds to constant across the neck (Figure 2—figure supplement 1).
Analytical solution for the free volume case
In the free volume case, and the force balance equations admits the solution and constant (Figure 2—figure supplement 1a). Considering a shape with the active and passive regions forming spheres with radii and separated by an infinitesimal neck, the condition of constant results in
with and the trace of the curvature tensor in the active and passive regions, and the difference in isotropic active bending moment between the passive and active regions. This results in Equation 24, taking into account that depending on whether the active region is curved towards the outside or the inside part of the surface. In addition, the condition results in conservation of area of the active and passive regions compared to the undeformed sphere:
where we have used Equation 85. When corresponding to the active region towards the outside, and , Equation 24 implies that which further requires to satisfy Equation 94.
Combining Equations 93 and 94 gives a condition defining a curve in the parameter space , :
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 . In an orthonormal frame , where and , it is characterised by an angle as
The director components in the basis are then
The traceless and symmetric nematic order parameter can be constructed from the director and a magnitude as
Its components read
We assume here that there is no azimuthal flow on the surface, . The -component of the tangential force balance then reads for constant , :
which requires excluding divergence of at the poles; therefore, the only possible orientations for the director, compatible with our assumption of vanishing azimuthal flows, are . Therefore, in the axisymmetric setup we consider only one non-zero component
which can take the values , corresponding to azimuthal or longitudinal orientation of the director , respectively.
The total free energy of the nematic Equation 9 reads in terms of
We minimise this energy with respect to on a given surface. The resulting Euler–Lagrange Equation 16 is obtained from
where
Here, we have used the definition from the functional derivative, . Equation 16 can be written as two first-order equations
The requirement that Equation 106 should be regular at the poles of the surface results in the two boundary conditions
and also implies
The limit of Equation 106 at is given by
and equivalently at . Therefore, Equation 106 does not provide a limit value for at the poles. When solving Equations 105 and 106 numerically, we use the analytical limits at the poles
with two free parameters and , 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 on a sphere of radius to axisymmetric perturbations; a more general analysis can be found in Napoli and Vergori, 2012. We note that for a spherical shape, with and at first order in :
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 is provided by taking derivatives of axisymmetric spherical harmonics, with with the Legendre polynomial of degree , for . One then finds
The isotropic state is stable if for , or for
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 (Figure 4b). We consider here the limit where .
We postulate that the limit shape reached as consists of two parallel flat central discs of radius , separated by a distance , and connected by a narrow curved region (Appendix 4—figure 1). One denotes the arclength of the extreme position of the shape where is maximal, and 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 . One looks for an asymptotic shape which satisfies ; we show later that this condition requires .
The force and torque balance Equations 63–65 can then be rewritten
with . The last equation implies that at (i.e. at the point of the surface with the extremal value of ), and therefore at this point . However, our definition of deformation implies that . Therefore, as , one must have (Appendix 4—figure 2), and the equilibrium equation for the nematic order parameter , Equations 105 and 106 can be linearised.
Introducing a renormalised order parameter and renormalised tension , the force and torque balance equation and the linearised equilibrium equation for the order parameter read
to which one can add the boundary conditions , , and a condition on which follows from Equation 143:
In the asymptotic regime where , 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 , and Equation 119 becomes
with solution, using , , with a constant to determine, and are the Bessel functions of the first kind. The condition then yields the expression for the radius of the disc:
with the smallest positive solution of ().
With this solution at hand, solving Equation 117 gives
using that . The constant can be determined from Equation 120:
which can be rewritten:
where one has used in the last line, following Equation 124. is a decreasing function from to and . In the limit , Equation 126 is satisfied provided that , which sets the constant in that limit. Because , this implies .
In the following, we denote the deformation reached at the end of the circular plate. Because of the arguments given above, when , . In general, implies that : we assume this is the case in the following.
Narrow curved connecting region
In the narrow curved region, at leading order in small , and are homogeneous, and . The force and torque balance Equations 117 and 118 now give
which can be combined to obtain, using because of the symmetry of the shape:
We look for a solution of this equation , with the boundary conditions , and, using ,
It is then helpful to introduce the following differential equation:
which admits a solution which can be found numerically. The solution of Equation 129 can then be written
and the constraint Equation 130 gives
where we have introduced , with the approximate numerical value . Together with Equation 122, this analysis indicates that the distance converges to a constant value as (Appendix 4—figure 2). Overall, this analysis predicts the limit value:
such that the condition implies that the analysis is self-consistent for .
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 (). 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
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
Because at low Reynolds number the total force acting on the surface must vanish, the value of should be set to zero by the numerical solver. Indeed, inspection of Equation 72 shows that the conditions also imply , and since is constant, . In practice, 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:
where , , , . Here, Equation 138 has been obtained by using Equation 72 to rewrite the normal force balance (Equation 64) as
with the integral of the force defined in Equation 136. The term in Equation 147 is regular at the poles by construction since it has the Taylor expansions at the south pole and at the north pole. Therefore, the limits are and .
The geometric boundary conditions for the shape are , . Further, we require , , and . The equations for the nematic require and , as discussed in Appendix 3. If , the unknown length is determined from the condition . If , the domain of integration is split into and . The additional condition sets L1. All unknown functions are matched at the internal boundary , except for the curvature and strain. These may acquire a jump, and , to ensure continuity in and .
If applicable, it is convenient to formulate the volume constraint as a boundary value problem for the partial volume defined in Equation 51,
Due to the geometric singularities appearing in several of the equations at the poles of the surface and , we derive the limits of these equations there. Denoting the solution vector by the limiting expressions of Equations 137–148 at the south and north poles, respectively, are
In summary, when considering the full interval and conserved volume, the boundary conditions are
and the free parameters are , , , and . Otherwise, if the volume is free to change, then and the condition is removed.
In the dimensionless equations, arc lengths are transformed to the unit interval by if , and to two unit intervals in the case of step-profiles. With and , the dimensionless variables are with in the first interval and with 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 and , leading to typical adaptive grid sizes of , 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 , where is the vector space spanned by all (or a subset of) parameters of the model, for example, and , 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 the current state and by the previous state of the system and approximate by the tangent vector at the current state. To find a new solution on the curve, a step of length 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,
which allows us to eliminate one (e.g. the first) component of p(i+1) via
provided that and where denotes the vector without the first component, etc. The remaining parameters are determined by the solver such that the new point is a continuation of the solution branch, which can be achieved for small enough . It suffices to include only a subset of the free parameters in this construction, for example, 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 , can be written as
where in Equation 156 one has to replace
The solution vector
is determined from
where the linear operator is constructed from Equations 155–158 together with two trivial relations relating to and to . The system of ode’s (Equation 162) is solved on the full interval with the boundary conditions
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
where
In Equations 165 and 166, we have used that and , defined as continuous functions on the closed surface, satisfy
Prior to solving system Equation 162, at every time step the nematic profiles and are determined on the shape as solutions of the Euler–Lagrange Equation 16,
with the boundary conditions , as discussed in Appendix 3.
Constraints
Volume
The volume of a closed surface reads
According to Equations 34, 36, and 38, for a Lagrangian surface update with we have
This allows to calculate from Equation 172 the rate of change in volume
where we have used two integrations by part and the relations , and .
If the volume is fixed then for all the dynamics has to satisfy
On the axisymmetric surface, this results in the integral constraint . We define a partial rate of volume change , such that the constraint Equation 179 can be written as
This is solved simultaneously with the boundary value problem (Equations 162–168), where the pressure 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:
In this case as the dynamics simulation approaches a steady state, in agreement with the direct steady-state calculations in which .
Rigid-body translation
We note that in the absence of external force, for any flow profile solution of the force and torque balance Equations 155–158, the flow
with an arbitrary constant , is again a solution. The addition of the uniform flow field corresponds to a rigid-body translation in the -direction which does not affect force balance and therefore makes the task of numerically determining 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
For a Lagrangian surface update one obtains
where we have used Equation 174 to obtain , and the divergence theorem on curved surfaces (Equation 33). It follows that the velocity of the centroid is given by
and we want to fix
for all . We note that tangential velocity components do not contribute to the centroid displacement.
On the axisymmetric surface and the constraint Equation 189 becomes
Analogously to the volume constraint, we introduce a partial centroid velocity and add the constraint Equation 190 to the boundary value problem (Equations 162–168) in the form of
Note that due to the free parameter 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 consists of two sub-steps: first, the material points are advected using a Lagrangian coordinate as given in Equation 23, then the surface is reparametrised in a new arc length coordinate , for which .
From Equation 22 we find the time updates for the shape descriptors:
where the last equation can be shown by taking the time derivative of the normal vector and using Equation 36. It is convenient to derive the time update for and its derivative from the constitutive Equation 6 and the torque balance (Equation 65),
because the trace of the corotational derivative (Equation 8) is the Lagrangian time evolution of , as can be seen by using Equation 39 with . From the variation (Equation 39) the circumferential curvature is updated as
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
with the trace of the strain rate tensor given in Equation 59.
The displaced surface is reparametrised in a new arc length coordinate ,
such that , or equivalently . From
we obtain
so the relationship between the two arc length parameters is given by the differential equation
We rewrite this equation in terms of a rate of change of arc length ,
and it is added to the linear system (Equations 162–164). The new arc length is obtained as and the perimeter length is updated as
For the active profiles defined via sigmoid functions, as given in the main text in Equation 18, the parameters chosen for the simulations are , , , and . In the time step the bending moment profile, as defined at time via Equation 19, is updated as
and an analogous relation holds for , , and . For the spatial dependence, we keep track of the arc length on the initial sphere 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 . For example,
where the Jacobian prefactor for the derivative is given by Equation 200.
The size of the adaptive time step is determined using a standard step doubling method (Press et al., 2007), where the shape after a full time step, denoted by , is compared to that after two half steps, denoted by . The relative error is calculated from the shape components and the curvature derivative as
on a uniform grid which is defined on and projected onto and , whereby the poles and values of and , which are too close to zero (), 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 , which defines . As an example we show in Appendix 7—figure 1 the convergence results for the active bending profile
with , , 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 ) is calculated as
on a dimensionless, uniform grid , , with , which is obtained through division by or , respectively, from the corresponding simulation. As discussed in Appendix 5, the deviation of the parameter (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 normalised by the pressure . We find good convergence of the dynamics to the steady-state result with decreasing tolerance on the time step. Based on these results, we take 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: , .
Appendix 8
Isotropic active tension
Consider first a shell with vanishing internal hydrostatic pressure, , and a step-like tension profile given on the reference surface by
One can verify that the spherical shape, given by , with , is a solution of the Equations 137–142. The strain has a jump , such that
Using this to solve Equation 143 for with yields
and similarly for
These two equations determine the unknown radius and the deformed active region size as
How the ratio changes with the tension jump is plotted in Appendix 8—figure 1 for . The above rescaling holds only for since the active region contracts to a point for .
For conserved volume, the spherical shape with , , is a solution of Equations 137–142. Here, only the relative size of the active patch changes from to . As before the strain is piecewise constant with a jump at , , and integrating Equation 143 with , gives the additional conditions
from which , , can be determined.
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 , with one control parameter . The line of steady state solutions is given by
where solutions are parametrised by . Taking the derivative one obtains . We consider a fold in the solution curve where and . In that case:
and the linear stability matrix at the fold thus has (at least) one eigenvalue that changes sign at the fold. Assuming that the solution branch is stable up to the fold, this indicates generically the appearance of an unstable mode; except in the special case where the 0 eigenvalue is a local maximum.
Data availability
The current manuscript is a computational study, so no data have been generated for this manuscript. Modelling code is available at the GitHub repository https://github.com/DianaKhoromskaia/EpithelialShell, (copy archived at swh:1:rev:0838f09f1b2228d8d7da5183fc68f3b49c6ee734).
References
-
Fission of a multiphase membrane tubePhysical Review Letters 93:158104.https://doi.org/10.1103/PhysRevLett.93.158104
-
Bending lipid membranes: experiments after W. helfrich’s modelAdvances in Colloid and Interface Science 208:47–57.https://doi.org/10.1016/j.cis.2014.02.002
-
Integer topological defects of cell monolayers: mechanics and flowsPhysical Review. E 103:012405.https://doi.org/10.1103/PhysRevE.103.012405
-
Stresses in lipid membranesJournal of Physics A 35:6233–6247.https://doi.org/10.1088/0305-4470/35/30/302
-
BookInternational Series of Monographs on PhysicsOxford, United Kingdom: Oxford University Press.
-
Spontaneous shear flow in confined cellular nematicsNature Physics 14:728–732.https://doi.org/10.1038/s41567-018-0099-7
-
Forces controlling organ growth and sizeMechanisms of Development 144:53–61.https://doi.org/10.1016/j.mod.2016.11.005
-
Primary invagination of the vegetal plate during sea urchin gastrulationAmerican Zoologist 24:571–588.https://doi.org/10.1093/icb/24.3.571
-
Defects in nematic membranes can buckle into pseudospheresPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 77:041705.https://doi.org/10.1103/PhysRevE.77.041705
-
Size control: the developmental physiology of body and organ size regulationWiley Interdisciplinary Reviews. Developmental Biology 4:335–356.https://doi.org/10.1002/wdev.181
-
Aspiration of biological viscoelastic dropsPhysical Review Letters 104:218101.https://doi.org/10.1103/PhysRevLett.104.218101
-
Intrinsic bending force in anisotropic membranes made of chiral moleculesPhysical Review. A, General Physics 38:3065–3068.https://doi.org/10.1103/physreva.38.3065
-
Dynamics of a Volvox embryo turning itself inside outPhysical Review Letters 114:178101.https://doi.org/10.1103/PhysRevLett.114.178101
-
Abbildungsklassenn-dimensionaler mannigfaltigkeitenMathematische Annalen 96:209–224.https://doi.org/10.1007/BF01209163
-
The hope and the hype of organoid researchDevelopment 144:938–941.https://doi.org/10.1242/dev.150201
-
Spontaneous symmetry breaking and pattern formation of organoidsCurrent Opinion in Systems Biology 11:123–128.https://doi.org/10.1016/j.coisb.2018.06.002
-
Vesicle shape, molecular tilt, and the suppression of necksPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 76:031908.https://doi.org/10.1103/PhysRevE.76.031908
-
Domain-Induced budding of vesiclesPhysical Review Letters 70:2964–2967.https://doi.org/10.1103/PhysRevLett.70.2964
-
Shape transformations of vesicles with intramembrane domainsPhysical Review E 53:2670–2683.https://doi.org/10.1103/PhysRevE.53.2670
-
Perspective: the promise of multi-cellular engineered living systemsAPL Bioengineering 2:040901.https://doi.org/10.1063/1.5038337
-
A BVP solver based on residual control and the maltab PSEACM Transactions on Mathematical Software 27:299–316.https://doi.org/10.1145/502800.502801
-
Buckling of spherical capsulesPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 84:046608.https://doi.org/10.1103/PhysRevE.84.046608
-
Gastrulation in the sea urchin embryo: a model system for analyzing the morphogenesis of a monolayered epitheliumDevelopment, Growth & Differentiation 46:309–326.https://doi.org/10.1111/j.1440-169x.2004.00755.x
-
Curvature control of valence on nematic shellsSoft Matter 7:670–683.https://doi.org/10.1039/C0SM00378F
-
Mechanics of nonplanar membranes with force-dipole activityPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 73:061913.https://doi.org/10.1103/PhysRevE.73.061913
-
Orientational order and vesicle shapeJournal de Physique II 2:371–382.https://doi.org/10.1051/jp2:1992133
-
Orientational order, topology, and vesicle shapesPhysical Review Letters 67:1169–1172.https://doi.org/10.1103/PhysRevLett.67.1169
-
Stability and foldsArchive for Rational Mechanics and Analysis 99:301–328.https://doi.org/10.1007/BF00282049
-
Activating membranesPhysical Review Letters 112:258101.https://doi.org/10.1103/PhysRevLett.112.258101
-
Hydrodynamics of soft active matterReviews of Modern Physics 85:1143–1189.https://doi.org/10.1103/RevModPhys.85.1143
-
Self-organized cytoskeletal alignment during Drosophila mesoderm invaginationPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 375:20190551.https://doi.org/10.1098/rstb.2019.0551
-
Topology and morphology of self-deforming active shellsPhysical Review Letters 123:208001.https://doi.org/10.1103/PhysRevLett.123.208001
-
Minimal model of cellular symmetry breakingPhysical Review Letters 123:188101.https://doi.org/10.1103/PhysRevLett.123.188101
-
Active morphogenesis of epithelial monolayersPhysical Review 100:022413.https://doi.org/10.1103/PhysRevE.100.022413
-
Emergence of active nematic behavior in monolayers of isotropic cellsPhysical Review Letters 122:048004.https://doi.org/10.1103/PhysRevLett.122.048004
-
Surface free energies for nematic shellsPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 85:061701.https://doi.org/10.1103/PhysRevE.85.061701
-
N -atic order and continuous shape changes of deformable surfaces of genus zeroEurophysics Letters 20:279–284.https://doi.org/10.1209/0295-5075/20/3/015
-
Geometrical control of active turbulence in curved topographiesPhysical Review Letters 122:168002.https://doi.org/10.1103/PhysRevLett.122.168002
-
Active dynamics of tissue shear flowNew Journal of Physics 19:033006.https://doi.org/10.1088/1367-2630/aa5756
-
BookNumerical Recipes 3rd Edition: The Art of Scientific ComputingCambridge University Press.
-
BookTheory and Analysis of Elastic Plates and ShellsFlorida, United States: CRC Press.https://doi.org/10.1201/9780849384165
-
Progress and potential in organoid researchNature Reviews. Genetics 19:671–687.https://doi.org/10.1038/s41576-018-0051-9
-
Irreversible thermodynamics of curved lipid membranesPhysical Review. E 96:042409.https://doi.org/10.1103/PhysRevE.96.042409
-
Hydrodynamics of cellular cortical flows and the formation of contractile ringsPhysical Review Letters 103:058102.https://doi.org/10.1103/PhysRevLett.103.058102
-
Mechanics of active surfacesPhysical Review. E 96:032404.https://doi.org/10.1103/PhysRevE.96.032404
-
Theory of nematic and polar active fluid surfacesPhysical Review Research 4:033158.https://doi.org/10.1103/PhysRevResearch.4.033158
-
Shape transformations of vesicles: phase diagram for spontaneous- curvature and bilayer-coupling modelsPhysical Review. A, Atomic, Molecular, and Optical Physics 44:1182–1202.https://doi.org/10.1103/physreva.44.1182
-
Configurations of fluid membranes and vesiclesAdvances in Physics 46:13–137.https://doi.org/10.1080/00018739700101488
-
Theory of chiral lipid tubulesPhysical Review Letters 71:4091–4094.https://doi.org/10.1103/PhysRevLett.71.4091
-
Theory of cylindrical tubules and helical ribbons of chiral lipid membranesPhysical Review E 53:3804–3818.https://doi.org/10.1103/PhysRevE.53.3804
-
Theoretical model for the formation of caveolae and similar membrane invaginationsBiophysical Journal 86:2049–2057.https://doi.org/10.1016/S0006-3495(04)74266-6
-
Confined bilayers passively regulate shape and stressPhysical Review Letters 110:028101.https://doi.org/10.1103/PhysRevLett.110.028101
-
Morphogenesis and compartmentalization of the intestinal cryptDevelopmental Cell 45:183–197.https://doi.org/10.1016/j.devcel.2018.03.024
-
Concise theory of chiral lipid membranesPhysical Review E 76:603.https://doi.org/10.1103/PhysRevE.76.031603
-
Furrow constriction in animal cell cytokinesisBiophysical Journal 106:114–123.https://doi.org/10.1016/j.bpj.2013.11.014
-
Cell fate coordinates mechano-osmotic forces in intestinal crypt formationNature Cell Biology 23:733–744.https://doi.org/10.1038/s41556-021-00700-2
Article and author information
Author details
Funding
Cancer Research UK (C55977/A23342)
- Guillaume Salbreux
Cancer Research UK (FC001317)
- Guillaume Salbreux
Medical Research Council (FC001317)
- Guillaume Salbreux
Wellcome Trust (FC001317)
- Guillaume Salbreux
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. For the purpose of Open Access, the authors have applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission.
Acknowledgements
We thank S Grigolon for useful discussions and N Cuny for comments on the manuscript. DK and GS acknowledge support from the Francis Crick Institute, which receives its core funding from Cancer Research UK (FC001317), the UK Medical Research Council (FC001317), and the Wellcome Trust (FC001317), and from the CRUK multidisciplinary project award C55977/A23342.
Copyright
© 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.
Metrics
-
- 1,731
- views
-
- 282
- downloads
-
- 12
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
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)
Further reading
-
- Physics of Living Systems
We propose the Self Returning Excluded Volume (SR-EV) model for the structure of chromatin based on stochastic rules and physical interactions. The SR-EV rules of return generate conformationally defined domains observed by single-cell imaging techniques. From nucleosome to chromosome scales, the model captures the overall chromatin organization as a corrugated system, with dense and dilute regions alternating in a manner that resembles the mixing of two disordered bi-continuous phases. This particular organizational topology is a consequence of the multiplicity of interactions and processes occurring in the nuclei, and mimicked by the proposed return rules. Single configuration properties and ensemble averages show a robust agreement between theoretical and experimental results including chromatin volume concentration, contact probability, packing domain identification and size characterization, and packing scaling behavior. Model and experimental results suggest that there is an inherent chromatin organization regardless of the cell character and resistant to an external forcing such as RAD21 degradation.
-
- Physics of Living Systems
A hallmark of biomolecular condensates formed via liquid-liquid phase separation is that they dynamically exchange material with their surroundings, and this process can be crucial to condensate function. Intuitively, the rate of exchange can be limited by the flux from the dilute phase or by the mixing speed in the dense phase. Surprisingly, a recent experiment suggests that exchange can also be limited by the dynamics at the droplet interface, implying the existence of an ‘interface resistance’. Here, we first derive an analytical expression for the timescale of condensate material exchange, which clearly conveys the physical factors controlling exchange dynamics. We then utilize sticker-spacer polymer models to show that interface resistance can arise when incident molecules transiently touch the interface without entering the dense phase, i.e., the molecules ‘bounce’ from the interface. Our work provides insight into condensate exchange dynamics, with implications for both natural and synthetic systems.