Introduction

Tissue elongation is one of the basic processes of morphogenesis. Tissues can elongate by oriented cell divisions and growth [1, 2], or via local rearrangements of cells. The latter happens during epithelial “convergent extension” – a common motif of early development and organogenesis in many organisms – where an epithelium elongates along one axis while contracting along the perpendicular direction [3]. Epithelial convergent extension exemplifies the important role of cell and tissue mechanics in morphogenetic processes, which must be studied alongside developmental genetics and signaling. The fundamental question of developmental mechanics is how force generation is coordinated across cells to produce a coherent morphogenetic outcome.

The Drosophila embryo is one of the best-established models of developmental mechanics as it is ideal for live imaging and offers an extensive set of genetic tools [4]. Progress in live imaging and computational image analysis has produced remarkably quantitative data [5, 6]. In this work, we will take advantage of a previously published dataset which we will reanalyze and use as test-bed for theory development. Figure 1 recapitulates the basic quantitative features of germ band elongation (GBE) during Drosophila embryogenesis based on the light-sheet imaging data from [6, 7] which, thanks surface extraction [8] and cell segmentation and tracking [7], provides a global picture of tissue dynamics with cellular resolution. During Drosophila gastrulation, the embryonic blastoderm (an epithelial monolayer of about 6000 cells on the surface of the embryo) undergoes dramatic deformation that changes tissue topology and gives rise to the three germ layers [9]. Gastrulation starts (see Fig. 1A) with the formation of the Ventral Furrow (VF) which initiates the internalization of mesoderm, followed immediately by Germ Band Extension (GBE), which involves convergent extension of the lateral ectoderm (or germ band) and the flow of tissue from the ventral onto the dorsal side of the embryo [10–12] (see Movie 1). Simultaneous to GBE, the Posterior MidGut (PMG) moves from the posterior pole towards the anterior on the dorsal side of the embryo and invaginates progressively as it moves.

Figure 1:

VF and PMG invagination and GBE have been extensively studied, leading to the identification of relevant developmental patterning genes [10–13]. Live imaging has also uncovered the pertinent cell behavior of GBE, namely intercalation of neighboring cells in the lateral ectoderm [13–15]. The elementary step of this cell rearrangement process is called a T1 transition. It involves a quartet of cells where two cells lose contact and the two other cells gain contact (see Fig. 1B, bottom). The role of intercalations is highlighted by a “tissue tectonics” [16] analysis, which decomposes the tissue strain rate into cell-level contributions: cell-shape deformation and cell rearrangement (see Fig. 1B and Fig. S1). During VF invagination, cells in the lateral ectoderm are stretched along the DV axis with little rearrangement (Fig. S1A) Subsequently, during GBE, the lateral ectoderm (i.e. the germ band) extends (Fig. 1D) by cell rearrangements with very little cell elongation (Fig. 1C, bottom; and Fig. S1B and D). By contrast, the dorsal tissue (which will become the amnioserosa) deforms in the opposite orthogonal as the germ band, by a combination of cell shape changes and rearrangements (Fig. 1C, top; and Fig. S1C). Overall, we find that the majority of cell rearrangements happen via T1s (Fig. 1E and E’) with only a small number of “rosettes” (which involve more than 4-cells transiently sharing a vertex) being formed. This prevalence of T1s is consistent with previous analysis ([7, 17]) with the high coordination of rearrangements that is apparent in the lateral ectoderm (Fig. 1C, bottom).

While the above tissue tectonics analysis reveals the cell scale “kinetics” during GBE, it does not address the fundamental question of the underlying driving forces. Intercalations are associated with localized activity of the force generating protein non-muscle myosin II [12, 14, 15] (henceforth simply “myosin”). However, the relative contribution of such locally generated forces vs pulling by adjacent tissues, such as the invaginating PMG remain debated. Previous studies addressing this key question on the tissue level have come to conflicting conclusions [13, 17–22]. A second key question is how myosin dynamics is controlled on a cellular level and how it is orchestrated across cells to create a coherent global morphogenetic flow. Myosin recruitment depends on the expression of key developmental genes [11, 12, 14, 23], and, in addition, is subject to positive and negative mechanical feedback that depends on the stress [24, 25] and the rate of strain (rate of cell deformation) [26]. Despite considerable understanding of the genetic and cell-biological components involved in GBE, the relative roles of the genetic pre-pattern (top down) vs local “self-organization” of myosin activity via mechanical feedback loops (bottom up) are still not clear. All of these questions call for a coherent theoretical framework that allows interpreting and reconciling existing experimental findings.

Our approach is based on the assumption of force balance of stresses concentrated in the cell cortices. Intercellular adhesion effectively links the cytoskeletons of neighboring cells, so that a two dimensional epithelial sheet constitutes a transcellular mechanical network put under tension by myosin motors. This internal tension is revealed by laser ablation of cell interfaces, which causes rapid recoil [27]. Cortical force balance provides a direct link between mechanics and geometry: It allows inference of tensions from the angles at the vertices where interfaces meet [28–30], which we utilize to identify from images the stereotyped local geometry and tension dynamics associated with internally driven (active) and externally driven (passive) cell rearrangements (T1s).

The capacity of epithelial tissue to support tension makes it markedly different from passive fluids [31]. How one can reconcile the solid-like capacity of an epithelial monolayer to support tension (as well as shear stress) with its ability to change shape and rearrange internally like a fluid is intimately linked to the fundamental question of internal vs external driving. Our findings suggest that morphogenetic flow of epithelia can be understood as adiabatic deformation of cell array geometry controlled by changes in the internal state of tension. In other words, the tissue behaves as a plastically deforming active solid rather than a fluid. Based on these insights, we provide evidence and a minimal model for self-organization of the internally driven cell rearrangements via a local mechanical feedback mechanism. This model reproduces the experimentally observed dynamics on the level of cell quartets and forms the basis for a predictive tissue-scale model that is formulated in a companion paper [32]. Finally, we address how cellular behaviors (shape changes and intercalations) are coordinated across the tissue to generate coherent tissue flow. We will show that coordination of active T1 events among neighboring cells involves a characteristic pattern of cortical tensions which we quantify by introducing the “Local Tension Configuration” (LTC) order parameter. In the companion paper [32], we employ this order parameter to quantitatively compare tissue-scale simulations with experimental data on the cell scale. Taken together, our findings identify the dominant role of active internal tension in the lateral ectoderm in driving the embryo scale flow and suggest mechanical feedback as the mechanism for self-organization on the cell-scale.

Results

Force balance and cell geometry in an epithelial monolayer

We begin by laying out the assumptions and concepts that underlie our framework. Epithelial tissues are under internally generated tension, which is revealed by recoil in response to laser ablation. The timescale of this recoil (∼ 10 s) on the scale of cells is at least 10-fold faster than the timescale of local tissue flow [33, 34], so that the tissue can be regarded as being in approximate mechanical equilibrium. This suggests that the apparent tissue flow can be understood in terms of adiabatic remodelling of a quasistatic force balance network [30]. This view contrasts with regular fluid flow, where externally or internally generated forces are balanced by viscosity or substrate friction.

We further assume that mechanical stress in the epithelium is primarily carried by the adherens-junctional cytoskeleton, which resides on cell interfaces. This is supported by the observation that cell-cell interfaces in the Drosophila blastoderm are mostly straight. Cells are attached to their neighbors via adherens junctions that are linked to the junctional cytoskeleton in each cell (Fig. 2A) [35]. The myosin motors exert a contractile force on the actin fibers and thereby keep the cortex under active tension (T), which we refer to as “cortical tension”[36]. Together, the adherens junctions define a tissue-wide mechanical structure, capable of (i) generating locally controlled internal tension and (ii) adaptively remodelling its architecture. The Drosophila blastoderm lies on top of a fluid yolk [37] which exerts negligible drag forces on tissue motion on the surface [38], suggesting that all forces are balanced within the epithelial layer.

Figure 2:

Concentration of tension in interfaces tightly links force balance with the readily observable geometry of cell-cell interfaces. Force balance requires that the forces T_ab exerted by the cortical tensions sum to zero at each vertex where three (or more) interfaces meet. (Here, the indices a, b label the two cells that meet at a interface; see Fig. 2A.) These force-balance constraints relate the relative tension on cell interfaces at a vertex to the angles at which the interfaces meet [28]. This link of hardto-observe cortical tensions to the readily observable local cell geometry enables tension inference methods which have been extensively validated by computational robustness checks [29], direct comparison with measured laser ablation recoils [39] and by correlation with local myosin abundance [30].

Crucially, the link between force balance and the geometry of the cell array goes beyond inference and to the very basis of our proposed mechanism of tissue flow. Force balance implies that tension vectors meeting at a vertex sum to zero and hence form a closed triangle [31], as illustrated in Fig. 4B. The angles of this triangle are complementary to the corresponding angles at which the interfaces meet. (This becomes evident by rotating each triangle edge by 90°. Two angles are complimentary if they sum to 180°) As adjacent vertices share an interface, the corresponding tension triangles must share an edge: the tension triangles form a triangulation, i.e. fit together to form a tiling of the plane. The triangulation is dual to the cell array: For each cell, there is a corresponding vertex in the tension triangulation (Fig. 4B’). It reflects the fact that the force-balance conditions at neighboring vertices are not independent because they share the interface that connects them. (Dual force tilings go back to Maxwell [40] and have been applied to the statics of beam assemblies [41] and granular materials [42].) This triangulation establishes a geometric structure in tension space. Force balance requires that the angles at vertices in the physical cell array are complementary to those in the tension triangulation. This intimately links tension space and epithelial geometry in real space. Myosin-driven local changes in tension change the shapes of tension-space triangles and hence remodel the tension triangulation. The induced changes in the geometry of the cell array drive both the local rearrangement of cells and the global tissue flow. The remainder of the Results section will provide quantitative data analysis and modelling that will work out the above ideas and their implications.

Cell scale analysis

To map out the relative tensions in the tissue, we perform local tension inference for each interface relative to its four neighbors. Geometrically in tension space, the local tension configuration is represented by two adjoined tension triangles forming a kite (see Fig. 2B) and the inference is a simple application of the law of sines on the two triangles (see SI Sec. 1.4). Figures 2C and C’ show snapshots from this local tension inference in the dorsal ectoderm (amnioserosa) and the lateral ectoderm (germ band) respectively. Initially, relative tensions are close to unity throughout the embryo, since the cell array is approximately a hexagonal lattice, with vertex angles close to 120°. As gastrulation progresses, the cortical tensions change and one starts to see characteristic differences between the dorsal ectoderm and the lateral ectoderm. In the former, interfaces that contract remain under approximately constant tension while interfaces oriented parallel to the direction of tissue stretching (i.e. along the DV axis) extend and are under increasing tension (Fig. 2C). A very different picture emerges in the lateral ectoderm, where one observes an alternating pattern of high and low tensions before intercalations start (22 min). The high tension interfaces contract, leading to coordinated T1 transitions (30 min, cf. Fig. 2C’ and Movie 2). As GBE transitions from the fast phase to the slow phase at around 40 min, the pattern of tensions becomes more disordered. We will return to the pattern of local tension configurations below in the discussion of tissue scale dynamics.

The differences in the patterns of inferred tension in the amnioserosa compared to the lateral ectoderm suggest very distinct mechanisms for cell intercalations in these two tissue regions, matching the fact that their levels of cortical myosin are very different (high in the lateral ectoderm, low in the amnioserosa [18]). In the following, we first focus on intercalating cell quartets to quantitatively analyze these different mechanisms. A quantitative understanding of the cell-scale dynamics will then form the basis for bridging to the tissue scale.

Relative tension dynamics distinguishes active and passive intercalations

We identify all cell quartets which undergo neighbor exchanges (T1 processes), calculate the length and relative tension of all collapsing and emerging interfaces and align the data to the time of the neighbor exchange [43]. Pooling the data for each of the bands of cells colored in Fig. 1A, we find two distinct scenarios for ventro-lateral quartets and dorsal quartets (Fig. 2D and E; for a breakdown by individual bands, see Fig. S3).

