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 “convergent extension”—a common motif of early development and organogenesis in many organisms—where tissue elongates along one axis while contracting along the perpendicular direction [3]. Convergent extension also exemplifies the important role of cell and tissue mechanics in morphogenetic processes, which must be studied alongside developmental genetics and signaling. Mechanics is what directly determines physical shape, and also plays a key regulatory role further upstream by facilitating self-organized mechanical pattern formation [4, 5], in diverse systems such as feather primordia [6], the amniote primitive streak [7], and hydra body axis formation [8]. Convergent extension is often conceptualized on the continuum scale, where the morphogenetic process resembles a shear flow [9, 10]. The key challenge is to relate the large scale flow of tissue to its behavior on the scale of cells, where the active, local forces that drive the flow are generated and controlled.

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 [11]. The capacity of epithelial tissue to support tension makes it markedly different from passive fluids [12]. How does one 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? By analysing empirical data on the dynamics of cells during tissue remodelling, we shall arrive at a model of morphogenetic flow as adiabatic deformation of cell array geometry controlled by the feedback-driven changes in the internal state of tension.

One of the best-established models of morphogenetic mechanics is provided by the Drosophila embryo, which is ideal for live imaging and offers an extensive set of genetic tools [13]. Progress in live imaging and computational image analysis has produced remarkably quantitative data [14, 15]. In this work, we will take advantage of previously published data [15] 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 [15] which, thanks to computational image processing (surface extraction [16], cell segmentation and tracking [15]), 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 [17]. 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 [18–20] (see Movie 1). GBE has been extensively studied, leading to the identification of relevant developmental patterning genes [18–21]. Live imaging has also uncovered the pertinent cell behavior, namely intercalation of neighboring cells in the lateral ectoderm [21–23]. The elementary process of cell rearrangement 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” [24] 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 so amnioserosa) deforms in the opposite direction 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 is consistent with the high degree of spatial coordination of rearrangements that is apparent in the lateral ectoderm (Fig. 1C, bottom).

Figure 1:

Intercalations are associated with localized activity of the force generating protein non-muscle myosin II [20, 22, 23] (henceforth, “myosin”). Myosin recruitment depends on the expression of key developmental genes [19, 20, 22, 25], and, in addition, is subject to positive and negative mechanical feedback that depends on the stress [26, 27] and the rate of strain (rate of cell deformation) [28]. Despite considerable understanding of the genetic and cell-biological components involved in GBE, the relative roles of the genetic pre-pattern versus mechanical feedback-driven “self-organization” of myosin activity are still not clear. Similarly, despite rich experimental data and great biological interest, no consensus has emerged regarding the cell-level mechanics of GBE [24, 26, 29–31]. Key questions remain open: How is myosin dynamics controlled on a cellular level? How are intercalations are driven mechanically [29, 32]? What is the relative strength of locally generated forces vs pulling by adjacent tissues? How is active force generation orchestrated across cells to create a coherent global morphogenetic flow? What sets the total amount of tissue elongation? These questions call for a coherent theoretical framework that allows interpreting and reconciling existing experimental findings.

The rich dataset shown in Fig. 1 exposes a wealth of information about the mechanical behavior of cells and thus allows us to address these problems. Our approach is based on the assumption of force balance of stresses concentrated in the cell cortices. Cortical force balance provides a direct link between mechanics and geometry: It allows inference of tensions from the angles at the vertices vertices where interfaces meet [33–35], 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). We will use these insights to formulate a model for an actively intercalating cell quartet where cortical tensions dynamically self-organize via a positive feedback mechanism while remaining in force balance. 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 [36]. 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 [36], we employ this order parameter to quantitatively compare tissue-scale simulations with experimental data 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 [37, 38], 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 [35]. 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) [39]. 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”[40]. 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 [41] which exerts negligible drag forces on tissue motion on the surface [42], 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 [33]. This link of hard-to-observe cortical tensions to the readily observable local cell geometry enables tension inference methods which have been extensively validated by computational robustness checks [34], direct comparison with measured laser ablation recoils [43] and by correlation with local myosin abundance [35].

Crucially, the link between force balance and the geometry of the cell tessellation 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 [12], as illustrated in Fig. 3B. 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 tessellation: For each cell, there is a corresponding vertex in the tension triangulation (Fig. 3B’). 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 [44] and have been applied to the statics of beam assemblies [45] and granular materials [46].) This triangulation establishes a geometric structure in tension space. Force balance requires that the angles at vertices in the physical cell tessellation 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 tessellation 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.

