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 (${\displaystyle {\displaystyle \sim 10\mathrm{s}}}$) on the scale of cells is at least 10-fold faster than the timescale of local tissue flow (Bambardekar et al., 2015; Munjal et al., 2015) 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 remodeling of a quasistatic force balance network (Noll et al., 2020). 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 (Figure 2A, Lecuit and Yap, 2015). The myosin motors exert a contractile force on the actin fibers and thereby keeping the cortex under active tension (${\displaystyle T}$), which we refer to as ‘cortical tension’ (Prost et al., 2015). Together, the adherens junctions define a tissue-wide mechanical structure, capable of (i) generating locally controlled internal tension and (ii) adaptively remodeling its architecture. The Drosophila blastoderm lies on top of a fluid yolk (Doubrovinski et al., 2017) which exerts negligible drag forces on tissue motion on the surface (Cheikh et al., 2022), suggesting that all forces are balanced within the epithelial layer.

Figure 2 with 2 supplements see all

Download asset Open asset

The concentration of tension in interfaces tightly links force balance with the readily observable geometry of cell-cell interfaces. Force balance requires that the forces ${\displaystyle {\mathbf{T}}_{ab}}$ exerted by the cortical tensions sum to zero at each vertex where three (or more) interfaces meet. (Here, the indices ${\displaystyle a,b}$ label the two cells that meet at an interface; see Figure 2A). These force-balance constraints relate to the relative tension on cell interfaces at a vertex to the angles at which the interfaces meet (Chiou et al., 2012). This relation between the hard-to-observe cortical tensions and the readily observable local cell geometry enables tension inference methods which have been extensively validated by computational robustness checks (Ishihara et al., 2013), direct comparison with measured laser ablation recoils (Kong et al., 2019) and by correlation with local myosin abundance (Noll et al., 2020).

Crucially, the link between force balance and the geometry of the cell array goes beyond inference and to the very basis of our proposed mechanism of tissue flow. Force balance implies that the tension vectors meeting at a vertex sum to zero and hence forming a a closed triangle (Noll et al., 2017), as illustrated in Figure 2B. 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 ${\displaystyle {90}^{\circ }}$. Two angles are complimentary if they sum to 180°). As adjacent vertices share an interface, the corresponding tension triangles must share an edge: the tension triangles form a triangulation, i.e., fit together to form a tiling of the plane. The triangulation is dual to the cell array: For each cell, there is a corresponding vertex in the tension triangulation (Figure 2B’). 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 Clerk Maxwell, 1864 and have been applied to the statics of beam assemblies Varignon, 1725 and granular materials Tighe et al., 2008). This triangulation establishes a geometric structure in tension space. Force balance requires that the angles at vertices in the physical cell array are complementary to those in the tension triangulation. This intimately links tension space and epithelial geometry in real space. Myosin-driven local changes in tension changes the shapes of tension-space triangles and hence remodel the tension triangulation. The induced changes in the geometry of the cell array drive both the local rearrangement of cells and the global tissue flow. The remainder of the Results section will provide quantitative data analysis and modeling that will work out of the above ideas and their implications.

To map out the relative tensions in the tissue, we perform local tension inference for each interface relative to its four neighbors. Geometrically in tension space, the local tension configuration is represented by two adjoined tension triangles forming a kite (see Figure 2B) and the inference is a simple application of the law of sines on the two triangles (see Appendix 2). Figure 2C and C’ show snapshots from this local tension inference in the dorsal ectoderm (amnioserosa) and the lateral ectoderm (germ band), respectively. Initially, relative tensions are close to unity throughout the embryo, since the cell array is approximately a hexagonal lattice, with vertex angles close to ${\displaystyle {120}^{\circ }}$. As gastrulation progresses, the cortical tensions change, and one starts to see characteristic differences between the dorsal ectoderm and the lateral ectoderm. In the former, interfaces that contract remain under approximately constant tension while interfaces oriented parallel to the direction of tissue stretching (i.e. along the DV axis) extend and are under increasing tension (Figure 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, Figure 2C’ and Video 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.

Video 2

Download asset

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 Streichan et al., 2018). 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.

We identify all cell quartets that 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 (Stern et al., 2020). Pooling the data for each of the bands of cells colored in Figure 1A, we find two distinct scenarios for ventrolateral quartets and dorsal quartets (Figure 2D and E; for a breakdown by individual bands, see Figure 2—figure supplement 1).