The length dynamics of collapsing and extending interfaces in the amnioserosa and the lateral ectoderm is qualitatively similar (Fig. 2D and E). In the lateral ectoderm, there is a slight increase in interface length preceding contraction. This transient stretching is caused by the VF invagination that precedes GBE [26].

By contrast, the tension dynamics is markedly different between the amnioserosa and the lateral ectoderm. In the lateral ectoderm, the tension on the contracting edge grows non-linearly and reaches its maximum just before the neighbor exchange. In terms of the local cell geometry, this increasing relative tension reflects the fact that the angles facing away from the interface decrease as the interface contracts (Fig. 2F). Notably, the non-linearly increasing tension, concomitant with an accelerating rate of interface contraction, is evidence that positive tension feedback plays a role in myosin recruitment [24, 25, 44] and is in excellent agreement with predictions from a recent model where such feedback is a key ingredient [45]. From the data shown in Fig. 2D we can read off the average relative tension threshold ⟨T_crit ⟩= 1.530 ± 0.005 for interface collapse. As we will see further below, this threshold can be predicted from simple geometric considerations.

The correlation between increasing tension and interface contraction during active T1s can be used to predict active T1s. Indeed, plotting the time to interface collapse against the relative tension for quartets in the lateral ectoderm shows a clear negative correlation (Fig. 2H). Conversely, relative tensions below 1 correlate with interfaces that never collapse (Fig. 2J).

After the neighbor exchange, the relative tension on the new interface starts at a lower value and then continues decreasing back to 1, corresponding to vertex angles of 120°. Because Figs. 2D and E show tensions on the collapsing edge for t < 0 and on the emerging edge for t > 0, there is no reason why the tension should be continuous at t = 0. The apparent jump in relative tension is a consequence of geometry: Because the angles facing away from the interface are < 90° before the neighbor exchange, they are necessarily > 90° afterwards. This implies that a relative tension on the collapsing interface is necessarily followed by a tension on the new interface.

Let us now turn to T1s in the amnioserosa. Here, the relative tension on the inner edge remains almost constant near 1 prior to the neighbor exchange, i.e. the vertex angles remain close to 120° (see Fig. 2G). As a consequence of this tension homeostasis on collapsing interfaces, there is no correlation between relative tension and the time until the interface collapses in the amnioserosa (see Fig. S4). On the new interface emerging after the neighbor exchange, tension is high and remains constant for an extended period. Again, the tension apparent jump across the neighbor exchange is a consequence of geometry. Just before the neighbor exchange, the vertex angles are close to 120° which implies that the angles facing away from the interface are 60° after the neighbor exchange (see Fig. 2G). This corresponds to a relative tension of on the central interface. Just before the neighbor exchange the vertex angles are close to 120° which implies that the angles facing away from the interface are 60° after the neighbor exchange (see Fig. 2G). This corresponds to a relative tension of on the emerging central interface.

To understand the high tension observed on emerging interfaces, recall that tension inference only yields the total tension, but not on how this tension is generated in the cytoskeleton. On an extending junction, tension carried by passive crosslinkers will add to the tension generated by myosin. The passive tension rapidly relaxes as crosslinkers turn over, giving an effective viscoelastic relaxation timescale on the order of minutes [46]. This passive tension relaxation is a crucial ingredient in the model presented below. In the amnioserosa, the high tension is sustained for a longer time because the tissue there is continually getting stretched as the germ band contracts along the DV axis. Indeed, increased tension is also found in interfaces that start out DV oriented (cf. Fig. 2C).

Minimal model based on tension feedback reproduces lengthand tension-dynamics of T1s

Tension inference and the pooled analysis of T1 events has revealed the cortical tension dynamics on cell-cell interfaces during active and passive cell intercalations. The behavior of cortical tensions differs significantly between distinct spatial regions of the embryo and, on first glance, appears quite counter-intuitive. In particular, it is very different from that of springs or rubber bands: the length and tension of a spring are tied to one another by a constitutive relationship. By contrast, both the length and tension of cell-cell interfaces can change independently and are actively regulated by turnover of junctional proteins (such as actin, myosin and E-cadherin).

These observations call for a new modeling approach, where tension are not governed by constitutive relationships such as the typical area–perimeter elasticity of the vertex model [47, 48]. Instead, our model directly builds on the same assumptions that underlie the tension inference: dominant cortical tensions that are in adiabatic force balance. The dynamics resides in the changes in tension governed by mechanical feedback.

To keep the number of parameters and variables (degrees of freedom) in the model to a minimum, we consider a quartet of cells with identical shapes. This describes a representative quartet in a periodic cell array. Such a cell array is characterized by the three angles at each vertex, ϕ_i, and the three interface lengths 𝓁_i, with i ∈ {0, 1, 2} (see Fig. 3A). The vertex angles are determined by the relative cortical tensions T_i via the condition of force balance (see Fig. 4A). Motivated by the nonlinearly increasing tension observed on contracting interfaces (cf. Fig. 2D), we equip each interface with self-reinforcing tension dynamics

Figure 3:

Figure 4:

which models tension-induced myosin recruitment [24] (with excitability exponent n set to n = 4 in our simulation). The dynamics conserves the total active tension Σ_i T_i, corresponding to a finite myosin pool in each cell. (Other choices for fixing the overall scale of tension, e.g. via the triangle area, as well as well as other values of n do not change the results qualitatively; see Ref. [32].) The timescale of tension dynamics τ_T is fitted to the experimental data. Positive feedback drives the relative tension on the quartet’s central interface towards the tension threshold discussed above (see Fig. 3B) and thus drives active T1s. The interface that collapses is the one that starts with the highest initial tension, which we label with i = 0 by convention such that T₀ > T_1,2.

The tension dynamics described above determines the angles at vertices via adiabatic force balance but leaves the interface lengths 𝓁 as independent degrees of freedom. Further below, we will identify these degrees of freedom as so-called isogonal (angle preserving) modes, which are an immediate consequence of having cortical tensions governed by feedback instead of a constitutive relation [31]. To find the lengths 𝓁_i, we need to account for subdominant mechanical contributions from internal structures of the cell, such as medial myosin contractility [19, 49] the nucleus [50] and microtubules [51–53]. We account for these contributions through a “cell elastic energy” in terms of the deviation of the cell “shape tensor” S_C = Σ_i(e_i ⊗ e_i)/𝓁_i from a target shape tensor S₀:

Here, e_i are the interface vectors with lengths | e_i | = 𝓁_i (see Fig. 3A) and ⊗ denotes the tensor product (We choose this shape tensor because it is linear in interface length, meaning that it does not change if an interface is subdivided by inserting an additional vertex). Importantly, we assume the scale of E_C to be much smaller than the elastic energy due to cortical tensions. Because of this separation of scales, the angles at vertices are fixed by tensions, and the relaxation of the cell elastic energy only affects isogonal modes.

The coefficients λ and µ control the cell’s resistance to isotropic compression/dilation and shear deformations, respectively (see SI Sec. 3 for details). Notably, this elastic energy does not engender an “energy barrier” for T1 transitions as is found in vertex models with area–perimeter elasticity [54] (see also Fig. S10D). Therefore, our simulations do not require stochastic fluctuations to drive T1s.

The model described thus far captures the dynamics up to the neighbor exchange (time < 0 in Fig. 3E and F) and qualitatively reproduces the dynamics of interface length and tension observed in the experimental data (cf. Fig. 2D). Interface contraction during an active T1 is driven by active remodeling of the force balance geometry: Myosin recruitment increases the cortical tension and thus drives it out of force balance with the external tension from adjacent interfaces. As a result, the interfaces remodel, changing the angles at vertices until force balance is reestablished. Positive tension feedback continually causes further myosin recruitment, which in turn drives further contraction until the interface has fully collapsed.

How does the subsequent resolution and interface elongation work? On the emerging interface, there is typically little myosin [15, 55]. Therefore, active contractility cannot balance the external tension exerted on the interface by its neighbors: The new interface is thrown out of active force balance. This force imbalance naturally leads to the extension of the new edge. As the interface extends, passive tension bearing elements of the cytoskeleton, such as crosslinkers, get loaded until force balance is reestablished. To account for this, we split the total cortical tension T into the myosin-borne, active tension T_myo and the passive tension T_p. Cortical remodeling due to turnover of crosslinkers will gradually relax this passive tension. We account for this by an exponential decay with a characteristic remodeling rate τ_p. This passive tension effectively represents Maxwell-type viscoelasticity of the extending interface, illustrated by the spring and dashpot in Fig. 3D. Such visco-elastic behavior on a scale of minutes has been reported in experimental data from the Drosophila embryo [46].

Completing the model requires to understand what sets the initial motor protein level on an emerging interface, which will set the initial condition for the tension dynamics on that interface. We propose a myosin “handover” mechanism based on the assumption that the myosin level along the cortex within each cell changes continuously along the interfaces and across vertices [56] (see Fig. 3D). (Force balance requires that the total tension, the sum of the tensions in the two abutting cortices, is uniform along an interface. However, the individual tensions on either side can be non-uniform, as the resulting traction forces are exchanged via E-cadherin molecules linking the cells along the interface [28].) Importantly, this “handover” model predicts that the myosin level (not the total tension) on the newly formed interface is always lower than that on the collapsing interface, in agreement with experimental observations [15, 55].

Simulating this model, we find good qualitative agreement between the singlequartet model and experimental data, both for the tension dynamics (compare Fig. 3E,F and Fig. 2D,E) and for the tension–isogonal decomposition (compare Fig. S10A,A’ and Fig. 4C,C’). Movie 3 shows the simulation of an active T1 transition for a symmetric quartet. In the absence of positive tension feedback, the model successfully describes passive intercalations when the displacements of cell centroids are prescribed to mimic external stresses (Fig. 3F; see SI Sec. 3 for details).

Notably, in our model there is no additional active mechanism for interface extension. Instead, interface extension is a consequence of the fundamental temporal asymmetry of T1 processes: high myosin levels on the collapsing interface versus low myosin levels on the emerging interface. No additional active ingredients are required.

Here, we considered a minimal model for an idealized regular lattice geometry here. In a real tissue, there will inevitably be disorder which will impact how cells can coordinate with one another. In the companion paper [32], we investigate the role of such disorder in tissue-scale simulations comprising many cells. For the remainder of this manuscript, we return to the analysis of experimental data, going from the cell scale to the tissue scale.

A tension–isogonal decomposition quantifies active vs passive contributions to tissue deformation

The above analysis and minimal model have revealed the distinct mechanical forces that drive active and passive cell intercalations. These forces manifested themselves through their effect on the geometry of the cell array. Dominant cortical tensions constrain the angles at vertices. These angle constraints are represented by the tension triangulation. Interface lengths can change collectively while keeping those angles fixed via isogonal modes [31], as illustrated in Fig. 4A. These isogonal modes therefore represent cell and tissue deformations that take place under constant cortical tensions. They can be caused, for instance, by external forces and cell-internal elasticity (medial myosin, stiffness of the nucleus, intermediate filaments). It is important to keep in mind that this association of isogonal deformations with nonjunctional stresses is based on the assumption that (active) cortical stresses are the dominant source of stress in the tissue (separation of scales).

These above considerations suggest that we can use the observed geometry of the cell array to decompose cell and tissue deformations into two distinct components: (i) deformations of the tension triangulation reflecting changing cortical tensions, and (ii) isogonal deformations reflecting all other mechanical forces. We will call this a tension–isogonal decomposition. This decomposition only assumes the dominance of cortical tensions over other stresses in the tissue and in particular is independent of the specific form of the cell elastic energy Eq. (2) we used in the minimal model above.

The key elements of the tension–isogonal decomposition are illustrated in Fig. 4B. Changes in cortical tension are reflected in deformations of the tension triangulation (i.e. changes of the angles at vertices). To each tension triangulation, we can construct a reference state for the polygonal cell array – compatible with the force balance constraints – from the tension triangulation using a Voronoi tessellation.

