Morphoelasticity of large bending deformations of cell sheets during development

Deformations of cell sheets during morphogenesis are driven by developmental processes such as cell division and cell shape changes. In morphoelastic shell theories of development, these processes appear as variations of the intrinsic geometry of a thin elastic shell. However, morphogenesis often involves large bending deformations that are outside the formal range of validity of these shell theories. Here, by asymptotic expansion of three-dimensional incompressible morphoelasticity in the limit of a thin shell, we derive a shell theory for large intrinsic bending deformations and emphasize the resulting geometric material anisotropy and the elastic role of cell constriction. Taking the invagination of the green alga Volvox as a model developmental event, we show how results for this theory differ from those for a classical shell theory that is not formally valid for these large bending deformations and reveal how these geometric effects stabilize invagination.


I. INTRODUCTION
During animal and plant development, cell division, cell shape changes, and related processes can drive deformations of cell sheets [1][2][3][4][5][6].In elastic continuum theories of the development of the green alga Volvox [7][8][9][10], of tissue folding in Drosophila [11,12], or of more abstract active surfaces [13], these driving processes appear as changes of the reference or intrinsic geometry of thin elastic shells.
Just as classical thin shell theories arise from an asymptotic expansion of bulk elasticity in the small thickness of the shell [14][15][16], these "morphoelastic" shell theories should be asymptotic limits of a bulk theory.While there is now a well-established framework of three-dimensional morphoelasticity [17,18], based on a multiplicative decomposition of the deformation gradient tensor into intrinsic and elastic deformations [19], studies of this asymptotic limit appear to have been restricted to the case of flat morphoelastic plates.Extensions of the classical Föppl-von Kármán equations [20,21] have been derived for this case, but we are not aware of any corresponding study for curved morphoelastic shells.Rather, such studies have remained more phenomenological: some models [7,8,[11][12][13] simply replaced the elastic strains in classical shell theories [15,22,23] with measures of the difference of the intrinsic and deformed geometries.Other studies [9,10] took a more geometric approach, mirroring geometric derivations of classical shell theories [22] based on the so-called Kirchhoff "hypothesis".This is the asymptotic result [15] that the normals of the midsurface of the undeformed shell remain, at leading order, normal to the deformed midsurface.
There is however one more serious limitation of these models: tissues in development undergo large bending deformations (Fig. 1) that are outside the formal range of validity of the underlying thin shell theories, which assume that the thickness of the shell is much smaller than all lengthscales of the midsurface of the shell [15,22,23].However, even if the thickness of the cell sheet is much smaller than its undeformed radius of curvature, this radius of curvature may become comparable, locally, to the thickness of the cell sheet as it deforms (Fig. 1).This is associated with cells contracting at one cell pole to splay and thereby bend the cell sheet [4].
Here, we derive a theory of thin incompressible morphoelastic shells undergoing large bending deformations by asymptotic expansion of three-dimensional elasticity.We reveal how, even in a constitutively isotropic material, this new scaling limit of large bending deformations induces, in the thin shell limit, a geometric anisotropy absent from classical shell theories: different deformation directions exhibit different deformation responses.We stress how this geometric effect is associated with the geometric singularity of cell constriction, i.e. the limit of wedged triangular cells [Fig.1(b), inset] associated with these large bending deformations.Specialising to the invagination of the green alga Volvox [25,26], we then FIG. 1.Large bending deformations during morphogenesis: even if the thickness of the cell sheet is small compared to the undeformed radius of curvature, the local radius of curvature need not remain small compared to the cell sheet thickness as the sheet deforms.(a) Cross section of ventral furrow formation in Drosophila, reproduced from Ref. [24].(b) Midsagittal cross section of invagination in the spherical alga Volvox globator, reproduced from Ref. [8].Inset: cartoon of constricted triangular cells in the bend region.Scale bars: 20 µm.
show how results for this theory differ from those for a classical theory that is not formally valid in this large bending limit, and reveal how invagination is stabilised by the geometry of large bending deformations.

II. ELASTIC MODEL
In this section, we describe large bending deformations of a thin incompressible morphoelastic shell, starting from threedimensional morphoelasticity.We shall have to distinguish between three configurations of the shell [Fig.2(a)]: (i) the undeformed configuration of the shell, (ii) the deformed configuration of the shell, and (iii) the intrinsic configuration of the shell that encodes the local, intrinsic deformations of the shell, i.e. the cell shape changes or cell division in the biological system.These intrinsic deformations are not in general compatible with the global geometry of the shell.Elasticity must therefore intervene to "glue" the intrinsically deformed patches of cell sheet back together, as illustrated in Fig. 2(a).Configurations (i) and (ii) are related by the geometric deformation gradient F. This tensor decomposes multiplicatively into an intrinsic contribution F 0 that relates configurations (i) and (iii), and an elastic contribution F = F F 0 −1 .This is the multiplicative decomposition of morphoelasticity [17,18].
In this section, we restrict to torsionless deformations of an axisymmetric shell.The analysis can be extended to more general deformations of the shell, and, for the sake of completeness, we do so in Appendix A, but the restriction to axisymmetric deformations eschews the mire of tensorial notation that arises in the general case.
The derivation of the shell theory for large bending deformations divides, like derivations of classical shell theories, into two steps: first, in subsection II A, we describe the kinematics of the deformations and derive expressions for the geometric, intrinsic, and elastic deformations gradients.Second, in subsection II B, we analyse the mechanics of the shell and expand the three-dimensional elastic energy and equilibrium conditions asymptotically.At the end of this section, in subsection II C, we discuss the limit of small bending deformations that gives rise to classical shell theories.

A. Axisymmetric deformations of an elastic shell
We consider an elastic shell of undeformed thickness εh, where ε 1 is a small asymptotic parameter expressing the thinness of the shell compared to other lengthscales associated with its midsurface.Large bending deformations will be introduced by allowing one of the intrinsic radii of curvature of the shell to be of order O(ε).We now derive an expression for the elastic deformation gradient F for torsionless deformations of an axisymmetric shell.

Undeformed configuration of the shell
We shall describe the undeformed configuration V of the shell with reference to a midsurface S that we shall choose later.With respect to the basis {u r , u φ , u z } of cylindrical coordinates, we define the position vector of a point on S , with s denoting arclength and φ being the azimuthal coordinate [Fig.2(b)].The tangent angle ψ(s) of S is defined by in which dashes denote differentiation with respect to s.The vectors e s (s, φ) = cos ψ(s)u r (φ) + sin ψ(s)u z , e φ (φ) = u φ (φ) (3) thus constitute a basis of the tangent space of S [Fig.2(c)], which we extend to an orthonormal basis B = e s , e φ , n for V by adjoining the normal to S , We complete the description of S by computing its curvatures, Now, the position of a point in V is where we have introduced the transverse coordinate ζ [Fig.2(c)].Noting the partial derivatives ∂ n/∂s = s e s and ∂ n/∂φ = φ e φ , we obtain

Deformed configuration of the shell
As the shell deforms into its deformed configuration Ṽ , the midsurface S maps to the deformed midsurface S [Fig.2(d)], with position vector where, in particular, s is again the undeformed arclength.Denoting by s the deformed arclength, we define the stretches which enable us to define the tangent angle ψ(s) of S by r (s) = f s cos ψ(s), z (s) = f s sin ψ(s), (10) where dashes still denote differentiation with respect to s.
Similarly to the analysis of the undeformed configuration, we introduce the tangent vectors and the normal vector to define an orthonormal basis B = ẽs , ẽφ , ñ to describe Ṽ [Fig.2(e)].The curvatures of the deformed shell are As the shell deforms, the normals to S need not remain normal to S , and so a point in V at a distance εζ from S will end up, in Ṽ , at a distance ε ζ from S , and displaced by a distance ε ς parallel to S [Fig.2(e)].By definition of the midsurface, Finally, using commata to denote partial differentiation, we note the partial derivatives and

Calculation of the deformation gradient tensors
Combining Eqs.(7) and (15), we obtain an expression for the geometric deformation gradient tensor and hence write down an expression for the intrinsic deformation gradient tensor, viz.
In the expression for F 0 , we have introduced the intrinsic stretches f 0 s , f 0 φ , the intrinsic curvatures κ 0 s , κ 0 φ , and the intrinsic displacement ζ 0 from the intrinsic midsurface S 0 of the intrinsic configuration V 0 of the shell.In writing down this expression for F 0 , we have assumed that there is no intrinsic displacement parallel to the midsurface, ς 0 = 0.In Eq. ( 16), F is expressed with respect to the mixed basis B ⊗ B, while F 0 is expressed with respect to another mixed basis, B 0 ⊗ B, where B 0 = {E s , E φ , N } is a (local) orthonormal basis for the intrinsically deformed configuration of the shell.The elastic deformation gradient is thus expressed here with respect to the natural mixed basis B ⊗ B 0 .Since B and B 0 are orthonormal, there exists a rotation, represented by an orthogonal matrix R, that maps B onto B 0 .With respect to B 0 ⊗ B 0 , the elastic deformation gradient tensor is then represented by the matrix RF.This rotation can, as explained below, be in fact ignored in the calculations that follow.This simplification is not however possible for general deformations, as discussed in Appendix A.

B. Thin shell theory for large bending deformations
In this subsection, we derive the effective elastic energy for the shell by asymptotic expansion of three-dimensional elasticity.We assume the simplest constitutive law, that the shell is made of an incompressible neo-Hookean material [17], so that its elastic energy is wherein C is a material parameter, and I 1 is the first invariant of the right Cauchy-Green tensor C = F F [17].We emphasise that the integration of the strain energy density e is over the intrinsic configuration V 0 of the shell, which has volume element dV 0 .The (first) Piola-Kirchhoff stress tensor for this neo-Hookean material [21] is, with respect to B 0 ⊗ B 0 , with Q = F − pF − , and where the Lagrange multiplier p is proportional to pressure and imposes the incompressibility condition det F = 1.The second equality follows, since R is orthogonal, from The configuration of the shell minimising the energy E in Eq. ( 18) is determined by Div 0 P = 0 [17], where the divergence is taken with respect to V 0 .Separating into components parallel and perpendicular to the midsurface and using Eq. ( 19), since R is independent of ζ 0 and orthogonal, where the nabla operator denotes the gradient on S 0 , and where J = I − N ⊗ N denotes projection onto the intrinsic midsurface, normal to N .

Scaling assumptions
We now introduce large intrinsic bending deformations explicitly by scaling the intrinsic curvatures so as to allow small radii of curvature in the meridional direction, viz.
wherein the rescaling by the prefactor f 0 s f 0 φ absorbs, as do similar rescalings below, the intrinsic stretching of the midsurface.This scaling regime in which the meridional intrinsic curvature becomes comparable to the thickness of the cell sheet is the one relevant for Volvox invagination, as shown in Fig. 1(b).Appendix A treats the general case in which all components of the curvature tensor are allowed to be large.
Next, we make the standard scaling assumption of shell theory, that the elastic strains are small, i.e. that the stretches and curvatures in the deformed configuration do not differ "too much" from the intrinsic stretches and curvatures.More formally, we introduce the shell strains E s , E φ by writing and the curvature strains L s , L φ by letting Finally, we introduce the scaled variables

Boundary and incompressibility conditions
We now impose the incompressibility condition det F = 1 and the boundary conditions, which require that there be no external forces on the surfaces of the shell.These boundary conditions are relevant for many problems in developmental biology, where deformations are, as discussed in the introduction, driven by changes of the intrinsic geometry only; including external forces does not pose any additional difficulty, though.
These force-free boundary conditions read P ± N ± = 0 [17], where, from Eq. ( 19), P ± = CRQ ± are the Piola-Kirchhoff stress tensors evaluated on the shell surfaces of the intrinsic configuration V 0 , and N ± are the unit outward normals to these shell surfaces.Since the rotation R is an isometry, the boundary conditions become To complete the derivation of the boundary conditions, we must obtain expressions for the normals N ± to the shell surfaces.At this stage, we complete specifying the midsurface by choosing S 0 so that the shell surfaces are at ζ 0 = ±h 0 /2.Similarly to the calculations leading up to Eqs. (15), we compute the unit tangent vectors to the shell surfaces in the intrinsic configuration, expressed here with respect to Introducing a normalisation factor and expanding while recalling that κ 0 s = O ε −1 , we obtain
b. Solution at orders O(ε) and O ε 2 .We now expand the incompressibility condition further, finding On solving the resulting differential equation for Z (1) by imposing Z (1) = 0 at Z 0 = 0, we obtain Next, we formally expand the entries of the deformation gradient in Eq. ( 17) further by writing In particular, we obtain ) , w (1) can be obtained in terms of the expansions ( 27), but will turn out to be of no consequence.Hence the incompressibility condition becomes Next, we notice that Q (0) = O from Eq. ( 36), and hence Eq. ( 20) from Eq. ( 36), we find Q (1) N = ε(v (1) + w (1) ), 0, O(ε) .From this and from Eq. ( 38), we infer