Figure 3:

Cell scale analysis

By applying local tension inference to each interface – using the angles it forms with its four neighboring interfaces – we can map out the relative tensions in the tissue, as shown in Fig. 2C and C’. Initially, relative tensions are close to unity throughout the embryo, since the cell tessellation 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 (amnioserosa) and the lateral ectoderm (germ band). 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 [10]). 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 [47]. 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 [28].

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 [26, 27] and is in excellent agreement with predictions from a recent model where such feedback is a key ingredient [48]. 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 jumps to a lower value and then continues decreasing back to 1, corresponding to vertex angles of 120^°. The jump in relative tension between the collapsing and the emerging interface 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 a jump of the relative tension from to .

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 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 central interface.

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

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 cortical tension behaviour is quite counter-intuitive and differs significantly between distinct spatial regions of the embryo. In particular, it is very differently 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 are can change independently and are actively regulated by turnover of junctional proteins (such as actin, myosin and E-cadherin).

This raises the question how the regulation of tensions controls the cell and tissue shape. The answer to this question lies in the force balance constraint that links tension and geometry which we have previously used to inferring the tensions from the observed epithelial geometry. Recall that force balance implies that the cortical tensions form a triangulation such that the angles at vertices in the physical cell tessellation are complementary to those in the tension triangulation (Fig. 3A, left). However, force balance does not uniquely determine the cell tessellation (i.e. the positions of the cell vertices) since interface lengths may collectively change while preserving angles via so-called isogonal modes [12], as illustrated in Fig. 3A, right. Thanks to the relative sliding of cytoskeletal elements (e.g. due to the turnover of passive crosslinkers or walking of myosin motors), cell-cell interfaces can extend or contract under constant tension, behaving similar to a muscle rather than a spring. These considerations motivate a decomposition of the tissue deformations into two distinct contributions: deformations of the tension triangulation reflecting changing cortical tensions, and isogonal (angle-preserving) deformations (see Fig. 3B). We will call this a tension–isogonal decomposition. The isogonal degrees of freedom allow for tissue deformations at fixed tensions (tension homeostasis) and therefore account for cell shape deformations driven by external forces or cell-internal stresses (hydrostatic pressure, medial myosin [29, 32], stiffness of the nucleus, intermediate filaments and microtubules [49–51]). It is important to keep in mind that this association of isogonal deformations with non-junctional (e.g. cell-elastic or external) stresses is based on the assumption that active cortical stresses are the dominant source of stress in the tissue (separation of scales).

Implicitly, the tension–isogonal decomposition uses the fact that we can construct a reference state for the cell tessellation—compatible with the force balance constraints—from the tension triangulation (e.g. using a Voronoi tessellation). Isogonal deformations then transform this reference state into the physical cell tessellation. In practice, we do not need to explicitly construct the reference cell tessellation. Instead, we utilize the fact that the physical cell tessellation can be quantified by the triangulation defined by the cells’ centroids [52]. 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. 3C). 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. 3C–C’. The quartet’s isogonal deformation tensor I_quartet transforms the cell centroids into the tension vertices . Given and c_i, I_quartet can be obtained using least squares. The constant, global scale factor η≈ 6 μm translates units from relative tensions (a.u.) to length (μm) and is fitted to minimize initial strain in the initial state. In Fig. 3C and C’, we plot the quartet shape aspect ratios and the elongation factor of the isogonal deformation against the length of the quartet’s inner edge, which serves as a pseudo-time parameterization of the T1 process. (Plots against real time are shown in Fig. S6B and B’.) In the lateral ectoderm, the tension-triangle deformation leads the quartet shape deformation, indicating that tension dynamics drives the intercalation (Fig. 3C). 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. 3C’). 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 [53] or due to intracellular actin cables [54], 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 [12].

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. 3C), 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. 3C’).

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. From the decomposition Eq. (1) it is 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. 3E). By a geometric construction (Fig. S6D) one finds that the corresponding isogonal deformation has an aspect ratio of 3, in good agreement with observations (Fig. 3C’).

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 (see Fig. 3D-D’). For simplicity, we consider the case of two identical, isosceles 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 [36], 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 (see Fig. 3D’). In general, the condition for the neighbor exchange defines the “T1 threshold” in the Δℓ_iso-T plane (see Fig. 3E). 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 lateral ectoderm, in particular near the VF (purple region in Fig. 1A), DV oriented interfaces are stretched due to VF invagination. The average isogonal contribution to interface length ⟨Δℓ_iso⟩≈ 1.9 μm of DV oriented interfaces in the lateral ectoderm can be estimated from the tissue-scale isogonal strain there (see Tissue Scale section below). One can read off the predicted relative tension threshold T_crit ≈ 1.53 in Fig. 3. This is in excellent agreement with the value found by tension inference (Fig. 2D). For twist and snail mutants where the VF is abolished [55] such that there is no isogonal stretching of the lateral ectoderm, we predict a T1 threshold for active T1s. This implies that there will be no tension jump between the contracting and the extending interface for active T1s.

