Structure and viscosity of non-confluent tissues:

(A) Bright-field single-plane images of an exemplary embryo of zebrafish before (t = −60 min), at the onset (t = 0 min), and after blastoderm spreading (t = 60 min). (B) Snapshot of 2D confocal sections at the 1st–2nd deep-cell layer of the blastoderm at t = 60 min. (A) and (B) are taken from [14]. (C) Viscosity η of zebrafish blastoderm as a function of ϕ in a log-linear scale using the data from [14]. The dashed line is the fit to VFT equation. Note that η does not change significantly beyond ϕ ≥ 0.87. (D) A typical snapshot taken from cell-based simulations for ϕ = 0.93. Cells are colored according to their radii (in μm) (color bar shown on the right). (E) The pair correlation function, g(r), as a function of r for ϕ = 0.93. The vertical dashed line is the position of the first peak (rmax = 17.0 μm). The pair correlation function does not exhibit signs of long-range order. Scale bars in (A) is 100 μm and in (B) is 50 μm. range order. Thus, the polydisperse cell system exhibits liquid-like structure even at the high ϕ.

Saturation in viscosity and relaxation time:

(A) Effective viscosity as a function of ϕ, with the solid line being the fit to VFT equation. The inset shows at high ϕ. The dashed line in the inset is the expected behavior assuming that the VFT relation holds at all ϕ. (B) The self-intermediate scattering function Fs(q, t) as a function of t for 0.70 ≤ ϕ ≤ 0.905. The dashed line corresponds to . (C) A similar plot for ϕ > 0.905. (D) The logarithm of the relaxation time τα(s) as a function of ϕ. The VFT fit is given by the dashed line. The inset shows a zoomed-in view for ϕϕS.

Spectrum of relaxation times:

(A) Scatter plot of relaxation times τα(s) as a function of cell radius. From top to bottom, the plot corresponds to decreasing ϕ. The vertical dashed line is for Ri = 4.25 μm, beyond which the τα changes sharply at high packing fractions. (B) Histogram P (ln(τα)) as a function of ln(τα). Beyond ϕ = 0.90 (ϕS), the histogram peaks do not shift substantially towards a high τα values. (C) For ϕϕS P (ln(τα)) (scaled by P max(ln(τα))) falls on a master curve, as described in the main text. (D) Same as (C) except the results are for ϕ > 0.90. The data deviates from the Gaussian fit, shown by the dashed line.

Density dependent cell-cell overlap :

(A) Probability of overlap (hij) between two cells, P (hij), for various ϕ values. The peak in the distribution function shifts to higher values as ϕ increases. (B) Mean ⟨hij⟩ = dhijP (hij) as a function of ϕ. Inset shows a pictorial illustration of h12 between two cells with radii R1 and R2 at a distance r12.

Changes in free area fraction with ϕ:

(A) Voronoi tessellation of cells for ϕ = 0.93 for a single realization. The orange circles represent actual cell sizes. The blue polygons show the Voronoi cell size. (B) Distribution of Voronoi cell size A as a function of ϕ. (C) Mean Voronoi cell size ⟨A⟩ as a function of ϕ. A zoomed-in view for ϕ > 0.860 is shown in the inset. (D) Distribution of free area P (Afree) for all ϕ. The vertical blue dashed line shows that the maximum in the distribution is at Afree ∼ 50μm2. (E) Free area fraction ϕfree as a function of ϕ. Note that ϕfree saturates beyond ϕ = 0.90. An expanded view of the saturated region is shown in the right panel of (E).

Relaxation in the polydisperse cell system is independent of the waiting time:

(A) Fs(q, t) for ϕ = 0.92 at different waiting times (τω = 106(s)). Regardless of the value of τω, all the Fs(q, t) curves collapse onto a master curve. (B) Relaxation time, ln(τα), as a function of τω. Over a 3 orders of magnitude change in tω, the variation in relaxation times is less than the sample-to-sample fluctuations, as shown by the error bar.

Parameters used in the simulation.

Fit of the stress-stress correlation functions to stretched exponential functions:

(A) The stress-stress correlation function ⟨Pμν(t)Pμν(0)⟩ divided by the value at t = 0 Pμν(0)2 as a function of t for ϕ ∈ (0.75 − 0.87). (B) Similar plot for ϕ ∈ (0.89 − 0.93). (C) The long time decay of ⟨Pμν(t)Pμν(0)⟩ is fit to , as shown by the dashed lines. The inset shows the dependence of β on ϕ. (D) The data that is fit using the stretched exponential function (black dashed line) is combined with the short time data (blue solid line), which is fit using the cubic spline function. The resulting fits produces a smooth curve⟨Pμν(t)Pμν(0)⟩combined, as shown in the inset.

Area distribution of the cells:

(A) Simulation snapshot for monodisperse cell system. The number of cells in the two-dimensional periodic box is N = 500. (B) Pair correlation function, g(r), as a function of r. There is clear evidence of order, as reflected in the sharp peaks at regular intervals, which reflects the packing in (A). (C) A schematic picture of polydisperse cell system from the simulations. Color bar on the right shows the scale of radii in μm. There is no discernible order. (D) Distribution of cell area extracted from experiment during morphogenesis of Zebrafish blastoderm (extracted from Fig. S2 (A)) [14]. (E) Same as (D) except, P (Ai), used in a typical simulation. Cell radii vary from 2μm to 15μm.

Structure and relaxation behavior for a binary mixture of cells:

(A) A typical simulation snapshot for binary mixture of cells at ϕ = 0.93. (B) The corresponding pair correlation function, g(r), between all the cells. The vertical dashed line is at the first peak position (rmax). (C) rmax Fs(q, t), with , where rmax is the location of the first peak in the g(r), as function of time at various ϕ values. (D) The logarithm of the relaxation time, τα, as a function of ϕ. Over the entire range of ϕ, the increase in τα is well fit by the Vogel-Fucher-Tamman (VFT) relation. Most importantly, the relaxation time does not saturate, which means the evolving tissue cannot be modeled using a 50:50 binary mixture. (E) Effective shear viscosity as a function of ϕ reflects the behavior of τ as a function of ϕ in (D).

Free area decreases monotonically for the binary mixture of cells:

(A) Mean voronoi cell size, ⟨A⟩, as a function of ϕ for the 50:50 binary system. (B) The free area fraction, ϕfree, as a function of ϕ shows that ϕfree decreases monotonically as ϕ increases.

Measure of ergodicity:

(A) Ergodic measure Ω(t) scaled by the value at t = 0 (Ω(0)) as a function of t for ϕ = 0.85. (B) and (C) Similar plots for ϕ = 0.90 and ϕ = 0.92 respectively. At long time Ω(t)/Ω(0) reach 0.01, 0.016 and 0.026 for ϕ = 0.85, ϕ = 0.90 and ϕ = 0.92 respectively. (D) Ω(0)/Ω(t) as a function of t for ϕ = 0.90. The dashed line shows a linear fit.

Cell size dependent structures and dynamics:

(A) Radial distribution function gSS(r) between small sized cells (RS ≤ 4.5μm) at ϕ = 0.905 (blue) and 0.92 (red). These values are greater than ϕs ≈ 0.90 (B) Same as (A) except the results are for gBB(r) between large cells (RB ≥ 12.0μm). (C) Fs(q, t) for cells with RS ≤ 4.5μm at ϕ = 0.905 and ϕ = 0.92. Note that even at these dense packings, the mobility of the smaller sized cells is substantial, which is reflected in the time dependence of Fs(q, t). (D) Fs(q, t) for cells with RB ≥ 12.0μm at ϕ = 0.905 and ϕ = 0.92. The black dashed lines are fits to stretched exponential functions, where τα is the relaxation time and β is the stretching exponent. The dotted lines correspond to the value .

Simulation snapshot and trajectories for a few smaller and bigger sized cells:

(A) Cells (ϕ = 0.91) are coloured according to their sizes (gray colors). A few small sized cells are shown in different colors (pink, blue, orange, purple, cyan, light purple and yellow). (B) The corresponding trajectories are shown over the entire simulation time. (C) Similar plot as (A) but for a few bigger sized cells shown in purple, yellow, light green, red, cyan and green colors. (D) Same as (B) expect the trajectories of the large sized cells are highlighted. Clearly, the large cells are jammed.

Dynamical rearrangement of jammed cells:

The changing local environment of a randomly selected cell(black) over time. (Top Panels: From left to right t = 9.41τα, 10.01τα and 25.39τα). The black coloured cell is completely jammed by other cells. (Bottom Panels: Form left to right t = 10.97τα, 25.44τα and 27.49τα). Dynamical facilitation, resulting in collective rearrangement of the cells surrounding the black cell, enables it to move in the dynamically created free volume.

Finite size effects:

Fs(q, t) for N = 200 (A) and N = 750 (B). Logarithm of τα as a function of ϕ for N = 200 (C) and for N = 750 (D). The dashed lines are the VFT fits.

Mean coordination number and cell area fraction:

(A), (B) and (C) shows the distribution of coordination number P (Nc) for ϕ =0.85, 0.90 and 0.93 respectively. The orange lines are Gaussian fits to the histograms. (D) Shows mean ⟨Nc⟩ as a function of ϕ. The dashed line shows the linear relationship between them.

Viscosity and coordination number:

(A) Shows ⟨C⟩ as a function of ⟨Nc⟩. Clearly they are linearly related as shown by the dashed line. Viscosity as a function of ⟨Nc⟩ (B) and ⟨C⟩ (C).

Connectivity profile:

Connectivity map for ϕ = 0.80, 0.85, 0.89, 0.90, 0.92 and 0.93 are shown in (A), (B), (C), (D), (E) and (F) respectively. For ϕ ≥ 0.89 there is a path that connects the cells in the entire sample. The percolation transition occurs over a very narrow range of ϕ (roughly at ϕ ≈ 0.89 -orange map), which also coincides with the sharp increase in η, thus linking equilibrium transition to geometric connectivity.