Asymptotic expansion of the constitutive relations
On computing the expansion of the eigenvalues of F from Eq. ( 36) and hence of I 1 , and simplifying using Eqs.(39), we obtain Hence, from Eqs. (37) and on introducing x = λ 0 s Z 0 and This determines the leading-order term of the asymptotic expansion of the energy density in Eq. ( 18).On defining, from Eq. ( 36), the (symmetric) effective two-dimensional deformation gradient and associated two-dimensional strain tensor, wherein I is the identity, we rewrite Eq. (40a) as This shows how, at leading order, the energy density depends only on the two invariants of the effective two-dimensional strain tensor.In the asymptotic limit of a thin shell, the constitutive relations have thus become effectively two-dimensional.

Intrinsic volume conservation
To complete the derivation of the shell energy, we must specify the intrinsic deformations, and in particular how ζ 0 and hence Z 0 depend on ζ.Volume conservation between the undeformed and intrinsic configurations of the shell requires equality of the corresponding volume elements, dV = dV 0 , which condition becomes (43) This, on recalling that κ 0 s = O ε −1 yields, at leading order, a differential equation for Z 0 (ζ), on imposing Z 0 = 0 at ζ = 0.The shell surfaces are at ζ 0 = ±h 0 /2 in the intrinsic configuration, and at ζ = ±h ± in the undeformed configuration, where h + + h − = h is the undeformed thickness of the cell sheet.Let H 0 = h 0 f 0 s f 0 φ .Eq. ( 44) yields Finally, we introduce η = λ 0 s h, so that the shell surfaces are at x = ±η/2.We note that Eq. ( 45) is a leading-order result only, since we have ignored O(ε) corrections in Eq. ( 44).

