Dense active matter model of motion patterns in confluent cell monolayers

TL;DR

This paper demonstrates that dense active matter physics can explain the complex motion patterns in confluent cell monolayers, combining analytical models, simulations, and experiments to match observed velocity correlations.

Contribution

It introduces a unified active matter framework to describe cell sheet dynamics, integrating continuum elasticity, normal modes, and agent-based models with experimental validation.

Findings

Analytical model accurately predicts velocity correlations.

Simulations match experimental data over long timescales.

Persistent uncoordinated motility drives swirl-like correlations.

Abstract

Epithelial cell monolayers show remarkable displacement and velocity correlations over distances of ten or more cell sizes that are reminiscent of supercooled liquids and active nematics. We show that many observed features can be described within the framework of dense active matter, and argue that persistent uncoordinated cell motility coupled to the collective elastic modes of the cell sheet is sufficient to produce swirl-like correlations. We obtain this result using both continuum active linear elasticity and a normal modes formalism, and validate analytical predictions with numerical simulations of two agent-based cell models, soft elastic particles and the self-propelled Voronoi model together with in-vitro experiments of confluent corneal epithelial cell sheets. Simulations and normal mode analysis perfectly match when tissue-level reorganisation occurs on times longer than the persistence time of cell motility. Our analytical model quantitatively matches measured velocity correlation functions over more than a decade with a single fitting parameter.

Figures

Equations

$$\zeta {\dot{\mathbf{r}}}_i={\mathbf{F}}_i^{\text{act}}+{\mathbf{F}}_i^{\text{int}},$$

$${\dot{\theta }}_i={\eta }_i,⟨{\eta }_i\left(t\right){\eta }_j(t^{\prime })⟩=\frac{1}{\tau }{\delta }_{ij}\delta \left(t-t^{\prime }\right),$$

$$\zeta \delta {\dot{\mathbf{r}}}_i=\zeta v_0{\hat{\mathbf{n}}}_i-\sum _j{\mathbf{K}}_{ij}\cdot \delta {\mathbf{r}}_j,$$

$$\zeta {\dot{a}}_{\nu }=-{\lambda }_{\nu }{a}_{\nu }+{\eta }_{\nu },$$

$$E_{\nu }=\frac{1}{2}{\lambda }_{\nu }⟨a_{\nu }^2⟩=\frac{\zeta v_0^2\tau }{4\left(1+{\lambda }_{\nu }\tau /\zeta \right)},$$

$$⟨|\mathbf{v}(\mathbf{q}){|}^2⟩=\sum _{\nu }\frac{v_0^2}{2(1+{\lambda }_{\nu }\tau /\zeta )}|{\mathbf{\xi }}_{\nu }(\mathbf{q}){|}^2,$$

$$\zeta \dot{\mathbf{u}}(\mathbf{R})=\zeta v_0\hat{\mathbf{n}}(\mathbf{R})-\sum _{{\mathbf{R}}^{\prime }}\mathbf{D}(\mathbf{R}-{\mathbf{R}}^{\prime })\mathbf{u}({\mathbf{R}}^{\prime }),$$

$$-\mathrm{i}\zeta \omega \mathbf{u}(\mathbf{q},\omega )={\mathbf{F}}^{\mathrm{act}}(\mathbf{q},\omega )-\mathbf{D}(\mathbf{q})\mathbf{u}(\mathbf{q},\omega ),$$

$$⟨{\mathbf{F}}^{\mathrm{act}}(\mathbf{q},\omega )\cdot {\mathbf{F}}^{\mathrm{act}}(-\mathbf{q},{\omega }^{\prime })⟩=2\pi N{\zeta }^2{v}_0^2\frac{2\tau }{1+{\left(\tau \omega \right)}^2}\delta (\omega +{\omega }^{\prime }).$$

$$⟨{\left|\mathbf{v}(\mathbf{q})\right|}^2⟩=\frac{N{v}_0^2}{2}\left[\frac{1}{1+({\xi }_{\mathrm{L}}q{)}^2}+\frac{1}{1+({\xi }_{\mathrm{T}}q{)}^2}\right],$$

$$⟨{\left|\mathbf{v}\right|}^2⟩=\frac{v_0^2{a}^2}{8\pi }\left[\frac{1}{{\xi }_{\mathrm{L}}^2}\log \left(1+{\xi }_{\mathrm{L}}^2{q}_{\text{m}}^2\right)+\frac{1}{{\xi }_{\mathrm{T}}^2}\log \left(1+{\xi }_{\mathrm{T}}^2{q}_{\text{m}}^2\right)\right],$$

$$C_{vv}(\mathbf{r})=\frac{1}{L^2}\int d^2{\mathbf{r}}_0⟨\mathbf{v}({\mathbf{r}}_0+\mathbf{r})\cdot \mathbf{v}({\mathbf{r}}_0)⟩.$$

$$C_{vv}(\mathbf{r})=\frac{a^2}{(2\pi {)}^{2}N}\int d^2\mathbf{q}⟨|\mathbf{v}(\mathbf{q}){|}^2⟩{\mathrm{e}}^{-\mathrm{i}\mathbf{q}\cdot \mathbf{r}}.$$

$$C_{vv}(\mathbf{r})=\frac{a^2{v}_0^2}{4\pi }\left[\frac{K_0(r/{\xi }_{\mathrm{L}})}{{\xi }_{\mathrm{L}}^2}+\frac{K_0(r/{\xi }_{\mathrm{T}})}{{\xi }_{\mathrm{T}}^2}\right],$$

$$C_{vv}(\mathbf{r})\approx \frac{a^2{v}_0^2}{4\pi }\sqrt{\frac{\pi }{2r}}\left(\frac{{\mathrm{e}}^{-r/{\xi }_{\mathrm{L}}}}{{\xi }_{\mathrm{L}}^{3/2}}+\frac{{\mathrm{e}}^{-r/{\xi }_{\mathrm{T}}}}{{\xi }_{\mathrm{T}}^{3/2}}\right),$$

$$⟨{\left|\bar{\mathbf{v}}(\mathbf{q})\right|}^2⟩=\frac{N}{2}{\left(\frac{v_0}{\bar{v}}\right)}^2\left[\frac{1}{1+{\xi }_{\mathrm{L}}^2{q}^2}+\frac{1}{1+{\xi }_{\mathrm{T}}^2{q}^2}\right],$$

$$\zeta \delta {\dot{\mathbf{r}}}_i=\zeta v_0{\hat{\mathbf{n}}}_i-\sum _j{\mathbf{K}}_{ij}\cdot \delta {\mathbf{r}}_j,$$

$${\dot{\varphi }}_i={\eta }_i\left(t\right),⟨{\eta }_i(t)⟩=0,⟨{\eta }_i\left(t\right){\eta }_j\left(t^{\prime }\right)⟩=\frac{2}{\tau }{\delta }_{ij}\delta \left(t-t^{\prime }\right),$$

$$\delta {\mathbf{r}}_i=\sum _{\nu =1}^{2N}{a}_{\nu }{\mathbf{\xi }}_i^{\nu }.$$

$$\zeta {\dot{a}}_{\nu }=-{\lambda }_{\nu }{a}_{\nu }+{\eta }_{\nu },\text{where}{\eta }_{\nu }=v_0\zeta \sum _{i=1}^{2N}{\hat{\mathbf{n}}}_i\cdot {\mathbf{\xi }}_i^{\nu },$$

$$⟨{\eta }_{\nu }\left(t\right){\eta }_{{\nu }^{\prime }}\left(t^{\prime }\right)⟩=C\left(t-t^{\prime }\right){\delta }_{\nu ,{\nu }^{\prime }}\text{with}C\left(t-t^{\prime }\right)=\frac{{\zeta }^2{v}_0^2}{2}⟨\cos \left[\varphi \left(t\right)-\varphi \left(t^{\prime }\right)\right]⟩,$$

$$p\left(\Delta \varphi ,t\right)=\frac{1}{\sqrt{4\pi |t|/\tau }}{\mathrm{e}}^{-(\Delta \varphi {)}^2\frac{\tau }{4|t|}}.$$

$$C\left(t\right)=\frac{{\zeta }^2{v}_0^2}{2}{\mathrm{e}}^{-|t|/\tau },$$

$$a_{\nu }\left(t\right)=a_{\nu }\left(0\right){\mathrm{e}}^{-\frac{{\lambda }_{\nu }}{\zeta }t}+{\int }_0^{t}d{t}^{\prime }\frac{{\eta }_{\nu }\left(t^{\prime }\right)}{\zeta }{\mathrm{e}}^{-\frac{{\lambda }_{\nu }}{\zeta }\left(t-t^{\prime }\right)}.$$

$$⟨a_{\nu }^2⟩=\frac{\zeta }{{\lambda }_{\nu }}{\int }_0^{\infty }dv\frac{1}{{\zeta }^2}C\left(v\right){\mathrm{e}}^{-\frac{{\lambda }_{\nu }}{\zeta }v}.$$

$$⟨a_{\nu }^2⟩=\frac{\zeta }{{\lambda }_{\nu }}{\int }_0^{\infty }dv\frac{v_0^2}{2}{\mathrm{e}}^{-v/\tau }{\mathrm{e}}^{-\frac{{\lambda }_{\nu }}{\zeta }v}=\frac{\zeta v_0^2}{2{\lambda }_{\nu }}{\int }_0^{\infty }dv{\mathrm{e}}^{-\left(\frac{1}{\tau }+\frac{{\lambda }_{\nu }}{\zeta }\right)v}=\frac{\zeta v_0^2\tau }{2{\lambda }_{\nu }\left(1+\frac{{\lambda }_{\nu }}{\zeta }\tau \right)}$$

