Understanding how forces and material properties give rise to tissue shapes is a fundamental issue in developmental biology. Although Drosophila gastrulation is a well-used system for investigating tissue morphogenesis, a consensus mechanical model that explains all the key features of this process does not exist. One key feature of Drosophila gastrulation is its anisotropy: the mesoderm constricts much more along one axis than along the other. Previous explanations have involved graded stress, anisotropic stresses or material properties, or mechanosensitive feedback. Here, we show that these mechanisms are not required to explain the anisotropy of constriction. Instead, constriction can be anisotropic if only two conditions are met: the tissue is elastic, as was demonstrated in our recent study; and the contractile domain is asymmetric. This conclusion is general and does not depend on the values of model parameters. Our model can explain results from classical tissue-grafting experiments and from more-recent laser ablation studies. Furthermore, our model may provide alternative explanations for experiments in other developmental systems, including C. elegans and zebrafish.

A fundamental process in animal development is gastrulation, whereby a single epithelial sheet gives rise to a multilayered structure (Sweeton et al., 1991; Leptin and Grunewald, 1990). In Drosophila, gastrulation initiates when the prospective mesoderm cells, which are located in a rectangular domain in the ventral part of the embryo, constrict apically (i.e. on the outer surface of the embryo, see schematic in Fig. 1). As the mesodermal cells constrict across their apical faces, they become wedge shaped and elongate along the apico-basal axis (Sweeton et al., 1991; Leptin and Grunewald, 1990). Apical constriction of the mesodermal cells is widely believed to be driven by myosin-generated active stresses in the apical domains of those cells (Dawes-Hoang et al., 2005). In accordance with this, the concentration of apically localized myosin increases dramatically in the constricting cells concomitantly with the onset of apical constriction (Martin et al., 2010). Following apical constriction, the surface of the embryo forms a furrow at the ventral midline. The furrow deepens and closes off; in this way, the mesoderm is brought into the interior of the embryo.

The focus of this study is the initial phase of gastrulation when mesodermal cells constrict apically, before the surface begins to fold inward. Notably, as mesodermal cells shrink, the rectangular mesodermal domain contracts strongly along its shorter axis and much less so along the longer axis. For brevity, the length of the mesodermal domain along the shorter mediolateral axis will be referred to as the ‘width’, and the length of its longer anteroposteior axis will be referred to as the ‘length’. Thus, the length-to-width ratio increases drastically in the course of the constriction of the mesoderm, see schematic in Fig. 1. Because mesoderm constriction in Drosophila melanogaster is not accompanied by appreciable cell rearrangements, the length-to-width ratio of individual mesodermal cells must increase in the course of tissue contraction (Martin et al., 2010; Spahn and Reuter, 2013).

This anisotropic constriction has been the subject of considerable interest. One key experiment was carried out in a classical paper by Maria Leptin and Siegfried Roth (Leptin and Roth, 1994), in which patches of mesodermal tissue (either nuclei or cytoplasm) were transplanted into mutant embryos that lacked mesoderm. If the contractile domain were symmetric (e.g. round or square) it would constrict symmetrically, whereas if the contractile domain were elongated, it would constrict anisotropically. Furthermore, the geometry of contraction of mesodermal patches was independent of the position and orientation of those patches within the embryo and was dependent only on the geometry of the patch. This showed that anisotropic constriction is a locally autonomous effect, and is not, for example, induced by maternal gradients (i.e. factors that localize in a patterned manner in a freshly laid egg, notably Bicoid, Nanos, Torsolike and Dorsal). The mechanism driving anistropic constriction remained unclear, however. Two more recent papers have proposed graded tension (Spahn and Reuter, 2013) or mechanosensitivity (in particular, feedback between stress and nematic order of actin filaments; Chanet et al., 2017) as the source of anisotropy. Here, we propose a simple mechanistic explanation for anisotropic constriction of the mesoderm that does not invoke any of these effects.

The model