Derivation of the thin shell theory
We are now set up to average out the transverse coordinate and thus obtain the thin shell theory.We note, from Eq. ( 43), the expression for the volume element in the intrinsic configuration, On substituting Eqs.(40b) and ( 46) into Eq.( 18), integrating with respect to x, and using axisymmetry, we obtain with the first integration over the undeformed axisymmetric midsurface S and the second over the curve C generating S .The effective two-dimensional energy density ê in Eq. ( 47a) is wherein only.All of these coefficient functions diverge as η → ±2; this divergence, absent from theories not valid for large bending deformations, is not surprising.Indeed, the limit η → ±2 corresponds to constricted cells, i.e. wedge-shaped, triangular cells [Fig.1(b), inset] for which the radius of curvature is half the intrinsic cell sheet thickness: either the apical or basal surface has become geometrically singular by contracting to a point.As the intrinsic configuration approaches this constricted limit somewhere, deviations from the intrinsic configuration become more and more expensive energetically there compared to other positions in the cell sheet.Fig. 3 makes a more general point by plotting the coefficient functions in Eqs.(48), arbitrarily scaled with α ss to absorb their divergence as η → ±2.This shows how the relative importance of different deformation modes depends on the amount of intrinsic bending.In other words, large bending deformations break the material isotropy, so that different directions of stretching have different effective stretching moduli; similarly, different effective bending moduli are associated with different directions of bending.This anisotropy is therefore geometric; as discussed in more detail below, this effect is absent from previous theories not valid for large bending deformations.
Moreover, rearranging Eqs. ( 22) and ( 23), the shell strains in Eq. (47b) are while the curvature strains are where we have defined At leading order (i.e. at the order to which the shell theory is valid), we may choose to replace the curvature strains L s , L φ in Eq. (47b) with the alternative curvature strains K s , K φ defined above.The latter choice may be more natural since K s , K φ vanish for pure stretching deformations, whereas L s , L φ do not: consider a shell, the undeformed (and intrinsic) configuration of which is a sphere of radius R, and which deforms into a sphere of radius R = f R, for example because of a pressure difference between the inside and outside.For this deformation, In other words, Eqs.(51a) show that the stretching deformations associated with changes in curvature can be neglected at leading order in the shell theory.Ref. [15] has also discussed this point.
This completes the derivation of the elastic energy (47a) of a thin shell undergoing large axisymmetric bending deformations.In Appendix B, we derive the governing equations associated with Eq. (47a).We solve these equations numerically using the boundary value problem solver bvp4c of M -(The MathWorks, Inc.) and the continuation software [28].

