Glassy Dynamics in Models of Confluent Tissue with Mitosis and Apoptosis

TL;DR

This study investigates how cell division and death influence the glassy, slow-moving behavior of confluent tissues, revealing that low cell cycling rates can lead to glass-like states despite ongoing cell turnover.

Contribution

The paper introduces a simulation-based analysis showing that low cell cycling rates can produce glassy tissue dynamics, reconciling experimental observations with theoretical models.

Findings

Glass-like tissue states occur at low cell cycling rates.

Cell birth and death rates contribute additively to tissue fluidity.

Subdiffusive cell displacements indicate caging behavior.

Abstract

Recent work on particle-based models of tissues has suggested that any finite rate of cell division and cell death is sufficient to fluidize an epithelial tissue. At the same time, experimental evidence has indicated the existence of glassy dynamics in some epithelial layers despite continued cell cycling. To address this discrepancy, we quantify the role of cell birth and death on glassy states in confluent tissues using simulations of an active vertex model that includes cell motility, cell division, and cell death. Our simulation data is consistent with a simple ansatz in which the rate of cell-life cycling and the rate of relaxation of the tissue in the absence of cell cycling contribute independently and additively to the overall rate of cell motion. Specifically, we find that a glass-like regime with caging behavior indicated by subdiffusive cell displacements can be achieved in systems with sufficiently low rates of cell cycling.

Equations

$$E_{Shape}=\sum _a{E}_a=\sum _a\left[{\kappa }_A{\left(A_a-A_{0a}\right)}^2+{\kappa }_P{\left(P_a-P_{0a}\right)}^2\right],$$

$$\hat{\mathbf{n}}({\theta }_a)=\cos ({\theta }_a)\hat{\mathbf{x}}+\sin ({\theta }_a)\hat{\mathbf{y}},$$

$${\partial }_t{\theta }_a={\eta }_a,$$

$$<{\eta }_a(t)>=0$$

$$<{\eta }_a(t_1){\eta }_{a^{\prime }}(t_2)>=2D_r\delta (t_1-t_2){\delta }_{a{a}^{\prime }}.$$

$${\partial }_t{\mathbf{r}}^{\mu }=-\frac{1}{\gamma }\frac{\partial }{\partial {\mathbf{r}}^{\mu }}{E}_{AVM},$$

$$E_{AVM}=E_{Shape}-\gamma v_0\sum _{cells-a}\hat{\mathbf{n}}({\theta }_a)\cdot {\mathbf{r}}_a,$$

$${\mathbf{r}}_a=\frac{1}{6A_a}\sum _{\nu }\left({\mathbf{r}}_a^{\nu }+{\mathbf{r}}_a^{\nu +1}\right)({\mathbf{r}}_a^{\nu }\times {\mathbf{r}}_a^{\nu +1})\cdot \hat{\mathbf{k}},$$

$${\partial }_t{\mathbf{r}}_a$$

$$=v_0\hat{\mathbf{n}}({\theta }_a)-\frac{1}{\gamma }\frac{\partial }{\partial {\mathbf{r}}_a}{E}_{shape},$$

$$\mathcal{O}(t)=\frac{1}{N_{cells}}\sum _a\Theta \left(b-|{\mathbf{r}}_a(t)-{\mathbf{r}}_a(0)|\right),$$

$$\frac{1}{{\tau }_{\alpha }}=\frac{1}{{\tau }_{\alpha 0}}+\frac{C_2}{{\tau }_{\delta }},$$

$$\delta u_{(m,a)}(r{)}^2=|{\mathbf{w}}_{(m,a)}{|}^2(r)-|{\mathbf{w}}_{(m,a)}{|}^2(\infty ),$$

$${\mathrm{D}}_{\delta }={\mathrm{D}}_{(a)}^{mean}+{\mathrm{D}}_{(m)}^{mean}+{\mathrm{D}}_{(a)}^{fluc}+{\mathrm{D}}_{(m)}^{fluc}.$$

$${\mathrm{D}}_i=\frac{N_{cells}}{4{\tau }_{\delta }}{⟨d^2⟩}_i,$$

$${⟨d^2⟩}_i=\frac{1}{A_T}\int {\mathrm{d}}^2\mathbf{x}{d}_i^2(r,\theta )$$

$${\partial }_t{\mathbf{r}}^{\mu }=-\frac{1}{\gamma }\frac{\partial }{\partial {\mathbf{r}}^{\mu }}E,$$

$$E=E_{shape}-\gamma v_0\sum _{cells-a}\hat{\mathbf{n}}({\theta }_a)\cdot {\mathbf{r}}_a,$$

$${\mathbf{r}}_a=\frac{1}{6A_a}\sum _{\nu }\left({\mathbf{r}}_a^{\nu }+{\mathbf{r}}_a^{\nu +1}\right)({\mathbf{r}}_a^{\nu }\times {\mathbf{r}}_a^{\nu +1})\cdot \hat{\mathbf{k}},$$

$$A_a=\frac{1}{2}\sum _{\nu }({\mathbf{r}}_a^{\nu }\times {\mathbf{r}}_a^{\nu +1})\cdot \hat{\mathbf{k}}.$$

$${\partial }_t{\mathbf{r}}^{\mu }$$

$$+v_0\left(\frac{\partial (\hat{\mathbf{n}}({\theta }_a)\cdot {\mathbf{r}}_a)}{\partial {\mathbf{r}}^{\mu }}+\frac{\partial (\hat{\mathbf{n}}({\theta }_b)\cdot {\mathbf{r}}_b)}{\partial {\mathbf{r}}^{\mu }}+\frac{\partial (\hat{\mathbf{n}}({\theta }_c)\cdot {\mathbf{r}}_c)}{\partial {\mathbf{r}}^{\mu }}\right).$$

$$\frac{\partial }{\partial {\mathbf{r}}^{\mu }}{E}_a=-\frac{{\Pi }_a}{2}({\hat{\mathbf{n}}}_{ab}{l}^{\mu \gamma }+{\hat{\mathbf{n}}}_{ac}{l}^{\mu \lambda })-T_a({\hat{\mathbf{l}}}^{\mu \gamma }+{\hat{\mathbf{l}}}^{\mu \lambda }),$$

$$T_a=\frac{\partial E_{shape}}{\partial P_a}{\Pi }_a=-\frac{\partial E_{shape}}{\partial A_a}.$$

$$\frac{\partial (\hat{\mathbf{n}}({\theta }_a)\cdot {\mathbf{r}}_a)}{\partial {\mathbf{r}}^{\mu }}=\frac{\partial {\mathbf{r}}_a}{\partial {\mathbf{r}}^{\mu }}\cdot \hat{\mathbf{n}}({\theta }_a)$$

$$\frac{\partial r_{a,j}}{\partial r_i^{\mu }}$$

$$+(r_j^{\lambda }+r_j^{\mu })R(\pi /2{)}_{ik}{r}_k^{\lambda }+(r_j^{\mu }+r_j^{\gamma })R(-\pi /2{)}_{ik}{r}_k^{\gamma }$$