$$E_{\nu }=⟨\frac{1}{2}{\lambda }_{\nu }{a}_{\nu }^2⟩=\frac{\zeta v_0^2\tau }{4\left(1+\frac{{\lambda }_{\nu }}{\zeta }\tau \right)}.$$

$$\hat{G}\left(\mathbf{q}\right)=⟨\mathbf{v}\left(\mathbf{q}\right)\cdot {\mathbf{v}}^{\ast }\left(\mathbf{q}\right)⟩\text{with}\mathbf{v}\left(\mathbf{q}\right)=\frac{1}{N}\sum _{j=1}^N{\mathrm{e}}^{\mathrm{i}\mathbf{q}\cdot {\mathbf{r}}_j^0}\delta {\dot{\mathbf{r}}}_j,$$

$$\hat{G}\left(\mathbf{q}\right)=\sum _{\nu ,{\nu }^{\prime }}⟨{\dot{a}}_{\nu }{\dot{a}}_{{\nu }^{\prime }}⟩{\mathbf{\xi }}_{\nu }\left(\mathbf{q}\right)\cdot {\mathbf{\xi }}_{{\nu }^{\prime }}^{\ast }\left(\mathbf{q}\right),\text{with}{\mathbf{\xi }}_{\nu }\left(\mathbf{q}\right)=\frac{1}{N}\sum _{j=1}^N{\mathrm{e}}^{\mathrm{i}\mathbf{q}\cdot {\mathbf{r}}_j^0}{\mathbf{\xi }}_j^{\nu },$$

$$⟨{\dot{a}}_{\nu }{\dot{a}}_{{\nu }^{\prime }}⟩=\frac{1}{{\zeta }^2}\left[{\lambda }_{\nu }{\lambda }_{{\nu }^{\prime }}⟨a_{\nu }{a}_{{\nu }^{\prime }}⟩-{\lambda }_{\nu }⟨a_{\nu }{\eta }_{{\nu }^{\prime }}⟩-{\lambda }_{{\nu }^{\prime }}⟨a_{{\nu }^{\prime }}{\eta }_{\nu }⟩+⟨{\eta }_{\nu }{\eta }_{{\nu }^{\prime }}⟩\right]=\frac{1}{{\zeta }^2}\left[{\lambda }_{\nu }^2⟨a_{\nu }^2⟩-2{\lambda }_{\nu }⟨a_{\nu }{\eta }_{\nu }⟩+⟨{\eta }_{\nu }^2⟩\right]{\delta }_{\nu ,{\nu }^{\prime }},$$

Full text

I Introduction

Collective cell migration is of fundamental importance in embryonic development Weijer (2009); Friedl and Gilmour (2009); Scarpa and Mayor (2016); Saw et al. (2017), organ regeneration and wound healing Safferling et al. (2013). During embryogenesis, robust regulation of collective cell migration is key for formation of complex tissues and organs. In adult tissues, a paradigmatic model of collective cell migration is the radial migration of corneal epithelial cells across the surface of the eye Collinson et al. (2002); Di Girolamo et al. (2015). Major advances in our understanding of collective cell migration have been obtained from in-vitro experiments on epithelial cell monolayers Poujade et al. (2007); Trepat et al. (2009); Tambe et al. (2011); Chepizhko et al. (2016); Saw et al. (2017). A key observation is that collective cell migration is an emergent, strongly correlated phenomenon that cannot be understood by studying the migration of individual cells Trepat and Sahai (2018). For example, forces in a monolayer are transmitted over long distances via a global tug-of-war mechanism Trepat et al. (2009). The landscape of mechanical stresses is rugged with local stresses that are correlated over distances spanning multiple cell sizes Tambe et al. (2011). These strong correlations lead to the tendency of individual cells to migrate along the local orientation of the maximal principal stress (plithotaxis Tambe et al. (2011); Trepat and Fredberg (2011)) and a tendency of a collection of migrating epithelial cells to move towards empty regions of space (kenotaxis Kim et al. (2013)). Furthermore, such coordination mechanisms lead to propagating waves in confined clusters Notbohm et al. (2016); Deforet et al. (2014), expanding colonies Serra-Picamal et al. (2012) and in colliding monolayers Rodríguez-Franco et al. (2017), which all occur in the absence of inertia.

Active matter physics Vicsek and Zafeiris (2012); Marchetti et al. (2013) offers a natural framework for describing subcellular, cellular and tissue-level processes. It studies the collective motion patterns of agents each internally able to convert energy into directed motion. In the dense limit, motility leads to a number of unexpected motion patterns, including flocking Kumar et al. (2014), oscillations Henkes et al. (2011), active liquid crystalline Marchetti et al. (2013), and arrested, glassy phases Berthier and Kurchan (2013). In-silico studies Grossman et al. (2008); Szabó et al. (2010); Henkes et al. (2011); Kabla (2012); Chepizhko et al. (2018) in the dense regime have been instrumental in describing and classifying experimentally observed collective active motion.

Continuum active gel theories Kruse et al. (2005); Prost et al. (2015) are able to capture many aspects of cell mechanics Banerjee et al. (2019), including spontaneous flow of cortical actin Joanny and Prost (2009) and contractile cell traction profiles with the substrate Kruse et al. (2006). In some cases, cell shapes form a nematic-like texture Saw et al. (2017); Kawaguchi et al. (2017) and topological defects present in such texture have been argued to assist in the extrusion of apoptotic cells Saw et al. (2017). To date, however, the cell-level origin of the heterogeneity in flow patterns and stress profiles in cell sheets is still poorly understood. Many epithelial tissues show little or no local nematic order or polarization, and even where order is present, the local flow and stress patterns only follow the continuum prediction on average, while individual patterns are dominated by fluctuations. This suggests that active nematic and active gel approaches capture only part of the picture.

Confluent cell monolayers exhibit similar dynamical behaviour to supercooled liquids approaching a glass transition. One observes spatio-temporally correlated heterogeneous patterns in cell displacements Angelini et al. (2011) known as dynamic heterogeneities Berthier (2011), a hallmark of the glass transition Berthier et al. (2011) between a slow, albeit flowing liquid phase and an arrested amorphous glassy state. The notion that collectives of cells reside in the vicinity of a liquid to solid transition provides profound biological insight into the mechanisms of collective cell migration. By tuning the motility and internal properties of individual cells, e.g. cell shape Bi et al. (2016); Barton et al. (2017); Merkel and Manning (2018) or cell-cell adhesion Garcia et al. (2015), a living system can drive itself across this transition and rather accurately control cell motion within the sheet. This establishes a picture in which tissue level patterning is not solely determined by biochemistry (e.g., the distribution of morphogens) but is also driven by mechanical cues.

In this paper, we show that the cell-level heterogeneity, that is variations in size, shape, mechanical properties or motility between individual cells of the same type inherent to any cell monolayer, together with individual, persistent, cell motility and soft elastic repulsion between neighbouring cells leads to correlation patterns in the cell motion, with correlation lengths exceeding ten or more cell sizes. Inspired by the theory of sheared granular materials Maloney and Lemaître (2004, 2006), we develop a normal modes formalism for the linear response of confluent cell sheets to active perturbations (see Fig.  1A), and derive a displacement correlation function with a characteristic length scale of flow patterns. Using numerical simulations of models for cell sheets, including a soft disk model as well as a self-propelled Voronoi model (SPV) Bi et al. (2016); Barton et al. (2017), we show that our analytical model provides an excellent match for both types of simulations up to a point where substantial flow in the sheet begins to subtly alter the correlation functions (Fig.  1B). At the level of linear elasticity, we are able to make an analytical prediction for the velocity correlation function and the mean velocity in a generic cell sheet. We test our theoretical predictions, which apply to any confluent epithelial cell sheet on a solid substrate dominated by uncoordinated migration, with time-lapse observations of corneal epithelial cells grown to confluence on a tissue culture plastic substrate. We find very good agreement between experimental velocity correlations and analytical predictions and are, thus, able to construct fully parametrized soft disk and SPV model simulations of the system (Fig.  1C) that quantitatively match the experiment. Garcia, et al. Garcia et al. (2015) observed similar correlations and proposed a scaling theory based on coherently moving cell clusters, and either cell-substrate or cell-cell dissipation. Our approach generalises their result for cell-substrate dissipation, and we recover both the scaling results and also find quantitative agreement with the experiments presented in ref.   Garcia et al. (2015).

II Results

II.1 Model overview

We model the monolayer as a dense packing of soft, self-propelled agents that move with overdamped dynamics, and where the main source of dissipation is cell-substrate friction. The equations of motion for cell centers are

[TABLE]

where $\zeta$ is the cell-substrate friction coefficient, $\mathbf{F}_{i}^{\text{act}}$ is the net motile force resulting from the cell-substrate stress transfer, and $\mathbf{F}_{i}^{\text{int}}$ is the interaction force between cell $i$ and its neighbours. Commonly used interaction models are short-ranged pair forces with attractive and repulsive components Szabó et al. (2006) and SPV models Bi et al. (2016); Barton et al. (2017). Here we only require that the inter-cell forces can be written as the gradients of a potential energy that depends on the positions of cell centres, $\mathbf{F}_{i}^{\text{int}}=-\nabla_{\mathbf{r}_{i}}V\left(\{\mathbf{r}_{j}\}\right)$. Furthermore, we neglect cell division and extrusion for now, but we will reconsider the issue when we match simulations to experiment below. The precise form and molecular origin of the active propulsion force $\mathbf{F}_{i}^{\text{act}}$ is a topic of ongoing debate, and interactions between cells through flocking, nematic alignment, plithotaxis and kenotaxis have all been proposed. What is clear, however, is that all alignment mechanisms occur over a substantial background of uncoordinated motility, and therefore, as a base model, we assume that the active cell forces undergo random, uncorrelated fluctuations in direction. With $\mathbf{F}_{i}^{\text{act}}=F^{\text{act}}\hat{\mathbf{n}}_{i}$, where $\hat{\mathbf{n}}_{i}$ is the unit vector that makes an angle $\theta_{i}$ with the $x-$axis of the laboratory frame, the angular dynamics is