We start by describing model assumptions, then present the corresponding equations below. We model the surface of the embryo as a flat (two-dimensional) elastic sheet. This elastic sheet may be pictured as a network of spheres connected by springs as shown in Fig. 1. The model sheet will contain a rectangular contractile domain with horizontal and vertical dimensions Lx (length) and Ly (width), respectively (see Fig. 1). The size of the contractile domain is taken to be much smaller than the size of the simulated sheet in order to avoid boundary effects. In order for this contractile domain to shrink it must carry active contractile stress, or ‘tension’, which is measured in units of force per unit length. In Fig. 1, the presence of active stress in the contractile domain is illustrated by red arrows, which are directed in parallel with the springs. It is assumed that the active stress within the contractile domain is uniform (having the same magnitude everywhere) and isotropic (same magnitude in every direction). In particular, in Fig. 1, all motors are the same regardless of the orientation of the corresponding spring or the position of the spring. Outside of the contractile domain, the actively generated stress is absent. In the real embryo, active contractile stress is generated by active myosin motors. As active contractile stresses do not depend on instantaneous displacement, there is no feedback to force generation in this model.

The above assumptions make up a complete description of the model. Combining these assumptions with standard equations of linear elasticity, assuming that the elastic sheet is in mechanical equilibrium [see Landau and Lifshitz (1970), in particular equation 13.4 therein], we have:
(1)
where u is the in-plane deformation of the sheet, E is Young's modulus (more accurately, this quantity represents Eh, where E is the three-dimensional Young's modulus and h is the thickness of the sheet), σ is the Poisson's ratio and μ is the active stress tensor within the contractile mesodermal domain (and does not, in general, correspond to a conservative force, e.g. Jülicher et al., 2007). We take active stress to be isotropic, μ=φI, with I being the identity matrix, and φ(r)=φ0Π(x/Lx)Π(y/Ly), where Π denotes the ‘top hat’ function that equals identity between −0.5 and 0.5 and zero otherwise. The magnitude of active stress is constant in the mesodermal domain and zero outside. (If some uniform active stress is present outside the contractile domain, φ0 may be considered the difference between the active stress inside and that outside. Since Eqn 1 takes a derivative of stress, the solution depends only on the difference in contractility between the domains.) We ignore viscous relaxation of stresses, as it was found experimentally that those stresses persist on a time-scale comparable with (but not significantly shorter than) the time-scale of tissue dynamics (Doubrovinski et al., 2017). Finally, we note that the assumption that the elastic sheet is in mechanical equilibrium is an idealization: it is not known whether developmental dynamics proceed adiabatically.

The model is schematically illustrated in Fig. 1. The elastic sheet is represented by a spring network with nodes illustrated as spherical beads. For simplicity, the illustration was made somewhat misleading: a network of squares is not isotropic, whereas Eqn 1 describes an isotropic elastic material. Direct rheological measurements have demonstrated that the membrane surface of embryonic tissue is highly elastic and that elasticity persists on a time-scale that is well comparable with the time-scale of gastrulation (Doubrovinski et al., 2017). Myosin-generated tension in the mesodermal domain is well established (Dawes-Hoang et al., 2005; Martin et al., 2009, 2010).

Additionally, the proposed model involves only a few dimensional parameters. The elastic sheet that represents the epithelium has Young's modulus E, which determines the force required to stretch material a given amount and may be thought of as a multi-dimensional counterpart of a spring constant. An elastic material is further characterized by its Poisson's ratio, σ, which determines how much a material shrinks in one direction when it is being stretched along a perpendicular direction. Finally, the active stress (or, equivalently, tension) inside the contractile domain has magnitude of φ0. As tension is assumed not to vary with either space or direction within the contractile domain, only one number is needed to describe this active stress.

Analysis of the model: constriction geometry

To analyze the behavior of the minimal model (Eqn 1) we performed numerical simulations; implementation details are given in the Materials and Methods. If the contractile domain is chosen to be asymmetric (i.e. Lx>> Ly, or length is larger than width), simulations reveal constriction to be greater along the shorter axis of the domain, see Fig. 2A-C and Fig. S1. In contrast, when the contractile domain is chosen to be symmetric (i.e. Lx=Ly), constriction is equal along both axes, Fig. 2F. Thus, this simple model, which involves only elastic material and an isotropic contractile force, replicates the key experimental finding that the anisotropy versus isotropy of constriction is dependent on geometry (Leptin and Roth, 1994).