Minimal local model based on tension feedback reproduces length- and tension-dynamics of T1s

We conclude our analysis of intercalating cell quartets by condensing the insights gained into a parsimonious model that combines tension dynamics with a constitutive relation for the isogonal modes. Importantly, we model the dynamics in adiabatic force balance, i.e. on a timescale where the relaxation of elastic energy has already come to its conclusion. That means that the positions of the vertices do not follow an explicit dynamical equation (like the overdamped energy relaxation dynamics in a classical vertex model) but are instead implicitly determined by the tension triangulation and the isogonal modes. The dynamics resides in the changes in tension that transform and remodel the tension triangulation.

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. 4A). The vertex angles are determined by the relative cortical tensions T_i via the force balance condition (see Fig. 3A). Motivated by the nonlinearly increasing tension observed on contracting interfaces (cf. Fig. 2D), we equip each interface with self-reinforcing tension dynamics

which models tension-induced myosin recruitment [26] (with excitability exponent n set to n = 4 in our simulation). The dynamics conserves the total active tension , 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. [36].) 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. 4B) 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.

Figure 4:

To find the physical cell shape for given cortical tensions, we need to determine the interface lengths ℓ_i, i.e. fix the isogonal mode (Fig. 4C). To this end, we introduce a cell elastic energy which is the isogonal mode must minimize. This elastic energy accounts for internal structures of the cell such as the nucleus and microtubuli. We assume the scale of this 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. We write the 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. 4A) 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. 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 [56] (see also Fig. S9D). 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. 4E 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 [23, 57]. 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. 4D.

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 [58] (see Fig. 4D). (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 [33].) 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 [23, 57].

Simulating this model, we find good qualitative agreement between the single-quartet model and experimental data, both for the tension dynamics (compare Fig. 4E,F and Fig. 2D,E) and for the tension–isogonal decomposition (compare Fig. S9A,A’ and Fig. 3C,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. 4F; 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. In passing, we note that we have ignored myosin fluctuations in the above arguments. If these fluctuations are strong compared to the average myosin levels, they may cause the collapse of the emerging interface, thereby reversing a T1, as observed in some genetic mutants [47].

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 [10, 57, 59, 60].

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 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 [28]. 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 [36] 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 [24] 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. 3C). 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 away from the VF also accumulates some isogonal extension along the DV axis (green shaded region, see time traces in Fig. 6E). This causes an increase of the tension threshold for active T1s (cf. Fig. 3E).

Figure 6:

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 [15] and low myosin density in the amnioserosa [10]. Taken together, isogonal strain identifies regions where passive tissue deformation is caused by mechanical coupling across the tissue.

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. 7A (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 [26, 61]. 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. 7B): 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. 7A, right), corresponding to an obtuse triangle (Fig. 7B, top right). In contrast, an acute tension triangle (Fig. 7F, 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. 7A, left). We can therefore use the tension triangle shape to define a local tension configuration (LTC) parameter. The LTC-parameter space (i.e. the space of triangle shapes), shown in Fig. 7B, is spanned by two axes, measuring how elongated and how acute or obtuse the tension triangle is (see Ref. [36], Fig. S7 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 7:

Figure 8:

LTC order choreographs T1s

The geometric condition for neighbor exchanges discussed above (cf. Fig. 3E) defines a threshold in the LTC-parameter space (dashed line in Fig. 7B, see [36] 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. 7C). The bridge fraction reaches its maximum around 25 min, just at the onset of GBE (see Fig. 7D). The corresponding alternating pattern of tensions is clearly visible in Fig. 7A (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. 7C). 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. 7G, 50 min; and Fig. 7H). In fact, the distribution of tension triangle shapes appears to approach that of a random Delaunay triangulation (see SI Fig. S8), 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.

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 [12] 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. Yet, the focus on the geometry of the force balance 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 [34, 35, 42, 43].