[TABLE]

where $\tau$ sets a persistence time scale, and different cells are not coupled (Fig.  1A). This dynamics is equivalent to active Brownian particles Marchetti et al. (2016), and in isolation, model cell motion is a persistent random walk. At sufficiently low driving, such models form active glasses Berthier and Kurchan (2013); Mandal et al. (2016); Nandi et al. (2018); Mandal et al. (2019), where the system moves through a series of local energy minima (i.e., spatial configurations of cells) on the time scale of the alpha-relaxation time $\tau_{\alpha}$, which diverges at dynamical arrest.

We now develop a linear response formalism. As shown in Fig.  1A, on time scales below $\tau_{\alpha}$, the self-propulsion reduces to a stochastic, time-correlated force fluctuating inside a local energy minimum. If the persistence time scale $\tau\ll\tau_{\alpha}$, the full dynamics can be described by a statistical average over long periods fluctuating around different energy minima, by assuming that the brief periods during which the system rearranges do not contribute appreciably (see also ref.   Mandal et al. (2019)). We linearize the interaction forces in the vicinity of an energy minimum, i.e., a mechanically stable or jammed configuration $\{\mathbf{r}_{i}^{0}\}$ by introducing $\delta\mathbf{r}_{i}=\mathbf{r}_{i}-\mathbf{r}_{i}^{0}$. After introducing the active velocity $v_{0}=F^{\text{act}}/\zeta$, Eq. (1) becomes

[TABLE]

where $\mathbf{K}_{ij}=\frac{\partial^{2}V\left(\{\mathbf{r}_{i}\}\right)}{\partial\mathbf{r}_{i}\partial\mathbf{r}_{j}}|_{\{\mathbf{r}_{j}^{0}\}}$ is the dynamical matrix Wyart et al. (2005), organised as $2\times 2$ blocks corresponding to cells $i$ and $j$. In this limit, we can solve the dynamics exactly, see Supplementary Note 1.

II.2 Normal mode formulation

Assuming that there are a sufficient number of inter-cell forces to constrain the tissue to be elastic at short time scales, the dynamical matrix has $2N$ independent normal modes ${\bm{\xi}}^{\nu}$ with positive eigenvalues $\lambda_{\nu}$. If we project Eq. (3) onto the normal modes, we obtain

[TABLE]

where $a_{\nu}=\sum_{i}\delta{\bf r}_{i}\cdot{\bm{\xi}}^{\nu}_{i}$ and the self-propulsion force has been projected onto the modes, $\eta_{\nu}=\zeta v_{0}\sum_{i}\hat{\bf n}_{i}\cdot{\bm{\xi}}^{\nu}_{i}$. The self-propulsion then acts like a time-correlated Ornstein-Uhlenbeck noise (see Supplementary Note 1), with $\langle\eta_{\nu}(t)\eta_{\nu^{\prime}}(t^{\prime})\rangle=\frac{1}{2}\zeta^{2}v_{0}^{2}\exp({-|t-t^{\prime}|/\tau})\delta_{\nu,\nu^{\prime}}$. We can integrate Eq. (4) and obtain the moments of $a_{\nu}$. In particular, the mean energy per mode is given by

[TABLE]

explicitly showing that equipartition is broken due to the mode-dependence induced by $\lambda_{\nu}$ in Eq. (5). In the limit $\tau\rightarrow 0$, we recover an effective thermal equilibrium, $E_{\nu}\rightarrow\zeta v_{0}^{2}\tau/4:=T_{\text{eff}}/2$, where $T_{\text{eff}}=\zeta v_{0}^{2}\tau/2$, consistent with previous work Bi et al. (2016); Mandal et al. (2016). In Fig.  2A, the left-most column is for a simulated thermal system, with properties that are nearly indistinguishable from the $\tau=0.2$ results. In the opposite, high persistence limit when $\tau\to\infty$, we obtain instead $E_{\nu}=\zeta^{2}v_{0}^{2}/4\lambda_{\nu}$, i.e., a divergence of the contribution of the lowest modes. A predominance of the lowest modes in active driven systems was also noted in ref.  Henkes et al. (2011); Ferrante et al. (2013); Bi et al. (2016).

It thus becomes clear that for large values of $\tau$, $T_{\text{eff}}$ can no longer be interpreted as temperature since the fluctuation-dissipation theorem is no longer valid. An analogous result to the $\tau\rightarrow\infty$ limit has been obtain in granular material with an externally applied shear Maloney and Lemaître (2004, 2006), showing that the mechanisms at play are generic, and that tuning $\tau$ allows active systems to bridge between features of thermal systems and (self-)sheared systems (see also ref.  Mandal et al. (2019); Bi et al. (2016)). However, the glass transition lies on a curve of constant $T_{\text{eff}}$ with moderate $\tau$ contributions (Fig.  2 and ref.  Nandi et al. (2018)), making $T_{\text{eff}}$ a convenient parameter, in spite of its lack of genuine thermodynamic interpretation.

In order to make connections to experiments on cells sheets, we compute several directly measurable quantities. One measure that is easily extracted from microscopy images is the velocity field, using particle image velocimetry (PIV) Raffel et al. (2007). We compute the Fourier space velocity correlation function, $\langle|\mathbf{v}(\mathbf{q})|^{2}\rangle=\langle{\bf v}({\bf q})\cdot{\bf v}^{*}({\bf q})\rangle$, with $\mathbf{v}(\mathbf{q})=1/N\sum_{j=1}^{N}\mathrm{e}^{\mathrm{i}{\bf q}\cdot\mathbf{r}_{j}^{0}}\,\delta\dot{\mathbf{r}}_{j}$, where the $\{{\bf r}_{j}^{0}\}$ are the positions of the cell centres at mechanical equilibrium. Expanding over the normal modes, and taking into account the statistical independence of the time derivatives $\dot{a}_{\nu}$ of the modes amplitudes for different modes, we first derive (see Supplementary Note 1) the mode correlations $\langle\dot{a}_{\nu}^{2}\rangle=v_{0}^{2}/\left[2\left(1+\lambda_{\nu}\tau/\zeta\right)\right]$, so that in Fourier space, we obtain

[TABLE]

where ${\bm{\xi}}_{\nu}({\bf q})$ is the Fourier transform of the vector ${\bm{\xi}}_{i}^{\nu}$.

II.3 Continuum elastic formulation

In most practical situations, it is impossible to extract either the normal modes or their eigenvalues. While it is possible to do so in, e.g., colloidal particle experiments Chen et al. (2010); Henkes et al. (2012), the current methods are strictly restricted to thermal equilibrium, and also require an extreme amount of data. Fortunately, the results above are easily recast into the language of solid state physics Ashcroft and Mermin (1976). We rewrite Eq. (3) as

[TABLE]

where $\mathbf{u}(\mathbf{R})$ denotes the elastic deformations from the equilibrium positions $\mathbf{R}$ in the solid, and $\mathbf{D}(\mathbf{R}-\mathbf{R}^{\prime})$ is the continuum dynamical matrix. The normal modes of the system are now simply Fourier modes with

[TABLE]

where $\mathbf{F}^{\mathrm{act}}(\mathbf{q},\omega)=\zeta v_{0}\int_{-\infty}^{\infty}\!\!\!dt\sum_{R}\hat{\mathbf{n}}(\mathbf{R},t)\mathrm{e}^{\mathrm{i}\omega t}\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot R}$ and $\mathbf{D}(\mathbf{q})$ are Fourier transforms of the active force and the dynamical matrix, respectively. Note that we assume that the system has a finite volume, so that Fourier modes are discrete. At scales above the cell size $a$, the noise $\hat{\mathbf{n}}(\mathbf{R},t)$ is spatially uncorrelated, and we find the noise correlators (see Supplementary Note 2)

[TABLE]

The dynamical matrix $\mathbf{D}(\mathbf{q})$ has two independent eigenmodes in two dimensions, one longitudinal $\hat{\mathbf{q}}$ with eigenvalue $B+\mu$ and one transverse one $\hat{\mathbf{q}}^{\perp}$ with eigenvalue $\mu$, where $B$ and $\mu$ are the bulk and shear moduli, respectively. We can then decompose our solution into longitudinal and transverse parts, $\mathbf{u}(\mathbf{q},\omega)=u_{\mathrm{L}}(\mathbf{q},\omega)\hat{\mathbf{q}}+u_{\mathrm{T}}(\mathbf{q},\omega)\hat{\mathbf{q}}^{\perp}$. We are interested in the equal time, Fourier transform of the velocity, which we find to be (see Supplementary Note 2)

[TABLE]

where we have introduced the longitudinal and transverse correlation lengths $\xi_{\mathrm{L}}^{2}=\left(B+\mu\right)\tau/\zeta$ and $\xi_{\mathrm{T}}^{2}=\mu\tau/\zeta$. Note that there are subtle differences in prefactors between expressions for velocity correlation function in Fourier space (cf., Eq. (S19), Eq. (10) and Eq. (52) in Supplementary Note 2), and that in two dimensions, $[\mu/\zeta]=[B/\zeta]=L^{2}T^{-1}$. As discussed in detail in Supplementary Note 2, these differences are due the use of discrete vs. continuum Fourier transforms and are important for comparison with simulations and experiments. Finally, the mean square velocity of the particles $\langle\left|\mathbf{v}\right|^{2}\rangle=\langle\frac{1}{N}\sum_{i}|\mathbf{v}_{i}|^{2}\rangle$ decreases with active correlation time as

