Mechanics and Variability of Cell Sheet Folding in the Embryonic Inversion of $Volvox$
Pierre A. Haas Stephanie S. M. H. Höhn Aurelia R. Honerkamp-Smith Julius B. Kirkegaard Raymond E. Goldstein


Integral Curvature Representation and Matching Algorithms
 for
Identification of Dolphins and Whales
    [
    December 30, 2023
====================================================================================================




Many embryonic deformations during development are the global result of local cell shape changes and other local active cell sheet deformations. Morphogenesis does not only therefore rely on the ability of the tissue to produce these active deformations, but also on the ability to regulate them in such a way as to overcome the intrinsic variability of and geometric constraints on the tissue. Here, we explore the interplay of regulation and variability in the green alga Volvox, whose spherical embryos turn themselves inside out to enable motility. Through a combination of light sheet microscopy and theoretical analysis, we quantify the variability of this inversion and analyse its mechanics in detail to show how shape variability arises from a combination of geometry, mechanics, and active regulation.




§ INTRODUCTION

Julian Huxley's pronouncement, “In some colony like [the green alga] Volvox, there once lay hidden the secret to the body and shape of [humans]” <cit.>, emphasises that morphogenesis across kingdoms relies on the fundamental ability of organisms to, firstly, produce active forces that drive the deformations of cell sheets underlying the development of many organs and tissues and, secondly, regulate these active deformations in such a way to complete morphogenesis. Unravelling the biomechanics of these processes is therefore of crucial importance to understand pathological errors and foster bioengineering to address these errors <cit.>. Local cellular changes can produce forces that are transmitted along the cell sheet to drive its global deformations <cit.>. Simple events of cell sheet folding such as ventral furrow formation in Drosophila can be driven primarily by cell shape changes <cit.>. In more complex metazoan developmental processes such as gastrulation <cit.>, optic cup formation <cit.>, neurulation <cit.> and related processes <cit.>, the effect of such cell shape changes is overlaid by that of other cellular changes such as cell migration, cell intercalation, cell differentiation, and cell division. 

In all of these processes however, these local cellular changes occur in specific regions of the cell sheet and at specific stages of morphogenesis. On the one hand, the spatio-temporal distribution of these local cellular changes affects the global tissue shape. On the other hand, a certain amount of noise is unavoidable in biological systems; indeed, it may even be necessary for robust development, as demonstrated for example by <cit.>, who showed that variability in cell growth is necessary for reproducible sepal size and shape in Arabidopsis. 
While some processes may be subject to less intrinsic variability than others, one must therefore ask: how are these processes orchestrated so that development can complete despite the intrinsic biological variability? Differences in the observed shapes of organisms at certain stages of development (i.e. what one might term their geometric variability) stem from a combination of mechanical variability (i.e. differences in mechanical properties or mechanical state) and active variability (i.e. differences in the active forces generated by individual cells). What experimental data there are suggest that the mechanical properties are subject to a large amount of variability <cit.>. Finally, differences in the mechanical stress state of the tissue are another facet of mechanical variability that is induced by active variability. 

The first mechanical models of morphogenesis <cit.> represented cells as discrete collections of springs and dashpots; they were soon followed by elastic continuum models <cit.>. Notable among this early modelling of morphogenesis is for example the work of  (, ), who combined models of several mechanisms of sea urchin gastrulation with measurements of mechanical properties to test the plausibility of these different mechanisms. These models heralded the emergence of a veritable plethora of mechanical modelling approaches over the subsequent decades <cit.>, though the choice of model must ultimately be informed by the questions one seeks to answer <cit.>. More recent endeavours were directed at deriving models that can represent the chemical and mechanical contributions to morphogenesis and their interactions <cit.> and at establishing the continuum laws that govern these out-of-equilibrium processes <cit.>. 

There is, however, a rather curious gap in the study of the variability of development: the importance of quantifying the morphogenesis and its variability has been recognised <cit.>, yet accounts of the variability of development, e.g. in the loach <cit.>, have often been merely descriptive. For this reason, the interplay between mechanics and active variability has seemingly received little attention and hence a question we believe to be fundamental appears to lie in uncharted waters: how does active variability lead to geometric variability? Conversely, what does geometric variability tell us about active variability?





This is the question that we explore in this paper in the context of the development of the multicellular green alga Volvox (Fig. <ref>a). Volvox and the related Volvocine algal genera have been recognised since the work of <cit.> as model organisms for the evolution of multicellularity <cit.>, spawning more recent investigations of kindred questions in fluid dynamics and biological physics <cit.>. The cells of Volvox (Fig. <ref>b) are differentiated into biflagellated somatic cells and a small number of germ cells, or gonidia, that will form daughter colonies <cit.>. The somatic cells in the adult are embedded in a glycoprotein-rich extracellular matrix <cit.>. The germ cells undergo several rounds of cell division, after which each embryo consists of several thousand cells arrayed to form a thin spherical sheet confined to a fluid-filled vesicle. Cells are connected to their neighbours by cytoplasmic bridges (Fig. <ref>b), thin membrane tubes resulting from incomplete cell division <cit.>. Those cell poles whence will emanate the flagella however point into the sphere at this stage, and so the embryos must turn themselves inside out through an opening at the anterior pole of the cell sheet (the phialopore), to enable motility and thus complete their development <cit.>. Because of this process of inversion, Volvox has become a model organism for the study of cell sheet deformations, too <cit.>.

Inversion in Volvox <cit.> and in related species <cit.> results from cell shape changes only, 
without the complicating additional processes found in metazoan development discussed above.  This simplification 
facilitates the study of morphogenesis.
While different species of Volvox have developed different ways of turning themselves inside out <cit.>, here, we focus on the so-called type-B inversion arising, for example, in Volvox globator <cit.>. This shares features such as invagination and involution with developmental events in metazoans <cit.>. This inversion scenario is distinct from type-A inversion, in which four lips open at the anterior of the shell and peel back to achieve inversion <cit.>. Type-B inversion begins with the appearance of a circular bend region at the equator of the embryo (Fig. <ref>c,d, Fig. <ref>a): cells there become wedge-shaped by developing narrow basal stalks <cit.>. At the same time, the cells move relative to the cytoplasmic bridges so as to be connected at their thin stalks, thus splaying the cells and bending and, eventually, invaginating the cell sheet <cit.>. <cit.> showed that inversion is arrested in the absence of analogous motion of cells relative to the cytoplasmic bridges in type-A inversion in Volvox carteri. The relative motion results from a kinesin associated to the microtubule cytoskeleton (Fig. <ref>, figure supplement 1); orthologues of this kinesin are found throughout the Volvocine algae <cit.>. After invagination, the posterior hemisphere moves into the anterior (Fig. <ref>e), the phialopore widens and the anterior hemisphere moves over the subjacent posterior (Fig. <ref>f) while `rolling' over a second circular bend region, the anterior cap <cit.>. Additional cell shape changes (Fig. <ref>d–f, Fig. <ref>b–d) in the anterior and posterior hemispheres are implicated in the relative contraction and expansion of either hemisphere with respect to the other <cit.>. This plethora of cell shape changes is possible as Volvox cells do not have a cell wall <cit.>.

In a previous study <cit.>, we combined light sheet microscopy and theory to analyse the early stages of inversion, showing that only a combination of active bending and active stretching (i.e. expansion or contraction) can account for the cell sheet deformations observed during invagination. The crucial role of active stretching was also highlighted by <cit.> who showed that type-A inversion in Volvox carteri cannot complete if acto-moysin mediated contraction is inhibited chemically. We later analysed the mechanics of this competition between bending and stretching in more detail <cit.>. Here, we analyse experimentally the variability of the shapes of inverting Volvox globator at consecutive stages of inversion. We refine our theoretical model to capture later stages of the inversion process, and finally combine theory and experiment to untangle the geometric, mechanical, and active contributions to the observed spatial structure of the shape variations.





§ RESULTS

We acquired three-dimensional time-lapse visualisations of inverting Volvox globator embryos (Video 1) using a selective-plane-illumination-microscopy setup (Methods) based on the OpenSPIM system <cit.>. Data were recorded for 13 parent spheroids containing, on average, 6 embryos. Summary statistics for 33 embryos were obtained from the recorded z-stacks and, for a more quantitative analysis of inversion, embryo outlines were traced on midsagittal sections of 11 of the recorded inversion processes, selected for optimal image quality (Methods).



In our previous work <cit.>, we discussed in detail three geometric descriptors of the traced embryo outlines, which we have reproduced for this dataset (Methods):

 
  * the distance e (Fig. <ref>a) from the posterior pole to the plane of the circular bend region; this serves as an indicator of the progress of the `upwards' movement of the posterior hemisphere;
 
  * the embryonic surface area A (Fig. <ref>b), which was computed by determining a surface of revolution from each half of the midsagittal slice and averaging the two values for each timepoint;
 
  * the minimal (most negative) value κ_∗ of the meridional curvature in the bend region (Fig. <ref>c).

We have computed three additional descriptors associated with the progress of later of inversion:

 
  * the diameter d of the phialopore (Fig. <ref>d) as an indicator of progress of inversion of the anterior hemisphere;
 
  * the width w of the bend region (Fig. <ref>e), where the bend region is defined as the region of negative curvature;
 
  * the position of the bend region (Fig. <ref>f), measured along the arclength of the deformed shell from the posterior pole to the midpoint of the bend region.

The computation of these descriptors is discussed in the Methods section. Each of these descriptors evolves in qualitatively similar ways in individual embryos, yet their evolution occurs over different timescales in different embryos and the local maxima in surface area (Fig. <ref>b) and in phialopore width (Fig. <ref>d) occurs at different relative times in different embryos. This initial impression of the variability of inversion is confirmed by the analysis of three summary statistics: (1) the duration of inversion, from appearance of the bend region to closure of the phialopore, (2) the diameter of the embryos post-inversion, (3) the relative time during inversion at which the phialopore starts to open. Histograms of these quantities in Fig. <ref>a–c reveal considerable variability, thus showing that the noise does not only affect the global duration of inversion, but also the relative timing of parts of it. We additionally note that there is no correlation between the size of an embryo and the duration of its inversion (Fig. <ref>d), not even between embryos from the same parent spheroid.



It is natural to ask to what extent the different deformations of inversion must arise in a particular order: while invagination occurs before phialopore opening in all our samples, analysis of characteristic `checkpoints' of inversion (Fig. <ref>, figure supplement 1) reveals that there is still considerable leeway in the timing of posterior inversion and phialopore opening. To further quantify the variability of inversion, we must define an average inversion sequence; our averaging approach must take into account these different types of variability.



 §.§ The Local Variability of Inversion

To define an average inversion sequence and analyse its mechanics, we compare the local geometry of the traced curves. The question of how to define an appropriate metric for this kind of comparison goes back at least to the work of D'Arcy Thompson <cit.>, and is altogether a rather philosophical one, to which there is no unique answer. Thompson showed for example how the outlines of fish of different species could be mapped onto one another by dilations, shears, and compositions thereof. In Volvox inversion, these shape differences are likely to arise from variations in cell shape and variations in the positions of cell shape changes. Our averaging method must therefore allow for these local variations as well as for differences in the timing of the cell shape changes (as suggested by the analysis of the summary statistics), while recognising that the posterior poles and the rims of the phialopores of the different embryos must correspond to each other. Our approach is therefore based on minimising the Euclidean distance between individual embryo shapes and their averages, with alignments obtained using dynamic time warping (Methods and Fig. <ref>, figure supplement 1). Results are shown in Fig. <ref>.

Averaging approaches that do not consider both stretching in time of individual inversions and local stretching of corresponding points of individual shapes tend to give unsatisfactory results: the simplest averaging approach is to align the inversion sequences by a single time point, say when the posterior-to-bend distance reaches half of its initial value (Methods and Fig. <ref>, figure supplement 2). The absence of time stretching, however, means that large variations arise at later stages of inversion. (Given the dramatic embryonic shape changes during inversion, it is not suprising that there should be no single parameter that could be used to align inversions of different embryos.) A better alignment is obtained if we allow stretching in time (Methods and Fig. <ref>, figure supplement 3), but this method, without local stretching of individual shapes relative to each other, produces unsatisfactory kinks in the bend region of the average shapes (Fig. <ref>, figure supplement 3).



The averages reveal that inversion seems to proceed at an approximately constant speed relative to the average inversion sequence (Fig. <ref>a,b). However, the alignment reveals that different stages of inversion take different times in different embryos (Fig. <ref>a), with some embryos seeming to linger in certain stages. This is the same non-linearity that we already saw earlier on the timelines in Fig. <ref>, figure supplement 1 obtained from the measurements in Fig. <ref>.
 
To analyse the local variations of the embryo shapes, we define, at each point of the average shapes, a covariance ellipse. The curves that are parallel to the average shape and tangent to the covariance ellipse define what we shall term the standard deviation shape. These standard deviation shapes measure the variability of the average shapes and are shown in Fig. <ref>. The variations they represent naturally divide into two components: first, those variations that are parallel to the average shape, and second, those perpendicular to the average shape. The former represent mere local stretches of the average shapes, while the latter correspond to actual variations of the shapes; we shall therefore refer to the thickness of the standard deviation shapes as `shape variation' in what follows. We report the mean shape variation and its standard error in Fig. <ref>c. This plot shows that the mean shape variation reaches a maximal value around the stages in Fig. <ref>g–i: different embryos start from the same shape and reach the same inverted shape after inversion (up to a scaling), but may take different inversion paths. Plotting the mean shape variation for different averaging methods (Fig. <ref>, figure supplement 1), we confirm that the present averaging method yields a better alignment than the alternative methods discussed earlier. 





It is intriguing, however, to note the spatial structure of the local shape variations. In particular, during the early stages of posterior inversion (Fig. <ref>d–f), the shape variation is smaller in the active bend region than in the adjacent anterior cap (Fig. <ref>e, the second bend region of increased positive curvature). As the phialopore opens and the anterior begins to peel back over the partially inverted posterior (Fig. <ref>h) the relative shape variations becomes smaller in the anterior cap. The initially small variation in the bend region is especially intriguing since this is where cells become wedge-shaped to drive invagination, while the anterior cap bends passively <cit.>. In other words, the shape variation is reduced in the part of the cell sheet where the active cell shape changes driving inversion arise. This correspondence characterises what one might term, from a teleological point of view, a `good' inversion. We shall focus on a less exalted question, the answer to which will be falsifiable, however: how is this spatial structure of the variability related to the mechanics of inversion? Before addressing this question, we need to analyse the mechanics of inversion in some more detail.



 §.§ Active Bending and Stretching during Inversion

Which active deformations are required for inversion? In our previous work <cit.>, we addressed 
this question for the early stages of inversion: at a mechanical level of description, invagination arises from an interplay of active bending and stretching <cit.> associated with different types of cell shape changes . A key role is played by the cells close to the equator of the cell sheet (Fig. <ref>d, Fig. <ref>a), which become wedge-shaped <cit.>, thus splaying the cells and hence imparting intrinsic curvature to the cell sheet <cit.>. Yet no such cell wedging has been reported at the anterior cap at later stages of inversion, when the anterior hemisphere peels back over the partly inverted posterior (Fig. <ref>f, Fig. <ref>d). 

To resolve this conundrum, we ask whether the additional cell shape changes observed during type-B inversion <cit.> are sufficient to explain anterior peeling: cells in the anterior hemisphere have flattened, ellipsoidal shapes, while cells on the posterior side of the anterior cap are pencil-shaped (Fig. <ref>e, Fig. <ref>d). We have previously described the early stages of inversion using a mathematical model <cit.> in which cell shape changes appear as local variations of the intrinsic (meridional and circumferential) curvatures κ_s^0,κ_ϕ^0 and stretches f_s^0,f_ϕ^0 of an elastic shell. We recall the difference between open, one-dimensional elastic filaments and two-dimensional elastic shells in this context: the former can simply adopt a shape in which the curvature and stretch are everywhere equal to their intrinsic values. For the latter, by contrast, the intrinsic curvatures and stretches may not be compatible with the global geometry, causing the shell to deform elastically and adopt actual (meridional and circumferential) curvatures κ_s,κ_ϕ and stretches f_s,f_ϕ different from the imposed intrinsic curvatures and stretches. In order to address these later stages of inversion, we must first generalise our previous mathematical model using ideas from morphoelasticity (Methods).
Indeed, that model was derived under the assumption of small strains. While the elastic strains are small indeed (since the metric tensor, which describes the deformed shape, is close to the intrinsic tensor defined by the cell shape changes), the geometric strains are large: both the metric tensor of the deformed shell and the intrinsic tensor differ considerably from the metric tensor of the undeformed sphere.



To test whether anterior peeling can be achieved by contraction of the cell sheet alone, we impose functional forms for the intrinsic stretches of the shell (Fig. <ref>a–c) representing these cell shape changes, but we do not modify the intrinsic curvatures in the anterior hemisphere (Fig. <ref>c). In particular, the linear variation of the circumferential stretch in the anterior hemisphere represents the different orientations of the ellipsoidal cells at the phialopore (Fig. <ref>c), where the long axis is the circumferential axis, and at the anterior cap, where the long axis is the meridional axis <cit.>. In our quasi-static simulation, we approximate the shape in Fig. <ref>h by a configuration with inverted posterior hemisphere (Fig. <ref>d), and displace the intrinsic `peeling front' (Fig. <ref>a,b). The shell responds by peeling (Fig. <ref>e), with shapes in qualitative agreement with the experimentally observed shapes. Since the peeling front is located at the anterior cap, where the shape variation is reduced during anterior peeling as discussed previously, we again see a correlation between reduced shape variations and the location of the active cell shape changes driving inversion. 

These considerations suggest that contraction is sufficient to drive the peeling stage of inversion, even without changes in intrinsic curvature. Although the position of the cytoplasmic bridges (Fig. <ref>d), on the inside end of the cells at the end of inversion <cit.>, suggests that the intrinsic curvature may change sign in the anterior hemisphere, too, this appears to be a seconday effect. Hence intrinsic bending complements intrinsic stretching. By contrast, our previous work <cit.> revealed that stretching complements bending during invagination. The roles of stretching and bending are thus interchanged during inversion of the posterior and anterior hemispheres, and the embryo uses these two different deformation modes for different tasks during inversion.



  §.§.§ Analysis of Cell Shape Changes

For a more quantitative analysis of the data and to validate our model, we proceed to fit the elastic model to the experimental average shapes (Methods). In the model, we impose a larger extent of the phialopore than in the biological system, where the phialopore is initially very small (Fig. <ref>a). This is an important simplification to deal with the discrete nature of the few cells that meet up at the phialopore. Nonetheless, using fifteen fitting parameters to represent previously observed cell shape changes <cit.> in terms of the intrinsic stretches and curvatures (Methods), the model captures the various stages of inversion (Fig. <ref>). This supports our interpretation of the observed cell shape changes (Fig. <ref>) and their functions. Comparing the geometric descriptors discussed previously (Fig. <ref>) for the experimental averages and the fitted shapes (Fig. <ref>, figure supplement 1), we notice that the fitted shapes underestimate the width of the bend region. Because curvature is a second derivative of shape, it is not surprising that larger differences arise in the minimal bend region curvature of the average and fitted shapes (Fig. <ref>, figure supplement 1).



Nonetheless, the fitted values of the intrinsic curvature of the cell sheet also resolve a cell shape conundrum: during invagination, the curvature in the bend region increases (Fig. <ref>), yet <cit.> reported similar wedge-shaped cells in the bend region at early and late invagination stages, although the number of wedge-shaped cells in the bend region increases as invagination progresses <cit.>. The fitted parameters indeed suggest a constant value of the intrinsic curvature at early stages of inversion, while the actual curvature in the bend region increases (Fig. <ref>a). This serves to illustrate that the intrinsic parameters cannot simply be read off the deformed shapes and confirms that there is but a single type of cell change, expanding in a wave to encompass more cells, and thus driving invagination. It is only at later stages of inversion, when the wedge-shaped cells in the bend region become pencil-shaped <cit.>, that both the intrinsic curvature and the actual curvature in the bend region decrease (Fig. <ref>a). The fitted shapes also yield the posterior and anterior limits of the bend region (Fig. <ref>b), i.e. the original positions, relative to the undeformed sphere, of the corresponding cells. Because of the varying spatial stretches of the shell, these positions cannot simply be read off the deformed shapes, but must be inferred from the fits. The fitted data suggest that invagination results from an intrinsic bend region of constant width, complemented by other cell shape changes (Fig. <ref>d, Fig. <ref>b,c). The 
region of wedge-shaped cells (and, by implication, of negative intrinsic curvature) starts to expand into the posterior at constant speed (i.e. at a constant number of cell shape changes per unit of time) between the stages in Fig. <ref>e,f. Anterior inversion starts about five minutes later when this region begins to expand into the anterior just after the stage in Fig. <ref>g. 



Fig. <ref> also shows the stretches f_s,f_ϕ in the fitted shapes. It is particularly interesting to relate the values of f_s,f_ϕ in the fitted shapes to the measurements of individual cells by <cit.>: before inversion starts, the cells are teardrop-shaped, and measure 3-5 µ m in the plane of the cell sheet. As invagination starts, the cells in the posterior hemisphere become spindle-shaped, measuring 2-3 µ m. This suggests values f_s,f_ϕ≈ 0.6-0.66 in the posterior hemisphere during invagination, in agreement with the fitted data (Fig. <ref>d). At later stages of inversion, the cells in the bend region become pencil-shaped, measuring 1.5-2 µ m in the meridional direction, suggesting smaller values f_s≈ 0.4-0.5 there, again in agreement with the fitted data (Fig. <ref>h). The large stretches f_s>2 seen in the anterior cap during inversion of the posterior hemisphere (Fig. <ref>f) cannot be accounted for by the disc-shaped cells in the anterior (which only measure 4-6 µ m) in the meridional direction). While examination of the thin sections of <cit.> does suggest, in qualitative agreement with the fits, that the largest meridional stretches arise in the anterior cap, the fact that the model overestimates the actual values of these stretches may stem from the simplified modelling of the phialopore. Further, at the very latest stages of inversion (Fig. <ref>j), the fitted shapes suggest very small values f_s<0.3 and corresponding values f_ϕ>3 that are not borne out by the cell measurements.



  §.§.§ Phialopore Opening and Cell Rearrangement