To examine the influence of model parameters on these predictions, we performed a parameter sweep (Fig. 2D,E). The non-dimensionalized version of Eqn 1 has the following dimensionless parameters: the ratio of Young's modulus to active stress magnitude E/φ0, Poisson's ratio σ and the length-to-width ratio of the contractile domain Lx/Ly. The value of Poisson's ratio σ was varied over the whole physically relevant range, which is between 0 and 1/2 (Greaves et al., 2011). It is seen that the qualitative nature of the solution does not vary appreciably with the value of the Poisson's ratio. Simulation results for three different choices of the active stress-to-stiffness ratio φ0/E are given in Fig. 2D. Since Eqn 1 is linear, the displacement is proportional to active stress and therefore the solutions (i.e. the ux and uy components of the displacement, shown for one particular choice of parameters in Fig. 2C) for different φ0/E are related by simple scaling. From this, it easily follows that anisotropy is zero in the limit of zero force and increases arbitrarily for sufficiently high values of φ0/E. Additionally, we systematically examined the influence of the initial aspect ratio of the contractile domain on the final equilibrium geometry (Fig. 2E); to achieve this, we used analytical expressions in the Materials and Methods. It was found that contraction is always asymmetric, although the degree of the asymmetry may vary. Symmetric domains contract symmetrically (Fig. 2E). In summary, our results are not sensitive to the choice of the elasticity parameters.

To examine the influence of the domain in which the contractile domain is embedded, we performed additional simulations. To achieve this, we simulated square as well as rectangular contractile domains embedded in either square or rectangular stress-free domains. The deformation of the contractile domain becomes essentially independent of the geometry of the embedding domain, provided the embedding domain was chosen to be sufficiently large (Fig. S2C). To further formally examine the influence of the geometry of the embedding domain, we systematically determined the final equilibrium aspect ratio of the contractile domain as a function of model parameters and the size of the embedding domain. Fig. S3 shows that for a fixed geometry of the contractile domain and fixed model parameters, the equilibrium ratio converges to a limit as the embedding domain is made larger and larger. Finally, to further study the influence of the boundary conditions, we compared simulations with both free and fixed boundaries using finite element simulations (Fig. S2B). Our results also indicate close quantitative agreement in this case, again confirming that the boundary effects are irrelevant as long as the embedding domain is chosen to be sufficiently large. As convergence with respect to the size of the embedding domain is reached for very modest sizes of the embedding domain (for the relevant aspect ratio of the contractile domain, see Fig. S3), we believe that our results are fully relevant to the biological situation.

Importantly, if the contractile domain is not embedded in an elastic continuum, the equilibrium configuration will not be anisotropic, as may be seen by directly solving Eqn 1 analytically. In this simpler case, a solution can be found by assuming that the solution is linear, i.e. ux=bxx, uy=byy (with bx,y being constants), and using the boundary conditions. Thus, the conclusion that the deformation remains anisotropic at long times appears to be a consequence of the particular geometry in our problem. The main effect described in our model is a special case of a more-general phenomenon when an (elastic) solid body is subjected to uniform isotropic active stress: when mechanical equilibrium sets in, total stress needs not remain either isotropic or uniform. A further discussion of a similar effect in the context of a different developmental model, C. elegans, may be found in Vuong-Brender et al. (2017). In addition, related effects for the case of a viscous model have been previously discussed by Salbreux et al. (2009).

Analysis of the model: laser ablations

Next, we tested whether the model correctly predicts the outcome of laser ablation experiments. Specifically, it has been shown that when the surface of the mesoderm is laser ablated to create an initially round hole, the hole will expand predominantly along the longer axis of the embryo. This anisotropy of the response to ablation appears to increase as gastrulation proceeds: immediately after the onset of gastrulation the response to ablation is approximately isotropic; later ablations become gradually more anisotropic with time (Martin et al., 2010).

To simulate these experiments, we implement the model with a rectangular contractile domain as before, but after the deformation field has evolved to mechanical equilibrium edges within a circular region are removed, to mimic a laser ablation (Fig. 3A). The initially round hole expands anisotropically, with its long axis aligned with the long axis of the contractile domain. This result qualitatively agrees well with the experiments.