[TABLE]

where $q_{\text{m}}=2\pi/a$ is the maximum wavenumber and the high-$q$ cutoff $a$ is of the order of the cell size. Eq. (10) shows that the correlation length of the system scales as $\sqrt{\tau}$. In the limit $\tau\rightarrow\infty$, $\langle|\mathbf{v}(\mathbf{q})|^{2}\rangle$ diverges at low $q$, as was found in ref.  Szamel (2016). The dominant scaling $\langle\left|\mathbf{v}\right|^{2}\rangle\sim 1/\xi^{2}$ is the same as results from the scaling Ansatz for cell-substrate dominated coordinated motion obtained by ref.  Garcia et al. (2015).

While Eq. (10) is elegant, correlations of cell velocities expressed in the Fourier space are not easy to interpret. Therefore, we derive a more intuitive, real space expression for the correlation of velocities of cells separated by $r$, defined as

[TABLE]

In the infinite size limit $L\to\infty$, the real space correlation function $C_{vv}(\mathbf{r})$ can be evaluated from the Fourier correlation $\langle|\mathbf{v}(\mathbf{q})|^{2}\rangle$ as

[TABLE]

Using Eq. (11), one finds the explicit result (see Supplementary Note 2)

[TABLE]

where $K_{0}$ is the modified Bessel function of the second kind. Note that this expression describes velocity correlations for $r>a$, with $a$ the cell size. For $r/\xi_{\mathrm{L},\mathrm{T}}\gg 1$, i.e. for distance much larger than the correlation lengths,

[TABLE]

i.e., as expected and consistent with the results of ref.   Garcia et al. (2015), $C_{vv}$ decays exponentially at large distances.

II.4 Comparison to simulations

We proceed to compare predictions made in the previous section to the correlation function measured in numerical simulations of an active Brownian soft disk model, as well as to an SPV model. The active Brownian model is defined by Eq. (1) and Eq. (2), with self-propulsion force ${\bf F}_{i}^{\rm act}=v_{0}\hat{\mathbf{n}}_{i}$ and pair interaction forces ${\bf F}_{ij}$ that are purely repulsive. We simulate a confluent sheet in this model by setting the packing fraction to $\phi=1$ in periodic boundary conditions. The SPV model is the same as introduced in refs.  Bi et al. (2016); Barton et al. (2017), and assumes that every cell is defined by the Voronoi tile corresponding to its centre. For this model, we choose the dimensionless shape factor $\bar{p}_{0}=3.6$, putting the passive system into the solid part of the phase diagram Bi et al. (2016); Sussman and Merkel (2018), and we employ open boundary conditions. Please see the method section for full details of the numerical models and simulation protocols.

The effective temperature $T_{\text{eff}}=\zeta v_{0}^{2}\tau/2$ has emerged as a good predictor of the active glass transition Szamel (2016); Mandal et al. (2016), at least at low $\tau$, and we use it together with $\tau$ itself as the axes of our phase diagram. The liquid or glassy behaviour of the model can be characterised by the alpha relaxation time $\tau_{\alpha}$. Fig.  2 provides a coarse-grained phase diagram where $\tau_{\alpha}$ is represented in gray scale as a function of persistence time $\tau$ and $T_{\text{eff}}$. For a fixed persistence time, the system is liquid at high enough temperature and glassy at low temperature, as expected. Now fixing the effective temperature, the system becomes more glassy when $\tau$ increases. This non-trivial result is consistent with the recent RFOT theory of the active glass transition Nandi et al. (2018) and related simulation results Nandi et al. (2018); Mandal et al. (2019). It can be partly understood from the fact that $v_{0}$ decreases when $\tau$ increases at fixed $T_{\text{eff}}$, meaning that the active force decreases and it becomes more difficult to cross energy barriers. In contrast, existing mode coupling theories of the active glass transition Szamel (2016) only apply in the small $\tau$ regime. We note that the features of the active glass transition of the soft disk model and the SPV model are very similar.

As is apparent from Fig.  1B, the growing correlation length with increasing $\tau$ is readily apparent as swirl-like motion (see also Supplementary Movies 1-4). Fig.  3A-C shows the Fourier velocity correlation $\langle|\mathbf{v}(\mathbf{q})|^{2}\rangle$ measured in the numerical simulations for different values of $v_{0}$, after normalizing $\langle|\mathbf{v}(\mathbf{q})|^{2}\rangle$ by $v_{0}^{2}N$. In panel A, we show that for soft disks at low $T_{\text{eff}}=0.005$, where the system is solid, the correlation function develops a dramatic $1/q^{2}$ slope as $\tau$ increases (dots), exactly in line with our modes predictions (lines). We can determine the bulk and shear moduli of the soft disk system ($B=1.684\pm 0.008$, $\mu=0.510\pm 0.004$, see Methods section) and then draw the predictions of Eq. (10) on the same plot (dashed lines). At low $q$, where the continuum elastic approximation is valid, we have excellent agreement, and at larger $q$, the peak associated with the static structure factor becomes apparent (in the limit $\tau\rightarrow 0$, the correlation function reduces to $S(q)$, see Supplementary Figure 2). In panel C, we show the same simulation results for the SPV model (dots), accompanied by the continuum predictions (dashed lines) using $B=7.0$ and $\mu=0.5$, as estimated from ref.   Sussman and Merkel (2018) for $\bar{p}_{0}=3.6$. We did not compute normal modes for the SPV model. Note that due to $\mu\ll\mu+B$, the contribution of the transverse correlations dominate the analytical results in both cases. In panel $B$, we show the soft disk simulation at $\tau=20$ when the transition to a liquid is crossed as a function of $T_{\text{eff}}$. Deviations from the normal mode predictions become apparent only at the two largest values of $T_{\text{eff}}$, when $\tau>\tau_{\alpha}$ (Fig.  2A), and even in these very liquid systems, a significant activity-induced correlation length persists. In Supplementary Figure 3, we show that for all $\tau$ and both soft disk and SPV models, our predictions remain in excellent agreement with the simulations for $T_{\text{eff}}=0.02$, where $\tau\lesssim\tau_{\alpha}$. In Fig.  3D, we show the mean square velocity normalized by $v_{0}$ as a function of the dimensionless transverse correlation length $\xi_{T}/a\sim\sqrt{\tau}$, for all our simulations, using $a=\sigma$, the particle radius. The dramatic drop corresponds to elastic energy being stored in distortions of the sheet, and it is in very good agreement with our analytical prediction in Eq. (11) (solid lines). In Supplementary Figure 4, we compare our numerical results for the spatial velocity correlations to the analytical prediction Eq. (14). The data and the predictions are in reasonable agreement.

II.5 Comparison to experiment

We now compare our theoretical predictions and numerical simulations with experimental data obtained from immortalised human corneal epithelial cells grown on a tissue culture plastic substrate (see Methods section and Supplementary Movie 5). We use PIV to extract the velocity fields corresponding to collective cell migration (Fig.  4A, Fig. 1C and Supplementary Movie 6). We first extract a mean velocity of $\bar{v}=\sqrt{\langle|\mathbf{v}(\mathbf{r},t)|^{2}\rangle}=12\pm 2$ µm  h*-1* ($n=5$ experiments, see Fig.  4E), consistent with the typical mean velocities of confluent epithelial cell lines grown on hard substrates. To reduce the effects of varying mean cell speed at different times and in different experiments, we use $\bar{\mathbf{v}}(\mathbf{r},t)=\mathbf{v}(\mathbf{r},t)/\sqrt{\langle|\mathbf{v}(\mathbf{r},t)|^{2}\rangle_{\mathbf{r}}}$, i.e., the velocity normalized by its mean-square spatial average at that moment in time. Using direct counting, we find an area per cell of $\langle A\rangle\approx 380$ µm2 corresponding to a particle radius of $\langle\sigma\rangle\approx 11$ µm, setting the microscopic length scale $a$. To compare the experimental result to our theoretical predictions, we perform a Fourier transform on the PIV velocity field and compute $\langle\left|\bar{\mathbf{v}}(\mathbf{q})\right|^{2}\rangle$, shown in Fig.  4B. Using the results from Eq.  (11) we can rewrite Eq.  (10) as

[TABLE]

where $\xi_{\mathrm{L}}^{2}$ and $\xi_{\mathrm{T}}^{2}$ are the longitudinal and transverse squared correlation lengths (with units of µm2) defined below Eq. (10). As can be seen from Eq. (11), the ratio $\frac{v_{0}}{\bar{v}}$ is a function only of the dimensionless ratios $\xi_{\mathrm{L}}/a$ and $\xi_{\mathrm{T}}/a$. If we further make the plausible assumption that the ratio of elastic moduli is the same as in the soft disk simulations, $(\mu+B)/\mu=4.3$, the correlation lengths $\xi_{\mathrm{L}}$ and $\xi_{\mathrm{T}}$ are not independent. Therefore, we are left with a single fitting parameter, $\xi_{\mathrm{T}}$. The best fit to the theory is obtained with $\xi_{\mathrm{T}}=100$ µm as indicated by the solid black line in Fig.  4A, with the interval of confidence denoted by dashed lines. The $q=0$ intercept of the correlation function gives a ratio $v_{0}/\bar{v}\approx 10$, corresponding to the high activity limit where most self propulsion is absorbed by the elastic deformation of the cells. Consistent with this, on the dimensionless plot Fig.  3d we are located at the point $\xi_{T}/a\approx 5,\>\langle v\rangle/v_{0}\approx 0.1$, on the right, strongly active side. The deviations between theory and experiment in the tail of the distribution are not due to loss of high-$q$ information in imaging, as far as we could determine, and the disappearance of the peak at high $q$ is particularly striking in this context.

