# Marangoni-like tissue flows enhance symmetry breaking of embryonic organoids

Simon Gsell, Sham Tlili, Matthias Merkel, Pierre-François Lenne

► To cite this version:

Simon Gsell, Sham Tlili, Matthias Merkel, Pierre-François Lenne. Marangoni-like tissue flows enhance symmetry breaking of embryonic organoids. Nature Physics, 2025, 21 (4), pp.644-653. 10.1038/s41567-025-02802-2. hal-04251716v2

HAL Id: hal-04251716

https://hal.science/hal-04251716v2

Submitted on 23 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Marangoni-like tissue flows enhance symmetry breaking of embryonic organoids

Simon Gsell,¹,²,³,* Sham Tlili,¹,† Matthias Merkel,²,‡ and Pierre-François Lenne¹,§

¹Aix Marseille Univ, CNRS, IBDM (UMR 7288),

Turing Centre for Living systems, Marseille, France

²Aix Marseille Univ, Université de Toulon, CNRS, CPT (UMR 7332),

Turing Centre for Living systems, Marseille, France

³Aix Marseille Univ, CNRS, Centrale Méditerranée, IRPHE (UMR 7342),

Turing Centre for Living Systems, Marseille, France

(Dated: May 23, 2025)

During early development of multi-cellular animals, cells self-organize to set up the body axes, such as the primary head-to-tail axis. Several signalling pathways are known to control body axis formation. Here, we show that also tissue mechanics play an important role. We focus on the emergence of a primary axis in initially spherical aggregates of mouse embryonic stem cells, which mirrors events in the development of the early mouse embryo. These aggregates break rotational symmetry by establishing domains of different expression profiles, for example of the transcription factor T/Bra and the adhesion molecule E-cadherin. Combining quantitative microscopy and physical modeling, we identify large-scale tissue flows with a recirculating component, which contribute substantially to the symmetry breaking. We show that the recirculating flows are – akin to Marangoni flows – driven by a difference in tissue surface tensions, whose existence we further confirm using aggregate fusion experiments. Our work highlights that body axis formation is not only driven by biochemical processes, but that it can also be amplified by tissue flows. We expect that this type of amplification may operate in many other organoid and in-vivo systems.

Deciphering the mechanisms underlying the formation of the mammalian body plan is an important challenge in biology. This process takes place in the early embryo and is characterized by the establishment of a multicellular spatial organization with reference to an emerging system of orthogonal axes: the head-to-tail axis, also called the anteroposterior (AP) axis, and an orthogonal dorsoventral axis [1]. There is evidence that a small number of conserved morphogens, including Wnt, BMP, Activin/Nodal, and FGF, control the emergence of the body axes [2]. These signals participate in tissue patterning and trigger differential behaviors, which ultimately partition an isotropic group of cells into different territories. Yet, the physical mechanisms underlying axis formation remain largely unknown.

Here we focus on the emergence of the head-to-tail body axis, which breaks the initial rotational symmetry. In vertebrates, the AP axis establishes first by locating the transcription factor T/Brachyury (T/Bra) to one end of the embryo. Yet, understanding the emergence of this T/Bra pole is a daunting task because of the experimental difficulty in accessing the embryo, in particular in mammals [3–5]. To address this issue, an excellent experimental system is provided by 3D aggregates of mouse and human embryonic stem cells (ESCs), called gastruloids. Gastruloids are embryonic organoids that undergo symmetry breaking and exhibit axial organization and gene expression patterns that mirror events

in the embryo [6–10]. Depending on the culture conditions and the signaling molecules present in the medium, gastruloids have the capacity to develop advanced structures strikingly similar to organs such as (1) somites and neural-tube [8, 11–13] or (2) gut and heart [11, 14, 15]. Past works on gastruloids have identified the temporal sequence of gene expression [7] and networks of genes that regulate the main steps of their early development, including during symmetry breaking [16].

Yet, there is growing evidence that tissue mechanics and flows are crucial for many stages of embryonic development [17–21]. Indeed, in avians it has recently been pointed out that tissue mechanics and flows could also contribute to embryonic axis formation [19, 22, 23], but their potential role for axis formation in mammals remains unclear. Recently, using a mouse gastruloid protocol derived from the standard protocol in the field (Figure 1a, Methods) [24], we have shown that symmetry breaking in mouse gastruloids is accompanied by large-scale tissue flows [15]. However, what role these flows play in gastruloid polarization and how they are created is so far unclear.

Here we quantitatively show that in our mouse gastruloids, large-scale tissue flows contribute substantially to symmetry breaking. We show that the dominating component of the tissue flows is a recirculating flow whose direction correlates strongly with the T/Bra and E-cad polarization of the aggregate. We then present a minimal effective model for the emergence of these flows. In particular, we show that effective interfacial tensions between low and high T/Bra or E-cad domains, as well as differential surface tensions of these domains, are sufficient to explain the recirculating flows. Finally, we further validate our hypothesis of differential tensions by measuring

<footnote>* S. Gsell and S. Tlili contributed equally.; simon.gsell@univ-amu.fr</footnote>

<footnote>† S. Gsell and S. Tlili contributed equally.; sham.tlili@univ-amu.fr</footnote>

<footnote>‡ matthias.merkel@univ-amu.fr</footnote>

<footnote>§ pierre-francois.lenne@univ-amu.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1. Experimental workflow for quantifying symmetry breaking in gastruloids. (a) Symmetry breaking in gastruloids: mouse embryonic stem cells are seeded in low-adherence microwells. After 48 h, a 24 h Chiron pulse activates Wnt signaling, promoting mesodermal differentiation and inducing T/Brachyury expression (green). From 72 h to 96 h, graded T/Brachyury expression aligns with gastruloid elongation, marking the anterior (A, low T/Brachyury) and posterior (P, high T/Brachyury) poles (see Methods). (b) Inverted two-photon microscopy is used to image sub-volumes of the gastruloids around their midplane, with each orange plane representing a scan at a specific depth (see Methods). (c) On the left, time-averaged bidimensional velocity field of tissue collective flows during symmetry breaking in the gastruloid midplane, derived from live imaging data using optic flow. Vector orientations are color-coded, with red indicating flow towards the posterior pole and blue indicating flow towards the anterior pole. On the right, the corresponding T/Brachyury expression pattern is shown with inverted grayscale for clarity. (See Methods for details.) (d) E-cadherin and Snail expression patterns in an 80-hour gastruloid similar as the one in (c), obtained via immunochemistry (see Methods and Extended Data Table 1). Grayscale levels are inverted for enhanced visibility and the scale bar corresponds to 50 µm.</figcaption></figure>

ing interface angles in fusion experiments between gastruloids. Taken together, our work quantitatively shows that not only genetic and biochemical interactions, but also tension-driven large-scale tissue flows promote symmetry breaking in gastruloids.

## TISSUE FLOWS SUBSTANTIALLY CONTRIBUTE TO GASTRULOID SYMMETRY BREAKING

To precisely quantify the contribution of large-scale tissue flows to the overall gastruloid symmetry breaking, we need to simultaneously measure tissue flows and T/Bra or E-cad expression patterns. We generate the gastruloids following our previously published protocol (Figure 1a, Methods) [15], where we focus on the time between 72 and 96 hAAF, when gastruloids start to elongate and display polarized T/Bra and E-cad levels (Figure 1c,d, Methods). We live imaged the gastruloids dur-

ing this stage with live reporters for T/Bra, E-cad, cell membranes, and extra-cellular space (Methods). We reduced light scattering (gastruloid diameters of a few hundred µm) while preventing phototoxic effects by using two-photon microscopy and imaging 60 µm-thick subvolumes around their mid-planes (Figure 1b, Methods). To quantify flow fields, we first applied a registration removing any global aggregate translation and rotation (Methods), and then applied an optical flow method [25] to the resulting movie (Methods, Supplementary Figure 6-7, Supplementary movie 6).

To quantify the degree of symmetry breaking of a gastruloid, i.e. its polarization, we use the dipole moment $P(t)$ of its T/Bra or E-cad distribution (Figure 2a-c, Methods) given by the 2D vector:

$$P(t) = \int_{\Omega} \phi(\mathbf{r}, t) [\mathbf{r} - \mathbf{c}(t)] \, dS. \quad (1)$$

Here, $\Omega$ represents the 2D cross section of the aggregate, $\phi$ denotes the T/Bra or E-cad fluorescence intensity field,
$r$ is the position vector that we integrate over, and $c$ is the barycenter of the aggregate. The direction of the vector $P$ points towards the part of the gastruloid that has the highest fluorescence intensity in T/Bra or E-cad, and its magnitude quantifies the degree of polarization.

We found that after registration, the direction of the dipole moment does not change much over time (Supplementary Figure 8a). This indicates that very early processes already define an axis that is kept later on. However, we also see that the magnitude of $P$ increases over time, reflecting a progressive build up of gastruloid polarization (Figure 2b,c).

We wondered how important different cellular processes are for the progressive increase in gastruloid polarization over time. There are different cell-scale processes that might contribute, which include (i) spatially heterogeneous cell differentiation, (ii) individual cell motion, (iii) collective cell motion (i.e. tissue flows), and differential (iv) cell division and/or (v) apoptosis (Figure 2d). To differentiate the contribution of all these processes, segmentation and tracking of the cell nuclei would be required, which is currently difficult due to light scattering and dense nuclear packing. Thus, we base our analysis here on coarse-grained optical flow and fluorescence data, which still allows us to differentiate three different contributions to gastruloid polarization: (1) the combined effect of cell differentiation and individual cell motion (contributions i and ii), which we call for brevity "reaction-diffusion" (RD) contribution, (2) advection of T/Bra or E-cad with tissue flows (contribution iii), and (3) differential volume growth, arising for instance from cell division and apoptosis (contributions iv and v), which we call for brevity just "growth" contribution (Figure 2d).

Using the optical flow and fluorescence data, we can quantify the contributions 1-3 to the overall increase in polarization. To this end, we decompose the time derivative of $P$ into four parts (Figure 2e, Methods, and section 1 of the SI):

$$ \frac{dP}{dt} = R + A + G + Q_{3D}. \quad (2) $$

The first three terms on the right-hand-side of Eq. (2) respectively represent the contributions 1-3 (Figure 2e), and $Q_{3D}$ is a geometric contribution related to 3D motion of the aggregate (section 1 of the SI). Besides this geometric contribution $Q_{3D}$, our mathematical decomposition also confirms that there is no other contribution to gastruloid polarization than contributions 1-3.

For each gastruloid, we quantify each term in Eq. (2) over time and plot their respective cumulative contributions to the overall dipole moment (Figure 2e and Supplementary Figures 8c-1). Figure 2f shows the final overall contributions across all gastruloids. We find that reaction-diffusion processes generally contribute most to the overall polarization. However, our analysis also reveals a significant contribution of advection, which is responsible for approximately 1/3 of the total polarization.

Taken together, our analysis shows that advection with

large-scale tissue flows significantly contributes to gastruloid symmetry breaking. To better understand how, we next analyzed the structure of the flow pattern across different gastruloids.

## TISSUE FLOWS EXHIBIT A DOMINANT RECIRCULATING COMPONENT.

Qualitatively, tissue flows emerging during the polarization process exhibit substantial fluctuations, both spatio-temporally and across different gastruloids (Supplementary Figures 7 and 9-12). In order to capture the dominant large-scale flow structures and compare them across gastruloids, we decomposed the tissue flow field into different modes (Figure 3, Figure 4, Methods, and section 2 of the SI).

We apply our mode decomposition analysis on 14 gastruloids imaged between 78hAAF and 88hAAF. The analysis reveals that about 50% of the tissue flow field is provided by only 5 low-index modes, i.e. large-scale features, which are robustly observed across different gastruloids (Figure 4a-c). In particular, the three modes of highest amplitude together contribute approximately 40% to the total flow field. The superposition of these three modes corresponds to a large-scale recirculating pattern (Figure 4d). Further analysis using the full roughly 16-hour observation window revealed that while the recirculation flows corresponded to roughly 5-10% of a full turn, individual cells moved substantially further during the same time (Supplementary figure 15).

The next most relevant mode identified in Figure 4c, $v_{2,0}^r(r)$ with amplitude $A_{2,0}^r$ (Figure 4a right), contributes to gastruloid elongation. It corresponds to a flow field that extends along one axis and contracts along the perpendicular axis, consistent with gastruloid elongation between 72 and 96hAAF.

We find that the direction of the large-scale recirculating pattern strongly correlates with the T/Bra polarization, where the direction of the flow in the aggregate center points towards high T/Bra concentrations (Figure 4e). This correlation between the directions of T/Bra polarization and recirculating flow could be created by the advection of T/Bra with the tissue flows. However, could it also be possible that, conversely, aggregate polarization drives the recirculating flows?

