Embryonic Tissues as Active Foams

The physical state of embryonic tissues emerges from non-equilibrium, collective interactions among constituent cells. Cellular jamming, rigidity transitions and characteristics of glassy dynamics have all been observed in multicellular systems, but there is no unifying framework to describe all these behaviors. Here we develop a general computational framework that enables the description of embryonic tissue dynamics, accounting for the presence of extracellular spaces, complex cell shapes and tension fluctuations. In addition to previously reported rigidity transitions, we find a distinct rigidity transition governed by the magnitude of tension fluctuations. Our results indicate that tissues are maximally rigid at the structural transition between confluent and non-confluent states, with actively-generated tension fluctuations controlling stress relaxation and tissue fluidization. Comparing simulation results to experimental data, we show that tension fluctuations do control rigidity transitions in embryonic tissues, highlighting a key role of non-equilibrium tension dynamics in developmental processes.

state of the tissue 13 .
Here we develop a computational framework to study the dynamics, structure and physical state of active multicellular systems, including embryonic tissues. The simulations consist of a fully dynamic vertex model that accounts for extracellular spaces, complex cell shapes and the dynamics of tensions at cell-cell contacts. We recapitulate known density-dependent and densityindependent rigidity transitions, and show the existence of a qualitatively distinct rigidity transition controlled by the magnitude of active tension fluctuations that occurs both in confluent and nonconfluent systems. Comparing these simulations to quantitative data, our results indicate that the observed rigidity transition guiding vertebrate axis elongation 13 is mainly controlled by tension fluctuations.

Dynamic Vertex Model with Extracellular Spaces
To study the dynamics of embryonic tissues, we generalize 2D vertex models by accounting for  (Fig. 1g), which follow from force balance, namely where R α is the position of the vertex α, η R is a friction coefficient characterizing the dissipation associated with moving a vertex, T ij is the effective tension at the contact between cell i and j (with F (α) representing the set of all cells sharing vertex α) and N ij are the normal forces acting on vertex α. Finally, Θ(·) represents the Heaviside step function and prevents unrealistic negative tensions 13,24 .
Normal forces arise from osmotic pressure differences in adjacent cells, namely N ij = (∆P i − ∆P j ) L ij /2, where L ij and ∆P i (∆P j ) are, respectively, the contour length of the contact between cells i and j and the osmotic pressure difference across cell i (j). While we use the experimentally measured relation between osmotic pressure and cell volume 25 (area in 2D), the specific functional form does not qualitatively affect our results as long as cell volume decreases with increasing osmotic pressure (Methods).
To capture key observed features of the dynamics of effective tensions, T ij , namely their fluctuating nature and their finite persistence time 13 , in a simple, yet generic manner, we write them as an Ornstein-Uhlenbeck process with a tension that fluctuates around a fixed point T 0 ij and has a persistence time τ T 26, 27 , specifically where ∆T is the amplitude of tension fluctuations and ξ is gaussian white noise ( Fig. 1f; Methods).
The fixed point effective tensions depend on both the average cortical tension, T 0 , and average strength of cell-cell adhesion, W , which, like ∆T , are different at cell-cell contacts and free cell boundaries ( Fig. 1b; Methods).
Scaling all quantities with characteristic scales, we obtain a dimensionless set of equations as well as the relevant dimensionless parameters (Methods). The dimensionless parameter ∆T /T 0 characterizes the magnitude of tension fluctuations (Fig. 1e,f), the ratio W/T 0 characterizes the relative strength of cell-cell adhesion and cortical tension (Fig. 1b), the parameter P 0 L 0 /T 0 measures the relative magnitude of normal to tensional forces, the ratio τ T /τ R of the persistence time of tension fluctuations and the characteristic timescale τ R ≡ η T L 0 /T 0 of dissipation at vertices and, finally, the normalized system density ρ ≡ N L 2 0 /A T , with A T and N being the total area of the system and the number of cells, respectively (Fig. 1d). Since P 0 L 0 /T 0 and τ T /τ R can be estimated from existing experimental data (Methods), we focus on the parameter space spanned by ∆T /T 0 , W/T 0 and ρ.

Structural Transitions at Equilibrium
In the absence of tension fluctuations (∆T /T 0 = 0), the amount of extracellular spaces in equilibrium configurations is determined by force balance and varies with both the cell density ρ and the relative cell adhesion strength W/T 0 (Fig. 2a). Increasing cell density results in larger cellular volume fraction φ (Fig. 2b) and cell contact number z (Fig. 2c), with the system eventually becoming confluent (φ = 1) at an adhesion-dependent critical density, ρ c (W/T 0 ). For any fixed cell density, the system undergoes a non-confluent to confluent transition as W/T 0 is increased, since higher adhesion promotes stronger cell-cell contacts (Fig. 2a,b).
Changes in relative cell adhesion not only affect the volume fraction but also the structural  In the limit of vanishing cell adhesion, the system should behave as foams/emulsions, which display a jamming transition at a critical value φ c . The isostatic condition (z = 4 in 2D) sets the critical volume fraction φ c ≈ 0.83 of the system 30 ( Fig. 2f). Both this φ c value and the power-law dependence of z (Fig. 2f), are in agreement with recent equilibrium simulations of deformable particles 31 , indicating that our description accurately describes the foam limit.
Previously reported density-independent rigidity transitions in confluent states are also captured within this framework. In these transitions, which are controlled by the cell shape factor s (defined as s = P/ √ A, with P and A being the cell perimeter and area, respectively), the system switches from a solid to a fluid state at s c 3.81, with cells transiting from isotropic to anisotropic shapes 15 . Importantly, the fluid state in these descriptions is characterized by vanishing effective tensions (T 0 ij = 0), allowing neighbor exchanges at no energetic cost 17 . In the present framework, density-independent transitions are recapitulated at W/T 0 = 2, which leads to vanishing effective tensions and associated fluid state, as well as a sharply larger average shape factor due to the emergence of anisotropic cell shapes (Fig. 2g).