In Fig.  4C, we show the real-space velocity correlations for the experiments, and the analytical prediction Eq. (14) with $\xi_{T}=100$ µm. We obtain a very good fit for experiments $1$, $2$, $3$ and $7$, but experiments $5$ and $6$ have significantly longer-ranged correlations. This indicates that the precise value of the correlation length is very sensitive to the exact experimental conditions that are not simple to accurately control. The qualitative features of the correlation function are, however, robust. Note that experiments $5$ and $6$ have the same mean density as experiments $1-3,7$.

To match experiments and simulations, we consider the temporal autocorrelation function $\langle\mathbf{v}(t)\cdot\mathbf{v}(0)\rangle$ in Fig.  4D. As Eq. (57) in Supplementary Note 2 shows, it is a complex function with a characteristic inverse S-shape that also depends on the moduli and $q_{\text{m}}$. Using the value of $\xi_{\mathrm{T}}^{2}=10^{4}$ µm2 extracted from fitting $\langle\left|\bar{\mathbf{v}}(\mathbf{q})\right|^{2}\rangle$ and the ratio $(\mu+B)/\mu=4.3$, we used different $\mu/\zeta$ and $\tau$ compatible with $\xi_{T}^{2}=\mu\sigma^{2}/\zeta\tau$ to obtain the best (numerically integrated) analytical fit to the experimental result. We settled on a best fit autocorrelation time of $\tau=2.5$ h and $\mu/(\sigma^{2}\zeta)=60.5~{}$h*-1* (black line in Fig.  4D). Experiments $5$ and $6$ have significantly longer autocorrelation times, and we achieve a good fit to experiment $5$ for $\tau_{5}=20$ h and the same $\mu/\zeta$ (light gray line). This is consistent with the longer spatial correlations observed in Fig.  4D, where the light grey line corresponds to $\xi_{\mathrm{T}}=100$ µm $\sqrt{\tau_{5}/\tau}\approx 283$ µm. There is also potentially weak local cell alignment in the experiment, not considered in the present theory.

We can now fully parametrise particle and SPV model simulations to the experiment as follows. Our results for $\bar{v}$ and the ratio $v_{0}/\bar{v}$ can be combined to give an initial estimate of $v_{0}=120$ µm h*-1*. Then, the normalised time autocorrelation function of the cell velocities is only a function of $\xi$ and $\tau$, and we can use it to determine the elastic moduli. Then, finally, we can determine the appropriate model parameters: In Fig.  4B, the red and blue dashed lines show Eq. (10) with $B$ and $\mu$ chosen with the same ratio as in the previous particle (respectively, SPV) simulations. From these values, we can approximate the parameter values $k/\zeta=\mu/\sigma^{2}$ for the particle model and $K/\zeta=\mu/\langle A\rangle^{2},\Gamma/\zeta=\mu/\langle A\rangle$ for the SPV model. The solid red and blue curves in Fig.  4B show the best fit simulations that we obtain this way, for $k/\zeta=\Gamma/\zeta=55$ h*-1*, $K/\zeta=0.454$ µm*-2* h*-1* and $v_{0}=90$ µm h*-1*, and snapshots are shown in Fig.  1C (see also Supplementary Movies 7, 9 and 10). The red and blue dashed lines in Fig.  4D show the autocorrelations of our matched simulations for soft disk and vertex simulations, respectively.

Our results are in quantitative agreement with ref.   Garcia et al. (2015) for confluent but still motile cells, with a reported maximum correlation length of $\xi=100$ µm, and a cell crawling speed that drops by a factor of $10$ from $\sim 90$ µm h*-1* at low density to this point of maximal correlation. Note also that we have an elastic time scale $(k/\zeta)^{-1}\approx 0.02h$, much shorter than our correlation time scale $\tau=2.5h$, confirming again that we are in the strongly active regime.

As can be seen in Supplementary Movie 5, a significant number of divisions take place in the epithelial sheet during the $48$h of the experiment. While it is difficult to adapt our theory to include divisions, we can simulate our particle model with a steady-state division and extrusion rate at confluence using the model developed in refs.   Matoz-Fernandez et al. (2017a, b). With a typical cell cycle time of $\tau_{\text{div}}=48$h, we obtain results (green line) that are very similar to the model without division (red line), suggesting that typical cell division rates do not change the velocity correlations noticeably (see also Supplementary Movie 8). This result is consistent with the observed separation of motility time scale $\tau$ and division time scale $\tau_{\text{div}}$. We have also considered the effect of weak polar alignment between cells, using the model from ref.  Henkes et al. (2011). For weak alignment with time scales $\tau_{v}\geq\tau$ that do not lead to global flocking of the sheet, which we do not observe, $\langle|\mathbf{v}(\mathbf{q})|^{2}\rangle$ does not change significantly, though we find somewhat longer autocorrelation times. Finally, the simulations can also give us information about the velocity distribution function, a quantity that is not accessible from our theory. In Fig.  4D, we show the experimental normalised velocity distribution (black line with grey confidence interval), together with the distribution we find from the best fit simulations (coloured lines). As can be seen, there is an excellent match in particular with the SPV model simulation. The particle model with division has additional weight in the tail due to the particular division algorithm implemented in the model (overlapping cells pushing away from each other).

III Discussion

In this study, we have developed a general theory of motion in dense epithelial cell sheets (or indeed other dense active assemblies Briand and Dauchot (2016)) that only relies on the interplay between persistent active driving and elastic response. We find an emerging correlation length that depends only on elastic moduli, the substrate friction coefficient and scales with persistence time as $\tau^{1/2}$. While we found an excellent match between theory and simulations, further experimental validations with different cell lines and on larger systems should be performed. Note that without a substrate, the mechanisms of cell activity are very different Rozbicki et al. (2015). More generally, including cell-cell dissipation in addition to cell-substrate dissipation could significantly modify scalings Garcia et al. (2015), a known result in continuum models of dry (substrate dissipation) vs. wet (internal dissipation) active materials. Due to the suppression of tangential slipping between cells, it could also be responsible for the disappearance of the high-$q$ peak in the velocity correlations. The assumption of uncoordinated activity between cells is a strong one, and it will be interesting to extend the theory by including different local mechanisms of alignment Henkes et al. (2011); Ferrante et al. (2013). From a fundamental point of view, our theoretical results (and also the results of ref.  Henkes et al. (2011)) are examples of a larger class of non-equilibrium steady-states that can be treated using a linear response formalism Liverpool (2018).

IV methods

IV.1 Experiment

Spontaneously immortalised, human corneal epithelial cells (HCE-S) (Notara & Daniels, 2010) were plated into a $12$-well plate using growth medium consisting of DMEM/F12 (Gibco) Glutamax, $10$% fetal bovine serum, and $1$% penicillin/streptomycin solution (Gibco). The medium was warmed to $37^{\circ}$C prior to plating and the cells were kept in a humidified incubator at $37^{\circ}$C and an atmosphere of $5$% CO2 overnight until the cells reached confluence. Before imaging the cells were washed with PBS and the medium was replaced with fresh medium buffered with HEPES. The cells were imaged using a phase contrast Leica DM IRB inverted microscope enclosed in a chamber which kept the temperature at $37^{\circ}$C. The automated time-lapse imaging setup took an image at $10$ minute intervals at a magnification of $10$x, corresponding to a field of view of $867$ µm $\times$ $662$ µm that was saved at a resolution of $1300\times 1000$ pixels. The total experimental run time for each culture averaged $48$ hours, or $288$ separate images. The collected data consists of $7$ experimental imaging runs, of which number $3$ and $4$ were consecutive on the same well-plate (number $4$ was not used in this article). Cell extrusions were counted at three time points during the experimental run by direct observation and counting from the still image (extruded cells detach from the surface and round up, appearing as white circles above the cell sheet in phase contrast). From this data, a typical cell number of $N=1400$, and a typical cell radius of $r=10.95$ µm were extracted. Cell movements were determined using Particle Image Velocimetry (PIV), using an iterative plugin for ImageJ PIV . At the finest resolution (level $3$), it provides displacement vectors on a $54\times 40$ grid, corresponding to a resolution of $16$ µm in the $x$ and the $y$-direction, i.e., slightly less than $1$ cell diameter. The numerous extruded cells and the nucleoli inside the nuclei acted as convenient tracer particles for the PIV allowing for accurate measurements.

IV.2 Simulations

The main simulations consist of $N=3183$ particles simulated with either soft repulsion, or the SPV model (in literature also refered to as the Active Vertex Model (AVM)) potential, using SAMoS SAM . The interaction potential for soft harmonic disks is $V_{i}=\sum_{j}\frac{k}{2}(\sigma_{i}+\sigma_{j}-|\mathbf{r}_{j}-\mathbf{r}_{i}|)^{2}$ if $|\mathbf{r}_{j}-\mathbf{r}_{i}|\leq\sigma_{i}+\sigma_{j}$ and [math] otherwise. To emulate a confluent cell sheet, we used periodic boundary conditions at packing fraction $\phi=1$, where $\phi=\sum_{i}\pi\sigma_{i}^{2}/L^{2}$ and thus double-counts overlaps. We also introduce $30\%$ polydispersity in radius, to emulate cell size heterogeneity. At this density, the model at zero activity is deep within the jammed region ($\phi>0.842$) and has a significant range of linear response.