## THE RECIRCULATING FLOWS COULD BE DRIVEN BY PATTERNING-INDUCED TENSION.

To understand how the aggregate polarization – in terms of gene expression patterns – could affect tissue mechanics and flows, we first noted that from a coarse-grained perspective, the expression patterns of many genes are qualitatively similar in showing a distinct AP gradient (Figure 1c,d and Supplementary Figures 1-5).
<figure><img src="image_3.png" /><figcaption>Figure 2. Quantitative analysis of the polarization process. (a) Example of T/Bra distribution and corresponding coarse-grained field $\phi$. We define the dipole moment by Eq. (1), where $\Omega$ represents the 2D cross section of the aggregate, $r$ is the position vector that we integrate over and $c$ is the barycenter of the aggregate. (b) Example of the time evolution of the T/Bra dipole moment and corresponding T/Bra distributions. We normalize the dipole moment $P$ by a reference value $P_0(t)$ that depends on the spatial standard deviation of $\phi$ and the total cross-sectional area of the aggregate (Methods and section 1 of the SI). (c) Time evolution of T/Bra and E-cad dipole moments for different samples. Each sample is indicated with a distinct color. (d) At the cellular level, the polarization process relies on the interplay between cell differentiation, cell rearrangements and cell division/apoptosis, resulting in three effective mechanisms at the coarse-grained level, namely reaction-diffusion, advection and growth. (e) Example of time evolution of the cumulative contributions to T/Bra polarization. For instance, $\bar{P}(t) := \int_{t_i}^t (dP(t')/dt') dt' = P(t) - P(t_i)$, $\bar{R} := \int_{t_i}^t R(t') dt'$, etc (SI). The grey line indicates geometric contributions to polarization related to 3D motion of the aggregate (section 1 of the SI). (f) Relative total contribution of each process at the end of the polarization. In (f), boxes extend from the first quartile to the third one, the horizontal line indicates the median value, and whiskers end at the most distant point whose value lies within 1.5 times the interquartile range.</figcaption></figure>

We thus studied a minimal model that represents these AP patterning genes by a single effective expression pattern $\phi(r)$, which ranges from -1 (anterior-most) to +1 (posterior-most). The field $\phi(r)$ reflects the combined effect of the AP patterning genes on tissue flows. In particular, the AP patterning genes include differentially expressed cadherins such as E-cad [15], which are known to lead to tension-driven cell separation [26–30].

To understand whether the observed gastruloid patterning can drive the recirculating flow, we study a sim-

ple continuum model that describes the gastruloid as an incompressible viscous fluid whose flows are driven by gradients of the field $\phi(r)$. Formally, this is described by Stokes’ equations:

$$0 = \nabla \cdot \mathbf{v} \qquad (3a)$$

$$0 = -\nabla \Pi + \eta \nabla^2 \mathbf{v} + \nabla \cdot (\boldsymbol{\sigma}_i + \boldsymbol{\sigma}_s), \qquad (3b)$$

where $\Pi$ is the hydrodynamic pressure ensuring incompressibility, $\eta$ is the effective tissue viscosity, and $\boldsymbol{\sigma}_i$ and $\boldsymbol{\sigma}_s$ are symmetric, traceless tensors that denote addi-
<figure><img src="image_4.png" /><figcaption>Figure 3. Mode decomposition of the tissue flow. (a) Coarse-grained tissue flows are analyzed in a circular disk (red dashed line), centered on the aggregate barycenter c. (b) Tissue velocity v is decomposed as the sum of a radial component v_r and an angular component v_θ. Similarly, positions within the aggregate can be expressed in polar coordinates (r, θ). (c) The velocity components v_r(r) and v_θ(r) are decomposed into a series of angular (index n) and radial (index p) modes (Methods and section 2 of the SI). These modes are plotted for n ≤ 4 and p ≤ 4. (d) The first few modes v_{n,p}^{r/θ} are illustrated by their spectrum, the corresponding angular and radial functions, and the resulting radial and angular velocity fields.</figcaption></figure>

tional mechanical stresses driving tissue flows. The stress $\sigma_i$ is internally created by the effective expression pattern $\phi(r)$, and $\sigma_s$ represents the tension at the aggregate surface.

There are different ways in which the internal stress $\sigma_i$ could depend on $\phi(r)$. However, as pointed out previously [31], the leading-order effect that drives incompressible flow is given by

$$ \sigma_i = -\kappa_i \left[ (\nabla\phi) \circ (\nabla\phi) - \frac{1}{2}(\nabla\phi)^2 I \right]. \quad (4) $$

Here, the coefficient $\kappa_i$ controls the magnitude of the internal mechanical stress $\sigma_i$, the symbol $\circ$ denotes the dyadic product, and $I$ denotes the identity tensor. The internal stress tensor $\sigma_i$ is locally aligned with the gradient of $\phi(r)$, and its magnitude increases monotonically with the magnitude of the gradient (Figure 5a). At an interface between regions of low and high $\phi(r)$, the tensor $\sigma_i$ determines an effective interface tension $\gamma_i$, which is obtained by integrating $\sigma_i$ across the interface width (Figure 5a,b).

To describe aggregate surface tension $\sigma_s$, we introduce a second field $\psi(r)$, which is equal to 1 inside the aggregate and 0 outside (Figure 5c). The stress $\sigma_s$ can be computed from gradients of $\psi(r)$ analogously to Eq. (4), but with a magnitude $\kappa_s$ that depends on local cell fate, i.e. on the $\phi(r)$ field. In the two limiting cases of $\phi = -1$ (anterior regions) and $\phi = +1$ (posterior regions), we

denote the surface stress magnitude by $\kappa_s^-$ and $\kappa_s^+$, respectively (Figure 5c, details in section 3 of the SI).

To test how far our simple effective model can account for the observed tissue flows, we solve Eqs. (3) and (4) using a lattice Boltzmann approach [32, 33] (section 3 of the SI), where we use for $\phi(r)$ the coarse-grained experimentally determined T/Bra or E-cad fields, as representative AP expression patterns. In Figure 5d, we plot for one of the T/Bra fields the resulting tissue flow field varying two dimensionless stress ratios, a relative surface stress difference, $\Delta\kappa_s := (\kappa_s^+ - \kappa_s^-)/\kappa_i$ and an interface-to-surface stress ratio, $r := 2\kappa_i/(\kappa_s^+ + \kappa_s^-)$. We find that the surface tension difference $\Delta\kappa_s$ indeed creates a recirculating flow pattern. This is because the surface tension gradient drives surface flows, a mechanism known as Marangoni effect [34]. For $\Delta\kappa_s < 0$ the flow on the aggregate surface is oriented towards the anterior pole, similarly to what we observe in our experiments, while for $\Delta\kappa_s > 0$ it is oriented towards the posterior pole (Figure 5d). This qualitative comparison suggests that in our gastruloids the surface tension of anterior regions may be higher than that of posterior regions. This conclusion is further corroborated by an analysis of the residuals between experimentally measured and predicted flow fields depending on the simulation parameters (Figure 5f-h, Methods).

Our results show that at the coarse-grained level, both the T/Bra and E-cad fields result in similar computational flow patterns (see also Figures S9, S10), as only
<figure><img src="image_5.png" /><figcaption>Figure 4. Analysis of tissue flow during polarization. (a) Example of a tissue velocity field expressed as a sum of the base modes $v_{n,p}^{r/\theta}$ shown in Figure 3d, where the modes are weighted by the amplitudes $A_{n,p}^r$ and $A_{n,p}^\theta$. Below each mode we indicate the approximate contribution to the total variance of the flow. (b) Spectra displaying the contributions of each mode to the overall flow (compare Figure 3d, leftmost column). Same sample as in (a). (c) Leading mode contributions over different samples. (d) We reconstruct the dominant recirculation mode over samples using median amplitudes of the three leading modes $\overline{A}_{1,1}^\theta$, $\overline{A}_{1,0}^\theta$ and $\overline{A}_{1,0}^r$. (e) Scatter plot showing the orientation angle of the recirculation mode versus the orientation angle of the dipole moment for each sample. The black line shows the linear fit and PCC indicates the Pearson correlation coefficient. In (c), boxes extend from the first quartile to the third one, the horizontal line indicates the median value, and whiskers end at the most distant point whose value lies within 1.5 times the interquartile range.</figcaption></figure>

the overall AP pattern is required to drive the main re-circulation flow. Taken together, the comparison between our model and experimental data suggests: (i) a higher surface tension in anterior (T/Bra−, E-cad−) than in posterior (T/Bra+, E-cad+) tissues and (ii) an interface tension between these tissues. To further test this, we performed fusion experiments between gastruloids.

## FUSION EXPERIMENTS CONFIRM TISSUE-TYPE-DEPENDENT TENSIONS.

To assess the existence of tissue-dependent tensions in the gastruloids, we performed fusion experiments to force posterior-like (T/Bra+) tissues to mechanically interact with anterior-like (T/Bra−) tissues (Figure 6, Supplementary Movie 7). To this end, we used our observations that gastruloids at 72hAAF display high T/Bra levels almost everywhere, while 84hAAF gastruloids are already polarized, exhibiting both T/Bra+ and T/Bra− regions.

We thus brought 72hAAF gastruloids in contact with 84hAAF gastruloids.

Different fusion scenarios occurred depending on the relative position of both aggregates (Figure 6, Supplementary Movie 8). We first focus on interactions between T/Bra+ and T/Bra− tissues by selecting experiments with only a heterotypical +− contact of the 72hAAF gastruloid with a T/Bra− region of the 84hAAF gastruloid (Figure 6a). As a control experiment, we also consider homotypical ++ fusion experiments between two 72hAAF T/Bra+ aggregates (Supplementary Movie 9). We track the fusion angles $\theta_1$, $\theta_2$ and $\theta_3$ over time during the fusion process (Figure 6a,c,d, Methods). In the case of homotypical ++ fusions, the outer fusion angle $\theta_1$ saturates to 180° (red solid curve in Figure 6c), while the inner angles $\theta_2$ and $\theta_3$ both saturate at approximately 90° (red and orange solid curves in Figure 6d). This behavior is expected and quantitatively reproduced by our model when we set $\kappa_s^− = \kappa_s^+$ and $r = 0$ (Figure 6b, simulation 1; red dashed curves in Figure 6c,d).
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

Figure 5. Computational modeling of the recirculating tissue flows. (a) Gradients of a scalar field $\phi$ control stresses (black crossed arrows), which results in an effective tension $\gamma_i$ (black double arrows) at the interface between domains of high and low $\phi$. The interface stress $\sigma_i$ is controlled by the stress coefficient $\kappa_i$. $I$ denotes the identity matrix. (b) We construct coarse-grained scalar fields from the experimentally observed T/Bra or Ecad fields and we compute the resulting active stresses within the aggregate. (c) We describe the aggregate surface as a smooth interface using a second scalar field $\psi$. The effective tension at the aggregate surface is controlled by a $\phi$-dependent stress coefficient, with two limit values $\kappa_s^-$ and $\kappa_s^+$ corresponding to surfaces of T/Bra- (or E-cad-) and T/Bra+ (or E-cad+) tissues, respectively. (d) We compute tissue flows driven by the T/Bra- or E-cad-dependent stresses over ranges of the interface ($\kappa_i$) and surface ($\kappa_s^+$ and $\kappa_s^-$) stress coefficients, varied through the non-dimensional parameters $r$ and $\Delta\kappa_s$. The simulated flow fields shown here were computed from the T/Bra field plotted in panel b. (e) Experimental velocity field corresponding to panels b and d. (f) Residual between simulated and experimentally observed flow field, averaged over all gastruloids, as a function of the model parameters $\Delta\kappa_s$ and $r$. (left) Residual between the full flow fields, (right) residual between the modes contributing to the recirculation flow only. (g) Dominant modes $v_{n,p}^{r/\theta}$ (Figure 3, Figure 4 and Methods) contributing to the simulated flow field for two parameter sets, respectively minimizing the total flow residual ($r = 2.5, \Delta\kappa_s = -1$) and the recirculating flow residual ($r = 0.9, \Delta\kappa_s = -0.2$). (h) Reconstructed recirculation flow in experiments and simulations. In (g), boxes extend from the first quartile to the third one, the horizontal line indicates the median value, and whiskers end at the most distant point whose value lies within 1.5 times the interquartile range.
<figure><img src="image_13.png" /></figure>