The length dynamics of collapsing and extending interfaces in the amnioserosa and the lateral ectoderm is qualitatively similar (Figure 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 (Gustafson et al., 2022).

By contrast, the tension dynamics are 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 (Figure 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 (Blankenship et al., 2006; Fernandez-Gonzalez et al., 2009; Duda et al., 2019) and is in excellent agreement with predictions from a recent model where such feedback is a key ingredient (Sknepnek et al., 2023). From the data shown in Figure 2D we can read off the average relative tension threshold ${\displaystyle ⟨T_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}⟩=1.530\pm 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 (Figure 2H). Conversely, relative tensions below 1 correlate with interfaces that never collapse (Figure 2J).

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

Let us now turn to T1s in the amnioserosa. Here, the relative tension in the inner edge remains almost constant near 1 prior to the neighbor exchange, i.e., the vertex angles remain close to ${\displaystyle {120}^{\circ }}$ (see Figure 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 Figure 2—figure supplement 2). On the new interface emerging after the neighbor exchange, tension is high and remains constant for an extended period. Again, the apparent tension jump across the neighbor exchange is a consequence of geometry. Just before the neighbor exchange, the vertex angles are close to ${\displaystyle {120}^{\circ }}$ which implies that the angles facing away from the interface are ${\displaystyle {60}^{\circ }}$ after the neighbor exchange (see Figure 2G). This corresponds to a relative tension of ${\displaystyle \sqrt{3}\approx 1.73}$ on the central interface.

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

Tension inference and the pooled analysis of T1 events have revealed the cortical tension dynamics on cell-cell interfaces during active and passive cell intercalations. The behavior of cortical tensions differs significantly between distinct spatial regions of the embryo and, at first glance, appears quite counter-intuitive. In particular, it is very different from that of springs or rubber bands. The length and tension of a springs are tied to one another by a constitutive relationship. By contrast, the length and tension of cell-cell interfaces can change independently and are actively regulated by the turnover of junctional proteins (such as actin, myosin, and E-cadherin). Experiments where actin turnover is decreased by Cytochalasin D treatment show reduced cell intercalations supporting the point of view that the decoupling of tension regulation from a spring-like constitutive relation is important for physiological behavior.

These observations call for a new modeling approach, where tensions are not governed by constitutive relationships such as the typical area–perimeter elasticity of the vertex model (Farhadifar et al., 2007; Hufnagel et al., 2007). Instead, our model directly builds on the same assumptions that underlie the tension inference: dominant cortical tensions that are in adiabatic force balance. The dynamics reside in the changes in tension governed by mechanical feedback.

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

which models tension-induced myosin recruitment (Blankenship et al., 2006) (with exponent ${\displaystyle n=4}$ in our simulation). The dynamics conserve the total active tension ${\displaystyle \sum _i{T}_i}$, corresponding to a finite myosin pool in each cell. (Other choices for fixing the overall scale of tension, e.g., via the triangle area, as well as other values of ${\displaystyle n}$ do not change the results qualitatively; see Claussen et al., 2024). The timescale of tension dynamics ${\displaystyle {\tau }_{\mathrm{T}}}$ is fitted to the experimental data. Positive feedback drives the relative tension on the quartet’s central interface towards the tension threshold discussed below (see Figure 4E) and thus drives active T1s. The interface that collapses is the one that starts with the highest initial tension, which we label with ${\displaystyle i=0}$ by convention, such that ${\displaystyle T_0>T_{1,2}}$.

Figure 3 with 1 supplement see all

Download asset Open asset

Figure 4 with 1 supplement see all

Download asset Open asset

The tension dynamics described above determine the angles at vertices via adiabatic force balance but leave the interface lengths ${\displaystyle \ell }$ as independent degrees of freedom. Further below, we will identify these degrees of freedom as the so-called isogonal (angle preserving) modes, which are an immediate consequence of having cortical tensions governed by feedback instead of a constitutive relation (Noll et al., 2017). To find the lengths ${\displaystyle {\ell }_i}$, we need to account for subdominant mechanical contributions from internal structures of the cell, such as medial myosin contractility (Collinet et al., 2015; Vanderleest et al., 2018) the nucleus (Kim et al., 2024) and microtubules (Ramms et al., 2013; Charrier et al., 2018; Pensalfini et al., 2023; Singh et al., 2024). We account for these contributions through a ‘cell elastic energy’ in terms of the deviation of the cell ‘shape tensor’ ${\displaystyle {\mathbf{S}}_{\mathrm{C}}=\sum _i{\ell }_i{\mathbf{e}}_i\otimes {\mathbf{e}}_i}$ from a target shape tensor ${\displaystyle {\mathbf{S}}_0}$.

Here, ${\displaystyle {\mathbf{e}}_i}$ are the unit vectors pointing along the interfaces with lengths ${\displaystyle \left|{\mathbf{e}}_i\right|={\ell }_i}$ (see Figure 3A) and ⊗ denotes the tensor product. We choose this shape tensor because it is linear in interface length, meaning that it does not change if an interface is subdivided by inserting an additional vertex.

Importantly, we assume the scale of ${\displaystyle E_{\mathcal{C}}}$ to be much smaller than the elastic energy due to cortical tensions. Because of this separation of scales, the angles at vertices are fixed by tensions, and the relaxation of the cell elastic energy only affects interface lengths. The coefficients ${\displaystyle \lambda }$ and µ control the cell’s resistance to isotropic compression/dilation and shear deformations, respectively (see Appendix 5 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 (Bi et al., 2015) (see also Figure 3—figure supplement 1D). 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 ${\displaystyle <0}$ in Figure 3E, F) and qualitatively reproduces the dynamics of interface length and tension observed in the experimental data (Figure 2D, E). 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 (Bertet et al., 2004; Rauzi et al., 2008). 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 ${\displaystyle T}$ into the myosin-borne, active tension ${\displaystyle T_{\mathrm{m}\mathrm{y}\mathrm{o}}}$ and the passive tension ${\displaystyle T_{\mathrm{p}}}$. Cortical remodeling due to turnover of crosslinkers will gradually relax this passive tension. We account for this by an exponential decay ${\displaystyle T_{\mathrm{p}}\propto e^{-t/{\tau }_{\mathrm{p}}}}$ with a characteristic remodeling rate ${\displaystyle {\tau }_{\mathrm{p}}}$. This passive tension effectively represents Maxwell-type viscoelasticity of the extending interface, illustrated by the spring and dashpot in Figure 3D. Such visco-elastic behavior on a scale of minutes has been reported in experimental data from the Drosophila embryo (Clément et al., 2017).

Completing the model requires understanding 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 (Kale et al., 2018; see Figure 3D). (Force balance requires that the total tension, the sum of the tensions in the two abutting cortices, is uniform along an interface. However, the individual tensions on either side can be non-uniform, as the resulting traction forces are exchanged via E-cadherin molecules linking the cells along the interface Chiou et al., 2012). 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 (Bertet et al., 2004; Rauzi et al., 2008).

Simulating this model, we find good qualitative agreement between the single-quartet model and experimental data, both for the tension dynamics (compare Figure 3E and F and Figure 2D and E) and for the tension–isogonal decomposition which we will introduce below (compare Figure 3—figure supplement 1A, A’ and Figure 4C and C’). Video 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 (Figure 3F; see Appendix 5 for details).

Video 3

Download asset

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

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

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

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

The key elements of the tension–isogonal decomposition are illustrated in Figure 4B. Changes in cortical tension are reflected in deformations of the tension triangulation (i.e. changes in the angles at vertices). For each tension triangulation, we can construct a reference state for the polygonal cell array – compatible with the force balance constraints – from the tension triangulation using a Voronoi tessellation. Isogonal deformations then transform this reference state into the physical cell array. In practice, we do not need to explicitly construct the reference cell array. Instead, we utilize the fact that the physical cell array can be quantified by the triangulation defined by the cells’ centroids (Merkel et al., 2017). 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 labeled ${\displaystyle i=1-4}$, we quantify the physical shape in terms of the centroid positions ${\displaystyle {\mathbf{c}}_i}$ by calculating the moment of inertia tensor. The tension reference shape is given by the tension vectors ${\displaystyle {\hat{\mathbf{T}}}_{ij}}$, which are obtained from the interface angles using tension inference and are rotated by ${\displaystyle {90}^{\circ }}$ relative to the interfaces (see Figure 2B). The ratios of eigenvalues of the moment-of-inertia tensors define the shape aspect ratios, illustrated by ellipses in Figure 4C and C’. The quartet’s isogonal deformation tensor ${\displaystyle {\mathbf{I}}_{\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{t}}}$ transforms the cell centroids into the tension vertices ${\displaystyle {\ell }_0{\mathbf{I}}_{\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{t}}.{\hat{\mathbf{T}}}_{ij}={\mathbf{c}}_i-{\mathbf{c}}_j}$ (see Figure 4—figure supplement 1A for a labeled sketch). Given ${\displaystyle {\hat{\mathbf{T}}}_{ij}}$ and ${\displaystyle {\mathbf{c}}_i}$, the equations can be (approximately) solved for the best fitting ${\displaystyle {\mathbf{I}}_{\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{t}}}$ using least squares (an exact solution is not possible in general because the equations are overdetermined; see Appendix 3). The constant global scale factor ${\displaystyle {\ell }_0\approx 4.2\mu \mathrm{m}}$ translates units from relative tensions (a.u.) to length (${\displaystyle \mu \mathrm{m}}$) and can be calculated from the average interface length ${\displaystyle ⟨\ell (0)⟩\approx 3.5\mu \mathrm{m}}$ at the first timepoint using the Voronoi–Delaunay duality as explained in Appendix 3.

In Figure 4C and C’, we plot the quartet shape aspect ratios and the eigenvalue ratio of the isogonal deformation tensor against the length of the quartet’s inner edge. The latter serves as a pseudo-time parameterization of the T1 process; plots against real-time are shown in Figure 4—figure supplement 1B, B’. In the lateral ectoderm, tension-triangle deformation accounts for quartet shape deformation, indicating that tension dynamics drives the intercalation (Figure 4C). This confirms our conclusion from the relative tension analysis (Figure 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 (Figure 4C’). Here, the tension mode is constant—indicating tension homeostasis—while the cell quartet deformation is entirely stored in the isogonal mode. While the behaviors observed in the lateral ectoderm and the amnioserosa represent the extreme cases, being purely tension-driven (active) in the former and purely isogonal (passive) in the latter, there is a 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 (Bocanegra-Moreno et al., 2023) or due to intracellular actin cables (Priess and Hirsh, 1986), the isogonal mode is actively controlled. Another example of this is apical constriction during VF invagination which is driven by medial myosin pulses that drive isogonal cell contractions (Noll et al., 2017).

A puzzle that has remained open so far is the threshold value of relative tension where a neighbor exchange takes place. In the lateral ectoderm where the isogonal strain is approximately constant (Figure 4C) the neighbor exchange takes place at a critical relative tension of ${\displaystyle 1.530\pm 0.005}$ on average (Figure 2D). By contrast, in the amnioserosa, where T1s take place under constant tension (Figure 2E), there is an average critical aspect ratio ${\displaystyle \sim 3}$ of the isogonal deformation at the point of neighbor exchange (see Figure 4C’).

We now explain these numbers in terms of the geometric tension–isogonal decomposition introduced in the previous section. We start by decomposing the length ${\displaystyle \ell }$ of a cell quartet’s central interface into a reference length ${\displaystyle {\ell }_{\mathrm{r}\mathrm{e}\mathrm{f}}(\{T_i\})}$, determined purely in terms of the local relative tensions, and an isogonal length ${\displaystyle \mathrm{\Delta }{\ell }_{\mathrm{i}\mathrm{s}\mathrm{o}}}$ depending on the quartet’s isogonal strain:

A neighbor exchange takes place when ${\displaystyle \ell }$ shrinks to zero. Thus, setting ${\displaystyle \ell =0}$ in Equation 3 defines the threshold values of the local tensions ${\displaystyle \{T_i\}}$ and the isogonal length ${\displaystyle \mathrm{\Delta }{\ell }_{\mathrm{i}\mathrm{s}\mathrm{o}}}$ where T1s take place (see green line in Figure 4E). It is immediately apparent that the interface collapse can be driven by the isogonal mode, by changes in tensions, or by a mixture of both.

Let us start with the passive T1s observed near the dorsal pole. In this case, the relative tensions are homeostatic (Figure 2E), so ${\displaystyle {\ell }_{\mathrm{r}\mathrm{e}\mathrm{f}}\approx ⟨\ell (0)⟩\approx 3.5\mu \mathrm{m}}$ remains constant. The neighbor exchange, therefore, happens when ${\displaystyle \mathrm{\Delta }{\ell }_{\mathrm{i}\mathrm{s}\mathrm{o}}=-⟨\ell (0)⟩}$ (see blue arrow in Figure 4E). By a geometric construction (Figure 4—figure supplement 1D) one finds that the corresponding isogonal deformation has an aspect ratio of 3, in good agreement with observations (Figure 4C’).

To find the tension threshold for active T1s, we calculate ${\displaystyle {\ell }_{\mathrm{r}\mathrm{e}\mathrm{f}}}$ using the Voronoi construction based on the vertices of the tension triangulation (see Appendix 3 for details). The length of a quartet’s central interface in the Voronoi reference state is determined by the pair of tension triangles corresponding to the cell quartet as illustrated in Figure 4D. For simplicity, we consider the case of two identical, isosceles tension triangles which are fully characterized by the relative tension ${\displaystyle T=2\sin (\varphi /2)}$ on the central interface, with ${\displaystyle \varphi }$ the angle in the tension triangle. (The general case of asymmetric tensions is discussed in the companion paper Claussen et al., 2024, 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 ${\displaystyle {\ell }_{\mathrm{r}\mathrm{e}\mathrm{f}}={\ell }_{0}T\cot (\varphi )}$. In the absence of isogonal deformations (${\displaystyle \mathrm{\Delta }{\ell }_{\mathrm{i}\mathrm{s}\mathrm{o}}=0}$), the interface length vanishes when ${\displaystyle \varphi =\pi /2}$, implying the critical relative tension ${\displaystyle T_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}^0=T_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}(\mathrm{\Delta }{\ell }_{\mathrm{i}\mathrm{s}\mathrm{o}}=0)=\sqrt{2}\approx 1.41}$. One can understand the value ${\displaystyle \sqrt{2}}$ for the critical relative tension directly from the circumcircle construction of the Voronoi tessellation as shown in Figure 4D’. ${\displaystyle {\ell }_{\mathrm{r}\mathrm{e}\mathrm{f}}}$ vanishes when the two adjacent tension triangles share the same circumcircle. In the case of isosceles triangles, this implies that they form a square such that we can read off the critical relative tension from the length of the diagonal.

In general, the condition ${\displaystyle {\displaystyle {\ell }_{\mathrm{r}\mathrm{e}\mathrm{f}}(T)=-\mathrm{\Delta }{\ell }_{\mathrm{i}\mathrm{s}\mathrm{o}}}}$ for the neighbor exchange defines the ‘T1 threshold’ in the ${\displaystyle \mathrm{\Delta }{\ell }_{\mathrm{i}\mathrm{s}\mathrm{o}}}$-${\displaystyle T}$ plane (see Figure 4E). This diagram quantifies how active and passive forces interact to drive cell intercalations. Active tension dynamics and passive isogonal strain appear as orthogonal ways to drive T1s, as illustrated by the red and blue arrows.

In the following section, we will generalize the tension–isogonal decomposition to the tissue scale. From this analysis, one can then estimate a local average value for ${\displaystyle \mathrm{\Delta }{\ell }_{\mathrm{i}\mathrm{s}\mathrm{o}}}$ and then read off the T1 threshold from Figure 4E. This will allow us to explain the value ${\displaystyle T_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}\approx 1.53}$.

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.

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 of the interfaces which will collapse, and hence, the direction of tissue flow. Therefore, tissue-scale anisotropy of active tension is central to driving and orienting convergent extension flow (Rauzi et al., 2008; Etournay et al., 2015; Lau et al., 2015; Streichan et al., 2018).

We assess tension anisotropy by locally averaging the anisotropy of inferred tensions ${\displaystyle {\mathbf{T}}_{ij}}$ at a given individual vertices as illustrated in Figure 5A (see Appendix 2 for details). The inferred tension anisotropy (Figure 5B) and its time course (Figure 5C) are consistent with previously published experimental observations: During VF invagination, DV-oriented interfaces in the lateral ectoderm are stretched causing myosin recruitment (Gustafson et al., 2022), increasing tension anisotropy. Remarkably, we find strong DV alignment of tension anisotropy in the trunk already before VF invagination (Figure 5D, 0 min). This supports our hypothesis that tension anisotropy is set up in the initial condition. As GBE progresses, the DV alignment of tension anisotropy decreases (Figure 5C and D, 37 min). Numerical simulations presented in the companion paper (Claussen et al., 2024) reproduce this loss of global tension anisotropy and show that it is responsible for the slowdown of GBE after a twofold extension.

Figure 5

Download asset Open asset

We began our investigation with a ‘tissue tectonics’ analysis (Blanchard et al., 2009) 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 (Figure 4B). Like 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 Figure 6A and Appendix 3 for details). Locally averaging over nearby vertices then yields the tissue-scale isogonal deformation (see Figure 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 the purple shaded region in Figure 6B). The lateral ectoderm further dorsal also accumulates some isogonal extension along the DV axis (green shaded region, see time traces in Figure 6C). Specifically, from the DV-DV component of the isogonal strain tensor, we can estimate the average isogonal contribution to the length of DV-oriented interfaces ${\displaystyle ⟨\mathrm{\Delta }{\ell }_{\mathrm{i}\mathrm{s}\mathrm{o}}⟩\approx 1.9\mu \mathrm{m}}$ in the lateral ectoderm. Using this value, one can read off the predicted relative tension threshold ${\displaystyle T_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}\approx 1.53}$ in from the green curve in Figure 4. This value is in excellent agreement with the value found by tension inference (Figure 2D), thus validating the model for active T1s. An immediate prediction from this model is that abolishing VF invagination will eliminate isogonal stretching of the lateral ectoderm and thereby shift the T1 threshold to ${\displaystyle T_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}=\sqrt{2}}$. To test this prediction, we analyzed a light-sheet recording of a snail mutant from the Drosophila morphodynamic atlas (Mitchell et al., 2022). This mutant lacks a VF as illustrated in Figure 7A (Gheisari et al., 2020). As expected, we find almost no isogonal strain in the lateral ectoderm of the snail mutant (see Figure 7—figure supplement 1). Analyzing the dynamics of relative tension in intercalating quartets, we find ${\displaystyle T_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}\approx \sqrt{2}}$ (see Figure 7B), confirming that tension–isogonal decomposition faithfully captures the interplay of internally generated stresses and external forces during T1 transitions.