To see this intuitively, let us again consider the model in Fig. 2B,C. After the onset of contraction, the domain will shrink more along the shorter vertical direction and less along the perpendicular horizontal direction. Now, consider what would happen if one were to remove a single vertical spring-edge in the interior of the contractile domain after contraction sets in (Fig. 3B). This vertical spring-edge has two adjacent vertical edges, colored blue in Fig. 3B. If the contractile domain had already undergone some contraction, those adjacent vertical neighbors would have shrunk to become shorter than their preferred rest length. Spring forces in those adjacent spring edges would tend to expand those adjacent springs, thus counteracting the expansion of the ‘hole’ along the vertical direction. Now, instead, consider what would happen if a horizontal spring-edge is removed. The adjacent horizontal spring-edges (colored green in Fig. 3B) are under relatively less compression (i.e. deviate less from their rest length). Thus, elastic forces in those springs will counteract the expansion of the ‘hole’ to a lesser extent than was the case with the vertical spring-edge. In this argument, it has been assumed that active forces on every spring in the interior of the contractile domain are the same on every edge of the network. In order to correctly interpret the outcome of the ablation experiments, it is absolutely necessary to distinguish between ‘total stress’, ‘elastic stress’ and ‘active stress’: elastic stress is due to the springs; active stress is due to the motors (myosin); and total stress is the sum of the two. Ablation experiments measure total stress, which becomes increasingly anisotropic over the course of the contraction. Active stress, however, remains isotropic throughout the course of the dynamics. Previously, the anisotropy in response to ablation was interpreted as evidence for anisotropy of the active stress (Martin et al., 2010). The present model provides a significantly simpler interpretation: total stress must become anisotropic if the domain is asymmetric (and if the tissue is elastic), whereas active stress need not be anisotropic. In addition, in the present model, the asymmetry of contraction is the cause, not a consequence, of anisotropic total stress as measured by laser ablation.

This paper introduces a simple model that can explain the asymmetry of contraction of the mesoderm in Drosophila melanogaster. In the proposed model, the asymmetry of contraction is a consequence of the asymmetry of the geometry of the contractile domain and not because of, for example, anisotropies in active stresses or graded distribution of those stresses within the contractile domain. Importantly, in our model, anisotropy arises from interactions of the mesoderm with the ambient tissue and is not an inherent property of the mesoderm itself. Let us now compare and contrast the proposed mechanism with a number of alternative models.

One alternative possibility is that active stresses in the contractile domain are anisotropic, i.e. the cells have a preferred axis along which to exert active stress. We believe that this scenario is unlikely. First, there is no evidence of planar cell polarity in prospective mesodermal cells of the ventral furrow. Moreover, this interpretation is at odds with the observations by Leptin and Roth (Leptin and Roth, 1994), where it was shown that contractile grafts of mesodermal tissue contract along their shorter dimension regardless of their position or orientation within the embryo. If active stresses in those contractile patches were anisotropic, a mechanism has to exist to align those active stresses with the shorter axis of the domain through some self-organized process. This general idea has been proposed, but a concrete mechanism has not been demonstrated (Chanet et al., 2017). Although not impossible, this appears a much more complicated explanation than the one we propose. Additionally, if active stresses were anisotropic at the onset of contraction, the anisotropy as measured by laser ablation experiments would be present immediately after the onset of contraction, rather than building up gradually, as has been shown experimentally (Martin et al., 2010).

A second model could propose that material properties of the contractile domain are anisotropic. For example, one can imagine that the cells of the contractile domain are easier to contract along the mediolateral axis than along the anteroposterior direction. This possibility may be refuted on the same basis as the possibility of anisotropic active stresses.

Other possible explanations include spatial nonuniformity of active stresses within the contractile domain, or spatial nonuniformity of material properties. For example, much like in the model proposed by Mayer et al. (2010) and Behrndt et al. (2012), one could imagine that active stresses are isotropic but form a gradient that peaks at the center of the contractile domain and decays gradually along the mediolateral direction. In this case, it can be shown that contraction would be directed along the mediolateral axis of the embryo (Spahn and Reuter, 2013). In fact, there are gradients of gene expression for two regulators of myosin, t48 and fog (Lim et al., 2017), and there may be a gradient of myosin density along the mediolateral axis of the embryo (Lim et al., 2017; Spahn and Reuter, 2013), although it is unclear from published data whether the myosin gradient appears early enough to be the cause rather than the consequence of anisotropic constriction. Although these data suggest that graded tension might contribute to anisotropic constriction, this potential mechanism is not sufficient to explain all experimental data. In particular, Leptin and Roth demonstrated that a patch of transplanted mesodermal cells will consistently constrict along its shorter axis (Leptin and Roth, 1994). In these experiments, the patch of transplanted mesodermal cells descended from a small number of transplanted nuclei and it seems highly unlikely that any gradient would be preserved in this process. Therefore, some alternative mechanism must drive anisotropic constriction in this case.

Anisotropy of mesoderm contraction has been the focus of a recent experimental paper (Chanet et al., 2017). Chanet et al. modified the geometry of the contractile domain by genetic means. The observations obtained (Chanet et al., 2017) are fully consistent with and complementary to the results from Leptin and Roth (1994): more symmetric domains constrict less anisotropically. According to Chanet et al. (2017), the anisotropy of constriction was attributed to a mechanosensitive response of mesodermal cells leading to re-orientation of actin filaments perpendicularly to the axis of contraction. Although our work does not rule out this possibility, it shows that such effects are not necessary to explain the key aspects of tissue dynamics.

Crucially, the phenomenology presented here may apply widely to other developmental systems. In a number of recent publications (Mayer et al., 2010; Behrndt et al., 2012), the response of embryonic tissue to laser ablation was examined in zebrafish or C. elegans. In these studies, the initially round hole generated through ablation expanded asymmetrically to become elliptical; as mentioned, the same results have been reported in gastrulating Drosophila tissue (Martin et al., 2010). In particular, the hole expands more in the direction perpendicular to the axis of tissue contraction. In the former papers, the observation was interpreted as an effect of anisotropic viscous shear (Mayer et al., 2010; Behrndt et al., 2012). In Chanet et al.’s study, this result was interpreted to mean that actively generated tension was anisotropic (Chanet et al., 2017). In the present study, it is shown that neither anisotropic viscous shear nor anisotropic active tension is required for an anisotropic response to laser ablation. Furthermore, it is shown that an anisotropic response to ablation is an inevitable consequence of tissue elasticity. In this way, the present study provides a simple parsimonious explanation of a series of important observations that we believe may have previously been interpreted incorrectly.

An earlier theoretical work by Spahn and Reuter (Spahn and Reuter, 2013) used a vertex model to study the mechanisms that underlie anisotropic tissue constriction during Drosophila gastrulation. Based on this modeling, the authors proposed two potential mechanisms to generate anisotropic constriction. One of these models relied on a gradient of contractility along the dorso-ventral axis; as discussed in the previous paragraph, work from Leptin and Roth (Leptin and Roth, 1994) indicates that such a gradient is dispensable for anisotropic constriction. In the second model proposed by Spahn and Reuter anisotropic constriction could be due to the contractile domain being rigidly anchored at the anterior and posterior domain boundaries (i.e. anisotropy being essentially a boundary effect). This solution seems biologically irrelevant: Drosophila mesoderm does not appear to be flanked by a row of cells that are being artificially held fixed at one end. Furthermore, our model explains anisotropic constriction without invoking such effects.

The model presented here qualitatively agrees with the available experimental data. In the future, this model could be tested more rigorously by systematically mapping the distribution of material properties of the tissue [similar to the approach used by Doubrovinski et al. (2017)] and the distribution of total stress (using, for example, laser ablation). With these data, the pattern of deformation can be determined uniquely without fitting.

Analytical treatment

In this section we solve Eqn 1 using the Fourier transform. We will be considering Eqn 1 on a rectangular domain with side lengths Λx and Λy. Rectangular contractile region x0<x<Lxx0, y0<y<Lyy0 is assumed to be under constant isotropic active stress of magnitude φ. We Fourier-expand u=[ux(x, y), uy(x, y)] as:
(A.1)
It is useful to introduce simplified notation αE/(2(1+σ)), βE/(2(1−σ)). Substituting (Eqn A.1) into Eqn 1 and collecting coefficients corresponding to the same modes we obtain:
(A.2)
where the constant coefficients , are given by:
(A.3)
With the obvious identifications, Eqn A.2 have the form:
(A.4)
which is readily solved as:
(A.5)
Combining Eqn A.1-A.5, one obtains an explicit expression for displacement u as a function of position. As the complete expression is rather lengthy and not especially telling, we do not provide it here. Instead, we present plots showing the solution u obtained by directly evaluating the (truncated) Fourier expansion in Fig. S1A.

Numerical simulations I