Isogonal deformations then transform this reference state into the physical cell array. In practice, we do not need to explicitly construct the reference cell array. Instead, we utilize the fact that the physical cell array can be quantified by the triangulation defined by the cells’ centroids [57]. The isogonal modes are calculated as the transformations that deform the tension triangulation into the centroidal triangulation.

We first apply the tension–isogonal decomposition to intercalating quartets. For a given quartet of cells labelled i = 1−4, we quantify the physical shape in terms of the centroid positions c_i by calculating the moment of inertia tensor (see Fig. 4C). The tension reference shape is given by the tension vectors , which are obtained from the interface angles using tension inference and are rotated by 90° relative to the interfaces (see Fig. 2B). The ratios of eigenvalues of the moment-of-inertia tensors defines the shape aspect ratios, illustrated by ellipses in Fig. 4C and C’. The quartet’s isogonal deformation tensor I_quartet transforms the cell centroids into the tension vertices (see Fig. S7A for a labelled sketch). Given and c_i, the equations can be (approximately) solved for the best fitting I_quartet using least squares (an exact solution is not possible in general because the equations are overdetermined; see SI Sec. 1.5.) The constant global scale factor η ≈ 4.2 µm translates units from relative tensions (a.u.) to length (µm) and can be calculated from the average interface length ⟨ 𝓁(0) ⟩ ≈ 3.5 µm at the first timepoint using the Voronoi–Delaunay duality as explained in SI Sec. 1.5.

In Figs. 4C and C’, we plot the quartet shape aspect ratios and the eigenvalue ratio of the isogonal deformation tensor against the length of the quartet’s inner edge The latter serves as a pseudo-time parameterization of the T1 process; plots against real time are shown in Figs. S7B and B’. In the lateral ectoderm, the tensiontriangle deformation leads the quartet shape deformation, indicating that tension dynamics drives the intercalation (Fig. 4C). This confirms our conclusion from the relative tension analysis (Fig. 2D). The isogonal strain is approximately constant and, in fact, slightly counteracts the deformation. This is a consequence of the ventral furrow invagination which pulls on the lateral ectoderm cells and passively stretches them along the DV axis. A very different picture emerges in the passively deforming amnioserosa (Fig. 4C’). Here, the tension mode is constant—indicating tension homeostasis—while the cell quartet deformation is entirely stored in the isogonal mode. While the behaviors observed in the lateral ectoderm and in the amnioserosa represent the extreme cases, being purely tension driven (active) in the former and purely isogonal (passive) in the latter, there is continuous spectrum in between. Intermediate scenarios where both cortical tensions and external stresses contribute to tissue deformation are conceivable and can be detected and analyzed using the tension–isogonal decomposition.

Taken together, the tension–isogonal decomposition quantitatively distinguishes between active and passive intercalations. Tensions are controlled locally, while the isogonal modes accommodate external, non-local forces, such as pulling by adjacent tissue. We associate isogonal deformations with passive deformations because active stresses in the lateral ectoderm are generated at interfaces. In systems where active stresses are generated in the cell’s interior, e.g. due to nuclear migration [58] or due to intracellular actin cables [59], the isogonal mode is actively controlled. Another example for this is apical constriction during VF invagination which is driven by medial myosin pulses that drive isogonal cell contractions [31].

Tension space geometry sets the threshold for interface collapse

A puzzle that has remained open so far is the threshold for interface collapse during a cell neighbor exchange. In the lateral ectoderm, where the isogonal strain is approximately constant (Fig. 4C), the neighbor exchange takes place at a critical relative tension of 1.530 ± 0.005 on average (cf. Fig. 2D). By contrast, in the amnioserosa, where T1s take place under constant tension (Fig. 2E), there is an average critical aspect ratio ∼ 3 of the isogonal deformation at the point of neighbor exchange (see Fig. 4C’).

We now explain these numbers in terms of the geometric tension–isogonal decomposition introduced in the previous section. We start by decomposing the length 𝓁 of a cell quartet’s central interface into a reference length 𝓁_ref ({T_i}), determined purely in terms of the local relative tensions, and an isogonal length Δ𝓁_iso depending on the quartet’s isogonal strain:

A neighbor exchange takes place when 𝓁 shrinks to zero. Thus, setting 𝓁 = 0 in Eq. (3) defines the threshold values of the local tensions {T_i} and the isogonal length Δ𝓁_iso where T1s take place (see green line in Fig. 4E). It is immediately apparent that the interface collapse can be driven by the isogonal mode, by changes in tensions, or by a mixture of both.

Let us start with the passive T1s observed near the dorsal pole. In this case, the relative tensions are homeostatic (cf. Fig. 2E), so 𝓁_ref ≈ ⟨𝓁(0)⟩ ≈ 3.5 µm remains constant. The neighbor exchange therefore happens when Δ𝓁_iso = ⟨ − 𝓁(0) ⟩ (see blue arrow in Fig. 4E). By a geometric construction (Fig. S7D) one finds that the corresponding isogonal deformation has an aspect ratio of 3, in good agreement with observations (Fig. 4C’).

To find the tension threshold for active T1s, we calculate 𝓁_ref using the Voronoi construction based on the vertices of the tension triangulation (see SI Sec. 1.5 for details). The length of a quartet’s central interface in the Voronoi reference state is determined by the pair of tension triangles corresponding to the cell quartet as illustrated in Fig. 4D. For simplicity, we consider the case of two identical, isosceles tension triangles which is fully characterized by the relative tension T = 2 sin(ϕ/2) on the central interface, with ϕ the angle in the tension triangle. (The general case of asymmetric tensions is discussed in the companion paper [32], where we show that the T1 threshold is shifted towards stronger tension anisotropy with increasing asymmetry of the tension configurations.) In the symmetric case, the Voronoi construction for the reference length yields 𝓁_ref = ηT cot(ϕ). In the absence of isogonal deformations (Δ𝓁_iso = 0), the interface length vanishes when ϕ = π/2, implying the critical relative tension . One can understand the value for the critical relative tension directly from the circumcircle construction of the Voronoi tessellation as shown in Fig. 4D’. 𝓁_ref vanishes when the two adjacent tension triangles share the same circumcircle. In the case of isosceles triangles, this implies that they form a square such that we can read off the critical relative tension from the length of the diagonal.

In general, the condition for the neighbor exchange defines the “T1 threshold” in the Δ𝓁_iso-T plane (see Fig. 4E). This diagram quantifies how active and passive forces interact to drive cell intercalations. Active tension dynamics and passive isogonal strain appear as orthogonal ways to drive T1s, as illustrated by the red and blue arrows.

In the following section, we will generalize the tension–isogonal decomposition to the tissue scale. From this analysis one can then estimate a local average value for Δ𝓁_iso and then read off the T1 threshold from Fig. 4E. This will allow us to explain the value T_crit ≈ 1.53.

Tissue scale

Our analysis so far has focused on the mechanism of individual T1 transitions in cell quartets. To bridge the gap from cell quartets to tissue-scale convergent–extension flow, we need to address three key questions: (i) How are active T1s oriented? (ii) Which regions of the tissue deform actively due to internally generated local tension and which ones yield passively to stresses created by adjacent regions? (iii) How are active T1s coordinated across cells, so that different interfaces do not “attempt” to execute incompatible T1s? In the following, we will address these three questions by building on the tools (tension inference and tension–isogonal decomposition) we employed above to analyze intercalating cell quartets.

Initial anisotropy of tension directs flow

Convergent–extension flow of tissue requires that the T1 transitions are oriented. If cortical tensions are regulated by positive feedback, the tissue-scale tension anisotropy sets the orientation the interfaces which will collapse, and hence, the direction of tissue flow. Therefore, tissue-scale anisotropy of active tension is central to drive and orient convergent–extension flow [18, 55, 60, 61].

We asses tension anisotropy by locally averaging the anisotropy of inferred tensions T_ij at a given individual vertices as illustrated in Fig. 5A (see SI Sec. 1.4 for details). The inferred tension anisotropy (Fig. 5B) and its time course (Fig. 5C) are consistent with previously published experimental observations: During VF invagination, DV oriented interfaces in the lateral ectoderm are stretched causing myosin recruitment [26]. Remarkably, we find strong DV alignment of tension anisotropy in the trunk already before VF invagination (Fig. 5D, 0 min). This confirms our hypothesis that tension anisotropy is set up in the initial condition. As GBE progresses, the DV alignment of tension anisotropy decreases (Fig. 5C and Fig. 5D, 37 min). Numerical simulations presented in the companion paper [32] reproduce this loss of global tension anisotropy and show that it is responsible for the slow down of GBE after two-fold extension.

Figure 5:

Isogonal strain identifies regions of passive tissue deformation

We started our investigation by a “tissue tectonics” analysis [16] that identifies the contributions of cell intercalations and cell shape changes to tissue deformations. This kinematic quantification in itself is not informative of the mechanical forces driving epithelial dynamics. To quantify the relative contributions of active (local) vs passive (non-local) forces acting on intercalating cell quartets we introduced the tension–isogonal decomposition (Fig. 4C). Just as a tissue tectonics analysis, this decomposition can be performed on the scale of entire tissue regions without the need to track cells and identify T1 events. Specifically, we exploit the fact that the isogonal strain tensor can be calculated for each individual triplet of cells that meet at a vertex (see Fig. 6A and SI Sec. 1.5 for details). Locally averaging over nearby vertices then yields the tissue-scale isogonal deformation (see Fig. 6B and B’). At the end of VF invagination (25 min), significant isogonal strain has built up adjacent to the VF indicating that the tissue there is passively stretched by the invaginating VF (see purple shaded region in Fig. 6B). The lateral ectoderm further dorsal also accumulates some isogonal extension along the DV axis (green shaded region, see time traces in Fig. 6E). Specifically, from the DV-DV component of the isogonal strain tensor, we can estimate the average isogonal contribution to the length of DV oriented interfaces ⟨Δ𝓁_iso⟩ ≈ 1.9 µm in the lateral ectoderm. Using this value, one can read off the predicted relative tension threshold T_crit ≈ 1.53 in from the green curve in Fig. 4. This value is in excellent agreement with the value found by tension inference (cf. Fig. 2D), thus validating the model for active T1s. An immediate prediction from this model is that abolishing VF invagination will eliminate isogonal stretching of the lateral ectoderm and thereby shift the T1 threshold to . To test this prediction, we analyzed a light-sheet recording of a snail mutant, which lacks the VF as illustrated in Fig. 7A [62]. This movie is taken from the Drosophila morphodynamic atlas [6]. As expected, we find almost no isogonal strain in the lateral ectoderm of the snail mutant (see Fig. S6). Analyzing the dynamics of relative tension in intercalating quartets, we find (see Fig. 7B), confirming that tension–isogonal decomposition faithfully captures the interplay of internally generated stresses and external forces during T1 transitions.

Figure 6:

Figure 7:

During GBE, the isogonal strain near the ventral furrow relaxes back to nearly isotropic, suggesting that the passive deformation stored in the isogonal mode is elastic rather than plastic (Fig. 6B’ and C). At the same time, the actively contracting lateral ectoderm stretches the dorsal tissue (amnioserosa) along the DV axis. This passive deformation of the amnioserosa is consistent with the transition of cells from columnar to squamous [7] and low myosin density in the amnioserosa [18]. There is also significant isogonal strain just anterior of the PMG (highlighted by the dashed outline in Fig. 6B’), as well as on the lateral side of the PMG, suggesting that the hyperbolic (convergent-extension) tissue flow near the dorsal pole [18] exerts a pulling force on the PMG. Quantitative re-analysis of previously published cauterization experiments [19] provides independent evidence of such a pulling force (see SI Sec.4).