$$\displaystyle-3(R(\pi/2)_{ik}r^{\lambda}_{k}+R(-\pi/2)_{ik}r^{\gamma}_{k})r_{a,j}\Bigg{]}\,,$$

$$R(\theta)_{ij}=\left[\begin{tabular}[]{cc}$\cos(\theta)$&$-\sin(\theta)$\\ $\sin(\theta)$&$\cos(\theta)$\end{tabular}\right]\,$$

$${\mathbf{f}}_i^{\mu (a)}$$

$$=\frac{v_0\hat{n}({\theta }_a{)}_i}{3A_a}[A_a^{\lambda \mu }+A_a^{\mu \gamma }]$$

$$\displaystyle+\frac{v_{0}}{6A_{a}}\Bigg{[}R(\pi/2)_{ik}\tilde{r}^{\lambda}_{k}\hat{n}(\theta_{a})_{j}(\tilde{r}_{j}^{\lambda}+\tilde{r}_{j}^{\mu})$$

$$\displaystyle+R(-\pi/2)_{ik}\tilde{r}^{\gamma}_{k}\hat{n}(\theta_{a})_{j}(\tilde{r}_{j}^{\mu}+\tilde{r}_{j}^{\gamma})\Bigg{]}\,,$$

$${\partial }_t{\mathbf{r}}^{\mu }=-\frac{1}{\gamma }\frac{\partial }{\partial {\mathbf{r}}^{\mu }}E+{\eta }^{\mu }$$

Full text

I Introduction

A number of experiments conducted over the past decade Angelini et al. (2011); Schotz et al. (2013); Park et al. (2015) have indicated the existence of glass-like states in confluent epithelial tissues. Evidence for glassy behavior includes the observation of dynamical heterogeneities characteristic of supercooled liquids Angelini et al. (2011) and caging effects in motile cell trajectories Park et al. (2015). In these living tissues, the constituent cells will often undergo mitosis (cell division) and apoptosis (programmed cell death) in a regulated cell life cycle. These cell life events are necessary for the survival and function of many tissues Renehan et al. (2001); Thompson (1995), and they clearly differentiate these dense cellular systems from the molecular and polymeric glasses that are so well studied. Introducing and removing material elements from any system should impact the dynamics, and it is therefore important to examine the effect of cell division and cell death on glassy tissue dynamics.

Recent simulations have suggested that the presence of cell division and death should generically fluidize a tissue Matoz-Fernandez et al. (2017a, b); Ranft et al. (2010); Malmi-Kakkada et al. (2018). Employing a particle-based model, for instance, Matoz-Fernandez et. al. Matoz-Fernandez et al. (2017a) have suggested that, even at small rates $k_{\delta}$ of cell division and death per cell, one will not observe the sub-diffusive behavior of cell trajectories, which we take to be a hallmark of glassy systems. This, therefore, poses a serious question for the experimental and theoretical studies of these systems: why are caging behavior and other signatures of glassy dynamics observed at all?

Earlier work by Ranft et. al. Ranft et al. (2010) at the continuum level considered the effect of cell division and death events on a model of 3-dimensional elastic tissue. Their analytical and numerical findings suggest that cell dynamics will be diffusive and controlled by $k_{\delta}$ in the limit of long times and infinitely large tissues. However, both this work and many other previous investigations have focused on cell division and death as the dominant driving force behind tissue dynamics. Epithelial cells can also escape the cages set by their neighbors through motility driven by traction forces the cells can exert on a substrate. Such forces play important roles in the development and function of organisms Puliafito et al. (2012); Martin et al. (2004); Trepat et al. (2009); Giavazzi et al. (2017, 2018); even in the absence of cell division and cell death, the timescale over which cells may escape their cages changes by orders of magnitude with variation of propulsive forces and intercellular tensions Henkes et al. (2011); Berthier and Kurchan (2013); Berthier (2014); Bi et al. (2016); Giavazzi et al. (2018). The full dynamics of tissues which experience motile forces in addition to mitosis and apoptosis should therefore have contributions from each of these sources of internal driving. In particular, it is natural to expect that the observed relaxation time of a tissue will be set by a competition between the timescale of cell cycling events and the timescale of cage escape due to the motile cells encountering mechanical energy barriers.

A complicating factor in previous investigations of model tissues with both active dynamics and cell cycle events is that they have focused on particle-based modes and they often couple the details of the cell cycle events to the current configurational state of the tissue Matoz-Fernandez et al. (2017a, b); Ranft et al. (2010) (for instance, triggering cell divisions when a particle is surrounded by a sufficient amount of space, or triggering cell deaths when the local density or contact number is too high). While such couplings may be biologically motivated, they also make it difficult to separate out the mechanically fluidizing effects of cell division and death on the glassy behavior of the model tissues.

In this work, we focus on models of confluent tissues, where there are no gaps or overlaps between cells. Such tissues are well-described by so-called “vertex models” that have been used extensively to understand patterning and rigidity in tissues Farhadifar et al. (2007); Hufnagel et al. (2007); Staple et al. (2010); Bi et al. (2014, 2015); Merkel et al. (2017). Although there has been substantial recent work on “Voronoi” versions of such models Bi et al. (2016); Barton et al. (2017), here we introduce an “Active Vertex Model” (AVM) which naturally incorporates cell motile forces similar to the previously studied Self-Propelled Voronoi (SPV) model, but without the restriction that cell shapes be described by Voronoi volumes.

The remainder of this work is organized as follows. In Sec. II we describe our Active Vertex Model and characterize its dynamical behavior in the absence of cell division and cell death. In Sec. III we extend this model to include a single controlled rate of mitosis and apoptosis, and we demonstrate that the phase diagram for this model has two dynamical regimes: one where the system is dominated by cell cycle events and another where the system is controlled by the inherent dynamics of the motility-only model. These combined regimes are well-described by a simple ansatz in which the overall relaxation time is an additive independent sum of the inherent rate of cell motility and the rate of cell cycling, revealing a universal crossover between cell-cycle-dominated and motility-dominated dynamics. In Sec. IV we investigate the displacement fields associated with individual mitosis and apoptosis events. Each event produces an Eshelby-like displacement field of tissue cells, in addition to a mechanical noise stemming from the disordered geometry of the tissue. Quantifying these displacements allows us to construct a prediction for the cell diffusion across a broad range of cell cycling rates, confirming that the effect of cell cycling on diffusion is similar to that predicted by previous continuum models Matoz-Fernandez et al. (2017a); Ranft et al. (2010), although the fluctuating motions generated by disordered geometries that are typically ignored by those models also contribute to overall cell motion. We close in Sec. V with a discussion of how our results compare to previous work, and we suggest experiments to guide future developments.

II A Vertex Model for Motile Tissues

