Dense active matter model of motion patterns in confluent cell monolayers

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.


I. INTRODUCTION
Collective cell migration is of fundamental importance in embryonic development [1][2][3][4], organ regeneration and wound healing [5]. 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 [6,7]. Major advances in our understanding of collective cell migration have been obtained from in-vitro experiments on epithelial cell monolayers [4,[8][9][10][11]. 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 [12]. For example, forces in a monolayer are transmitted over long distances via a global tug-of-war mechanism [9]. The landscape of mechanical stresses is rugged with local stresses that are correlated over distances spanning multiple cell sizes [10]. These strong correlations lead to the tendency of individual cells to migrate along the local orientation of the maximal principal stress (plithotaxis [10,13]) and a tendency of a collection of migrating epithelial cells to move towards empty regions of space (kenotaxis [14]). Furthermore, such coordination mechanisms lead to propagating waves in confined clusters [15,16], expanding colonies [17] and in colliding monolayers [18], which all occur in the absence of inertia.
Active matter physics [19,20] 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 [21], oscillations [22], active liquid crystalline [20], and arrested, glassy phases [23]. In-silico studies [22,[24][25][26][27] in the dense regime have been instrumental in describing and classifying experimentally observed collective active motion.
Continuum active gel theories [28,29] are able to capture many aspects of cell mechanics [30], including spontaneous flow of cortical actin [31] and contractile cell traction profiles with the substrate [32]. In some cases, cell shapes form a nematic-like texture [4,33] and topological defects present in such texture have been argued to assist in the extrusion of apoptotic cells [4]. 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 [34] known as dynamic heterogeneities [35], a hallmark of the glass transition [36] 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 [37][38][39] or cell-cell adhesion [40], 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 [41,42], 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 selfpropelled Voronoi model (SPV) [37,38], 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. [40] 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. [40].

A. 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 cellsubstrate friction. The equations of motion for cell centers are where ζ is the cell-substrate friction coefficient, F act i is the net motile force resulting from the cell-substrate stress transfer, and F int i is the interaction force between cell i and its neighbours. Commonly used interaction models are short-ranged pair forces with attractive and repulsive components [43] and SPV models [37,38]. 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, F int i = −∇ ri V ({r j }). 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 F act i 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 F act i = F actn i , wheren i is the unit vector that makes an angle θ i with the x−axis of the laboratory frame, the angular dynamics isθ where τ sets a persistence time scale, and different cells are not coupled (Fig. 1A). This dynamics is equivalent to active Brownian particles [44], and in isolation, model cell motion is a persistent random walk. At sufficiently low driving, such models form active glasses [23,[45][46][47], 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 τ α , which diverges at dynamical arrest. We now develop a linear response formalism. As shown in Fig. 1A, on time scales below τ α , the self-propulsion reduces to a stochastic, time-correlated force fluctuating inside a local energy minimum. If the persistence time scale τ τ α , 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. [47]). We linearize the interaction forces in the vicinity of an energy minimum, i.e., a mechanically stable or jammed configuration {r 0 i } by introducing δr i = r i − r 0 i . After introducing the active velocity v 0 = F act /ζ, Eq. (1) becomes [48], organised as 2 × 2 blocks corresponding to cells i and j. In this limit, we can solve the dynamics exactly, see Supplementary Note 1.

B. Normal mode formulation
Assuming that there are a sufficient number of intercell forces to constrain the tissue to be elastic at short  time scales, the dynamical matrix has 2N independent normal modes ξ ν with positive eigenvalues λ ν . If we project Eq. (3) onto the normal modes, we obtain where a ν = i δr i · ξ ν i and the self-propulsion force has been projected onto the modes, η ν = ζv 0 in i · ξ ν i . The self-propulsion then acts like a time-correlated Ornstein-Uhlenbeck noise (see Supplementary Note 1), with η ν (t)η ν (t ) = 1 2 ζ 2 v 2 0 exp(−|t − t |/τ )δ ν,ν . We can integrate Eq. (4) and obtain the moments of a ν . In particular, the mean energy per mode is given by explicitly showing that equipartition is broken due to the mode-dependence induced by λ ν in Eq. (5). In the limit τ → 0, we recover an effective thermal equilibrium, E ν → ζv 2 0 τ /4 := T eff /2, where T eff = ζv 2 0 τ /2, consistent with previous work [37,45]. In Fig. 2A, the left-most column is for a simulated thermal system, with properties that are nearly indistinguishable from the τ = 0.2 results. In the opposite, high persistence limit when τ → ∞, we obtain instead E ν = ζ 2 v 2 0 /4λ ν , 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. [22,37,49].
It thus becomes clear that for large values of τ , T eff can no longer be interpreted as temperature since the fluctuation-dissipation theorem is no longer valid. An analogous result to the τ → ∞ limit has been obtain in granular material with an externally applied shear [41,42], showing that the mechanisms at play are generic, and that tuning τ allows active systems to bridge between features of thermal systems and (self-)sheared systems (see also ref. [37,47]). However, the glass transition lies on a curve of constant T eff with moderate τ contributions ( Fig. 2 and ref. [46]), making T 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) [50]. We compute the Fourier space velocity correlation function, |v(q)| 2 = v(q) · v * (q) , with v(q) = 1/N N j=1 e iq·r 0 j δṙ j , where the {r 0 j } 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 derivativeṡ a ν of the modes amplitudes for different modes, we first derive (see Supplementary Note 1) the mode correlations ȧ 2 ν = v 2 0 / [2 (1 + λ ν τ /ζ)], so that in Fourier space, we obtain where ξ ν (q) is the Fourier transform of the vector ξ ν i .

C. 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 [51,52], 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 [53]. We rewrite Eq. (3) as where u(R) denotes the elastic deformations from the equilibrium positions R in the solid, and D(R − R ) is the continuum dynamical matrix. The normal modes of the system are now simply Fourier modes with where F act (q, ω) = ζv 0 ∞ −∞ dt Rn (R, t)e iωt e iq·R and D(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 noisen(R, t) is spatially uncorrelated, and we find the noise correlators (see Supplementary Note 2) The dynamical matrix D(q) has two independent eigenmodes in two dimensions, one longitudinalq with eigenvalue B + µ and one transverse oneq ⊥ with eigenvalue µ, where B and µ are the bulk and shear moduli, respectively. We can then decompose our solution into longitudinal and transverse parts, u(q, ω) = u L (q, ω)q + u T (q, ω)q ⊥ . We are interested in the equal time, Fourier transform of the velocity, which we find to be (see Supplementary Note 2) where we have introduced the longitudinal and transverse correlation lengths ξ 2 L = (B + µ) τ /ζ and ξ 2 T = µτ /ζ. 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, [µ/ζ] = [B/ζ] = 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 |v| 2 = 1 N i |v i | 2 decreases with active correlation time as where q m = 2π/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 √ τ . In the limit τ → ∞, |v(q)| 2 diverges at low q, as was found in ref. [54]. The dominant scaling |v| 2 ∼ 1/ξ 2 is the same as results from the scaling Ansatz for cell-substrate dominated coordinated motion obtained by ref. [40].
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 In the infinite size limit L → ∞, the real space correlation function C vv (r) can be evaluated from the Fourier correlation |v(q)| 2 as Using Eq. (11), one finds the explicit result (see Supplementary Note 2) 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/ξ L,T 1, i.e. for distance much larger than the correlation lengths, i.e., as expected and consistent with the results of ref. [40], C vv decays exponentially at large distances.

D. 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 = v 0ni and pair interaction forces F ij that are purely repulsive. We simulate a confluent sheet in this model by setting the packing fraction to φ = 1 in periodic boundary conditions. The SPV model is the same as introduced in refs. [37,38], and assumes that every cell is defined by the Voronoi tile corresponding to its centre. For this model, we choose the dimensionless shape factorp 0 = 3.6, putting the passive system into the solid part of the phase diagram [37,55], 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 eff = ζv 2 0 τ /2 has emerged as a good predictor of the active glass transition [45,54], at least at low τ , and we use it together with τ itself as the axes of our phase diagram. The liquid or glassy behaviour of the model can be characterised by the alpha relaxation time τ α . Fig. 2 provides a coarse-grained phase diagram where τ α is represented in gray scale as a function of persistence time τ and T 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 τ increases. This non-trivial result is consistent with the recent RFOT theory of the active glass transition [46] and related simulation results [46,47]. It can be partly understood from the fact that v 0 decreases when τ increases at fixed T 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 [54] only apply in the small τ 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 τ is readily apparent as swirl-like motion (see also Supplementary Movies 1-4). Fig. 3A-C shows the Fourier velocity correlation |v(q)| 2 measured in the numerical simulations for different values of v 0 , after normalizing |v(q)| 2 by v 2 0 N . In panel A, we show that for soft disks at low T eff = 0.005, where the system is solid, the correlation function develops a dramatic 1/q 2 slope as τ 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 ± 0.008, µ = 0.510 ± 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 τ → 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 µ = 0.5, as estimated from ref. [55] forp 0 = 3.6. We did not compute normal modes for the SPV model. Note that due to µ µ+B, the contribution of the transverse correlations dominate the analytical results in both cases. In panel B, we show the soft disk simulation at τ = 20 when the transition to a liquid is crossed as a function of T eff . Deviations from the normal mode predictions become apparent only at the two largest values of T eff , when τ > τ α ( 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 τ and both soft disk and SPV models, our predictions remain in excellent agreement with the simulations for T eff = 0.02, where τ τ α . In Fig. 3D, we show the mean square velocity normalized by v 0 as a function of the dimensionless transverse correlation length ξ T /a ∼ √ τ , for all our simulations, using a = σ, 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.

E. 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 ofv = |v(r, t)| 2 = 12 ± 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 usev(r, t) = v(r, t)/ |v(r, t)| 2 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 A ≈ 380 µm 2 corresponding to a particle radius of σ ≈ 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 |v(q)| 2 , shown in Fig. 4B. Using the results from Eq. (11) we can rewrite Eq. (10) as where ξ 2 L and ξ 2 T are the longitudinal and transverse squared correlation lengths (with units of µm 2 ) defined below Eq. (10). As can be seen from Eq. (11), the ratio v0 v is a function only of the dimensionless ratios ξ L /a and ξ T /a. If we further make the plausible assumption that the ratio of elastic moduli is the same as in the soft disk simulations, (µ + B)/µ = 4.3, the correlation lengths ξ L and ξ T are not independent. Therefore, we are left with a single fitting parameter, ξ T . The best fit to the theory is obtained with ξ 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 /v ≈ 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 ξ T /a ≈ 5, v /v 0 ≈ 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 ξ 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 v(t)·v(0) 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 m . Using the value of ξ 2 T = 10 4 µm 2 extracted from fitting |v(q)| 2 and the ratio (µ + B)/µ = 4.3, we used different µ/ζ and τ compatible with ξ 2 T = µσ 2 /ζτ to obtain the best (numerically integrated) analytical fit to the experimental result. We settled on a best fit autocorrelation time of τ = 2.5 h and µ/(σ 2 ζ) = 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 τ 5 = 20 h and the same µ/ζ (light gray line). This is consistent with the longer spatial correlations observed in Fig. 4D, where the light grey line corresponds to ξ T = 100 µm τ 5 /τ ≈ 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 forv and the ratio v 0 /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 ξ and τ , 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 µ chosen with the same ratio as in the previous particle (respectively, SPV) simulations. From these values, we can approximate the parameter values k/ζ = µ/σ 2 for the particle model and K/ζ = µ/ A 2 , Γ/ζ = µ/ A 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/ζ = Γ/ζ = 55 h −1 , K/ζ = 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. [40] for confluent but still motile cells, with a reported maximum correlation length of ξ = 100 µm, and a cell crawling speed that drops by a factor of 10 from ∼ 90 µm h −1 at low density to this point of maximal correlation. Note also that we have an elastic time scale (k/ζ) −1 ≈ 0.02h, much shorter than our correlation time scale τ = 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 48h 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. [56,57]. With a typical cell cycle time of τ div = 48h, 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 τ and division time scale τ div . We have also considered the effect of weak polar alignment between cells, using the model from ref. [22]. For weak alignment with time scales τ v ≥ τ that do not lead to global flocking of the sheet, which we do not observe, |v(q)| 2 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 [58]) 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 τ 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 [59]. More generally, including cellcell dissipation in addition to cell-substrate dissipation could significantly modify scalings [40], 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 [22,49]. From a fundamental point of view, our theoretical results (and also the results of ref. [22]) are examples of a larger class of non-equilibrium steady-states that can be treated using a linear response formalism [60].

A. 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 • C prior to plating and the cells were kept in a humidified incubator at 37 • C and an atmosphere of 5% CO 2 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 • C. The automated time-lapse imaging setup took an image at 10 minute intervals at a magnification of 10x, corresponding to a field of view of 867 µm × 662 µm that was saved at a resolution of 1300 × 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 [61]. At the finest resolution (level 3), it provides displacement vectors on a 54 × 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.

B. 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 [62]. The interaction potential for soft harmonic disks is V i = j k 2 (σ i + σ j − |r j −r i |) 2 if |r j −r i | ≤ σ i +σ j and 0 otherwise. To emulate a confluent cell sheet, we used periodic boundary conditions at packing fraction φ = 1, where φ = i πσ 2 i /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 (φ > 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 where A i is the area of the tile, and P i is its perimeter, K and Γ are the area and perimeter stiffness coefficients and A 0 and P 0 are reference area and perimeter, respec-tively. SPV is confluent by construction, and its effective rigidity is set by the dimensionless shape parameter p 0 = P 0 / √ A 0 , with a transition from a solid to a fluid that occurs forp 0 ≈ 3.812. We simulate the model at p 0 = 3.6, well within the solid region at zero activity [37,55], and also introduce 30% variability in A 0 . AVM was implemented with open boundary conditions, and we use a boundary line tension λ = 0.3 to avoid a fingering instability at the border that appears especially at large τ .
Both models are simulated with overdamped active Brownian dynamics Equations of motions are integrated using a first order scheme with time step δt = 0.01. Simulations are 5 × 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.

C. 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 q and configuration, the correlation function is |v(q)| 2 = v(q) · v * (q), of which we then take a radial 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 × 40 grid points as normalization. We compute the α-relaxation time from the self-intermediate scattering function S(q, t) = 1 N j e iq·(rj (t0+t)−rj (t0) t0,|q|=q , where the angle brackets indicate time and radial averages. At q = 2π/σ, we determine τ α as the first time point where S(q, t) < 0.5, bounded from above by the simulation time.

D. Normal mode analysis
The normal modes are the eigenvalues and eigenvectors of the Hessian matrix K ij = ∂ 2 V ({r i })/∂r i ∂r j , evaluated at mechanical equilibrium. We first equilibrate the t = 2500 snapshot with v 0 = 0 for 2 × 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 τ = 2000). For the soft disk model, each individual ij contact with contact normaln ij and tangentialt ij vectors contributes a term K ij = −kn ij ×n ij + |f ij |t ij ×t ij to the ij 2 × 2 off-diagonal element of the matrix and −K ij is added to the ii diagonal element [48]. 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 ν through ξ ν (q) = 1 N j e iq·rj ξ j ν and then |ξ ν (q)| 2 = ξ ν (q) · (ξ ν (q)) * . The 2 × 2 continuum Fourier space dynamical matrix D(q) has one longitudinal eigenmode alongq with eigenvalue (B + µ)q 2 and one transverse eigenmode alongq ⊥ with eigenvalue µ. We compute D(q) αβ = 1 N j l e iqαrj,α H jl,αβ e −iq β r j,β , where the greek indices α, β 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 ontoq, and from there the longitudinal and transverse eigenvalues λ L (q) and λ T (q) after a radial q average. We fit B + µ as the slope of λ L (q) vs q 2 up to q = 1.5, and the same for µ and λ T (q) (Supplementary Figure 3).

DATA AVAILABILITY
Data supporting the findings of this manuscript are available from the corresponding authors upon reasonable request.  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 r 0 i , corresponding to a local minimum of the elastic energy. If the interaction potential is linearized around the energy minimum in terms of the displacement δr i = r i − r 0 i , the dynamics is described by the equation where the K ij 's are the 2 × 2 blocks of the 2N × 2N dynamical matrix, v 0ni is the self-propulsion term witĥ n i = cos φ i e x + sin φ i e y (i.e., direction ofn is given by the angle φ i with the x axis of a laboratory reference frame) and ζ is the friction coefficient. In the absence of inter-particle alignment, the angle φ i obeys a simple rotational diffusive dynamics with white noise η i (t): where we have expressed the inverse rotational diffusion constant as a time scale, τ = 1/D r . We note that in general, the system is far out of thermodynamic equilibrium and D r and ζ 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 δr i . To this aim, we first expand δ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 twodimensional vectors (ξ ν 1 , . . . , ξ ν N ), where the index ν = 1, . . . , 2N labels the mode; the associated eigenvalue is denoted as λ ν . This form of the normal modes is useful as it allows the decomposition of δr i to be written in the simple form Projecting Eq. (S1) on the normal modes, we find the uncoupled set of equations is the projection of the self-propulsion force onto the normal mode ν.

B. Self-propulsion force as a persistent noise
We consider the projection η ν of the self-propulsion force on normal mode ν as a correlated noise, which we now characterize. Since η ν 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 φ i , that η ν (t) = 0. We thus simply need to evaluate the two-time correlation function of η ν (t). Using the fact that the eigenvectors of the dynamical matrix form an orthonormal basis, we have where φ (t) 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 . We can thus set t = 0 without loss of generality. Solving Eq. (S2), the quantity ∆φ = φ (t) − φ (0) is distributed according to One then finds, using Eqs. (S5) and (S6), i.e., the time correlation of the noise η ν decays exponentially with the correlation (or persistence) time τ . It is worth emphasizing that the statistical properties of the noise η ν are independent of the mode ν.

C. 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 η ν (t), one finds From this expression, one can compute the average value a 2 ν (t) , leading for t → ∞ to Using Eq. (S7), we obtain or, in terms of average energy per mode (S11) For very short correlation time τ (i.e., large diffusion coefficient D r ), one recovers an effective equipartition of energy over the modes, E ν ≈ ζv 2 0 τ 4 even though the system is out-of-equilibrium. For finite correlation time, this result remains valid in the range of modes ν such that τ ζλ −1 ν , if such a range exists. However, for large correlation time τ , that is, as soon as there is a wide range of modes such that τ ζλ −1 ν , equipartition is broken, and the energy spectrum is given by E ν ≈ ζ 2 v 2 0 4λν .

D. Velocity correlation
Following [S3], we consider the velocity-velocity correlation functionĜ(q) in Fourier space, where one can express the (discrete) Fourier transform v (q) as a function of the particles reference positions r 0 i : where the star denotes the complex conjugate. Expanding over the normal modes, one findŝ where ξ ν (q) is the Fourier transform of the vectors ξ ν . From Eq. (S4), the quantity ȧ νȧν is expressed as where the last equality is due to the modes being uncorrelated. The cross-correlation is in fact not 0, but crucial: (S15) To sum up, one has according to Eq. (S14) with (S17) Further, using Eqs. (S5), (S7), (S10) and (S15), one obtains (S18) Combining Eqs. (S13), (S16) and (S18), we derive the final expression for the velocity correlation function: (S19) Note that we can compute the equal-time, spatial mean square velocity through Parseval's theorem as

A. 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 µ can be written as [S1, S2] whereû is the strain tensor with components u αβ = 1 2 [∂ α u β + ∂ β u α ] written as spatial derivatives of the components α, β ∈ {x, y} of the displacement vectors u (r) = r (r) − r from a reference state r to the deformed state r (r). The stress tensor σ αβ = δF el δu αβ can then be written as where summation over pairs of repeated indices is assumed. Hence, its divergence is given by We can then write the overdamped equations of motion for the displacement field This last equation can be rewritten in vectorial notation as In Fourier space, we can write this relation as where D (q) is the Fourier space dynamical matrix, and q 2 = q 2 x + q 2 y . The two eigenvalues of the dynamical matrix are with normalized eigenvectors In other words, for each q, we obtain one longitudinal and one transverse eigenmode, with diffusive equations of motion, where the diffusion coefficients are the two elastic moduli:

B. Overdamped dynamics with activity
Now including the self-propulsion force, the continuum version of the active equations of motion is given by where we have included an active force F act (r, t) = ζv 0n (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 → ∞ and a → 0, i.e., using the continuous Fourier transform. For consistency between the two approaches, we use the following space continuous Fourier transform When the finite system and particle sizes need to be taken into account, we discretize the integrals into where N = L 2 /a 2 is the number of particles, at unity packing fraction. In the sum, q takes discrete values defined by the geometry of the problem. For instance, for a square lattice of linear size L, q = (2πm/L, 2πn/L) where (m, n) are integers satisfying 0 ≤ m, n ≤ L/a − 1. From this discretization, we get that the discrete space Fourier transform u(q, t) is consistently related to the continuous Fourier transformũ(q, t) through u(q, t) = a 2 u(q, t).
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 With these definitions, the active equation of motion (S28) can be rewritten in Fourier space as where we have defined the continuous Fourier transformF act (q, ω) of the random active force F act (r, t) in Fourier space asF (S36)

C. 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 haven , and the noise remains spatially uncorrelated. We thus have 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 replacen i by a continuous field, and we substitute δ i,j by its Dirac counterpart, namely We then have that, in the continuum limit, n(r, t) ·n(r , t ) = a 2 δ(r − r ) e −|t−t |/τ .
In view of Eq. (S36), it is clear that F act (q, ω) = 0, as cos φ = sin φ = 0. The second order correlations are simply Using Eqs. (S39) and (S40), a straightforward calculation then yields Note that Eq. (S41) is obtained in the continuum formulation, where δ(q + q ) is a Dirac delta distribution, which is infinite if one sets q = −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 We are thus led to define the space-discrete Fourier transform F act (q, ω) =F act (q, ω)/a 2 [see Eq. (S32)] for discrete wavevectors q (note that ω remains a continuous variable). The correlation of the discrete Fourier transform F act (q, ω) of the active noise then reads in agreement with Eq. (9) of the main text.

D. Fourier modes properties
We decompose equation (S35) into longitudinal and transverse modes:ũ (q, ω) =ũ L (q, ω)q +ũ T (q, ω)q ⊥ along and perpendicular to the eigenvectors of the dynamical matrix, Eq. (S24). We obtain two equations whereF act L (q, ω) =F act (q, ω) ·q andF act T (q, ω) =F act (q, ω) ·q ⊥ . We can use these expressions to obtain velocity correlation functions that can be directly measured in experiments and simulations. Asṽ (q, ω) = −iωũ (q, ω), we can simply write It is easy to show that the longitudinal and transverse components of the active force contribute equally to the correlation, namely 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 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 Using the decomposition a straightforward integration leads to ṽ L (q, t)ṽ L (q , t) = 2π 2 a 2 ζv 2 0 (B + µ)τ q 2 + ζ δ(q + q ).
It is important to note that Eq. (S52) is obtained in the continuum formulation, where δ(q + q ) is a Dirac delta distribution. Hence ṽ (q, t) ·ṽ (q , t) is infinite if one sets q = −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ṽ (q, t) with the discrete one, v (q, t), according toṽ (q, t) = a 2 v (q, t) [see Eq. (S32)]. We thus end up with, 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 ṽ (q, t) ·ṽ (q , t ) . This two-time correlation function is found to decay with the time lag |t − t | over three different characteristic times, the persistence time τ of the noise and two elastic time scales τ L = ζ (B+µ)q 2 and τ T = ζ µq 2 associated with longitudinal and transverse modes respectively.

E. Mean-square velocity and velocity autocorrelation function
We conclude by computing the real-space mean-square velocity |v(r, t) 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 m = 2π a . Therefore, using d 2 q = 2π q dq = π d(q 2 ) when integrating a function of q 2 , one obtains Note that |v(r, t)| 2 is independent of position (and time) and is thus also equal to Finally, generalizing the above calculation one can also compute the autocorrelation function of the velocity field, yielding F. 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 δ(q + q ) in Eq. (S52). We define the real space correlation function of the velocity field as C vv (r) = 1 L 2 d 2 r 0 v(r 0 + r) · v(r 0 ) (S58) as well as its Fourier transform C vv (q) = d 2 r C vv (r) e iq·r .
Note that the space integration is done on the finite volume L 2 , so that the wavevector q is discretized. A straightforward calculation leads to where v(q) is the discrete Fourier transform of the velocity field, and |v(q)| 2 is given in Eq. (S53) as well as in Eq. (10) of the main text. From Eq. (S60), one can evaluate C vv (r) by computing the inverse Fourier transform of C vv (q). The inverse discrete Fourier transform of C vv (q) can be turned into an integral by taking the limit L → ∞, yielding Using Eqs. (S60) and (S53), we obtain with r = |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 [S4] 1 2π π −π dθ e ix sin θ = J 0 (x), that is, an exponential decay of C vv (r) at large distances, with algebraic corrections.
SUPPLEMENTARY NOTE 3

A. Fitting to experiment and simulations
To compare simulations to our continuum predictions, we need to determine B and µ. As detailed in the Supplementary Note 2, we determine D(q) by Fourier-transforming the dynamical matrix on the q grid appropriate to the simulations box. The longitudinal and transverse eigenvalues of the resulting 2 × 2 matrix are then (B + µ)q 2 and µq 2 , respectively. In Supplementary Figure S1A, we show the radially 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 atp 0 = 3.6.