Figure 6 with 1 supplement see all

Download asset Open asset

Figure 7 with 1 supplement see all

Download asset Open asset

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 (Figure 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 (Stern et al., 2022) and low myosin density in the amnioserosa (Streichan et al., 2018). There is also significant isogonal strain just anterior of the PMG (highlighted by the dashed outline in Figure 6B’), as well as on the lateral side of the PMG, suggesting that the convergent extension tissue flow near the dorsal pole (Streichan et al., 2018) exerts a pulling force on the PMG. In this picture, the tissue on the dorsal pole is stretched along the DV axis by active convergent extension in the germ band and contracts along the AP axis to conserve total cell area, leading to an anterior-ward pulling force on the PMG. Quantitative re-analysis of previously published cauterization experiments (Collinet et al., 2015) provides independent evidence of such a pulling force (see Appendix 6).

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

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 Figure 8A, center (Figure 2C’). In contrast, when several connected interfaces are under high tension, they form a tension cable, contraction of which leads to rosette formation (Blankenship et al., 2006; Harding et al., 2014). 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 Figure 8B): The elementary motif of the alternating tension pattern is a ‘tension bridge’ where a high-tension interface is surrounded by low-tension interfaces (see Figure 8A, right), corresponding to an obtuse triangle (Figure 8B, top right). In contrast, an acute tension triangle (Figure 8F, bottom right) corresponds to the local motif of a tension cable, where neighboring interfaces along the cable are under high tension while those transverse to the cable are under low tension (see Figure 8A, left). We can, therefore, use the tension triangle shape to define a LTC parameter. The LTC-parameter space (i.e. the space of triangle shapes), shown in Figure 8B, is spanned by two axes, measuring how elongated and how acute or obtuse the tension triangle is (see Claussen et al., 2024, Figure 8—figure supplement 1 and Appendix 4 for details). These axes correspond to the magnitude of local tension anisotropy and the cable vs bridge character of the local tension configuration, respectively.

Figure 8 with 2 supplements see all

Download asset Open asset

The geometric condition for neighbor exchanges discussed above defines a threshold in the LTC-parameter space (dashed line in Figure 8B, see Claussen et al., 2024 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 Figure 8C). The bridge fraction reaches its maximum around 25 min, just at the onset of GBE (see Figure 8D). The corresponding alternating pattern of tensions is clearly visible in Figure 8A (center). This suggests that the pattern of tensions on the (sub-)cellular scale is biologically relevant to choreographing parallel T1s. As tension anisotropy increases (Figure 5C), the local tension configurations eventually hit the T1 threshold where a cell neighbor exchange occurs, causing a reconfiguration of local tensions. Therefore, tension configurations beyond the T1 threshold are strongly suppressed (hatched region in Figure 8C). We observe that as cell intercalations take place in the lateral ectoderm, the fraction of tension bridges decreases and the fraction of tension cables increases (see Figure 8C, 50 min; and Figure 8D). In fact, the distribution of tension triangle shapes appears to approach that of a random Delaunay triangulation (see Figure 8—figure supplement 2), a completely disordered distribution. This suggests that the triangle edge flips can be statistically understood as a random ‘mixing’ process that generically causes a loss of LTC order and global alignment of tension anisotropy. This result, as well as the findings of the tissue-scale analysis more generally, are reproduced by tissue-scale simulations of our minimal model, as shown in the companion paper Claussen et al., 2024.

Recall that in a snail mutant, the lack of DV-stretching of the lateral ectoderm implies that T1s occur at a reduced tension anisotropy compared to the wild type. Figure 8B shows that the shifted T1 threshold affects tension bridges more than tension cables. As a result, we expect that tension-bridge configurations undergo T1s more rapidly in a snail mutant and, therefore, contribute less to the overall distribution of local tension configurations. Indeed, no increase in the fraction of tension bridges is observed in this mutant (orange line in Figure 8D).

We have developed a novel perspective on tissue dynamics based on the 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 (Noll et al., 2017) 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 remodeled, resulting in spatial displacement of material points. However, the focus on the force balance geometry associated with the internal structure is more useful as it holds the key to understanding possible mechanisms of active control. We find two mechanically distinct modes of tissue deformation: a cortical-tension-driven mode and an isogonal (angle preserving) mode that is unconstrained by the cortical force-balance geometry. The former represents active remodeling of interfaces, e.g., by myosin recruitment, while the latter accounts for passive tissue deformations that are controlled by external (non-local) forces and cell shape elasticity. The tension–isogonal decomposition quantifies the contributions of these two modes in experimental data based on cell geometry alone. This allows us to disentangle whether deformations are due to locally or non-locally generated forces, i.e., whether they are active or passive. Importantly, tension–isogonal decomposition is based on the same assumptions as tension inference which has been validated extensively in previous literature (Ishihara et al., 2013; Kong et al., 2019; Noll et al., 2020; Cheikh et al., 2022).

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

This predicts that the local relative tension threshold will be lower (namely ${\displaystyle \sqrt{2}}$) in mutants that lack a ventral furrow (twist and snail Gheisari et al., 2020). Indeed, analysis of a snail mutant embryo confirmed this prediction and thus validates that tension–isogonal decomposition captures the interplay of internal and external forces acting on a tissue.

Where the forces driving tissue deformations during morphogenesis originate is a fundamental question of developmental biology. In the context of Drosophila gastrulation, it has been has been intensively debated whether the germ band elongates due to internally generated stresses (Irvine and Wieschaus, 1994; Bertet et al., 2004; Rauzi et al., 2008; Butler et al., 2009; Gustafson et al., 2022; Streichan et al., 2018) or due to external pulling by the posterior midgut invagination (Collinet et al., 2015; Lye et al., 2015; Bailles et al., 2019; Gehrels et al., 2023). Experimental evidence shows that both processes contribute, making a nuanced, quantitative analysis necessary. We, therefore, re-analyzed microscopy data from previously published cauterization and mutant experiments (Collinet et al., 2015) to quantify tissue flow (see Appendix 6). The results from this quantitative analysis indicate that forces generated in the germ band contribute significantly to tissue flow. This conclusion is further supported by the observations from mutants where posterior midgut invagination is disrupted (fog, torso-like, scab, corkscrew, ksr). In these mutants, the germ band buckles forming ectopic folds (Zusman and Wieschaus, 1985; Parks and Wieschaus, 1991) or twists into a corkscrew shape (Perkins et al., 1992; Smits et al., 2023) as it extends, pointing towards a buckling instability characteristic of internally driven extensile flows (Liu et al., 2006; Senoussi et al., 2019). This suggests that the main effect of PMG invagination on the germ band lies not in creating pulling forces, but rather in ‘making room’ to allow for its orderly extension.

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