Numerical studies of the influence of motility on tissue dynamics have been conducted in particulate models Henkes et al. (2011), the Potts model Kabla (2012); Chiang and Marenduzzo (2016), and Voronoi / vertex models Bi et al. (2016); Giavazzi et al. (2018); Sussman and Merkel (2018); Sussman et al. (2018). These investigations have generally indicated that motile forces tend to promote fluidity and enhance diffusion. For instance, using the SPV model – in which cells are described as polygons obtained from a Voronoi tiling of the plane that self-propel analogously with self-propelled particle models – Bi et. al. Bi et al. (2016) obtained a transition line separating glassy states from fluid-like states where cells frequently exchange neighbors.

While the Voronoi description is theoretically appealing, the restriction of cell shapes to local Voronoi volumes is not always desirable. More seriously, Voronoi models are not easily extendable to describe tissues with free boundaries, as would be needed to describe wound-healing geometries in 2D or isolated spheroids in 3D. To avoid these issues, while simultaneously preserving the spirit of the self-propelled Voronoi model, we have constructed a model for a motile tissue that preserves the freedom of polygonal cell shapes. In this work we restrict our attention to the generic case where every vertex is three-fold coordinated, although this restriction could be easily relaxed to study the sorts of open boundaries described above. To our knowledge, the dynamics of the Active Vertex Model presented herein has not been explored in previous work. Therefore, we include a more detailed description of the model, as well as quantification of its dynamical states in the absence of cell cycle events.

II.1 Shape forces in the Vertex Model

In the Vertex Model Farhadifar et al. (2007); Staple et al. (2010); Hufnagel et al. (2007); Bi et al. (2014) cells are modeled as irregular polygons tiling the plane, but in contrast to Voronoi models, the degrees of freedom are the positions of the vertices of the spatial tiling. Previous studies using this model have typically involved searching for geometric states which minimize the tissue energy,

[TABLE]

where the sum runs over all cells $a$ and $\kappa_{A}$ and $\kappa_{P}$ are elastic moduli. This energy drives the area $A_{a}$ and perimeter $P_{a}$ of cell $a$ toward a target area $A_{0a}$ and a target perimeter $P_{0a}$ respectively. The area term in Eq. 1 comes from the resistance to fluctuations in cell height. The perimeter term may be seen as arising from the competition between cell-cell adhesive interactions and the non-linear tension created by the cortical actin network. While the target shape properties may vary from cell to cell, in the absence of cell division and death we presume cells identical with uniform $A_{0}$ and $P_{0}$ values. An important parameter in such models is the dimensionless target “shape index” $s_{0}=P_{0}/\sqrt{A_{0}}$, which controls rigidity in both static Bi et al. (2015) and motile Bi et al. (2016) tissues.

II.2 The Active Vertex Model

Although the simplest implementation of activity in a vertex model allows each vertex to be self-propelled, we would like to model the behavior of a polarized motile cell moving persistently along a specific direction. Therefore, we extend this model to include self propulsion of cell $a$ in the direction

[TABLE]

where the cell propulsion angle $\theta_{a}$ is governed by

[TABLE]

with $\eta_{a}$ being a white gaussian noise defined by

[TABLE]

and

[TABLE]

Assuming that cell dynamics take place in the overdamped limit, vertex $\mu$ will follow an equation of motion,

[TABLE]

where $\gamma$ is the substrate friction. Here and throughout we use superscripted greek letters to refer to vertex positions, and subscripted latin letters to refer to cell positions. The total effective energy,

[TABLE]

now captures the cellular self-propulsion forces (of magnitude $v_{0}$) in addition to the standard shape energy from Eq. 1. In our AVM simulations we define $\mathbf{r}^{a}$ as the cell centroid, or geometric center,

[TABLE]

where $\nu$ here indexes the $N_{a}$ vertices on cell $a$ in counterclockwise fashion and $\mathbf{r}_{a}^{N_{a}+1}=\mathbf{r}_{a}^{1}$.

With these equations of motion, each cell attempts to carry out a persistent random walk, as is commonly observed for cells isolated from the surrounding tissue and in the absence of external signals Metzner et al. (2015); Passucci et al. (2019). While it is uncommon to think of such a self-propulsion force as coming from an energy Sknepnek and Henkes (2015), writing things in this way permits an analogy with the Voronoi model. In the SPV, the degrees of freedom are the Voronoi centers and we might imagine constructing the same energy Eq. 7, using these Voronoi centers as the $\{\mathbf{r}_{a}\}$. It is easy to see that using these as the degrees of freedom in the overdamped equation of motion,

[TABLE]

indeed leads to the standard SPV dynamics Bi et al. (2016); Yang et al. (2017); Barton et al. (2017). We are therefore using the closest energetic analog that implements SPV dynamics in a vertex model. In the AVM we will instead apply shape-based and motility-based forces to the vertices. These forces may be precisely defined by taking the derivatives in Eq. 6, as explained in greater depth in Appendix A.1.

II.3 Glassy dynamics in the Active Vertex Model

Since this model for confluent tissue dynamics has not been explored in previous work, we first examine its behavior and search for glassy states. We quantify the dynamical state of the system by using two familiar metrics: the mean square displacement and the self-overlap function. The self-overlap function is defined as

[TABLE]

where $\Theta$ is the heaviside function and $b$ represents the size of a typical cage in natural units, which we set to $0.5$. All lengths are measured in units of $1/\sqrt{\rho}$, where $\rho$ is the system number density. The function $\mathcal{O}$ has a value of $1$ at $t=0$ and decays towards zero as cells move beyond the caging distance $b$. The structural relaxation time $\tau_{\alpha 0}$ is defined as the point at which the self-overlap decays below $1/e$.

As indicated in the introduction, in this work we use the term “glassy” to primarily refer to the presence of a subdiffusive regime in the mean squared displacements, and to the observation of uncaging times $\tau_{\alpha 0}$ that are orders of magnitude longer than expected from free diffusion. The uncaging time averaged over 10 systems of $N_{cells}=300$ is displayed in Fig. 1. Its behavior is consistent with previous work on the Self-Propelled Bi et al. (2016) and Thermal Sussman et al. (2018) Voronoi models. Less motile tissues at lower $s_{0}$ (higher cortical tension) are indeed glassy and have long relaxation times. The tissue may be effectively fluidized by either increasing $v_{0}$ (higher effective temperature) or increasing $s_{0}$ (higher cell-cell adhesion).

Although the overlap function shown in Fig. 1 appears qualitatively similar to those observed in other glassy materials, both it and the mean-squared displacement shows an unusual behavior. In most glassy materials, the mean-squared displacement levels off to a subdiffusive regime where $\Delta r^{2}\sim t^{0}$, followed later by a diffusive regime $\Delta r^{2}\sim t^{1}$. Here we find anomalous sub-diffusive behavior with a non-standard exponent, $\Delta r^{2}\sim t^{\alpha}$, where $0<\alpha<1$, is found to persist across a wide range of timescales. In this same parameter regime, we also see multiple timescales in the shear stress autocorrelation function. These results are highlighted in Appendix A.2. Such properties are highly atypical in the context of standard glasses, and may be related to the unusual zero-temperature rigidity transition in such models Sussman and Merkel (2018). Understanding the origins of this behavior will be an interesting avenue for future work.