To understand how these values of the stretches at odds with the observed cell shape changes arise in the fitted shapes, we must analyse the opening of the phialopore in more detail. The observations of <cit.> show that the cytoplasmic bridges stretch considerably, to many times their initial length, as the phialopore opens. Circumferential elongation of cells as a means to increase effective radius was discussed in some detail by <cit.>, but is not sufficient to explain the circumferential stretches observed at the phialopore. Additional elongation of cytoplasmic bridges as a means to further increase the effective radius (Fig. <ref>) may suffice to produce the large circumferential stretches, but does not explain the small values of meridional stretch at the phialopore in the fitted shapes. For this reason, we additionally imaged the opening of the phialopore using confocal laser scanning microscopy (Methods) to resolve single cells close to the phialopore (Video 2). 

L6.3cm

    < g r a p h i c s >

Mechanisms of phialopore stretching: cell shape changes, stretching of cytoplasmic bridges, and cell rearrangements. Red lines represent cytoplasmic bridges; fainter colours signify other, out-of-plane cells. CSC: cell shape changes, CB: cytoplasmic bridge.
 


The data reveal that cells rearrange near the phialopore, suggesting an additional mechanism to stretch the phialopore sufficiently for the anterior to be able to peel over the inverted posterior (Fig. <ref>). Video 2 shows how, initially, only a small number of cells form a ring at the anterior pole. When the phialopore widens, cells that were initially located away from this initial ring come to be positioned at the rim of the phialopore. It is unclear whether the cytoplasmic bridges between these cells stretch or break, or whether these cells were not connected by cytoplasmic bridges in the first place. While such cell rearrangement is beyond the scope of the current model, it is nevertheless captured qualitatively by the small values of f_s near the phialopore. <cit.> observed elongation of cytoplasmic bridges near the phialopore of Volvox aureus, but not in small fragments of broken-up embryos, and concluded that the elongation of cytoplasmic bridges was the result of passive mechanical forces. By contrast, in our model, the opening of the phialopore is the result of active cell shape changes there. This discrepancy may herald a breakdown of the approximations made to represent the phialopore.
The data also hint that there may be a different mechanical contribution at later stages of inversion (Fig. <ref>i), where the rim of the phialopore may be in contact with the inverted posterior. Since the model does not resolve the rim of the phialpore in the first place, we do not pursue this further here. For completeness of the mechanical analysis, we analyse such a contact configuration in Appendix <ref>, where we also discuss a toy problem to highlight the intricate interplay of mechanics and geometry in the contact configuration. 





 §.§ Mechanics and Regulation of Local Shape Variations