Taken together, isogonal strain identifies regions where passive tissue deformation is caused by mechanical coupling across the tissue. Importantly, the absence of isogonal strain along the AP axis in the germ band provides strong evidence that the main driver of GBE is internally generated stress, rather than an external pulling force from the PMG. This conclusion from cell-scale analysis is in line with quantitative tissue-scale analysis of previously published experimental data which we address in the Discussion.

A geometric order parameter quantifies local tension configurations

Tissue-scale tension anisotropy is not enough to drive coherent “parallel” T1s. Such coordination requires an alternating pattern of high and low tensions as seen in Fig. 8A (center; cf. Fig. 2C’). In contrast, when several connected interfaces are under high tension, they form a tension cable, contraction of which leads to rosette formation [24, 63]. Tension patterns do not appear globally with perfect regularity, but rather as local motifs that can be distinguished based on the tension configuration at individual vertices. These configurations correspond to characteristic triangle shapes in the tension triangulation (see Fig. 8B): The elementary motif of the alternating tension pattern is a “tension bridge” where a high tension interface is surrounded by low tension interfaces (see Fig. 8A, right), corresponding to an obtuse triangle (Fig. 8B, top right). In contrast, an acute tension triangle (Fig. 8F, bottom right) corresponds to the local motif of a tension cable, where neighboring interfaces along the cable are under high tension while those transverse to the cable are under low tension (see Fig. 8A, left). We can therefore use the tension triangle shape to define a local tension configuration (LTC) parameter. The LTCparameter space (i.e. the space of triangle shapes), shown in Fig. 8B, is spanned by two axes, measuring how elongated and how acute or obtuse the tension triangle is (see Ref. [32], Fig. S8 and SI Sec. 1.7 for details). These axes correspond to the magnitude of local tension anisotropy and the cable vs bridge character of the local tension configuration, respectively.

Figure 8:

Figure 9:

LTC order choreographs T1s

The geometric condition for neighbor exchanges discussed above (cf. Fig. 4E) defines a threshold in the LTC-parameter space (dashed line in Fig. 8B, see [32] for details). When the shapes of a pair of adjacent tension triangles cross this threshold, the interface corresponding to the triangles’ shared edge collapses, giving rise to a T1 transition. From the position of the T1 threshold in the LTC-parameter space, it is immediately evident that for tension cables (acute tension triangles), the tension anisotropy required to cross the T1 threshold is higher compared to tension bridges. Tension cables are therefore less efficient at driving T1s: they require stronger anisotropy of tension to cause an interface collapse.

Two-dimensional histograms of the LTC-parameter distribution in the lateral ectoderm show that the fraction tension of bridges increases before the onset of GBE (see Fig. 8C). The bridge fraction reaches its maximum around 25 min, just at the onset of GBE (see Fig. 8D). The corresponding alternating pattern of tensions is clearly visible in Fig. 8A (center). This suggests that the pattern of tensions on the (sub-)cellular scale is biologically relevant to choreograph parallel T1s. As the tension anisotropy increases (cf. Fig. 5C), the local tension configurations eventually hit the T1 threshold where a cell neighbor exchange happens, causing a reconfiguration of local tensions. Therefore, tension configurations beyond the T1 threshold are strongly suppressed (hatched region in Fig. 8C). We observe that as cell intercalations take place in the lateral ectoderm, the fraction of tension bridges decreases and the fraction of tension cables increases (see Fig. 8G, 50 min; and Fig. 8H). In fact, the distribution of tension triangle shapes appears to approach that of a random Delaunay triangulation (see SI Fig. S9), a completely disordered distribution. This suggests that the triangle edge flips can be statistically understood as a random “mixing” process that generically causes a loss of LTC order and global alignment of tension anisotropy. This result, as well as the results of the tissue–scale analysis more generally, can be reproduced by tissue–scale simulations of our minimal model, as shown in the companion paper Ref. [32].

Discussion

Cortical tensions drive and constrain tissue deformations on the cellular scale

We have developed a novel perspective on tissue dynamics based on the underlying principle of force balance. Central to our approach is the balanced network of active forces generated and transmitted by cells. In epithelia where contractility in the adherens-junctional cortex is the strongest source of stress, the force network takes the form of a triangulation in tension space [31] which fixes the angles at tri-cellular vertices in the tissue. Tissue deformations are driven by the adiabatic transformation of this force-balance geometry, allowing the epithelium to rearrange like a fluid, while supporting internal tension like a solid. Alternatively, one can think of this behavior as a form of plasticity where the internal “reference” structure of the material can be actively remodelled, resulting in spatial displacement of material points. However, the focus on the force balance geometry associated with the internal structure is more useful as it holds the key to understanding possible mechanisms of active control. We find two mechanically distinct modes of tissue deformation: a cortical-tension driven mode and an isogonal (angle preserving) mode that is unconstrained by the cortical force-balance geometry. The former represents active remodeling of interfaces, e.g. by myosin recruitment, while the latter accounts for passive tissue deformations that are controlled by external (non-local) forces and cell shape elasticity. The tension–isogonal decomposition quantifies the contributions of these two modes in experimental data based on cell geometry alone. This allows us to disentangle whether deformations are due to locally or non-locally generated forces, i.e. whether they are active or passive. Importantly, tension–isogonal decomposition is based on the same assumptions as tension inference which has been validated extensively in previous literature [29, 30, 38, 39].

In the early Drosophila embryo, we observe both active and passive cell rearrangements (T1 events). Generally, both modes—tension-driven and isogonal— interact and contribute to tissue deformations at the same time. We have found that ventral furrow invagination isogonally stretches the lateral ectoderm which increases the tension threshold for active T1s there. This predicts that the local relative tension threshold will be lower (namely ) in mutants which lack a ventral furrow (twist and snail [62]). Indeed, analysis of a snail mutant embryo confirmed this prediction and thus validates that tension–isogonal decomposition captures the interplay of internal and external forces acting on a tissue.

Internal and external contributions to germ band extension

Where the forces driving tissue deformations during morphogenesis originate is a fundamental question of developmental biology. In the context of Drosophila gastrulation, the question is whether the germ band elongates due to internally generated stresses [13, 15, 18, 20, 26, 55] or due to external pulling by the posterior midgut invagination [19, 21, 22, 64] has been intensively debated. Experimental evidence shows that both processes contribute, making a nuanced, quantitative analysis necessary. We therefore re-analyzed microscopy data from previously published cauterization and mutant experiments [19] to quantify tissue flow (see SI Sec. 4). The results from this quantitative analysis indicate that forces generated in the germ band contribute significantly to tissue flow. This conclusion is further supported by the observations from mutants where posterior midgut invagination is disrupted (fog, torso-like, scab, corkscrew, ksr). In these mutants, the germ band buckles forming ectopic folds [65, 66] or twists into a corkscrew shape [67, 68] as it extends, pointing towards a buckling instability characteristic of internally driven extensile flows [69, 70]. This suggests that the main effect of PMG invagination on the germ band lies not not in creating pulling forces, but rather in “making room” to allow for its orderly extension.

However, these tissue-scale observations only provide circumstantial evidence for internally driven GBE. To conclusively settle the debate, evidence from the cell scale – where the forces are generated – is needed. Such evidence is provided by our tension–isogonal analysis, which yields a fine grained picture with various regions of active and passive deformation (see Figs. 6B and B’). It clearly shows that tissue deformation in the germ band is driven by internal remodeling of tensions and therefore active. In stark contrast, the amnioserosa, and the tissue just anterior and lateral of the invaginating posterior midgut deforms passively.

Cells orchestrate tissue flow by self-organizing in tension space

The internally driven nature of germ-band elongation flow immediately raises the question of how force generation is coordinated across cells in order to drive coherent tissue flow. In other words, how can local behavior of cells orchestrate global tissue flow? Cells exert active stresses on each other and, at the same time, constantly sense their mechanical environment [71–73]. Mechanical stresses and strains can propagate over long distances and contain information about the tissue geometry (e.g. in the form of hoop stresses [74]). In an epithelium dominated by cortical tensions, this mechanical environment forms a “tension space”, linked to physical space via force balance. Tension space takes the form of a triangulation which allows an intuitive visualization. The angles in the cell array are fixed by complementary angles in the tension triangulation. Thus, the tensions are geometric dials cells can directly sense and control.

We found that cells can control their configuration in tension space by defining local dynamics of cortical tension. In experimental data from gastrulating Drosophila embryos, we identified two distinct behaviors: (i) amplification of tension on the interfaces that are already most tense, suggestive of positive tension feedback (observed in lateral ectoderm) and (ii) apparent tension homeostasis (observed in the amnioserosa). Mechanical homeostasis is found in various systems and organisms [75–77]. We hypothesize that tension homeostasis allows the amnioserosa to undergo large cell deformations while maintaining tissue integrity [76, 78]. Tension feedback, on the other hand, continuously modifies local force balance, driving the change in cell and tissue geometry. To the extent that cortical responses can be controlled by spatially modulated gene expression, evolution thus has the means to define a program of non-trivial spatio-temporal dynamics of tissue during morphogenesis. Evidence for positive feedback is provided by the nonlinearly increasing inferred tension on contracting interfaces (Fig. 2D) and by laser ablation experiments in earlier work [24]. The underlying molecular mechanisms are yet to be identified and might rely on the catch-bond behavior of myosin [79, 80] or mechano-sensitive binding of α-catenin [81, 82].

We have found that T1 events can be explained quantitatively a simple geometric criterion which defines a “T1 threshold” in terms of the local tension configuration and the local isogonal strain. A minimal mathematical model of an intercalating cell quartet demonstrates how positive tension feedback can drive the tensions towards this T1 threshold and thus generate active T1s. In contrast, passive T1s result from external forces changing the isogonal modes until a cell neighbor exchange occurs, while tensions (and thus interface angles) remain constant.

Our mathematical model for an intercalating cell quartet also sheds light on the long-standing puzzle of how interfaces extend after a neighbor exchange, even though cortical myosin can only exert contractile forces. We show that no additional active mechanism is necessary. Our minimal model captures interface extension as a purely passive process. In our model, interface extension during both active and passive intercalations is driven by the tension on adjacent interfaces which exceeds the low active contractility of the new edge. The low myosin level on the emerging interface is predicted by a simple model for the myosin “handover” during a neighbor exchange. In tissue-scale simulations [32], the model reproduces experimental observations (interface extension even if tissue extension is blocked, generation of irregular cell shapes) that were previously taken as evidence for actively driven interface extension [19, 83].

Our findings on the dynamics of interface length and tension during active T1s (Fig. 2D) are similar to another recently proposed model [45]. This model extends the classical “vertex model” with perimeter–area elasticity by active cortical contractility. Like in our model, active T1s are driven by a positive feedback loop of cortical tension. Since our model is formulated in adiabatic force balance, it requires fewer parameters than the one of Ref. [45] and might be obtained from the latter in the limit of fast mechanical relaxation and dominant active tensions.

Other recent theoretical work has interpreted Drosophila germ band extension as a result of tissue fluidization [84], i.e. the loss of resistance to shear stress due to vanishing of cortical tensions [47, 85]. However, laser ablation experiments [19, 27] show that cortical tensions remain finite during the morphogenetic process. In fact, interfaces under vanishing would generically buckle into a wiggly shape (as observed, for instance, in the amnioserosa during dorsal closure [86]). This is not observed in the Drosophila ectoderm during gastrulation where interfaces instead appears taught straight, serving as a simple visual indicator of finite tension. We have therefore investigated tissue mechanics in the regime where the tissue internally generated tension at all times, and hence is a solid rather than a fluid. Tissue flow occurs as a result of internally driven plastic rearrangements, and does not require or imply tissue fluidization. The capacity of an epithelial tissue to internally rearrange in a controlled mannerwhich may be thought of as active plasticity falls beyond the fluidization paradigm. The results of our analyses on the cell and tissue scale corroborate this point of view. We return to the question of fluid vs solid behavior in the companion paper [32].

Order in local tension configurations coordinates active T1s