As we are interested in characterizing the simple long-time diffusion of cells in our model, this unusual behavior is problematic. To acknowledge this, throughout the manuscript we highlight with distinct symbols any data presented from simulations that remain sub-diffusive on long timescales. In Section III we will see that the parameter regimes of the active vertex model that are diffusive at long times have dynamics which are easily interpreted as a competition between glassy relaxation and relaxation driven by cell division and death.

III Simulation of mitosis and apoptosis events

We would now like to incorporate the influence of cell cycle events in these vertex model simulations. To identify an appropriate procedure, we first note that division and death events can change the total number of cells and the cell areal density, with potentially drastic changes in tissue dynamics. To minimize such effects, we choose to work in a constant-number ensemble. In our simple implementation, each instance of apoptosis will be accompanied by a mitosis event somewhere else in the tissue. In addition, the implementation of the individual cell cycle events is chosen to preserve the sum of the cell target areas.

In apoptosis, a cell will abruptly contract to a small size and then extrude itself, effectively disappearing from the 2d monolayer Rosenblatt et al. (2001). Shown in Fig.2(d-f), our simple realization of apoptosis on cell $a$ begins with setting the target area and perimeter to zero. This change induces the rapid contraction of cell $a$ to small size and generally will lead to a final triangular shape. The simulation then detects triangular cells that are smaller than some threshold area $A_{min}$ and deletes them.

Similarly, in the process of mitosis a cell will expand, eventually reaching a threshold size and dividing into two complete cells Conlon and Raff (2003); Dolznig et al. (2004); Streichan et al. (2014), all the while maintaining tissue cohesion Baker and Garrod (1993). While this process typically spans a much longer time period than the event of a cell death, we model this growth process as similarly instantaneous to suppress fluctuations of tissue density. Shown in Fig 2, in the division of cell $a$ two non-adjacent edges of the cell are chosen at random, and their midpoints determine the axis of the division. A vertex is then added to the center of each edge and these new vertices are connected with a new edge. Naturally the product is two closed polygonal cells out of one, and the shape parameters of these new cells are then set to the values of the parent cell. These cells are then allowed to expand dynamically in the simulation. Similar division dynamics were studied in a quasi-static system by others, including Farhadifar et. al. Farhadifar et al. (2007).

Such constant number simulations allow us to enforce a global cell cycle event rate $k_{\Delta}$, which is implemented as a Poisson process. This sets the effective cycle rate per cell as $k_{\delta}=k_{\Delta}/N_{cells}$. Note that $k_{\delta}$ is a more appropriate parameter for real biological systems whose cell cycle timing does not depend strongly on the extent of the surrounding tissue. We therefore use $k_{\delta}=1/\tau_{\delta}$ as our tuning parameter.

In simulations with cell division and death, complications arise when trying to extract cell trajectory data. This is because many trajectories will start and end during the course of the simulation. Rather than parse this network of trajectories, we follow the method used in Matoz-Fernandez et al. (2017a) and exclude a small (%10) subset of the cells from the cell cycle events. Using these “tracer” cells to obtain our dynamical data is both convenient and affords the closest comparison with the previous literature Matoz-Fernandez et al. (2017a).

We note that fast rates of cell division tend to cause tissues to reach rare topological states that can generate a breakdown in our simulation framework; to avoid this we restrict the data presented here to $\tau_{\delta}>3000$.

III.1 Tissue dynamics in the presence of mitosis and apoptosis

Given this simulation protocol, we search for signatures of glassy behavior. Tissue dynamics is quantified using the uncaging time $\tau_{\alpha}$, defined as the time required for the observed self-overlap to decay below $1/e$. The notation intentionally differs from that in Section II.3, to distinguish $\tau_{\alpha 0}$ as the uncaging time in the limit of the model without cell cycle events (i.e., the limit $k_{\delta}=0$).

Our search for glassy behavior is guided by the simple expectation that the dynamics is determined by a competition between the timescale $\tau_{\delta}$ of division and death events and the timescale of motility driven cage escape $\tau_{\alpha 0}$. Naturally, we recover the dynamics of the “bare” (free of cell division and death) AVM in the limit $\tau_{\delta}\gg\tau_{\alpha 0}$ where mitosis and apoptosis will play a negligible role. One may then expect to find the opposite behavior in the limit $\tau_{\delta}\ll\tau_{\alpha 0}$, where motility-based dynamics should become negligible. To check this, we increase $s_{0}$ to decrease $\tau_{\alpha 0}$. In the top-left region of Fig. 3-a where $\tau_{\alpha 0}$ is large and $\tau_{\delta}$ is small ($k_{\delta}$ is high) we observe structural relaxation times $\tau_{\alpha}$ which are approximately independent of $s_{0}$ and by proxy, $\tau_{\alpha 0}$.

III.2 Simple model for interaction between $\tau_{\delta}$ and $\tau_{\alpha 0}$

Building on the limiting behaviors described above we develop a very simple ansatz for the interplay between cell death and division and glassy dynamics. Specifically, we assume that the overlap decay rate $\frac{1}{\tau_{\alpha}}$ of a dividing tissue is determined by the weighted sum of the bare cage escape rate $\frac{1}{\tau_{\alpha 0}}$ and the rate of division and death events according to

[TABLE]

where $C_{2}$ captures the displacements of surrounding cells resulting from the division and death events. The strong assumption we have made here is that the two rates add in series and are not strongly correlated. Our model predicts that the quantity $\frac{\tau_{\alpha}}{\tau_{\alpha 0}}$ will be a function only of $\frac{\tau_{\delta}}{\tau_{\alpha 0}}$. As shown in Fig. 4, this works reasonably well, indicating the generic presence of a regime where fast divisions dominate the dynamics and $\tau_{\alpha}$ is proportional to $\tau_{\delta}$ (left side of Fig. 4), along with a regime of slow divisions where $\tau_{\alpha}$ becomes independent of $\tau_{\delta}$ (right side of Fig. 4). To understand what controls the crossover, we quantify the impact of cell division and cell death events on the motion of surrounding cells in the next section.

IV Flow and fluidization from individual mitosis and apoptosis events

To quantify the effect of a single cell cycling (division or death) event on the motion of the surrounding tissue, we perform a special set of “single-event” simulations. In these, the standard AVM is run for a short equilibration time before a single cell is chosen at random to undergo either apoptosis or mitosis in a randomly chosen direction. We then monitor the motion of the surrounding tissue cells and average over 1000 realizations.

We measure the individual cell displacement vectors $\{\mathbf{u}_{i}\}$ and use them to construct a few useful quantities. In order to resolve coherent spatial data, the cells are first binned based on their distance from the event (and in the case of cell division based on their angle relative to the division axis orientation). Within each bin, we calculate two quantities: the vector averaged displacement $\mathbf{u}_{(m,a)}(r,\theta)$ and the vector standard deviation of this average, $\bm{w}_{(m,a)}(r,\theta)$, where the subscripts $m$ and $a$ refer to mitosis and apoptosis events, respectively.