Figure 6. Quantitative analysis of aggregate fusion experiments. (a) Fusion experiments with either a homotypical ++ contact between two 72hAAF gastruloids or a heterotypical +− contact between a 72hAAF gastruloid and the T/Bra− region of a 84hAAF gastruloid. We define the fusion angles $θ_1$, $θ_2$ and $θ_3$ as indicated. $t_f$ corresponds to the end of the fusion experiment, which lasts 600 minutes. (b) Experiments are compared to model simulations of fusing droplets. we consider two particular cases: simulation 1 is performed in the absence of tension between the droplets, while simulation 2 is performed in the presence of interface tension and surface tension difference, namely $r = 0.8$ and $\Delta\kappa_s = -0.2$. (c,d) Time evolution of the (c) outer angle $θ_1$ and (d) inner angles $θ_2$ and $θ_3$, averaged over samples and compared to the angles in the respective simulation. The shaded areas indicate the 95% confidence interval. (e) Example of hybrid fusion, where a 72hAAF gastruloid is in contact with both a T/Bra+ and a T/Bra− region of a 84hAAF gastruloid. (f) Time evolution of the fusion interface length, decomposed as the sum of homotypic contact length between two T/Bra+ tissue regions (++), and heterotypic contact length between a T/Bra+ and a T/Bra− region (−−).

For heterotypical +− fusions, the outer fusion angle saturates to a lower value of $θ_1 \approx 140°$ (blue solid curve in Figure 6c). Moreover, the inner angle $θ_2$ in the T/Bra+ tissue saturates to a lower value than the inner angle $θ_3$ in the T/Bra− tissue (blue and green solid curves, respectively, in Figure 6d). This is consistent with our model, where simulations performed with $\Delta\kappa_s = -0.2$ and $r = 0.8$ quantitatively reproduce the observed fusion angles (Figure 6b, simulation 2; blue and green dashed curves in Figure 6c,d). Taken together, our fusion observations are consistent with our previous conclusions that (i) T/Bra− tissues have a higher surface tension than T/Bra+ tissues ($θ_2 < θ_3$), and (ii) there is a finite interface tension between T/Bra+ and T/Bra− tissues ($θ_1 < 180°$).

The existence of finite interface tension between T/Bra+ and T/Bra− tissues is further supported by hybrid fusions, when a 72hAAF gastruloid was in contact with both T/Bra+ and T/Bra− regions of a 84hAAF gastruloid (Figure 6e, Supplementary Movie 8). In this case, the homotypic (++) interface extends in length while the heterotypic (+−) interface slightly shrinks (Figure 6f). This suggests that the affinity between two T/Bra+ regions is larger than that between T/Bra+ and T/Bra− regions, i.e. there is an effective interface tension between T/Bra+ and T/Bra− regions.

While in the previous section the velocity amplitude was essentially a fit parameter (section 3 of the SI),
the fusion experiments allow us to independently predict the order of magnitude of the expected recirculating flow speed. The time scale on which two droplets with viscosity $\eta$, radius $R$, and surface tension $\gamma_s$ fuse is given by $t_F \sim R\eta/\gamma_s$ [35, 36]. According to Figure 6e, we estimate $t_F \sim 1$ h. Moreover, in a droplet whose surface tension varies between $\gamma_s^-$ and $\gamma_s^+$, recirculating Marangoni flows traverse the size of the droplet on the time scale $t_M \sim R\eta/|\gamma_s^+ - \gamma_s^-| \simeq t_F/(r|\Delta\kappa_s|)$ [37]. Using $\Delta\kappa_s = -0.2$ and $r = 0.8$, we thus obtain $t_M \sim 6$ h, which is close to the time scale of gastruloid symmetry breaking $t_{SB} \sim 10$ h. This confirms that the magnitude of the differential tension observed in our fusion experiments is sufficient to drive the observed recirculation flows.

Taken together, our results suggest that large-scale recirculating flows are created by a small, early gastruloid polarization. In a positive feedback loop, these flows then further amplify and thus stabilize gastruloid polarization.

## DISCUSSION

Since the mid-20th century, it has been proposed that effective tissue surface tensions can drive the separation of cell population [26, 28, 29, 38–41], often as a consequence of differential cadherin expression between cells [29, 30]. Indeed, in our gastruloids, cells progressively lose E-cadherin and show spatial anti-correlation between E-cadherin and N-cadherin/OB-cadherin [15] (Supplementary Figure 16). However, past research on differential tissue tensions driving cell population separation considered only a diffusive unmixing process, including in stem cell aggregates [42, 43]. In contrast, our work demonstrates that differential tissue surface tensions can generate long-range hydrodynamic flows sustained over extended periods (>$\sim$ 10 hours), which have been shown to significantly accelerate cellular unmixing [44].

There are several possible future extensions of this work. First, most past work on cellular unmixing has focused on differentiated cells [45]. Yet, if cells keep their initial differentiation states, flows driven by tissue tensions may eventually stop [44]. In contrast, our gastruloids maintain recirculating flows over several hours, raising the question of whether these flows are sustained by patterned cell differentiation. Indeed, in previous work, we have shown that T/Bra-negative cells in the posterior pole are pluripotent, differentiating into T/Bra-positive cells, while T/Bra-positive and E-Cad-negative cells evolve into T/Bra-negative and E-Cad-negative cells [15]. Future research will explore how the existence and maintenance of the pluripotent cell population in the posterior pole may sustain the observed flows.

Second, here we studied a 2D section of the 3D flows in the aggregates, which was due to technical limitations in gastruloid *in toto* imaging, including photo-toxicity and line scanning speed constraints. Consequently, we had to filter out aggregates where the inferred out-of-plane polarization component was too high (section 1 of the SI,

Supplementary movie 10). It remains a future challenge to overcome light scattering for *in toto* live imaging of the aggregates in order to obtain the full 3D expression patterns and flows [46]. Similarly, in future work, it will be important to further differentiate the reaction-diffusion term of the polarization into contributions by gene expression dynamics and small-scale cellular motion. Yet, this will require tracking individual cells or nuclei, which is still challenging in gastruloids.

Third, it will be interesting to examine the cellular mechanisms underlying differences in effective tissue surface tension [28, 41, 47], considering factors like adhesion, cell contractility, and motility (Supplementary figures 16-17). These contributions are not mutually exclusive, and determining their exact roles will require generating specific knock-out mutant cell lines or optogenetic constructs, which is beyond the scope of this work. Modifying the tissue surface tension using these tools would also provide an independent test for whether and how gastruloids polarize with reduced large-scale recirculating flows.

Fourth, while we focused on recirculating flows in the aggregate, we did not address the mechanisms behind aggregate elongation, which require active forces and aligned orientational information [48]. Research on fish embryonic stem cell aggregates shows that cell polarity proteins are crucial for elongation [49, 50]. Heterogeneous mechanical properties could also play a role [51]. Identifying the precise mechanisms driving gastruloid elongation will be important for future work.

We expect Marangoni-effect-driven recirculating tissue flows to operate in various organoid systems formed from pluripotent cells. These flows can be induced artificially [52] or driven internally, arising from distinct cell populations with different surface tensions, amplifying symmetry breaking. They may also be relevant in the gastrulating mouse embryo, where epithelial cells undergo EMT and delaminate from the epiblast, migrating as mesenchymal cells [2, 53]. Recirculating flows have also been observed during the formation of the cardiac crescent in gastrulation [3]. Thus, the feedback loop seen in our mouse embryonic stem cell aggregates may be significant across various *in vitro* and *in vivo* systems.

## ACKNOWLEDGMENTS

M.M. thanks the Centre Interdisciplinaire de Nanoscience de Marseille (CINaM) for providing office space. This work is supported by the French National Research Agency (“France 2030”, ANR-16-CONV-0001 from Excellence Initiative of Aix-Marseille University - A*MIDEX and generic grant to P.-F.L. ANR-19-CE13-0022), the Fondation de la Recherche Médicale (to P.-F.L. EQU202003010407), a generic grant (to S.T. ANR-22-CE30-0021) and an ERC grant (to M.M. and P.-F.L., ERC SyG 101072123). We also acknowledge the France-Bioimaging Infrastructure
(ANR-10-INBS-04). The Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources. S.T. and P.F.L. received financial support from Inserm to the Booster Program Mecacell3D. S.T. thanks Daniel Sapede for his invaluable help on implementing multiposition timelapse.

## AUTHOR CONTRIBUTIONS STATEMENT

S.G., S.T., M.M. and P.-F.L. designed the study, interpreted the results, and wrote the paper. S.T. performed the experiments, the imaging and the image and optical flow analysis. S.G. performed the data analysis and simulations.

## COMPETING INTERESTS STATEMENT

The authors declare no competing interest.

## REFERENCES

[1] L. Solnica-Krezel and D. S. Sepich, Gastrulation: Making and Shaping Germ Layers, Annual Review of Cell and Developmental Biology 28, 687-717 (2012).

[2] S. J. Arnold and E. J. Robertson, Making a commitment: cell lineage allocation and axis patterning in the early mouse embryo., Nature reviews. Molecular cell biology 10, 91-103 (2009).

[3] M. H. Dominguez, A. L. Krup, J. M. Muncie, and B. G. Bruneau, Graded mesoderm assembly governs cell fate and morphogenesis of the early mammalian heart, Cell 186, 479-496 (2023).

[4] K. McDole, L. Guignard, F. Amat, A. Berger, G. Marylandin, L. A. Royer, S. C. Turaga, K. Branson, and P. J. Keller, In toto imaging and reconstruction of post-implantation mouse development at the single-cell level, Cell 175, 859-876 (2018).

[5] B. Saykali, N. Mathiah, W. Nahaboo, M.-L. Racu, L. Hammou, M. Defrance, and I. Migeotte, Distinct mesoderm migration phenotypes in extra-embryonic and embryonic regions of the early mouse embryo, eLife 8, e42434 (2019).

[6] S. C. van den Brink, P. Baillie-Johnson, T. Balayo, A.-K. Hadjantonakis, S. Nowotschin, D. A. Turner, and A. Martinez Arias, Symmetry breaking, germ layer specification and axial organisation in aggregates of mouse embryonic stem cells, Development 141, 4231-4242 (2014).

[7] L. Beccari, N. Moris, M. Girgin, D. A. Turner, P. Baillie-Johnson, A.-C. Cossy, M. P. Lutolf, D. Duboule, and A. M. Arias, Multi-axial self-organization properties of mouse embryonic stem cells into gastruloids, Nature 562, 272-276 (2018).

[8] S. C. van den Brink, A. Alemany, V. van Batenburg, N. Moris, M. Blotenburg, J. Vivié, P. Baillie-Johnson, J. Nichols, K. F. Sonnen, A. Martinez Arias, and A. van Oudenaarden, Single-cell and spatial transcriptomics reveal somitogenesis in gastruloids, Nature 582, 405-409 (2020).

[9] N. Moris, K. Anlas, S. C. van den Brink, A. Alemany, J. Schröder, S. Ghimire, T. Balayo, A. van Oudenaarden, and A. Martinez Arias, An in vitro model of early anteroposterior organization during human development,

Nature 582, 410-415 (2020).

[10] N. Moris, A. Martinez Arias, and B. Steventon, Experimental embryology of gastrulation: Pluripotent stem cells as a new model system, Current Opinion in Genetics & Development Cell Reprogramming, Regeneration and Repair, 64, 78-83 (2020).

[11] S. C. van den Brink and A. van Oudenaarden, 3d gastruloids: A novel frontier in stem cell-based in vitro modeling of mammalian gastrulation, Trends in Cell Biology 31, 747-759 (2021).

[12] J. V. Veenvliet, A. Bolondi, H. Kretzmer, L. Haut, M. Scholze-Wittler, D. Schifferl, F. Koch, L. Guignard, A. S. Kumar, M. Pustet, S. Heimann, R. Buschow, L. Wittler, B. Timmermann, A. Meissner, and B. G. Herrmann, Mouse embryonic stem cells self-organize into trunk-like structures with neural tube and somites, Science 370, eaba4937 (2020).

[13] M. Matsuda, Y. Yamanaka, M. Uemura, M. Osawa, M. K. Saito, A. Nagahashi, M. Nishio, L. Guo, S. Ikegawa, S. Sakurai, S. Kihara, T. L. Maurissen, M. Nakamura, T. Matsumoto, H. Yoshitomi, M. Ikeya, N. Kawakami, T. Yamamoto, K. Woltjen, M. Ebisuya, J. Toguchida, and C. Alev, Recapitulating the human segmentation clock with pluripotent stem cells, Nature 580, 124-129 (2020).

[14] G. Rossi, N. Broguiere, M. Miyamoto, A. Boni, R. Guiet, M. Girgin, R. G. Kelly, C. Kwon, and M. P. Lutolf, Capturing cardiogenesis in gastruloids, Cell stem cell 28, 230-240 (2021).