We found a striking alternating pattern of cortical tensions that explains how T1s are coordinated between neighboring cells. The prevalence of this organizing motif is quantified by an order parameter for local tension configurations (LTC). LTC order precedes intercalations in lateral ectoderm and the loss of LTC order correlates with the slowdown of tissue flow at the end of GBE. This raises two important questions: First, how does LTC order emerge? Second, is loss of LTC order the cause for the termination GBE? These questions are addressed in a companion paper [32] using a tissue-scale mathematical model based on the quartet-scale model introduced here. We show that LTC order is driven by positive tension feedback and requires an initial regular hexagonal packing of cells. Intercalations disrupt this regular hexagonal packing and thus lead to the loss of LTC order as GBE progresses. Together with a degradation of coherent tension orientation, loss of LTC order causes T1s to become incoherent and therefore ineffective at driving tissue flow. Taken together, this implies that convergent–extension flow is selflimiting and that the final tissue shape is encoded in the initial degree hexagonal order and tension anisotropy. Moreover, we predict that disrupting the hexagonal packing of nuclei prior to cellularization will prevent the emergence of LTC order and thus cause slower GBE.

Genetic patterning and initial anisotropy

Because Drosophila embryo is a closed surface, convergent extension of the germ band must be compensated by an orthogonal deformation of the dorsal tissue, the amnioserosa (see Fig. 6B’). To achieve that, mechanical properties of the tissue must be modulated along the dorso-ventral axis: Positive tension feedback driving active rearrangements in the lateral ectoderm, and tension homeostasis in the amnioserosa allowing it to yield to external forces while maintaining tissue integrity. This highlights the importance of the DV-axis patterning. Recent work [26] has shown that mechanical feedback loops in the embryo are patterned along the DV axis. Notably, the DV patterning system is conserved across species [87]. Our work also highlights a second interaction between tissue flow and DV patterning: because the coherent T1s we observe do not mix cells in the tissue (Fig. 1A), distinct DV fates (e.g. neural and surface ectoderm) remain clearly demarcated, a clear biological necessity.

In addition to the DV gradient, tension anisotropy along the DV axis is required to orient the GBE flow along the AP axis. This anisotropy might be due to mechanics (pulling of the ventral furrow [26] and hoop stresses resulting from ellipsoidal embryo geometry [74], since we observe anisotropy already prior to VF formation), or indeed due to AP-striped genetic patterning [88, 89]. Thanks to the positive tension feedback, a weak initial anisotropy is sufficient to bias the direction of tissue elongation. This explains why twist and snail mutants, which have significantly reduced initial myosin anisotropy, still extend their germ band, albeit at a slower rate [26].

Our findings show how genetic patterning and self-organization are interconnected during GBE. Genetic patterning provides tissue scale input by demarcating distinct tissue regions with different mechanical properties. Self-organization via positive tension feedback orchestrates myosin contractility on the cellular scale and drives T1s. In the companion paper [32], we show that mechanical self-organization via tension feedback also allows cells to coordinate their behaviors (such as active T1s) on the tissue scale to give rise to coherent morphogenetic flow. As of yet, no genetic mechanism has been found that produces cell scale coordination – as manifested, for instance, in the alternating pattern of high and low tensions (Fig. 8A, center) – from predetermined genetic patterns alone. Moreover, many interfaces in the lateral ectoderm that undergo T1 transitions at late stages of GBE are initially oriented along the AP, rather than the DV axis and only rotate into DV alignment briefly before they collapse (see SI Fig. S2). This shows that an initially set up genetic pattern is not sufficient to explain the coordination of cell rearrangements, in line with evidence from tissue scale analysis [74]. Instead, the emerging overarching picture is that of morphogenesis driven by an interplay of bottom-up local self-organization controlled by top-down genetic patterning that sets initial conditions and modulates parameters on the tissue scale [90, 91]. Indeed, this interplay has recently been shown to underlie posterior midgut invagination, where genetic patterning initiates and channels a wave that mechanically propagates due to a local feedback mechanism acting on myosin [64].

Going forward, we expect that geometric insight into tension space will provide a new perspective on many different systems [92], for example Xenopus neural tube formation [93], amniote primitive streak formation [94–96], wing disk elongation [97, 98], the sea urchin archenteron [99], or kidney tubule elongation [100]. Further, it will be interesting to study intercalations driven by actin “protrusions” [101]: they too are governed by cortical force balance on the cellular level, making tension inference and the tension–isogonal decomposition applicable.

Acknowledgements