As noted by Puosi et. al. Puosi et al. (2014) for the similar case of shear transformations in thermal sphere packings, the vector averaged $\mathbf{u}_{(m,a)}(r,\theta)$ captures the mean elastic response of the surrounding medium. A straightforward calculation of the response of a homogeneous elastic medium to localized strains Puosi et al. (2014) suggests that the mean deformations associated with either a cell death or cell division event should fall off as $r^{-1}$ in two dimensions. This is consistent with the numerical observations in our model, shown in Fig. 5 (a), (b) and (c). In particular, the dashed line in Fig. 5(c) shows the expected scaling of $r^{-1}$, which is in reasonable agreement with the data.

Additional cell displacements are expected to arise in the vicinity of the event due to disorder in the tissue structure. These additional fluctuating displacements are captured by $\bm{w}_{(m,a)}(r,\theta)$. It is important to note that this quantity captures both local contributions from the mitosis (apoptosis) event, as well as contributions everywhere from active, motility driven, cell motions. To separate the local contribution of these fluctuating elastic displacements from the noise generated by the active forces in the AVM, we define

[TABLE]

where $|w_{(m,a)}|^{2}(\infty)$ is the far-field plateau in the fluctuating displacement field generated by active noise, shown in the inset to Fig. 5 (d). The resulting rotationally averaged fluctuating displacement field capturing only the motion caused by the cell cycling event, $\delta u_{(m,a)}(r)$, is plotted in Fig. 5(d) for a mitosis event.

We note immediately from the data in Fig. 5(c) that the mean deformations show only minor variation with the change in $s_{0}$ and therefore depend little on the quiescent dynamics. In contrast, in Fig. 5(d) we see that the fluctuating motions decrease as we go from more rigid tissue at low $s_{0}$ to more fluid tissue at $s_{0}\sim 3.9$. The degree of variation is, however, still minor and in general the dynamical response to cell division and death events is not expected to be a good indicator of the underlying dynamical state of the tissue.

IV.1 Predicting Cell Displacements from Mitosis and Apoptosis

In the simple model in Eq. 11 it is assumed that motility and cell cycling contribute independently and linearly to the rate of cell motion. This same assumption implies that the diffusion $\mathrm{D}$ will be the sum of independent parts $\mathrm{D}=\mathrm{D}_{0}+\mathrm{D}_{\delta}$, where $\mathrm{D}_{0}$ is the contribution from motile forces and $\mathrm{D}_{\delta}$ is the contribution from the cell cycling. Below, we estimate $\mathrm{D}_{\delta}$ from the displacements measured for single cell-death and cell-division events.

As pointed out by Ranft et. al. Ranft et al. (2010), the mean displacements identified above produce a constantly changing reference state in an elastic material. Therefore, while cells may oscillate in their cages, the cages themselves move as a result of each event via the displacement fields $\mathbf{u}_{(m,a)}(r,\theta)$. We would like to estimate the effect of this changing reference state on the mean square displacement of our tracer cells. These events will simultaneously create fluctuating displacements, which will also contribute to the net cell motion.

To estimate the total diffusion produced by these different contributions (mitosis and apoptosis, mean and fluctuating displacements), we again make the simplifying assumption that they are uncorrelated. We can then write $\mathrm{D}_{\delta}$ as a sum of the decoupled parts

[TABLE]

Here, $\mathrm{D}^{mean}_{(m,a)}$ captures the diffusion due to the mean $\mathbf{u}_{(m,a)}(r,\theta)$, while $\mathrm{D}^{fluc}_{(m,a)}$ captures diffusion due to the fluctuating $\delta u_{(m,a)}(r,\theta)$. Quantifying each of these diffusion constants will rely on the simplifying assumptions (1) that we may ignore randomness in the timing of the divisions and deaths and (2) that the effect from each event is felt instantaneously by the surrounding tissue. This enables us to simply sum up the contributions from the $n(t)=tN_{cells}/\tau_{\delta}$ events that will have taken place after a time $t$. Therefore, each term in Equation 13 may be estimated in terms of an average magnitude of motion per cell per event $d_{i}$ as

[TABLE]

where the four contributions from the rhs of Eq. 13 are now indexed by $i$. As shown in Appendix C, the values of $\left\langle d^{2}\right\rangle_{i}$ may be then estimated from summation of the single-event data shown in Fig. 5(b) and (c). In practice, $\left\langle d^{2}\right\rangle_{i}$ is computed as a spatial average

[TABLE]

where $d^{2}_{i}(r,\theta)$ represents the square displacements captured in our previously averaged values $\mathbf{u}^{2}_{(a,m)}(r,\theta)$ and $\delta u^{2}_{(a,m)}(r,\theta)$. As these simulation data are measured on a radial grid in discrete bins, the integral in Eq. 15 will turn into a sum over bins, and will have a cutoff radius at $r\sim L/2$. While this excludes cell displacements in the corners of the simulation box, this contribution is small and the predictions constructed without these corners validate their exclusion. As shown in Fig. 6-(a), each of these contributions to the diffusion $\mathrm{D}_{i}$ is essential to form a complete prediction of the observed diffusion. In Fig. 6 this prediction provides a reasonable estimate, deviating most significantly for tissues which display glassy signatures in the ($\tau_{\delta}\rightarrow\infty$ limit) bare-AVM. This breakdown is expected because the corresponding component $\mathrm{D}_{0}$ in our prediction is no longer truly a diffusion, since the AVM itself displays anomalous subdiffusion in this regime as discussed in Appendix II.

This prediction for the diffusion constant may be extended to a prediction for the uncaging time $\tau_{\alpha}$, which performs similarly well. Relating these two quantities relies on the simple result $\mathrm{D}\cdot\tau_{\alpha}=c_{B}$ which applies to free overdamped Brownian motion in 2D. The constant $c_{B}\sim 0.1363$ is straightforwardly derived as shown in Appendix. B. Using this relation, our prediction of $\tau_{\alpha}$ then relies only on the assumption that, on long time scales, the motion of cells is Brownian.

This is an assumption which, again, begins to break down as motion becomes glassy and cage-escapes become intermittent. In Fig. 4, the black dashed line prediction matches well with (red) data. However, even where anomalous subdiffusive glassy signatures arise, our prediction captures this crossover between motility-dominated and cell division/death dominated dynamics reasonably well.

V Discussion

To study the influence of propulsive forces on the dynamical state of tissues, we have developed an Active Vertex Model that naturally generalizes the Vertex Model framework to incorporate cell motility analogously with self-propelled Voronoi models Bi et al. (2016). In Section II, we have shown that this model qualitatively reproduces some of the dynamical results of the SPV model. As found previously by Bi et. al. Bi et al. (2016), we locate “glassy” states which exhibit sub-diffusive mean square displacements. The transition from a fluid-like state to a solid-like state is achieved either by decreasing the propulsive forces, or by increasing the effective cell-cell interfacial tension.