The internally driven nature of germ-band elongation flow immediately raises the question of how force generation is coordinated across cells to drive coherent tissue flow. In other words, how can the 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 (Heisenberg and Bellaïche, 2013; Collinet and Lecuit, 2021; Pinheiro and Bellaïche, 2018). Mechanical stresses and strains can propagate over long distances and contain information about the tissue geometry (e.g. in the form of hoop stresses Lefebvre et al., 2023). In an epithelium dominated by cortical tensions, this mechanical environment forms a ‘tension space,’ linked to physical space via force balance. Tension space takes the form of a triangulation which allows an intuitive visualization. The angles in the cell array are fixed by complementary angles in the tension triangulation. Thus, the tensions are geometric dials cells can directly sense and control.

We found that cells can control their configuration in tension space by defining the local dynamics of cortical tension. In experimental data from gastrulating Drosophila embryos, we identified two distinct behaviors: (i) amplification of tension on the interfaces that are already most tense, suggestive of positive tension feedback (observed in lateral ectoderm) and (ii) apparent tension homeostasis (observed in the amnioserosa). Mechanical homeostasis is found in various systems and organisms (Humphrey, 2008; Latorre et al., 2018; Stamenović and Smith, 2020). We hypothesize that tension homeostasis allows the amnioserosa to undergo large cell deformations while maintaining tissue integrity (Latorre et al., 2018; Jodoin et al., 2015). 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 spatiotemporal dynamics of tissue during morphogenesis. Evidence for positive feedback is provided by the nonlinearly increasing inferred tension on contracting interfaces (Figure 2D) and by laser ablation experiments in earlier work (Blankenship et al., 2006). The underlying molecular mechanisms are yet to be identified and might rely on the catch-bond behavior of myosin Veigel et al., 2003; Laakso et al., 2008 or mechano-sensitive binding of ${\displaystyle \alpha }$-catenin (Barry et al., 2014; Mei et al., 2020).

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

Our mathematical model for an intercalating cell quartet also sheds light on the long-standing puzzle of how interfaces extend after a neighbor exchange, even though cortical myosin can only exert contractile forces. We show that no additional active mechanism is necessary. Our minimal model captures interface extension as a purely passive process. In our model, interface extension during both active and passive intercalations is driven by the tension on adjacent interfaces which exceeds the low active contractility of the new edge. The low myosin level on the emerging interface is predicted by a simple model for the myosin ‘handover’ during a neighbor exchange. In tissue-scale simulations (Claussen et al., 2024), 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 (Collinet et al., 2015; Vanderleest et al., 2022).

Our findings on the dynamics of interface length and tension during active T1s (Figure 2D) are similar to another recently proposed model (Sknepnek et al., 2023). This model extends the classical ‘vertex model’ with area–perimeter 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 Sknepnek et al., 2023 and might be obtained from the latter in the limit of fast mechanical relaxation and dominant active tensions.

In other recent theoretical work, Drosophila extension germ band has been interpreted as resulting from tissue fluidization (Wang et al., 2020), i.e., the loss of resistance to shear stress due to the vanishing of cortical tensions (Farhadifar et al., 2007; Yan and Bi, 2019). However, laser ablation experiments (Ma et al., 2009; Collinet et al., 2015) show that cortical tensions remain finite during the morphogenetic process. In fact, interfaces under vanishing would generically buckle into a wiggly shape (as observed, for instance, in the amnioserosa during dorsal closure Tah et al., 2023 and in mutants with increased medial myosin Latorre et al., 2018). This is not observed in the Drosophila ectoderm during gastrulation where interfaces instead appear taught straight, serving as a simple visual indicator of finite tension. We have, therefore, investigated tissue mechanics in the regime where the tissue internally generated tension at all times and hence is a solid rather than a fluid. Tissue flow occurs as a result of internally driven plastic rearrangements, and does not require or imply tissue fluidization. The capacity of epithelial tissue to internally rearrange in a controlled manner – which may be thought of as active plasticity – falls beyond the fluidization paradigm. The results of our analyses on the cell and tissue scale corroborate this point of view. We return to the question of fluid vs solid behavior in the companion paper (Claussen et al., 2024).

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 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 the loss of LTC order the cause for the termination GBE? These questions are addressed in a companion paper (Claussen et al., 2024) 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 leads to the loss of LTC order as GBE progresses (see illustration in Figure 9B). Together with the 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.

Figure 9

Download asset Open asset

In addition to facilitating rapid and efficient convergent extension flow, the coordination of T1s also helps preserve the genetic patterns that are set up before germ band extension (Nüsslein-Volhard and Wieschaus, 1980) (see Appendix 1.4 and Figure 1—figure supplement 3). These patterns have been found to have sufficient fidelity to resolve the positions of individual cells along the AP axis (Dubuis et al., 2013) and are, therefore, highly susceptible to disorder resulting from incoherent cell rearrangements. This suggests that high-fidelity patterning and highly coordinated cell rearrangements might have co-evolved to maximize the speed of development.

Because the 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 Figure 9A). To achieve that, the mechanical properties of the tissue must be modulated along the dorsoventral axis: Positive tension feedback drives 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 (Gustafson et al., 2022) has shown that mechanical feedback loops in the embryo are patterned along the DV axis. Notably, the DV patterning system is conserved across species (Holley and Ferguson, 1997). Our work also highlights a second interaction between tissue flow and DV patterning: because the coherent T1s we observe do not mix cells in the tissue (Figure 1A), distinct DV fates (e.g. neural and surface ectoderm) remain clearly demarcated, a clear biological necessity.

In addition to the DV gradient, tension anisotropy along the DV axis is required to orient the GBE flow along the AP axis. This anisotropy might be due to mechanics (pulling of the ventral furrow [Gustafson et al., 2022] and hoop stresses resulting from ellipsoidal embryo geometry [Lefebvre et al., 2023] since we observe anisotropy already prior to VF formation), or indeed due to AP-striped genetic patterning (Butler and Wallingford, 2017; Lavalou et al., 2021). 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 (Gustafson et al., 2022).

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 (Claussen et al., 2024), we show that mechanical self-organization via tension feedback also allows cells to coordinate their behaviors (such as active T1s) on the tissue scale to give rise to coherent morphogenetic flow. As of yet, no genetic mechanism has been found that produces cell scale coordination – as manifested, for instance, in the alternating pattern of high and low tensions (Figure 8A, center) – from predetermined genetic patterns alone. Moreover, many interfaces in the lateral ectoderm that undergo T1 transitions at late stages of GBE are initially oriented along the AP, rather than the DV axis and only rotate into DV alignment briefly before they collapse (see Figure 1—figure supplement 2). This shows that an initially set-up genetic pattern is not sufficient to explain the coordination of cell rearrangements, in line with evidence from tissue scale analysis (Lefebvre et al., 2023). Instead, the emerging overarching picture is that of morphogenesis driven by an interplay of bottom-up local self-organization controlled by top-down genetic patterning that sets initial conditions and modulates parameters on the tissue scale (Schweisguth and Corson, 2019; Lenne et al., 2021). Indeed, this interplay has recently been shown to underlie posterior midgut invagination, where genetic patterning initiates and channels a wave that mechanically propagates due to a local feedback mechanism acting on myosin (Bailles et al., 2019).

Going forward, we expect that geometric insight into tension space will provide a new perspective on many different systems (Shindo, 2018), for example, Xenopus neural tube formation (Butler and Wallingford, 2018), amniote primitive streak formation (Voiculescu et al., 2007; Saadaoui et al., 2020; Rozbicki et al., 2015), wing disk elongation (Dye et al., 2021; Dye et al., 2017), the sea urchin archenteron (Hardin and Weliky, 2019), or kidney tubule elongation (Lienkamp et al., 2012). Furthermore, it will be interesting to study intercalations driven by actin ‘protrusions’ (Weng et al., 2022): they too are governed by cortical force balance on the cellular level, making tension inference and the tension–isogonal decomposition applicable.

Based on segmented and tracked cells, the tissue strain rate and cell deformation rate can be calculated via the methods described in Blanchard et al., 2009; see also Guirao et al., 2015 and Merkel et al., 2017. In brief, the local tissue strain rate is calculated from the centroid positions ${\displaystyle {\mathbf{c}}_i}$ of the cells in a small neighborhood ${\displaystyle N}$ of a cell (we use nearest neighbors). One looks for a linear transformation that maps the positions of the cell centroids at one-time point 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 1—figure supplement 1A–B’ show the tissue and cell deformation rates during VF invagination (Figure 1—figure supplement 1A, A’) and GBE (Figure 1—figure supplement 1B, B’). Figure 1—figure supplement 1C shows the cell elongation as measured by the Beltrami coefficient ${\displaystyle \mu =|a-b|/(a+b)}$, where ${\displaystyle a}$ and ${\displaystyle b}$ are the half-axes of a best-fit ellipse to the cell.

