Introduction

There is great interest in characterizing the mechanical and dynamical properties of embryonic tissues because they regulate embryo development [15]. Measurements of bulk properties, such as viscosity and elastic modulus, and the dynamics of individual cells through imaging techniques, have been interpreted by adapting concepts developed to describe phase transitions, glass transition, and active matter [69]. Several experiments have shown that during embryo morphogenesis, material properties of the tissues change dramatically [10– 14]. Of relevance to our study is a remarkable finding that provided evidence that a phase transition (PT) occurs during zebrafish blastoderm morphogenesis, which was analyzed using rigidity percolation theory [1417]. The authors also estimated the viscosity (η) of the blastoderm tissue using the micropipette aspiration technique [13, 18]. It was found that change in η is correlated with cell connectivity (⟨C⟩), rising sharply over a narrow range of ⟨C⟩. Surprisingly, a single geometrical quantity, the cell-cell contact topology controls both the rigidity PT as well as changes in η in this non-confluent tissue, thus linking equilibrium and transport properties.

Here, we focus on two pertinent questions that arise from the experiments on zebrafish blastoderm. First, experiments [14] showed that η increases as a function of the cell packing fraction (ϕ) till ϕ ≤ 0.87. The dependence of η on ϕ follows the well-known Vogel-Fulcher-Tamman (VFT) law [19], which predicts that η grows monotonically with ϕ. The VFT law, which is commonly used to analyze the viscosity of a class of glass forming materials [20], is given by η expwhere ϕ0 is a constant. Surprisingly, for packing fractions, ϕ ≥ ϕS ≈ 0.90, η deviates from the VFT law and is independent of ϕ, which cannot be explained using conventional theories for glasses [8, 21]. Second, the experimental data [14] was interpreted using equilibrium rigidity percolation theory [1517] for an embryonic tissue in which cells undergo random cell divisions. A priori it is unclear why equilibrium concepts should hold in zebrafish morphogenesis, which one would expect is controlled by non-equilibrium processes such as self-propulsion, growth and cell division.

We show that the two conundrums (saturation of η at high packing fractions and the use of equilibrium statistical mechanics in a growing system to explain phase transition) may be rationalized by (i) assuming that the interactions between the cells are soft, (ii) the cell sizes are highly heterogeneous (polydisperse), which is the case in zebrafish blastoderm. Using an agent-based (particle) simulation model of a two-dimensional non-confluent tissue, we explore the consequences of varying ϕ (see section Materials and Methods for the definition) of interacting self-propelled polydisperse soft cells, on η. The central results of our study are: (i) The calculated effective viscosity (for details see Appendix F), for the polydisperse cell system, shows that for ϕ ≤ ϕS ≈ 0.90 the increase in viscosity follows the VFT law. Just as in experiments, η is essentially independent of ϕ at high (≥ ϕS) packing fractions. The unusual dependence of η at ϕ ≥ ϕS is quantitatively explained using the notion of available free area fraction (ϕfree), which is the net void space that can be explored by the cells when they are jammed. At high densities, a given cell requires free space in order to move. The free area is created by movement of the neighboring cells into the available void space. One would intuitively expect that the ϕfree should decrease with increasing packing fractions, due to cell jamming, which should slow down the overall dynamics. Indeed, we find that ϕfree decreases with increasing packing fraction (ϕ) until ϕS. The simulations show that when ϕ exceeds ϕS, the free area ϕfree saturates because the soft cells (characterized by “soft deformable disks”) can overlap with each other, resulting in the collective dynamics of cells becoming independent of ϕ for ϕ ≥ ϕS. As a consequence η saturates at high ϕ. (iii) Cells whose sizes are comparable to the available free area move almost like a particle in a liquid. The motility of small sized cells facilitates adjacent cells to move through multi-cell rearrangement even in a highly jammed environment. The facilitation mechanism, invoked in glassy systems [22] allows large cells to move with low mobility. A cascade of such facilitation processes enable all the cells to remain dynamic even above the onset packing fraction of the phase transition. (iv) We find that the relaxation time does not depend on the waiting time for measurements even in the regime where viscosity saturates. In other words, there is no evidence of aging even in the regime where viscosity saturates. Strikingly, the tissue exhibits ergodic [23] behavior at all densities. The cell-based simulations, which reproduce the salient experimental features, may be described using equilibrium statistical mechanics, thus providing credence to the use of cell contact mechanics to describe both rigidity PT and dynamics in an evolving non-confluent tissue [14].

Results

A. Experimental results

We first describe the experimental observations, which serve as the basis for carrying out the agent-based simulations. Fig. 1(A) shows the bright-field images of distinct stages during zebrafish morphogenesis. A 2D section of zebrafish blastoderm (Fig. 1 (B)) shows that there is considerable dispersion in cell sizes. The statistical properties of the cell sizes are shown in Fig. 8 (D). Fig. 1 (C) shows that η increases sharply over a narrow ϕ range, and saturates when ϕ exceeds ϕS ≈ 0.90.

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