We then extended the AVM to investigate the influence of cell death and cell division on tissue dynamics. Our investigation was guided by the expectation that the rate of cell division and death will compete with the “bare” uncaging rate from the AVM to determine the tissue dynamics. This simple picture allows us to identify evidence of glassy states even at finite rates of apoptosis and mitosis. The glassy behavior found here appears to be in a different regime from the one explored in Matoz-Fernandez et al. (2017a), where sources of activity beyond cell cycling were not the focus. We again emphasize that another pertinent difference is the choice of whether or not to couple cell cycle events directly to the local structural state of the individual cells. This may resolve a standing discrepancy between existing experiments and theory. While cell division and death events will generally lead to diffusive fluid-like behavior on the longest timescales, subdiffusive behavior can be observed on intermediate timescales when the bare uncaging rate is slow enough to allow it and the cell division rate is low enough to not obscure it. The subdiffusive dynamics also suggests, by analogy with supercooled fluids, that these tissues will behave elastically on these intermediate timescales.

We further characterized the spatial distribution of cell displacements in the vicinity of each individual mitosis and apoptosis event. As expected from previous work for the case of area-preserving deformations Puosi et al. (2014), we have identified an average displacement field that matches the one obtained from elasticity theory, as well as “fluctuating” displacements that arise from the disordered geometry. The mean elastic response decays spatially as $1/r$, with a magnitude that does not depend significantly on tissue dynamics. The fluctuating displacements, in contrast, vary in magnitude with model parameters, such as cell interfacial tension. Specifically, mitosis and apoptosis events tend to generate more rearrangements in rigid tissues than in fluid ones.

Using the calculated mean and fluctuating displacements resulting from these division and death events, we assess the relative importance of motile forces versus division and death in the tissue dynamics. From this data, we are able to construct predictions for both the diffusion coefficient and the uncaging time. These predictions are very accurate for systems where true diffusive exponents are measured and begin to break down as subdiffusive glassy signatures emerge. Even in these glassy states, however, the prediction gives a good qualitative description of the crossover from motility-dominated to cell cycle-dominated dynamics.

Taken together, our results provide a general perspective on the role of cell division and cell death in tissues. To reasonably good approximation, cell cycle events appear to produce long-time diffusivity with adds linearly with the existing cell dynamics. This linear behavior appears independent of the approximate cell lifecycle $2\tau_{\delta}$ for the values probed here, but is expected to break down as this lifetime gets faster than the speed of sound in the tissue and events begin to interact.

In addition, we have not addressed the question of whether glassy signatures persist as the tissue size is increased. In Appendix C.1 we present evidence that the total diffusion due to mean displacements scales as $\log(N)$ for our two-dimensional simulations. In disordered systems we expect, however, that the signal generated by the mean displacements will destructively interfere with an effective noise generated by motions from other events at a characteristic lengthscale, as has been previously highlighted for sheared particulate packings Lemaître and Caroli (2009). Beyond this lengthscale, we no longer expect the $\log(N)$ scaling to hold, and so the system will remain glassy in the thermodynamic limit. While we do not identify such a lengthscale for the simulations performed here, testing this hypothesis is an important avenue for future investigations.

So far, we have focused on dynamics that explicitly keeps the tissue at a fixed number density, in order to pinpoint the effect of cell divisions on tissue dynamics. However, in many experiments on epithelial monolayers, a common and interesting regime is where cell divisions outnumber cell deaths and so the number density increases rapidly. An important open question is how Active Vertex Model parameters, such as $s_{0}$ or $v_{0}$, change during epithelial densification. Several experiments on pairs or triplets of cells indicate that many epithelial cells exhibit “contact inhibition of locomotion”, where cells reduce traction forces and slow down upon contact Puliafito et al. (2012); Abercrombie (1970). Such behavior could be modeled by decreasing $v_{0}$ as a function of time in densifying monolayers, and additional experiments to measure traction forces and fluctuations during such processes would be very valuable in constraining model predictions. Perhaps even more interestingly, it remains unclear how cells regulate their shape $s_{0}$ during such processes. One possibility is that cells attempt to maintain the same shape despite changes to density, meaning that they must change their perimeter to match their decreasing area as a function of time. Alternatively, one could postulate that cells keep their preferred perimeter fixed while their area decreases dramatically, resulting in an increasing $s_{0}$ that would fluidize a tissue. This latter hypothesis seems inconsistent with existing experimental data, where tissues are generically found to solidify as they become more dense. Careful studies of cell shape coupled with careful studies of interfacial tensions Saraswathibhatla and Notbohm (2019) would help to constrain models and better test model predictions.

Finally, it is interesting to think about how cell division and death rates are themselves affected by tissue dynamics. For example, it is known that in many cell types cell division rates are governed by the magnitude of local stresses that build up in the tissue Shraiman (2005); Streichan et al. (2014). In addition, the orientation of a cell division is also controlled by local stresses. This creates the possibility for interesting feedbacks, where tissue dynamics is controlled by the rate of cell division and the rate of cell division is controlled by tissue dynamics. Understanding precisely how division and death affect dynamics is therefore very important for predicting how such feedback loops can control tissue growth and patterning. It would be interesting to explore how such feedback loops generate patterns in AVMs, and compare to experiment.

Acknowledgements

This work was supported primarily by Simons Foundation Grants Number 446222 and 342354 (MC and MCM). Partial support also came from: NSF-DGE-1068780 (MC), NSF-DMR-1609208 (MCM and MC), NSF-DMR-1352184 (MC and MLM), NSF-PHY-1607416 (MLM and DMS) as well as Simons Foundation Grant Number 454947 (MLM and DMS).

VI Appendices

Appendix A The Active Vertex Model

Here we include a more detailed description of the Active Vertex Model implemented in this work.

A.1 Forces in the Active Vertex Model

In this Active Vertex Model (AVM), the degrees of freedom are the positions of the vertices of the spatial tiling. The vertex labeled $\mu$ follows an overdamped equation of motion

[TABLE]

where $\gamma$ is the substrate friction. The effective energy

[TABLE]

captures the self-propulsion force (of magnitude $v_{0}$) of each cell $a$ in direction $\theta_{a}$ in addition to the standard shape energy terms from Eq. 1. The geometric center (centroid) $\mathbf{r}^{a}$ of cell $a$ defined by

[TABLE]

captures the center of mass of a polygon of uniform mass density. Here, $\nu$ indexes the $N_{a}$ vertices on cell $a$ in counterclockwise fashion and $\mathbf{r}_{a}^{N_{a}+1}=\mathbf{r}_{a}^{1}$. The polygon area may be expressed in similar terms as

[TABLE]

The force (and therefore motion) on each vertex can be calculated by carrying out the derivatives in Eq. 6. Each vertex in this model is connected to 3 cells and each cell energy will contribute separate terms to the net motion of the vertex. If vertex $\mu$ is connected to cells $a$, $b$ and $c$, then the motion breaks down into

[TABLE]

For simplicity, we may focus on the contributions from cell $a$. As identified in previous work Yang et al. (2017), the shape energy produces tension-based and pressure-based forces on each vertex. The shape-based force on vertex $\mu$ from cell $a$ reads

[TABLE]

where $\gamma$ and $\lambda$ index the vertices of cell $a$ which are adjacent to $\mu$, $l^{\mu\gamma}$ ($l^{\mu\lambda}$) and $\mathbf{\hat{l}}^{\mu\gamma}$ ($\mathbf{\hat{l}}^{\mu\lambda}$) are the length and direction of the edge connecting vertex $\mu$ to vertex $\gamma$ ($\lambda$) and the unit vector $\mathbf{\hat{n}}_{ab}$ ($\mathbf{\hat{n}}_{ac}$) points across the edge shared by cell $a$ and cell $b$ ($c$) as in Fig 7. The tension and pressure of cell $a$ are respectively

[TABLE]

Similar to the above, we now evaluate derivatives to understand the self propulsion forces which act on the vertices. Again, we will consider only the contributions from cell $a$. The derivative will take the form

[TABLE]

With these expressions, taking the derivative in Eq. 23 is tedious but straightforward. This becomes

[TABLE]

where

[TABLE]

is a vector rotation by an angle $\theta$, $(i,j)$ index cartesian vector components in the $x-y$ plane and Einstein summation convention is assumed for repeated indices. While the expression in Eq. 24 is a bit unwieldy, we note that it is translationally invariant and may consider it in a more convenient coordinate system $\tilde{\mathbf{r}}=\mathbf{r}-\mathbf{r}_{a}$ which has its origin at $\mathbf{r}_{a}$. Using this, we find the force on vertex $\mu$ due to the motility of cell $a$ as

[TABLE]

where $A_{a}^{\lambda\mu}$ and $A_{a}^{\mu\gamma}$ are the area of triangles with vertices $\{\mathbf{r}_{a},\mathbf{r}^{\lambda},\mathbf{r}^{\mu}\}$ and $\{\mathbf{r}_{a},\mathbf{r}^{\mu},\mathbf{r}^{\gamma}\}$ respectively, as shown in Fig. 7. While this expression does not lend itself to much insight, we can see on inspection that the first terms will move the vertex in the direction of $\mathbf{\hat{n}}(\theta_{a})$. The second terms in the large square brackets will roughly serve to shrink the length of interfaces at the back of the cell, while expanding the length of interfaces at the front.

While these forces determine the continuous time-evolution of the vertices, there can be no large-scale migration of cells through the tissue until topological rearrangements are allowed. We therefore include a protocol for T1 rearrangements (as shown in Fig 8) which are sufficient to explore the space of cellular topological configurations at constant density. In practice, the edge lengths are periodically checked after a time $t_{T1}$ for values lower than a threshold $\ell_{c}$. These edges are then topologically rearranged so that the two cells which initially share the short edge are no longer in contact. The edge is then rotated by $\pi/2$ and the length extended by a factor $\lambda_{T1}$. The value of $t_{T1}=0.05$ used here is chosen for speed of simulation, while the rescaling factor $\lambda_{T1}=2$ is chosen to avoid “T1-traps” whereby the same transition may repeat itself regardless of energetic favorability. The cutoff length $\ell_{c}=0.04$ is chosen small enough to make the transition appear continuous but also large enough so that a vertex may “find” and activate the desired T1. The parameters $\ell_{c}$, $\lambda_{T1}$ and $t_{T1}$ have been separately varied within reasonable ranges and the impact on the dynamics appears insubstantial. These simulations are implemented using the “cellGPU” codebase previously developed by one of us Sussman (2017).

We note again that this form of the Active Vertex Model is chosen to minimize the differences with the Self-Propelled Voronoi model. The differences are limited to the following: (1) the Active Vertex Model has more degrees of freedom, avoiding the shape constraints of Voronoi tesselations, (2) T1 rearrangements in the AVM must be done by hand, while in the SPV they come about naturally and (3) motility in the AVM is designed to propel the centroid of the cell, whereas for the SPV this is replaced by the Voronoi center.

A.2 Subdiffusive Exponents and Exotic Dynamical States in the Active Vertex Model

While the main text describes the glassy behavior present in the AVM, Vertex Models governed by an energy of the type in Eq. 1 may also be subject to anomalous energy landscapes and non-trivial inherent stresses Sussman and Merkel (2018); Sussman et al. (2018). As a reflection of these features, dynamical signatures such as the mean-square displacement (MSD) and the self-overlap (as defined in Eq. 1) may appear anomalous in the vicinity of the glassy phase. Here, we explore these anomalies in greater detail.

The mean-square displacement in thermal particulate systems approaching a glass transition will display the signatures of “caging”. This manifests as a flat plateau in the MSD as particles reach this cage size (radius) and are temporarily trapped before diffusing away on a longer timescale Berthier and Biroli (2011); Weeks and Weitz (2002). The AVM MSD, in contrast, near the onset of glassiness does not appear to completely flatten at any discernible length scale. Instead, we observe an extended subdiffusive regime which emerges as the system approaches the glassy state, implying that there may not be a well defined cage-size and that the cells are exploring an unusual distribution of local metastable states as they attempt to realize their target shapes. It may therefore not be useful to think of the travel between these metastable states as ”cage-breaking” which come paired with intermittent jumps in displacement on the order of the particle size. Instead, it appears that a distribution of effective cages may be traversed by the cells in a more continuous manner.

The above points to the unconventional geometric mechanics of the Vertex Model as the cause of this anomalous behavior, but we must also consider the details of the non-equilibrium dynamics chosen here. As we have employed a non-zero persistence time of these forces as $\tau_{p}=1/\mathrm{D}_{r}=1$, we also perform simulations with the motility replaced by thermal forces on the vertices at effective temperature $T$. These thermal simulations are run in the AVM with $v_{0}$ set to zero, and the equation of motion of a vertex instead obeying

[TABLE]

where

[TABLE]

and

[TABLE]

Results for this thermal vertex model are summarized in Fig. 10 and indicate qualitatively the same behavior as observed for cell motile forces. We also note that previous studies in the Self-Propelled Voronoi and Thermal Voronoi models have not explicitly identified this subdiffusive dynamical feature. This indicates that this strange extended subdiffusion is associated with the Vertex Model geometric energy along with the vertex degrees of freedom that it acts on.

Finally, we investigate temporal fluctuations in the mechanical stress as defined instantaneously by

[TABLE]

as in Ref. Yang et al. (2017). From this global tissue stress, we compute the stress autocorrelation function via

[TABLE]

Where the angular brackets denote the average over both noise and over the time $t_{0}$. In practice, the computation of this correlation is accomplished using the multiple-time algorithm described by Ramírez et. al. in Ref. Ramírez et al. (2010). This stress auto-correlation function, shown in Fig. 10 for a range of $s_{0}$ values, indicates additional unusual behavior. Specifically, we do not resolve the traditionally expected “beta” (fast) and “alpha”(slow) decorrelation times. Instead, the stress appears to decorrelate in three or more stages across a range of timescales. This may also indicate that the process of cage breaking can occur across a broad range of different length scales (and therefore different time scales) in epithelial tissues.