For the SPV, cells are defined as Voronoi polygonal tiles around cell centers, and the multiparticle interaction potential is given by $V_{i}=\frac{K}{2}(A_{i}-A_{0})^{2}+\frac{\Gamma}{2}(P_{i}-P_{0})^{2}$, where $A_{i}$ is the area of the tile, and $P_{i}$ is its perimeter, $K$ and $\Gamma$ are the area and perimeter stiffness coefficients and $A_{0}$ and $P_{0}$ are reference area and perimeter, respectively. SPV is confluent by construction, and its effective rigidity is set by the dimensionless shape parameter $\bar{p}_{0}=P_{0}/\sqrt{A_{0}}$, with a transition from a solid to a fluid that occurs for $\bar{p}_{0}\approx 3.812$. We simulate the model at $\bar{p}_{0}=3.6$, well within the solid region at zero activity Bi et al. (2016); Sussman and Merkel (2018), and also introduce $30\%$ variability in $A_{0}$. AVM was implemented with open boundary conditions, and we use a boundary line tension $\lambda=0.3$ to avoid a fingering instability at the border that appears especially at large $\tau$.

Both models are simulated with overdamped active Brownian dynamics $\zeta\dot{\mathbf{r}}_{i}=v_{0}\hat{\mathbf{n}}_{i}-\nabla_{\mathbf{r}_{i}}V_{i}$, where the orientation vector $\hat{\mathbf{n}}_{i}=(\cos\theta_{i},\sin\theta_{i})$ follows $\dot{\theta}_{i}=\eta_{i}$, $\langle\eta_{i}(t)\eta_{j}(t^{\prime})\rangle=\frac{1}{\tau}\delta_{ij}\delta(t-t^{\prime})$. Equations of motions are integrated using a first order scheme with time step $\delta t=0.01$. Simulations are $5\times 10^{4}$ time units long, with snapshots saved every $50$ time units, and the first $1250$ time units of data are discarded in the data analysis.

IV.3 Velocity correlations and glassy dynamics

We compute the velocity correlation function for a given simulation directly from particle positions and velocities by first computing the Fourier transform. Then for a given $\mathbf{q}$ and configuration, the correlation function is $|\mathbf{v}(\mathbf{q})|^{2}=\mathbf{v}(\mathbf{q})\cdot\mathbf{v}^{*}(\mathbf{q})$, of which we then take a radial $\mathbf{q}$ average, followed by a time average. The procedure is identical for the experimental PIV fields using the grid positions and velocities, with $N=54\times 40$ grid points as normalization. We compute the $\alpha$-relaxation time from the self-intermediate scattering function $S(q,t)=\langle\frac{1}{N}\sum_{j}\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot(\mathbf{r}_{j}(t_{0}+t)-\mathbf{r}_{j}(t_{0})}\rangle_{t_{0},|\mathbf{q}|=q}$, where the angle brackets indicate time and radial averages. At $q=2\pi/\sigma$, we determine $\tau_{\alpha}$ as the first time point where $S(q,t)<0.5$, bounded from above by the simulation time.

IV.4 Normal mode analysis

The normal modes are the eigenvalues and eigenvectors of the Hessian matrix

${\bf K}_{ij}=\partial^{2}V(\{\mathbf{r}_{i}\})/\partial\mathbf{r}_{i}\partial\mathbf{r}_{j}$, evaluated at mechanical equilibrium. We first equilibrate the $t=2500$ snapshot with $v_{0}=0$ for $2\times 10^{5}$ time steps, equivalent to a steepest descent energy minimization. We made sure that results are not sensitive to the choice of snapshot as equilibration starting point (with the exception of the deviations apparent in Fig. fig:fouriervel at $\tau=2000$). For the soft disk model, each individual $ij$ contact with contact normal $\hat{\mathbf{n}}_{ij}$ and tangential $\hat{\mathbf{t}}_{ij}$ vectors contributes a term $K_{ij}=-k\hat{\mathbf{n}}_{ij}\times\hat{\mathbf{n}}_{ij}+|\mathbf{f}_{ij}|\hat{\mathbf{t}}_{ij}\times\hat{\mathbf{t}}_{ij}$ to the $ij$ $2\times 2$ off-diagonal element of the matrix and $-K_{ij}$ is added to the $ii$ diagonal element Wyart et al. (2005). We use the NumPy eigh function (numpy.linalg.eigh). As the system is deep in the jammed phase, with the exception of two translation modes, all eigenvalues of the Hessian are positive. We compute the Fourier spectrum of mode $\nu$ through $\bm{\xi}_{\nu}(\mathbf{q})=\frac{1}{N}\sum_{j}\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{r}_{j}}\bm{\xi}_{\nu}^{j}$ and then $|\bm{\xi}_{\nu}(\mathbf{q})|^{2}=\xi_{\nu}(\mathbf{q})\cdot(\xi_{\nu}(\mathbf{q}))^{*}$. The $2\times 2$ continuum Fourier space dynamical matrix $\mathbf{D}(\mathbf{q})$ has one longitudinal eigenmode along $\hat{\mathbf{q}}$ with eigenvalue $(B+\mu)q^{2}$ and one transverse eigenmode along $\hat{\mathbf{q}}^{\perp}$ with eigenvalue $\mu$. We compute $\mathbf{D}(\mathbf{q})_{\alpha\beta}=\frac{1}{N}\sum_{j}\sum_{l}\mathrm{e}^{\mathrm{i}q_{\alpha}r_{j,\alpha}}H_{jl,\alpha\beta}\mathrm{e}^{-\mathrm{i}q_{\beta}r_{j,\beta}}$, where the greek indices $\alpha,\beta$ correspond to $x$ or $y$ and there is no sum implied. We diagonalise the resulting matrix, and choose the longitudinal eigenvector as the one with the larger projection onto $\hat{\mathbf{q}}$, and from there the longitudinal and transverse eigenvalues $\lambda_{\mathrm{L}}(q)$ and $\lambda_{\mathrm{T}}(q)$ after a radial $\mathbf{q}$ average. We fit $B+\mu$ as the slope of $\lambda_{\mathrm{L}}(q)$ vs $q^{2}$ up to $q=1.5$, and the same for $\mu$ and $\lambda_{\mathrm{T}}(q)$ (Supplementary Figure 3).

Data availability

Data supporting the findings of this manuscript are available from the corresponding authors upon reasonable request.

Code availability

Simulation and analysis code used in this study are available under an open source (GNU GPL v3.0) licence at: https://github.com/sknepneklab/SAMoS (http://doi.org/10.5281/zenodo.3616475).

Author contributions

SH, EB and RS developed the theory. KK performed the experiment and JMC coordinated it. RS developed the SAMoS code used for both particle and SPV model simulations. SH and KK performed the numerical simulations and analysed the numerical and experimental data. SH, EB, RS and JMC wrote the paper.

Acknowledgements.

We acknowledge many helpful discussions with C. Huepe, D. Matoz Fernandez, K. Martens, I. Näthke, R. Sunyer, X. Trepat and C. J. Weijer. SH acknowledges support by the UK BBSRC (grant number BB/N009150/1-2). RS acknowledges support by the UK BBSRC (grant numbers BB/N009789/1-2). JMC was funded by BBSRC Research Grant BB/J015237/1. KK is funded by a BBSRC EASTBIO PhD studentship. The authors declare no competing interests.

Supplementary Note 1

IV.1 Normal mode formulation

We present here a more detailed methodological account of the normal modes formalism used. We consider a model system of $N$ self-propelled soft interacting particles with overdamped dynamics, in the jammed state. In the absence of self-propulsion, the particles have an equilibrium position ${\bf r}_{i}^{0}$, corresponding to a local minimum of the elastic energy. If the interaction potential is linearized around the energy minimum in terms of the displacement $\delta{\bf r}_{i}={\bf r}_{i}-{\bf r}_{i}^{0}$, the dynamics is described by the equation

[TABLE]

where the ${\bf K}_{ij}$’s are the $2\times 2$ blocks of the $2N\times 2N$ dynamical matrix, $v_{0}\hat{{\bf n}}_{i}$ is the self-propulsion term with $\hat{{\bf n}}_{i}=\cos\phi_{i}\mathbf{e}_{x}+\sin\phi_{i}\mathbf{e}_{y}$ (i.e., direction of $\hat{\mathbf{n}}$ is given by the angle $\phi_{i}$ with the $x$ axis of a laboratory reference frame) and $\zeta$ is the friction coefficient. In the absence of inter-particle alignment, the angle $\phi_{i}$ obeys a simple rotational diffusive dynamics with white noise $\eta_{i}\left(t\right)$:

[TABLE]

where we have expressed the inverse rotational diffusion constant as a time scale, $\tau=1/D_{r}$. We note that in general, the system is far out of thermodynamic equilibrium and $D_{r}$ and $\zeta$ are not simply related to each other. In the following, we consider the self-propulsion noise as a (vectorial) colored noise, and characterize its statistics as well as the statistics of the displacements $\delta{\bf r}_{i}$. To this aim, we first expand $\delta{\bf r}_{i}$ over the normal modes, i.e., the eigenvectors of the dynamical matrix. Each normal mode is a $2N$-dimensional vector that can be written as a list of $N$ two-dimensional vectors $({\bm{\xi}}_{1}^{\nu},\dots,{\bm{\xi}}_{N}^{\nu})$, where the index $\nu=1,\dots,2N$ labels the mode; the associated eigenvalue is denoted as $\lambda_{\nu}$. This form of the normal modes is useful as it allows the decomposition of $\delta{\bf r}_{i}$ to be written in the simple form

[TABLE]

Projecting Eq. (S1) on the normal modes, we find the uncoupled set of equations

[TABLE]

is the projection of the self-propulsion force onto the normal mode $\nu$.

IV.2 Self-propulsion force as a persistent noise