[15] A. Hashmi, S. Tlili, P. Perrin, M. Lowndes, H. Peradziryi, J. M. Brickman, A. Martínez Arias, and P.-F. Lenne, Cell-state transitions and collective cell movement generate an endoderm-like region in gastruloids, eLife 11, e59371 (2022).

[16] S. Suppinger, M. Zinner, N. Aizarani, I. Lukonin, R. Ortiz, C. Azzi, M. B. Stadler, S. Vianello, G. Palla, H. Kohler, A. Mayran, M. P. Lutolf, and P. Liberali, Multimodal characterization of murine gastruloid development, Cell Stem Cell 30, 867-884 (2023).

[17] R. Etournay, M. Popović, M. Merkel, A. Nandi, C. Blasse, B. Aigouy, H. Brandl, G. Myers, G. Salbreux,
F. Jülicher, and S. Eaton, Interplay of cell dynamics and
epithelial tension during morphogenesis of the Drosophila
pupal wing, eLife 4, e07090 (2015).

[18] E. Rozbicki, M. Chuai, A. I. Karjalainen, F. Song, H. M. Sang, R. Martin, H. J. Knölker, M. P. Macdonald, and C. J. Weijer, Myosin-II-mediated cell shape changes and cell intercalation contribute to primitive streak formation, *Nature Cell Biology* **17**, 397–408 (2015).

[19] M. Saadaoui, D. Rocancourt, J. Roussel, F. Corson, and J. Gros, A tensile ring drives tissue flows to shape the gas-trulating amniote embryo, *Science* **367**, 453–458 (2020).

[20] E. W. Gehrels, B. Chakrabortty, M.-E. Perrin, M. Merkel, and T. Lecuit, Curvature gradient drives po-larized tissue flow in the Drosophila embryo, *Proceedings of the National Academy of Sciences* **120**, e2214205120 (2023).

[21] S. Tlili, J. Yin, J.-F. Rupprecht, M. Mendieta-Serrano, G. Weissbart, N. Verma, X. Teng, Y. Toyama, J. Prost, and T. Saunders, Shaping the zebrafish myotome by in-tertissue friction and active stress, *Proceedings of the Na-tional Academy of Sciences* **116**, 25430–25439 (2019).

[22] O. Voiculescu, F. Bertocchini, L. Wolpert, R. E. Keller, and C. D. Stern, The amniote primitive streak is defined by epithelial cell intercalation before gastrulation, *Nature* **449**, 1049–1052 (2007).

[23] P. Caldarelli, A. Chamolly, A. Villedieu, O. Alegria-Prévot, C. Phan, J. Gros, and F. Corson, Self-organized tissue mechanics underlie embryonic regulation, *Nature* **633**, 887–894 (2024).

[24] D. A. Turner, M. Girgin, L. Alonso-Crisostomo, V. Trivedi, P. Baillie-Johnson, C. R. Glodowski, P. C. Hayward, J. Collignon, C. Gustavsen, P. Serup, B. Steventon, M. P. Lutolf, and A. M. Arias, Anteroposterior polarity and elongation in the absence of extra-embryonic tissues and of spatially localised signalling in gastruloids: Mammalian embryonic organoids, *Development* **144**, 3894–3906 (2017).

[25] B. D. Lucas and T. Kanade, An iterative image registration technique with an application to stereo vision, in *Proceedings of the 7th international joint conference on Artificial intelligence - Volume 2*, IJCAI'81 (Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1981) pp. 674–679.

[26] M. S. Steinberg, Reconstruction of Tissues by Dissociated Cells, *Science* **141**, 401–408 (1963).

[27] R. A. Foty and M. S. Steinberg, The differential adhesion hypothesis: A direct evaluation, *Developmental Biology* **278**, 255–263 (2005).

[28] J. D. Amack and M. L. Manning, Knowing the Bound-aries: Extending the differential adhesion hypothesis, *Science* **338**, 212–215 (2012).

[29] J.-L. Maitre, H. Berthoumieux, S. F. G. Krens, G. Salbreux, F. Julicher, E. Paluch, and C.-P. Heisenberg, Adhesion Functions in Cell Sorting by Mechanically Coupling the Cortices of Adhering Cells, *Science* **338**, 253–256 (2012).

[30] E. Cachat, W. Liu, K. C. Martin, X. Yuan, H. Yin, P. Hohenstein, and J. A. Davies, 2- and 3-dimensional synthetic large-scale of novo patterning by mammalian cells through phase separation, *Sci Rep* **6**, 20664 (2016).

[31] A. Tiribocchi, R. Wittkowski, D. Marenduzzo, and M. E. Cates, Active Model H: Scalar Active Matter in a Momentum-Conserving Fluid, *Physical Review Letters* **115**, 188302 (2015).

[32] S. Gsell, U. D'Ortona, and J. Favier, Multigrid dual-time-stepping lattice boltzmann method, *Physical Review E* **101**, 023309 (2020).

[33] S. Gsell, U. D'Ortona, and J. Favier, Lattice-Boltzmann simulation of creeping generalized newtonian flows: theory and guidelines, *Journal of Computational Physics* **429**, 109943 (2021).

[34] M. G. Velarde, R. K. Zeytounian, M. G. Velarde, M. Sayir, W. Schneider, B. Schrefler, G. Bianchi, and C. Tasso, eds., *Interfacial Phenomena and the Marangoni Effect*, CISM International Centre for Mechanical Sciences, Vol. 428 (Springer, Vienna, 2002).

[35] E. Flenner, L. Janosi, B. Barz, A. Neagu, G. Forgacs, and I. Kosztin, Kinetic monte carlo and cellular particle dynamics simulations of multicellular systems, *Physical Review E* **85**, 031907 (2012).

[36] D. Oriola, M. Marin-Riera, K. Anlaş, N. Gritti, M. Sanaki-Matsumiya, G. Aalderink, M. Ebisuya, J. Sharpe, and V. Trivedi, Arrested coalescence of multicellular aggregates, *Soft Matter* **18**, 3771–3780 (2022), arxiv:2012.01455.

[37] M. Schmitt and H. Stark, Marangoni flow at droplet interfaces: Three-dimensional solution and applications, *Physics of Fluids* **28**, 012106 (2016).

[38] J. Holtfreter, Gewebeaffinität, ein mittel der embryonalen formbildung, Arch. Zellforsch. **23**, 169–209 (1939).

[39] P. L. Townes and J. Holtfreter, Directed movements and selective adhesion of embryonic amphibian cells, *Journal of Experimental Zoology* **128**, 53–120 (1955).

[40] A. K. Harris, Is Cell sorting caused by differences in the work of intercellular adhesion? A critique of the Steinberg hypothesis., *Journal of theoretical biology* **61**, 267–285 (1976).

[41] M. L. Manning, R. A. Foty, M. S. Steinberg, and E.-M. Schoetz, Coaction of intercellular adhesion and cortical tension specifies tissue surface tension, *Proceedings of the National Academy of Sciences* **107**, 12517–12522 (2010).

[42] M. Bao, J. Cornwall-Scoones, E. Sanchez-Vasquez, A. L. Cox, D.-Y. Chen, J. De Jonghe, S. Shadkhoo, F. Hollfelder, M. Thomson, D. M. Glover, and M. Zernicka-Goetz, Stem cell-derived synthetic embryos self-assemble by exploiting cadherin codes and cortical tension, *Nat Cell Biol* **24**, 1341–1349 (2022).

[43] M. A. de Jong, E. Adegeest, N. M. L. P. Bérenger-Currias, M. Mircea, R. M. H. Merks, and S. Semrau, A combination of convergent extension and differential adhesion explains the shapes of elongating gastruloids, *bioRxiv* 10.1101/2023.05.24.541949 (2023).

[44] S. Gsell and M. Merkel, Phase separation dynamics in deformable droplets, *Soft Matter* **18**, 2672–2683 (2022).

[45] E.-M. Schötz, R. D. Burdine, F. Jülicher, M. S. Steinberg, C.-P. Heisenberg, and R. A. Foty, Quantitative differences in tissue surface tension influence zebrafish germ layer positioning., *HFSP J* **2**, 42–56 (2008).

[46] A. Gros, J. Vanaret, V. Dunsing-Eichenauer, A. Ros-tan, P. Roudot, P.-F. Lenne, L. Guignard, and S. Tlili, A quantitative pipeline for whole-mount deep imaging and multiscale analysis of gastruloids, *bioRxiv* 10.1101/2024.08.13.607832 (2024).

[47] T. V. Stirbat, A. Mgharbel, S. Bodennec, K. Ferri, H. C. Mertani, J.-P. Rieu, and H. Delanoë-Ayari, Fine tuning of tissues’ viscosity and surface tension through contractility suggests a new role for α-catenin, *PloS one* **8**, e52554 (2013).
[48] M. Ibrahimi and M. Merkel, Deforming polar active matter in a scalar field gradient, New Journal of Physics 25, 013022 (2023).

[49] T. Fulton, V. Trivedi, A. Attardi, K. Anlas, C. Dingare, A. M. Arias, and B. Steventon, Axis Specification in Zebrafish Is Robust to Cell Mixing and Reveals a Regulation of Pattern Formation by Morphogenesis, Current Biology 30, 2984-2994 (2020).

[50] M. L. Williams and L. Solnica-Krezel, Nodal and planar cell polarity signaling cooperate to regulate zebrafish convergence and extension gastrulation movements, eLife 9, e54445 (2020).

[51] A. Mongera, P. Rowghanian, H. J. Gustafson, E. Shelton, D. A. Kealhofer, E. K. Carn, F. Serwane, A. A. Lucio, J. Giammona, and O. Campàs, A fluid-to-solid jamming transition underlies vertebrate body axis elongation, Nature 561, 401-405 (2018).

[52] V. Yadav, M. S. Yousafzai, S. Amiri, R. W. Style, E. R. Dufresne, and M. Murrell, Gradients in solid surface tension drive marangoni-like motions in cell aggregates, Phys. Rev. Fluids 7, L031101 (2022).

[53] A. E. Sutherland, Tissue morphodynamics shaping the early mouse embryo, in *Seminars in cell & developmental biology*, Vol. 55 (Elsevier, 2016) pp. 89-98.

[54] E. J. Underhill and J. E. Toettcher, Control of gastruloid patterning and morphogenesis by the Erk and Akt signaling pathways, Development 150, dev201663 (2023).

## METHODS

**Cell culture and gastruloid formation.** In our protocol derived from the standard protocol in the field [24], we first seed a controlled number of pluripotent mouse ESCs into non-adherent micro-wells (Figure 1a, Methods). After a few hours, the cells form 3D aggregates. 48 hours after aggregate formation (hAAF), differentiation is stimulated by adding the Wnt-agonist Chiron, which is washed out 24 hours later. At the end of Chiron exposure (72 hAAF), gastruloids are spherical and express T/Bra throughout. Between 72 and 96 hAAF, gastruloids start to elongate and display AP-graded T/Bra and E-cadherin (E-cad) expression patterns with higher expression posteriorly (Figure 1a-d). These heterogeneousities are linked to an epithelial-to-mesenchymal transition (as evidenced by the differential expression of Snail, Figure 1d, Supplementary Figure 1), which has since been identified in other differentiation protocols [54]. The main tissue types produced after symmetry breaking in our differentiation protocol are endoderm and mesoderm ([15], Supplementary Figures 2-5).

We used the two cell lines, E-cad-GFP/Oct4-mCherry and T-Bra-GFP/NE-mKate2 (Histone GFP under the T-Bra promoter). In brief, between 100 and 200 cells were seeded in low adhesive 96 wells plates in a neural differentiation medium. Note that, following the usual convention in the field, when indicating time points in hAAF, the time of aggregate formation refers to the time when the cells were seeded. In order to promote the formation

of endoderm formation and to reduce symmetry breaking variability compared to the standard gastruloids protocol, we enriched the medium with Fibroblast Growth Factor and Activin. A complete description of the cell lines, the culture conditions and the protocol for making gastruloids is presented in [15]. All cell lines were tested to be free from mycoplasma contamination using qPCR.

**Two-photon imaging.** The gastruloids were imaged in a chamber maintained at 37°C, 5% CO2 with a humidifier. We performed two-photon imaging using a Zeiss 510 NLO (Inverse - LSM) with a femtosecond laser (Laser Mai Tai DeepSee HP) with a 40 x/1.2 C Apochromat objective. 72hAAF gastruloids were transferred from 96-well plates to either MatTek dishes (MatTek corporation, ref: P35G-1.5-14 C) or to microwell plates (500 µm wells, SUNBIOSCIENCE ref: Gri3D), which enable to image in parallel up to four gastruloids thanks to the small distance between two neighbouring wells, Supplementary Figure 18. To visualize inter-cellular space, we added 10 µL/mL of a Sulforhodamin B solution (1 mg/mL of Sulforhodamin B powder, 230,162 Aldrich) to the medium and excited both GFP and Sulforhodamin B at 900 nm. To image cell nuclei, we added 1 µL of a SPY555-DNA (spirochrome) 1000x stock solution per mL of medium and excited both GFP and SPY555-DNA at 880 nm. Nuclei staining progressively improved over time, likely due to a progressive decrease in cell packing density and to nuclei accessibility (Supplementary movies 2-5).