To account for the results in Fig. 1 (C), we first simulated a mono-disperse system in which all the cells have identical radius (R = 8.5 μm). Because the system crystallizes (Fig. 8 (A) & (B)), we concluded that the dynamics observed in experiments cannot be explained using this model. A 1:1 binary mixture of cells, with different radii gives glass-like behavior for all ϕ, with the relaxation time τα as well as the effective viscosity (defined in Eq. 1) following the VFT behavior (see Appendix B).

B. Polydispersity and cell-cell interactions

In typical cell tissues, and zebrafish in particular, there is a great dispersion in the cell sizes, which vary in a single tissue by a factor of ∼ 5 − 6 [14] (Fig. 1 (B) and Fig. 8 (D)). In addition, the elastic forces characterizing cell-cell interactions are soft, which implies that the cells can overlap, with rij −(Ri +Rj) < 0 when they are jammed ((Fig. 1 (B) and (Fig. 1 (D)). Thus, both polydispersity as well as soft interactions between the cells must control the relaxation dynamics. To test this proposition, we simulated a highly polydisperse system (PDs) in which the cell sizes vary by a factor of ∼ 8 (Fig. 1 (D) and Fig. 8 (E)).

A simulation snapshot (Fig. 1 (D)) for ϕ = 0.93 shows that different sized cells are well mixed. In other words, the cells do not phase separate. The structure of the tissue can be described using the pair correlation function, , where is the number density, δ is the Dirac delta function, is the position of the ith cell, and the angular bracket ⟨⟩ denotes an average over different ensembles. The g(r) function (Fig. 1 (E)) has a peak around r ∼ 17μm, which is approximately the average diameter of the cells. The absence of peaks in g(r) beyond the second one, suggests there is no long

C. Effective shear viscosity() as a function of ϕ

A fit of the experimental data for η using the VFT [24, 25] relation in the range ϕ ≤ 0.87 (Fig. 1 (C)) yields ϕ0 ≈ 0.95 and D ≈ 0.51 [19]. The VFT equation for cells, which is related to the Doolittle equation [26] for fluidity that is based on free space available for motion in an amorphous system [27, 28], is η = η exp, where D is the apparent activation energy. In order to compare with experiments, we calculated an effective shear viscosity() for the polydisperse system using a Green-Kubo-type relation [29],

The stress tensor Pμν(t) in the above equation is,

where μ, ν ∈ (x, y) are the Cartesian components of coordinates, is the force between ith and jth cells and A is the area of the simulation box. Note that should be viewed as a proxy for shear viscosity because it does not contain the kinetic term and the factor is not included in Eqn. (1) because temperature is not a relevant variable in the highly over-damped model for cells considered here.

Plot of as a function of ϕ in Fig. 2 (A), shows qualitatively the same behavior as the estimate of viscosity (using dimensional arguments) made in experiments. Two features about Fig. 1 (C) and Fig. 2 (A) are worth noting. (i) Both simulations and experiments show that up to ϕ ≈ 0.90, follows the VFT relation with ϕ0 ∼ 0.94 and D ∼ 0.5. More importantly is independent of ϕ when ϕ > 0.90. (ii) The values of ϕ0 and D obtained by fitting the experimental estimate of η to the VFT equation and simulation results are almost identical. Moreover, the onset of the plateau packing fraction in simulations and experiments occur at the same value (ϕS ∼ 0.90). The overall agreement with experiments is remarkable given that the model was not created to mimic the zebrafish tissue.

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.

To provide additional insights into the dynamics, we calculated the isotropic self-intermediate scattering function, Fs(q, t),

Where is the wave vector, and is the position of a cell at time t. The degree of dynamic correlation between two cells can be inferred from the decay of Fs(q, t). The angle bracket ⟨…⟩ is an average over different time origins and different trajectories. We chose where rmax is the position of the first peak in g(r) between all cells (see Fig. 1 (E)). The relaxation time τα is calculated using .

From Figs. 2 (B) & (C), which show Fs(q, t) as a function of t for various ϕ, it is clear that the dynamics become sluggish as ϕ increases. The relaxation profiles exhibit a two step decay with a plateau in the intermediate time scales. The dynamics continues to slow down dramatically until ϕ ≤ 0.90. Surprisingly, the increase in the duration of the plateau in Fs(q, t) ceases when ϕ exceeds ≈ 0.90 (Fig. 2 (C)), a puzzling finding that is also reflected in the dependence of τα on ϕ in Fig. 2 (D). The relaxation time increases dramatically, following the VFT relation, till ϕ ≈ 0.90, and subsequently reaches a plateau (see the inset in Fig. 2 (D)). If the VFT relation continued to hold for all ϕ, as in glasses or in binary mixture of 2D cells (see Appendix B), then the fit yields ϕ0 ≈ 0.95 and D ≈ 0.50. However, the simulations show that τα is nearly a constant when ϕ exceeds 0.90. We should note that the behavior in Fig. 2 (D) differs from the dependence of τα on ϕ for 2D monodisperse polymer rings, used as a model for soft colloids. Simulations [30] showed τα increases till a critical ϕS but it decreases substantially beyond ϕS with no saturation.

D. Relaxation dynamics of individual cells

Plot of τα as a function of the radius of cells Ri (Fig. 3 (A)) shows nearly eight orders of magnitude change. The size dependence of τα on ϕ is striking. That τα should increase for large sized cells (see the data beyond the vertical dashed line in Fig. 3 (A)) is not unexpected. However, even when cell sizes increase beyond Ri = 4.25 μm, the dispersion in τα is substantial, especially when ϕ exceeds ϕS. The relaxation times for cells with Ri < 4.25 μm are relatively short even though the system as a whole is jammed. For ϕ ≥ 0.90, τα for small sized cells have a weak dependence on ϕ. Although τα for cells with radius < 4 μm is short, it is clear that for a given ϕ (for example ϕ = 0.93) the variations in τα is substantial. In contrast, τα’s for larger cells (R ≥ 7μm) are substantially large, possibly exceeding the typical cell division time in experiments. In what follows, we interpret these results in terms of available free area ⟨Afree⟩ for cells. The smaller sized cells have the largest (RS is the radius of the small cell).

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.

The effect of jamming on the dramatic increase in τα occurs near Ri ≈ 4.5 μm, which is comparable to the length scale of short range interactions. For ϕ ≤ 0.90, τα increases as the cell size increases. However, at higher packing fractions, even cells of similar sizes show substantial variations in τα, which change by almost 3 − 4 orders of magnitude (see the data around the vertical dashed line for ϕ ≥ 0.915 in Fig. 3 (A)). This is a consequence of large variations in the local density (Fig. 14). Some of the similar-sized cells are trapped in the jammed environment, whereas others are in less crowded regions (see Fig. 14). The spread in τα increases dramatically for ϕ > ϕS (≈ 0.90) and effectively overlap with each other. This is vividly illustrated in the the histogram, P(log(τα)) shown in Fig. 3 (B). For ϕ < ϕs, the peak in P(log(τα)) monotonically shifts to higher log(τα) values. In contrast, when ϕ exceeds ϕS there is overlap in P(log(τα)), which is reflected in the saturation of and τα. There are cells (typically with small sizes) that move faster even in a highly jammed environment (see Fig. 12 (C) and Fig. 13). The motions of the fast-moving cells change the local environment, which effectively facilitates the bigger cells to move in a crowded environment (see Fig. 12 (D) and Fig. 13 and and Movie 1 (ϕ = 0.92 > ϕS) and Movie 2 (ϕ = 0.90 = ϕS)). In contrast, for ϕ = 0.85 < ϕS small and large sized cells move without hindrance because of adequate availability of free area (Movie 3). The videos vividly illustrate the large scale facilitated rearrangements that enable the large sized cells to move.

The dependence of τα on ϕ for ϕ ≤ ϕS (Fig. 2 (D)) implies that the polydisperse cell systems behave as a soft glass in this regime. On theoretical grounds, it was predicted that in glass-forming systems [8]. Remarkably, we found that this prediction is valid in the polydisperse cell system (Fig. 3 (C)). However, above ϕS the predicted relation is not satisfied (see Fig. 3 (D)).

E. Available free area explains viscosity saturation at high ϕ

We explain the saturation in the viscosity by calculating the available free area per cell, as ϕ increases. In a hard disk system, one would expect that the free area would decrease monotonically with ϕ until it is fully jammed at the close packing fraction (∼ 0.84 [31, 32]). Because the cells are modeled as soft deformable disks, they could overlap with each other even when fully jammed. Therefore, the region where cells overlap create free area in the immediate neighborhood. The extent of overlap (hij) is reflected in distribution P(hij). The width in P(hij) increases with ϕ, and the peak shifts to higher values of hij (Fig. 4 (A)). The mean, ⟨hij⟩ increases with ϕ (Fig. 4 (B)). Thus, even if the cells are highly jammed at ϕ ≈ ϕS, free area is available because of an increase in the overlap between cells (see Fig. 5). When ϕ exceeds ϕS, the mobility of small sized cells facilitates the larger cells to move, as is assumed in the free volume theory of polymer glasses [28, 3335]. As a result of the motion of small cells, a void is temporarily created, which allows other (possibly large) cells to move. In addition to the release of space, the cells can also interpenetrate (Figs. 4 (A) & (B)). If hij increases, as is the case when the extent of compression increases (Figs. 4 (A)), the available space for nearby cells would also increase. This effect is expected to occur with high probability at ϕS and beyond, resulting in high overlap between the cells. These arguments suggest that the combined effect of polydispersity and cell-cell overlap creates, via the self-propulsion of cells, additional free area that drives larger cells to move even under jammed conditions.

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

In order to quantify the physical picture given above, we calculated an effective area for each cell by first calculating Voronoi cell area A. A plot for Voronoi tessellation is presented in Fig. 5 (A) for ϕ = 0.93, and the histogram of A is shown in Fig. 5 (B). As ϕ increases, the distribution shifts towards lower Voronoi cell size ⟨A⟩. The mean Voronoi cell size ⟨A⟩ as a function of ϕ in Fig. 5 (C) shows ⟨A⟩ decreases as ϕ is increased. As cells interpenetrate, the Voronoi cell size will be smaller than the actual cell size in many instances (Fig. 5 (A)). To demonstrate this quantitatively, we calculated . The value of Afree could be negative if the overlap between neighbouring cells is substantial; Afree is positive only when the Voronoi cell size is greater than the actual cell size. Positive Afree is an estimate of the available free area. The histograms of Afree for all the packing fractions in Fig. 5 (D) show that the distributions saturate beyond ϕ = 0.90. All the distributions have a substantial region in which Afree is negative. The negative value of Afree increases with increasing ϕ, which implies that the amount of interpenetration between cells increases.

Because of the overlap between the cells, the available free area fraction ϕfree, is higher than the expected free area fraction (1.0 − ϕ) for all ϕ. We define an effective free area fraction ϕfree as,

where Np is the number of positive free area in jth snapshots, Nt is the total number of snapshots, Abox is the simulation box area and is the positive free area of ith cell in jth snapshot. The calculated ϕfree, plotted as a function of ϕ in Fig. 5 (E), shows that ϕfree decreases with ϕ until ϕ = 0.90, and then it saturates near a value ϕfree ≈ 0.22 (see the right panel in Fig. 5 (E)). Thus, the saturation in as a function of ϕ is explained by the free area picture, which arises due to combined effect of the size variations and the ability of cells to overlap.

F. Aging does not explain viscosity saturation

Our main result, which we explain by adopting the free volume theory developed in the context of glasses [28, 3335], is that above a critical packing fraction ϕS ∼ 0.90 the viscosity saturates. Relaxation time, τα, measured using dynamic light scattering, in nearly monodisperse microgel poly(N-isopropylacrylamide) (PNiPAM) [36] was found to depend only weakly on the volume fraction (3D), if ϕV, exceeds a critical value. It was suggested that the near saturation of τα at high ϕV is due to aging, which is a non-equilibrium effect. If saturation in viscosity and relaxation time in the embryonic tissue at high ϕ is due to aging then τα should increase sharply as the waiting time, τω, is lengthened. We wondered if aging could explain the observed saturation of η in the embryonic tissue above ϕS. If aging causes the plateau in the tissue dynamics, then η or τα should be an increasing function of the waiting time, τω. To test the effect of τω on τα, we calculated the self-intermediate scattering function Fs(q, t+τω) as a function of t by varying τω over three orders of magnitude at ϕ = 0.92 (Fig. 6 (A)). There is literally no change in Fs(q, t + τω) over the entire range of tω. We conclude that, τα, extracted from Fs(q, t+τω) is independent of τω. The variations in τα (Fig. 6 (B)), with respect to τω, is significantly smaller than the errors in the simulation. Thus, the saturation in η or τα when ϕ > ϕS is not a consequence of aging.

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.

There are two implications related to the absence of aging in the dynamics of the non-confluent embryonic tissues. (i) Although active forces drive the dynamics of the cells, as they presumably do in reality, the cell collectives can be treated as being near equilibrium, justifying the use of Green-Kubo relation to calculate η. (ii) Parenthetically, we note that the absence of significant non-equilibrium effects, even though Zebrafish is a living system, further justifies the use of equilibrium rigidity percolation theory to analyze the experimental data [14].

Discussion

Extensive computer simulations of a two-dimensional dense tissue using a particle-based model of soft deformable cells with active self-propulsion have successfully reproduced the dynamical behavior observed in the blastoderm tissue of zebrafish. The dependence of viscosity (η) and relaxation time (τα) (before the saturation) is well fit by the VFT equation. The value of ϕ0 obtained from simulations, ϕ0 ∼ 0.95, is close to ϕ0 ∼ 0.94 extracted by fitting the experimental data to the VFT equation. Thus, the dynamics for ϕ ≤ ϕS resembles the behavior expected for glass forming systems. Remarkably, the dependence of η on ϕ over the entire range (VFT regime followed by a plateau) may be understood using available free area picture with essentially a single parameter, an idea that was proposed nearly 70 years ago. We discovered that polydispersity as well as the ease of deformation of the cells that creates free area under high jamming conditions, is the mechanism that explains viscosity saturation at high cell densities. The mechanism suggested here is an important step that links equilibrium phase transition to dynamics during zebrafish development [37].

One could legitimately wonder if the extent of polydispersity (PD) used in the soft discs simulations, which seems substantial, is needed to recapitulate the observed dependence of η on ϕ. Furthermore, such large values of PD may not represent biological tissues. Although the choice of PD was made in part by the two-dimensional projection of area reported in experiments [14], it is expected that PD values have to be less in three dimensions. We performed preliminary simulations in three dimensions with considerably reduced PD, and calculated the dependence of relaxation time (τα) as a function of ϕ. The results show that τα does indeed saturate at high volume fractions.

The proposed model neglects adhesive interactions between cells, which of course is not unimportant. It is crucial to wonder if the proposed mechanism would change if cell-cell adhesion is taken into account. We wanted to create the simplest model to explain the ex-perimental data. We do think that realistic values of adhesion strength does not significantly alter the forces between cells [38]. Thus, we expect a similar mechanism. Furthermore, the physics of the dynamics in glass forming materials does not change in systems with and without attractive forces [8]. Universal behavior, such as Vogel-Fucher-Tammann relation, is valid for a broad class of unrelated materials (see Fig.1 in [20]). Needless to say, non-universal quantities such as glass transition temperature Tg or effective free energy barriers for relaxation will change. In our case, we expect that changing the adhesion strength, within a reasonable range, would change ϕS without qualitatively altering the dependence of η on ϕ. For these reasons, in the first pass we neglected adhesion, whose effects have to be investigated in the future.

In the physical considerations leading to Eq. 6, the random activity term (μ) plays an important role. Is it possible to create a passive model by maintaining the system at a finite temperature using stochastic noise with μ = 0, which would show the observed viscosity behavior? First, in such a system of stochastic equations, the coefficient of noise (a diffusion constant) would be related to γi in Eq.6 through fluctuation dissipation theorem (FDT). Thus, only γi can be varied. In contrast, in Eq.6 the two parameters (γi and μ) maybe independently changed, which implies that the two sets of stochastic equations of motion are not equivalent. Second, the passive system describes particles that interact by soft Hertz potential. In analogy with systems in which the particles interact with harmonic potential [39], we expect that the passive model would form a glass in which the viscosity would follow the VFT law.

We find it surprising that the calculation of viscosity using linear response theory (valid for systems close to equilibrium), and the link to free area quantitatively explain the simulation results and by implication the experimental data for a living and growing tissue. The calculation of free area of the cells is based on the geometrical effects of packing, which in turn is determined by cell to cell contact topology. These considerations, that are firmly established here, explain why equilibrium phase transitions are related to steep increase in viscosity [8] as the packing fraction changes over a narrow range. The absence of aging suggests that, although a large number of cell divisions occur, they must be essentially independent, thus allowing the cells to reach local equilibrium.

Acknowledgements

We acknowledge Anne D. Bowen at the Visualization Laboratory (Vislab), Texas Advanced Computing Center, for help with video visualizations. We are grateful to Mauro Mugnai and Nicoletta Petridou for useful discussions. This work was supported by grants from the National Science Foundation (PHY 23-10639 and CHE 23-20256, and the Welch Foundation (F-0019).

Materials and Methods

Two dimensional cell model

Following our earlier studies [38, 40], we simulated a two dimensional (2D) version of a particle-based cell model. We did not explicitly include cell division in the simulations. This is physically reasonable because in the experiments [14] the time scales over which cell division induced local stresses relax are short compared to cell division time. Thus, local equilibrium is established in between random cell division events. We performed simulations in 2D because experiments reported the dependence of viscosity as a function of area fraction.

In our model, cells are modeled as soft deformable disks [38, 4143] interacting via short ranged forces. The elastic (repulsive) force between two cells with radii Ri and Rj is Hertzian, which is given by,

where . The repulsive force acts along the unit vector , which points from the center of the jth cell to the centre of ith cell. The total force on the ith cell is

where NN(i) is the number of near neighbor cells that is in contact with the ith cell. The jth cell is the nearest neighbor of the ith cell, if hij > 0. The near neighbor condition ensures that the cells interpenetrate each other to some extent, thus mimicking the cell softness. For simplicity, we assume that the elastic moduli (E) and the Poisson ratios (ν) for all the cells are identical. Polydispersity in the cell sizes is important in recovering the plateau in the viscosity as a function of packing fraction. Thus, the distribution of cell areas (Ai = πR2) is assumed to have a distribution that mimics the broad area distribution discovered in experiments.

Self-propulsion and equations of motion

In addition to the repulsive Hertz force, we include an active force arising from self-propulsion mobility (μ), which is a proxy for the intrinsically generated forces within a cell. For illustration purposes, we take μ to independent of the cells, although this can be relaxed readily. We assume that the dynamics of each cell obeys the phenomenological equation,

where γi is the friction coefficient of ith cell, and 𝒲i(t) is a noise term. The friction coefficient γi is taken to be γ0Ri [44]. By scaling t by the characteristic time scale, in Eqn. (6), one can show that the results should be insensitive to the exact value of μ. The noise term 𝒲i(t) is chosen such that ⟨𝒲i(t)⟩ = 0 and . In our model there is no dynamics with only systematic forces because the temperature is zero. The observed dynamics arises solely due to the self-propulsion (Eqn. (6)).

We place N cells in a square box that is periodically replicated. The size of the box is L so that the packing fraction (in our two dimensional system it is the area fraction) is . We performed extensive simulations by varying ϕ in the range 0.700 ≤ ϕ ≤ 0.950. The results reported in main text are obtained with N = 500. Finite size effects are discussed in Appendix G.

To mimic the variations in the area of cells in a tissue [14], we use a broad distribution of cell radii (see Appendix A for details). The parameters for the model are given in Table I. In the present study, we do not consider the growth and division of cells. Thus, our simulations describe steady-state dynamics of the tissue. For each ϕ, we performed simulations for at least (5−10)τα before storing the data. For each ϕ, we performed 24 independent simulations. The calculation of viscosity was performed by averaging over 40 independent simulations at each ϕ.

Parameters used in the simulation.

Calculation of viscosity

We calculated the effective viscosity for various values of ϕ by integrating the off-diagonal part of the stress-stress correlation function ⟨Pμν(t)Pμν(0)⟩, using the Green-Kubo relation [46] (without the pre-factor),

where μ and ν denote Cartesian components (x and y) of the stress tensor Pμν(t) (see main text for the definition of Pμν(t)). The definition of , which relates the decay of stresses as a function of times in the non-confluent tissue, is akin to the methods used to calculate viscosity in simple fluids. (Eq.7) The time dependence of ⟨Pμν(t)Pμν(0)⟩, normalized by ⟨Pμν(0)2⟩, for different values of ϕ (Fig. 7 (A) & (B)) shows that the stress relaxation is clearly non-exponential, decaying to zero in two steps. After an initial rapid decay followed by a plateau at intermediate times (clearly visible for ϕ ≥ 0.91) the normalized ⟨Pμν(t)Pμν(0)⟩ decays to zero as a stretched exponential. The black dashed lines in Fig. 7 (C)) show that a stretched exponential function, , where τ is the characteristic time in which stress relax and β is the stretching exponent, provides an excellent fit to the long time decay of ⟨Pμν(t)Pμν(0)⟩ (from theI plateau1region to zero) as a function of t. Therefore, we utilized the fit function, , to replace the noisy long time part of ⟨Pμν(t)Pμν(0)⟩ by a smooth fit data before evaluating the integral in Eqn. (7). The details of the procedure to compute is described below.

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.

We divided ⟨Pμν(t)Pμν(0)⟩ in two parts. (i) The short time part (⟨Pμν(t)Pμν(0)⟩short) – the smooth initial rapid decay until the plateau is reached (for example see the blue circles in Fig. 7(D) for ϕ = 0.93). For the n data points at short times, t1, ⟨Pμν(t1)Pμν(0)⟩short,…,(tn, ⟨Pμν(tn)Pμν(0)⟩short), we constructed a spline S(t) using a set of cubic polynomials:

The polynomials satisfy the following properties. (a) Si(ti) = ⟨Pμν(ti)Pμν(0)⟩short and Si(ti+1) = ⟨Pμν(ti+1)Pμν(0)⟩short for i = 1, .., n − 1 which guarantees that the spline function S(t) interpolates between the data points. for i = 2, .., n − 1 so that S(t) is continuous in the interval [t1, tn]. so that S′′(t) is continuous in the interval [t1, tn]. By solving for the unknown parameters, bi, ci and di using the above mentioned properties, we constructed the function S(t). We used S(t) to fit ⟨Pμν(t)Pμν(0)⟩short to get an evenly spaced (δt = 10s) smooth data (solid blue line in Fig. 7 (D)). The fitting was done using the software “Xmgrace”.

(ii) The long time part (⟨Pμν(t)Pμν(0)⟩long) – from the plateau until it decays to zero-is shown by the red circles in Fig. 7 (D). The long time part was fit using the analytical function (black dashed line in Fig. 7 (D)). We refer to the fit data (δt = 10s) as . We then combined ⟨Pμν(t)Pμν(0)⟩short and to obtain ⟨Pμν(t)Pμν(0)⟩combined (see inset of Fig. 7 (D)). Finally, we calculated using the equation,

where t1δt is the end point of ⟨Pμν(t)Pμν(0)⟩short and T δt is the end point of ⟨Pμν(iδt)Pμν(0)⟩combined.

Appendix A Cell polydispersity is needed to account for viscosity saturation

To explain the observed saturation in viscosity (Fig. 1 (C)) when cell area fraction, ϕ, exceeds ϕS (≈ 0.90), we first simulated a monodisperse cell system using the model described in the Method section. The monodisperse cell system (R = 8.5μm) crystallizes (see Fig. 8 (A) and (B)), which excludes it from being a viable model for explaining the experimental findings. We also find that a system consisting of 50:50 binary mixture of cells can not account for the experimental data even though crystallization is avoided (see the next section). These findings forced us to take polydispersity into account.

In order to develop an agent based model that accounts for polydispersity effects, we first extracted the distribution, P(Ai), of the area (Ai) in the Zebrafish cells from the experimental data [14]. Figure 8 (D), shows that P(Ai) is broad, implying that cell sizes are highly heterogeneous. Based on this finding, we sought a model of the non-confluent tissue that approximately mimics polydispersity found in experiments. In other words, (ΔA is the dispersion in P(Ai) and ⟨A⟩ is the mean value) should be similar to the data in Fig. 8 (D). The radii (Ri) of the cells in the simulations are sampled from a Gaussian distribution ∼ exp (−(Ri − ⟨R⟩)2/2ΔR2) with ⟨R⟩ = 8.5μm and ΔR = 4.5μm. The resulting P(Ai) for one of the realizations (see Fig. 8 (C)) is shown in Fig. 8 (E). The value of is 0.5, which compares favourably with the experimental estimate (∼ 0.4). The unusual dependence of the viscosity (η) as a function of packing fraction (ϕ) cannot be reproduced in the absence of polydispersity.

Appendix B: Relaxation time in the binary mixture of cells does not saturate at high ϕ

We showed that the saturation in the relaxation time above a critical area fraction, ϕS, is related to two factors. One is that there ought to be dispersion in the cell sizes (Fig. 8). The extent of dispersion is likely to less in three rather in two dimensions. The second criterion is that the cells should be soft, allowing them to overlap at high area fraction. In other words, the cell diameters are explicitly non-additive. The Hertz potential captures the squishy nature of the cells.

In order to reveal the importance of polydispersity, we first investigated if a binary mixture of cells (see Fig. 9 (A)) would reproduce the observed dependence of τα on ϕ. We created a 50 : 50 mixture with a cell size ratio ∼ 1.4 (RB = 8.5μm and RS = 6.1μm). All other parameters are the same as in the polydisperse system (see Table I). The radial distribution function, g(r), calculated by considering both the cell types, exhibits intermediate but not long range translational order (Fig. 9 (B)). A peak in g(r) appears at rmax = 14.0μ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).

In order to determine the dependence of the relaxation time, τα, on ϕ, we first calculated Fs(q, t) at q = 2π/rmax (Fig. 9 (C)). The decay of Fs(q, t) is similar to what one finds in typical glass forming systems. As ϕ increases, the decay of Fs(q, t) slows down dramatically. When ϕ exceeds 0.93 there is a visible plateau at intermediate times followed by a slow decay. By fitting the long time decay of Fs(q, t) to exp(−(t/τα)β) (β is the stretching exponent), we find that τα(ϕ) as a function of ϕ is well characterized by the VFT relation (Fig. 9 (D)). There is no evidence of saturation in τα(ϕ) at high ϕ. The ϕ dependence of the effective shear viscosity (Fig. 9 (E)), calculated using equation (8), also shows no sign of saturation. The VFT behavior, with ϕ0 ≈ 1 and D ≈ 1.2, shows that the two component cell system behaves as a “fragile” glass [20]. A comment regarding the binary system of cells is in order. The variables that characterize this system are, λ = RB/RS (RB (RS) is the radius of the big (small) cells), ΦB = NB/(NA + NB) (NB (NA) is the number of big (small) cells) with ΦA = 1 − ΦB, and the packing fraction, AS where AS is the area of the sample. The results in Fig. 9 were obtained using λ = 1.4, ΦB = 0.5. With this choice the value of polydispersity is (λ − 1)/(λ + 1) ≈ 0.17, which is smaller than the experimental value. It is possible that by thoroughly exploring the parameter space, λ and ϕ, one could find regions in which the τα in the two component cell system would saturate beyond ϕS.

However, in light of the results in Fig. 8 (D)), we choose to simulate systems that also have high degree of polydispersity (Fig. 8 (E)).