Dynamic Simulations with Finite Tension Fluctuations
In Comparing MSD values at long timescales (t = 10 2 τ T τ T > τ R ) shows that increasing tension fluctuations always leads to higher MSD values, regardless of cell adhesion strength.
However, for a fixed level of activity, MSD values vary non-monotonically as adhesion increases, displaying very reduced cell movements at the transition between non-confluent to confluent states shape factors in the confluent regime (Fig. 2g).

Physical state of multicellular systems
To directly assess the physical state of the system, we monitor shear stress relaxation after imposing an affine (shear) deformation ( Fig. 4a; Methods). Since tissue fluidization in biological tissues is associated with the non-linear mechanical response 13 , we impose large strains xy of 150% ( xy = 1.5). In the absence of tension fluctuations, the initial stress jump is largest with no adhesion (Fig. 4b), and vanishes when effective tensions vanish (W/T 0 = 2; Supplementary Information).
Subsequently, shear stress relaxes with a characteristic timescale τ R towards a constant value at long timescales, namely the yield stress σ Y 32 , which depends non-monotonically on cell adhesion ( Fig. 4c; Supplementary movie 5). For low relative adhesion W/T 0 , the system is non-confluent and the yield stress increases with increasing adhesion as extracellular spaces close down. In contrast, for adhesion values leading to confluence, increasing W/T 0 leads to a decreasing yield stress due to lower effective tensions. The vanishing yield stress at W/T 0 = 2 indicates a fluid state, in agreement with density-independent rigidity transitions ( Fig. 2g). At equilibrium, the system is maximally rigid (maximal yield stress) at the structural transition between confluent and non-confluent states (Fig. 4c,d).
Introducing finite tension fluctuations qualitatively changes stress relaxation at long timescales, which displays a slow stress decay rather than plateauing ( As in colloidal glasses 35 , whether the system behaves as a fluid or a solid depends on observation timescales. In embryonic tissues, if the characteristic timescale τ d of developmental processes is larger (smaller) than the stress relaxation timescale τ SR , namely τ d /τ SR the tissue behaves as a fluid (solid). Using typical developmental timescales (τ d ∼ 1-2 hours; τ d /τ T ∼ 10 2 ), we obtain the tissue phase diagram (Fig. 4h). Over a critical value of tension fluctuations, the tissue is always fluid regardless of cell adhesion levels. Below that critical activity value, the tissue is fluid at both low and high adhesion levels, but solid in between, in the region of the phase diagram surrounding the structural transition between non-confluent to confluent states.

Rigidity transition during vertebrate axis elongation
Recent experiments have shown that during body axis elongation in zebrafish, posterior tissues transit from a fluid state in the mesodermal progenitor zone (MPZ) to a solid state in the presomitic mesoderm (PSM) (Fig. 5a,b). To study the nature of rigidity transitions in these embryonic tissues, we compare our simulation results to both new and existing data 13 of the structure, dynamics and mechanics of the MPZ and PSM tissues (Fig. 5).
The presence of extracellular spaces in both the MPZ and PSM ( Fig. 5b; Methods) restrict the values of the relative adhesion parameter to W/T 0 0.6 ( Fig. 5c), with the cell density being the same in both regions 13 . Simultaneously fitting the distribution of cell-cell contact length fluctuations (Fig. 5d), the measured MSD (Fig. 5e), the cell volume fraction (Fig. 5f) and the yield stress (Fig. 5g), it is possible to reproduce all observations, provided that tensional fluctuations and relative cell adhesion are larger in the MPZ than in the PSM (∆T /T 0 = 0.5, 1.5 and W/T 0 = 0.3, 0.2 in PSM and MPZ, respectively). It is also possible to reproduce the observed higher cell boundary dynamics in the MPZ than in the PSM (Fig. 5h) and the relative changes in cell shape factor from PSM to MPZ (Fig. 5i) Indeed, cell movements in the PSM of N-cadherin (cdh2 −/− ) mutants (Methods), which lack Ncadherin adhesion 36 and are characterized by reduced cell volume fraction, lower yield stress but similar cell-cell contact fluctuations as wild type embryos 13 (Fig. 5d,f,g), show the same MSD behavior as wild type embryos and a solid-like tissue state (Fig. 5e). Placing the obtained parameter values for MPZ and PSM in the phase diagram indicates that the magnitude of tension fluctuations plays a key role in the control of the observed fluid-to-solid transition (Fig. 5c).

Discussion
We developed a general framework to study the structure, dynamics and mechanics of active (nonequilibrium) multicellular systems in a wide range of conditions. While jamming transitions and density-independent transitions are controlled by cell density and cell shape, respectively, we show

Methods
Pressure and volume relation. The (osmotic) pressure P of a cell has been experimentally shown to vary with its volume according to P = K/V 25 , meaning that changes in the osmotic pressure difference between the inside (P ) and outside (P 0 ) of the cell lead to changes in the cell volume V , with K characterizing the cell compressibility (we assumed for simplicity that cells are never compressed close to their dry mass limit). In our 2D description, the cell area A plays the role of the volume, so that P = K/A, with the characteristic cell size L 0 (or preferred area A 0 = L 2 0 ) being set by this relation, namely L 0 ≡ K/P 0 . However, the specific P-V (or P-A) functional form does not change our results qualitatively as long as the cell volume decreases with increasing applied osmotic pressure, namely P − P 0 = −Γ(A − A 0 ) (with Γ being a coefficient characterizing the cell's compressibility), and the deviations from the preferred cell area in 2D (A 0 ) are mild. Parameter estimation. Since the relaxation timescale τ R has been measured to be much smaller than the persistence timescale of tension variations (τ R 20s 37, 38 ; τ T 90s 13 ) we set the ratio to τ T /τ R = 10. While the values of osmotic pressure are unknown in vivo, they are expected to be larger than cortical stresses 25 , so that P 0 L 0 /T 0 1. Consequently, we fix P 0 L 0 /T 0 = 10 for all simulations. This ensures relatively mild cell size variations, as observed experimentally 13 MSD values were averaged over all cells and also all initial time points t 0 for any given set of model parameters.
Application of a step strain. A large step strain was applied by deforming the simulation box from a square to a parallelogram of the target shear strain (150%). An affine deformation is applied to all vertex positions when imposing the step strain. Lees-Edwards periodic boundary conditions were imposed throughout the simulation to avoid mismatch of cell geometry across system boundaries.
where ∆p i ≡ p i − p 0 , m, n are indices indicating the spatial direction (m = x, y and n = x, y), a i is the dimensionless cell area and l ij is the vector form of the edge length between cells i and j.
The shear stress term can be written as Images of extracellular spaces were obtained using an outcross of Tg(hsp70:secP-mCherry) p142 and Tg(bact2:mem-neonGreen-neonGreen) hm40 , to label the extracellular spaces and cell membranes, respectively. To trigger the expression of mCherry in the extracellular space, a 1 hour heat shock in a water bath at 39 • C was performed when the embryos were at 75% epiboly stage.
Confocal stacks through PSM and MPZ tissues were acquired with a z-step size of 0.55µm. Cell membranes were segmented using CellPose 43 . The resulting segmented image is overlaid on the confocal section showing both the secreted fluorescent reporter in the extracellular spaces and cell membranes (Fig. 5b).
Tracking of cell movements in the PSM of N-cadherin mutants. Both mutant (cdh2 −/− ) and sibling (cdh2 +/− ) embryos were injected with 80-100 pg each of H2B-mRFP mRNA and membrane-GFP mRNA at 1-2 cell stage to label nuclei and membranes, respectively. Confocal stacks through the PSM were acquired with a z-step size of 1µm and time interval of 1 minute for 30 mins, and processed using Imaris (Bitplane). Data were smoothed using a 1-pixel gaussian filter, to correct for photobleaching over time the normalize timepoints function was used, then attenuation correction was applied to correct for z-attenuation. After processing, data were cropped to the tissue of interest, then nuclei were detected using the spots function, and tracked using the Brownian motion algorithm. Nuclei positions were output for further analyses. The Mean Squared Relative Displacement (MSRD) is calculated from the cell trajectories 13 . MSRD in 3D corresponds to 3 times the MSD in 2D.
Measurement of cell shapes and shape factor. To measure cell shape and the shape factor we used Tg(h2afva:eGFP) kca6 xTg(actb2:memCherry2) hm29 double transgenic embryos (wild type phenotype) and N-cadherin (cdh2 −/− ) mutants injected with 80-100pg H2B-mRFP mRNA and 80-100 pg membrane-GFP mRNA at 1-2 cell stage. In all cases, confocal stacks were acquired with a z-step size of 0.55 µm. Segmentation analysis of membrane labeled images were performed in FIJI 44 . To measure cell shape, we chose the measurement plane of the confocal section when the nucleus area was maximized. The shape factor was calculated by fitting a polygon to the shape of the cell membrane signal. To visualize dynamics of cell shape contour, the template matching plugin 45 was first applied to measured timelapses of membrane-labeled images to remove drift and an average intensity projection image was obtained from the drift corrected timelapse images.
Data availability. All data is available from the authors upon reasonable request.
Code availability. The code developed for this manuscript is available from the authors upon reasonable request.