We now return to the spatial structure of the shape variations discussed previously. It is clear that some of this structure is geometrical: since the shapes are aligned so that the positions of their centres of mass along the axis coincide, the shape variations accumulate, and are thus expected to e.g. increase in the anterior hemisphere, towards the phialopore, as at the stage in Fig. <ref>c. At the same stage however, the shape variation is smaller in the bend region than in the adjacent anterior cap. Both of these regions are, however, close to the centre of mass, and so we do not expect this difference to arise from mere geometric accumulation of shape variations. We must therefore ask: can this structure arise purely mechanically (i.e. from a uniform distribution of the intrinsic parameters), but possibly as a statistical fluke, or must there be some regulation (i.e. non-uniform variation of the intrinsic parameters)?

To answer this question, we analyse random perturbations of the fitted intrinsic parameters of the inversion stage in Fig. <ref>c. We observe that, if the relative size of perturbations (the `noise level') exceeds about 4% at this stage of inversion, computation of the perturbed shapes fails for some parameter choices. This mechanical effect is not surprising: our previous analysis of invagination <cit.> revealed strong shape non-linearities and the possibility of bifurcations as the magnitude of the intrinsic curvature in the bend region is increased. While we may therefore expect more leeway in some parameters than in others, we shall simply discard those perturbations for which the computation fails; further estimation of the distribution of possible perturbations is beyond the scope of the present discussion. We now estimate, for each noise level, the mean shape variation from 1000 perturbations of the fitted shape. By comparing this to the mean shape variation estimated from the N=22 embryo halves in Fig. <ref>c, we roughly estimate a noise level of 7.5% (Fig. <ref>a). At this noise level, about 15% of perturbations fail; while the non-uniformities are small, they are statistically significant (Methods).

With this noise level, we obtain 10000 samples of N=22 perturbations to the fitted shape each (Fig. <ref>b), and we compute their averages in the same way as for the experimental samples. While these samples qualitatively capture the spatial structure of the shape variation, they overestimate the shape variation at the poles. More strikingly, they feature a local maximum of the shape variation in the bend region, rather than in the anterior cap. From the sample distribution of the position of these local maxima (Fig. <ref>c), it is clear that the experimental distribution with the local maximum in the anterior cap, is very unlikely to arise under this model. We make this statement more precise statistically in the Methods section. To explain the observed structure of the shape variation, we therefore allow more variability in the meridional stretch in the anterior cap (with a noise level of 80%, compared to 2.5% for the remaining parameters to reproduce the mean shape variation). The resulting distribution is consistent with the experimentally observed position of the local maximum of shape variation in the anterior cap (Fig. <ref>b,c). While still overestimating the variability near the posterior pole, this modified distribution of the parameter variability captures the magnitude of the variability in the anterior cap much better than the original one. 

Thus, at this early stage of inversion (Fig. <ref>c), the observed embryo shapes are consistent with an increased variation of the intrinsic meridional stretch in the anterior cap. We can take the interpretation of this active regulation (or lack thereof) further by relating it to the observed cell shape changes: at the stage of Fig. <ref>c, the variations of the meridional stretch in the anterior cap correspond to the formation of disc-shaped cells there (Fig. <ref>d). This indicates that invagination and initiation of the expansion of the anterior hemisphere (via the formation of disc-shaped cells) are really two separate processes the relative timing of which is not crucial. (The formation of disc-shaped cells starting at different times also explains the large noise level in the meridional stretch under the modified model, although there is no fitting involved, here.) This adds to our earlier point, that these processes rely on different deformation modes (active bending for invagination and active contraction and stretching for inversion of the anterior hemisphere). These considerations also rationalise our second observation concerning the spatial structure of shape variations, that the variation in the anterior cap is reduced as inversion of the posterior hemisphere ends (Fig. <ref>h): there are no longer two separate processes at work. We finally point to a purely mechanical aspect of the structure of the shape variations: despite the increased variability in the anterior cap, the mechanics ensure that the variability is lowest in the bend region, where the main cell shape changes driving invagination take place.





§ DISCUSSION

In this paper, we have combined experiment and theory to analyse the variability of Volvox inversion and obtain a detailed mechanical description of this process. From observations of the structure of the variability of the shapes of inverting Volvox embryos, we showed, using our mathematical model, that this structure results from a combination of geometry, mechanics, and active regulation. The simplest scenario with which the observed shape variations are consistent is that type-B inversion in Volvox globator results from two separate processes, with most of the variability at the invagination stage attributed to the relative timing of these processes in individual embryos. The difference between these processes is mirrored, at a mechanical level, by the different types of deformations driving them: the first process, to invert the posterior hemisphere, mainly relies on active bending, whereas the second process, to invert the anterior hemisphere, is mainly driven by active expansion and contraction. We anticipate that these ideas and methods can be applied to other morphogenetic events in other model organisms to add to our understanding of the regulation of morphogenesis: what amount of regulation, be it spatial or temporal, of the cell-level processes is there, and how does it relate to the amount required mechanically for the processes to be able to complete? Additionally, <cit.> showed that diffusion of two morphogens with inhibition à la <cit.> has error-correcting properties that can explain the precise domain specification that is observed in Drosophila embryos in spite of the huge variability of morphogen gradients <cit.>. Does the interplay of geometry and mechanics yield analogous error-correcting properties?

While we have begun to analyse the mechanical regulation of development in the context of Volvox inversion, our answers this far have been either negative (excluding certain mechanisms of regulation) or of what one might term the Occam's razor variety (invoking the law of parsimony to find the simplest modification of the model that can explain the observations). This approach of testing falsifiable hypotheses mitigates the risk of drawing conclusions that are mere teleology <cit.>. Nonetheless, a fuller answer to the questions above requires estimation of the variability of the model parameters from the experimental data, yet that endeavour entails significant statistical, computational, and experimental difficulties: to estimate the variability with statistical signficance we need a large number of experimental samples to estimate the experimental distribution; for each step of the optimisation algorithm used to estimate the large number of variability parameters, a large number of computational samples must be computed to estimate the distribution under the model. Similar difficulties arise when estimating the variability allowed mechanically. While we have previously noted <cit.> that the dynamic data for type-B inversion suggest that invagination proceeds without a `snap-through' bifurcation, there is no general requirement for individual developmental paths to lie on one and the same side of a mechanical bifurcation boundary. This poses an additional challenge for modelling approaches.

After this discussion of general challenges for a mechanobiological analysis of morphogenesis and its regulation, we mention some of the remaining questions specific to Volvox inversion: our model does not resolve the details of the phialopore, and hence does not describe the closure of the phialopore at the end of inversion, which remains a combined challenge for experiment and theory: as discussed above,  the cytoplasmic bridges elongate drastically at the phialopore <cit.>, and confocal imaging has revealed the possibility of rearrangements within the cell sheet at the phialopore. Do some cytoplasmic bridges rend to make such rearrangements possible? Understanding the details of the opening of the phialopore may also require answering a more fundamental question the answer to which has remained elusive <cit.>: what subcellular structures are located within the cytoplasmic bridges and how is it possible for them to stretch to such an extent? At the theoretical level, rearrangements of cells near the phialopore raise more fundamental questions of morphoelasticity <cit.>: in particular, how does one describe the evolution of the boundary of the manifold underlying the elastic description? Cytoplasmic bridges rending next to the phialopore would lead to the formation of lips similar to those seen in type-A inversion <cit.>. Is there a simple theory to describe the elasticity of this non-axisymmetric setup?

At the close of this discussion, it is meet to briefly dwell on a question of more evolutionary flavour: how did different species of Volvox evolve different ways of turning themselves inside out? Mapping inversion types to a phylogenetic tree of Volvocine algae shows shows that different inversion types evolved several times independently in different lineages <cit.>. Additionally, <cit.> reported that in Volvox rousseletii and Volvox capensis, inversion type depends on the (sexual or asexual) reproduction mode. This may be a manifestation of the poorly understood role of environmental and evolutionary cues in morphogenesis <cit.>, but it is natural to wonder whether there is a mechanical side to this issue. Ultimately, this is another incentive to study the mechanics of type-A inversion in more detail.  



§ METHODS AND MATERIALS




 §.§ Acquisition of Experimental Data

Wild-type strain Volvox globator Linné (SAG 199.80) was obtained from the Culture Collection of Algae at the University of Göttingen, Germany <cit.>, and cultured as previously described <cit.> with a cycle of 16 h light at 24°C and 8 h dark at 22°C.



  §.§.§ OpenSPIM Imaging

A selective plane illumination microscope (SPIM) was assembled based on the OpenSPIM setup <cit.>, with modifications to accommodate a Stradus® Versalase™ laser system with multiple wavelengths (Vortran Laser Technology, Inc., Sacramento, CA, USA) and a CoolSnap Myo CCD camera (1940× 1460 pixels; Photometrics, AZ, USA). Moreover, to decrease the loss of data due to shadowing a second illumination arm was added to the setup (Fig. <ref>). Illumination from both sides improved the image quality and enabled re-slicing of the z-stacks when embryos began to spin during anterior inversion. 

Volvox globator parent spheroids were mounted in a column of low-melting-point agarose and suspended in fluid medium in the sample chamber. To visualise the cell sheet deformations of inverting Volvox globator embryos, chlorophyll-autofluorescence was excited at λ=561 nm and detected at λ>570 nm. Z-stacks were recorded at intervals of 60 s over 4-6 hours to capture inversion of all embryos in a parent spheroid. We acquired time-lapse data of 13 different parent spheroids each containing 4-7 embryos. 

L.46

    < g r a p h i c s >

SPIM imaging setup. a: beamsplitter cube, b: mirror, c: beam expander, d: cylindrical lens, e: telescope, f: illumination objective, g: detection objective, h: emission filter, i: camera.





  §.§.§ Confocal Laser Scanning Microscopy

Samples were immobilised on glass-bottom dishes by embedding them in low-melting-point agarose and covered with fluid medium. Chlorophyll-autofluorescence was excited at  and detected at λ>647 nm. Z-stacks were recorded at intervals of 30 s over 1-2 hours to capture inversion of a single embryo. Trajectories of individual cells close to the phialopore were obtained using Fiji <cit.>. Experiments were carried out using a Observer Z1 spinning-disk microscope (Zeiss, Germany).



  §.§.§ Image Tracing

To ensure optimal image quality (traceability) for the quantitative analyses of inversion, from the inversion processes recorded with the SPIM, we selected 11 inversions (in 6 different parent spheroids) in which the acquisition plane was initially approximately parallel to the midsagittal plane of the embryos. Midsagittal cross-sections were obtained using Fiji <cit.> and Amira (FEI, OR, USA).

Splines were fitted to these cross-sections using the following semi-automated approach implemented in Python/C++: in a preprocessing step, images were bandpass-filtered to remove short-range noise and large-range intensity correlations. Low-variance Gaussian filters were applied to smooth out the images slightly. Splines were obtained from the pre-processed images I(x⃗) using the active contour model <cit.>, with modifications to deal with intensity variations and noise in the image: the spline x⃗_⃗s⃗(s), where s is arclength, minimises an energy

    ℰ[x⃗_s]=ℰ_image[x⃗_s]+ℰ_spline[x⃗_s]+ℰ_skel[x⃗_s],

where 


    ℰ_image[x⃗_s]    =-α∫I(x⃗_⃗s⃗(s)) ds,
    ℰ_spline[x⃗_s]    = β∫∂^2x⃗_⃗s⃗∂ s^2^2ds+γ(∫ds-L_0)^2,
    ℰ_skel[x⃗_s]   =δ∫I_skel(x⃗_⃗s⃗(s)) ds,


wherein α,β,γ,δ are parameters, L_0 is the estimated length of the shape outline, and I_skel is obtained by skeletonising I using the algorithm of <cit.> to minimise the number of branches. 

The energy ℰ was minimised using stochastic gradient descent. Initial guesses for the splines were obtained by manually initialising about 15 timepoints for each inversion using a few guidepoints and polynomial interpolation. An initial guess for other frames was obtained from these frames by interpolation; these interpolated shapes were used to estimate L_0.

With δ = 0, the standard active contour model of <cit.> is recovered. We found that this model was not sufficient to yield fits of sufficient quality, because of the existence of local minima at small values of α, while larger values of α lead to noisy splines. Thresholding methods on their own were not sufficient either, because of branching and, in particular, since they failed to capture the bend region properly.  Dynamic thresholding methods <cit.> are not applicable either because of the fast variations of the brightness of the images. The modified active contour model did however produce good fits when we progressively reduced δ to zero with increasing iteration number of the minimisation scheme, yielding smooth splines, while overcoming the local minima (or, from the point of view of the skeletonisation method, choosing the correct, branchless part of the skeleton). All outlines obtained from this algorithm were manually checked and corrected.



 §.§ Analysis of Traced Embryo Shapes

From the traced cell sheet outlines, anterior-posterior axes of the embryos were determined as follows: for shapes for which the bend region was visible on either side of the cross-section, the embryo axis was defined to be the line through the centre of mass of the shape that is perpendicular to the common tangent to the two bend regions (the apex line). Shapes were then rotated and translated manually so that their axes coincided. Since embryos do not rotate much before the flagella grow, the orientation of the axes of the earliest traces (for which the bend regions are not apparent) were taken to be the same as that of the earliest timepoint for which two bend regions were visible. The intersection of the embryo trace and axis defines the posterior pole. After manually recentring some embryos with more pronounced asymmetry, embryos were halved to obtain N=22 embryo halves. 



  §.§.§ Computation of Inversion Descriptors

From the aligned shapes, the geometric descriptors of inversion reported in Fig. <ref> were computed as follows: the posterior-to-bend distance e was computed as the distance from the apex line to the posterior pole. The maximal surface area A_max and the most negative value of curvature κ_∗ in the bend region were computed as described previously <cit.>; traces were smoothed before computing the curvature. The phialopore width d was computed as the absolute distance between the two ends of a complete embryo trace. The bend region was defined as the region of negative curvature; the distance between the first and last points of negative curvature defined the bend region width. The bend region position is defined by the distance, along the embryo trace, between the posterior pole and the midpoint of the bend region. (The latter may differ from the point where the most negative value of curvature is attained.) The values of bend region width and position obtained for each embryo half were averaged to yield the reported values.



  §.§.§ Aligning and Averaging Embryo Shapes

To align embryos to each other, one embryo half was arbitrarily taken as the reference shape, and T=10 regularly spaced timepoints were chosen for fitting. (These timepoints were chosen to be well after invagination had started and before the phialopore had closed, so that defining the start and end of inversion was not required.) For each of the remaining N-1 embryo halves, a scale and T corresponding timepoints were then sought, with shapes being (linearly) interpolated at intermediate timepoints. The interpolated and scaled shapes were centred so that the centres of mass of the cross-sections coincided. This fixes the degree of freedom of translation parallel to the embryo axis; the position perpendicular to the axis is fixed by requiring that the embryo axes coincide (Fig. <ref>, figure supplement 1a). The motivation for using the centres of mass of the cross-sections (rather than of that of the embryos, which assigns the same mass to each cell by assigning more mass to those points of the cross-section that are farther away from the embryo axis) is a biological one: because of the cylindrical symmetry of the cell shape changes, this average assigns the same mass to each cell shape change. 

For aligning embryo shapes, we distribute M=100 averaging points uniformly along the (possibly different) arclength of each of the embryo halves. Corresponding points were determined using dynamic time warping (DTW) as described by e.g. <cit.>, and the distances between these shapes and their averages were minimised as explained in what follows. The parameters describing the alignment are thus the scale factors S_1 = 1, S_2,… S_N and the averaging time points τ⃗_1 = (τ_11,τ_12,…,τ_1T),τ⃗_2,…,τ⃗_N, where τ⃗_1 is fixed. Each choice of these parameters yields a set of shapes X⃗_1=(x_11,…,x_1M),X⃗_2,…,X⃗_N with points matched up by maps σ_1,σ_2,…,σ_N obtained from the DTW algorithm. The effect of the local stretching allowed by the DTW algorithm is illustrated in Fig. <ref>, figure supplement 1b,c. The mean shapes having been determined, the sum of Euclidean distances between shapes of individual embryos and the mean,

    ∑_t=1^T{∑_n=1^N∑_m=1^M(x_nσ_n(m)-x_m)^2}^1/2,    x_m=1N∑_n=1^Nx_nσ_n(m),

was minimised over the space of all these alignment parameters using the Matlab (The MathWorks, Inc.) routine , modified to incorporate the variant of the Nelder–Mead algorithm suggested by <cit.> for problems with a large number of parameters. After the algorithm had converged, each of the alignment parameters was modified randomly, and the algorithm was run again. This was repeated until the alignment score defined by (<ref>) did not decrease further. The means x_1,x_2,…,x_M for the alignment minimising (<ref>) define the average embryo shapes.

Aligning shapes in this way using dynamic time warping requires a considerable amount of computer time. To make the problem computationally tractable, we invoked the usual heuristics of only computing pairwise DTW distances, and reducing the size of the DTW matrix by only computing a band centred on the diagonal. To verify the algorithm, we also ran several instantiations of the alignment algorithm without DTW (i.e. with σ_n=id) and with larger parameter randomisations, confirming that the modified Nelder–Mead algorithm finds an appropriate alignment. This also enabled us to verify that results do not change qualitatively if the centres of mass of the cross-sections are replaced with those of the embryo halves (even though, as noted in the main text, the shapes without DTW are unsatisfactory since they have kinks in the bend region that are not seen in individual embryo shapes).

For the simple alternative averaging method in Fig. <ref>, figure supplement 2, different numbers of averaging points were distributed at equal arclength spacing along all individual shapes. Differences in arclengths of individual embryos mean that the rims of the phialopores of individual embryo halves are not necessarily matched up (Fig. <ref>, figure supplement 1c). No time stretching was applied. The averaging method in Fig. <ref>, figure supplement 3, is the method discussed above, without DTW (i.e. with σ_n=id).



 §.§ Elastic Model




We consider a spherical shell of radius R and uniform thickness h≪ R (Fig. <ref>a), characterised by its arclength s and distance from the axis of revolution ρ(s), to which correspond arclength S(s) and distance from the axis of revolution r(s) in the axisymmetric deformed configuration (Fig. <ref>b). We define the meridional and circumferential stretches

    f_s(s)=dSds,       f_ϕ(s)=r(s)ρ(s).

The position vector of a point on the midsurface of the deformed shell is thus

    r⃗(s,ϕ)=r(s)u⃗_⃗r⃗(ϕ)+z(s)u⃗_⃗z⃗,

in a right-handed set of axes (u⃗_⃗r⃗,u⃗_⃗ϕ⃗,u⃗_⃗z⃗) and so the tangent vectors to the deformed midsurface are

    e⃗_⃗s⃗   =r'u⃗_⃗r⃗+z'u⃗_⃗z⃗,      e⃗_⃗ϕ⃗=ru⃗_⃗ϕ⃗,

where dashes denote differentiation with respect to s. By definition, r'^2+z'^2=f_s^2, and so we may write 

    r'=f_scosβ,       z'=f_ssinβ.

Hence the normal to the deformed midsurface is

    n⃗= r'u⃗_⃗z⃗-z'u⃗_⃗r⃗f_s=cosβ u⃗_⃗z⃗-sinβ u⃗_⃗r⃗.

We now make the Kirchhoff `hypothesis' <cit.>, that the normals to the undeformed midsurface remain normal to the deformed midsurface (Fig. <ref>c). Taking a coordinate ζ across the thickness h of the undeformed shell, the position vector of a general point in the shell is

    r⃗(s,ϕ,ζ)=ru⃗_⃗r⃗+zu⃗_⃗z⃗+ζn⃗=(r-ζsinβ)u⃗_⃗r⃗+(z+ζcosβ)u⃗_⃗z⃗.

The tangent vectors to the shell are thus

    e⃗_⃗s⃗   = f_s(1-κ_sζ)(cosβ u⃗_⃗r⃗+sinβ u⃗_⃗z⃗),   e⃗_⃗ϕ⃗   = ρ f_ϕ(1-κ_ϕζ)u⃗_⃗ϕ⃗,

where κ_s=β'/f_s and κ_ϕ=sinβ/r are the curvatures of the deformed midsurface. The metric of the deformed shell under the Kirchhoff hypothesis accordingly takes the form

    dr^2 = f_s^2(1-κ_sζ)^2ds^2+f_ϕ^2(1-κ_ϕζ)^2ρ^2dϕ^2.

The geometric and intrinsic deformation gradient tensors are thus

    F^g=([ f_s(1-κ_sζ)           0;           0 f_ϕ(1-κ_ϕζ) ]),       F^0=([ f_s^0(1-κ_s^0ζ)               0;               0 f_ϕ^0(1-κ_ϕ^0ζ) ]),

where f_s^0,f_ϕ^0 and κ_s^0,κ_ϕ^0 are the intrinsic stretches and curvatures of the shell. Thence, invoking the standard multiplicative decomposition of morphoelasticity <cit.>, the elastic deformation gradient tensor is

    F=F^g(F^0)^-1=([ f_s(1-κ_sζ)f_s^0(1-κ_s^0ζ)                          0;                          0 f_ϕ(1-κ_ϕζ)f_ϕ^0(1-κ_ϕ^0ζ) ]).

While we do not make any assumption about the geometric or intrinsic strains derived from F^g and F^0, respectively, we assume that the elastic strains derived from F remain small; we may thus approximate 

    ε_ss≈f_s(1-κ_sζ)f_s^0(1-κ_s^0ζ)-1,       ε_ϕϕ≈f_s(1-κ_sζ)f_s^0(1-κ_s^0ζ)-1,

with the off-diagonal elements vanishing, ε_sϕ=ε_ϕ s=0. For a Hookean material with elastic modulus E and Poisson's ratio ν <cit.>, the elastic energy density (per unit extent in the meridional direction) is found by integrating across the thickness of the shell:

    ℰ2πρ   =E2(1-ν^2)∫_-h/2^h/2(ε_ss^2+ε_ϕϕ^2+2νε_ssε_ϕϕ)dζ
       =Eh2(1-ν^2){(1+h^24κ_s^0^2)E_s^2+(1+h^24κ_ϕ^0^2)E_ϕ^2+2ν(1+h^212(κ_s^0^2+κ_s^0κ_ϕ^0+κ_ϕ^0^2))E_sE_ϕ}
           +Eh^324(1-ν^2){K_s^2+K_ϕ^2+2ν K_sK_ϕ-4κ_s^0E_sK_s-4κ_ϕ^0E_ϕ K_ϕ-2ν(κ_s^0+κ_ϕ^0)(E_ϕ K_s+E_sK_ϕ)},

where we have expanded the energy up to third order in the thickness, and where we have defined the shell strains and curvature strains

    E_s=f_s-f_s^0f_s^0,      E_ϕ=f_ϕ-f_ϕ^0f_ϕ^0,        K_s=f_sκ_s-f_s^0κ_s^0f_s^0,      K_ϕ=f_ϕκ_ϕ-f_ϕ^0κ_ϕ^0f_ϕ^0.

As in our previous work <cit.>, the elastic modulus is an overall constant that ensures that ℰ has units of energy, but does not otherwise affect the shapes. We shall also assume that ν=1/2 for incompressible biological material; the cell size measurements of <cit.> for type-A inversion in Volvox carteri support this assumption qualitatively. (These considerations also explain why we do not perturb these mechanical parameters in our analysis of the shape variations.) We finally set h/R=0.15 as in our previous work.



  §.§.§ Derivation of the Governing Equations

The derivation of the governing equations proceeds similarly to standard shell theories <cit.>. In fact, the resulting equations turn out to have a form very similar to those of standard shell theories, but a host of extra terms arise in the expressions for the shell stresses and moments due to the assumptions of morphoelasticity. The variation of the elastic energy takes the form

    δℰ2πρ=n_sδ E_s+n_ϕδ E_ϕ+m_sδ K_ϕ+m_ϕδ K_ϕ,

with

    δ E_s   =δ f_sf_s^0=1f_s^0(β δ r'+f_stanβ δβ),   δ E_ϕ   =δ f_ϕf_ϕ^0=δ rf_ϕ^0ρ,
    δ K_s   =δ(f_sκ_s)f_s^0=δβ'f_s^0,   δ K_ϕ   =δ(f_ϕκ_ϕ)f_ϕ^0=cosβf_ϕ^0ρδβ,


wherein dashes again denote differentiation with respect to s, and where the shell stresses and moments are defined by

    n_s   =Eh1-ν^2{E_s+ν E_ϕ+h^212(3κ_s^0^2E_s+ν(κ_s^0^2+κ_s^0κ_ϕ^0+κ_ϕ^0^2)E_ϕ-2κ_s^0K_s-ν(κ_s^0+κ_ϕ^0)K_ϕ)},
    
    n_ϕ   =Eh1-ν^2{E_ϕ+ν E_s+h^212(3κ_ϕ^0^2E_ϕ+ν(κ_s^0^2+κ_s^0κ_ϕ^0+κ_ϕ^0^2) E_s-2κ_ϕ^0K_ϕ-ν(κ_s^0+κ_ϕ^0)K_s)},


and


    m_s   =Eh^312(1-ν^2){K_s+ν K_ϕ-2κ_s^0 E_s-ν(κ_s^0+κ_ϕ^0)E_ϕ},
    
    m_ϕ   =Eh^312(1-ν^2){K_ϕ+ν K_s-2κ_ϕ^0 E_ϕ-ν(κ_s^0+κ_ϕ^0)E_s}.


Defining

    N_s=n_sf_s^0f_ϕ,       N_ϕ=n_ϕf_ϕ^0f_s,      M_s=m_sf_s^0f_ϕ,       M_ϕ=m_ϕf_ϕ^0f_s,

the variation becomes

    δℰ2π   = r N_sβ δ r+r M_s δβ
       -∫{(dds(r N_sβ)-f_sN_ϕ)δ r - (r f_sN_stanβ+f_sM_ϕcosβ-dds(r M_s))δβ} ds.

The Euler–Lagrange equations of (<ref>) are thus

    dds(rN_sβ)-f_sN_ϕ=0,       dds(rM_s)-f_sM_ϕcosβ-rf_sN_stanβ=0.

To remove the singularity that arises in the second of (<ref>) when β=π/2, we define the transverse shear tension T=-N_stanβ as in standard shell theories. The governing equations can then be rearranged to give

    dN_sds   =f_s(N_ϕ-N_srcosβ+κ_s T),   dM_sds   =f_s(M_ϕ-M_srcosβ-T).

By differentiating the definition of T and using the first of (<ref>), one finds that

    dTds=-f_s(κ_sN_s+κ_ϕ N_ϕ+Trcosβ).

Together with the geometrical equations r'=f_scosβ and β'=f_sκ_s, equations (<ref>) and (<ref>) describe the deformed shell. The five required boundary conditions can be read off the variation (<ref>) and the definition of T,


    β   = 0,    r   = 0,       T = 0      at the posterior pole,
    
    N_s    = 0,    M_s    = 0            at the phialopore.


We solve these equations numerically using the boundary value-problem solver  of Matlab (The MathWorks, Inc.). 

For completeness, we note that if external forces are applied to the shell, and δ𝒲 is the variation of the work done by these forces, then the variational condition is δℰ+δ𝒲=0. In that case, it is useful to write the variation (<ref>) in terms of δ r and δ z. We note that δ r'=-f_ssinβ δβ and δ z'=f_scosβ δβ, and so

    f_sδβ=cosβ δ z'-sinβ δ r'.

Using this geometric relation and integrating by parts, we obtain

    δℰ2π   = r M_s δβ+{rN_scosβ-sinβf_s(M_ϕcosβ-dds(rM_s))}δ r
       +{rN_ssinβ+cosβf_s(M_ϕcosβ-dds(rM_s))}δ z
       +∫{f_sN_ϕ-dds(rN_scosβ-sinβf_s(M_ϕcosβ-dds(rM_s)))} δ r ds
       -∫dds(rN_ssinβ+cosβf_s{M_ϕcosβ-dds(rM_s)})δ z ds.




  §.§.§ Limitations of the Theory

The theory presented here has a singularity in a biologically relevant limit: the intrinsic deformation gradient F^0 becomes singular at |κ_s^0|=(h/2)^-1 or |κ_ϕ^0|=(h/2)^-1. This value corresponds precisely to the case of cells that are constricted to a point at one cell pole. 

The way around this issue would presumably involve writing down an energy directly relative to the (possibly incompatible) intrinsic configuration of the shell. Working in the intrinsic configuration of the shell raises another issue to contend with, however: intrinsic volume conservation, which implies that the thickness H of the intrinsically deformed shell, which is close to the thickness of the deformed shell by assumption, differs from the thickness h of the undeformed shell. For a doubly curved shell, the relative thickness η=H/h is a function of both the intrinsic stretches f_s^0,f_ϕ^0 and the intrinsic curvatures κ_s^0,κ_ϕ^0. The volume of an element of shell is

    ∫_-H/2^H/2f_s^0f_ϕ^0(1-κ_s^0ζ)(1-κ_ϕ^0ζ)ρ ds dϕ dζ= f_s^0f_ϕ^0H(1+H^212κ_s^0κ_ϕ^0)ρ ds dϕ.

It follows that η satisfies satisfies the cubic equation

    (h^212f_s^0f_ϕ^0κ_s^0κ_ϕ^0)η^3+f_s^0f_ϕ^0η-(1+h^212R^2)=0,

the solution of which can be expressed in closed form. It is clear that this equation always has a solution if κ_s^0κ_ϕ^0>0. If κ_s^0κ_ϕ^0<0, there is a solution if and only if

    |κ_s^0κ_ϕ^0|<(4f_s^0f_ϕ^03h)^2(1+h^212R^2)^-2.

Since 16/9<4, this condition may fail before the intrinsic geometry becomes singular, so this additional condition is not vacuous. This brief discussion therefore points to some interesting, more fundamental problems in the theory of morphoelastic shells.

There is an additional subtlety associated with the geometric and intrinsic deformation gradient tensors in Eq. (<ref>): the components of F^g are expressed in (<ref>) relative to the (natural) mixed basis {ê⃗_⃗s⃗,ê⃗_⃗ϕ⃗}⊗{Ê⃗_⃗s⃗,Ê⃗_⃗ϕ⃗}, where ê⃗_⃗s⃗,ê⃗_⃗ϕ⃗ are the unit vectors tangent to the deformed configuration of the shell and Ê⃗_⃗s⃗,Ê⃗_⃗ϕ⃗ are defined analogously for the undeformed configuration. We have implicitly written down the components of F^0 relative to the same basis. In general however, the components of F^0 in (<ref>) are those relative to the basis {ê⃗^⃗0⃗_⃗s⃗,ê⃗^⃗0⃗_⃗ϕ⃗}⊗{Ê⃗_⃗s⃗,Ê⃗_⃗ϕ⃗}, where the unit basis {ê⃗^⃗0⃗_⃗s⃗,ê⃗^⃗0⃗_⃗ϕ⃗} can a priori be specified freely. We have neglected these additional degrees of freedom in the above derivation; the question of how to define a natural intrinsic tangent basis {ê⃗^⃗0⃗_⃗s⃗,ê⃗^⃗0⃗_⃗ϕ⃗} is however an interesting one, since the intrinsic stretches and curvatures need not be compatible.



 §.§ Fitting Embryo Shapes

For the purpose of fitting the model to the observed averages shapes, we define a family of piecewise constant or  linear functional forms for the intrinsic stretches and curvatures, shown in Fig. <ref>. This family of intrinsic stretches and curvatures is defined in terms of fifteen parameters, which are to be fitted for. Their functional forms are based on observations of cell shape changes by <cit.> summarised below: 


  * The intrinsic stretches f_s^0,f_ϕ^0 vary in the both hemispheres (Fig. <ref>a): in the posterior hemisphere, the initially teardrop-shaped cells thin into spindle-shaped cells (Fig. <ref>c,d, Fig. <ref>b), while, in the anterior hemisphere, they flatten into disc-shaped (`pancake-shaped') cells (Fig. <ref>d,e, Fig. <ref>c). While the evolution towards spindle-shaped cells appears to occur at the same time all over the posterior hemisphere, the data from thin sections suggest that the transition to disc-shaped cells starts at the bend region and progresses towards the phialopore (Fig. <ref>d,e). Moreover, the spindle-shaped cells are isotropic, f_s^0≈ f_ϕ^0, while the pancake-shaped cells are markedly anisotropic: next to the bend region, the long axis of their elliptical cross-section is the meridional one; next to the phialopore, it is the circumferential axis (Fig. <ref>c). 

  * The meridional intrinsic curvature κ_s^0 (Fig. <ref>b) is expected to vary most drastically in the region where paddle-shaped cells with thin wedge ends form (Fig. <ref>d, Fig. <ref>a). Because of the motion of cytoplasmic bridges relative to the cells, some additional, yet slighter, variation may be expected.

  * The variations of the circumferential intrinsic curvature κ_ϕ^0 are less clear: on the one hand, κ_ϕ^0 does not vary as drastically as the meridional one, because of the anisotropy of the paddle-shaped cells. This is a marked difference to type-A inversion, where the flasks-shaped cells are isotropic <cit.>, and both intrinsic curvatures therefore vary more dramatically in the bend region. On the other hand, some variation of the circumferential intrinsic curvature may be expected because of the motion of cytoplasmic bridges (Fig. <ref>c). We impose a continuous functional form for κ_ϕ^0, regularising a step function over a distance Δ s in arclength (Fig. <ref>c), but we do not fit for Δ s since we lack detailed information about the cell shape changes that define it.

The other geometrical parameter of the shell, the angular extent P of the phialopore, is not fitted for. We arbitrarily set P=0.3. The reasons for this simplification are discussed in the main text.

Numerical shapes were fitted to the average shapes by distributing M=100 points uniformly along the arclength of the numerical and average shapes, and minimising a Euclidean distance between them using Matlab (The MathWorks, Inc.) routine , modified as discussed above. A custom-written adaptive stepper was used to move about in parameter space and select the initial guess for the Nelder–Mead simplex. For each shape, the fit for the previous stage of inversion was used as the initial guess for the optimisation.





 §.§ Shape Perturbations and Statistical Statements

To define perturbations for the F=15 fitted model parameters P⃗_⃗0⃗∈ℝ^F at noise level δ, we draw independent N uniform random samples X⃗∼𝒰[0,1]^F on the unit interval and define the perturbed parameters P⃗=P⃗_⃗0⃗(1+2δ(X⃗-1)).


  §.§.§ Uniformity of the Distribution of Perturbations

As discussed in the main text, some of these perturbed parameters must be discarded. As a result, the samples that are retained are uniform on an unknown set 𝒜⊆[0,1]^F with means μ⃗. To establish that these means are not all the same, we derive confidence intervals for μ_i-μ_j. Since |X_i-X_j|⩽ 1, we may bound the variance of these differences by Var(X_i-X_j)⩽ 1, and hence, by the central limit theorem, a 100(1-p)% confidence interval is

    ⟨ X_i⟩-⟨ X_j⟩±z√(N),  z=Φ^-1(1-p/2F2),

wherein Φ^-1 is the inverse of the cumulative distribution function of the 𝒩(0,1) distribution, and where have included a multiple-testing correction. At noise level δ=0.075, we have run 10000 perturbations, finding M=max⟨X⃗⟩≈ 0.526 and m=min⟨X⃗⟩≈ 0.485. With M-m≈ 0.041 and , we infer that the 99% confidence interval for the maximum difference of the means does not contain zero, and hence that the means are not all the same. We notice however that these deviations of the means are small, in that they are not statistically signficantly different from 0.5. 


  §.§.§ Position of the Maxima of Shape Variation

We now make quantitative our statement, based on the cumulative distributions in Fig. <ref>c, that the experimental distribution of shape variation (with a maximum in the anterior cap) is very unlikely to arise under the uniform model. We ask: what is the probability p, under the uniform model, for the maximum in shape variation to lie in the anterior cap (Fig. <ref>c)? For 10000 perturbations, we found that 757 had a maximum in the anterior cap. Among these perturbations, 2345 yielded a single maximum in shape variation, with 60 of these maxima in the anterior cap. With 99% confidence, we therefore have upper bounds p<0.0757+0.0129<0.09 from all perturbations, and p<0.0256+0.0266<0.06 if we restrict to shape variations with a single maximum. 



§ ACKNOWLEDGEMENTS

We are grateful to D. Page-Croft and C. Hitch for instrument fabrication. We thank S. Hilgenfeldt for asking a question at the right time, T. B. Berrett for a conversation on matters statistical, and the Engineering and Physical Sciences Research Council, the Schlumberger Chair Fund, and the Wellcome Trust for partially funding this work.





In this appendix, we analyse the configuration where the rim of the phialopore is in contact with the inverted posterior for completeness of the mechanical analysis. We also analyse a toy problem to illustrate the intricate interplay of geometry and mechanics during contact.





 §.§ Elastic Model in the Contact Configuration

Let P be the angular extent of the axisymmetric phialopore at the anterior pole of the shell. Here, we discuss the contact problem where the shell has deformed in such a way that the rim of the phialopore (at θ=π-P=Q, where θ=s/R is the polar angle) is in contact with the shell at some as yet unknown position θ=C, as shown in Fig. <ref>a,b. 




    < g r a p h i c s >

figureAnalysis of the contact problem. (a) Undeformed configuration and (b) contact configuration. The phialopore at θ=Q=π-P touches the shell at θ=C, where θ is the polar angle. (c) Increasing circumferential stretch f_ϕ with advancing position θ_p of the peeling front, at constant intrinsic stretch f_ϕ^0. Insets: configuration with inverted posterior (as in Fig. <ref>d), at beginning of contact, and at a later stage. (d) Advancing contact position with advancing peeling front. Insets: configurations at beginning of contact and at later stage, as in (c).




As in the derivation of the governing equations without contact (Methods), we shall express the variations in terms of δ r and δβ. The third variation, δ z, is not independent of the former two, and so the condition that the vertical positions of the shell at the point of contact and at the phialopore match must be incorporated via a Lagrange multiplier, U. <cit.> raised a related issue in the derivation of the shape equations for vesicles. The Lagrangian for the problem is therefore

    ℒ=ℰ-2π U∫_C^Qf_ssinβ dθ,

where the prefactor has been introduced for mere convenience. We note the variation of Eq. (<ref>),

    δℒ2π=δℰ2π- Utanβ δ r_C_+^Q+Uf_s(C_+)sinβ(C) δ C+U∫_C^Q{f_sκ_s^2β δ r-f_sβ δβ} dθ.

Next, expanding the condition β(C_-)=β(C_+) of geometric continuity that we have already implicitly applied in the above, we note that

    δβ(C_-)+f_s(C_-)κ_s(C_-) δ C=δβ(C_+)+f_s(C_+)κ_s(C_+) δ C.

Since the outer part of the shell can rotate freely with respect to the inner part at the points of contact, the variations δβ(C_±) and δβ(Q) are, by contrast, independent. This is not true of the variations δ r(C_±) and δ r(Q), however:

    δ r(Q) = δ r(C_-)+f_s(C_-)cosβ(C) δ C=δ r(C_+)+f_s(C_+)cosβ(C) δ C.

Analogous expansions were used by <cit.> for discussing an adhesion problem for vesicles. Next, a straightforward calculation reveals that the governing equations (<ref>) and (<ref>) remain unchanged if we define T=-N_stanβ+Uβ/r for C⩽θ⩽ Q. For convenience, we adjoin the equation dz/ds=f_ssinβ to the system (thereby fixing the degree of freedom of vertical translation). The system thus becomes a system of six first-order differential equations on two regions, with two unknown parameters (the contact position C and the Lagrange multiplier U). We thus have to impose fourteen boundary conditions:


    r(0)   =0,   z(0)   =0,   β(0)   =0,   T(0)   =0,
    
    r(Q)   = r(C),   z(Q)   = z(C),   N_s(Q)   =0,   M_s(Q)   =0,

as well as the continuity conditions at θ=C,

    β   =0,    r   =0,    z   =0,    M_s   =0,

and

    r(C) N_sβ(C)-Utanβ(C)=r(Q)N_s(Q)β(Q)-Utanβ(Q),
        T = - N_stanβ(C) + Uβ(C)r(C),    ℰ2π=r(C) f_sN_s+r(C)M_s(C) f_sκ_s.


We also note that the conditions r(Q)=r(C) and z(C)=z(Q) do not take into account the finite, but small, thickness of the shell. A more detailed condition would require knowledge of the nature of the contact (and is anyway beyond the remit of a thin shell theory). 

We briefly explore shapes in the contact configuration in what follows. We start from a configuration where the posterior hemisphere has inverted, as in Fig. <ref>d, and advance the peeling front, but now without increasing the intrinsic circumferential stretch f_ϕ^0 at the phialopore. As the peeling front advances, the circumferential stretch at the phialopore increases (Fig. <ref>c) at constant f_ϕ^0, showing how the phialopore is pushed open by the posterior hemisphere. The procession of the point of contact between the posterior and the phialopore along the inverted posterior speeds up with advancing peeling front position (Fig. <ref>d) because the closer the point of contact is to the posterior, the more the latter resists the progression of the contact point because of the changing tangent angle.

The inset configurations in Fig. <ref>c,d also suggest that, as the peeling front advances, the regime of contact at a point discussed here gives way to a second contact regime, where the contact is over a finite extent of the meridian of the shell. We do not pursue this further.



 §.§ Asymptotic Analysis of a Toy Problem

Some analytic progress can be made and additional insight into the contact configuration can be gained by asymptotic analysis of a toy problem: two elastic spherical shells, an inner shell of radius R_1 and an outer, open shell of radius R_2>R_1, touch at the respective angular positions Θ_1 and Θ_2<Θ_1 (Fig. <ref>a), so that R_2/R_1=sinΘ_1/sinΘ_2. The intrinsic stretches and curvatures are those of the undeformed shells. For the remainder of this section, we non-dimensionalise distances with respect to the radius R_1 of the inner shell; stresses we non-dimensionalise with Eh. 

If the outer shell is moved relative to the inner shell by a distance d (Fig. <ref>b), the two shells deform in asymptotically small regions near the point of contact. This point of contact moves a distance dΞ down along the inner shell, determined by matching the displacements of the contact point and the forces exerted by one shell on the other. We assume in particular that the nature of the contact is such that the shells do not exert  torques on each other. Since we have non-dimensionalised distances with R_1, our asymptotic small parameter is

    ε^2=112(1-ν^2)h^2R_1^2≪ 1.




    < g r a p h i c s >

figureAsymptotic toy contact problem. (a) Two shells of radii R_1 and R_2 are in contact at angular positions Θ_1 and Θ_2, respectively. (b) Relative motion of one shell with respect to the other by a distance d induces deformations of the shell in an asymptotic inner layer of size δ, and causes the point of contact to move by a distance dΞ along the inner shell. (c) Contours of Ξ in the (Θ_1,Θ_2) plane.



The classical leading-order scalings for this problem are discussed by <cit.>, for example: deformations are localised to asymptotic inner regions of width δ∼ε^1/2, in which deviations of the tangent angle from its equilibrium value are of order d/δ, and we assume that d≪δ. We introduce an inner coordinate ξ, and write the polar angles as θ_1=Θ_1+δξ+𝒪(d), θ_2=Θ_2+δξ. We thus expand

    β_1(θ_1)   =Θ_1+(d/δ)b_1(ξ),    β_2(θ_2)   =Θ_2+(d/δ)b_2(ξ).

Assuming that δ^2≪ d≪δ, we then have the leading-order expansions


    N_s^(1)   =Eh δ d σ_1(ξ),    N_s^(2)   =Eh δ d σ_2(ξ),
    
    N_ϕ^(1)   (∗)=Eh E_ϕ^(1)+ν N_s^(1)=Eh δ a_1(ξ),   N_ϕ^(2)   (∗)=Eh E_ϕ^(2)+ν N_s^(2)=Eh δ a_2(ξ),


where a_1,a_2 are hoop strains. We note that the relations marked (∗) are only valid at leading order, where we may approximate f_s≈ f_ϕ≈ 1. Let F_r and F_z denote the (suitably scaled) radial and vertical forces exerted by the outer shell on the inner shell. We obtain the leading-order force balances from the energy variation (<ref>): using dashes to denote differentiation with respect to ξ,

    σ_1'sin^2Θ_1-b_1”'cosΘ_1sinΘ_1=F_zδ(ξ),      σ_1'sinΘ_1cosΘ_1-a_1+b_1”'sin^2Θ_1=F_rδ(ξ).

This system is closed, at leading order, by the geometric relation a_1'=-b_1, as in <cit.>. Eliminating σ_1, we obtain

    b_1””+b_1=(F_r-F_zΘ_1)δ'(ξ).

The matching conditions b_1→ 0 as ξ→±∞ reduce the number of undetermined constants to four, which are determined by the jump conditions at the contact point ξ=0.

The asymptotic balance for the outer shell is of course the same, but we must remember that the system has been non-dimensionalised with the radius of the inner shell, for which reason a geometric factor arises in the equations. Thus

    b_2””+(sinΘ_1sinΘ_2)^4b_2=0,

with the matching condition b_2→ 0 as ξ→∞, leaving two boundary conditions to be imposed on this equation. Since the shells do not exert any moments on each other, b_2'(0)=0. The second condition is obtained from the force balance: the vertical force balance can be integrated once to yield

    sinΘ_1sinΘ_2{σ_2-Θ_2(sinΘ_2sinΘ_1)^4b_2”}=F_z.

Matching to the undeformed, unstressed shell as ξ→∞ implies F_z=0. The radial force boundary condition resulting from (<ref>) is

    sinΘ_1cosΘ_2{σ_2(0)+tanΘ_2(sinΘ_2sinΘ_1)^4b_2”(0)}=F_r,

which, upon imposing (<ref>), reduces to

    b_2”(0)=(sinΘ_1sinΘ_2)^3F_r.

Let U_r^(1),U_z^(1) and U_r^(2),U_z^(2) denote the respective (non-dimensional) displacements of the contact point ξ=0, scaled with d. Then


    U_r^(1)   =sinΘ_1∫_0^∞b_1 dξ=-F_r2√(2)sinΘ_1,   U_z^(1)   =-cosΘ_1∫_0^∞b_1 dξ=F_r2√(2)cosΘ_1,
    
    U_r^(2)   =sinΘ_1∫_0^∞b_2 dξ=-√(2)F_rsinΘ_1,   U_z^(2)   =-sinΘ_1Θ_2∫_0^∞b_2 dξ=√(2)F_rsinΘ_1Θ_2.


In particular, these expressions once again contain additional geometric factors resulting from the non-dimensionalisation. 

The values of the two remaining undetermined constants, F_r and Ξ, are finally obtained by imposing continuity of the displacement of the contact point, i.e.

    U_r^(1)+ΞcosΘ_1   =U_r^(2),   U_z^(1)+ΞsinΘ_1   =U_z^(2)+1.

Notice that arclength is computed here from the anterior pole of the shell to match the asymptotic setup of <cit.>, and so the `vertical' axis is pointing downwards in Fig. <ref>, giving rise to some sign changes. In particular, we obtain

    Ξ=3sinΘ_11+2cosecΘ_2sin(2Θ_1-Θ_2).

The contours of this expression are plotted in Fig. <ref>c. The very non-linear nature of this expression illustrates that the contact geometry is quite intricate; in particular, Θ_2(Θ_1) at fixed Ξ is not a monotonic function, but, as expected (since it is easier for the the contact point to slide along the inner shell the more parallel it is to the axis of symmetry), at fixed Θ_1, Ξ increases with Θ_2.