Appendix C Absence of saturation in the free area in binary mixture of cells

In the main text, we established that the effective viscosity, , which should be a proxy for the true η and the relaxation time (τα) in the polydisperse cell system saturates beyond ϕS. The dynamics saturates beyond ϕS because the available area per cells (quantified as ϕfree) is roughly a constant in this region (see Fig. 5 (E) in main text). In the binary system, however, we do not find any saturation in the τα. Therefore, if ϕfree controls the dynamics, one would expect that ϕfree should decrease monotonically with ϕ for the binary cells system. We calculated the average Vornoi cell size ⟨A⟩ (Fig. 10 (A)) and ϕfree (Fig. 10 (B)) for the binary system. Both ⟨A⟩ and ϕfree decrease with ϕ, which is consistent with the free volume picture proposed in the context of glasses.

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.

Appendix D Absence of broken ergodicity

We have shown in the earlier section that, the saturation of τα above ϕS is not a consequence of aging. In the context of glasses and supercooled liquids [23], it has been shown that if ergodicity is broken then a variety of observables would depend on initial conditions. In order to test if ergodicity is broken in our model of non-confluent tissues, we define a measure Ω(t) given by,

where, k and l represent systems with different initial conditions, and with Pμν(s) being the value of stress (see Eq. (2) in the main text) at time s. If the system is ergodic, implying the system has explored the whole phase space on the simulation time scale, the values of ωk(t) and ωl(t) would be independent of k and l. Therefore, Ω(t) should vanish at very long time for ergodic systems. On the other hand, if ergodicity is broken then Ω(t) would be a constant whose value would depend on k and l.