We thank Noah Mitchell, Sebastian Streichan, and Arthur Hernandez for stimulating discussions and careful reading of the manuscript. We would also like to thank the anonymous reviewers for excellent feedback on the manuscript. FB acknowledges support of the GBMF post-doctoral fellowship (under grant #2919). NHC was supported by NIGMS R35-GM138203 and NSF PHY:1707973. MFL was support by NSF Career grant No. PHY-2047140. EFW acknowledges support from Howard Hughes Medical Institute. BIS acknowledges support via NSF PHY:1707973 and NSF PHY:2210612.

Supplementary Material

1 Analysis of experimental data

1.1 Data source

Whole embryo cell segmentation and tracking data for a WT embryo was obtained from the repository DOI:10.6084/m9.figshare.18551420.v3 deposited with Ref. [7]. We analyzed the dataset with id number 1620 since it had the highest time resolution (15 s) and covered the longest time period (50 min), starting ca. 7 min before the onset of ventral furrow invagination.

For the analysis of the snail mutant, we used data from the Drosophila morphodynamic atlas [6]. This data is available on the data repository Ref. [102]. We used recording No. 202207111700 of a halo snail[IIG05]; CAAX-mCherry fly. For cell segmentation and tracking, we used the Fiji plugin TissueAnalyzer [103] and custom Mathematica code. For analysis of T1 transitions and isogonal strain, we excluded regions that form folds. Note that in recording No. 202207111700 the one dorsal fold ectopically extends all the way to the ventral side.

1.2 Tissue tectonics analysis

Based on segmented and tracked cells, the tissue strain rate and cell deformation rate can be calculated via the methods described in Ref. [16]; see also [104] and [57]. In brief, the local tissue strain rate is calculated from the centroid positions c_i of the cells in a small neighborhood N of a cell (we use nearest neighbors). One looks for a linear transformation that maps the positions of the cell centroids at one timepoint to those at the next one. Because the resulting equations are overdetermined, an approximate solution minimizing the error is obtained using the least-squares algorithm. The linear transformation is then further decomposed into a symmetric part, describing local shear and area changes, and an antisymmetric part, describing the tissue rotation. Because the Drosophila ectoderm is a surface in three dimensional space, we need to calculate these quantities in the local tangent plane. The cell deformation rate is calculated from the rate of change of the cell’s moment of inertia tensor in the frame that co-rotates with the tissue.

Figure S1A–B’ show the tissue and cell deformation rates during VF invagination (A and A’) and during GBE (B and B’). Figure S1C shows the cell elongation as measured by the Beltrami coefficient µ = |a − b| /(a + b), where a and b are the half-axes of a best-fit ellipse to the cell.

Local tangent plane and coordinate frames

The analysis of local tissue geometry (including tissue tectonics, tension inference, and tension–isogonal decomposition) is performed in the local tangent plane. The local tangent plane is determined by the local normal vector n which we find by a singular value decomposition of the point set of interest (vertex coordinates for a single cell, cell centroids for a patch of cells). The normal vector is read-off as the basis vector corresponding to the smallest singular value.

Depending on the context, we use three different coordinate bases within the tangent plane: (i) The “Eulerian” (or “lab frame”) basis formed by the azimuthal vector along the DV axis and its orthogonal complement , along the AP axis. The Eulerian frame is used to plot orientation and strain fields in the flat map projections (“pullbacks”) of the tissue (e.g. Fig. 4D and Fig. 8B). (ii) A “Lagrangian” (or “tissue frame”) basis which co-rotates with the tissue is used to quantify the DV component of the isogonal deformation tensor (Fig. 4E). This co-rotating frame is defined by tracking a patch of cells and fitting a linear transformation A(t) to the displacement of cell centroids relative to the patch’s center of mass. The linear transformation is constrained to transform the initial into the final normal vector A.n(0) = n(t). A polar decomposition A(t) = S(t).R(t) then yields the local rotation matrix R(t) for the tissue patch. The co-rotated basis vectors are obtained by applying R to the Eulerian basis vectors at time zero: . (iii) For the tension–isogonal decomposition of intercalating cell quartets (Fig. 4C and C’ and Fig. S7), we define a co-moving orthogonal basis based on the cell centroids: The basis vectors are defined as the pair of orthogonal vectors that best fit the displacement vectors that connect opposite cells in the quartet. Before the neighbor exchange, the first basis vector is aligned with the two cells that touch and the second basis vector is aligned with the cells that don’t touch. After the neighbor exchange, the roles are flipped. This ensures that the basis vectors vary smoothly through the neighbor exchange.

1.3 Identification and classification of intercalation events

Cell interface collapse and emergence events can be identified by tracking cell contacts. We start by identifying all pairs of cells that share a interface at some timepoint. For each pair, one can then identify the time periods where the cells are in contact. This allows one to identify interface collapse and emergence events. We exclude cases where a cell pair loses contact because one of the cells invaginates or divides and cases where cells come into contact because cells in between them invaginate.

Some cell pairs transiently lose or gain contact. For the interface collapse rate shown in Fig. 1E, we only account for the first collapse event for cell pairs that are not in contact at the last frame, i.e. for those that permanently lose contact. Similarly, for interface emergence events, we only account for the last event where a pair of cells not in contact at the first frame comes into contact.

Spatial distribution of interface collapses

From the tracked interface collapses, we can count the number of interface collapses that each cell is involved in, giving us a spatially resolved map of intercalations (Fig. S1D). In the lateral ectoderm (regions 4–8 in Fig. S3A), we find 2.58 interface collapses per cell. From this, we can estimate the total tissue extension by cell rearrangements. In a hexagonal lattice, synchronous T1s cause two interface collapses per cell and elongate the tissue by a factor , assuming that the cells do not elongate. (The tissue also contracts by a factor along the orthogonal axis.) For 2.58 interface collapses per cell, therefore cause tissue elongation by a factor , in good agreement with the total elongation of the germ band (cf. Fig. 1D).

In contrast to the lateral ectoderm where cells don’t elongate coherently, cell elongation contributes significantly to tissue deformation near the dorsal pole (see Fig. S1C). This is a generic feature of passive T1s driven by isogonal deformations (see Fig. S7D).

Orientation of collapsing interfaces

The convergent–extension flows in the lateral ectoderm and near the dorsal pole (amnioserosa) are oriented perpendicular to each other (see Fig. S1). We therefore also expect perpendicular orientation of the collapsing interfaces in these two regions. We quantify interface orientation of collapsing edges at the blastoderm stage (i.e. before the onset of tissue flow, t = 0 min) as shown in Fig. S2A and C. In the ectoderm, interfaces that collapse before t = 42 min, i.e. within the first ca. 15 min of GBE are predominantly DV oriented at t = 0 min (see Fig. S2B). By contrast, interfaces near the dorsal pole predominantly oriented along AP. The orientation of interfaces in the lateral ectoderm that collapse after t = 42 min is initially not strongly biased along the DV axis Fig. S2C and D. However, during GBE, these interfaces progressively align with the DV co-rotated DV axis before they eventually start contracting (see Fig. S2C’ and D; the coordinate frame co-rotating with the tissue is described in Sec. 1.2).

Tracking quartets and identifying rosettes

The above procedure for identifying intercalation events is very fast and easy to implement. However, it does not distinguish between T1s and higher order intercalation events (rosettes). Moreover, it cannot easily be determined which interface collapse corresponds to which interface emergence.

Figure S1:

Figure S2:

We therefore implemented an alternative algorithm to identify intercalations, by tracking quartets of cells. This is possible by identifying cycles of length four in the neighborhood graph of the cell array at each time point. We exclude cases where the four cells surround one or more cells that are not part of the cycle and cases where one of the cells is surrounded by the three others. Notably, the quartets don’t necessarily exist for the entire time series. Namely, when one of the “inner” interfaces in the quartet, other than the central interface, collapses, the quartet is destroyed.

A T1 event happens when the opposite pairs of cells in the quartet exchange neighbors (i.e. one pair loses contact and the other one gains it). In some cases, the neighbor exchange is not instantaneous, but instead there is a fourfold vertex in the center of the quartet for an extended period. Cases where the central interface and one of its neighbors are simultaneously shorter than a threshold length (0.5 µm) are classified as rosettes and not counted as T1s. The threshold value 0.5 µm corresponds to the average vertex-positions fluctuations on the timescale of 1 min. Higher order rosettes can be “decomposed” into overlapping 5-fold rosettes and therefore can be identified by finding such overlaps.

Using this procedure, we identified 2231 T1 events, 498 five-fold rosettes, 40 six-fold rosettes, 4 seven-fold rosettes and no higher order rosettes (see Fig. 1E’). Because the number of rosettes is small compared to the number of T1s, we have focused on T1 events for the subsequent analysis. Expressed in terms of a participation ratio of interface contractions in T1s vs rosettes, we find a ratio 2231 : (2 · 498 + 3 · 40 + 4 · 4) ≈ 2 : 1, matching observations from previous studies [7, 17].

1.4 Local tension inference

Consider the triplet of cells that meet at a tri-cellular vertex and label them i = 1, 2, 3. Denote the unit vectors pointing along the three interfaces that meet at the vertex by e_ij, where i, j are the labels of the pair of cells that share the respective interface. Force balance implies that

where T_ij are the cortical tensions. Multiplying the above equation with ê₁₂ and ê₂₃ (where ê denotes the vector orthogonal to e and with the same magnitude) yields two equations

Here, α_ijk denotes the angle subtended by the interfaces (ij) and (jk). An alternative, more geometric way to arrive at these equations is by applying the law of sines to the tension triangle, using that its angles are complementary to those at the vertex and sin(α − π) = sin(α).

Eqs. (5) equations do not constrain the overall level of tension. We fix this by setting the average tension to 1 which gives the solution These equations are solved by

with the normalization factor a = Σ _(ijk) sin α_ijk.

The above calculation can be generalized to an arbitrary number of vertices. Notably, the equations become overconstrained once the set of vertices contains entire cells [31]. This is because force balance implies that the tension vectors T_ij = T_ije_ij form a triangulation (see below).

Here, we limit ourselves to local tension inference for the central interface in a cell quartet (see Fig. 2A). In this case, the equations are not overconstrained. Labelling the cells 1–4, the equations read

We solve these equations under the normalization condition that the average of the “outer tensions” T₁₂, T₂₃, T₃₄ and T₁₄ is 1. The inferred value for “inner” cortical tension T₁₃ then is that relative to the average tension of its neighbors. We call this the relative tension of a interface throughout the manuscript.

The use of local tension inference is justified because we are interested in local cell behaviors, namely rearrangements. Tension inference is also most robust on the local level, since this is where force balance, the underlying physical determinant of the link between mechanics and geometry, resides. In global tension inference, spurious large scale gradients can appear when small deviations from local force balance accumulate over large distances. For this reason we do not try to infer tension gradients in the magnitude of tension across the embryo, e.g. between the amnioserosa and the lateral ectoderm.

Analysis of collapsing and emerging interfaces

For the analysis presented in Fig. 2D and E, collapsing and emerging interfaces were identified separately as described in Sec. 1.3. All events before 12 min were excluded from the analysis.

When a interface is short, the interface angles become highly susceptible to fluctuations of the vertex positions. We therefore excluded all timepoints where the central interface or one of its neighbors is shorter than a threshold length (0.75 µm).

Figure S3:

This value was chosen because the vertex positions fluctuate about 0.5 µm between timepoints. These fluctuations are partially due to intrinsic fluctuations in cortical tension and partially result from limitations in vertex localization during segmentation because of limited optical resolution.

Time series of interface length and relative tension are then aligned such that the interface collapse (respectively emergence) is at t = 0. Figure S3 shows the pooled data resolved by the position of the quartet along the DV axis.

1.5 Tension-isogonal decomposition

Because the angles in the tension triangulation are complementary to the angle at vertices in the cell array, one can find a unique tension triangulation (up to a global scale factor) for a given physical cell array in force balance (see dashed gray arrow in Fig. S7A). However, the reverse is not true. For a given tension triangulation, there is a continuous family of cell arrays that share the same tension triangulation and that are related by isogonal deformations which change the interface lengths collectively while preserving the angles. (Mathematically, there are fewer triangulation vertices than cell vertices, so that the triangulation leaves one isogonal mode degree of freedom per cell undetermined. We will use this in Sec. 2 to define an “isogonal potential” that parametrizes the isogonal modes.)

Figure S4:

The isogonal degrees of freedom reflect tissue deformations that are not constrained by force balance of cortical tensions. This means that weaker contributions to stress, e.g. from intermediate filaments, microtubules and the nucleus in the cell interior, as well as externally applied stresses will drive the isogonal modes. This has two important consequences: (i) For the analysis of experimental data, it implies that we can decompose the tissue deformations into active, tension-driven and passive, isogonal contributions. (ii) For modeling, it implies that we need to specify a cell shape elasticity to fix the isogonal modes.

In the following, we discuss the details of the tension–isogonal decomposition. To achieve a well-defined decomposition of tissue deformation, we need to specify a reference tessellation that is uniquely determined by the tension triangulation. There will then be a unique isogonal deformation that deforms the reference configuration into the physical cell array. Fortunately, we do not need to explicitly construct the reference configuration. Rather, we use the triangulation formed by the cell centroids to quantify the physical tissue shape [57]. (This implies to some degree of coarse-graining, as information about the non-affine displacements of individual vertices is lost.) We can then obtain the isogonal deformation by comparing the tension triangulation to the centroidal triangulation.

Single-vertex level

To make this construction concrete, we first consider a single tri-cellular vertex, with the associated tension triangle and the centroids c_i. The indices i, j = 1, 2, 3 label the three cells that meet at the vertex (see Fig. S5A’). The tension triangle vectors are obtained by local tension inference based on the angles at the vertex (see Sec. 1.4). The orientation of the tension vectors is determined by demanding that they are at right angles to the cell-cell interfaces e_ij. To fix the tension scale, we normalize the area of the tension triangle to unity. Fixing the tension triangle area is justified because it maintains the total area of the tension triangulation such that cell dilation/contraction are purely isogonal. The numerical value of the tension triangle area is immaterial, since it can be absorbed in the scale factor η.

We define the local isogonal deformation tensor I_vertex such that it transforms the tension triangle into the centroidal triangle. Mathematically, this definition is written as the equations . Since there are two independent displacement vectors in tension space and physical space each, these equations represent four independent constraints which uniquely determine the four components of the 2×2 matrix ηI_vertex. The scale factor η converts units of tension (a.u.) into units of length. How we fix η is described below.

Figure S5:

The isogonal deformation tensor I_vertex can be decomposed into a symmetric strain tensor U and a rotation R using polar decomposition UR = I_vertex. The appearance of a rotational component in the isogonal deformation might appear surprising at first glance, since isogonal modes leave angles at vertices invariant. However, the centroidal triangle can be sheared relative to the tension triangle by a simple shear which is a composition of pure shear with a rotation.

Fixing the scale factor η

To fix η we use that the tension triangles are initially close to equilateral (as evidenced by the angle distribution peaked around π/3, see Fig. 8L). Recall that we have normalized the area of the tension triangle to unity by convention. For an equilateral tension triangle, the edge length T is uniquely determined by the area, giving T = 2/3^1/4 ≈ 1.52. To find η, we compare this value to the average edge length of the centroidal triangles, i.e. the average centroid distance between adjacent cells, which we find to be 6.38 µm. The ratio of these two numbers gives a scale factor η ≈ 4.2 µm. Note that the scale factor η≈ 4.2 µm minimizes the isogonal strain at t = 0 (see Fig. S5B), indicating that at early times, the cell array geometry closely matches Voronoi dual of the tension triangulation. Indeed, we can also calculate η from the initial average cell-cell interface length ⟨𝓁(0)⟩ ≈ 3.6 µm. To this end, we use that the Voronoi dual to an equilateral triangulation with edge length T = 2/3^1/4 is a hexagonal lattice with edge length , giving η ≈ 4.1. Finally, minimizing ⟨||I||⟩ at a reference time is provides an alternative way to determine η in cases where the initial triangulation is not as uniform and regular as in the Drosophila blastoderm.

Coarse-grained isogonal strain

To find the tissue scale isogonal strain, we coarse grain the isogonal deformation tensor on a square grid in the flat map projection (“pullback”). The grid spacing is set to 20 µm, such that each grid facet contains ca 10–15 cells (grid facets near the poles contain less cells because of the distortion caused by the pullback projection). We calculate I_vertex for all tri-cellular vertices in a grid facet and take the average to find I_facet.

Figures 4D and S5B–D (top panels) show the isogonal strain tensor I_facet − 𝕀, where 𝕀 is the identity matrix. The bottom panels in Fig. S5B–D show the individual, single-vertex strain tensors I_facet − 𝕀 at each vertex. Crosses indicate the principal eigenvectors of the strain tensors. The bar length is proportional to the eigenvalue. Green (magenta) indicates a positive (negative) eigenvalue, i.e. contraction (elongation).

To calculate the average isogonal deformation in different tissue regions (Fig. 4E), we first transform the single-vertex isogonal deformation tensor to the local corotating (“Lagrangian”) frame (see Sec. 1.2). This is important because the germ band rotates as it moves over the posterior pole during elongation.

To estimate the critical tension for T1s in the lateral ectoderm, we need the average isogonal contribution to the length of DV-oriented interfaces, ⟨Δ𝓁_iso⟩ (see SI Sec. 1.5). At the end of VF invagination (t = 25 min), we find ⟨ (I_vertex)₂₂ ⟩ ≈ 1.538 ± 0.006, where the index 2 denotes the DV component and the average is taken over cells in the lateral ectoderm (region shaded in blue in Fig. S5C). We thus have ⟨ Δ𝓁_iso ⟩ = (I₂₂ − 1) ⟨ 𝓁(0) ⟩ ≈ 1.9 µm for DV-oriented interfaces, where ⟨ 𝓁(0) ⟩ ≈ 3.5 µm is the average initial interface length.

Figure S6 shows the coarse-grained isogonal strain in a snail mutant, where the ventral furrow is abolished. Without the ventral furrow invagination pulling on the lateral ectoderm, no isogonal strain is found in the lateral ectoderm before and during germ-band extension. This validates that isogonal strain reliably indicates tissue deformation driven by external forces, and, vice versa, the absence of isogonal strain indicates the lack of external driving forces.

Filtering of degenerate tension triangles

When one of the angles at a vertex is close to π, the tension triangle becomes highly elongated. In fact, sometimes one of the angles at a vertex exceeds π due to fluctuations and inaccuracies in detection. In the tension inference, this leads to negative tensions. The (nearly) degenerate tension triangles, in turn lead to extreme apparent local isogonal deformations. We therefore exclude nearly degenerate tension triangles whose perimeter is larger than 8 after normalizing the area to 1 (ca. 2% of the tension triangles are filtered out this way).

Figure S6:

Tension–isogonal decomposition for cell quartets

The calculation for cell quartets (vertex pairs) is analogous to the case of single vertices. The two tension triangles associated to the inner vertices of the cell quartet form a “kite” (see Fig. S7A). We normalize tension such that the average area of the tension triangles is one. The quartet’s isogonal deformation tensor I_quartet is defined such that it transforms the vectors T_ij into the corresponding centroidal displacement vectors c_j − c_i. Here, we only use the vectors around the perimeter of the tension and centroidal kite. This ensures that the isogonal deformation tensor varies continuously through the cell neighbor exchange where the inner edge in the kite is flipped. As in the single-vertex case, we find I_quartet using least squares.

To quantify the deformations of the tension kite and the physical cell quartet, we define the shape tensors

where the index pairs (ij) run over cell pairs going around the quartet, i.e. (ij) = (12), (23), (34), (41). Notably, with these definitions, one has the relation

Remarks on the Voronoi construction for the reference cell array

For the tension–isogonal decomposition, we have circumvented the explicit construction of a reference cell array by directly comparing the tension triangulation to the centroidal triangulation of the physical cells. This approach is particularly useful because a local isogonal deformation tensor can be obtained on the level of a single tri-cellular vertex.

The explicit construction of the reference cell array allows us understand the geometry of cell-neighbor exchanges and find a criterion for when they occur (see main text for a discussion simple symmetric case and the companion paper [32] for the general case). We define the reference tesselation using a generalized Voronoi construction based on the tension triangulation. The vertices of this generalized Voronoi tessellation are given by the circumcircle centers of the triangles in the triangulation. The generalized Voronoi tessellation is only free of self-interesections when the triangulation fulfills the Delaunay criterion that opposite angles in adjacent triangles must sum to a value smaller than π. As we will see, in the following section, this provides a purely geometric criterion for neighbor exchanges in the absence of isogonal strain. In the presence of isogonal strain, the reference can have self intersections (i.e. negative interface lengths). This is similar to how the reference state encoding plastic deformations in finite strain theory isn’t necessarily geometrically compatible [105].

Figure S7:

In general, the centroids of the generalized Voronoi cells will not coincide with the vertices of the tension triangulation. This means that the definition of the isogonal deformation tensor based on the tension triangulation vertices is not exactly identical to that based on an explicit Voronoi construction. However, the two agree for a periodic lattice of identical tension triangles. In this case, the vertices of the tension triangulation coincide exactly with the centroids of the generalized Voronoi tessellation. In other words, using the tension triangulation vertices to define isogonal deformations amounts to approximating tissue as locally lattice-like, which is valid when spatial gradients are small on the cell scale.

In passing, we note that the generalized Voronoi tessellation is not the only possible reference state. Any tessellation that is uniquely defined by a given tension triangulation is in principle a valid reference state. The Voronoi tessellation is a useful choice for the following analysis because the Voronoi cell shapes behave similar to the actual cell shapes and because it is defined purely geometrically. Other choices of the reference state might be of use depending on the particular question and when one has more detailed information about the mechanical properties of the cells.

1.6 Local tension anisotropy

From the tension vectors T_ij at a vertex, we calculate the “tension anisotropy” tensor

where the sum runs over pairs from the set of the three cells that meet a the vertex. The normalization factor is chosen such that Tr 𝕋 = 1. The difference of the eigenvalues of 𝕋 quantifies the magnitude of tension anisotropy. When the interface angles at the vertex are all identical to 120° the tension anisotropy vanishes. For non-identical vertex angles the tension is locally anisotropic and the dominant eigenvector of 𝕋 points along the principal axis of stress. Anisotropic tension also implies that the tension triangle is non-equilateral. The tension triangle’s axis of elongation is perpendicular to the principal axis of tension because the tension triangle vectors are rotated by 90° relative to the interfaces (cf. Fig. 2B).

In passing, we note that the tension anisotropy tensor is closely related, but not exactly equal, to the stress tensor σ. Instead it has the form of a metric tensor and thus describes the local geometry of tension space. Specifically, the 𝕋 measures the extrinsic anisotropy of tension and is complementary to the intrinsic shape tensor that we will introduce in the next section.

1.7 Local tension configuration

In the main text, we argued that local configuration of tensions at a vertex is characterized by the shape of the tension triangle. Specifically, acute triangles correspond to tension cables, while obtuse triangles correspond to tension bridges (see Fig. S8A).

Figure S8:

Whether a triangle is obtuse or acute is an intrinsic property, i.e. it is independent of the position, orientation and size of the triangle in the embedding space. The intrinsic shape of a triangle is given by its angles. Since these angles sum to π, we can represent the space of triangle shapes in barycentric coordinates (see Fig. S8B). The center of the barycentric shape space corresponds to an equilateral triangle, while the edges correspond to triangles where one angle is zero and the corners to triangles where two angles are zero. Because the order of labeling the angles in the triangle is arbitrary, the shape space is invariant under permutations of the angles and a single fundamental domain (highlighted in Fig. S8B) is representative.

In the triangle shape space, we have colored each point according to how anisotropic and how acute-vs-obtuse it is. We already know that we can read off the magnitude of anisotropy from the eigenvalues of the tensor 𝕋, defined in the previous subsection. To quantify the shape (acute vs obtuse), we define a second tensor

where the indices I, J ∈ {1, 2, 3} are short hand labels for the three interfaces that meet a tri-cellular vertex. Because Σ_I T_J = 0 (force balance), this tensor has the null vector (1, 1, 1)^T. To get rid of this nullspace, we go to barycentric coordinates by defining the pair of vectors, τ_1,2 that we pack into a matrix

These vectors have the defining property that they are orthogonal to one another and to the null vector (1, 1, 1)^T of 𝒮. Now we we can define the shape tensor

Observe that

which implies that 𝕋 and 𝕊 share the same eigenvalues. In fact we can perform a singular value decomposition of 𝔗

where R(ϕ) is the rotation matrix with angle ϕ and Λ = diag (λ₁, λ₂) is a diagonal matrix of singular values whose squares are the eigenvalues of 𝕋 and 𝕊. Because we normalized the tension vectors, . ϕ is the angle of the dominant eigenvector in physical space, which gives the orientation of the principal axis of stress. As we will see below, the angle ψ indicates whether the triangle is obtuse or acute. The reasoning to define the rotation based on ψ/2 will become apparent below when we represent the shape tensor by a complex number.

The shape tensor 𝕊 has only two degrees of freedom, the difference of eigenvalues (which it shares with the tensor 𝕋) and the angle ψ that sets the orientation of its two orthogonal eigenvectors. We can represent 𝕊 by a complex number

where is the magnitude of tension anisotropy and the phase ψ is the “shape space angle” defined by Eq. (16).

Permutations of the edge indices α correspond to sign changes and rotations of ψ by integer multiples of 2π/3. To mod out these group actions (which correspond to the symmetries of the barycentric shape space Fig. S8B), we define

Figure S9:

where x mod_d n = (x+d mod n) − d ∈ [− d, n − d] is the modulo operation with offset d. This shape parameter informs how acute or obtuse a tension triangle is. For an acute isosceles triangle, , whereas for an obtuse isosceles triangle . For non-isosceles triangles, takes values in [0, π] quantifying the relative degree of acuteness vs obtuseness. For an equilateral triangle, |Ψ| = 0 and is not defined. The axes |Ψ| and span the space of triangle shapes shown in Fig. S8.

As an alternative to the above mathematical derivation of the shape parameter, we briefly discuss a more intuitive parameter as follows: For a given tension triangle, label its edges in ascending order of length (i.e. T₂ and T₃ will second highest and highest tension, respectively. We can now define a the parameter

which quantifies the relative tension difference between the two largest tensions at a vertex. For a tension cable, there are two large tensions which are similar in magnitude, implying that p is small. For a tension bridge, there is one high tension edge meeting two low tension edges so p will be closer to 1. The maximal value p = 1 is reached for T₁ = T₂ = T₃/2. Remarkably, to very good approximation .

Random Delaunay triangulation

As a reference for the triangle shape distributions, we generated random Delaunay triangulations based on the Ginibre random point process [106]. A small sample of such a random triangulation is shown in Fig. S9A. The histograms in the triangle shape space Fig. S9B shows good agreement between the Ginibre-based random Delaunay triangulation and the distribution of tension triangle shapes at the end of GBE. Moreover, the marginalized angle distribution (black dashed line), closely matches the experimentally observed distribution (Fig. S9D).

These findings suggest that the local coordination of cortical tensions is lost during GBE and that the local tension configurations approach a maximally disordered state.

2 Isogonal shear

As discussed in the main text (and further elaborated in Ref. [32]) our simulations show that a non-zero cell-level shear modulus is required for extension via active T1s. Here, we analyze the relation of cell-level and tissue-level shear modulus. Crucially, because of the dominance of cortical tensions, a tissue patch in our model will respond to externally applied forces via an isogonal deformation. We first show that an isogonal deformation can create tissue-scale pure shear (on the level of cell centroid displacement – the transformation of cell vertices is necessarily non-affine to preserve angles). Then we show that these pure shear modes correspond to the response of the tissue to external force by computing the Hessian of our cell elastic energy in the subspace spanned by isogonal deformations, and measure the shear modulus.

To show that isogonal modes can create pure shear, not just dilation/contraction of cells, we make use of the isogonal mode parametrization introduced in Ref. [31]. It assigns an isogonal “potential” Θ_i to each cell, and calculates the cell displacements from the Θ_i and the edge tension vectors. In the following, we will show that a constant gradient in the isogonal potential generates a uniform translation in real space. By integration, this implies that a quadratic spatial profile of the isogonal potential creates a pure shear.

Let us identify the real space edge unit vectors by the two adjacent cells e_ij = −e_ji and denote the corresponding tensions as T_ij. Then the tension vectors form a triangulation, where â denotes the normal vector to a, i.e. â.a = 0 and ||â|| = ||a||.

The isogonal displacement u_ijk of the real space vertices r_ijk (identified by the three adjacent cells) is given by

where is the area of the tension triangle (ijk).

First, observe that the uniform isogonal mode Θ_i = const. has no effect on the vertex positions because

Now we aim to show that a constant gradient in Θ_i drives a uniform displacement of the r_ijk. Specifically, by uniform gradient we mean Θ_i = t_i.a, i.e. a linear gradient in the tension space (t_i is the position of the tension triangulation vertex corresponding to cell i, such that ). To show that the displacement in real space is uniform, it is enough to show that two adjacent vertices are displaced identically. By induction, this implies that all displacements are identical. It is therefore sufficient to consider a quartet of cells (i = 1−−4), corresponding to a “kite” in tension space (note that â.b is identical to the wedge product a ∧ b.)

Because a constant can be arbitrarily added to all Θ_i, we can set Θ₁ = 0 and thus have for i = 2, 3, 4. The displacements now read

To show that these displacements are identical, we project them onto two conveniently chosen, linearly independent vectors, namely and . For the latter we find

where we used that and applied the definition of S_ijk.

Projecting (22) and (23) onto gives

To show equality of these two right-hand sides, we use that given and we can find a and substitute the result into . We start by “expanding the identity”

Explicitly writing out the inverse matrix then gives

With this, we find the relation

where we used T_ij = − T_jk to flip the indices on T₁₄. Comparing to (26) and (27) now shows the identity of their RHSs.

Taken together, we have shown that

Because and are linearly independent, it follows that u₁₂₃ = u₁₃₄. QED. We just showed that a constant gradient in Θ_i corresponds to a uniform displacement of the real space vertices r_ijk. We can therefore think of Θ(t) as a “potential” for the isogonal displacement field: u ≈ ∇_tΘ, where the approximation is valid for slowly varying gradients and exact for constant gradients. The gradient ∇_t is taken in tension space because the function Θ(t_i) = Θ_i is defined on the vertices of the tension triangulation t_i.

A pure shear aligned with the coordinate axes is given by a displacement field u(r) = ε diag (1, −1).r and is therefore generated (approximately) by a quadratic isogonal potential Θ = ε t^T.diag (−1, 1).t = ε [(t¹)² − (t²)²].

3 Simulation methods

In the following, we explain the simulation details for Fig. 3, where we present a model of a single intercalating quartet. The code for the simulations shown is available online: https://github.com/nikolas-claussen/geometric_basis_of_convergent_extension_simulation.

Our model is based on two assumptions: adiabatic force balance, and dominance of cortical tensions. The first assumption means that we obtain the instantaneous cell geometry by solving force-balance equations. The second assumption means the we do so in two steps: first, we determine the interface angles from the cortical tension force balance. The remaining forces in the system are much weaker than cortical tensions and only affect the residual degrees of freedom, the isogonal modes.

In this work, we only consider a symmetric lattice of identical tension triangles, so a number of simplifications occur (the general case of a disordered cell array is the subject of the companion paper [107]). In this case, the cell shape is completely specified by the three vertex angles ϕ_i and edge lengths 𝓁_i at each vertex, with i ∈ {0, 1, 2 }. Because the cell array is a symmetric lattice, the edge lengths 𝓁_i can be changed without affecting the angles and therefore parameterize the isogonal degrees of freedom.

We assume that the dynamics of the tensions are determined by a mechanosensitive feedback loop implementing positive tension feedback. To account for the dynamics of a junction post T1, each junction at a vertex is characterized by the total tensions T_i, and the passive tension T_p,i. The tensions at a vertex evolve according to:

where n > 1 and τ_T is a time scale converting simulation time into minutes (fit to the data). The passive tension is 0 on all junctions except those that were newly created by a T1 transition. The initial value of T_p on those interfaces is discussed below. Given the tensions, we first calculate the angles from the tensions using the law of sines as , where R is the tension triangle circumradius.

Cell shape energy

To determine the 𝓁_i, we minimize the an elastic energy based on the cell shape tensor:

where i runs over the edges of the cell 𝒞, and e_i is the vector pointing from one vertex of the interface to the other. Note that in this shape tensor, each edge contributes linearly in length to the cell shape. This means that artificially subdividing an edge has no effect on the cell shape tensor. This makes sense if we assume that the elasticity we aim to model using S resides in the cell interior (incompressibility, microtubules, intermediate filaments, nucleus). The shape tensor can also be defined using vectors from the cell centroid to its vertices (in the lattice, the two definitions are equivalent). An alternative definition of the elastic energy using instead gives broadly similar results (although interface collapse happens more abruptly, because of the higher order non-linearity).

We assign a reference shape S₀. An isotropic S₀ ∝ Id favours equilateral hexagons. We chose a slightly anisotropic reference shape S₀ (15% anisotropy) to model the isogonal stretching caused by the ventral furrow before onset of GBE. This anisotropy sets the angle at which the inner interface collapses, 𝓁₀ = 0. The experimental value of the collapse tension was therefore used to fit the S₀-anisotropy.

We define the cell elastic energy via

with bulk modulus λ and shear modulus µ. We used a shear/bulk ration of µ/λ = 0.1. Because of the separation of scales between cortical tensions and elastic energy baked into the model, the absolute values of λ, µ are irrelevant. Further, for the case of active T1s, the results do not depend on µ/λ, as long as µ > 0.

The edge lengths 𝓁_i can now be determined by minimizing the elastic energy Eq. (34) w.r.t. the 𝓁_i, for the angles determined by the tension dynamics.

Comparison to area-perimeter elastic energy

Strikingly, within the singlequartet model, the “shape strain” S_C −S₀ can always be set to 0 by choice of the edge lengths, and the energy E_C = 0 throughout. This can already be seen from a degreeof-freedom count (three 𝓁_i for the three independent components of S_C − S₀). The elastic energy therefore acts only on the isogonal modes (i.e. ) and there is no energy barrier for intercalations. Consequently, there is no need for noise to drive intercalations in our model. By contrast, for the widely-used areaperimeter elastic energy E = (A−A₀)² +(P −P₀)² (where A, P are the cell area and perimeter, and A₀, P₀ their target values) [54], there exists an energy barrier, and the inner interface 𝓁₀ only collapses when ϕ₀ = π. This is shown in Fig. S10B. Note that the area-perimeter energy is a special case because of the geometric incompatibility of area and perimeter constraints. Combining area elasticity with shear elasticity based on the shape tensor, E = (A − A₀)² + µTr[(S_C − S₀)²], (or perimeter elasticity with cell shape bulk elasticity) leads to similar results as Eq. (34). Note also that because of the degree-of-freedom count, the system is under-specified if the shear modulus is 0, foreshadowing the fact that without shear modulus, no convergenceextension takes place [107].

Myosin handover and passive tension

When an interface collapses, 𝓁_i = 0, the tension triangulation is modified topologically: the edge corresponding to the collapsed interfaces is replaced by one corresponding to the new interface (triangulation “flip”). To complete our model, we must specify the initial conditions of this new edge.

As shown in Fig. 3D, we propose a myosin handover mechanism to explain the extension of the new interface post intercalation. A interface with cortical tension T is comprised of the adherens-junctional actomyosin cortex of the two adjacent cells, which are coupled mechanically via adherens junctions. Under force balance, the total tension T has to be constant along the cortex, but the individual tensions on either side can be non-uniform, as the resulting traction forces are exchanged via adherens junctions. In the following, we assume as a first order approximation that the level of active tension (i.e. myosin concentration) varies linearly along an interface (similar calculations have been performed in Ref. [56] to calculate interfacial shear stress). This allows us to geometrically obtain the myosin concentration at the individual cortices that will form the two juxtaposed sides of the new interface (see Fig 3D).

Consider the tensions T₀, T₁, T₂ at a vertex, and let 0 be the interface about to collapse. The three cells that meet at the vertex will be referred to by (01), (12), (21) (where e.g. (01) is the cell abutting interfaces 0 and 1). Let m₀₁, m₁₂, m₂₁ be the motor molecule concentrations (in units of tension) at the vertex in the junctional cortices of the three cells. Then, the tensions are related to the motor molecule concentrations as:

This uses the assumption of myosin continuity at vertices, and the fact that the tension on an interface is the sum of the tensions of the two cortices that make it up. The motor molecule concentration on the cortex belonging to the new interfaces, post collapse, will be equal to m₁₂. Solving for this in terms of the tensions:

The new interface consists of two cortices, coming from the two vertices of the old junction. Let the tensions at the two triangles be T₀, T₁, T₂ and . Let T_a be the active tension on the new junction immediately after the T1. It is the equal to

Note however that the total tension T_n on the new junction is not necessary equal to T_a. The total tension is defined geometrically from the angles at the new junction (or, equivalently, the tension triangle vertices). Indeed, generally, T_n > T_a, i.e. the active tension on the new junction is not enough to balance the tension due to the adjacent edges. We introduce a passive tension T_p on the new edge which balances this deficit

Figure S10:

For example, if a perfectly symmetric quartet collapses when the vertex angle facing the collapsing edge is 90°, and . Therefore, and . Note that by the triangle inequality, for any convex quadrilateral with perimeter P and diagonals D₁, D₂, one has P/2 ≤ D₁ + D₂ ≤ P . Applying this to the quadrilateral formed by the two tension triangles at the collapsing interface, we get T_a,n ≥ 0 and T_p,n ≥ 0: the handover formula always results in positive active and passive tensions. Further, the “handover” mechanism robustly generates irreversible T1s: if a junction were to collapse back after a T1, the newly formed junction would inherit high myosin levels and therefore be likely to collapse again.

The passive tension subsequently relaxes visco-elastically with rate . Combining this with the feedback equation yields Eq. (32). The relaxation rate τ_p = τ_T/6 was fit to the measured tension decay post intercalation.

Numerical implementation

We integrated Eq. (32) using a 4th order RungeKutta method as implemented in the SciPy software package [108], using as initial conditions the vertex angles in the experimental data at time t = 5 min (the vertex angles from the data were temporally smoothed with a window of 2 min to reduce noise).

In our simulations, we triggered an intercalation event once the collapsing edge length, as determined by energy minimization, reached 0. We then initialized the passive tension as given by Eq. (35) and simulated the combined active/passive tension dynamics. The overall timeand length-scale were likewise fit to the data. The feedback exponent was n = 4. Numerically, energy minimization was carried out using the scipy.optimize package, specifically its Nelder-Mead optimizer.

Passive intercalation

Within the same framework, we can also model passive intercalations, shown in Fig. 3F. We can either specify cortical tension dynamics (in this case, tension homeostasis, ) and then minimize an elastic energy that includes the influence of externally imposed shear. This will trivially reproduce the isogonal behavior shown in Fig. 2E. Alternatively one can reason that on the dorsal side there is very little active tension (cortical myosin), so that cortical tensions are expected to be low. Then, the passive elastic energy, together with the imposed shear, should determine the also angles. Both models are consistent with the data.

In Fig. 3F, we show simulations based on the latter option. We implement the externally imposed shear causing the passive intercalations as follows. We calculate the four cell centroids c_i of the cell quartets, and demand that at simulation time t

where S(t) is a diagonal matrix representing (area-preserving) shear, increasing over time with constant shear rate γ.

We then minimize the elastic energy under the constraint (numerically, this is done by adding a penalty term to the elastic energy), but this time with respect to both the lengths 𝓁_i and the angles ϕ_i. As long as the shear modulus µ ≪ λ, this procedure results in approximately isogonal behavior. We chose µ = 0.01λ. To obtain continuity across the intercalation, the rest shape S₀ needs to be adjusted. Otherwise, post intercalation, the cell quartet jumps back to a hexagonal lattice configuration. This can be done by introducing a viscous relaxation scale for the rest shape; for simplicity, we simply set the rest shape to the current shape tensor at the time of intercalation. The initial conditions are obtained by varying the angle between the shear axis and the collapsing interface orientation from 0° − 25°.

4 PIV and tissue tectonics analysis of experiments assessing the role of the posterior midgut

To obtain additional evidence on the influence of posterior midgut (PMG) invagination on GBE, we carried out a quantitative re-analysis of data published in Ref. [19]. In brief, in Ref. [19] – among other experiments and analyses – the authors used an IR-laser to fuse germ band tissue to the overlaying vitelline membrane by cauterization. This creates a “cauterization fence” along the the posterior end of the germ band which mechanically uncouples the germ band and the region of PMG invagination [19]. The authors then analyse tissue between the cauterization fence and the PMG. In wild-type embryos, they find that this tissue rips off from the cauterization fence, indicating pulling forces. This tissue ripping is also observed in mutants where force generation in the germ band is impaired (eve/rnt and twi/snl). In contrast, for tor mutants where the PMG does not invaginate, no tissue ripping is observed [19]. This suggests that the PMG exerts a pulling force on the germ band.

However, tissue ripping is only a qualitative measure for such pulling forces. To obtain a more quantitative estimate of the relative contributions of internally generated forces in the germ band and “external” pulling forces from the PMG we re-analyzed the movies from [19] using Particle Image Velocimetry (PIVlab [109]) to measure the average tissue flow in the region posterior to the cauterization fence. The results are shown in Fig. S11. In tor embryos, where no PMG invagination occurs, there is essentially no tissue flow, in line with the qualitative assessment from tissue ripping. Notably, in the eve/rnt and the twi/snl embryo, tissue flow velocity posterior of the cauterized region is reduced by about two-fold compared to wild-type. This shows that impairing active force generation in the germ band has a strong influence on the tissue flow anterior of the PMG. This flow originates from the dorso-ventral contraction of the germ band that pulls on the dorsal amnioserosa and drives the hyperbolic flow there [18]. These quantitative observations based on the data from [19] are in line with the much reduced tissue flow in eve and twist embryos imaged in-toto using light-sheet microscopy [74].

This above tissue flow quantification, is complemented by the following “anatomical” considerations. As shown in Ref. [64], the tissue that will invaginate to form the posterior midgut consists of two distinct regions, the “inititation” and the “propagation” region, both have an AP-length of approximately 7–8 cell diameters. On the other hand, the germ band has a length of approximately 40 cell diameters (Fig. 1A). Invagination of the initiation domain can hence create at most a strain ≈ 8/40 = 0.2 (since the initiation region is anchored on its anterior side). Moroever, the PMG forms a deep fold separating its anterior front, where the force driving invagination originates [64], from the germ band. It is unclear how a pulling force could be transmitted across this fold. Instead, the pulling force originating from the hyperbolic fixed point on the dorsal pole bypasses the PMG laterally where it manifests as isogonal stretching along the AP axis (see Fig. 6B’).

Figure S11: