A simplified mechanism for anisotropic constriction in Drosophila mesoderm

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.

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 L_x (length) and L_y (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 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 φ₀. 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.

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. L_x>> L_y, 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. L_x=L_y), 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/φ₀, Poisson's ratio σ and the length-to-width ratio of the contractile domain L_x/L_y. 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 φ₀/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 u_x and u_y components of the displacement, shown for one particular choice of parameters in Fig. 2C) for different φ₀/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 φ₀/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. u_x=b_xx, u_y=b_yy (with b_x,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).

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.

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⁻⁵. 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 (−r²/ξ²), 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⁻⁴.

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

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

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′.

In the case of fixed boundary conditions, we use a penalty technique to enforce u_j=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 A_pp=10³⁰ as well as F_p=0×10³⁰. Consequently, on the row of the stiffness matrix corresponding to that node, we have ⁠, thus leading to u_p=0.

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 u_s is a solution of Equation A.10, and if u_r is an eigenmode of the matrix A_ij with zero eigenvalue, then clearly u_h=u_s+u_r 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 (u_h). 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 thank Miriam Osterfield for helpful discussions.