The long time values of Ω(t), normalized by Ω(0) (at t = 0), are 0.01, 0.016 and 0.026 for ϕ = 0.85, 0.90 and 0.92 respectively (Fig. 11 (A), (B) & (C)). Because these values are sufficiently small, we surmise that effectively ergodicity is established. Therefore, our conclusion in the earlier section that the polydisperse cell system is in near equilibrium is justified, and also explains the absence of aging. Furthermore, it was also predicted previously [23] that in long time the ergodic measure (Ω(t) in our case) should decay as ≈ 1/t. Fig. 11 (D) shows that this is indeed the case -at long time Ω(t)/Ω(0) decays approximately as 1/t.

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.

Appendix E Dynamics of small and large cells are dramatically different

The structural and the dynamical behavior of the small (RS ≤ 4.5μm) and large (RB ≥ 12.0μm) cells are dramatically different in the non-confluent tissue. The pair correlation functions between small cells (gSS(r)) and between large cells (gBB(r)) (Fig. 12 (A) & (B)) for ϕ = 0.905 and ϕ = 0.92 show that, both small and large cells exhibit liquid like disordered structures. However, it is important to note that, gSS(r) has only one prominent peak and a modest second peak. In contrast, gBB(r) has three prominent peaks. Thus, the smaller sized cells exhibit liquid-like behavior whereas the large cells are jammed. This structural feature is reflected in the decay of Fs(q, t) with , where rmax is the position of the first peak in the g(r) (Fig. 12 (C) & (D)). There is a clear difference in the decay of Fs(q, t), spanning nearly four orders of magnitude, between small and large cells (compare Figs. 12 (C) and (D)). At the highly jammed packing fraction (ϕ = 0.92), small sized cells have the characteristics of fluid-like behavior.