We imaged sub-volumes of gastruloids around their midplane (z-stacks of 60 µm, with 12 to 30 slices at intervals between 5 µm and 2 µm depending on the movie) for a total duration between 10 h and 24 h (time interval between two frames between 3 and 10 min depending on the experiment). We imaged the GFP and the red fluorophores (either Sulforhodamin B or SPY555-DNA) channels using two non-descanned detectors. We ensured that the T/Bra and E-cad polarization direction, as well as the flow field directionality was approximately aligned with the imaging plane.

**Immunostaining.** Samples were immunostained using the same protocol than in [15] apart from the fact that primary antibodies were incubated during 3 consecutive days. Secondary antibodies (AF488, A568, AF647) and Hoechst 33342 were chosen to be compatible with four-color two-photon imaging (Antibodies are in Extended Data Table 1).
<table><thead><tr><th colspan="5">Antibodies</th></tr><tr><th>Antibody</th><th>Species</th><th>Reference</th><th>Provider</th><th>Dilution</th></tr></thead><tbody><tr><td>Brachyury</td><td>Goat</td><td>AF2085</td><td>RD system</td><td>1:40</td></tr><tr><td>E-cadherin</td><td>Rat</td><td>M108</td><td>Takara</td><td>1:200</td></tr><tr><td>FoxA2</td><td>Goat</td><td>sc-6554</td><td>Santa-Cruz</td><td>1:50</td></tr><tr><td>Snail</td><td>Goat</td><td>AF3639</td><td>RD system</td><td>1:50</td></tr><tr><td>OB-cadherin</td><td>Rabbit</td><td>71-7600</td><td>Invitrogen</td><td>1:500</td></tr><tr><td>Flk1</td><td>Rat</td><td>14-5821-85</td><td>Ebioscience</td><td>1:100</td></tr><tr><td>VEGFR2</td><td>Goat</td><td>AF644</td><td>BioTechnne</td><td>1:200</td></tr><tr><td>Snail</td><td>Rabbit</td><td>C15D3</td><td>Cell Signaling</td><td>1:100</td></tr><tr><td>Anti-Rat IgG H&amp;L (Alexa Fluor® 568)</td><td>Donkey</td><td>ab175475</td><td>Abcam</td><td>1/500</td></tr><tr><td>Anti-Goat IgG H&amp;L (Alexa Fluor® 647)</td><td>Donkey</td><td>A-21447</td><td>Invitrogen</td><td>1/500</td></tr><tr><td>Anti-Goat IgG H&amp;L (Alexa Fluor® 488)</td><td>Donkey</td><td>ab150129</td><td>Abcam</td><td>1/500</td></tr><tr><td>Anti-Rabbit IgG H&amp;L (Alexa Fluor® 546)</td><td>Donkey</td><td>A10040</td><td>Invitrogen</td><td>1/500</td></tr><tr><td>Anti-Rat IgG H&amp;L (Alexa Fluor® 647)</td><td>Donkey</td><td>A-21247</td><td>Invitrogen</td><td>1/500</td></tr><tr><td>Anti-Goat IgG H&amp;L (Alexa Fluor® 568)</td><td>Donkey</td><td>ab175704</td><td>Abcam</td><td>1/500</td></tr></tbody></table>

**Optical flow measurements.** To measure coarse-grained velocity fields, we use the signal collected in the red channel (either inter-cellular space with the Sulforhodamin B (SRB) or nuclear signal with the SPY555-DNA probe). The SRB signal corresponds to (1) cell contours, (2) dead cells, and (3) bright spots of subcellular size moving rapidly, corresponding probably to vesicles. As dead cells and bright spots are much brighter than cell contours, we filter them with a thresholding operation. The SPY555-DNA signal also accumulates more in bright spots and dead cells which are filtered by the same thresholding process. The mid-plane of the z-stack is determined (maximal area of the gastruloid) and an average projection is made between the mid-plane and the two planes under and above it (separated usually by 2 microns). The same averaging is performed on the GFP channel. We then perform 2D timelapse registration using the ImageJ Plugin MultiStackReg to remove rapid rotations and translations of the aggregate (the transformation is calculated on the red channel and applied to both red and green channels). We then use a custom-made Matlab optical flow code based on the Kanade Lu-

cas Tomasi (KLT) algorithm with a level 2 pyramid and a window size corresponding to a square of 2 cells by 2 cells [25]. The velocity fields are averaged over 12 consecutive frames (1 hour) to smooth both measurement noise and temporal fluctuations. In Supplementary Figure 7, we compare the velocity maps obtained by applying the optical flow method on two distinct signals (SRB versus nuclei or raw nuclei signal versus locally equalized nuclei signal).

**Dipole moment decomposition.** In our plots, we normalize the polarization vector $P(t)$ at each time point by some reference value $P_0(t)$ that depends on the instantaneous spatial standard deviation of the T/Bra or E-cad expression level and the total cross-sectional area of the gastruloid (Figure 2b,c, section 1 of the SI). This allows us to compensate for both time dependence and sample variability of overall expression level and size when comparing different gastruloids. Compared to a time-independent normalization, this resulted in a better collapse of the polarization dynamics across gastruloids (Figure 2c and Supplementary Figure 8b).

We decompose the time derivative of $P$ into four parts (Figure 2e, Methods, section 1 of the SI):

$$ \frac{d\mathbf{P}}{dt} = \mathbf{R} + \mathbf{A} + \mathbf{G} + \mathbf{Q}_{3D}. \quad (5) $$

with

$$ \mathbf{R}(t) = \int_{\Omega} \left( \frac{\partial \phi}{\partial t} + \mathbf{v} \cdot \nabla \phi \right) (\mathbf{r} - \mathbf{c}) \, dS \quad (\text{RD}) \quad (6) $$

$$ \mathbf{A}(t) = \int_{\Omega} \phi(\mathbf{v} - \bar{\mathbf{v}}) \, dS \quad (\text{advection}) \quad (7) $$

$$ \mathbf{G}(t) = \int_{\Omega} (\phi - \bar{\phi})(\nabla \cdot \mathbf{v})(\mathbf{r} - \mathbf{c}) \, dS \quad (\text{growth}). \quad (8) $$

Here, $\mathbf{v}$ is the local tissue velocity (Figure 1c) and a line over a symbol, $\bar{\cdot}$, denotes averaging over the aggregate cross section.

**Mode decomposition of the tissue flow.** We define a disk that encompasses most of the aggregate projection (Figure 3a), and express the coarse-grained tissue velocity $\mathbf{v}(\mathbf{r})$ in polar form, i.e. with radial and angular velocity components $v_r(\mathbf{r})$ and $v_{\theta}(\mathbf{r})$ (Figure 3b). We then decompose both spatial functions $v_r(\mathbf{r})$ and $v_{\theta}(\mathbf{r})$ into different modes, again using polar coordinates. To capture the variation of both functions with the angular coordinate $\theta$, we decompose them into angular Fourier modes (mode index $n$, Figure 3c). Similarly, we decompose the radius-dependent part of both functions into polynomial modes (mode index $p$, Figure 3c). As a result, both $v_r(\mathbf{r})$ and $v_{\theta}(\mathbf{r})$ are decomposed into a sum of base modes $v_{n,p}^r(\mathbf{r})$ and $v_{n,p}^{\theta}(\mathbf{r})$, each weighted by a
complex amplitude, denoted $A_{n,p}^r$ and $A_{n,p}^\theta$:

$$v_r(\mathbf{r}) = \sum_{n,p} A_{n,p}^r v_{n,p}^r(\mathbf{r}) \quad (9)$$

$$v_{\theta}(\mathbf{r}) = \sum_{n,p} A_{n,p}^{\theta} v_{n,p}^{\theta}(\mathbf{r}). \quad (10)$$

We display a few of the base modes $v_{n,p}^r(\mathbf{r})$ and $v_{n,p}^\theta(\mathbf{r})$ in Figure 3d. Generally, the modes with smaller indices $(n, p)$ correspond to larger-scale features of the flow field. Each of the complex mode amplitudes $A_{n,p}^{r/\theta}$ provides information about an angle and a magnitude. The angle corresponds to a rotation of the respective base mode around the barycenter, and the squared magnitude reflects the contribution of the respective base mode to the spatial variance of the flow field (Figure 4a,b; section 2 of the SI). In other words, the modes with the largest amplitudes correspond to the dominant features of the flow field.

**Residuals between measured and simulated flows.** We further quantitatively compared the flow fields obtained from our model to the experimentally observed one. We quantify the deviation between simulated and experimentally observed flow fields using a residual, which we define as the fraction of the spatial variance in the experimental flow field that is not reproduced by our simulations (section 3 of the SI). In Figure 5f left, we plot the residual depending on the two model parameters $\Delta\kappa_s$ and $r$. We find indeed that the residual is smallest for negative $\Delta\kappa_s$. However, we also find that the deviation becomes larger for a small interface-to-surface stress ratio $r$. This is because, for small $r$, internally created tension are smaller than surface tension, which induces flows that would lead to a more circular aggregate shape (e.g. $(\Delta\kappa_s, r) = (0, 0.5)$ in Figure 5d). However, we observe that even at negative $\Delta\kappa_s$ and large $r$, the residual is still very large, on the order of 85%.

To dissect the possible reasons for the relatively large residual between measured and predicted flow field (Figure 5f left), we suspected that small-scale spatio-temporal fluctuations contribute to the experimental flow fields. We thus compared between simulation and experiment only the flow contributions provided by the modes $v_{1,1}^\theta$, $v_{1,0}^\theta$ and $v_{1,0}^r$, which together create the recirculating flow pattern. The residual between just these contributions between simulation and experiment is shown in Figure 5f right. We find that our simulations predict the observed recirculation flow in a wide parameter region defined by $\Delta\kappa_s < 0$ with a residual that is reduced to ~40–50%, independent of the precise values of $\Delta\kappa_s$ and $r$. However, varying $\Delta\kappa_s$ and $r$ does affect the amplitude of other large-scale flow components that do not contribute to the recirculating flow. For example, reducing the values of $r$ and $|\Delta\kappa_s|$ results in an important flow contribution by the mode $v_{2,0}^r$ (Figure 5g), which corresponds to a rounding up of the gastruloid (e.g. $(\Delta\kappa_s, r) = (0, 0.5)$ in

Figure 5d).

**Fusion experiments.** Two gastruloid batches were prepared from the T/Bra GFP line with a 14 hours delay from one another. Both batches were incubated during the 24h long Chiron pulse in different CellTracker colors (CMTMR orange CellTracker, Thermofischer) for the older batch and blue (CMAC blue CellTracker, Thermofischer) for the youngest batch. They were rinsed in transparent medium (Neural differentiation medium enriched in FGF and Activin) at the time of fusion ($t = 72hAAF$ for the youngest gastruloid and $t = 86hAAF$ for the oldest one) and were put in contact in pairs (either a young/young pair or a young/old pair) in an individual well from a 96 wells plate. Pairs of fusing gastruloids were imaged with epifluorescence and brightfield at 10x magnification with a Zeiss microscope every 20 minutes in the red, blue, green and brightfield channels using a multiposition configuration.

**Fusion interface/angle analysis.** We used the T/Bra and CellTracker color channels to define three scalar fields: (i) a field $\psi(\mathbf{r})$ equal to 1 inside any one of the two gastruloids and 0 outside, (ii) a field $\chi(\mathbf{r})$ equal to 1 in the first aggregate (e.g. blue in Figure 6c) and -1 in the second aggregate (e.g. red in Figure 6c), and (iii) a field $\phi(\mathbf{r}) \in [0, 1]$ representing the T/Bra signal. We defined the fusion interface as the iso-line $\chi = 0$ with curvilinear coordinate $s_\chi$. The fusion interface length was then defined as

$$L = \int_{\chi=0} \psi \, ds_{\chi}. \quad (11)$$

In addition, the length of the homotypic ++ interface between T/Bra+ tissues was computed as

$$L_{++} = \int_{\chi=0} \phi \, ds_{\chi}. \quad (12)$$

The heterotypic +− interface length was then computed as $L_{+−} = L − L_{++}$.

We defined the gastruloids’ surface as the iso-line $\psi = 0.5$ with curvilinear coordinate $s_\psi$. The triple points where both gastruloids and the outer medium meet are then characterized by the coordinate $s_\psi^t$ on the aggregate surface that satisfies $\chi(s_\psi^t) = 0$. At a given triple point, we computed the outer fusion angle $\theta_1$ using