In the early Drosophila embryo, we observe both active and passive T1s. 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. We predict that the local relative tension threshold will be lower (namely ) in mutants which lack a ventral furrow (twist and snail [55]).

Cells orchestrate tissue flow by self-organizing in tension space

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 [62–64]. Mechanical stresses and strains can propagate over long distances and contain information about the tissue geometry (e.g. in the form of hoop stresses [65]). In an epithelium dominated by cortical tensions, this mechanical environment is 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 tessellation 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 area 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 [66–68]. We hypothesize that tension homeostasis allows the amnioserosa to undergo large cell deformations while maintaining tissue integrity [67, 69]. 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 [26]. The underlying molecular mechanisms are yet to be identified and might rely on the catch-bond behavior of myosin [70, 71] or mechano-sensitive binding of α-catenin [72, 73].

We have found that T1 events can be explained quantitatively a simple geometric criterion which defines a “T1 threshold” in terms of on 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 [36], 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 [29, 74].

Our findings on the dynamics of interface length and tension during active T1s (Fig. 2D) are similar to another recently proposed model [48]. 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. [48] and might be obtained in the limit of fast mechanical relaxation from the latter.

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 [36] 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 self-limiting 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 [28] has shown that mechanical feedback loops in the embryo are patterned along the DV axis. Notably, the DV patterning system is conserved across species [75]. Our work also highlights a second interaction between tissue flow and DV patterning: because the coherent T1s we observe do not scramble cells in the tissue (Fig. 1A), distinct DV fates (e.g. neural and surface ectoderm) do not mix, 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 [28] and hoop stresses resulting from ellipsoidal embryo geometry [65], since we observe anisotropy already prior to VF formation), or indeed due to AP-striped genetic patterning [76, 77]. 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 [28].

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 [36], we show that mechanical self-organization via tension feedback also allows cells to coordinate their behaviors (such as active T1s) with their neighbors. 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. 7A, center)—from predetermined genetic patterns alone. Instead, the emerging overarching picture is that of morphogenesis driven by an interplay of bottom-up self-organization and top-down genetic patterning [5, 78].

Going forward, we expect that geometric insight into tension space will provide a new perspective on many different systems [79], for example Xenopus neural tube formation [80], amniote primitive streak formation [9, 81, 82], wing disk elongation [83, 84], the sea urchin archenteron [85], or kidney tubule elongation [86]. Further, it will be interesting to study intercalations driven by actin “protrusions” [87]: they too are governed by cortical force balance on the cellular level, making tension inference and the tension–isogonal decomposition applicable.

Acknowledgements

We thank Matthew Lefebvre, Noah Mitchell, Sebastian Streichan, and Arthur Hernandez for stimulating discussions and careful reading of the manuscript. FB acknowledges support of the GBMF post-doctoral fellowship (GBMF award #2919). NHC was supported by NIGMS R35-GM138203 and NSF PHY:1707973. 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 was obtained from the repository DOI:10.6084/m9.figshare.18551420.v3 deposited with Ref. [15]. 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.

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. [24]; see also [88] and [52]. 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. 3D and Fig. 7B). (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. 3E). 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. 3C and C’ and Fig. S6), 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. S6D).

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.

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

Figure S1:

Figure S2:

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.

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). These 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 [12]. 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 with the additional constraint that the average of the T₁₂, T₂₃, T₃₄ and T₁₄ is 1. The solution 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 in the following.

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). 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 interfacial angles in the cell tessellation, one can find a unique tension triangulation (up to a global scale factor) for a given physical cell tessellation in force balance (see dashed gray arrow in Fig. S6A). However, the reverse is not true. For a given tension triangulation, there is a continuous family of cell tessellations 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 S3:

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 (see Sec. ?? and Sec. 3 below).

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 tessellation. 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 [52]. (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 interface unit vectors e_ij (see Sec. 1.4) and a subsequent rotation of the tension vectors by π/2. To fix the tension scale, we normalize the area of the tension triangle to one. 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.¹

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 η I_vertex. . The scale factor η converts units of tension (a.u.) into units of length. How we fix η is described below.

Because the isogonal modes must preserve angles, they cannot include rotations. In other words, I_vertex must be a symmetric tensor, and therefore has only 3 degrees of freedom. However, the above equations impose 4 constraints. We therefore determine I_vertex using least squares, i.e. by solving the minimization problem

where Sym₂ is the set of symmetric 2 × 2 matrices.

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. 7L). For an equilateral triangle, the edge length is determined by the area, which we have normalized to one, giving an edge length 2/3^1/4≈ 1.52. The average edge length of the centroidal triangles is given by 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. This scale factor minimizes the isogonal strain at t = 0 (see Fig. S5B). Indeed, minimizing ⟨||I||⟩ at a reference time is provides a more general way to determine η in cases where the initial triangulation is not as uniform and regular as in the Drosolphila 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.