A typical snapshot from one of the simulations (ϕ = 0.91) and the trajectories for a few small and big sized cells are displayed in Fig. 13 (A), (B), (C) and (D). The figures reflect the decay in Fs(q, t). Not only are the mobilities heterogeneous it is also clear that the the Mean Square Displacements (MSDs) of the small cells are greater than the large cells.

Appendix F: Dynamical changes in local packing fraction cause jammed cells to move

The mobilities of all the cells, even under highly jammed conditions (ϕ ≥ ϕS), is only possible if the local area fraction changes dynamically. This is a collective effect, which is difficult to quantify because it would involve multi-cell correlation function. As explained in the main text, the movement of a jammed cell can only occur if several neighboring cells move to create space. This picture is not that dissimilar to kinetic facilitation in glass-forming materials [47, 48]. However in glass like systems facilitation is due to thermal excitation but in the active system, it is self-propulsion that causes the cells to move.

The creation of free space may be visualized by tracking the positions of the nearest neighbour cells. In Fig. 14, we display the local free area of a black coloured cell at different times. The top panels show the configurations where the black cell is completely jammed by other cells. The cell, coloured in black, can move if the neighboring cells rearrange (caging effect in glass forming systems) in order to increase the available free space. The bottom panels show that upon rearrangement of the cells surrounding the black cell its mobility increases. Such rearrangement occurs continuously, which qualitatively explains the saturation in viscosity in the multicomponent cell system.

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.