$$\cos(\theta_1) = n_{\psi}(s_{\psi}^{t} + \delta_{s}) \cdot n_{\psi}(s_{\psi}^{t} - \delta_{s}), \quad (13)$$

where $n_\psi = -\nabla\psi/|\nabla\psi|$ is the outward-pointing unit vector normal to the aggregate surface and $\delta_s$ is a length scale parameter. We checked that the choice of $\delta_s$ does not affect the results in Figure 6.

Inner angles $\theta_2$ and $\theta_3$ were determined following a similar approach, but using the gastruloid-gastruloid interface $\chi = 0$ with curvilinear coordinate $s_\chi$. At a triple point, the coordinate $s_\chi^t$ on the gastruloid-gastruloid in-
terface satisfies $\psi(s_{\chi}^t) = 0.5$. The inner angles were computed using

request.

$$ \cos(\theta_2) = n_{\psi}(s_{\psi}^t + \delta_s) \cdot n_{\chi}(s_{\chi}^t - \delta_s) \quad (14) $$

and

$$ \cos(\theta_3) = n_{\psi}(s_{\psi}^t - \delta_s) \cdot n_{\chi}(s_{\chi}^t - \delta_s), \quad (15) $$

where $n_{\chi} = \nabla\chi/|\nabla\chi|$ is the unit vector normal to the fusion interface on the iso-line $\chi = 0$.

## DATA AVAILABILITY

All data that support the findings of this study are available from the corresponding author upon reasonable request.

## CODE AVAILABILITY

All code that support the findings of this study are available from the corresponding author upon reasonable

## STATISTICS AND REPRODUCIBILITY

No statistical method was used to predetermine sample size. The experiments were not randomized, and the investigators were not blinded to allocation during experiments and outcome assessment. The following data were excluded from the analyses: instances where tissue flows had an excessively strong component out of the imaging plane (orthogonal to the optical axis), preventing a bidimensional flow analysis (see Methods and Supplementary Information for further details).
# Supplementary Information

## Marangoni-like tissue flows enhance symmetry breaking of embryonic organoids

Simon Gsell<sup>a,b,c</sup>, Sham Tlili<sup>a</sup>, Matthias Merkel<sup>b</sup>, Pierre-François Lenne<sup>a</sup>

<sup>a</sup> Aix Marseille Univ, CNRS, IBDM, Turing Centre for Living systems, Marseille, France

<sup>b</sup> Aix Marseille Univ, Université de Toulon, CNRS, CPT (UMR 7332), Turing Centre for Living systems, Marseille, France

<sup>c</sup> Aix Marseille Univ, CNRS, Centrale Med, IRPHE (UMR 7342), Turing Centre for Living systems, Marseille, France

## 1 Dipole moment decomposition

### Definition

We define the dipole moment of the 2D cross section $\Omega$ of the aggregate as

$$ \mathbf{P}(t) = \int_{\Omega} \phi(\mathbf{r}, t) [\mathbf{r} - \mathbf{c}(t)] \, dS. \qquad (S1) $$

Here and below, $\mathbf{r}$ is the integration variable, $\phi(\mathbf{r}, t)$ is the protein fluorescence intensity (of T/Bra or E-cad), $\mathbf{c}(t) = (\int_{\Omega} \mathbf{r} \, dS)/S(t)$ is the aggregate cross section's barycenter, and $S(t) = \int_{\Omega} dS$ is its area.

### Normalization

As a reference dipole, we consider a Janus particle, i.e. a disk with the same area $S$ as the aggregate cross section, and composed of two semi-circular surfaces with homogeneous concentrations $\tilde{\phi}$ and $-\tilde{\phi}$, respectively. The resulting reference dipole moment reads

$$ P_0 = \frac{4}{3} \tilde{\phi} R^3, \qquad (S2) $$

with $R = \sqrt{S/\pi}$ the radius of the Janus particle. For the characteristic concentration amplitude $\tilde{\phi}$, we choose the root-mean-square of the $\phi$-field within the aggregate cross section.

Since both $\tilde{\phi}$ and $R$ vary over time, we compared two different strategies for the normalization of the dipole moment:

1. Defining a time-independent $P_0$ according to Eq. (S2) based on $\tilde{\phi}$ and $R$ at the end of the polarization process (Figure S8b), or

2. defining a time-dependent $P_0(t)$ according to Eq. (S2) using the instantaneous values of $\tilde{\phi}(t)$ and $R(t)$ (Figure 2c in main text).

In the main text, we have decided to use the second strategy throughout, because we observed that it allows for a better collapse of the polarization curves of different aggregates (compare Figure 2c in main text to Figure S8b).

### Decomposition

Using the Leibniz integral theorem, the time derivative of $P$ reads

$$ \frac{d\mathbf{P}}{dt} = \int_{\Omega} \frac{\partial}{\partial t} (\phi[\mathbf{r} - \mathbf{c}]) \, dS + \int_{\partial\Omega} \phi[\mathbf{r} - \mathbf{c}](\mathbf{u}_b \cdot \mathbf{n}_{2D}) \, dl, \qquad (S3) $$
where $\partial\Omega$ is the boundary of the aggregate cross section, $\mathbf{n}_{2D}$ is a unit vector within the cross-sectional plane normal to the boundary $\partial\Omega$, and $\mathbf{u}_b$ is the apparent boundary velocity in the cross-sectional plane. Since the aggregate can also move and deform in the direction normal to the cross-sectional plane, the apparent velocity $\mathbf{u}_b$ is the sum of (1) the components of the actual tissue velocity that lie in the cross-sectional plane, $\mathbf{v}_{2D}$, and (2) an apparent planar velocity resulting from motion of the boundary normal to the cross-sectional plane:

$$ \mathbf{u}_b = \mathbf{v}_{2D} + v_z n_z^* \mathbf{n}_{2D}, \qquad (S4) $$

where $v_z$ is the velocity normal to the cross-sectional plane, and $n_z^* = n_z/\sqrt{n_x^2 + n_y^2}$ with $\mathbf{n} = (n_x, n_y, n_z)$ being the 3D unit normal vector at the aggregate boundary.

Using the divergence theorem and denoting the gradient operator in the cross-sectional plane by $\nabla$, Eq. (S3) becomes

$$ \frac{d\mathbf{P}}{dt} = \int_{\Omega} \frac{\partial}{\partial t} (\phi[\mathbf{r} - \mathbf{c}]) \, dS + \int_{\Omega} \nabla \cdot (\phi \mathbf{v}_{2D}[\mathbf{r} - \mathbf{c}]) \, dS + \mathbf{W}_{3D} \qquad (S5) $$

with

$$ \mathbf{W}_{3D} := \int_{\partial\Omega} v_z n_z^* \phi [\mathbf{r} - \mathbf{c}] \, dl. \qquad (S6) $$

This can be written as

$$ \frac{d\mathbf{P}}{dt} = \int_{\Omega} \left[ \frac{\partial\phi}{\partial t} + \mathbf{v}_{2D} \cdot \nabla\phi \right] [\mathbf{r} - \mathbf{c}] \, dS - \frac{d\mathbf{c}}{dt} S \bar{\phi} + \int_{\Omega} \phi (\mathbf{v}_{2D} + (\nabla \cdot \mathbf{v}_{2D})[\mathbf{r} - \mathbf{c}]) \, dS + \mathbf{W}_{3D}, \qquad (S7) $$

where $\bar{\phi} := (\int_{\Omega} \phi \, dS)/S$. The first term on the right-hand side of Eq. (S7) represents the reaction-diffusion contributions to dipole moment variations, the second and third correspond to advective contributions, and the last term, $\mathbf{W}_{3D}$, represents 3D effects.

To further reduce the 3D effect contribution, we apply the Leibniz rule to $d\mathbf{c}/dt$ in Eq. (S7):

$$ \frac{d\mathbf{c}}{dt} = -\frac{\mathbf{c}}{S} \frac{dS}{dt} + \frac{1}{S} \int_{\Omega} (\nabla \cdot \mathbf{v}_{2D}) \mathbf{r} \, dS + \frac{1}{S} \int_{\partial\Omega} v_z n_z^* \mathbf{r} \, dl, \qquad (S8) $$

where $\overline{\mathbf{v}_{2D}} := (\int_{\Omega} \mathbf{v}_{2D} \, dS)/S$. Substituting Eq. (S8) into Eq. (S7), we obtain

$$ \frac{d\mathbf{P}}{dt} = \mathbf{R} + \mathbf{A} + \mathbf{G} + \mathbf{Q}_{3D}, \qquad (S9) $$

with

$$ \mathbf{R}(t) := \int_{\Omega} \left[ \frac{\partial \phi}{\partial t} + \mathbf{v}_{2D} \cdot \nabla \phi \right] [\mathbf{r} - \mathbf{c}] \, dS \qquad (\text{reaction-diffusion}) \qquad (S10) $$

$$ \mathbf{A}(t) := \int_{\Omega} \phi [\mathbf{v}_{2D} - \overline{\mathbf{v}_{2D}}] \, dS \qquad (\text{advection}) \qquad (S11) $$

$$ \mathbf{G}(t) := \int_{\Omega} (\phi - \bar{\phi})(\nabla \cdot \mathbf{v}_{2D})[\mathbf{r} - \mathbf{c}] \, dS, \qquad (\text{growth}) \qquad (S12) $$

and where we introduced the new 3D correction:

$$ \mathbf{Q}_{3D}(t) := \int_{\partial\Omega} v_z n_z^* (\phi - \bar{\phi}) [\mathbf{r} - \mathbf{c}] \, dl. \qquad (S13) $$

The 3D effects contribution $\mathbf{Q}_{3D}$ is generally expected to be smaller than $\mathbf{W}_{3D}$ defined in Eq. (S6), because the definition of $\mathbf{Q}_{3D}$ is reduced by the “mean-field 3D effects $\mathbf{W}_{3D}$” that assume a homogeneous
φ field. Therefore, in the paper we use Eq. (S9) to quantify contributions to dipole moment variations in the cross-sectional plane, where for notational simplicity $v_{2D}$ is replaced by $v$.

## Analysis procedure

We use our live imaging data to quantify each of the terms in Eq. (S9). The velocity field $v_{2D}$ is obtained from an optical flow analysis (Figure 1c left in the main text, Figure S6), and the scalar field $\phi$ corresponds to T/Bra or Ecad fluorescence intensity (Figure 1c right in the main text, Figure S1 and Figure S4). We study the dynamics of 10 gastruloids, including 6 for the T/Bra reporter line and 4 for the E-cad reporter line. The cumulative polarization change during the acquisition time interval $[t_i, t_f]$ is

$$ \tilde{P}(t) := \int_{t_i}^{t} \frac{d\mathbf{P}}{dt} dt = \tilde{R}(t) + \tilde{A}(t) + \tilde{G}(t) + \tilde{Q}_{3D}(t), \qquad (S14) $$

where $\tilde{R}$, $\tilde{A}$, and $\tilde{G}$ are computed as the time integrals of Eqs. (S10), (S11), and (S12), respectively. Because we cannot extract the precise values of $v_z n_z^*$ along the boundary, we determined the 3D contribution $\tilde{Q}_{3D}$ as

$$ \tilde{Q}_{3D}(t) = \tilde{P}(t) - \tilde{R}(t) - \tilde{A}(t) - \tilde{G}(t), \qquad (S15) $$

where $\tilde{P}(t)$ is computed based on equation (S1).

We discarded samples that showed too strong 3D effects; more precisely, we discarded the samples where $|\tilde{Q}_{3D}(t_f)|/|\tilde{P}(t_f)| > 0.2$ at the end of the acquisition time, $t_f$. The decomposition of the cumulative polarization change for the 5 resulting samples is plotted in Figure S8c-e,g,i. Two of the five discarded samples are shown in Figure S8f,h.

# 2 Mode decomposition of tissue flows

## Flow decomposition

In a disk of radius $R$, we express tissue velocity in polar coordinates $v(r, \theta) = [v_r(r, \theta), v_\theta(r, \theta)]$, which we decompose into base modes $v_{n,p}^r(r, \theta)$ and $v_{n,p}^\theta(r, \theta)$ as follows:

$$ v_r(r, \theta) = \sum_{n=0}^{+\infty} \sum_{p=0}^{+\infty} A_{n,p}^r v_{n,p}^r(r, \theta), \qquad (S16a) $$

$$ v_\theta(r, \theta) = \sum_{n=0}^{+\infty} \sum_{p=0}^{+\infty} A_{n,p}^\theta v_{n,p}^\theta(r, \theta), \qquad (S16b) $$

