Free volume theory explains the unusual behavior of viscosity in a nonconfluent tissue during morphogenesis
Abstract
A recent experiment on zebrafish blastoderm morphogenesis showed that the viscosity (η) of a nonconfluent embryonic tissue grows sharply until a critical cell packing fraction (ϕ_{S}). The increase in η up to ϕ_{S} is similar to the behavior observed in several glassforming materials, which suggests that the cell dynamics is sluggish or glasslike. Surprisingly, η is a constant above ϕ_{S}. To determine the mechanism of this unusual dependence of η on ϕ, we performed extensive simulations using an agentbased model of a dense nonconfluent twodimensional tissue. We show that polydispersity in the cell size, and the propensity of the cells to deform, results in the saturation of the available free area per cell beyond a critical packing fraction. Saturation in the free space not only explains the viscosity plateau above ϕ_{S} but also provides a relationship between equilibrium geometrical packing to the dramatic increase in the relaxation dynamics.
eLife assessment
This fundamental study substantially advances our physical understanding of the sharp increase and saturation of the viscosity of nonconfluent tissues with increasing cell density. Through the analysis of a simplified model, this study provides compelling evidence that polydispersity in cell size and the softness of cells together can lead to this phenomenon. The work will be of general interest to biologists and biophysicists working on development.
https://doi.org/10.7554/eLife.87966.4.sa0Introduction
There is great interest in characterizing the mechanical and dynamical properties of embryonic tissues because they regulate embryo development (Kimmel et al., 1995; Keller et al., 2008; Petridou and Heisenberg, 2019; Hannezo and Heisenberg, 2019; Autorino and Petridou, 2022). 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 (PTs), glass transition, and active matter (Shaebani et al., 2020; Marchetti et al., 2013; Kirkpatrick and Thirumalai, 2015; Bär et al., 2020).
Several experiments have shown that during embryo morphogenesis, material properties of the tissues change dramatically (Morita et al., 2017; Mongera et al., 2018; Barriga et al., 2018; Petridou et al., 2019; Petridou et al., 2021). Of relevance to our study is a remarkable finding that provided evidence that a PT occurs during zebrafish blastoderm morphogenesis, which was analyzed using rigidity percolation theory (Petridou et al., 2021; Jacobs and Thorpe, 1995; Jacobs and Thorpe, 1996; Jacobs and Hendrickson, 1997). The authors also estimated the viscosity ($\eta$) of the blastoderm tissue using the micropipette aspiration technique (Guevorkian et al., 2010; Petridou et al., 2019). It was found that change in $\eta$ is correlated with cell connectivity ($\u27e8C\u27e9$), rising sharply over a narrow range of $\u27e8C\u27e9$. Surprisingly, a single geometrical quantity, the cell–cell contact topology controls both the rigidity PT and changes in $\eta$ in this nonconfluent tissue, thus linking equilibrium and transport properties.
Here, we focus on two pertinent questions that arise from the experiments on zebrafish blastoderm. First, experiments (Sinha and Thirumalai, 2021) showed that $\eta$ increases as a function of the cell packing fraction ($\varphi$) till $\varphi \le 0.87$. The dependence of $\eta$ on $\varphi$ follows the wellknown Vogel–Fulcher–Tammann (VFT) law (Sinha and Thirumalai, 2021), which predicts that $\eta$ grows monotonically with $\varphi$. The VFT law, which is commonly used to analyze the viscosity of a class of glassforming materials (Angell, 1991), is given by $\eta \sim \mathrm{exp}\left[\frac{1}{{\varphi}_{0}/\varphi 1}\right]$, where $\varphi}_{0$ is a constant. Surprisingly, for packing fractions, $\varphi \ge {\varphi}_{S}\approx 0.90$, $\eta$ deviates from the VFT law and is independent of $\varphi$, which cannot be explained using conventional theories for glasses (Berthier and Biroli, 2011; Kirkpatrick and Thirumalai, 2015). Second, the experimental data (Petridou et al., 2021) was interpreted using equilibrium rigidity percolation theory (Jacobs and Thorpe, 1995; Jacobs and Thorpe, 1996; Jacobs and Hendrickson, 1997) 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 nonequilibrium processes such as selfpropulsion, growth, and cell division.
We show that the two conundrums (saturation of $\eta$ at high packing fractions and the use of equilibrium statistical mechanics in a growing system to explain PT) may be rationalized by (i) assuming that the interactions between the cells are soft, and (ii) the cell sizes are highly heterogeneous (polydisperse), which is the case in zebrafish blastoderm. Using an agentbased (particle) simulation model of a twodimensional (2D) nonconfluent tissue, we explore the consequences of varying $\varphi$ (see ‘Materials and methods’ for the definition) of interacting selfpropelled polydisperse soft cells, on $\eta$. The central results of our study are (i) the calculated effective viscosity $\overline{\eta}$ (for details see, Appendix 6, ‘Dynamical changes in local packing fraction cause jammed cells to move’), for the polydisperse cell system, shows that for $\varphi \le {\varphi}_{S}\approx 0.90$ the increase in viscosity follows the VFT law. Just as in experiments, $\eta$ is essentially independent of $\varphi$ at high ($\ge {\varphi}_{S}$) packing fractions. (ii) The unusual dependence of $\eta$ at $\varphi \ge {\varphi}_{S}$ is quantitatively explained using the notion of available free area fraction ($\varphi}_{\text{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 $\varphi}_{\text{free}$ should decrease with increasing packing fractions due to cell jamming, which should slow down the overall dynamics. Indeed, we find that $\varphi}_{\text{free}$ decreases with increasing packing fraction ($\varphi$) until $\varphi}_{S$. The simulations show that when $\varphi$ exceeds $\varphi}_{S$, the free area $\varphi}_{\text{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 $\varphi$ for $\varphi \ge {\varphi}_{S}$. As a consequence, $\eta$ saturates at high $\varphi$. (iii) Cells whose sizes are comparable to the available free area move almost like a particle in a liquid. The motility of smallsized cells facilitates adjacent cells to move through multicell rearrangement even in a highly jammed environment. The facilitation mechanism, invoked in glassy systems (Biroli and Garrahan, 2013), 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 PT. (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 (Thirumalai et al., 1989) behavior at all densities. The cellbased 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 nonconfluent tissue (Petridou et al., 2021).
Results
Experimental results
We first describe the experimental observations, which serve as the basis for carrying out the agentbased simulations. Figure 1A shows the brightfield images of distinct stages during zebrafish morphogenesis. A 2D section of zebrafish blastoderm (Figure 1B) shows that there is considerable dispersion in cell sizes. The statistical properties of the cell sizes are shown in Appendix 1—figure 1D. Figure 1C shows that $\eta$ increases sharply over a narrow $\varphi$ range and saturates when $\varphi$ exceeds ${\varphi}_{S}\approx 0.90$.
To account for the results in Figure 1C, we first simulated a monodisperse system in which all the cells have identical radius ($R=8.5\text{}\mu m$). Because the system crystallizes (Appendix 1—figure 1A and 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 glasslike behavior for all $\varphi$, with the relaxation time $\tau}_{\alpha$ as well as the effective viscosity $\overline{\eta}$ (defined in Equation 1) following the VFT behavior (see Appendix 2).
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 (Petridou et al., 2021; Figure 1B, Appendix 1—figure 1D). In addition, the elastic forces characterizing cell–cell interactions are soft, which implies that the cells can overlap, with $\displaystyle {r}_{ij}({R}_{i}+{R}_{j})<0$ when they are jammed (Figure 1B, D). Thus, both polydispersity (PD) and 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 (Figure 1D , Appendix 1—figure 1E).
A simulation snapshot (Figure 1D) for $\varphi =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, $g(r)=\frac{1}{\rho}\u27e8\frac{1}{N}\sum _{i}^{N}\sum _{j\ne i}^{N}\delta (r{\overrightarrow{r}}_{i}{\overrightarrow{r}}_{j})\u27e9$, where $\rho ={\textstyle \frac{N}{{L}^{2}}}$ is the number density, $\delta$ is the Dirac delta function, $\overrightarrow{r}}_{i$ is the position of the ith cell, and the angular bracket $\u27e8\u27e9$ denotes an average over different ensembles. The $g(r)$ function (Figure 1E) has a peak around $r\sim 17\mu 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 longrange order. Thus, the polydisperse cell system exhibits liquidlike structure even at the high $\varphi$.
Effective shear viscosity($\overline{\eta}$) as a function of $\varphi$
A fit of the experimental data for $\eta$ using the VFT (Tammann and Hesse, 1926; Fulcher, 1925) relation in the range $\varphi \le 0.87$ (Figure 1C) yields ${\varphi}_{0}\approx 0.95$ and $D\approx 0.51$ (Sinha and Thirumalai, 2021). The VFT equation for cells, which is related to the Doolittle equation (White and Lipson, 2016) for fluidity ($\frac{1}{\eta}$) that is based on free space available for motion in an amorphous system (Doolittle and Doolittle, 1957; Cohen and Turnbull, 1959), is $\eta ={\eta}_{0}\mathrm{exp}\left[\frac{D}{{\varphi}_{0}/\varphi 1}\right]$, where $D$ is the apparent activation energy. In order to compare with experiments, we calculated an effective shear viscosity ($\overline{\eta}$) for the polydisperse system using a Green–Kubotype relation (Hansen and McDonald, 2013)
The stress tensor ${P}_{\mu \nu}(t)$ in the above equation is
where $\mu ,\nu \in (x,y)$ are the Cartesian components of coordinates, $\overrightarrow{r}}_{ij}={\overrightarrow{r}}_{i}{\overrightarrow{r}}_{j$, $\overrightarrow{f}}_{ij$ is the force between ith and jth cells, and $A$ is the area of the simulation box. Note that $\overline{\eta}$ should be viewed as a proxy for shear viscosity because it does not contain the kinetic term and the factor $\textstyle \frac{A}{{k}_{B}T}$ is not included in Equation (1) because temperature is not a relevant variable in the highly overdamped model for cells considered here.
Plot of $\overline{\eta}$ as a function of $\varphi$ in Figure 2A shows qualitatively the same behavior as the estimate of viscosity (using dimensional arguments) made in experiments. Two features about Figures 1C and 2A are worth noting. (i) Both simulations and experiments show that up to $\varphi \approx 0.90$, $\overline{\eta}(\varphi )$ follows the VFT relation with ${\varphi}_{0}\sim 0.94$ and $D\sim 0.5$. More importantly, $\overline{\eta}$ is independent of $\varphi$ when $\displaystyle \varphi >0.90$. (ii) The values of $\varphi}_{0$ and $D$ obtained by fitting the experimental estimate of $\eta$ to the VFT equation and simulation results are almost identical. Moreover, the onset of the plateau packing fraction in simulations and experiments occurs at the same value (${\varphi}_{S}\sim 0.90$). The overall agreement with experiments is remarkable given that the model was not created to mimic the zebrafish tissue.
To provide additional insights into the dynamics, we calculated the isotropic selfintermediate scattering function, ${F}_{s}(q,t)$,
where $\overrightarrow{q}$ is the wave vector, and ${\overrightarrow{r}}_{j}(t)$ is the position of a cell at time $t$. The degree of dynamic correlation between two cells can be inferred from the decay of ${F}_{s}(q,t)$. The angle bracket $\u27e8...\u27e9$ is an average over different time origins and different trajectories. We chose $q={\textstyle \frac{2\pi}{{r}_{\text{max}}}}$, where $r}_{\text{max}$ is the position of the first peak in $g(r)$ between all cells (see Figure 1E). The relaxation time $\tau}_{\alpha$ is calculated using $F}_{s}(q,t={\tau}_{\alpha})={\textstyle \frac{1}{e}$.
From Figure 2B and C, which show ${F}_{s}(q,t)$ as a function of $t$ for various $\varphi$, it is clear that the dynamics become sluggish as $\varphi$ increases. The relaxation profiles exhibit a twostep decay with a plateau in the intermediate time scales. The dynamics continues to slow down dramatically until $\varphi \le 0.90$. Surprisingly, the increase in the duration of the plateau in ${F}_{s}(q,t)$ ceases when $\varphi$ exceeds $\approx 0.90$ (Figure 2C), a puzzling finding that is also reflected in the dependence of $\tau}_{\alpha$ on $\varphi$ in Figure 2D. The relaxation time increases dramatically, following the VFT relation, till $\varphi \approx 0.90$, and subsequently reaches a plateau (see the inset in Figure 2D).
If the VFT relation continued to hold for all $\varphi$, as in glasses or in binary mixture of 2D cells (see Appendix 2), then the fit yields ${\varphi}_{0}\approx 0.95$ and $D\approx 0.50$. However, the simulations show that $\tau}_{\alpha$ is nearly a constant when $\varphi$ exceeds $0.90$. We should note that the behavior in Figure 2D differs from the dependence of $\tau}_{\alpha$ on $\varphi$ for 2D monodisperse polymer rings, used as a model for soft colloids. Simulations (Gnan and Zaccarelli, 2019) showed $\tau}_{\alpha$ increases till a critical $\varphi}_{S$ but it decreases substantially beyond $\varphi}_{S$ with no saturation.
Relaxation dynamics of individual cells
Plot of $\tau}_{\alpha$ as a function of the radius of cells $R}_{i$ (Figure 3A) shows nearly eight orders of magnitude change. The size dependence of $\tau}_{\alpha$ on $\varphi$ is striking. That $\tau}_{\alpha$ should increase for largesized cells (see the data beyond the vertical dashed line in Figure 3A) is not unexpected. However, even when cell sizes increase beyond ${R}_{i}=4.25\text{}\mu m$, the dispersion in $\tau}_{\alpha$ is substantial, especially when $\varphi$ exceeds $\varphi}_{S$. The relaxation times for cells with ${R}_{i}<4.25\mu m$ are relatively short even though the system as a whole is jammed. For $\varphi \ge 0.90$, $\tau}_{\alpha$ for smallsized cells have a weak dependence on $\varphi$. Although $\tau}_{\alpha$ for cells with radius <4 µm is short, it is clear that for a given $\varphi$ (e.g., $\varphi =0.93$) the variations in $\tau}_{\alpha$ are substantial. In contrast, $\displaystyle {\tau}_{\alpha}\mathrm{s}$ for larger cells ($R\ge 7\mu 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 $\u27e8{A}_{\text{free}}\u27e9$ for cells. The smallersized cells have the largest $\u27e8{A}_{\text{free}}\u27e9\approx 50\text{}\mu {m}^{2}\approx \pi {R}_{S}^{2}({R}_{S}\approx 4\text{}\mu m)$ ($R}_{S$ is the radius of the small cell).
The effect of jamming on the dramatic increase in $\tau}_{\alpha$ occurs near ${R}_{i}\approx 4.5\text{}\mu m$, which is comparable to the length scale of shortrange interactions. For $\varphi \le 0.90$, $\tau}_{\alpha$ increases as the cell size increases. However, at higher packing fractions, even cells of similar sizes show substantial variations in $\tau}_{\alpha$, which change by almost 3–4 orders of magnitude (see the data around the vertical dashed line for $\varphi \ge 0.915$ in Figure 3A).
This is a consequence of large variations in the local density (Appendix 6—figure 1). Some of the similarsized cells are trapped in the jammed environment, whereas others are in less crowded regions (see Appendix 6—figure 1). The spread in $\tau}_{\alpha$ increases dramatically for $\displaystyle \varphi >{\varphi}_{S}$ ($\approx 0.90$) and effectively overlap with each other. This is vividly illustrated in the histogram, $P(\mathrm{log}({\tau}_{\alpha}))$, shown in Figure 3B. For $\displaystyle \varphi <{\varphi}_{s}$, the peak in $P(\mathrm{log}({\tau}_{\alpha}))$ monotonically shifts to higher $\mathrm{log}({\tau}_{\alpha})$ values. In contrast, when $\varphi$ exceeds $\varphi}_{S$ there is overlap in $P(\mathrm{log}({\tau}_{\alpha}))$, which is reflected in the saturation of $\overline{\eta}$ and $\tau}_{\alpha$.
There are cells (typically with small sizes) that move faster even in a highly jammed environment (see Appendix 5—figures 1C and 2). The motions of the fastmoving cells change the local environment, which effectively facilitates the bigger cells to move in a crowded environment (see Appendix 5—figures 1D and 2, Video 1 ($\displaystyle \varphi =0.92>{\varphi}_{S}$) and Video 2 ($\varphi =0.90={\varphi}_{S}$)). In contrast, for $\displaystyle \varphi =0.85<{\varphi}_{S}$, small and largesized cells move without hindrance because of adequate availability of free area (Video 3). The videos vividly illustrate the largescale facilitated rearrangements that enable the largesized cells to move.
The dependence of $\tau}_{\alpha$ on $\varphi$ for $\varphi \le {\varphi}_{S}$ (Figure 2D) implies that the polydisperse cell systems behave as a soft glass in this regime. On theoretical grounds, it was predicted that $P(\mathrm{ln}({\tau}_{\alpha}))\sim \mathrm{exp}[c(\mathrm{ln}(\frac{{\tau}_{\alpha}}{{\tau}_{0}}){)}^{2}]$ in glassforming systems (Kirkpatrick and Thirumalai, 2015). Remarkably, we found that this prediction is valid in the polydisperse cell system (Figure 3C). However, above $\varphi}_{S$ the predicted relation is not satisfied (see Figure 3D).
Available free area explains viscosity saturation at high $\varphi$
We explain the saturation in the viscosity by calculating the available free area per cell, as $\varphi$ increases. In a hard disk system, one would expect that the free area would decrease monotonically with $\varphi$ until it is fully jammed at the close packing fraction (∼0.84; Drocco et al., 2005; Reichhardt and Reichhardt, 2014). 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 creates free area in the immediate neighborhood.
The extent of overlap ($h}_{ij$) is reflected in distribution $P({h}_{ij})$. The width in $P({h}_{ij})$ increases with $\varphi$, and the peak shifts to higher values of $h}_{ij$ (Figure 4A). The mean, $\u27e8{h}_{ij}\u27e9$, increases with $\varphi$ (Figure 4B). Thus, even if the cells are highly jammed at $\varphi \approx {\varphi}_{S}$, free area is available because of an increase in the overlap between cells (see Figure 5).
When $\varphi$ exceeds $\varphi}_{S$, the mobility of smallsized cells facilitates the larger cells to move, as is assumed in the free volume theory of polymer glasses (Cohen and Turnbull, 1959; Turnbull and Cohen, 1961; Turnbull and Cohen, 1970; Falk et al., 2020). 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 (Figure 4A and B). If $h}_{ij$ increases, as is the case when the extent of compression increases (Figure 4A), the available space for nearby cells would also increase. This effect is expected to occur with high probability at $\varphi}_{S$ and beyond, resulting in high overlap between the cells. These arguments suggest that the combined effect of PD and cell–cell overlap creates, via the selfpropulsion of cells, additional free area that drives larger cells to move even under jammed conditions.
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 Figure 5A for $\varphi =0.93$, and the histogram of $A$ is shown in Figure 5B. As $\varphi$ increases, the distribution shifts toward lower Voronoi cell size $\u27e8A\u27e9$. The mean Voronoi cell size $\u27e8A\u27e9$ as a function of $\varphi$ in Figure 5C shows $\u27e8A\u27e9$ decreases as $\varphi$ is increased. As cells interpenetrate, the Voronoi cell size will be smaller than the actual cell size ($\pi {R}_{i}^{2}$) in many instances (Figure 5A). To demonstrate this quantitatively, we calculated $A}_{\text{free,i}}={A}_{i}\pi {R}_{i}^{2$. The value of $A}_{\text{free}$ could be negative if the overlap between neighboring cells is substantial; $A}_{\text{free}$ is positive only when the Voronoi cell size is greater than the actual cell size. Positive $A}_{\text{free}$ is an estimate of the available free area. The histograms of $A}_{\text{free}$ for all the packing fractions in Figure 5D show that the distributions saturate beyond $\varphi =0.90$. All the distributions have a substantial region in which $A}_{\text{free}$ is negative. The negative value of $A}_{\text{free}$ increases with increasing $\varphi$, which implies that the amount of interpenetration between cells increases.
Because of the overlap between the cells, the available free area fraction $\varphi}_{\text{free}$ is higher than the expected free area fraction ($1.0\varphi$) for all $\varphi$. We define an effective free area fraction $\varphi}_{\text{free}$ as
where $N}_{p$ is the number of positive free area in jth snapshots, $N}_{t$ is the total number of snapshots, $A}_{\text{box}$ is the simulation box area, and $A}_{{\text{free}}_{+,i}}^{j$ is the positive free area of ith cell in jth snapshot.
The calculated $\varphi}_{\text{free}$, plotted as a function of $\varphi$ in Figure 5E, shows that $\displaystyle {\varphi}_{\text{free}}$ decreases with $\varphi$ until $\varphi =0.90$, and then it saturates near a value ${\varphi}_{\text{free}}\approx 0.22$ (see the right panel in Figure 5E). Thus, the saturation in $\overline{\eta}$ as a function of $\varphi$ is explained by the free area picture, which arises due to combined effect of the size variations and the ability of cells to overlap.
Aging does not explain viscosity saturation
Our main result, which we explain by adopting the free volume theory developed in the context of glasses (Cohen and Turnbull, 1959; Turnbull and Cohen, 1961; Turnbull and Cohen, 1970; Falk et al., 2020), is that above a critical packing fraction ${\varphi}_{S}\sim 0.90$ the viscosity saturates. Relaxation time, $\tau}_{\alpha$, measured using dynamic light scattering, in nearly monodisperse microgel poly(N isopropylacrylamide) (PNiPAM) (Philippe et al., 2018) was found to depend only weakly on the volume fraction (3D), if $\varphi}_{V$ exceeds a critical value. It was suggested that the near saturation of $\tau}_{\alpha$ at high $\varphi}_{V$ is due to aging, which is a nonequilibrium effect. If saturation in viscosity and relaxation time in the embryonic tissue at high $\varphi$ is due to aging, then $\tau}_{\alpha$ should increase sharply as the waiting time, $\tau}_{\omega$, is lengthened. We wondered if aging could explain the observed saturation of $\eta$ in the embryonic tissue above $\varphi}_{S$. If aging causes the plateau in the tissue dynamics, then $\eta$ or $\tau}_{\alpha$ should be an increasing function of the waiting time, $\tau}_{\omega$. To test the effect of $\tau}_{\omega$ on $\tau}_{\alpha$, we calculated the selfintermediate scattering function ${F}_{s}(q,t+{\tau}_{\omega})$ as a function of $t$ by varying $\tau}_{\omega$ over three orders of magnitude at $\varphi =0.92$ (Figure 6A). There is literally no change in ${F}_{s}(q,t+{\tau}_{\omega})$ over the entire range of $\tau}_{\omega$. We conclude that, $\tau}_{\alpha$, extracted from ${F}_{s}(q,t+{\tau}_{\omega})$ is independent of $\tau}_{\omega$. The variations in $\tau}_{\alpha$ (Figure 6B), with respect to $\tau}_{\omega$, are significantly smaller than the errors in the simulation. Thus, the saturation in $\eta$ or $\tau}_{\alpha$ when $\displaystyle \varphi >{\varphi}_{S}$ is not a consequence of aging.
There are two implications related to the absence of aging in the dynamics of the nonconfluent 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 $\eta$. (ii) Parenthetically, we note that the absence of significant nonequilibrium effects, even though zebrafish is a living system, further justifies the use of equilibrium rigidity percolation theory to analyze the experimental data (Petridou et al., 2021).
Discussion
Extensive computer simulations of a 2D dense tissue using a particlebased model of soft deformable cells with active selfpropulsion have successfully reproduced the dynamical behavior observed in the blastoderm tissue of zebrafish.
The dependence of viscosity ($\eta$) and relaxation time ($\tau}_{\alpha$) (before the saturation) is well fit by the VFT equation. The value of $\varphi}_{0$ obtained from simulations, ${\varphi}_{0}\sim 0.95$, is close to ${\varphi}_{0}\sim 0.94$ extracted by fitting the experimental data to the VFT equation. Thus, the dynamics for $\varphi \le {\varphi}_{S}$ resembles the behavior expected for glassforming systems. Remarkably, the dependence of $\eta$ on $\varphi$ 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 y ago. We discovered that PD 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 PT to dynamics during zebrafish development (Hannezo and Heisenberg, 2022).
One could legitimately wonder if the extent of PD used in the soft discs simulations, which seems substantial, is needed to recapitulate the observed dependence of $\eta$ on $\varphi$. Furthermore, such large values of PD may not represent biological tissues. Although the choice of PD was made in part by the 2D projection of area reported in experiments (Petridou et al., 2021), 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 ($\tau}_{\alpha$) as a function of $\varphi$. The results show that $\tau}_{\alpha$ does indeed saturate at highvolume 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 experimental data. We do think that realistic values of adhesion strength would not significantly alter the forces between cells (MalmiKakkada et al., 2018). Thus, we expect a similar mechanism. Furthermore, the physics of the dynamics in glassforming materials does not change in systems with and without attractive forces (Kirkpatrick and Thirumalai, 2015). Universal behavior, such as VFT relation, is valid for a broad class of unrelated materials (see Figure 1 in Angell, 1991). Needless to say, nonuniversal quantities such as glass transition temperature $T}_{g$ 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 $\varphi}_{S$ without qualitatively altering the dependence of $\eta$ on $\varphi$. 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 Equation (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 $\gamma}_{i$ in Equation (6) through fluctuation dissipation theorem (FDT). Thus, only $\gamma}_{i$ can be varied. In contrast, in Equation (6) the two parameters ($\gamma}_{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 (Ikeda et al., 2012), 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 celltocell contact topology. These considerations, which are firmly established here, explain why equilibrium PTs are related to a steep increase in viscosity (Kirkpatrick and Thirumalai, 2015) 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.
Materials and methods
Twodimensional cell model
Request a detailed protocolFollowing our earlier studies (MalmiKakkada et al., 2018; Sinha et al., 2020), we simulated a 2D version of a particlebased cell model. We did not explicitly include cell division in the simulations. This is physically reasonable because in the experiments (Petridou et al., 2021) 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 (MatozFernandez et al., 2017; Drasdo and Höhme, 2005; Schaller and MeyerHermann, 2005; MalmiKakkada et al., 2018) interacting via shortranged forces. The elastic (repulsive) force between two cells with radii $R}_{i$ and $R}_{j$ is Hertzian, which is given by
where ${h}_{ij}=\text{max}[0,{R}_{i}+{R}_{j}{\overrightarrow{r}}_{i}{\overrightarrow{r}}_{j}]$. The repulsive force acts along the unit vector $\overrightarrow{n}}_{ij$, which points from the center of the jth cell to the center of the ith cell. The total force on the ith cell is
where $NN(i)$ is the number of nearneighbor cells that are in contact with the ith cell. The jth cell is the nearest neighbor of the ith cell, if $\displaystyle {h}_{ij}>0$. The nearneighbor 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 ($\nu$) for all the cells are identical. PD 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 ($A}_{i}=\pi {R}_{i}^{2$) is assumed to have a distribution that mimics the broad area distribution discovered in experiments.
Selfpropulsion and equations of motion
Request a detailed protocolIn addition to the repulsive Hertz force, we include an active force arising from selfpropulsion 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 $\gamma}_{i$ is the friction coefficient of ith cell, and ${\mathcal{W}}_{i}(t)$ is a noise term. The friction coefficient $\gamma}_{i$ is taken to be $\gamma}_{0}{R}_{i$ (Sinha et al., 2022). By scaling $t$ by the characteristic time scale, $\tau =\frac{{\u27e8R\u27e9}^{2}}{{\mu}^{2}}$ in Equation (6), one can show that the results should be insensitive to the exact value of µ. The noise term ${\mathcal{W}}_{i}(t)$ is chosen such that $\u27e8{\mathcal{W}}_{i}(t)\u27e9=0$ and $\u27e8{\mathcal{W}}_{i}^{\alpha}(t){\mathcal{W}}_{j}^{\beta}({t}^{\mathrm{\prime}})\u27e9=\delta (t{t}^{\mathrm{\prime}}){\delta}_{i,j}{\delta}^{\alpha ,\beta}$. In our model, there is no dynamics with only systematic forces because the temperature is zero. The observed dynamics arises solely due to the selfpropulsion (Equation 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 2D system it is the area fraction) is $\varphi ={\textstyle \frac{\sum _{i=1}^{N}\pi {R}_{i}^{2}}{{L}^{2}}}$. We performed extensive simulations by varying $\varphi$ in the range $0.700\le \varphi \le 0.950$. The results reported in main text are obtained with $N=500$. Finite size effects are discussed in Appendix 7.
To mimic the variations in the area of cells in a tissue (Petridou et al., 2021), we use a broad distribution of cell radii (see Appendix 1 for details). The parameters for the model are given in Table 1. In this study, we do not consider the growth and division of cells. Thus, our simulations describe steadystate dynamics of the tissue. For each $\varphi$, we performed simulations for at least (5–10)$\tau}_{\alpha$ before storing the data. For each $\varphi$, we performed 24 independent simulations. The calculation of viscosity was performed by averaging over 40 independent simulations at each $\varphi$.
Calculation of viscosity
Request a detailed protocolWe calculated the effective viscosity ($\overline{\eta}$) for various values of $\varphi$ by integrating the offdiagonal part of the stress–stress correlation function $\u27e8{P}_{\mu \nu}(t){P}_{\mu \nu}(0)\u27e9$ using the Green–Kubo relation (Hansen and McDonald, 2013) (without the prefactor $\textstyle \frac{A}{{k}_{B}T}$)
where µ and $\nu$ denote Cartesian components ($x$ and $y$) of the stress tensor ${P}_{\mu \nu}(t)$ (see main text for the definition of ${P}_{\mu \nu}(t)$). The definition of $\overline{\eta}$, which relates the decay of stresses as a function of times in the nonconfluent tissue, is akin to the methods used to calculate viscosity in simple fluids (Equation 7). The time dependence of $\u27e8{P}_{\mu \nu}(t){P}_{\mu \nu}(0)\u27e9$, normalized by $\u27e8{P}_{\mu \nu}(0{)}^{2}\u27e9$, for different values of $\varphi$ (Figure 7A and B) shows that the stress relaxation is clearly nonexponential, decaying to zero in two steps. After an initial rapid decay followed by a plateau at intermediate times (clearly visible for $\varphi \ge 0.91$), the normalized $\u27e8{P}_{\mu \nu}(t){P}_{\mu \nu}(0)\u27e9$ decays to zero as a stretched exponential. The black dashed lines in Figure 7C show that a stretched exponential function, ${C}_{s}\mathrm{exp}[{\left({\textstyle \frac{t}{{\tau}_{\eta}}}\right)}^{\beta}]$, where $\tau}_{\eta$ is the characteristic time in which stress relax and $\beta$ is the stretching exponent, provides an excellent fit to the long time decay of $\displaystyle \u27e8{P}_{\mu \nu}(t){P}_{\mu \nu}(0)\u27e9$ (from the plateau region to zero) as a function of $t$. Therefore, we utilized the fit function, ${C}_{s}\mathrm{exp}[{\left({\textstyle \frac{t}{{\tau}_{\eta}}}\right)}^{\beta}]$, to replace the noisy long time part of $\u27e8{P}_{\mu \nu}(t){P}_{\mu \nu}(0)\u27e9$ by a smooth fit data before evaluating the integral in Equation (7). The details of the procedure to compute $\overline{\eta}$ are described below.
We divided $\u27e8{P}_{\mu \nu}(t){P}_{\mu \nu}(0)\u27e9$ in two parts. (i) The short time part ($\u27e8{P}_{\mu \nu}(t){P}_{\mu \nu}(0)\u27e9}_{\text{short}$) – the smooth initial rapid decay until the plateau is reached (e.g., see the blue circles in Figure 7D for $\varphi =0.93$). For the $n$ data points at short times, $({t}_{1},{\u27e8{P}_{\mu \nu}({t}_{1}){P}_{\mu \nu}(0)\u27e9}_{\text{\% short}})$,…, $({t}_{n},{\u27e8{P}_{\mu \nu}({t}_{n}){P}_{\mu \nu}(0)\u27e9}_{\text{\% short}})$, we constructed a spline $S(t)$ using a set of cubic polynomials:
The polynomials satisfy the following properties. (a) $S}_{i}({t}_{i})={\u27e8{P}_{\mu \nu}({t}_{i}){P}_{\mu \nu}(0)\u27e9}_{\text{\% short}$ and $S}_{i}({t}_{i+1})={\u27e8{P}_{\mu \nu}({t}_{i+1}){P}_{\mu \nu}(0)\u27e9}_{\text{short}$ for $\displaystyle i=1,...,n1$ which guarantees that the spline function $S(t)$ interpolates between the data points. (b) ${S}_{i1}^{\mathrm{\prime}}(t)={S}_{i}^{\mathrm{\prime}}(t)$ for $\displaystyle i=2,...,n1$ so that ${S}^{\mathrm{\prime}}(t)$ is continuous in the interval $[{t}_{1},{t}_{n}]$. (c) ${S}_{i1}^{\mathrm{\prime}\mathrm{\prime}}(t)={S}_{i}^{\mathrm{\prime}\mathrm{\prime}}(t)$ for $\displaystyle i=2,...,n1$ so that ${S}^{\mathrm{\prime}\mathrm{\prime}}(t)$ is continuous in the interval $[{t}_{1},{t}_{n}]$. By solving for the unknown parameters, $b}_{i},{c}_{i$, and $d}_{i$, using the abovementioned properties, we constructed the function S(t). We used $S(t)$ to fit $\u27e8{P}_{\mu \nu}(t){P}_{\mu \nu}(0)\u27e9}_{\text{short}$ to get an evenly spaced ($\delta t=10s$) smooth data (solid blue line in Figure 7D). The fitting was done using the software ‘Xmgrace’.
(ii) The long time part ($\u27e8{P}_{\mu \nu}(t){P}_{\mu \nu}(0)\u27e9}_{\text{long}$) – from the plateau until it decays to zero – is shown by the red circles in Figure 7D. The long time part was fit using the analytical function ${C}_{s}\mathrm{exp}[{\left({\textstyle \frac{t}{{\tau}_{\eta}}}\right)}^{\beta}]$ (black dashed line in Figure 7D). We refer to the fit data ($\delta t=10s$) as $\u27e8{P}_{\mu \nu}(t){P}_{\mu \nu}(0)\u27e9}_{\text{long}}^{\text{fit}$.
We then combined $\u27e8{P}_{\mu \nu}(t){P}_{\mu \nu}(0)\u27e9}_{\text{short}$ and $\u27e8{P}_{\mu \nu}(t){P}_{\mu \nu}(0)\u27e9}_{\text{long}}^{\text{\% fitted}$ to obtain $\u27e8{P}_{\mu \nu}(t){P}_{\mu \nu}(0)\u27e9}_{\text{combined}$ (see inset of Figure 7D). Finally, we calculated $\overline{\eta}$ using the equation,
where ${t}_{1}\delta t$ is the end point of $\u27e8{P}_{\mu \nu}(t){P}_{\mu \nu}(0)\u27e9}_{\text{short}$ and $T\delta t$ is the end point of $\u27e8{P}_{\mu \nu}(i\delta t){P}_{\mu \nu}(0)\u27e9}_{\text{combined}$.
Appendix 1
Cell polydispersity is needed to account for viscosity saturation
To explain the observed saturation in viscosity (Figure 1C) when cell area fraction, $\varphi$, exceeds $\varphi}_{S$ ($\approx 0.90$), we first simulated a monodisperse cell system using the model described in the ‘Materials and methods’ section. The monodisperse cell system ($R=8.5\mu m$) crystallizes (see Appendix 1—figure 1A 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 cannot account for the experimental data even though crystallization is avoided (see the next section). These findings forced us to take PD into account.
In order to develop an agentbased model that accounts for PD effects, we first extracted the distribution, $P({A}_{i})$, of the area ($A}_{i$) in the zebrafish cells from the experimental data (Petridou et al., 2021). Appendix 1—figure 1D shows that $P({A}_{i})$ is broad, implying that cell sizes are highly heterogeneous. Based on this finding, we sought a model of the nonconfluent tissue that approximately mimics PD found in experiments. In other words, $\textstyle \frac{\mathrm{\Delta}A}{\u27e8A\u27e9}$ ($\mathrm{\Delta}A$ is the dispersion in $P({A}_{i})$ and $\u27e8A\u27e9$ is the mean value) should be similar to the data in Appendix 1—figure 1D. The radii ($R}_{i$) of the cells in the simulations are sampled from a Gaussian distribution $\sim \mathrm{exp}(({R}_{i}\u27e8R\u27e9{)}^{2}/2\mathrm{\Delta}{R}^{2})$, with $\u27e8R\u27e9=8.5\mu m$ and $\mathrm{\Delta}R=4.5\mu m$. The resulting $P({A}_{i})$ for one of the realizations (see Appendix 1—figure 1C) is shown in Appendix 1—figure 1E. The value of $\textstyle \frac{\mathrm{\Delta}A}{\u27e8A\u27e9}$ is $\sim 0.5$, which compares favorably with the experimental estimate (∼0.4). The unusual dependence of the viscosity ($\eta$) as a function of packing fraction ($\varphi$) cannot be reproduced in the absence of PD.
Appendix 2
Relaxation time in the binary mixture of cells does not saturate at high $\varphi$
We showed that the saturation in the relaxation time above a critical area fraction, $\varphi}_{S$, is related to two factors. One is that there ought to be dispersion in the cell sizes (Appendix 1—figure 1). 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 nonadditive. The Hertz potential captures the squishy nature of the cells.
In order to reveal the importance of PD, we first investigated if a binary mixture of cells (see Appendix 2—figure 1A) would reproduce the observed dependence of $\tau}_{\alpha$ on $\varphi$. We created a 50:50 mixture with a cell size ratio $\sim 1.4$ (${R}_{B}=8.5\mu m$ and ${R}_{S}=6.1\mu m$). All other parameters are the same as in the polydisperse system (see Table 1). The radial distribution function, $g(r)$, calculated by considering both the cell types, exhibits intermediate but not longrange translational order (Appendix 2—figure 1B). A peak in $g(r)$ appears at ${r}_{max}=14.0\mu m$.
In order to determine the dependence of the relaxation time, $\tau}_{\alpha$, on $\varphi$, we first calculated ${F}_{s}(q,t)$ at $q=2\pi /{r}_{max}$ (Appendix 2—figure 1C). The decay of ${F}_{s}(q,t)$ is similar to what one finds in typical glassforming systems.
As $\varphi$ increases, the decay of ${F}_{s}(q,t)$ slows down dramatically. When $\varphi$ exceeds $0.93$, there is a visible plateau at intermediate times followed by a slow decay. By fitting the long time decay of ${F}_{s}(q,t)$ to $\mathrm{exp}((t/{\tau}_{\alpha}{)}^{\beta})$ ($\beta$ is the stretching exponent), we find that ${\tau}_{\alpha}(\varphi )$ as a function of $\varphi$ is well characterized by the VFT relation (Appendix 2—figure 1D). There is no evidence of saturation in ${\tau}_{\alpha}(\varphi )$ at high $\varphi$. The $\varphi$ dependence of the effective shear viscosity $\overline{\eta}$ (Appendix 2—figure 1E), calculated using Equation (8), also shows no sign of saturation. The VFT behavior, with ${\varphi}_{0}\approx 1$ and $D\approx 1.2$, shows that the twocomponent cell system behaves as a ‘fragile’ glass (Angell, 1991).
A comment regarding the binary system of cells is in order. The variables that characterize this system are $\lambda ={R}_{B}/{R}_{S}$ ($R}_{B$ ($R}_{S$) is the radius of the big (small) cells), ${\mathrm{\Phi}}_{B}={N}_{B}/({N}_{A}+{N}_{B})$ ($N}_{B$ ($N}_{A$) is the number of big (small) cells) with $\mathrm{\Phi}}_{A}=1{\mathrm{\Phi}}_{B$, and the packing fraction, $\varphi =\pi ({N}_{B}{R}_{B}^{2}+{N}_{A}{R}_{A}^{2})/{A}_{S}$, where $A}_{S$ is the area of the sample. The results in Appendix 2—figure 1 were obtained using $\lambda =1.4$, ${\mathrm{\Phi}}_{B}=0.5$. With this choice, the value of PD is $(\lambda 1)/(\lambda +1)\approx 0.17$, which is smaller than the experimental value. It is possible that by thoroughly exploring the parameter space, $\lambda$ and $\varphi$, one could find regions in which the $\tau}_{\alpha$ in the twocomponent cell system would saturate beyond $\varphi}_{S$. However, in light of the results in Appendix 1—figure 1D, we choose to simulate systems that also have high degree of PD (Appendix 1—figure 1E).
Appendix 3
Absence of saturation in the free area in binary mixture of cells
In the main text, we established that the effective viscosity, $\overline{\eta}$, which should be a proxy for the true $\eta$ and the relaxation time ($\tau}_{\alpha$) in the polydisperse cell system, saturates beyond $\varphi}_{S$. The dynamics saturates beyond $\varphi}_{S$ because the available area per cells (quantified as $\varphi}_{\text{free}$) is roughly a constant in this region (see Figure 5E). In the binary system, however, we do not find any saturation in the $\tau}_{\alpha$. Therefore, if $\displaystyle {\varphi}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}$ controls the dynamics, one would expect that $\varphi}_{\text{free}$ should decrease monotonically with $\varphi$ for the binary cells system. We calculated the average Vornoi cell size $\u27e8A\u27e9$ (Appendix 3—figure 1A) and $\varphi}_{\text{free}$ (Appendix 3—figure 1B) for the binary system. Both $\u27e8A\u27e9$ and $\varphi}_{\text{free}$ decrease with $\varphi$, which is consistent with the free volume picture proposed in the context of glasses.
Appendix 4
Absence of broken ergodicity
We have shown in the earlier section that the saturation of $\tau}_{\alpha$ above $\varphi}_{S$ is not a consequence of aging. In the context of glasses and supercooled liquids (Thirumalai et al., 1989), 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 nonconfluent tissues, we define a measure $\mathrm{\Omega}(t)$ given by
where $k$ and $l$ represent systems with different initial conditions, and ${\omega}^{k}(t)=\frac{1}{t}{\int}_{0}^{t}{P}_{\mu \nu}^{k}(s)ds$ with ${P}_{\mu \nu}(s)$ being the value of stress (see Equation 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 ${\omega}^{k}(t)$ and ${\omega}^{l}(t)$ would be independent of $k$ and $l$. Therefore, $\mathrm{\Omega}(t)$ should vanish at very long time for ergodic systems. On the other hand, if ergodicity is broken, then $\mathrm{\Omega}(t)$ would be a constant whose value would depend on $k$ and $l$.
The long time values of $\mathrm{\Omega}(t)$, normalized by $\mathrm{\Omega}(0)$ (at $t=0$), are 0.01, 0.016, and 0.026 for $\varphi =0.85,0.90$, and 0.92, respectively (Appendix 4—figure 1A–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 (Thirumalai et al., 1989) that in long time the ergodic measure ($\mathrm{\Omega}(t)$ in our case) should decay as $\approx 1/t$. Appendix 4—figure 1D shows that this is indeed the case – at long time $\mathrm{\Omega}(t)/\mathrm{\Omega}(0)$ decays approximately as $1/t$.
Appendix 5
Dynamics of small and large cells are dramatically different
The structural and the dynamical behavior of the small (${R}_{S}\le 4.5\mu m$) and large (${R}_{B}\ge 12.0\mu m$) cells is dramatically different in the nonconfluent tissue. The pair correlation functions between small cells (${g}_{SS}(r)$) and between large cells (${g}_{BB}(r)$) (Appendix 5—figure 1A and B) for $\varphi =0.905$ and $\varphi =0.92$ show that both small and large cells exhibit liquidlike disordered structures. However, it is important to note that ${g}_{SS}(r)$ has only one prominent peak and a modest second peak. In contrast, ${g}_{BB}(r)$ has three prominent peaks. Thus, the smallersized cells exhibit liquidlike behavior, whereas the large cells are jammed. This structural feature is reflected in the decay of ${F}_{s}(q,t)$ with $q={\textstyle \frac{2\pi}{{r}_{max}}}$, where $r}_{max$ is the position of the first peak in the $g(r)$ (Appendix 5—figure 1C and D).
There is a clear difference in the decay of ${F}_{s}(q,t)$, spanning nearly eight orders of magnitude, between small and large cells (compare Appendix 5—figure 1C and D). At the highly jammed packing fraction ($\varphi =0.92$), smallsized cells have the characteristics of fluidlike behavior.
A typical snapshot from one of the simulations ($\varphi =0.91$) and the trajectories for a few small and bigsized cells are displayed in Appendix 5—figure 2A–D. The figures reflect the decay in ${F}_{s}(q,t)$. Not only are the mobilities heterogeneous, it is also clear that the displacements of the small cells are greater than the large cells.
Appendix 6
Dynamical changes in local packing fraction cause jammed cells to move
The mobilities of all the cells, even under highly jammed conditions ($\varphi \ge {\varphi}_{S}$), are only possible if the local area fraction changes dynamically. This is a collective effect, which is difficult to quantify because it would involve multicell 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 glassforming materials (Hedges et al., 2009; Biroli and Garrahan, 2013). However, in glasslike systems, facilitation is due to thermal excitation, but in the active system, it is selfpropulsion that causes the cells to move.
The creation of free space may be visualized by tracking the positions of the nearestneighbor cells. In Appendix 6—figure 1, we display the local free area of a blackcolored cell at different times. The top panels show the configurations where the black cell is completely jammed by other cells. The cell, colored in black, can move if the neighboring cells rearrange (caging effect in glassforming 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.
Appendix 7
Finite system size effects
In the main text, we report results for $N=500$. To asses if the unusual dynamics is not an effect of finite system size, we performed additional simulations with $N=200$ and $N=750$. As shown in Appendix 7—figure 1A and B, ${F}_{s}(q,t)$ saturates at $\varphi \ge {\varphi}_{S}$, which is reflected in the logarithm of $\tau}_{\alpha$ as a function of $\varphi$ (Appendix 7—figure 1C and D). The saturation value ${\varphi}_{S}\sim 0.90$ is independent of the system size. The value of $\varphi}_{0$ ($\approx 0.95$) is also nearly independent of system size. Therefore, the observed dynamics, reflected in the plateau in the viscosity at high $\varphi$, is likely not an effect of finite system size.
Appendix 8
Dependence of viscosity on average coordination number and average connectivity
In the zebrafish blastoderm experiment (Petridou et al., 2021), the change in viscosity ($\eta$) was shown as a function of mean connectivity $\u27e8C\u27e9$. To test whether our 2D tissue simulations also exhibit similar behavior, we calculated viscosity as a function of mean coordination number $\u27e8{N}_{c}\u27e9$ (defined below), which is equivalent to $\u27e8C\u27e9$. We define the coordination number, $N}_{c$, as the number of cells that are in contact with a given cell. Two cells with indices $i$ and $j$ are in contact if $\displaystyle {h}_{ij}={R}_{i}+{R}_{j}{r}_{ij}>0$. We calculate $N}_{c$ for all the cells for each $\varphi$ and calculate the histogram $P({N}_{c})$. The distributions of $\displaystyle P({N}_{c})$ are well fit by a Gaussian distribution function $A\mathrm{exp}[(\frac{x\mu}{\sigma}{)}^{2}]$ (Appendix 8—figure 1A–C). The calculated mean, $\u27e8{N}_{c}\u27e9$, from the fit is linearly related to the cell area fraction $\varphi$ (Appendix 8—figure 1D). We also calculate the average connectivity $\u27e8C\u27e9$ 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 $C=\frac{2m}{n}$ (Petridou et al., 2021). We calculate $C$ for all the snapshots for each $\varphi$ and estimated its mean value $\u27e8C\u27e9$. We find that $\u27e8C\u27e9$ and $\u27e8{N}_{c}\u27e9$ are of similar values (Appendix 8—figure 2A), and the dependence of viscosity $\overline{\eta}$ on $\u27e8{N}_{c}\u27e9$ and $\u27e8C\u27e9$ is similar to experimental results (Appendix 8—figure 2B and C; Petridou et al., 2021).
Appendix 9
Connectivity map
To pictorially observe the percolation transition in the simulations, we plot the connectivity map for various values $\varphi$ in Appendix 9—figure 1. Note that for smaller $\varphi \le 0.89$ the map shows that cells are loosely connected, suggesting a fluidlike behavior. For $\varphi \ge 0.89$, the connectivity between the cells spans the entire system, which was noted and analyzed elsewhere thoroughly (Petridou et al., 2021). 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 nonpercolated state to a percolated state occurs over a very narrow range of $\varphi$, which corresponds to the onset of rigidity percolation transition (Petridou et al., 2021). It is gratifying that the simple simulations reproduce the experimental observations.
Data availability
The experimental data was extracted from Petridou et al., 2021, and the data is available in GitHub link: https://github.com/rajsekhardas88/eLife (copy archived at Rajsekhar, 2023).
References

Relaxation in liquids, polymers and plastic crystals — strong/fragile patterns and problemsJournal of NonCrystalline Solids 131–133:13–31.https://doi.org/10.1016/00223093(91)902669

Critical phenomena in embryonic organizationCurrent Opinion in Systems Biology 31:100433.https://doi.org/10.1016/j.coisb.2022.100433

SelfPropelled rods: Insights and perspectives for active matterAnnual Review of Condensed Matter Physics 11:441–466.https://doi.org/10.1146/annurevconmatphys031119050611

Theoretical perspective on the glass transition and amorphous materialsReviews of Modern Physics 83:587–645.https://doi.org/10.1103/RevModPhys.83.587

Perspective: The glass transitionThe Journal of Chemical Physics 138:12A301.https://doi.org/10.1063/1.4795539

Molecular transport in liquids and glassesThe Journal of Chemical Physics 31:1164–1169.https://doi.org/10.1063/1.1730566

Studies in Newtonian Flow. V. Further Verification of the FreeSpace Viscosity EquationJournal of Applied Physics 28:901–905.https://doi.org/10.1063/1.1722884

Multiscaling at Point J: jamming is a critical phenomenonPhysical Review Letters 95:088001.https://doi.org/10.1103/PhysRevLett.95.088001

Nonempirical free volume viscosity model for alkane lubricants under severe pressuresPhysical Review Letters 124:105501.https://doi.org/10.1103/PhysRevLett.124.105501

Analysis of recent measurements of the viscosity of glassesJournal of the American Ceramic Society 8:339–355.https://doi.org/10.1111/j.11512916.1925.tb16731.x

Aspiration of biological viscoelastic dropsPhysical Review Letters 104:218101.https://doi.org/10.1103/PhysRevLett.104.218101

Rigidity transitions in development and diseaseTrends in Cell Biology 32:433–444.https://doi.org/10.1016/j.tcb.2021.12.006

Unified study of glass and jamming rheology in soft particle systemsPhysical Review Letters 109:018301.https://doi.org/10.1103/PhysRevLett.109.018301

Generic rigidity percolation: The pebble gamePhysical Review Letters 75:4051–4054.https://doi.org/10.1103/PhysRevLett.75.4051

Generic rigidity percolation in two dimensionsPhysical Review E 53:3682–3693.https://doi.org/10.1103/PhysRevE.53.3682

An Algorithm for TwoDimensional Rigidity Percolation: The Pebble GameJournal of Computational Physics 137:346–365.https://doi.org/10.1006/jcph.1997.5809

Stages of embryonic development of the zebrafishDevelopmental Dynamics 203:253–310.https://doi.org/10.1002/aja.1002030302

Colloquium : Random first order transition theory concepts in biology and physicsReviews of Modern Physics 87:183–209.https://doi.org/10.1103/RevModPhys.87.183

Hydrodynamics of soft active matterReviews of Modern Physics 85:1143–1189.https://doi.org/10.1103/RevModPhys.85.1143

Tissue rheology in embryonic organizationThe EMBO Journal 38:e102497.https://doi.org/10.15252/embj.2019102497

Glass transition of soft colloidsPhysical Review. E 97:040601.https://doi.org/10.1103/PhysRevE.97.040601

Multicellular tumor spheroid in an offlattice VoronoiDelaunay cell modelPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 71:051910.https://doi.org/10.1103/PhysRevE.71.051910

Computational models for active matterNature Reviews Physics 2:181–199.https://doi.org/10.1038/s4225402001521

Rigidity percolation in embryo morphogenesis: Physics meets biology (again)Journal Club for Condensed Matter Physics 7:1914–1928.https://doi.org/10.36471/JCCM_June_2021_03

Mechanical feedback controls the emergence of dynamical memory in growing tissue monolayersThe Journal of Chemical Physics 156:245101.https://doi.org/10.1063/5.0087815

Die Abhängigkeit der Viscosität von der Temperatur bie unterkühlten FlüssigkeitenZeitschrift Für Anorganische Und Allgemeine Chemie 156:245–257.https://doi.org/10.1002/zaac.19261560121

Ergodic behavior in supercooled liquids and in glassesPhysical Review A 39:3563–3574.https://doi.org/10.1103/PhysRevA.39.3563

FreeVolume model of the amorphous phase: glass transitionThe Journal of Chemical Physics 34:120–125.https://doi.org/10.1063/1.1731549

On the FreeVolume Model of the LiquidGlass TransitionThe Journal of Chemical Physics 52:3038–3041.https://doi.org/10.1063/1.1673434

Polymer free volume and its connection to the glass transitionMacromolecules 49:3987–4007.https://doi.org/10.1021/acs.macromol.6b00215
Peer review
Joint Public Review:
This paper explores how minimal active matter simulations can model tissue rheology, with applications to the in vivo situation of zebrafish morphogenesis. The authors explore the idea of active noise, particle softness and size heterogeneity cooperating to give rise to surprising features of experimental tissue rheologies (in particular an increase and then a plateau in viscosity with fluid fraction). In general, the paper is interesting from a theoretical standpoint, by providing a bridge between concepts from jamming of particulate systems and experiments in developmental biology. The idea of exploring a free space picture in this context is also interesting. It will be interesting in the future to see whether and how the findings change when considering 3D tissues with less size heterogeneity or how viscosity is impacted by the time scale of measurements.
https://doi.org/10.7554/eLife.87966.4.sa1Author response
The following is the authors’ response to the previous reviews.
We are pleased that Reviewers 1 and 3 have recommended that the revised paper be published.
Reviewer #2
For point A: Their preliminary simulation in 3D looks also nice, although it’s referenced in the discussion but not actually included in manuscript  I would advise adding it even under the mention of preliminary.
We appreciate the reviewer for liking our 3D results and suggesting to include them in the manuscript. However, these are preliminary results of our ongoing work. We are yet to establish the corresponding viscosity results quantitatively in the 3D simulations. Because the relationship between viscosity and relaxation time is not (always) linear in glass forming systems, we hesitate to report our results for publication. We hope to report the new results as part of a separate work.
For point B/C: I see some of the points of the authors  although not all of it made it in the main text. I still have some points that puzzle me. For instance, the authors mention that a single value of viscosity (from GreenKubo) is ”valid for all time scales and amplitude”. This sounds very surprising to me for a complex fluid even at equilibrium: doesn’t it for instance assume linear response (hence small amplitudes)? Fast vs slow probing of a complex medium should also matter (see refs previously mentioned). Related to this, it’s not clear how can selfpropulsion not matter if one would shear the system at a finite time scale, given past work on motilitydriven unjamming and the mechanism of the authors from facilitation ( wouldn’t shearing at time scales larger vs smaller than the typical time for given cells to spontaneously rearrange from selfpropulsion change drastically the effective complex modulus of the system?)
There might be a slight misunderstanding between the reviewer and us when
we say ‘single value of viscosity is valid for all timescales and amplitude’. Let us explain this point more carefully. In our problem, we are studying the dynamics of a many body system which is undergoing Brownian dynamics where the fluctuationdissipation theorem need not be valid (as the friction and the selfpropulsion noise strength are not related via FluctuationDissipation Theorem). Now, for us to use the concepts of linearresponse (which in the present study are the GreenKubo relations for the transport coefficients in terms of timecorrelations functions), we need to show that the within the simulation time, the system has reached state that could be described using an “equilibrium” probability measure. This is the precise reason we calculated the ergodicity measure, which is a way to show that all the phasespace have been sampled uniformly under the given Brownian dynamics. This suggests (does not prove) that the system has attained a stationary probability measure (i.e, near equilibrium) for the value of selfpropulsion used. Now for this value of selfpropulsion, the GreenKubo relations hold for ‘any timescale of the simulations’ so that we can perform a time average over the trajectories of the particles (which is an alias of the stationary probability measure under the values of selfpropulsion used). If we change the amplitude of the selfpropulsion, we need to again compute the ergodicity measure and show the stationarity of the probability measure. If the system is ergodic with respect to the new selfpropulsion, we can again use GreenKubo for the simulations. Note that we will definitely get a different value of viscosity under the new selfpropulsion as the shearstresses generated will be different but the GreenKubo holds. If the system is not ergodic, for the selfpropulsion with the new amplitude, we cannot use GreenKubo relations. Also a priori, one cannot say what is a large/small amplitude of selfpropulsion because it has to be compared with the intrinsic energy scale, which is encoded in the energy function, which is difficult to say without explicit calculations.
This is what we meant when we said, ‘single value of viscosity is valid for all timescales and amplitude’. It is valid for timescales of the simulations for a given amplitude of selfpropulsion only if the system is ergodic. Note that if the system is not ergodic, then the results of Ref. [14] (in the main text) could be questioned on theoretical grounds, because they were analyzed using 3 the equilibrium rigidity percolation theory. Nevertheless, the authors of Ref. [14] showed that equilibrium phase transition theory works in tissues. For these reasons, we have been, just like the Reviewer, puzzled that equilibrium ideas appear to be valid in the cell system. Additional theoretical work has to be done to clarify these links in tissues. Although this is not the last word, we hope this clarifies our view point.
For point D: I agree with the simplicity argument, although the added sentence from the discussion “Furthermore, the physics of the dynamics in glass forming materials does not change in systems with and without attractive forces” seems a bit strong given works like Lois et al., PRL, 2008 or Koeze et al, PRL, 2018 finding fundamentally different physics of jamming with or without adhesion.In the two cited papers the authors only consider equilibrium transitions in systems with attraction using computer simulations. Apparently, jamming properties depend on the strength of attraction. There are no attempts to characterize the dynamics, the focus of our work.
What we meant is that any universal relations, such as the VogelFulcherTammann relation, would still be valid. Of course, nonuniversal quantities such as glass transition temperature Tg or fragility will change. In our case, changing the adhesion strength would change ϕS, and the parameters in the VFT. However, our contention is that the overall finding that increase in viscosity followed by saturation is unlikely to change. We have added some clarifying statements in the manuscript to make this clear.
https://doi.org/10.7554/eLife.87966.4.sa2Article and author information
Author details
Funding
National Science Foundation (2345242)
 D Thirumalai
Welch Foundation (F0019)
 D Thirumalai
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
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 2310639 and CHE 2320256) and the Welch Foundation (F0019).
Senior Editor
 Aleksandra M Walczak, École Normale Supérieure  PSL, France
Reviewing Editor
 Karsten Kruse, University of Geneva, Switzerland
Version history
 Preprint posted: November 26, 2022 (view preprint)
 Sent for peer review: April 17, 2023
 Preprint posted: July 3, 2023 (view preprint)
 Preprint posted: September 28, 2023 (view preprint)
 Preprint posted: December 18, 2023 (view preprint)
 Version of Record published: January 19, 2024 (version 1)
Cite all versions
You can cite all versions using the DOI https://doi.org/10.7554/eLife.87966. This DOI represents all versions, and will always resolve to the latest one.
Copyright
© 2023, Das et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 185
 Page views

 31
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Physics of Living Systems
Many animals moving through fluids exhibit highly coordinated group movement that is thought to reduce the cost of locomotion. However, direct energetic measurements demonstrating the energysaving benefits of fluidmediated collective movements remain elusive. By characterizing both aerobic and anaerobic metabolic energy contributions in schools of giant danio (Devario aequipinnatus), we discovered that fish schools have a concave upward shaped metabolism–speed curve, with a minimum metabolic cost at ~1 body length s^{1}. We demonstrate that fish schools reduce total energy expenditure (TEE) per tail beat by up to 56% compared to solitary fish. When reaching their maximum sustained swimming speed, fish swimming in schools had a 44% higher maximum aerobic performance and used 65% less nonaerobic energy compared to solitary individuals, which lowered the TEE and total cost of transport by up to 53%, near the lowest recorded for any aquatic organism. Fish in schools also recovered from exercise 43% faster than solitary fish. The nonaerobic energetic savings that occur when fish in schools actively swim at high speed can considerably improve both peak and repeated performance which is likely to be beneficial for evading predators. These energetic savings may underlie the prevalence of coordinated group locomotion in fishes.

 Computational and Systems Biology
 Physics of Living Systems
The adaptive dynamics of evolving microbial populations takes place on a complex fitness landscape generated by epistatic interactions. The population generically consists of multiple competing strains, a phenomenon known as clonal interference. Microscopic epistasis and clonal interference are central aspects of evolution in microbes, but their combined effects on the functional form of the population’s mean fitness are poorly understood. Here, we develop a computational method that resolves the full microscopic complexity of a simulated evolving population subject to a standard serial dilution protocol. Through extensive numerical experimentation, we find that stronger microscopic epistasis gives rise to fitness trajectories with slower growth independent of the number of competing strains, which we quantify with powerlaw fits and understand mechanistically via a random walk model that neglects dynamical correlations between genes. We show that increasing the level of clonal interference leads to fitness trajectories with faster growth (in functional form) without microscopic epistasis, but leaves the rate of growth invariant when epistasis is sufficiently strong, indicating that the role of clonal interference depends intimately on the underlying fitness landscape. The simulation package for this work may be found at https://github.com/nmboffi/spin_glass_evodyn.