The code for all simulations is included in the supplementary information. Simulations in Fig. 2 were implemented as follows. We solve boundary problem 1 by introducing an additional term ∂tu on the right-hand side and evolving the resulting equation to equilibrium. All spatial derivatives were approximated as second-order central differences. The displacement of the (outermost) boundary was initially set to zero and not updated (fixed boundary conditions). As the deformation remained localized away from the boundary throughout the course of the simulation, the particular choice of the boundary conditions is expected to not influence the result. Time integration was carried out using the Euler-forward (explicit) scheme. The domain was discretized into a regular grid of 800×800; the Euler forward step was dt=10−5. In order to avoid development of singularities, we ‘smeared out’ the forces at the boundaries of the contractile domain over a region of finite size. Specifically, force distribution had Gaussian profile exp (−r2/ξ2), with ξ=0.3 (this was, however, found not to be necessary as no singularities developed in the limit of much smaller ξ). Parameters are listed in the corresponding figure caption.

Fig. 2D,E was generated by evaluating the Fourier-expanded expressions given in the next section. The 600 lowest Fourier modes were summed. We checked explicitly that analytical results closely matched the results of the numerical simulations. Although our treatment focuses on the description of a flat elastic sheet, it may readily be modified to describe a semi-infinite three-dimensional elastic continuum by replacing the coefficient in front of the last term of Eqn 1 with E/(1+σ)(1−2σ). The major conclusions about the anisotropy of the asymptotic state hold in this case as well.

Simulations in Fig. 3 were carried out using a numerical scheme that differed from simulations in Fig. 2 in order to simplify the implementation of ablation. Specifically, instead of using a finite difference scheme, we discretized the domain as an (unstructured) grid of equilateral triangles whose edges are linear springs. It has been shown that this approximation reduces to the equations of linear elasticity (i.e. Eqn 1) in the limit of small strains, see, for example, Seung and Nelson (1988) and Liang and Mahadevan (2011). Simulation parameters were chosen as follows. The simulated domain had a size of 2×2, the length of an individual spring edge was 0.03. Stiffness of an individual spring edge was set to 50. The contractile domain had a length of 1.2 and a width of 0.6, and was positioned in the center of the simulated domain. Edges within the contractile domain were subjected to compressive force of 2.5. We used fixed boundary conditions (as in Fig. 2). Ablation was simulated after mechanical equilibrium was reached (t=5) by effectively removing all nodes and edges within a circle of radius 0.025 in the center of the domain. Integration was carried out using Euler-forward scheme with a time step of 5 · 10−4.

For the parameter sweep of (1), we non-dimensionalized the equation. Dividing throughout by φ0, and rescaling length by the width of the contractile domain Ly one obtains the re-scaled equation in the same form as Eqn 1, except with φ0 set to 1, and Young's modulus E replaced by the dimensionless ratio E/φ0. The non-dimensionalized version of Eqn 1 has the following three dimensionless parameters: the ratio of Young's modulus to active stress magnitude E/φ0, Poisson's ratio σ and the length-to-width ratio of the contractile domain Lx/Ly. The plots in Figs S2C and S3 were constructed using analytical expressions for the solution in terms of Fourier-expansion, see Eqn A.1.

Numerical simulations II: finite element implementation

The displacement profiles shown in Fig. S2A,B were obtained by solving Eqn 1 with the finite element method (FEM) (Zienkiewicz and Taylor, 1991). To that end, we first rewrite Eqn 1 in its compact form as:
(A.6)
where φ(x, y) is the distribution of the active stress introduced in the main text, I is the identity tensor and s is the Cauchy stress tensor in the entire domain. This tensor is related to the displacement field u through the constitutive relation:
(A.7)
where is the strain tensor and λsμs are elastic constants that are, in the plane stress approximation, given by:
Here, E and σ are, respectively, the Young's modulus and Poisson's ratio introduced in the main text.

Weak formulation

In order to obtain the displacement field using FEM, we use a test function v and rewrite Eqn A.6 in the following weak form:
(A.8)

In Eqn A.8, Ω and ∂Ω are respectively the bulk and boundaries of the numerical domain, while n is the outward normal vector to ∂Ω.