We consider the projection $\eta_{\nu}$ of the self-propulsion force on normal mode $\nu$ as a correlated noise, which we now characterize. Since $\eta_{\nu}$ is the sum of many statistically independent contributions with bounded moments, using the Central Limit theorem, we can assume its statistics to be Gaussian. It is also clear, by averaging over the realizations of the stochastic angles $\phi_{i}$, that $\left\langle\eta_{\nu}\left(t\right)\right\rangle=0$. We thus simply need to evaluate the two-time correlation function of $\eta_{\nu}\left(t\right)$. Using the fact that the eigenvectors of the dynamical matrix form an orthonormal basis, we have $\sum_{i=1}^{N}{\bm{\xi}}^{\nu}_{i}\cdot{\bm{\xi}}^{\nu^{\prime}}_{i}=\delta_{\nu,\nu^{\prime}}$. We find

[TABLE]

where $\phi\left(t\right)$ obeys the diffusive dynamics of Eq. (S2). Note that we have used time translation invariance by assuming that the correlation function depends only on the time difference $t-t^{\prime}$. We can thus set $t^{\prime}=0$ without loss of generality. Solving Eq. (S2), the quantity $\Delta\phi=\phi\left(t\right)-\phi\left(0\right)$ is distributed according to

[TABLE]

One then finds, using Eqs. (S5) and (S6),

[TABLE]

i.e., the time correlation of the noise $\eta_{\nu}$ decays exponentially with the correlation (or persistence) time $\tau$. It is worth emphasizing that the statistical properties of the noise $\eta_{\nu}$ are independent of the mode $\nu$.

IV.3 Potential energy spectrum

We now turn to the computation of the average potential energy per mode. Solving Eq. (S4) explicitly for a given realization of the noise $\eta_{\nu}\left(t\right)$, one finds

[TABLE]

From this expression, one can compute the average value $\left\langle a_{\nu}^{2}\left(t\right)\right\rangle$, leading for $t\to\infty$ to

[TABLE]

Using Eq. (S7), we obtain

[TABLE]

or, in terms of average energy per mode

[TABLE]

For very short correlation time $\tau$ (i.e., large diffusion coefficient $D_{r}$), one recovers an effective equipartition of energy over the modes, $E_{\nu}\approx\frac{\zeta v_{0}^{2}\tau}{4}$ even though the system is out-of-equilibrium. For finite correlation time, this result remains valid in the range of modes $\nu$ such that $\tau\ll\zeta\lambda_{\nu}^{-1}$, if such a range exists. However, for large correlation time $\tau$, that is, as soon as there is a wide range of modes such that $\tau\gg\zeta\lambda_{\nu}^{-1}$, equipartition is broken, and the energy spectrum is given by $E_{\nu}\approx\frac{\zeta^{2}v_{0}^{2}}{4\lambda_{\nu}}.$

IV.4 Velocity correlation

Following Maloney06 , we consider the velocity-velocity correlation function $\hat{G}({\bf q})$ in Fourier space, where one can express the (discrete) Fourier transform ${\bf v}\left({\bf q}\right)$ as a function of the particles reference positions ${\bf r}_{i}^{0}$:

[TABLE]

where the star denotes the complex conjugate. Expanding over the normal modes, one finds

[TABLE]

where ${\bm{\xi}}_{\nu}\left({\bf q}\right)$ is the Fourier transform of the vectors $\bm{\xi}_{\nu}$. From Eq. (S4), the quantity $\left\langle\dot{a}_{\nu}\dot{a}_{\nu^{\prime}}\right\rangle$ is expressed as

[TABLE]

where the last equality is due to the modes being uncorrelated. The cross-correlation is in fact not 0, but crucial:

[TABLE]

To sum up, one has according to Eq. (S14)

[TABLE]

with

[TABLE]

Further, using Eqs. (S5), (S7), (S10) and (S15), one obtains

[TABLE]

Combining Eqs. (S13), (S16) and (S18), we derive the final expression for the velocity correlation function:

[TABLE]

Note that we can compute the equal-time, spatial mean square velocity through Parseval’s theorem as

[TABLE]

Supplementary Note 2

IV.1 Continuum elastic formulation

We now turn to the study of the overdamped equations of motion derived from the elastic energy, in the framework of continuum elastic. In two dimensions, the elastic energy of an isotropic elastic solid with bulk modulus $B$ and shear modulus $\mu$ can be written as Chaikin-Lubensky ; Landau

[TABLE]

where $\hat{u}$ is the strain tensor with components $u_{\alpha\beta}=\frac{1}{2}\left[\partial_{\alpha}u_{\beta}+\partial_{\beta}u_{\alpha}\right]$ written as spatial derivatives of the components $\alpha,\beta\in\left\{x,y\right\}$ of the displacement vectors $\mathbf{u}\left(\mathbf{r}\right)=\mathbf{r}^{\prime}\left(\mathbf{r}\right)-\mathbf{r}$ from a reference state $\mathbf{r}$ to the deformed state $\mathbf{r}^{\prime}\left(\mathbf{r}\right)$. The stress tensor $\sigma_{\alpha\beta}=\frac{\delta F_{\rm el}}{\delta u_{\alpha\beta}}$ can then be written as

[TABLE]

where summation over pairs of repeated indices is assumed. Hence, its divergence is given by

[TABLE]

We can then write the overdamped equations of motion for the displacement field

[TABLE]

This last equation can be rewritten in vectorial notation as

[TABLE]

In Fourier space, we can write this relation as

[TABLE]

where $\mathbf{D}\left(\mathbf{q}\right)$ is the Fourier space dynamical matrix, and $q^{2}=q_{x}^{2}+q_{y}^{2}$. The two eigenvalues of the dynamical matrix are

[TABLE]

with normalized eigenvectors

[TABLE]

In other words, for each $\mathbf{q}$, we obtain one longitudinal and one transverse eigenmode, with diffusive equations of motion, where the diffusion coefficients are the two elastic moduli:

[TABLE]

IV.2 Overdamped dynamics with activity

Now including the self-propulsion force, the continuum version of the active equations of motion is given by

[TABLE]

where we have included an active force $\mathbf{F}^{\rm act}({\bf r},t)=\zeta v_{0}\mathbf{\hat{n}}({\bf r},t)$, whose statistical properties will be discussed below. At this stage, we need a brief aside to properly define our conventions for the Fourier transform. This is particularly important because we wish to compare results from numerical simulations and from continuum theory. Numerical simulations are done in a system of relatively large, but finite linear size $L$, and with a minimal length scale given by the particle size $a$, which leads to the use of a discrete space Fourier transform. On the other hand, analytical calculations are made much easier by assuming whenever possible that $L\to\infty$ and $a\to 0$, i.e., using the continuous Fourier transform. For consistency between the two approaches, we use the following space continuous Fourier transform

[TABLE]

When the finite system and particle sizes need to be taken into account, we discretize the integrals into

[TABLE]

where $N=L^{2}/a^{2}$ is the number of particles, at unity packing fraction. In the sum, $\mathbf{q}$ takes discrete values defined by the geometry of the problem. For instance, for a square lattice of linear size $L$, $\mathbf{q}=(2\pi m/L,2\pi n/L)$ where $(m,n)$ are integers satisfying $0\leq m,n\leq L/a-1$. From this discretization, we get that the discrete space Fourier transform $\mathbf{u}(\mathbf{q},t)$ is consistently related to the continuous Fourier transform $\tilde{\mathbf{u}}(\mathbf{q},t)$ through

[TABLE]

This relation will be useful for comparison to the results of numerical simulations. In the following, we generically use the tilde notation for continuous Fourier transform, and drop the tilde when dealing with the discrete Fourier transform.

To proceed with the computations in the framework of the continuum theory, we now introduce the space and time Fourier transform

[TABLE]

With these definitions, the active equation of motion (S28) can be rewritten in Fourier space as

[TABLE]

where we have defined the continuous Fourier transform $\tilde{\mathbf{F}}^{\rm act}(\mathbf{q},\omega)$ of the random active force $\mathbf{F}^{\rm act}(\mathbf{r},t)$ in Fourier space as

[TABLE]

IV.3 Active noise correlations

To determine the correlation of the active noise, we need to start from a spatially discretized version of the model. For definiteness, we assume a square grid with lattice spacing $a$. Then for each grid node $i$ we have $\hat{\mathbf{n}}_{i}=(\cos\phi_{i},\sin\phi_{i})$ with dynamics $\dot{\phi}_{i}=\eta_{i}$, $\left\langle\eta_{i}\left(t\right)\eta_{j}\left(t^{\prime}\right)\right\rangle=\frac{2}{\tau}\delta_{ij}\delta\left(t-t^{\prime}\right)$, and the noise remains spatially uncorrelated. We thus have

[TABLE]

The exponential time dependence has been obtained using the same reasoning as in Eqs. (S5) to (S7). In order to take a continuum limit, we replace $\hat{\mathbf{n}}_{i}$ by a continuous field, and we substitute $\delta_{i,j}$ by its Dirac counterpart, namely

[TABLE]

We then have that, in the continuum limit,

[TABLE]

In view of Eq. (S36), it is clear that $\left\langle\tilde{\mathbf{F}}^{\rm act}\left(\mathbf{q},\omega\right)\right\rangle=0$, as $\left\langle\cos\phi\right\rangle=\left\langle\sin\phi\right\rangle=0$. The second order correlations are simply

[TABLE]

Using Eqs. (S39) and (S40), a straightforward calculation then yields

[TABLE]