C. Limit of small bending deformations
We conclude our calculations by taking the limit η → 0, in which the bending deformations become small compared to the thickness of the shell.The energy density in Eq. (47b) then limits to the form familiar from classical shell theories, up to corrections of order O ε 4 .This is the energy density of a thin Hookean shell [15,22,23] with Poisson's ratio ν = 1/2, implying incompressibility, and elastic modulus E = 6C.In particular, our analysis also provides a formal derivation of the morphoelastic version of this classical shell theory.The first term in the sum in square brackets is the stretching energy, while the second term is the bending energy.To pick up on a point made earlier, we note that, in this theory, the same stretching modulus E h 1 − ν 2 = 8Ch and the same bending modulus E h 3 12 1 − ν 2 = 2Ch 3 /3 are associated with all directions of stretching or bending; it is this isotropy resulting from the constitutively assumed isotropy of the material that is broken by the geometry of large bending deformations.Of course, Eq. ( 52) could be derived directly by imposing different scalings, of small intrinsic bending, replacing those for large bending deformations in Eq. ( 21); these scalings would considerably simplify the solutions of Eqs. ( 29), (30), and (34).Indeed, the structure of these calculations would be broadly similar to the earlier asymptotic derivation of the classical shell theories in Ref. [16].We are not however aware of previous work pointing out, as we did above, that the terms at order O ε 2 in the expansion (36) of the deformation gradient need not be computed explicitly.
We note that Eq. ( 52) has the same structure as the elastic energies used in previous models referenced in the introduction, but the morphoelastic definitions of the shell and curvature strains in Eqs. ( 50) and (51b) differ from those in these previous models: in models not based on morphoelasticity and its multiplicative decomposition of the deformation gradient [7,8,[11][12][13], the shell and curvature strains are simply differences of stretches or curvatures, missing the scaling factors of f 0 s , f 0 φ that appear in Eqs. ( 50) and (51b).We also note that the expressions for the curvature strains in Eqs.(51b) differ, by a factor of f 0 s f 0 φ , from those in Refs.[9,10], which, as discussed earlier, derived morphoelastic shell theories based on a geometric approach.The geometric role of this factor has been noticed previously in the context of uniform growth of an elastic shell [29].
Moreover, the geometric approach in Refs.[9,10] leads to additional terms in the energy density.The present analysis proves that these terms are not leading-order terms in the thin shell limit.However, there is no reason to expect this geometric approach to yield all terms at next order in the asymptotics.A complete expansion could in principle be obtained by continuing the asymptotic analysis presented here.This would permit asymptotic justification of the so-called shear deformation theories [30] in which the normals to the undeformed midsurface need not remain normals in the deformed configurations, but we do not pursue this further here.

A. Biological background
The green algal genus Volvox [31] has become a model for the study of the evolution of multicellularity [32,33], for biological fluid dynamics [34], and for problems in developmental biology [35,36].Adult Volvox colonies [Fig.4(a)] are spheroidal, consisting of several thousand biflagellated somatic cells that enclose a small number of germ cells [31].Each germ cell undergoes several rounds of cell division to form a spherical embryonic cell sheet [Figs.4(b) and 4(e)], at which stage those cell poles whence will emanate the flagella point into the sphere [31].To acquire motility, the embryo turns itself inside out in a process called inversion [25,37].
In some species of Volvox [25,26], inversion starts with the formation of a circular invagination [Figs.4(c) and 4(f)], reminiscent of the cell sheet folds associated with processes such as gastrulation or neurulation in higher organisms.At the cell level, this invagination results from two types of cell shape changes [7,26]: (1) cells near the equator become wedgeshaped [Fig.4(d)], while the cytoplasmic bridges (cell-cell connections resulting from incomplete division) rearrange to connect the cells at their thin wedge ends, and (2) cells in the posterior hemisphere narrow in the meridional direction.These cell shape changes arise simultaneously, with (1) splaying the cells and thereby bending the cell sheet [Fig.4(d)] and (2) contracting the posterior hemisphere to facilitate the subsequent inversion of the posterior hemisphere inside the as yet uninverted anterior hemisphere.In later stages of inversion, other cell shape changes arise in different parts of the cell sheet [9,26] to facilitate the peeling of the anterior hemisphere over the inverted posterior and thus complete inversion.

B. Results
Following our earlier work [7][8][9][10], we model Volvox inversion by considering the deformations of an incompressible elastic spherical shell under quasi-static axisymmetric variations of its intrinsic stretches and curvatures representing the cell shape changes driving inversion.The slow speed of inversion-it takes about an hour for a Volvox embryo to turn itself inside out [25,26]-justifies this quasi-static approximation.In more detail, Figs.4(g) and 4(h) show functional forms of the intrinsic stretches and curvatures encoding the cell shape changes driving invagination and define the model parameters κ p , κ b , κ a , f p , f a , s 0 , and w that encode the intrinsic curvatures and intrinsic stretches of different regions of the cell sheet and the extent of these regions.In numerical calculations, we regularise the step discontinuities in the definitions of the intrinsic stretches and curvatures in Figs.4(g) and 4(h), we non-dimensionalise all lengths with the pre-inversion radius R of the embryo, and we take εh = 0.15, appropriate for Volvox globator [7,9].During the invagination stage, the radius of curvature in the bend region of wedge-shaped cells [Fig.4(f)] becomes comparable to the thickness of the cell sheet: this is the scaling limit of large bending deformations studied in Section II.We therefore compare the resulting elastic model, with energy density (47b), to the classical theory, in which the energy density is given by Eq. (52).For weakly invaginated stages of Volvox inversion (corresponding to small values of η in the large bending theory), the two models yield, unsurprisingly, very similar shapes [Fig.5(a)], mirrored by very similar profiles of meridional curvature [Fig.5(b)].However, the more the intrinsic configuration of the cell sheet approaches the limit of cell constriction, the more the shapes resulting from the two models differ [Fig.5(c)].The curvature profile of these shapes differ in particular in the bend region of nearly constricted cells [Fig.5(d)], as expected: since the coefficients defined in Eqs.(48) diverge in the constriction limit, deviations from the intrinsic configuration become increasingly expensive there in the large bending theory, but not in the classical theory.For this reason, the cell sheet is more invaginated in the large bending theory than in the classical theory [Fig.5(c)].In this way, the geometric effects of large bending stabilise the invagination process.
These examples indicate that the results of the two models differ at a quantitative, if not at a qualitative level.We extend this observation by plotting k = −κ b against the displacement d of the posterior pole [Fig.5(e), inset] for different values of the width w of the bend region in Fig. 5(e).Again, the solution curves show similar behaviour in the two models, but differ at a quantitative level.There is a critical bend region width w * such that the solution curves in the (k, d) diagram are single-valued for w < w * , but become multivalued for w > w * , leading to discontinuous jumps in d as k is varied.Where multiple solutions exist for a given value of k, the one with the lowest value of d has the lowest energy (not shown).For the classical theory, we have discussed this bifurcation behaviour in Ref. [8], and rationalised it by constructing an effective energy that estimates different elastic contributions.It is therefore not surprising that, here, we find qualitatively identical bifurcation behaviour in the two models, but, again, there are quantitative differences in the bifurcation behaviour of the two theories.In particular, we note that w * is larger in the large bending theory than in the classical theory.In other words, continuous invagination is possible in a larger region of parameter space in the large bending theory than in the classical theory, so, again, the geometry of large bending deformations is stabilising.

IV. CONCLUSION
In this paper, we have derived a morphoelastic shell theory valid for the large bending deformations that are commonly observed in developmental biology (Fig. 1), and have shown how this new scaling limit of large bending deformations induces a purely geometric effective material anisotropy absent from classical theories.Taking the invagination of the green alga Volvox as an example, we have compared this theory to a simpler, classical theory not formally valid for these large bending deformations.Since the classical theory does not account for the geometric singularity of cell constriction, it differs, for strongly invaginated shapes as in Figs.1(b), 4(c), and 4(f), from the theory for large bending deformation at a quantitative, if not at a qualitative level.In particular, we have argued that this geometric limit stabilises invagination.
This and the growing interest in quantitative rather than merely qualitative analyses of morphogenesis [38,39] emphasise the importance of this new scaling limit of large bending deformations for studies of the mechanics of developmental biology.The theory we have derived here is not however the most general theory of these large bending deformations.Indeed, when writing down the expression for the intrinsic deformation gradient in Eq. ( 16), we assumed that there is no intrinsic displacement parallel to the midsurface, ς 0 = 0.For ς 0 0, we would replace the expansion of S posited in Eq. ( 27) with S = S 0 + S (0) + O(ε), where S 0 = f 0 s f 0 φ ς 0 .The nonlinear differential equations extending Eqs. ( 29) and (30) that arise in the expansions of the boundary and incompressibility conditions in this more general case still admit a trivial solution S (0) ≡ 0, Z (0) ≡ 0, p (0) = 1, but we were unable to extend our calculations in Section II to prove that this solution is unique; we raise a similar issue in Appendix C. It therefore remains unclear what form the extension of the Kirchhoff "hypothesis" [15] to this case takes.
In this paper, we assumed the simplest, incompressible Neo-Hookean constitutive relations when deriving our shell theory for large bending deformations.The restriction to incompressible elastic materials is justified by the biological context of our analysis, in which the models derived here describe sheets of fluid-filled cells that are therefore indeed incompressible to a first approximation.However, the bulk elastic response of biological materials such as brain tissue is not linear [40][41][42].The restriction to linear Neo-Hookean relations may therefore appear to be a limitation of the analysis, but that turns out not to be the case: in the thin shell limit, general hyperelastic constitutive relations reduce to Neo-Hookean relations.This result has been established previously for thin plates [20,43], and, in Appendix C, we (partially) extend it to the large bending deformations of thin shells considered here.In the context of shell theories, the problem of specifying the nonlinear constitutive relations of biological tissues does not therefore arise.However, we have recently shown that the continuum limit of a class of discrete models of cell sheets involves not only nonlinear elastic, but also nonlocal, nonelastic terms [44].Moreover, adding the geometric singularity of apical constriction (corresponding to triangular cells in the underlying discrete model) as a constraint to the variational problem that arises in this continuum limit remains an important open problem [44], solving which may provide a regularisation of the singularity that arises in the theory derived here.Meanwhile, this suggests that the journey towards understanding the continuum mechanics of biological materials, on which we have taken another step with the present analysis of large bending deformations of thin elastic shells, will continue to abound with new problems in nonlinear mechanics.

Non-axisymmetric deformations of an elastic shell
As in Section II, we begin by deriving an expression for the deformation gradient tensors of an elastic shell of thickness εh, where ε is, again, a small asymptotic parameter that expresses the thinness of the shell.

a. Undeformed configuration of the shell
We parameterise the undeformed midsurface S of the shell in terms of generalised, not necessarily orthogonal coordinates; we shall use Greek letters to denote these coordinates.Thus, if ρ is the position of a point on S , the tangent vectors there are e α = ∂ρ/∂α.The metric g of the midsurface thus has components g αβ = e α • e β , and we set g = det g.
Next, we define a basis B for the shell by adjoining the normal vector n to this tangent basis.The Weingarten and Gauß equations [45] in which commata denote partial differentiation, express the derivatives of the normal and tangent vectors in terms of the curvature tensor αβ = e α • n ,β and the Christoffel symbols associated with the metric of the midsurface [45].
The position of a point in the undeformed configuration V of the shell is r = ρ + εζ n, where ζ denotes the transverse coordinate, as defined in Fig. 2(c).Hence wherein we have used the Weingarten equation (A1), and where δ is the Kronecker delta.The metric tensor G of the undeformed configuration therefore has components

Deformed configuration of the shell
We take the same generalised coordinates to parameterise the deformed midsurface of the shell, and define the tangent vectors ẽα at a point ρ on its midsurface S .The metric g of the midsurface has components gαβ = ẽα • ẽβ .We extend the tangent basis of the midsurface to a basis B for the deformed shell by adding the unit normal ñ, and introduce the corresponding curvature tensor καβ = ẽα • ñ,β .Extending the setup in Fig. 2(e) to non-axisymmetric deformations, the position of a point in the deformed configuration Ṽ of the shell is where ζ and ςα are the transverse and parallel displacements of this point relative to the midsurface.In particular, for these non-axisymmetric deformations, the displacement parallel to the midsurface is no longer a scalar.Using the Weingarten and Gauß equations (A1), we find in which ςβ ;α = ςβ ,α + Γ αγ β ςγ is a covariant derivative.

c. Intrinsic configuration of the shell
We define the intrinsic configuration of the shell by specifying the intrinsic basis B 0 of the shell, containing the vectors E α , which we think of as the intrinsic tangent vectors, and a unit vector N normal to these tangent vectors.The metric g 0 of the intrinsic midsurface has components g 0 αβ = E α • E β , and we write g 0 = det g 0 .Further, we define the (symmetric) intrinsic curvature tensor κ 0 αβ = E α • N ,β , and the metric tensor G 0 of the intrinsic configuration.Denoting by ζ 0 the intrinsic transverse displacement, G 0 has, by analogy with Eqs.(A2) and (A3), components and When writing down these expressions, we have assumed, as we did in Section II, that there is no intrinsic displacement parallel to the midsurface, ς 0 α = 0.

d. Calculation of the deformation gradient tensors
The deformation gradient tensor relating the undeformed and deformed configurations of the shell is F = r,α ⊗ r ,α by definition.Now, from Eq. (A3), It follows that or, in block matrix notation, in which the asterisk denotes a dual basis, and where we have introduced and we write By analogy with Eqs.(A8), the intrinsic deformation gradient tensor is or, in block matrix notation, Here we have again assumed that there is no intrinsic displacement parallel to the midsurface, ς 0 α = 0, and we have introduced On inverting the block-lower triangular matrix in Eq. (A11b), we find The elastic deformation gradient is F = F F 0 −1 .From Eqs. (A8b) and (A13), we obtain We replace B with an intermediate basis B0 by replacing the tangent vectors using the transformation ẽα = U β α Ẽβ .Then The bases B or B0 and B 0 are not necessarily orthogonal or normalised, so the map that maps the first onto the second is not now in general a rotation, but we can around work around this issue using the following observation: since the unit normals ñ and N are orthogonal to Ẽα and E α respectively, we can choose Ẽα in such a way that B0 is a rotation of B 0 , represented by a matrix R. In particular, if (abusing notation) F denotes the matrix in Eq. (A14b), then the deformation gradient is represented, with respect to B 0 ⊗ B 0 * , by the matrix RF.Further Ẽα • Ẽβ = E α • E β since rotations preserve distances and angles.Hence, on recalling the definitions of the metrics gαβ and g 0 αβ of the deformed and intrinsic midsurfaces of the shell,

Thin shell theory for large bending deformations
As in Section II, we assume that the shell is made of an incompressible neo-Hookean material, with energy given by Eq. ( 18).Given our above remark on modifying the bases, Eq. ( 19) still provides the expression for the (first) Piola-Kirchhoff stress tensor P, now with respect to B 0 ⊗ B 0 * : we still have P = RQ, with Q = F − pF − , in which F is now given by Eq. (A14b).Moreover, Eq. ( 20) still applies.

a. Scaling assumptions
Again as in Section II, we now introduce large bending deformations explicitly by scaling the intrinsic curvature tensor so as to absorb the intrinsic stretching of the midsurface, Next, we make the standard scaling assumptions of shell theory, that the elastic strains remain small.First, we therefore require that the metrics g and g 0 be not "too different".Introducing the shell strain tensor E, and from Eq. (A15), we therefore impose For axisymmetric deformations, these tensors are symmetric, and we recover the definition of the shell strains in Eq. ( 22).Similarly, we require the curvature tensors of the intrinsic and deformed configurations to be not "too different".We therefore introduce the curvature strain tensor K by writing Finally, we introduce scaled variables

b. Intrinsic volume conservation
We now impose volume conservation of the intrinsic configuration of the shell compared to the undeformed configuration.We need one preliminary result: Lemma 1.Let M be a 2 × 2 matrix, and x be a scalar.Then Moreover, from Eqs. (A6) with the scalings introduced above and invoking Lemma 1, we find wherein H 0 = 1 2 tr λ 0 and K 0 = det λ 0 are the (scaled) mean and Gaussian intrinsic curvatures [45].The eigenvalues λ 1 , λ 2 of λ 0 are real [45], and H 0 = 1 2 (λ 1 + λ 2 ), K 0 = λ 1 λ 2 , from which follows the well-known inequality H 0 2 K 0 .Now, integrating the differential equation for Z 0 (ζ) resulting from Eqs. (A20) and imposing Z 0 = 0 at ζ = 0, we find Since Eq. (A20a) neglects O ε 2 corrections, this result holds at leading order only.Now, as discussed in Section II, the shell surfaces are at ζ 0 = ±h 0 /2 in the intrinsic configuration, and at ζ = ±h ± in the undeformed configuration, where h + + h − = h is the undeformed thickness of the cell sheet.On defining H 0 = h 0 g 0 /g, so that the shell surfaces are at Z 0 = ±H 0 /2 in the intrinsic configuration, Eq. (A21) yields and hence which is a depressed cubic equation for H 0 (h) that can be solved in closed form.In particular, Eq. (A22b) has a unique positive real solution if K 0 > 0, but has no positive real solution if h 2 K 0 < −16/9.If 0 > h 2 K 0 > −16/9, two positive real solutions exist; by continuity, the smaller must be chosen.More generally, we require that ζ increases with Z 0 , for Z 0 H 0 /2.As H 0 2 K 0 , the cubic in Eq. (A21) has The positions of the turning points at Z 0 = Z 0 ± are indicated, and ζ Z 0 must increase monotonically for Z 0 < H 0 /2.This condition excludes the dotted parts of the graphs.(b) Intrinsic volume conservation in K 0 h 2 , H 0 h space: conservation of intrinsic volume is only possible within the region of parameter space enclosed by the solid curve, in which −16/9 < h 2 K 0 < c = 64/9 and h H 0 < √ c = 8/3.The dashed lines delimit the regions of parameter space excluded by the inequality H 0 2 K 0 and the condition that Eq. (A22b) have a positive real solution.
two turning points [Fig.6(a)], at Z 0 = Z 0 ± , where explicit expressions for Z 0 − Z 0 + in terms of K 0 , H 0 can be found by solving a quadratic equation.The condition that ζ increases with Z 0 translates to inequalities Z 0 ± ≷ H 0 (h)/2 depending on the signs of K 0 , H 0 [Fig.6(a)].These inequalities involving h, H 0 , K 0 only depend on H 0 h and K 0 h 2 , since the curvatures can be nondimensionalised with h.The inequalities can then be solved numerically to determine the region in K 0 h 2 , H 0 h parameter space for which intrinsic volume conservation is possible [Fig.6(b)].In particular, Fig. 6(b) shows that intrinsic volume conservation requires −16/9 K 0 h 2 c and H 0 h √ c , where c is a numerical constant.We note that an expression for the boundary of this region can also be determined in closed form using M (Wolfram, Inc.); although rather complicated, this expression can be used to show that c = 64/9.
Finally, to relate this discussion to the results of Section II, we note that if K 0 = 0, the condition is H 0 1, which is equivalent to |η| 2, as expected.

c. Boundary conditions and condition of incompressibility
As in Section II, the incompressibility and boundary conditions read det F = 1 and Q ± N ± = 0.
To derive an expression for the normals N ± , we compute the tangent vectors E ± α to the surfaces of the intrinsic configuration of the shell, at ζ 0 = ±h 0 /2 ⇐⇒ Z 0 = ±H 0 /2.By analogy with Eq. (A5a) and using the scalings in Eqs.(A16) and (A19), We now order B 0 = {E 1 , E 2 , N } so that E 1 × E 2 = g 0 N , and hence E 1 × N = −E 2 , E 2 × N = E 1 .By expanding in components, we find On normalising this vector, we deduce (A25)

d. Expansion of the boundary and incompressibility conditions
To expand the boundary and incompressibility conditions, we extend Eqs. ( 27 and thence, from Eq. (A14b), in which dashes now denote differentiation with respect to Z 0 , and where, writing Ã = Ã(0) + O(ε), with Ã(0) given by Eq. (A27a), On computing the determinant of the block matrix [46] in Eq. (A28a), the incompressibility condition becomes Using further properties of block matrices [46] and from Eq. (A28a), we obtain + O(ε), (A30) cumbersome and therefore not presented here.For this reason, the general theory presented here is likely most useful to describe deformations with some additional symmetry, such as the axisymmetric deformations discussed in Section II.As in standard shell theories [23], the apparent singularity in the resulting equations is removed by introducing the transverse shear tension, T = −N s tan ψ, and we obtain, using Eqs.( 10) and ( 13 10) and ( 13), Eqs.(B6) determine the deformed configuration of the shell.Having solved these equations, integrating the otherwise redundant shape equation z = f s sin ψ from Eqs. (10) yields the shape of the shell.

Numerical solution of Eqs. (B6)
We conclude the derivation of the governing equations for axisymmetric deformations with two remarks on the numerical solution of Eqs.(B6).
First, we note that Eqs.(B6) are singular where r = 0.At such a point, geometric continuity implies ψ = 0. Hence T = 0 there by definition, and N φ = N s for regularity in Eq. (B6a).Moreover, by applying l'Hôpital's rule to the definitions in Eqs. ( 9) and ( 13 of which the first two follow by reflection across the axis of symmetry, and the last follows by applying l'Hôpital's rule to Eq. (B6c) and using the previous observations and Eqs.(B7).Second, as discussed also in Ref. [9,10], at each stage of the numerical solution, f s , f φ , κ s , κ φ must be determined from r, ψ, M s , N s .To begin with, f φ , κ φ and hence E φ , K φ are computed directly from r, ψ using the definitions (50) and (51b).We can then compute f s , κ s by noting that, once f φ , E φ , K φ are known, the definitions of N s , M s in Eqs.(B2a), (B2c), and (B4) define a system of linear equations for E s , K s .Their definitions (50) and (51b) then yield f s and finally κ s .We can then compute N φ , M φ using Eqs.(B2b), (B2d), and (B4), and thus continue the numerical integration.(Moreover, if r = 0, we similarly obtain two linear equations for f = f s = f φ and k = f s κ s = f φ κ φ , from the solution of which the numerical integration can be continued.)

FIG. 2 .
FIG. 2. Morphoelasticity of an axisymmetric shell.(a) The undeformed (top), deformed (left), and intrinsic (right) configurations of the shell are related by the three tensors F, F 0 , and F = F F 0 −1 .(b) Undeformed geometry of an axisymmetric shell of thickness εh, described by coordinates r(s), z(s), where s is arclength, with respect to the basis {u r , u φ , u z } of cylindrical polars.(c) Cross section of the undeformed shell, defining a local basis B = {e s , e φ , n} and the transverse coordinate ζ.(d) Deformed geometry of the shell: after a torsionless deformation, the shell has thickness ε h, arclength s, and is described by coordinates r(s), z(s) with respect to cylindrical polars.(e) Cross section of the deformed shell, defining a local basis B = { ẽs , ẽφ , ñ}.Normals to the midsurface rotate so that a point at a distance εζ from the undeformed midsurface S is at a distance ε ζ(s, ζ) from the deformed midsurface S , and displaced by a distance ε ς(s, ζ) parallel to S .

FIG. 4 .
FIG. 4. Invagination in Volvox.(a) Volvox colony, with somatic cells and one embryo labelled.(b) Light-sheet microscopy image of a spherical Volvox embryo before inversion.(c) Corresponding image at an early stage of inversion, when a circular invagination (I) has formed.(d) Splaying of cells and bending of the cell sheet result from the formation of wedge-shaped cells and the rearrangement of the cytoplasmic bridges (CBs); red lines indicate position of CBs.(e) Midsagittal cross-section of a Volvox embryo before inversion.(f) Corresponding cross-section during invagination, with the regions where wedge-shaped cells (W) and contracted spindle-shaped cells (C) have formed labelled.(g) Plot of the intrinsic curvature κ 0 s against arclength s, defined in the inset.The plot defines the model parameters κ p , κ b , κ a , s 0 , and w.Regions of cell shape changes (W, C) as in (f) are also indicated.(h) Corresponding plot of the intrinsic stretches f 0 s , f 0 φ , defining additional model parameters f p , f a .Panels (a)-(f) include microscopy images by Stephanie Höhn and have been redrawn from Ref. [8].Scale bars: (a) 50 µm; (e), (f) 20 µm.

FIG. 5 .
FIG.5.Comparison of the elastic model for large bending deformations (solid lines), with energy density given by Eq. (47b), and the classical model (dashed lines), with energy density given by Eq. (52).(a) Early invagination stage: the two models yield very similar shapes.Dotted line: undeformed configuration of the spherical shell.Parameter values:κ p = κ a = 1, κ b = −2, f p = 0.8, f a = 1, s 0 = 1.5, w = 0.2.(b) Corresponding plot of the meridional curvature and intrinsic curvature.(c) Later invagination stage: as the cells in the bend region approach the constriction limit, the shapes resulting from the two models start to differ.Parameter values are as in (a), except κ b = −8.5, w = 0.4.(d) Corresponding plot of the meridional curvature and intrinsic curvature.(e) Bifurcation diagram, for different values of w, in (k, d) space, with k = −κ b and with d the posterior displacement defined in the axis inset.Different lines correspond to parameter values w = 0.3, 0.5, 0.6, 0.7, 0.8, 0.9.Other parameter values are as in (a).The vertical line |η| = 2 corresponding to the constriction limit is also shown.Vertical arrows indicate discontinuous jumps in d as k is increased.