At this stage, it is straightforward to account for the boundary conditions on the boundaries on the embedding domain. Both in the case of traction-free conditions and zero displacement conditions, the first term of Eqn A.8 vanishes. In the former case, this is easily obtained, as s · n=0 on ∂Ω. In the latter case, that term can be reduced to zero by choosing test functions v that vanish on ∂Ω. Therefore, Eqn A.8 can be simplified as:
(A.9)
Next, we discretize Eqn A.9 by expanding the displacement field in the basis of shape functions (chosen to be the same as the test functions): , where the uj is a scalar unknown to be found and N is the number of nodes of the computational domain. Inserting this decomposition into the weak form and letting v=vi, the displacement field is obtained by inverting the following linear system:
(A.10)
where are components of the stiffness matrix and are components of the forcing arising from the active stress.

We discretize the domain Ω with a triangular mesh, choose P2 elements as shape functions and solve the system of Eqn A.10 using the finite element open software FreeFem++ (Hecht, 2012). In building the mesh, we ensure that the contractile (or inner) domain is much more refined than the embedding (or outer) domain. In order to ensure a smooth transition of the mesh between the contractile and embedding parts, we include an intermediate zone where the mesh progresses from being very refined near the inner domain to being less refined as one approaches the outer domain. Fig. S4 shows the mesh used to obtain the displacement field that we plot in Fig. S2A′ and Fig. S2B′.

Boundary conditions treatment

Case of zero displacement condition

In the case of fixed boundary conditions, we use a penalty technique to enforce uj=0 on the nodes corresponding to the domain boundary. That is, if p is the index of a node located on the boundary, we ascribe to it a very large number and write App=1030 as well as Fp=0×1030. Consequently, on the row of the stiffness matrix corresponding to that node, we have , thus leading to up=0.

Case of zero traction condition

In the case of traction-free boundary conditions, we simply solve Eqn A.10, which was obtained from Eqn A.8 with the traction term (s · n) put to zero on ∂Ω, the boundary of the embedding domain.

However, unlike in the case of zero displacement conditions on ∂Ω, the embedding domain here is kinematically unconstrained. As a result, we find numerical solutions of Eqn A.10 that contain contributions from rigid modes of the system, i.e. the eigenmodes of the stiffness matrix with zero eigenvalue. Indeed, if us is a solution of Equation A.10, and if ur is an eigenmode of the matrix Aij with zero eigenvalue, then clearly uh=us+ur is also a solution of Eqn A.10.

In order to compare the displacement field obtained under the conditions of zero displacement with that obtained under the condition of zero traction on ∂Ω, we must first remove the rigid modes contributions from the numerical solution (uh). This ensures that both displacement fields correspond to the elastic deformations due the active force only and, in absence of the latter force, the displacement fields would uniformly zero. To that end, we proceed as follows.

We first rewrite the numerical solution in the basis of the eigenmodes mj of the stiffness matrix, i.e. , where the coefficients of the expansion are given by the dot products . Then, by identifying the rigid modes as those of the stiffness matrix with zero eigenvalue, we find three such modes (two translations and one rotation) that we note mk with k=1, 2, 3. Last, by removing these rigid displacements from the previously computed solution, we obtain the following displacement field due to the active force only:
(A.11)
It is the displacement field us that we show in Fig. S2B and Fig. S2B.

We thank Miriam Osterfield for helpful discussions.

Author contributions

Conceptualization: K.D., J.T., K.M.; Methodology: K.D., J.T., K.M.; Software: K.D., J.T., K.M.; Validation: K.D., J.T., K.M.; Formal analysis: K.D., J.T., K.M.; Investigation: K.D., J.T., K.M.; Writing - original draft: K.D.; Writing - review & editing: K.D., J.T., K.M.; Supervision: K.D., K.M.

Funding

We acknowledge support from the Robert A. Welch Foundation (I-1950-20180324) and the National Institutes of Health (R01-GM110066). Deposited in PMC for release after 12 months.