Appendix G: Finite system size effects

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.

Appendix H Dependence of viscosity on average coordination number and average connectivity

In the Zebrafish blastoderm experiment [14], the change in viscosity (η) was shown as a function of mean connectivity ⟨C⟩. To test whether our two dimensional tissue simulations also exhibit similar behavior, we calculated viscosity as a function of mean coordination number ⟨Nc⟩ (defined below), which is equivalent to ⟨C⟩. We define the coordination number, Nc, as the number of cells that are in contact with a given cell. Two cells with indices i and j are in contact if hij = Ri + Rj − rij > 0. We calculate Nc for all the cells for each ϕ and calculate the histogram P(Nc). The distributions of P(Nc)s are well fit by a Gaussian distribution function A (Figs. 16 (A), (B) & (C)). The calculated mean, ⟨Nc⟩, from the fit is linearly related to the cell area fraction ϕ (Fig. 16 (D)). We also calculate the average connectivity ⟨C⟩ defined in the following way. Each cell is defined as a node and an edge is defined as the line connecting two nodes. If a snapshot has n nodes and m edges then the connectivity is defined as [14]. We calculate C for all the snapshots for each ϕ and estimated its mean value ⟨C⟩. We find that ⟨C⟩ and ⟨Nc⟩ are of similar values (Fig. 17 (A)), and the dependence of viscosity on ⟨Nc⟩ and ⟨C⟩ is similar to experimental results (Fig. 17 (B) & (C)) [14].

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

Appendix I Connectivity Map

To pictorially observe the percolation transition in the simulations, we plot the connectivity map for various values ϕ in Fig. 18. Note that for smaller ϕ ≤ 0.89 the map shows that cells are loosely connected, suggesting a fluid-like behavior. For ϕ ≥ 0.89 the connectivity between the cells spans the entire system, which was noted and analysed elsewhere thoroughly [14]. Cells at one side of simulation box are connected to cells at the side. The cell connectivity extends throughout the sample. The transition from a non-percolated state to a percolated state occurs over a very narrow range of ϕ, which corresponds to the onset of rigidity percolation transition occurs [14]. It is gratifying that the simple simulations reproduce the experimental observations.

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.