Appendix B Relating the Uncaging Time to the Diffusion Constant for Brownian Particles

This work relies largely on two scalar measures of tissue dynamics, the uncaging time $\tau_{\alpha}$ and the diffusion coefficient $\mathrm{D}$. Here we show that, for the case of Brownian motion in two dimensions, these quantities may be related explicitly as

[TABLE]

where $b$ is the radius defining the cage size in the overlap function (Eq. 10). This equation (which is related to the Stokes-Einstein relation) breaks down as a material approaches a glass transition, and therefore provides important context for interpreting dynamical parameters.

To understand Eq. 32, we note that free Brownian motion of an ensemble of particles may be described by a smooth probability density function $P(\mathbf{x},t)$ which is governed by Fick’s Law

[TABLE]

To understand the measured uncaging time, we must first solve for the $P(\mathbf{x},t)$ for a single particle starting at an arbitrary point $\mathbf{x}_{0}$ at a time $t_{0}$. This corresponds to solving with an initial condition $P(\mathbf{x}_{0},t_{0})=\delta(\mathbf{x}-\mathbf{x}_{0})$. Preserving the normalization of $P$, this is solved for $t>t_{0}$ by

[TABLE]

This probability distribution may be used to take averages and we may therefore rewrite the self-overlap,

[TABLE]

as

[TABLE]

where, for convenience, we have set $t_{0}=0$ and $r=|\mathbf{x}-\mathbf{x}_{0}|$. The uncaging time is then defined by the condition $\mathcal{O}(\tau_{\alpha})=1/e$ and we may insert Eq. 34 into Eq. 36, integrate, and manipulate to find our desired result in Eq. 32. Evaluating this constant for the case $b=0.5$, we find that $\tau_{\alpha}\mathrm{D}=c_{B}\sim 0.1363$.

Appendix C Displacements due to single division or death events

Here we describe our method of converting displacement data for single apoptosis and mitosis events into averaged motions per event $\left\langle d^{2}\right\rangle_{i}$ as well as a predicted diffusion coefficient $\mathrm{D}_{\delta}$ per Equations 13 and 15.

C.1 Estimating mean contributions to displacement

To estimate the first and second diffusion coefficients on the right hand side of Eq. 13, we will use the ”single-event” mean motion data $\mathbf{u}_{(m,a)i}(r,\theta)$ obtained in Section IV. On the long timescales which we are interested in, each event effectively instantaneously moves the cells via the mean displacement $\mathbf{u}_{(m,a)i}(r,\theta)$ and the cells then follow an equation of motion

[TABLE]

in which the $n(t)$ displacements due to previous events are simply summed up. With this EOM, the MSD may be computed as

[TABLE]

where we have implicitly defined the average mean motion per event $\left\langle d^{2}\right\rangle^{mean}_{(m,a)}$ and $\langle\rangle$ represents the average over the realizations of these apoptosis and mitosis events. As the stochasticity here comes from the spatial positioning (and orientation) of the mitosis (apoptosis) event with respect to our tracer, finding this average for $\left\langle d^{2}\right\rangle^{mean}_{(m,a)}=\langle\mathbf{u}^{2}_{(m,a)i}\rangle$ requires integrating over the possible positions and orientations of this mitosis (apoptosis) event. This amounts to a spatial integral

[TABLE]

over the mean displacement data where $\rho$ is the global number density of cells. To directly produce a numeric value of $\left\langle d^{2}\right\rangle^{mean}_{(m,a)}$ from the data, the integral is turned into a sum over all the bins in $r$ and $\theta$. The diffusion coefficient produced by these motions is then defined from Eq. 38 as

[TABLE]

The estimate produced above leads to useful and accurate predictions of the dynamics as shown in Figures 6 and 4 in the main text. However, this prediction provides little insight into the functional form of these mean diffusions. Some insight can instead be gained by estimating the integral in Eq. 39 using some assumptions based on the data. In Fig. 12(a) we see that averaging these mean displacements over the angular bins reveals a $1/r$ trend which may be captured by $u_{(m,a)}(r)=u_{(m,a)}/(\rho r)$. This trend persists out to a distance of $\sim 0.3L$ where the event begins to interfere with its periodic image and drops off quickly (here $L=\sqrt{N_{cells}/\rho}$ is the periodic box length). A good approximation of the average displacement is given by

[TABLE]

where the small $r_{small}$ and large $r_{large}\sim 0.3L$ cutoffs of integration respectively capture the typical cell neighbor spacing and the extent of the mean field. Inserting the above form of $u_{(m,a)}(r)$ and integrating, we find an expression for the diffusion

[TABLE]

Inserting $r_{large}=0.3\sqrt{N_{cells}/\rho}$ we will find that this contribution to the diffusion scales as $\mathrm{D}\sim\log(N)$. However, this predicted scaling will only continue as long as the appropriate long-distance cutoff $r_{large}$ is set by the box size. As we consider increasingly large tissues, these mean fields will interfere with that from other events rather than periodic images of the same event. Therefore we expect that the cutoff length (and by proxy, the predicted $\mathrm{D}^{mean}_{(m,a)}$) will become independent of the tissue size in large enough tissues.

C.2 Estimating fluctuating contributions to displacement

Similar to the previous section, we now consider the contribution to diffusion from the fluctuating (third and fourth) terms in Eq. 13. Similar to Eq. 37, the contribution to the tracer dynamics generated by these fluctuating displacements is written as

[TABLE]

where $\delta\mathbf{u}_{i}$ is the fluctuating displacement produced by the $i$-th mitosis(apoptosis) event. This displacement is assumed to have a random direction and magnitude determined by the distribution of $\delta u_{(m,a)}(r)$ from Eq. 12. The mean square displacement again relies on the quantity $\left\langle d^{2}\right\rangle^{fluc}_{(m,a)}$, which again takes the form of a spatial average

[TABLE]

Turning this integral into a sum over bins, the fluctuating displacement data from Section IV can be used to estimate these contributions to diffusion

[TABLE]

Similar to the case of the mean displacements, we may estimate the form of the diffusion Eq. 44 based on some more tangible parameters. This is accomplished by approximating the fluctuating displacements as a uniform region of fluidization with $\delta u_{(m,a)}(r)=w_{(m,a)}\Theta(R_{w}-r)$. The diffusion is then

[TABLE]

Fig. 5 (b) indicates that all of these parameters above are independent of system size as long as $L$ is sufficiently larger than $R_{w}$. Therefore the contribution of these displacements to cell dynamics is also independent of system size as long as $\sqrt{N_{cells}/\rho}>2R_{w}$.

These fluctuating components may also be checked for scaling with increasing system size. As shown in Fig. 12-b, the mitosis event creates a finite region of fluidization which is consistent across system sizes.