Behrndt
,
M.
,
Salbreux
,
G.
,
Campinho
,
P.
,
Hauschild
,
R.
,
Oswald
,
F.
,
Roensch
,
J.
,
Grill
,
S. W.
and
Heisenberg
,
C.-P.
(
2012
).
Forces driving epithelial spreading in zebrafish gastrulation
.
Science
338
,
257
.
Chanet
,
S.
,
Miller
,
J. C.
,
Vaishnav
,
E. D.
,
Ermentrout
,
B.
,
Davidson
,
L. A.
and
Adam Martin
,
A. C.
(
2017
).
Actomyosin meshwork mechanosensing enables tissue shape to orient cell force
.
Nature
8
,
15014
.
Dawes-Hoang
,
R. E.
,
Parmar
,
K. M.
,
Christiansen
,
E. A.
,
Phelps
,
C. B.
,
Brand
,
A. H.
and
Wieschaus
,
E. F.
(
2005
).
folded gastrulation, cell shape change and the control of myosin localization
.
Development
132
,
4165
-
4178
.
Doubrovinski
,
K.
,
Swan
,
M.
,
Polyakov
,
O.
and
Wieschaus
,
E. F.
(
2017
).
Measurement of cortical elasticity in Drosophila melanogaster embryos using ferrofluids
.
Proc. Natl. Acad. Sci. USA
114
,
1051
-
1056
.
Greaves
,
G. N.
,
Greer
,
A. L.
,
Lakes
,
R. S.
and
Rouxel
,
T.
(
2011
).
Poisson's ratio and modern materials
.
Nat. Mater.
10
,
823
-
837
.
Hecht
,
F.
(
2012
).
New development in freefem++
.
J. Num. Math.
20
,
251
.
Jülicher
,
F.
,
Kruse
,
K.
,
Prost
,
J.
and
Joanny
,
J.-F.
(
2007
).
Active behavior of the cytoskeleton
.
Phys. Rep.
449
,
3
-
28
.
Landau
,
L. D.
and
Lifshitz
,
E. M.
(
1970
).
Theory of Elasticity (Volume 7 of A Course of Theoretical Physics)
.
Pergamon Press
.
Leptin
,
M.
and
Grunewald
,
B.
(
1990
).
Cell shape changes during gastrulation in Drosophila
.
Development
110
,
73
-
84
.
Leptin
,
M.
and
Roth
,
S.
(
1994
).
Autonomy and non-autonomy in Drosophila mesoderm determination and morphogenesis
.
Develoment
120
,
853
-
859
.
Liang
,
H.
and
Mahadevan
,
L.
(
2011
).
Growth, geometry, and mechanics of a blooming lily
.
Proc. Natl. Acad. Sci. USA
108
,
5516
-
5521
.
Lim
,
B.
,
Levine
,
M.
and
Yamazaki
,
Y.
(
2017
).
Transcriptional pre-patterning of Drosophila gastrulation
.
Curr. Biol.
27
,
286
-
290
.
Martin
,
A.
,
Kaschube
,
M.
and
Wieschaus
,
E.
(
2009
).
Pulsed contractions of an actin-myosin network drive apical constriction
.
Nature
457
,
495
-
499
.
Martin
,
A. C.
,
Gelbart
,
M.
,
Fernandez-Gonzalez
,
R.
,
Kaschube
,
M.
and
Wieschaus
,
E.
(
2010
).
Integration of contractile forces during tissue invagination
.
J. Cell Biol.
188
,
735
.
Mayer
,
M.
,
Depken
,
M.
,
Bois
,
J. S.
,
Jülicher
,
F.
and
Grill
,
S. W.
(
2010
).
Anisotropies in cortical tension reveal the physical basis of polarizing cortical flows
.
Science
467
,
617
-
621
.
Salbreux
,
G.
,
Prost
,
J.
and
Joanny
,
J. F.
(
2009
).
Hydrodynamics of cellular cortical flows and the formation of contractile rings
.
Phys. Rev. Lett.
103
,
058102
.
Seung
,
H. S.
and
Nelson
,
D. R.
(
1988
).
Defects in flexible membranes with crystalline order
.
Phys. Rev. A.
38
,
1005
.
Spahn
,
P.
and
Reuter
,
R.
(
2013
).
A vertex model of Drosophila ventral furrow formation
.
PLoS ONE
8
,
e75051
.
Sweeton
,
D.
,
Parks
,
S.
,
Costa
,
M.
and
Wieschaus
,
E.
(
1991
).
Gastrulation in Drosophila: the formation of the ventral furrow and posterior midgut invaginations
.
Development
112
,
775
-
789
.
Vuong-Brender
,
T. T. K.
,
Amar
,
B. A.
,
Pontabry
,
J.
and
Labouesse
,
M.
(
2017
).
The interplay of stiffness and force anisotropies drive embryo elongation
.
eLife
6
,
e23866
.
Zienkiewicz
,
O. C.
and
Taylor
,
R. L.
(
1991
).
The Finite Element Method
.
McGraw-Hill
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information