Figure S5:

Figures 3D 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. 3E), 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)22⟩≈ 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.

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

Quartet shape tensors

The tension–isogonal decomposition for cell quartets is calculated analogously to the procedure for single vertices. The two tension triangles associated to the inner vertices of the cell quartet form a “kite” (see Fig. S6A). 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 tessellation

For the tension–isogonal decomposition, we have circumvented the explicit construction of a reference cell tessellation 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.

Figure S6:

The explicit construction of the reference cell tessellation 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 [36] 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 [89].

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.

Figure S7:

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. S7A).

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. S7B). 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. S7B) 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 performa 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. (17).

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. S7B), we define

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

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

Figure S8:

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 [90]. A small sample of such a random triangulation is shown in Fig. S8A. The histograms in the triangle shape space Fig. S8B 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. S8D).

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. [36]) 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. [12]. 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. and .

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 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 .a 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 (23) and (24) 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 (27) and (28) 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. 4, 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 [91]). 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 mechano-sensitive 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. (35) w.r.t. the ℓ_i, for the angles determined by the tension dynamics.

Comparison to area-perimeter elastic energy

Strikingly, within the single-quartet model, the “shape strain” S_𝒞 −S₀ can always be set to 0 by choice of the edge lengths, and the energy E_𝒞 = 0 throughout. This can already be seen from a degree-of-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 area-perimeter elastic energy E = (A−A₀)² +(P −P₀)² (where A, P are the cell area and perimeter, and A₀, P₀ their target values) [56], there exists an energy barrier, and the inner interface ℓ₀ only collapses when ϕ₀ = π. This is shown in Fig. S9B. 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_𝒞 − S₀)²], (or perimeter elasticity with cell shape bulk elasticity) leads to similar results as Eq. (35). 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 convergence-extension takes place [91].

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. 4D, 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. [58] 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 4D).

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

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. (33). The relaxation rate τ_p = τ_T/6 was fit to the measured tension decay post intercalation.

Numerical implementation

We integrated Eq. (33) using a 4th order Runge-Kutta method as implemented in the SciPy software package [92], 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. (36) and simulated the combined active/passive tension dynamics. The overall time- and 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. 4F. 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. 4F, 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^°.

Figure S9:

4 Movie captions

Movie 1: Side view and flat projection of one side of the embryo surface. Ventral furrow invagination (gray cells on the dorsal side) is followed by convergence–extension of the germ band (lateral ectoderm cells colored purple, blue and green). As the germ band elongates along the AP axis, the cells move over the posterior pole. The amnioserosa (orange and red cells) undergoes convergence– extension in the opposite direction of the germ band and exhibits significant cell shape elongation while cells in the germ band remain mostly isotropic in shape. Near the end of germ band extension (ca 35 min) cell divisions start. (Corresponds to Fig. 1A; invaginating cells are colored in gray; cells in the head are colored white; cells after division events are colored in light gray.)

Movie 2: Relative tension dynamics in the lateral ectoderm. Relative junctional tensions inferred from cell geometry reveals show the emergence of an alternating pattern of high and low tensions that organizes cell intercalations (T1 transitions). Coherent intercalations drive convergent–extension tissue flow which slows down significantly as cell scale order is lost.

Movie 3: Simulation of single intercalating cell quartet. Simulation of an intercalating cell quartet driven by positive tension feedback and myosin handover mechanism, corresponding to Fig. 4E. Out of the ensemble from Fig. 4E, the movie and shows a simulation for symmetric initial tensions (i.e. equal initial tensions on the two non-collapsing interfaces T₁ = T₂). After the edge flip, the blue parts of the inner tension triangulation edge indicate the passive tension on the newly formed interface. The passive tension rapidly relaxes while the active tension grows due to positive tension feedback.