Figure 5Figure 6

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 ${\displaystyle \mathbf{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 the basis vector corresponding to the smallest singular value.

Depending on the context, we use three different coordinate bases in the local tangent plane: (i) The ‘Eulerian’ (or ‘lab frame’) basis for a cell ${\displaystyle i}$ is formed by the azimuthal vector ${\displaystyle {\mathbf{e}}_{i,\mathrm{D}\mathrm{V}}^{\mathrm{E}}}$ along the DV axis and its orthogonal complement ${\displaystyle {\mathbf{e}}_{i,\mathrm{A}\mathrm{P}}^{\mathrm{E}}=\mathbf{n}\times {\mathbf{e}}_{i,\mathrm{D}\mathrm{V}}^{\mathrm{E}}}$, along the AP axis. The Eulerian basis vectors are illustrated in Appendix 1—figure 1 on the left (0 min).

(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 (Figure 5C) and the orientation of collapsing interfaces (Figure 1—figure supplement 2). This co-rotating frame is defined by tracking a patch of cells neighboring cell ${\displaystyle i}$ and fitting a linear transformation ${\displaystyle A_i(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 ${\displaystyle A_i.{\mathbf{n}}_i(0)={\mathbf{n}}_i(t)}$. A polar decomposition ${\displaystyle A_i(t)=S_i(t).R_i(t)}$ then yields the local rotation matrix ${\displaystyle R_i(t)}$ for the tissue patch. The co-rotated basis vectors for cell ${\displaystyle i}$ are obtained by applying ${\displaystyle R_i}$ to the cell’s Eulerian basis vectors at time zero: ${\displaystyle {\mathbf{e}}_{i,\mathrm{A}\mathrm{P}}^{\mathrm{L}}(t)=R(t).{\mathbf{e}}_{i,\mathrm{A}\mathrm{P}}^{\mathrm{E}}(0)}$, ${\displaystyle {\mathbf{e}}_{i,\mathrm{D}\mathrm{V}}^{\mathrm{L}}(t)=R(t).{\mathbf{e}}_{i,\mathrm{D}\mathrm{V}}^{\mathrm{E}}(0)}$. The Lagrangian basis vectors are illustrated in Appendix 1—figure 1 on the right (35 min). (iii) For the tension–isogonal decomposition of intercalating cell quartets (Figure 4C and C’ and Figure 4—figure supplement 1), 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 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.

Cell interface collapse and emergence events can be identified by tracking cell contacts. We start by identifying all pairs of cells that share an interface at some time point. 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 Figure 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.

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 (Figure 1—figure supplement 1D). In the lateral ectoderm (regions 4–8 in Figure 2—figure supplement 1A), we find an average of 1.29 interface collapses per cell. From this, we can estimate the total tissue extension by cell rearrangements. In a hexagonal lattice, synchronous T1s cause one interface collapses per cell and elongate the tissue by a factor ${\displaystyle \sqrt{3}}$, assuming that the cells do not elongate (The tissue also contracts by a factor ${\displaystyle 1/\sqrt{3}}$ along the orthogonal axis). For 1.29 interface collapses per cell, therefore, cause tissue elongation by a factor ${\displaystyle 1.29\times \sqrt{3}\approx 2.23}$, in good agreement with the total elongation of the germ band (Figure 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 Figure 1—figure supplement 1C). This is a generic feature of passive T1s driven by isogonal deformations (see Figure 4—figure supplement 1D).

The convergent extension flows in the lateral ectoderm and near the dorsal pole (amnioserosa) are oriented perpendicular to each other (see Figure 1—figure supplement 1). We, therefore, also expect perpendicular orientation of the collapsing interfaces in these two regions. We quantify the orientation of collapsing edges at the blastoderm stage (i.e. before the onset of tissue flow, ${\displaystyle t=0\mathrm{m}\mathrm{i}\mathrm{n}}$) as shown in Figure 1—figure supplement 2A and C. In the ectoderm, interfaces that collapse before ${\displaystyle {\displaystyle t=42\mathrm{m}\mathrm{i}\mathrm{n}}}$, i.e., within the first ca. 15 min of GBE are predominantly DV oriented at ${\displaystyle {\displaystyle t=0\mathrm{m}\mathrm{i}\mathrm{n}}}$ (see Figure 1—figure supplement 2B). By contrast, interfaces near the dorsal pole are predominantly oriented along AP. The orientation of interfaces in the lateral ectoderm that collapse after ${\displaystyle {\displaystyle t=42\mathrm{m}\mathrm{i}\mathrm{n}}}$ is initially not strongly biased along the DV axis Figure 1—figure supplement 2C and D. However, during GBE, these interfaces progressively align with the DV co-rotated DV axis before they eventually start contracting (see Figure 1—figure supplement 2C’ and D; the coordinate frame co-rotating with the tissue is described in Appendix 1.1).

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 array at each time point. We exclude cases where the four cells surround one or more cells that are not part of the cycle and cases where one of the cells is surrounded by the three others. Notably, the quartets don’t necessarily exist for the entire time series. Namely, when one of the ‘inner’ interfaces in the quartet, other than the central interface, collapses, the quartet is destroyed.

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

Using this procedure, we identified 2231 T1 events, 498 fivefold rosettes, 40 sixfold rosettes, 4 sevenfold rosettes, and no higher order rosettes (see Figure 1E’). Because the number of rosettes is small compared to the number of T1s, we have focused on T1 events for the subsequent analysis. Expressed in terms of a participation ratio of interface contractions in T1s vs rosettes, we find a ratio ${\displaystyle 2231:(2\cdot 498+3\cdot 40+4\cdot 4)\approx 2:1}$, matching observations from previous studies (Farrell et al., 2017; Stern et al., 2022).

During the early development of the Drosophila embryo a hierarchy of genetic patterns is set up before morphogenetic flow starts. This hierarchy culminates in the striped expression of the pair-rule genes which lay out the segmented body plan of the larva (Nüsslein-Volhard and Wieschaus, 1980). Detailed quantification of the gene expression profiles has revealed a remarkable fidelity of this patterning system which is sufficient to identify each cell’s position along the AP axis (Dubuis et al., 2013). Germ band extension starts after this pattern has been established, raising the question of how cell rearrangement affects the pattern fidelity. Cell tracking allows us to test this question computationally. To visualize the effect of cell rearrangements on an AP pattern, we assign colors to cells according to their initial AP position, coarse-grained into ca 2 cell-wide stripes (Figure 1—figure supplement 3A, 0 min). Rendering the cells with these colors at a later time shows sharp boundaries between stripes with few outliers, revealing that the stripe pattern is locally preserved with high fidelity (Figure 1—figure supplement 3A, 40 min). On large scales, the pattern is distorted because the posterior germband curves over the posterior pole as it elongates (Figure 1A).

The local preservation of pattern fidelity is in line with the highly choreographed active T1s in the germ band (Figure 1C). As a reference, we simulated purely passive cell rearrangement driven by externally imposed strain (Figure 1—figure supplement 3B). In brief, cells were simulated as foam bubbles with constant interfacial tension and area-elasticity. The externally imposed strain was implemented by forces applied to the top and bottom boundaries. Notably, the stripe boundaries after twofold elongation appear significantly more disrupted than in the embryo, suggesting that coordination of T1s in the germ band helps preserve pattern fidelity. Simulations were carried out as described in Claussen et al., 2024.

To quantify the local pattern fidelity we define the relative AP position for cell pairs ${\displaystyle (ij)}$ as

where ${\displaystyle {\mathbf{c}}_i(t)}$ is the centroid position of cell ${\displaystyle i}$ at time ${\displaystyle t}$ we project onto the advected (Lagrangian) local basis vectors ${\displaystyle {\mathbf{e}}_{i,\mathrm{A}\mathrm{P}}^{\mathrm{L}}(t)}$ defined inAppendix 1.2 to measure relative AP positions in relation to the original blastoderm stage coordinate frame. From the ${\displaystyle d_{ij}^{\mathrm{A}\mathrm{P}}(t)}$ we define the correlation coefficient

where ${\displaystyle ⟨\dots ⟩}$ denotes the average over all pairs that are neighbors at time ${\displaystyle t}$. The correlation ${\displaystyle C(t)}$ plotted against the number of interfaces lost per cell measures how well the local gradient of an AP pattern is preserved during convergent extension (Figure 1—figure supplement 3C). For the anterior germ band, we find significantly less degradation of fidelity compared to the passive stretching simulation. The posterior germ band, which is distorted as it curves over the posterior pole, exhibits a similar degradation rate as the passive reference case. Notably, the rate of pattern degradation (i.e. the slope of the curve) increases towards the end in the anterior germ band, reaching the same degradation rate as the posterior germ band and passive reference. This is consistent with the loss of coherence of T1s towards the end of germ band extension which likely results from the loss of large-scale tension anisotropy (Figure 5C and D) and LTC order (Figure 8C and D).

Appendix 1—figure 1

Download asset Open asset

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

where ${\displaystyle T_{ij}}$ are the cortical tensions. Multiplying the above equation with ${\displaystyle {\hat{\mathbf{e}}}_{12}}$ and ${\displaystyle {\hat{\mathbf{e}}}_{23}}$ (where ${\displaystyle \hat{\mathbf{e}}}$ denotes the vector orthogonal to ${\displaystyle \mathbf{e}}$ and with the same magnitude) yields two equations

Here, ${\displaystyle {\alpha }_{ijk}}$ denotes the angle subtended by the interfaces ${\displaystyle (ij)}$ and ${\displaystyle (jk)}$. An alternative, more geometric way to arrive at these equations is by applying the law of sines to the tension triangle, using that its angles are complementary to those at the vertex and ${\displaystyle \sin (\alpha -\pi )=\sin (\alpha )}$.

Eqs. (7) 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 ${\displaystyle \mathcal{N}=\sum _{(ijk)}\sin {\alpha }_{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 (Noll et al., 2017). This is because force balance implies that the tension vectors ${\displaystyle {\mathbf{T}}_{ij}=T_{ij}{\mathbf{e}}_{ij}}$ form a triangulation (see below).

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

We solve these equations under the normalization condition that the average of the ‘outer tensions’ ${\displaystyle T_{12},T_{23},T_{34}}$ and ${\displaystyle T_{14}}$ is 1. The inferred value for ‘inner’ cortical tension ${\displaystyle T_{13}}$ then is that relative to the average tension of its neighbors. We call this the relative tension of an interface throughout the manuscript.

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

For the analysis presented in Figure 2D and Figure 2E, collapsing and emerging interfaces were identified separately as described in Appendix 1.3. All events before ${\displaystyle {\displaystyle 12\mathrm{m}\mathrm{i}\mathrm{n}}}$ were excluded from the analysis.

When an interface is short, the interface angles become highly susceptible to fluctuations in the vertex positions. We, therefore, excluded all time points where the central interface or one of its neighbors is shorter than a threshold length (${\displaystyle 0.75\mu \mathrm{m}}$). This value was chosen because the vertex positions fluctuate about ${\displaystyle 0.5\mu \mathrm{m}}$ between time points. 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.

The time series of interface length and relative tension are then aligned such that the interface collapse (respectively emergence) is at ${\displaystyle t=0}$. Figure 2—figure supplement 1 shows the pooled data resolved by the position of the quartet along the DV axis.

Because the angles in the tension triangulation are complementary to the angle at vertices in the cell array, one can find a unique tension triangulation (up to a global scale factor) for a given physical cell array in force balance (see dashed gray arrow in Figure 4—figure supplement 1A). However, the reverse is not true. For a given tension triangulation, there is a continuous family of cell arrays that share the same tension triangulation and that are related by isogonal deformations which change the interface lengths collectively while preserving the angles. (Mathematically, there are fewer triangulation vertices than cell vertices, so the triangulation leaves one isogonal mode degree of freedom per cell undetermined. At the end of this section, will define an ‘isogonal potential’ that parametrizes the isogonal modes. Specifically, we find that a quadratic potential corresponds to a uniform pure shear of the cell centroids that preserves the angles at vertices).

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

In the following, we discuss the details of the tension–isogonal decomposition. To achieve a well-defined decomposition of tissue deformation, we need to specify a reference tessellation that is uniquely determined by the tension triangulation. There will then be a unique isogonal deformation that deforms the reference configuration into the physical cell array. Fortunately, we do not need to explicitly construct the reference configuration. Rather, we use the triangulation formed by the cell centroids to quantify the physical tissue shape (Merkel et al., 2017). (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.

To make this construction concrete, we first consider a single tri-cellular vertex, with the associated tension triangle ${\displaystyle {\hat{\mathbf{T}}}_{ij}}$ and the centroids ${\displaystyle {\mathbf{c}}_i}$. The indices ${\displaystyle i,j=1,2,3}$ label the three cells that meet at the vertex (see Figure 6—figure supplement 1A’). The tension triangle vectors ${\displaystyle {\hat{\mathbf{T}}}_{ij}}$ are obtained by local tension inference based on the angles at the vertex (see Appendix 2). The orientation of the tension vectors is determined by demanding that they are at right angles to the cell-cell interfaces ${\displaystyle {\mathbf{e}}_{ij}}$. To fix the tension scale, we normalize the area of the tension triangle to unity. Fixing the tension triangle area is justified because it maintains the total area of the tension triangulation such that cell dilation/contraction is purely isogonal. The numerical value of the tension triangle area is immaterial since it can be absorbed in the scale factor ${\displaystyle {\ell }_0}$.

We define the local isogonal deformation tensor ${\displaystyle {\mathbf{I}}_{\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{x}}}$ such that it transforms the tension triangle into the centroidal triangle. Mathematically, this definition is written as the equations ${\displaystyle {\ell }_0{\mathbf{I}}_{\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{x}}.{\hat{\mathbf{T}}}_{ij}={\mathbf{c}}_j-{\mathbf{c}}_i}$, ${\displaystyle i\ne j\in \{1,2,3\}}$. Since there are two independent displacement vectors in tension space and physical space each, these equations represent four independent constraints that uniquely determine the four components of the ${\displaystyle 2\times 2}$ matrix ${\displaystyle {\ell }_0{\mathbf{I}}_{\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{x}}}$. The scale factor ${\displaystyle {\ell }_0}$ converts units of tension (a.u.) into units of length. How we fix ${\displaystyle {\ell }_0}$ is described below.

The isogonal deformation tensor ${\displaystyle {\mathbf{I}}_{\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{x}}}$ can be decomposed into a symmetric tensor ${\displaystyle \mathbf{U}}$ and a rotation ${\displaystyle \mathbf{R}}$ using polar decomposition ${\displaystyle \mathbf{U}\mathbf{R}={\mathbf{I}}_{\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{x}}}$. The appearance of a rotational component in the isogonal deformation might appear surprising at first glance since isogonal modes leave angles at vertices invariant. However, the centroidal triangle can be sheared relative to the tension triangle by a simple shear – composition of pure shear with a rotation – while the angles at the vertex remain invariant.

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 ${\displaystyle 20\mu \mathrm{m}}$, such that each grid facet contains ca 10–15 cells (grid facets near the poles contain fewer cells because of the distortion caused by the pullback projection). We calculate ${\displaystyle {\mathbf{I}}_{\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{x}}}$ for all tri-cellular vertices in a grid facet and take the average to find ${\displaystyle {\mathbf{I}}_{\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}\mathrm{t}}}$.

Figures Figure 4D and Figure 6—figure supplement 1B–D (top panels) show the isogonal strain tensor ${\displaystyle {\mathbf{I}}_{\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}\mathrm{t}}-\mathbb{I}}$, where ${\displaystyle \mathbb{I}}$ is the identity matrix. The bottom panels in Figure 6—figure supplement 1B–D show the individual, single-vertex strain tensors ${\displaystyle {\mathbf{I}}_{\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}\mathrm{t}}-\mathbb{I}}$ 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 (Figure 4E), we first transform the single-vertex isogonal deformation tensor to the local co-rotating (‘Lagrangian’) frame (see Appendix 1.1). 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, ${\displaystyle ⟨\mathrm{\Delta }{\ell }_{\mathrm{i}\mathrm{s}\mathrm{o}}⟩}$ (see Appendix 3). At the end of VF invagination (${\displaystyle {\displaystyle t=25\mathrm{m}\mathrm{i}\mathrm{n}}}$), we find ${\displaystyle ⟨({\mathbf{I}}_{\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{x}}{)}_{22}⟩\approx 1.538\pm 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 Figure 6—figure supplement 1C). We thus have ${\displaystyle ⟨\mathrm{\Delta }{\ell }_{\mathrm{i}\mathrm{s}\mathrm{o}}⟩=⟨({\mathbf{I}}_{\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{x}}{)}_{22}-1⟩⟨\ell (0)⟩\approx 1.9\mu \mathrm{m}}$ for DV-oriented interfaces, where ${\displaystyle ⟨\ell (0)⟩\approx 3.5\mu \mathrm{m}}$ is the average initial interface length.

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

When one of the angles at a vertex is close to ${\displaystyle \pi }$, the tension triangle becomes highly elongated. In fact, sometimes one of the angles at a vertex exceeds ${\displaystyle \pi }$ 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 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).

From the tension vectors ${\displaystyle {\mathbf{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 the vertex. The normalization factor ${\displaystyle {\mathcal{N}}^2=\sum _{(ij)}{T}_{ij}^2}$ is chosen such that ${\displaystyle \mathrm{Tr}\mathbb{T}=1}$. The difference of the eigenvalues of ${\displaystyle \mathbb{T}}$ quantifies the magnitude of tension anisotropy. When the interface angles at the vertex are all identical to ${\displaystyle {120}^{\circ }}$ the tension anisotropy vanishes. For non-identical vertex angles the tension is locally anisotropic and the dominant eigenvector of ${\displaystyle \mathbb{T}}$ 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 ${\displaystyle {90}^{\circ }}$ relative to the interfaces (Figure 2B).

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

The calculation for cell quartets (vertex pairs) is analogous to the case of single vertices. The two tension triangles associated with the inner vertices of the cell quartet form a ‘kite’ (see Figure 4—figure supplement 1A). We normalize tension such that the average area of the tension triangles is one. The quartet’s isogonal deformation tensor ${\displaystyle {\mathbf{I}}_{\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{t}}}$ is defined such that it transforms the vectors ${\displaystyle {\mathbf{T}}_{ij}}$ into the corresponding centroidal displacement vectors ${\displaystyle {\mathbf{c}}_j-{\mathbf{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 ${\displaystyle {\mathbf{I}}_{\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{t}}}$ 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 ${\displaystyle (ij)}$ run over cell pairs going around the quartet, i.e., ${\displaystyle (ij)=(12),(23),(34),(41)}$. Notably, with these definitions, one has the relation

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

The explicit construction of the reference cell array allows us to understand the geometry of cell-neighbor exchanges and find a criterion for when they occur (see main text for a discussion of a simple symmetric case and the companion paper [Claussen et al., 2024] 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-intersections when the triangulation fulfills the Delaunay criterion that opposite angles in adjacent triangles must sum to a value smaller than ${\displaystyle \pi }$. 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 (Efrati et al., 2009).

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

As discussed in the main text (and further elaborated in Claussen et al., 2024) our simulations show that a non-zero cell-level shear modulus is required for extension via active T1s. Here, we analyze the relationship 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 deformations can shear and not just dilate/contract cells, we make use of the’isogonal mode parametrization introduced in Noll et al., 2017. It assigns an isogonal ‘potential’ ${\displaystyle {\mathrm{\Theta }}_i}$ to each cell, and calculates the cell displacements from the ${\displaystyle {\mathrm{\Theta }}_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 ${\displaystyle {\mathbf{e}}_{ij}=-{\mathbf{e}}_{ji}}$ and denote the corresponding tensions as ${\displaystyle T_{ij}}$. We denote by ${\displaystyle \hat{\mathbf{a}}}$ the vector ${\displaystyle \mathbf{a}}$ rotated by ${\displaystyle \pi /2}$ in counterclockwise direction, such that ${\displaystyle \hat{\mathbf{a}}}$ fulfills the relations ${\displaystyle \hat{\mathbf{a}}.\mathbf{a}=0}$ and ${\displaystyle ||\hat{\mathbf{a}}||=||\mathbf{a}||}$. In force balance, the rotated force vectors ${\displaystyle {\hat{\mathbf{T}}}_{ij}=T_{ij}{\hat{\mathbf{e}}}_{ij}}$ form a triangulation (Noll et al., 2017).

The isogonal displacement ${\displaystyle {\mathbf{u}}_{ijk}}$ of the real space vertices ${\displaystyle {\mathbf{r}}_{ijk}}$ (identified by the three adjacent cells) is given by

where ${\displaystyle S_{ijk}={\hat{\mathbf{T}}}_{ij}.{\mathbf{T}}_{ik}}$ is the area of the tension triangle ${\displaystyle (ijk)}$ and cyc. denotes cyclic permutations of ${\displaystyle (ijk)}$.

First, observe that the uniform isogonal mode ${\displaystyle {\displaystyle {\mathrm{\Theta }}_i=}}$ const. has no effect on the vertex positions because

Now we aim to show that a constant gradient in ${\displaystyle {\mathrm{\Theta }}_i}$ drives a uniform displacement of the ${\displaystyle {\mathbf{r}}_{ijk}}$. Specifically, by uniform gradient, we mean ${\displaystyle {\mathrm{\Theta }}_i={\mathbf{t}}_i.\mathbf{a}}$, i.e., a linear gradient in the tension space (${\displaystyle {\mathbf{t}}_i}$ is the position of the tension triangulation vertex corresponding to cell ${\displaystyle i}$, such that ${\displaystyle {\hat{\mathbf{T}}}_{ij}={\mathbf{t}}_j-{\mathbf{t}}_i}$). 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 (${\displaystyle i=1-4}$), corresponding to a ‘kite’ in tension space (note that ${\displaystyle \hat{\mathbf{a}}.\mathbf{b}}$ is identical to the wedge product ${\displaystyle \mathbf{a}\wedge \mathbf{b}}$).

Because a constant can be arbitrarily added to all ${\displaystyle {\mathrm{\Theta }}_i}$, we can set ${\displaystyle {\mathrm{\Theta }}_1=0}$ and thus have ${\displaystyle {\mathrm{\Theta }}_i={\hat{\mathbf{T}}}_{ij}.\mathbf{a}}$ for ${\displaystyle 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 ${\displaystyle {\hat{\mathbf{T}}}_{12}}$ and ${\displaystyle {\hat{\mathbf{T}}}_{13}}$. For the latter, we find

where we used that ${\displaystyle {\hat{\mathbf{T}}}_{ij}.{\mathbf{T}}_{ij}=0}$ and applied the definition of ${\displaystyle S_{ijk}}$.

Projecting (Equation 20) and (Equation 21) onto ${\displaystyle {\hat{\mathbf{T}}}_{12}}$ gives

To show the equality of these two right-hand sides, we use that given ${\displaystyle {\hat{\mathbf{T}}}_{13}.\mathbf{a}}$ and ${\displaystyle {\hat{\mathbf{T}}}_{14}.\mathbf{a}}$ we can find ${\displaystyle \mathbf{a}}$ and substitute the result into ${\displaystyle {\hat{\mathbf{T}}}_{12}.\mathbf{a}}$. We start by ‘expanding the identity’

Explicitly writing out the inverse matrix then gives

With this, we find the relation

where we used ${\displaystyle {\mathbf{T}}_{ij}=-{\mathbf{T}}_{jk}}$ to flip the indices on ${\displaystyle {\mathbf{T}}_{14}}$. Comparing to (Equation 24) and (Equation 25) now shows the identity of their RHSs.

Taken together, we have shown that

Because ${\displaystyle {\hat{\mathbf{T}}}_{12}}$ and ${\displaystyle {\hat{\mathbf{T}}}_{13}}$ are linearly independent, it follows that ${\displaystyle {\mathbf{u}}_{123}={\mathbf{u}}_{134}}$. QED.

We just showed that a constant gradient in ${\displaystyle {\mathrm{\Theta }}_i}$ corresponds to a uniform displacement of the real space vertices ${\displaystyle {\mathbf{r}}_{ijk}}$. We can, therefore, think of ${\displaystyle \mathrm{\Theta }(\mathbf{t})}$ as a ‘potential’ for the isogonal displacement field: ${\displaystyle \mathbf{u}\approx {\mathrm{\nabla }}_{\mathbf{t}}\mathrm{\Theta }}$, where the approximation is valid for slowly varying gradients and exact for constant gradients. The gradient ${\displaystyle {\mathrm{\nabla }}_{\mathbf{t}}}$ is taken in tension space because the function ${\displaystyle \mathrm{\Theta }({\mathbf{t}}_i)={\mathrm{\Theta }}_i}$ is defined on the vertices of the tension triangulation ${\displaystyle {\mathbf{t}}_i}$.

A pure shear aligned with the coordinate axes is given by a displacement field ${\displaystyle \mathbf{u}(\mathbf{r})=\epsilon \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(1,-1).\mathbf{r}}$ and is, therefore, generated (approximately) by a isogonal potential forming a hyperbolic paraboloid ${\displaystyle \mathrm{\Theta }=\epsilon {\mathbf{t}}^{\mathrm{T}}.\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(-1,1).\mathbf{t}=\epsilon \left(t_1^2-t_2^2\right)}$.

In the main text, we argued that the 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 Figure 8—figure supplement 1A).

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 ${\displaystyle \pi }$, we can represent the space of triangle shapes in barycentric coordinates (see Figure 8—figure supplement 1B). 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 Figure 8—figure supplement 1B) 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 ${\displaystyle \mathbb{T}}$, defined in the previous subsection. To quantify the shape (acute vs obtuse), we define a second tensor

where the indices ${\displaystyle I,J\in \{1,2,3\}}$ are short-hand labels for the three interfaces that meet a tri-cellular vertex. Because ${\displaystyle \sum _I{\mathbf{T}}_J=0}$ (force balance), this tensor has the null vector ${\displaystyle (1,1,1{)}^{\mathrm{T}}}$. To get rid of this nullspace, we go to barycentric coordinates by defining the pair of vectors, ${\mathit{𝝉}}_{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 ${\displaystyle (1,1,1{)}^{\mathrm{T}}}$ of ${\displaystyle \mathcal{S}}$. Now we can define the shape tensor

Observe that

which implies that ${\displaystyle \mathbb{T}}$ and ${\displaystyle \mathbb{S}}$ share the same eigenvalues. In fact, we can perform a singular value decomposition of ${\displaystyle \mathfrak{T}}$

where ${\displaystyle R(\varphi )}$ is the rotation matrix with angle ${\displaystyle \varphi }$ and $\mathbf{𝚲}=\mathrm{diag}({\lambda }_1,{\lambda }_2)$ is a diagonal matrix of singular values whose squares ${\displaystyle {\lambda }_1^2,{\lambda }_2^2}$ are the eigenvalues of ${\displaystyle \mathbb{T}}$ and ${\displaystyle \mathbb{S}}$. Because we normalized the tension vectors, ${\displaystyle {\displaystyle \mathrm{Tr}\mathbb{T}=\mathrm{Tr}\mathbb{S}={\lambda }_1^2+{\lambda }_2^2=1}}$ 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 ${\displaystyle \psi }$ indicates whether the triangle is obtuse or acute. The reasoning to define the rotation based on ${\displaystyle \psi /2}$ will become apparent below when we represent the shape tensor by a complex number.

The shape tensor ${\displaystyle \mathbb{S}}$ has only two degrees of freedom, the difference of eigenvalues ${\displaystyle {\lambda }_1^2-{\lambda }_2^2}$ (which it shares with the tensor ${\displaystyle \mathbb{T}}$) and the angle ${\displaystyle \psi }$ that sets the orientation of its two orthogonal eigenvectors. We can represent ${\displaystyle \mathbb{S}}$ by a complex number

where ${\displaystyle |\mathrm{\Psi }|=|{\lambda }_1^2-{\lambda }_2^2|}$ is the magnitude of tension anisotropy and the phase ${\displaystyle \psi }$ is the ‘shape space angle’ defined by Equation 30.

Permutations of the edge indices ${\displaystyle \alpha }$ correspond to sign changes and rotations of ${\displaystyle \psi }$ by integer multiples of ${\displaystyle 2\pi /3}$. To mod out these group actions (which correspond to the symmetries of the barycentric shape space Figure 8—figure supplement 1B), we define

where ${\displaystyle x{\mathrm{m}\mathrm{o}\mathrm{d}}_{d}n=(x+d\mathrm{m}\mathrm{o}\mathrm{d}n)-d\in [-d,n-d]}$ is the modulo operation with offset ${\displaystyle d}$. This shape parameter informs how acute or obtuse a tension triangle is. For an acute isosceles triangle, ${\displaystyle \tilde{\psi }=0}$, whereas for an obtuse isosceles triangle ${\displaystyle \tilde{\psi }=\pi }$. For non-isosceles triangles, ${\displaystyle \tilde{\psi }}$ takes values in ${\displaystyle [0,\pi ]}$ quantifying the relative degree of acuteness vs obtuseness. For an equilateral triangle, ${\displaystyle |\mathrm{\Psi }|=0}$ and ${\displaystyle \tilde{\psi }}$ is not defined. The axes ${\displaystyle |\mathrm{\Psi }|}$ and ${\displaystyle \tilde{\psi }/\pi }$ span the space of triangle shapes shown in Figure 8—figure supplement 1.

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. ${\displaystyle T_2}$ and ${\displaystyle T_3}$ will be second highest and highest tension, respectively). We can now define a parameter

which quantifies the relative tension difference between the two largest tensions at a vertex. For a tension cable, there are two large tensions that are similar in magnitude, implying that ${\displaystyle p}$ is small. For a tension bridge, there is one high-tension edge meeting two low-tension edges so ${\displaystyle p}$ will be closer to 1. The maximal value ${\displaystyle p=1}$ is reached for ${\displaystyle T_1=T_2=T_3/2}$. Remarkably, to a very good approximation ${\displaystyle p\approx |\mathrm{\Psi }|\tilde{\psi }/\pi }$.

As a reference for the triangle shape distributions, we generated random Delaunay triangulations based on the Ginibre random point process (Ginibre, 1965; the Ginibre process models the distribution of particles with repulsion and thus generates slightly more regular point distribution than the Poisson process where points are randomly positioned independently from one another). A small sample of such a Ginibre-based random triangulation is shown in Figure 8—figure supplement 2A. The histograms in the triangle shape space Figure 8—figure supplement 2B show 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 (Figure 8—figure supplement 2D).

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.

In the following, we explain the simulation details for Figure 3, where we present a model of a single intercalating quartet. The code for the simulations shown is available online: https://github.com/nikolas-claussen/geometric_basis_of_convergent_extension_simulation (copy archived at Claussen, 2024b).

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 that 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 several simplifications occur (the general case of a disordered cell array is the subject of the companion paper Claussen et al., 2024). In this case, the cell shape is completely specified by the three vertex angles ${\displaystyle {\varphi }_i}$ and edge lengths ${\displaystyle {\ell }_i}$ at each vertex, with ${\displaystyle i\in \{0,1,2\}}$. Because the cell array is a symmetric lattice, the edge lengths ${\displaystyle {\ell }_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 ${\displaystyle T_i}$, and the passive tension ${\displaystyle T_{p,i}}$. The tensions at a vertex evolve according to:

where ${\displaystyle n>1}$ and ${\displaystyle {\tau }_{\mathrm{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 ${\displaystyle 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 ${\displaystyle {\varphi }_i={\sin }^{-1}\left(\frac{T_i}{2R}\right)}$, where ${\displaystyle R}$ is the tension triangle circumradius.

To determine the ${\displaystyle {\ell }_i}$, we minimize an elastic energy based on the cell shape tensor:

where ${\displaystyle i}$ runs over the edges of the cell ${\displaystyle \mathcal{C}}$, and ${\displaystyle {\mathbf{e}}_i}$ are the unit vectors pointing along the interfaces. 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 ${\displaystyle S}$ resides in the cell interior (cytoplasm 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 ${\displaystyle {\tilde{S}}_{\mathcal{C}}=\sum _i{\ell }_i^2{\mathbf{e}}_i\otimes {\mathbf{e}}_i}$ instead gives broadly similar results (although interface collapse happens more abruptly, because of the higher order non-linearity).

We assign a reference shape ${\displaystyle S_0}$. An isotropic ${\displaystyle S_0\propto \mathrm{I}\mathrm{d}}$ favors equilateral hexagons. We chose a slightly anisotropic reference shape ${\displaystyle S_0}$ (15% anisotropy) to model the isogonal stretching caused by the ventral furrow before the onset of GBE. This anisotropy sets the angle at which the inner interface collapses, ${\displaystyle {\ell }_1=0}$. The experimental value of the collapse tension was, therefore, used to fit the ${\displaystyle S_0}$-anisotropy.

We define the cell’s elastic energy via

with bulk modulus ${\displaystyle \lambda }$ and shear modulus µ. We used a shear/bulk ratio of ${\displaystyle \mu /\lambda =0.1}$. Because of the separation of scales between cortical tensions and elastic energy baked into the model, the absolute values of ${\displaystyle \lambda ,\mu }$ are irrelevant. Furthermore, for the case of active T1s, the results do not depend on ${\displaystyle \mu /\lambda }$, as long as ${\displaystyle \mu >0}$.

The edge lengths ${\displaystyle {\ell }_i}$ can now be determined by minimizing the elastic energy Equation 36 w.r.t. the ${\displaystyle {\ell }_i}$, for the angles determined by the tension dynamics.

Strikingly, within the single-quartet model, the ‘shape strain’ ${\displaystyle S_{\mathcal{C}}-S_0}$ can always be set to 0 by choice of the edge lengths, and the energy ${\displaystyle E_{\mathcal{C}}=0}$ throughout. This can already be seen from a degree-of-freedom count (three ${\displaystyle {\ell }_i}$ for the three independent components of ${\displaystyle S_{\mathcal{C}}-S_0}$). The elastic energy, therefore, acts only on the isogonal modes (i.e. ${\displaystyle {\mathrm{\partial }}_{{\varphi }_i}{E}_{\mathcal{C}}{|}_{{\ell }_i=\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}\mathrm{r}\mathrm{s}}=0}$) 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 ${\displaystyle E=(A-A_0{)}^2+(P-P_0{)}^2}$ (where ${\displaystyle A,P}$ are the cell area and perimeter, and ${\displaystyle A_0,P_0}$ their target values) (Bi et al., 2015), there exists an energy barrier, and the inner interface ${\displaystyle {\ell }_1}$ only collapses when ${\displaystyle {\varphi }_0=\pi }$. This is shown in Figure 3—figure supplement 1B. 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, ${\displaystyle E=(A-A_0{)}^2+\mu \mathrm{T}\mathrm{r}[(S_{\mathcal{C}}-S_0{)}^2]}$, (or perimeter elasticity with cell shape bulk elasticity) leads to similar results as Equation 36. 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 (Claussen et al., 2024).

When an interface collapses, ${\displaystyle {\ell }_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 illustrated in Figure 3D, we propose a myosin handover mechanism to explain the extension of the new interface post-intercalation. An interface with cortical tension ${\displaystyle 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 ${\displaystyle 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 Kale et al., 2018 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 Figure 3D).

Consider the tensions ${\displaystyle T_0,T_1,T_2}$ 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 ${\displaystyle (01),(12),(21)}$ (where e.g. ${\displaystyle (01)}$ is the cell abutting interfaces 0 and 1). Let ${\displaystyle m_{01},m_{12},m_{21}}$ 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:

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 ${\displaystyle m_{12}}$. 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 ${\displaystyle T_0,T_1,T_2}$ and ${\displaystyle T_0,T_1^{\mathrm{\prime }},T_2^{\mathrm{\prime }}}$. Let ${\displaystyle T_a}$ be the active tension on the new junction immediately after the T1. It is equal to

Note, however, that the total tension ${\displaystyle T_n}$ on the new junction is not necessarily equal to ${\displaystyle T_a}$. The total tension is defined geometrically from the angles at the new junction (or, equivalently, the tension triangle vertices). Indeed, generally, ${\displaystyle 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 ${\displaystyle 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 ${\displaystyle {90}^{\circ }}$, ${\displaystyle T_1=T_2=T_1^{\mathrm{\prime }}=T_2^{\mathrm{\prime }}=1}$ and ${\displaystyle T_0=T_n=\sqrt{2}}$. Therefore, ${\displaystyle T_{a,n}=2-\sqrt{2}\approx 0.6}$ and ${\displaystyle T_{p,n}=\sqrt{2}-(2-\sqrt{2})\approx 0.8}$. Note that by the triangle inequality, for any convex quadrilateral with perimeter ${\displaystyle P}$ and diagonals ${\displaystyle D_1,D_2}$, one has ${\displaystyle P/2\le D_1+D_2\le P}$. Applying this to the quadrilateral formed by the two tension triangles at the collapsing interface, we get ${\displaystyle T_{a,n}\ge 0}$ and ${\displaystyle T_{p,n}\ge 0}$: the handover formula always results in positive active and passive tensions. Furthermore, 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 ${\displaystyle {\displaystyle {\tau }_p^{-1}}}$. Combining this with the feedback equation yields Equation 34. The relaxation rate ${\displaystyle {\tau }_p={\tau }_{\mathrm{T}}/6}$ was fit to the measured tension decay post intercalation.

We integrated Equation 34 using a 4th order Runge-Kutta method as implemented in the SciPy software package (Virtanen et al., 2020), using as initial conditions the vertex angles in the experimental data at time ${\displaystyle t=5\text{min}}$ (the vertex angles from the data were temporally smoothed with a window of ${\displaystyle 2\text{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 Equation 37 and simulated the combined active/passive tension dynamics. The overall time and length scales were likewise fit to the data. The feedback exponent was ${\displaystyle n=4}$. Numerically, energy minimization was carried out using the scipy.optimize package, specifically its Nelder-Mead optimizer.

Within the same framework, we can also model passive intercalations, shown in Figure 3F. We can either specify cortical tension dynamics (in this case, tension homeostasis, ${\displaystyle \dot{T}\sim -(T-1)}$) and then minimize an elastic energy that includes the influence of externally imposed shear. This will trivially reproduce the isogonal behavior shown in Figure 2E. Alternatively one can reason that on the dorsal side, there is very little active tension (cortical myosin), so 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 Figure 3F, we show simulations based on the latter option. We implement the externally imposed shear causing the passive intercalations as follows. We calculate the four cell centroids ${\displaystyle {\mathbf{c}}_i}$ of the cell quartets and demand that at simulation time ${\displaystyle t}$

where ${\displaystyle S(t)}$ is a diagonal matrix representing (area-preserving) shear, increasing over time with constant shear rate ${\displaystyle \gamma }$.

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 ${\displaystyle {\ell }_i}$ and the angles ${\displaystyle {\varphi }_i}$. As long as the shear modulus ${\displaystyle \mu \ll \lambda }$, this procedure results in approximately isogonal behavior. We chose ${\displaystyle \mu =0.01\lambda }$. To obtain continuity across the intercalation, the rest shape ${\displaystyle S_0}$ 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 ${\displaystyle 0^{\circ }-{25}^{\circ }}$.

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

However, tissue ripping is only a qualitative measure for such pulling forces. To obtain a more quantitative estimate of the relative contributions of internally generated forces in the germ band and ‘external’ pulling forces from the PMG we re-analyzed the movies from Collinet et al., 2015 using Particle Image Velocimetry (PIVlab Thielicke and Sonntag, 2021) to measure the average tissue flow in the region posterior to the cauterization fence. The results are shown in Appendix 6—figure 1. In tor embryos, where no PMG invagination occurs, there is essentially no tissue flow, in line with the qualitative assessment from tissue ripping. Since the PMG is attached to the vitelline membrane in the scab domain (Münster et al., 2019), it acts as a barrier in tor embryos, preventing the germ band from elongating. Notably, in the eve/rnt and the twi/snl embryo, tissue flow velocity posterior of the cauterized region is reduced by about two-fold compared to wild-type. This shows that impairing active force generation in the germ band has a strong influence on the tissue flow anterior of the PMG. This flow originates from the dorso-ventral contraction of the germ band that pulls on the dorsal amnioserosa and drives the hyperbolic flow there (Streichan et al., 2018). These quantitative observations based on the data from Collinet et al., 2015 are in line with the much-reduced tissue flow in eve and twist embryos imaged in toto using light-sheet microscopy (Lefebvre et al., 2023; Gustafson et al., 2022).

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

Going forward, it would be interesting to image tor and eve embryos with a membrane-marker in toto and analyze them with the methods presented here. In toto imaging is crucial to resolve the contributions of different tissue regions as it allows to precisely time-align the recordings (Mitchell et al., 2022).

Appendix 6—figure 1

Download asset Open asset