where $A_{n,p}^r$ and $A_{n,p}^\theta$ are complex amplitudes.

We construct the base modes from a radial part given by polynomials with real coefficients, $P_p(r)$, and an angular part given by Fourier modes with index $n$:

$$ v_{n,p}^r(r, \theta) = v_{n,p}^\theta(r, \theta) = P_p(r) \frac{e^{in\theta}}{\sqrt{2\pi}}. \qquad (S17) $$

Thus, both radial and angular modes, $v_{n,p}^r(r, \theta)$ and $v_{n,p}^\theta(r, \theta)$, are formally defined in the same way. We still use a different notation only to differentiate whether they are used to construct the radial or the angular part of the velocity field $v(r)$.

For convenience, we want these base modes to fulfill the following orthonormality condition:

$$ \int_0^R \int_0^{2\pi} v_{n,p}^r(r, \theta) [v_{n',p'}^r(r, \theta)]^* r \, d\theta dr = \delta_{nn'} \delta_{pp'}, \qquad (S18) $$
where the asterisk * denotes complex conjugation. This implies that

$$ \int_{0}^{R} r \mathcal{P}_{p}(r) \mathcal{P}_{p'}(r) dr = \delta_{pp'}. \qquad (S19) $$

We thus constructed the polynomials $\mathcal{P}_p(r)$ following a standard modified Gram-Schmidt (MGS) procedure with the following inner product

$$ \langle f, g \rangle = \int_{0}^{R} rf(r)g(r) dr, \qquad (S20) $$

with $f$ and $g$ two continuous functions on $[0, R]$. This procedure leads to a variant of the so-called Jacobi polynomials and ensures validity of Eqs. (S18) and (S19). In our construction, we ensured that each polynomial $\mathcal{P}_p(r)$ is of degree $p$, respectively (compare Figure 3c in the main text).

As a consequence of having a complete orthonormal basis, we can extract the mode amplitudes $A_{n,p}^r$ and $A_{n,p}^\theta$ as follows:

$$ A_{n,p}^{r} = \int_{0}^{R} \int_{0}^{2\pi} v_{r}(r, \theta) [v_{n,p}^{r}(r, \theta)]^{*} r d\theta dr, \qquad (S21a) $$

$$ A_{n,p}^{\theta} = \int_{0}^{R} \int_{0}^{2\pi} v_{\theta}(r, \theta) [v_{n,p}^{\theta}(r, \theta)]^{*} r d\theta dr. \qquad (S21b) $$