Note that Eq. (S41) is obtained in the continuum formulation, where $\delta(\mathbf{q}+\mathbf{q}^{\prime})$ is a Dirac delta distribution, which is infinite if one sets $\mathbf{q}^{\prime}=-\mathbf{q}$. To compare with the numerics, one has to come back to the discrete formulation, corresponding to a finite system size $L$. The Dirac delta is then replaced by a Kronecker delta according to the substitution rule

[TABLE]

We are thus led to define the space-discrete Fourier transform $\mathbf{F}^{\rm act}(\mathbf{q},\omega)=\tilde{\mathbf{F}}^{\rm act}(\mathbf{q},\omega)/a^{2}$ [see Eq. (S32)] for discrete wavevectors $\mathbf{q}$ (note that $\omega$ remains a continuous variable). The correlation of the discrete Fourier transform $\mathbf{F}^{\rm act}(\mathbf{q},\omega)$ of the active noise then reads

[TABLE]

in agreement with Eq. (9) of the main text.

IV.4 Fourier modes properties

We decompose equation (S35) into longitudinal and transverse modes: $\tilde{\mathbf{u}}\left(\mathbf{q},\omega\right)=\tilde{u}_{\mathrm{L}}\left(\mathbf{q},\omega\right)\hat{\mathbf{q}}+\tilde{u}_{\mathrm{T}}\left(\mathbf{q},\omega\right)\hat{\mathbf{q}}^{\perp}$ along and perpendicular to the eigenvectors of the dynamical matrix, Eq. (S24). We obtain two equations

[TABLE]

with solution

[TABLE]

where $\tilde{F}_{\mathrm{L}}^{\rm act}\left(\mathbf{q},\omega\right)=\tilde{\mathbf{F}}^{\rm act}\left(\mathbf{q},\omega\right)\cdot\hat{\mathbf{q}}$ and $\tilde{F}_{\mathrm{T}}^{\rm act}\left(\mathbf{q},\omega\right)=\tilde{\mathbf{F}}^{\rm act}\left(\mathbf{q},\omega\right)\cdot\hat{\mathbf{q}}^{\perp}$.

We can use these expressions to obtain velocity correlation functions that can be directly measured in experiments and simulations. As $\tilde{\mathbf{v}}\left(\mathbf{q},\omega\right)=-\mathrm{i}\omega\tilde{\mathbf{u}}\left(\mathbf{q},\omega\right)$, we can simply write

[TABLE]

It is easy to show that the longitudinal and transverse components of the active force contribute equally to the correlation, namely

[TABLE]

Using Eqs. (S41), (S44) and (S45), the correlation functions of the longitudinal and transverse components of the Fourier velocity field are then straightforward to compute, leading to

[TABLE]

Of particular interest is the equal-time Fourier transform of the velocity. In other words, we need to integrate over frequency. E.g., for the longitudinal velocity, we find

[TABLE]

Using the decomposition

[TABLE]

a straightforward integration leads to

[TABLE]

A similar calculation for the transverse component of the Fourier velocity field yields

[TABLE]

Introducing the longitudinal and transverse characteristic length scales

[TABLE]

the equal-time (continuous) Fourier velocity correlation can be expressed as

[TABLE]

The length scales $\xi_{\mathrm{L}}$ and $\xi_{\mathrm{T}}$ can be interpreted as the longitudinal and transverse correlations lengths that both diverge $\sim\tau^{1/2}$ for $\tau\to\infty$ (i.e., a fully persistent self-propulsion). The existence of those correlation lengths is a direct consequence of activity. In the “passive” limit $\tau\to 0$, these length scales vanish.

It is important to note that Eq. (S52) is obtained in the continuum formulation, where $\delta(\mathbf{q}+\mathbf{q}^{\prime})$ is a Dirac delta distribution. Hence $\left\langle\tilde{\mathbf{v}}\left(\mathbf{q},t\right)\cdot\tilde{\mathbf{v}}\left(\mathbf{q}^{\prime},t\right)\right\rangle$ is infinite if one sets $\mathbf{q}^{\prime}=-\mathbf{q}$. To compare with the numerics, one has to come back to the discrete formulation, corresponding to a finite system size $L$. The Dirac delta is then replaced by a Kronecker delta according to the substitution rule given in Eq. (S42). One also needs to replace the continuum Fourier transform $\tilde{\mathbf{v}}\left(\mathbf{q},t\right)$ with the discrete one, $\mathbf{v}\left(\mathbf{q},t\right)$, according to $\tilde{\mathbf{v}}\left(\mathbf{q},t\right)=a^{2}\mathbf{v}\left(\mathbf{q},t\right)$ [see Eq. (S32)]. We thus end up with, using $N=L^{2}/a^{2}$,

[TABLE]

which is precisely Eq. (10) of the main text.

In addition, one can also compute (using integration techniques in the complex plane) the two-time Fourier velocity correlation $\left\langle\tilde{\mathbf{v}}\left(\mathbf{q},t\right)\cdot\tilde{\mathbf{v}}\left(\mathbf{q}^{\prime},t^{\prime}\right)\right\rangle$. This two-time correlation function is found to decay with the time lag $|t-t^{\prime}|$ over three different characteristic times, the persistence time $\tau$ of the noise and two elastic time scales $\tau_{\mathrm{L}}=\frac{\zeta}{\left(B+\mu\right)q^{2}}$ and $\tau_{\mathrm{T}}=\frac{\zeta}{\mu q^{2}}$ associated with longitudinal and transverse modes respectively.

IV.5 Mean-square velocity and velocity autocorrelation function

We conclude by computing the real-space mean-square velocity $\langle\left|\mathbf{v}(\mathbf{r},t)\right|^{2}\rangle$. One has

[TABLE]

Using Eq. (S52) we get

[TABLE]

This integral diverges at the upper boundary. This divergence can be regularized if we note that the physical upper limit to this integral is set by the inverse particle size, i.e., by $q_{\rm m}=\frac{2\pi}{a}$. Therefore, using $\int d^{2}q=2\pi\int q\,dq=\pi\int d(q^{2})$ when integrating a function of $q^{2}$, one obtains

[TABLE]

Note that $\left\langle\left|\mathbf{v}(\mathbf{r},t)\right|^{2}\right\rangle$ is independent of position (and time) and is thus also equal to

[TABLE]

Finally, generalizing the above calculation one can also compute the autocorrelation function of the velocity field, yielding

[TABLE]

IV.6 Real space expression of the velocity autocorrelation function

We derive here the real space expression for the correlation of velocities of cells separated by $r$. This is analytically tractable only if the continuum inverse Fourier transform is used, i.e., in the limit of infinite system size. However, this calculation has to be done with care, using the discrete Fourier transform and eventually taking the infinite volume limit to evaluate sums as integrals. Using instead the continuum Fourier transform of the velocity field would lead to difficulties because of the delta function $\delta(\mathbf{q}+\mathbf{q}^{\prime})$ in Eq. (S52). We define the real space correlation function of the velocity field as

[TABLE]

as well as its Fourier transform

[TABLE]

Note that the space integration is done on the finite volume $L^{2}$, so that the wavevector $\mathbf{q}$ is discretized. A straightforward calculation leads to

[TABLE]

where $\mathbf{v}(\mathbf{q})$ is the discrete Fourier transform of the velocity field, and $\langle|\mathbf{v}(\mathbf{q})|^{2}\rangle$ is given in Eq. (S53) as well as in Eq. (10) of the main text. From Eq. (S60), one can evaluate $C_{vv}(\mathbf{r})$ by computing the inverse Fourier transform of $C_{vv}(\mathbf{q})$. The inverse discrete Fourier transform of $C_{vv}(\mathbf{q})$ can be turned into an integral by taking the limit $L\to\infty$, yielding

[TABLE]

Using Eqs. (S60) and (S53), we obtain

[TABLE]

with $r=|\mathbf{r}|$, and $K_{0}$ the modified Bessel function of the second kind. To obtain Eq. (S62), we have made use of the following identities involving Bessel functions Gradshteyn

[TABLE]

where $J_{0}$ is the Bessel function of the first kind. An asymptotic expansion of Eq. (S62) for $r\gg\xi_{\mathrm{L},\mathrm{T}}$ yields

[TABLE]

that is, an exponential decay of $C_{vv}(\mathbf{r})$ at large distances, with algebraic corrections.

Supplementary Note 3

IV.1 Fitting to experiment and simulations

To compare simulations to our continuum predictions, we need to determine $B$ and $\mu$. As detailed in the Supplementary Note 2, we determine $D(\mathbf{q})$ by Fourier-transforming the dynamical matrix on the $\mathbf{q}$ grid appropriate to the simulations box. The longitudinal and transverse eigenvalues of the resulting $2\times 2$ matrix are then $(B+\mu)q^{2}$ and $\mu q^{2}$, respectively. In Supplementary Figure S1A, we show the radially $\mathbf{q}$-averaged eigenvalues (dots) as a function of $q^{2}$, and the linear fit of the $15$ first points we use to extract the moduli.

In Supplementary Figure S1B, we show the Self-Intermediate function as a function of time for the experiments and all three fitted simulations. For the experiment, we numerically integrated the PIV field to obtain approximate trajectories for the regions belonging to each individual PIV arrow at $t=0$. Significant local non-affine motion and distortions emerged, and we stopped before $t=10$ hours and at motions of a couple of cell diameters. The match between experiment and simulation is good for the soft disk simulations; the much slower dynamics of the vertex model is due to its much higher bulk modulus for a given shear modulus at $\bar{p}_{0}=3.6$.

Supplementary Figures

Additional numerical simulation results on the Fourier velocity correlations and their dependence on the persistence time $\tau$ are shown in Supplementary Figures S2 and S3. In addition, Supplementary Figure S4 shows the real-space correlation functions for soft disks (left) and the vertex model (right), together with the predictions of Eq. (14) of the main text.