Moreover, the overall variance of the flow field can be expressed as follows (Parseval's identity):

$$ V := \int_C \mathbf{v}^2(\mathbf{r}) dA = \sum_{n=0}^{+\infty} \sum_{p=0}^{+\infty} (|A_{n,p}^r|^2 + |A_{n,p}^\theta|^2). \qquad (S22) $$

Here, $C$ denotes the circular disk of radius $R$. Thus, the spatial variance of the velocity field (up to the area factor) corresponds to the sum of the absolute squares of the complex mode amplitudes. This allows us to quantify the contribution of a given mode, e.g. $v_{n,p}^r$, to the overall spatial flow field fluctuations, which is given by $|A_{n,p}^r|^2/V$. Correspondingly, in Figures 3f,g and 4g in the main text, we normalized the squared amplitudes $|A_{n,p}^r|^2$ and $|A_{n,p}^\theta|^2$ by $V$.

### Recirculation mode

We define the recirculation mode as

$$ \mathbf{v}_R(r, \theta) = A_{1,0}^r v_{1,0}^r e_r + [A_{1,1}^\theta v_{1,1}^\theta + A_{1,0}^\theta v_{1,0}^\theta] e_\theta, \qquad (S23) $$

where the complex amplitudes $A_{1,0}^r$, $A_{1,1}^\theta$ and $A_{1,0}^\theta$ are determined from the measured flow field $\mathbf{v}(\mathbf{r})$ through Eqs. (S21). In order to determine a statistical representation of the recirculation mode over samples, we compute these complex amplitudes for all samples and then define median amplitudes $\overline{A}_{1,0}^r$, $\overline{A}_{1,1}^\theta$ and $\overline{A}_{1,0}^\theta$, obtained by taking the median of the norms and the median of the phases. In Figure 3i in the main text, we use the phase of $\overline{A}_{1,0}^r$ to define the orientation of the recirculation mode.

## 3 Computational model

### Phase field model

We model gastruloids as incompressible viscous fluids using Eqs. (5) and (6) in the main text. The
internally generated stress $\sigma_i$ is computed based on $\phi$ as in Eq. (6) in the main text with coefficient $\kappa_i$, and the surface stress $\sigma_s$ is computed analogously, replacing $\phi$ by $\psi$ and $\kappa_i$ by $\kappa_s$. However, the aggregate surface stress also depends on the field $\phi(\mathbf{r})$ representing T/Bra or E-cad concentration. To account for this $\phi$-dependence of the surface stress, we use a single free-energy functional defining the thermodynamics of the ternary mixture composed of T/Bra+ tissues, T/Bra− tissues and external medium, from which the driving force can be derived as discussed in [2, 3]. The effective stresses only depend on 3 parameters controlling the respective interface tensions between T/Bra+ and T/Bra− tissues ($\kappa_i$), the surface tension between T/Bra+ tissues and the medium ($\kappa_s^+$) and the surface tension between T/Bra− tissues and the medium ($\kappa_s^-$).

## Computational approach

We use experimental T/Bra fields to compute the stresses $\sigma_i$ and $\sigma_s$. The tissue velocity $\mathbf{v}$ and pressure $\Pi$ are then computed using a lattice-Boltzmann method [3, 4] accelerated by a multigrid method [5]. The computational procedure can be summarized as follows:

* We mapped the T/Bra fluorescence signal in the cross-sectional plane to the phase field $\phi$ by rescaling its magnitude by a global factor and adding a global offset, such that the maximum intensity is at $\phi = 1$ and the minimum at $\phi = -1$. We further defined the phase field $\psi$ equal to 1 inside the aggregate and 0 outside, which we also did using the T/Bra channel.

* We smoothed both $\phi$ and $\psi$ using convolutions with a Gaussian with a standard deviation of ~ 20 µm.

* We projected the $\phi$ and $\psi$ fields onto a 256 × 256 grid that is suitable for computational simulations. In both the $x$ and $y$ directions, the grid extends over approximately twice the gastruloid size, allowing a suitable resolution of the interface between the gastruloid and the outer medium.

* Based on the stress coefficients $\kappa_i$, $\kappa_s^+$ and $\kappa_s^-$, we compute the force driving tissue flows using the approach from [2, 3].

* We perform lattice-Boltzmann iterations until a steady flow solution is obtained.

An overview of the computational results for all samples is provided in Figure S9 and Figure S10.

## Computational residual

We perform simulations over ranges of the stress coefficients $\kappa_s^+$, $\kappa_s^-$ and $\kappa_i$. In particular, tissue flows depend on two non-dimensional ratios,

$$r = \frac{2\kappa_i}{\kappa_s^+ + \kappa_s^-} \qquad (S24)$$

and

$$\Delta\kappa_s = \frac{\kappa_s^+ - \kappa_s^-}{\kappa_i}. \qquad (S25)$$

A third parameter $\kappa_i/\eta$ sets the tissue velocity magnitude, but does not affect the flow pattern. We thus define the computational residual as

$$R(c) = \frac{\langle (\mathbf{u} - c\mathbf{v})^2 \rangle}{\langle \mathbf{u}^2 \rangle}, \qquad (S26)$$

where $\langle \cdot \rangle$ corresponds to an average over the whole aggregate cross-section, $\mathbf{u}$ and $\mathbf{v}$ are the experimental and computational residuals, respectively, and $c$ is a free coefficient scaling the computational velocity
magnitude. We can show analytically that the residual is minimal for $c = \langle \mathbf{u} \cdot \mathbf{v} \rangle / \langle \mathbf{v}^2 \rangle$. The minimal residual thus reads

$$R = 1 - \frac{\langle \mathbf{u} \cdot \mathbf{v} \rangle^2}{\langle \mathbf{u}^2 \rangle \langle \mathbf{v}^2 \rangle}. \qquad (S27)$$

In the paper, we used Eq. (S27) to compute the computational residual (Figure 4f in the main text).

## 4 Contact angles between T/Bra+ and T/Bra− tissues

Our mechanical model assumes the existence of an interface tension between T/Bra+ and T/Bra− tissues and a surface tension between these tissues and the outer medium. In this framework, T/Bra+ and T/Bra− tissues thus behave as two immiscible fluid droplets and we can thus use classical theory for wetting mechanics to predict contact angles between both tissues at mechanical equilibrium.

We consider the ternary fluid configuration represented in Figure S13. Interfaces between fluids 1, 2 and 3 are governed by surface tensions $\sigma_{12}$, $\sigma_{13}$ and $\sigma_{23}$. At mechanical equilibrium, the balance between tension forces at one triple point results in the Neumann triangle relation,

$$\frac{\sigma_{12}}{\sin(\theta_3)} = \frac{\sigma_{23}}{\sin(\theta_1)} = \frac{\sigma_{13}}{\sin(\theta_2)}. \qquad (S28)$$

As a result, the equilibrium angles $\theta_1$, $\theta_2$ and $\theta_3$ are entirely determined by the surface tensions $\sigma_{12}$, $\sigma_{13}$ and $\sigma_{23}$.

The full geometry of the equilibrium configuration can be computed considering the conserved area of droplets 2 and 3. At equilibrium, fluid 2 forms a lens characterized by one arc of radius $R_2$ and angle $\alpha_{21}$ and one arc of radius $R_{23}$ and angle $\alpha_{23}$. Similarly, fluid 3 forms a lens characterized by two arcs of radii $R_3$ and $R_{23}$ and angles $\alpha_{31}$ and $\alpha_{23}$. Angles $\theta_1$, $\theta_2$, $\theta_3$ and $\alpha_{21}$, $\alpha_{31}$, $\alpha_{23}$ are related through

$$\theta_3 = \frac{1}{2}(\alpha_{31} - \alpha_{23}), \qquad (S29)$$

$$\theta_2 = \frac{1}{2}(\alpha_{21} - (2\pi - \alpha_{31})) \qquad (S30)$$

and $\theta_1 = 2\pi - \theta_2 - \theta_3$.

The area $A_2$ of fluid 2 can be expressed as

$$A_2 = R_2^2 \left[ \alpha_{21} - \pi - \frac{1}{2} \sin(\alpha_{21}) \right] + R_{23}^2 \left[ \alpha_{23} - \frac{1}{2} \sin(\alpha_{23}) \right]. \qquad (S31)$$

Similarly, the area of fluid 3 is

$$A_3 = R_3^2 \left[ \alpha_{31} - \pi - \frac{1}{2} \sin(\alpha_{31}) \right] - R_{23}^2 \left[ \alpha_{23} - \frac{1}{2} \sin(\alpha_{23}) \right]. \qquad (S32)$$

Together with the extra condition required to ensure the connectivity of the three arc lengths, equations (S31) and (S32) allow to determine the unknown radii $R_2$, $R_3$ and $R_{23}$ from the conserved fluid area $A_2$ and $A_3$.

## 5 Supplementary movies

### Supplementary movie 1

One hour long timelapse of 3 gastruloids used to compute recirculation velocity fields presented in Supplementary Figure 5 (top 78 hAAF, middle 80 hAAF and bottom 86 hAAF). On the left, T/Bra signal
and on the right Sulforhodamin B signal. Both signals are inverted for better visibility.

**Supplementary movie 2**

z-stack of a 86 hAAF T/Bra GFP gastruloid. On the left in green, T/Bra signal and on the right in purple, SPY555-DNA signal.

**Supplementary movie 3**

Timelapse of a T/Bra GFP gastruloid between 76 hAAF and 91 hAAF. On the left in green, T/Bra signal and on the right in purple, SPY555-DNA signal.

**Supplementary movie 4**

z-stack of a 76 hAAF E-cadherin GFP gastruloid. On the left in green, E-cadherin signal and on the right in purple, SPY555-DNA signal.

**Supplementary movie 5**

Timelapse of a E-cadherin GFP gastruloid between 76 hAAF and 97 hAAF. On the left in green, E-cadherin signal and on the right in purple, SPY555-DNA signal.

**Supplementary movie 6**

Timelapse of a recirculation velocity field (moving average during one hour) corresponding to the Supplementary movie 3.

**Supplementary movie 7**

3 hours long epifluorescence timelapse of fusions between 72 hAAF gastruloids (in blue Cell Tracker) with 84 hAAF gastruloids (in orange Cell Tracker).

**Supplementary movie 8**

Zoom on few representative cases of different fusion scenarios that depend on the relative position of both aggregates extracted from Supplementary movie 8.

**Supplementary movie 9**

16 hours long epifluorescence timelapses of homogeneous fusions between 72 hAAF gastruloids (in blue Cell Tracker) and 72 hAAF gastruloids (in orange Cell Tracker).

**Supplementary movie 10**

Example of timelapse that has been discared because the inferred out-of-plane flow component was too high.

**References**

[1] Steven J. Kattman, Tara L. Huber, and Gordon M. Keller. Multipotent flk-1+ cardiovascular progenitor cells give rise to the cardiomyocyte, endothelial, and vascular smooth muscle lineages. *Developmental Cell*, 11(5):723-732, 2006.

[2] Ciro Semprebon, Timm Krüger, and Halim Kusumaatmaja. Ternary free-energy lattice Boltzmann model with tunable surface tensions and contact angles. *Physical Review E*, 93(3):033305, 2016.

[3] S. Gsell and M. Merkel. Phase separation dynamics in deformable droplets. *Soft Matter*, 18(13):2672-2683, 2022.

[4] S. Gsell, U. D'Ortona, and J. Favier. Lattice-Boltzmann simulation of creeping generalized newtonian flows: theory and guidelines. *Journal of Computational Physics*, 429:109943, 2021.

[5] S. Gsell, U. D'Ortona, and J. Favier. Multigrid dual-time-stepping lattice Boltzmann method. *Physical Review E*, 101(2):023309, 2020.
<figure><img src="image_14.png" /></figure>

E-cadherin

<figure><img src="image_15.png" /></figure>

Snail

<figure><img src="image_16.png" /></figure>

E-cadherin

<figure><img src="image_17.png" /></figure>

Snail

<figure><img src="image_18.png" /></figure>

E-cadherin

<figure><img src="image_19.png" /></figure>

Snail

<figure><img src="image_20.png" /></figure>

E-cadherin

<figure><img src="image_21.png" /></figure>

Snail

<figure><img src="image_22.png" /></figure>

E-cadherin

<figure><img src="image_23.png" /></figure>

Snail

Figure S1: Epithelial-mesenchymal transition (EMT). (a-e) Immunostaining of E-cadherin and Snail for 5 different 80hAAF gastruloids (see Methods, Immunostainings). Each image corresponds to an average projection of a 20 µm-thick stack (20 slices), centered around the gastruloid mid-plane. Images are inverted (high-intensity pixels appear in black) for better visibility.
<figure><img src="image_24.png" /><figcaption>Figure S2: Epithelial-mesenchymal transition evolution before and during symmetry breaking. Immunostainings of gastruloids between 48 hAAF and 80 hAAF with E-Cadherin, Snail, VEGFR2 and Dapi labeled (see Methods, Immunostainings). At 48 hAAF (just before the Chiron pulse), cells are all E-Cadherin positive and both Snail and Vegfr2 negative. At the middle of the 24h-long Chiron pulse (60 hAAF), Snail is expressed by a large proportion of cells and cells at the periphery of the gastruloid start to lose their E-Cadherin expression. At the end of the Chiron pulse (72 hAAF), a higher proportion of cells express Snail near the gastruloid surface and are both E-Cadherin negative and VEGFR2 positive. The center of the gastruloid is composed of a globally E-cadherin positive population with sparse cells which are Snail positive and E-cadherin negative. At 80 hAAF, the E-Cadherin negative and VEGFR2 positive cell population appears both at the anterior part of the gastruloid and at the surface of its posterior part. The letter A denotes for the anterior pole and P for the posterior one.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure S3: Schematics illustrating (a) the different cell populations within our gastruloids and (b) differential cadherin expression within these populations. (a) At 80hAAF, the gastruloids exhibit a stereotypical spatial separation between different cell populations. One pole expresses T/Bra, which indicates nascent mesoderm, mixed in a salt-and-pepper manner with cells expressing FoxA2, which indicates endoderm. The complementary pole is characterized by the absence of T/Bra and FoxA2 positive cells and the expression of Flk-1, which indicates more differentiated mesoderm, since Flk-1 is a hemato-cardiovascular cell lineages marker [1]. (b) At the coarse-grained level, the pole with many T/Bra-positive cells expresses high levels of E-cadherin and low levels of OB-cadherin, while the complementary population (Flk1 positive) expresses high levels of OB-cadherin and low levels of E-cadherin.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure S4: The E-cadherin positive population is a mix of endoderm (FoxA2+ cells) and nascent mesoderm (T/Bra+ cells). (a-d) Immunostaining of T/Bra, FoxA2 and E-cadherin for 4 different 80hAAF gastruloids (see Methods, Immunostainings). Each image corresponds to an average projection of a 20 µm-thick stack (20 slices) centered around the gastruloid mid-plane. Images are inverted (high-intensity pixels appear in black) for better visibility.</figcaption></figure>
(a)

<figure><img src="image_27.png" /><figcaption>Figure S5: (a) Immunostaining of T/Bra and Flk-1 for 4 different 80hAAF gastruloids (see Methods, Immunostainings).</figcaption></figure>
<figure><img src="image_28.png" /></figure>

Figure S6: Quantified velocity fields associated to the symmetry breaking superimposed on the T/Bra
field (inverted for better visibility). (a-f) Velocity fields are obtained by averaging over one hour (12
consecutive frames) the velocity fields were obtained by optical flow, performed on an average projection
of a 10 µm-thick stack (5 slices, Sulforhodamin B) centered around the gastruloid mid-plane (see Methods,
Optic Flow).
<figure><img src="image_29.png" /></figure>

Figure S7: Variations in velocity fields depending on the nature of the signal used in the optical flow pipeline. (a-b) Raw velocity fields obtained between two consecutive time points separated by 5 minutes during symmetry breaking: (a) on the T/Bra GFP signal (with enhanced contrast to improve visibility of even low-T/Bra cells), and (b) on the Sulforhodamin B signal. (c) Comparison between the velocity components obtained in (a) and (b) with the x component of the velocity in red and the y component of the velocity in brown. The black line corresponds to equality between results from a and b. (d-e) Raw velocity fields obtained between two consecutive time points separated by 5 minutes at early symmetry breaking: (d) on the raw nuclei signal (SPY555-DNA), and (e) on the signal with a locally equalized contrast in order to improve visibility of the inner nuclei. (f) Comparison between the velocity components obtained in (d) and (e) with the x component of the velocity in red and the y component of the velocity in brown. The black line corresponds to equality between results from d and e.
<figure><img src="image_30.png" /><figcaption>Figure S8: Dipole moment analysis. (a) Trajectories of the *x* and *y* components of the dipole moment *P* = (*P*<sub>x</sub>, *P*<sub>y</sub>), for different samples. The *y* axis is defined as the main axis of elongation of the aggregate at the last time point. (b) *P*<sub>y</sub> component of the dipole moment over time, for different samples. (c-l) Cumulative contributions to the polarization for 10 gastruloids. These include the 6 samples considered in Figure 2 of the main text, while 4 gastruloids were discarded (red frames) because of a too large 3D contribution (see SI section Analysis procedure).</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>Figure S9: Simulations of tissue flows during symmetry breaking of 9 gastruloids. The first column indicates the sample name and the relative contribution of the recirculation mode to the experimental flow fields. For each sample, we show (a) the computational residual as a function of the surface stress difference $\Delta\kappa_s$ and the interface-to-surface stress ratio $r$, (b) the residual comparing the recirculation flow, Eq. (S23), of simulated and experimentally observed velocity fields (columns f and g), (c) the experimental T/Bra field, (d) the experimental tissue velocity field, (e) the computational velocity field for $r = 2.5$ and $\Delta\kappa_s = -1$, (f) the reconstructed experimental recirculation flow and (g) the reconstructed computational recirculation flow for $r = 2.5$ and $\Delta\kappa_s = -1$.</figcaption></figure>
<figure><img src="image_32.png" /></figure>

Discarded samples - Weak polarization - Weak recirculation

<caption>Figure S10: Simulations of tissue flows during symmetry breaking of 7 gastruloids. The first column indicates the sample name and the relative contribution of the recirculation mode to the experimental flow fields. For each sample, we show (a) the computational residual as a function of the surface stress difference $\Delta\kappa_s$ and the interface-to-surface stress ratio $r$, (b) the residual comparing the recirculation flow, Eq. (S23), of simulated and experimentally observed velocity fields (columns f and g), (c) the experimental T/Bra or Ecad field, (d) the experimental tissue velocity field, (e) the computational velocity field for $r = 2.5$ and $\Delta\kappa_s = -1$, (f) the reconstructed experimental recirculation flow and (g) the reconstructed computational recirculation flow for $r = 2.5$ and $\Delta\kappa_s = -1$.</caption>
<figure><img src="image_33.png" /><figcaption>Figure S11: Tissue flows during symmetry breaking for one gastruloid with Ecad cell line, at different times. (a) Dipole moment decomposition. (b) Relative contribution of the recirculation flow in the experimental flow fields. (c) Ecad field, experimental and computational flow fields. Computational flows are computed for $r = 2.5$ and $\Delta\kappa_s = -1$.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>Figure S12: Tissue flows during symmetry breaking for one gastruloid with T/Bra cell line, at different times. (a) Dipole moment decomposition. (b) Relative contribution of the recirculation flow in the experimental flow fields. (c) T/Bra field, experimental and computational flow fields. Computational flows are computed for $r = 2.5$ and $\Delta\kappa_s = -1$.</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>Figure S13: Schematic representation of the equilibrium between two immiscible droplets.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Figure S14: Simulation of fusing droplets with different diameter ratios $D_1/D_2$, for $r = 0.8$ and $\Delta\kappa_s = -0.2$. (a) $D_1/D_2 = 1$. (b) $D_1/D_2 = 1/2$.</figcaption></figure>
<figure><img src="image_37.png" /><figcaption>Figure S15: Quantification of the overall cell motion during symmetry breaking. (a) We use the coarse-grained tissue velocity **u** to determine the average tissue rotation rate *N*(*t*) on both sides of the aggregate as well as the cumulative tissue rotation *R*(*t*). (b) Over 6 samples, the average cumulative rotation is between 1/16 and 1/8 complete rotation by the end of the symmetry breaking process. The symbols and colors are the same as in Figure 2(c) in the main manuscript. (c) Even when the computed average rotation is small, in some regions cells can move over approximately half of the aggregate length during polarization, as shown by selected trajectories computed using the coarse-grained tissue flow. This sample corresponds to the white circle symbols in (b). On trajectories, blue and red indicate initial and final position, respectively.</figcaption></figure>
<figure><img src="image_38.png" /></figure>

Figure S16: Spatial anti-correlation of E-cadherin and OB-cadherin (Cadherin 11). (a) Immunostaining of E-cadherin and OB-cadherin were performed for 3 different 76hAAF gastruloids (see Methods, Immunostainings, and Supp Movie OB-Ecad). Each image corresponds to a coarse-grained intensity (CGI) map of cadherins (each pixel in the map is spaced by 7 µm and its intensity corresponds to the raw image intensity averaged in a disk of radius 14 µm). Each raw image is obtained by averaging 10 µm-thick stacks (10 slices) at different z locations in the gastruloid. CGI values of OB-cadherin are plotted as a function of CGI values of E-cadherin in the same place and fitted linearly. (b) Slopes of the linear fits in panel (a) as a function of their respective *R*² values.
<figure><img src="image_39.png" /></figure>

Figure S17: At a coarse-grained level, T/Bra expression is positively correlated with tissue packing density and negatively correlated with cell-scale motion. (a) For a gastruloid mid-plane at 80hAAF, the T/Bra signal (inverted) is used to manually divide the cross-section into a T/Bra positive region (red) and a T/Bra negative region (green). (b,c) Regions as defined in panel (a) superimposed on the Sulforhodamin B signal, which indicates inter-cellular space. (b) Average ("SRB") over one hour (12 frames), which we use as a proxy for cell packing density. (c) Variance ("Fl") over one hour (12 frames), which we use as a proxy for cell-scale motion. (d) Ratios of both SRB and Fl as defined in panels (b) and (c) between the T/Bra negative region (green) and the T/Bra positive region (red). Both SRB and Fl are higher in the T/Bra negative region indicating a lower cell packing density and a higher cell motility.
<figure><img src="image_40.png" /><figcaption>Figure S18: Gastruloids random seeding in the microwells system (SUNBIOSCIENCE ref: Gri3D). 72 hAAF gastruloids are seeded in the system. Each well has zero, one or multiple gastruloids after the random seeding. A series of 3 or 4 neighbouring gastruloids is chosen (as the ones numbered from 1 to 4 in the schema) for parallel imaging, in order to minimize the distance covered with the water objective to avoid bubbles emergence over time.</figcaption></figure>