Energetics of mesoscale cell turbulence in two-dimensional monolayers

Investigation of energy mechanisms at the collective cell scale is a challenge for understanding various biological processes, such as embryonic development and tumor metastasis. Here we investigate the energetics of self-sustained mesoscale turbulence in confluent two-dimensional (2D) cell monolayers. We find that the kinetic energy and enstrophy of collective cell flows in both epithelial and non-epithelial cell monolayers collapse to a family of probability density functions, which follow the q-Gaussian distribution rather than the Maxwell–Boltzmann distribution. The enstrophy scales linearly with the kinetic energy as the monolayer matures. The energy spectra exhibit a power-decaying law at large wavenumbers, with a scaling exponent markedly different from that in the classical 2D Kolmogorov–Kraichnan turbulence. These energetic features are demonstrated to be common for all cell types on various substrates with a wide range of stiffness. This study provides unique clues to understand active natures of cell population and tissues. Understanding the mechanisms underlying collective motion in cells is central to understanding tissue development as well as the emergence of several pathologies. Here, the authors study the mesoscale turbulence of confluent cell monolayers through a combined experimental and numerical approach, finding that energetic statistics are independent of cell types and substrate stiffness.

Due to the active motility of cells and complex intercellular interactions, many nonequilibrium dynamic features emerge spontaneously in living cell collectives. For example, migrating cells can self-organize into mesoscale cell turbulence [16][17][18][19][20][21][22] , which is crucial in multicellular organization 23 , nutrient mixing 24 , and thrombopoiesis 25 . Though the mesoscale cell turbulence apparently resembles the classical turbulent flows, it takes place at an ultralow Reynolds number where the inertial effect can be negligible 26,27 , profoundly different from the passive turbulence resulting from excessive kinetic energy and inertia. However, the energetics of mesoscale cell turbulence remain unclear. Studying the energetic features in multicellular systems could advance our understanding of tissue morphodynamics, as well as the nonequilibrium statistical properties in active matter systems.
In this article, we combine large-scale experiments and theoretical analysis to investigate the energetic statistics of mesoscale cell turbulence in two-dimensional (2D) multicellular monolayers. By live-cell imaging on diverse cell monolayer systems, we reveal some unique energetic features of mesoscale cell turbulence. The enstrophy scales proportionally to the kinetic energy, which generalizes previous findings 23, 28 to monolayers of diverse cell lines, suggesting a common feature in living multicellular systems. The kinetic energy as well as the enstrophy over time collapse to a family of probability density functions (PDFs), which follow the q-Gaussian statistics. The energy spectra of mesoscale cell turbulence are profoundly different from that in the classical 2D Kolmogorov-Kraichnan turbulence. These features are not only common across different cell types, but also insensitive to the substrate stiffness. An active vertex model is used to reveal the physical mechanisms underlying these energetic statistics.

Results
Experiment. In our experiments, 2D cell monolayers consisting of cell lines including Madin Darby canine kidney (MDCK), human umbilical vein endothelial cells (HUVECs), C2C12 mouse myoblast cells (C2C12), and NIH-3T3 mouse embryo fibroblast cells (NIH-3T3) are investigated. Cells are seeded on 35 mm Petri dishes and cultured for enough time (typically 1 day) to grow into gap-free monolayers with a density of $ 1:4 10 5 cm À2 . Subsequently, live-cell imaging is performed to obtain consecutive phase contrast images ( Fig. 1a; "Methods"). Hereafter, the concept of time in experiments refers to time after live-cell imaging start. The field of view (FOV) has size 1.33 mm × 1.33 mm, which is much larger than the spatial correlation length (~200-300 µm) of mesoscale cell turbulence 16,18,19 , and thus large enough to access its energetic statistics.
Mesoscale cell turbulence in various cell monolayer systems. To quantitatively investigate the collective cell flows in a confluent cell monolayer, we calculate the velocity field via the particle image velocimetry (PIV) analysis ("Methods"; for the accuracy validation, see Supplementary Note 1 and Supplementary Fig. 1). Figure 1b, c shows representative velocity fields and corresponding vorticity fields in the four kinds of cell monolayers. Interestingly, turbulent motion patterns are observed for all the four kinds of cell monolayer systems throughout the live-cell imaging experimental duration of 12 h, in spite of the distinct difference of cell motility in these cell lines. This suggests that turbulent flows are ubiquitous in 2D cell monolayer systems.
Energetics of mesoscale cell turbulence in MDCK cell monolayers. We first study the energetics of the mesoscale turbulence in a confluent MDCK cell monolayer. Based on the cell velocity and vorticity fields, we calculate the kinetic energy E x; t ð Þ ¼ where v x; t ð Þ is the velocity vector and ω x; t ð Þ is the vorticity 23,28 . Figure 2a, b shows the representative kinetic energy and the enstrophy field in a MDCK cell monolayer, demonstrating distinctively spatial heterogeneity. We thus examine the PDF of the kinetic energy.
where Á h i x;t stands for an average over the entire FOV throughout 1 h period of time, during which cell motility keeps nearly invariant. The PDFs of the rescaled kinetic energy at different time points collapse to a family of PDFs which vary little over time (Fig. 2c), though the average kinetic energy decreases considerably with time. This suggests timeinvariant statistics of collective cell flows during the jamming transition of cellular motions as the monolayer matures 18,29 . These rescaled PDFs display a fat tail distribution characteristics, which deviates distinctly from the classical Maxwell-Boltzmann distribution of energy, p BẼ À Á ¼ exp ÀẼ À Á (Fig. 2c), evidencing that the examined cell monolayer system is nonequilibrium. The nonequilibrium feature of living cell systems stems from the active motility of cells.
In our previous study, we have explored the distribution statistics of the velocity field. It was found that cell velocities satisfy the 2D q-Gaussian distribution 30 (Supplementary Fig. 2a), depend on the Tsallis index q, i.e., the entropic index [31][32][33] ; Γ Á ð Þ is the Legendre gamma function. In statistical mechanics, the q-Gaussian distribution of velocities results from optimizing the Tsallis entropy 31,34 . This indicates that the collective cell dynamics is governed by the Tsallis entropy rather than the classical Boltzmann-Gibbs entropy. Based on the q-Gaussian distribution feature of velocities, we can deduce that the PDF of the kinetic energy should obey the following distribution function (also referred to as the q-Gaussian distribution hereafter for simplicity), As expected, our experimental data clearly demonstrate that the PDF of the kinetic energy can be well captured by the q-Gaussian distribution, i.e., Eq. (2), as shown in Fig. 2c. Similarly, we find that cell vorticities satisfy the 1D q-Gaussian distribution 30 (Supplementary Fig. 2b), whereω ¼ ω=ω rms is the rescaled vorticity with ω rms ¼ ffiffiffiffiffiffiffiffiffi ω 2 h i p being the root-mean-square vorticity; q Γðλ q À 1=2Þ and D q ¼ 1=ð2λ q À 3Þ. Based on the q-Gaussian distribution feature of vorticities, we can obtain the distribution law of the enstrophy Ω (Fig. 2d), which obeys the following q-Gaussian statistics: whereΩ ¼ Ω= Ω h i is the rescaled enstrophy. In addition, to verify whether the landscape of the energetic PDFs is affected by the PIV calculation, we have varied the window size of PIV analysis (from 1.5 cell size to 3 cell size) and found that it has no marked effect on the PDFs of the rescaled kinetic energy and the rescaled enstrophy ( Supplementary Fig. 3).
The kinetic energy and the enstrophy of the whole cell monolayer can be quantified by E t where Á h i x stands for an average over the entire FOV. Figure 2e plots the kinetic energy E(t) versus the enstrophy Ω(t) at different time points during the experimental duration. It shows that both E(t) and Ω(t) reduce over one order of magnitude within 10 h, indicating a jamming transition of cellular motions. Intriguingly, the enstrophy Ω(t) scales linearly with the kinetic energy E(t) during the jamming transition. This scaling behavior defines a conserved characteristic length scale ξ a ¼ ffiffiffiffiffiffiffiffiffi E=Ω p , which stems from the force balance between active and passive forces and roughly reflects the size of vortices 23,28,35 . For the MDCK cell monolayer, our measurement gives ξ a % 44:1 μm, approximately corresponding to a length scale of two cells. We note that a similar scaling behavior of E versus Ω has been found in other biological systems, such as human bronchial epithelial cell monolayers 23 and bacterial suspensions 28 , revealing a potentially conserved energetic property in these active matter systems.
We noted that previous theoretical studies have predicted some universal scaling laws for active nematic turbulence 36,37 . For example, there exists power scaling laws of the kinetic energy spectrum and the enstrophy spectrum. We thus seek to examine whether these energetic statistics hold for diverse cell monolayer systems. To gain further insights into the energetic landscape of mesoscale cell turbulence, we examine the energy spectrum E(k) of collective cell flows. E(k) is related to the velocity field by v r; t ð Þ j j 2 r;t ¼ 2 Þdk, and thus reflects the accumulation of kinetic energy over different spatial scales. In 2D space, E(k) can be calculated by where k ¼ k j j is the wavenumber. Figure 2f shows that the energy spectra E(k) at different time points exhibit a similar structure in the wavenumber space, suggesting a high robustness of energetic statistics over time. Our experimental data suggest two regimes that exhibit different behaviors. Let l c denote the average cell size. Then k c ¼ 2π=l c represents the wavenumber at the scale of one cell and, thus, sets the upper bound for the spectral range of the cell turbulence. At small wavenumbers, E(k) slightly increases with k. Previous theoretical studies predicted a scaling of E k ð Þ $ k À1 at small wavenumbers 36,37 . However, it seems it's not the case here. More data points (i.e., larger FOV size) are needed to identify the specific scaling law here. At large wavenumbers k=k c > 0:1, in contrast, E(k) displays an overall decay with increasing k. We observe three energy plateaus for wavenumbers (i) 0:06 < k=k c < 0:15, (ii) 0:4 < k=k c < 0:6, and (iii) 0:8 < k=k c < 1:0 (Fig. 2f). The first plateau (i) occurs at length scales ranging from 10 cells to 17 cells, which corresponds to the scale of swirls 16,18,19 . The second (ii) and the third plateau (iii) stem from the filtering and averaging effects of PIV calculation (see Supplementary Note 2). When the PIV window size is varied, these two plateaus shift in the wavenumber space, but the scaling behavior of the energy spectrum does not vary significantly ( Supplementary  Fig. 4). Therefore, the two plateaus will be ignored hereafter. The energy spectra obey a power scaling as E k ð Þ $ k Àβ for large wavenumbers 0:2 < k=k c < 0:4 and 0:6 < k=k c < 0:8. These two decay scaling regimes possess the same exponent β and are reminiscent of the "inertial subrange" for classical passive turbulence. The power-law exponent is fitted to be β ¼ 4:55 ± 0:17 ( Fig. 2f), which, nevertheless, strikingly differs from the classical k À5=3 -decay of 2D Kolmogorov-Kraichnan turbulence. This difference may be attributed to the different mechanisms of energy injection and dissipation. Compared to classical turbulence, mesoscale cell turbulence is powered by energy production of individual cells due to intracellular biochemical reactions and coordinated between neighbors by intercellular mechanical linkages, such as cadherin-mediated cell-cell junctions 38,39 . Similar energy-scaling laws but with different power-law exponents have been found in other active matter systems. For instance, a k À8=3 -decay scaling law was found in dense bacterial suspensions 40 , whereas a k À4 -decay scaling law was predicted for active nematics 36,37 . In addition, we have also examined the enstrophy spectra Ω(k) and observed similar scaling behavior (Fig. 2g). However, our experiments show that the enstrophy spectrum does not follow the law Ω k ð Þ $ k 2 E k ð Þ, which was predicted by active nematic theory 36,37 .
Energetics of mesoscale cell turbulence across diverse cell monolayers. Are these energetic characteristics of mesoscale cell turbulence specified to MDCK cell monolayers? To address this question, we examine the energetic landscape in cell monolayers of HUVECs, C2C12 cells, and NIH-3T3 cells. Across these systems, both the enstrophy Ω and the kinetic energy E decrease over time, and the linear relation Ω ¼ E=ξ 2 a still holds (Fig. 3a). This results in a characteristic length scale ranging from ξ a; HUVEC % 30:7 μm to ξ a; MDCK % 44:1 μm, corresponding to a narrow length range of 1.5-2 cells. Besides, both the PDFs of the kinetic energy (Fig. 3b) and the enstrophy (Fig. 3c) in these different cell monolayer systems collapse to a common q-Gaussian (2) for c and Eq. (4) for d) with the Tsallis index q ¼ 1:2. The cell densities corresponding to time 0, 3, 6, 9, and 12 h are manually counted and estimated approximately as 1:4 10 5 cm À2 , 1:6 10 5 cm À2 , 1:9 10 5 cm À2 , 2:1 10 5 cm À2 , and 2:3 10 5 cm À2 , respectively. e The enstrophy Ω versus the kinetic energy E over time. Different symbols represent different data set (eight independent FOVs and four independent samples); the dash line refers to the fitting Ω ¼ E=ξ 2 a with the characteristic length scale ξ a % 44:1 μm and the correlation coefficient The energy spectra E(K) (f) and the enstrophy spectrum Ω(k) (g) at different time points. The arrows indicate the plateaus. Here, k c represents the wavenumber at the scale of one cell. distribution, which deviates remarkably from the classical Maxwell-Boltzmann distribution. Further, Fig. 3d (Supplementary Fig. 5a, respectively) shows the energy spectra (enstrophy spectra, respectively) in these different cell monolayer systems, which demonstrate similar landscapes and scaling behaviors. We linearly fit these curves at large wavenumbers in the log-log plots, and extract the scaling exponent β. For these different cell monolayers, β is bounded in a narrow range from 4.05 to 4.55 (Fig. 3d). These results reveal nearly common statistics of 2D mesoscale cell turbulence.
Insensitivity of the energetics of mesoscale cell turbulence to the substrate stiffness. Since the stiffness of the underlying substrate plays a significant role in regulating cell-matrix adhesion [41][42][43][44] and cell migration 16,45,46 , we next explore how it affects the energetic statistics of 2D mesoscale cell turbulence. We then culture MDCK cells on various substrates of different Young's moduli, including polyacrylamide (PA) gel ($ 10 kPa), polydimethylsiloxane (PDMS; $ 2 MPa), plastic ($ 1 GPa), and glass ($ 70 GPa). Mesoscale cell turbulence can be observed in all MDCK cell monolayers adhering to these substrates. Interestingly, the energetic statistics of cell turbulence, including the scaling behavior between kinetic energy and enstrophy, energy distribution, and energy spectrum, are all insensitive to substrate stiffness (Fig. 4). Although the substrate stiffness varies over six orders, the characteristic length scale ξ a just ranges from 38:5 μm to 44:1 μm (Fig. 4a and Supplementary Fig. 6), which remains 2 cell length. The PDFs of rescaled kinetic energy (Fig. 4b), as well as rescaled enstrophy (Fig. 4c), collapse to a common q-Gaussian distribution. Besides, the scaling exponent β at large wavenumbers also remains within a narrow range of 4.43-4.55 ( Fig. 4d and Supplementary Fig. 5b). These findings reveal that the statistics of cell turbulence is almost independent of substrate stiffness.
Numerical simulations. We use an active vertex model 47 to further probe the energetics of 2D mesoscale cell turbulence. In this model, a cell monolayer is characterized by a polygonal network, and cells are described by interconnected polygons   Supplementary Fig. 7a). To account for the active motility of cells, each cell is assigned with a self-propelled force 48,49 . The direction of the self-propelled force, i.e., cell polarity, can be affected by intercellular social interactions and random fluctuations. We consider two typical intercellular social interactions, including local alignment (LA) and contact inhibition of locomotion (CIL). Both LA and CIL are crucial in many morphogenetic processes, such as wound healing and collective chemotaxis 13,14,50 . LA describes the motion alignment effect among neighboring cells 51,52 , similar to the Vicsek interaction 53,54 . In contrast, CIL describes the repulsive interaction among neighboring cells upon contact, as present in many cell types 13,14,50 , resembling the repulsion effect that occurs in crowded animal groups 55,56 .
The collective cell dynamics is determined by the evolution of vertex positions, r i t ð Þ, dictated by where f U i ¼ À∂U=∂r i stands for the potential force acting on the vertex i, with U being the potential energy of the system; γ is the friction coefficient; v sp denotes the self-propelled velocity, and p J ¼ cos θ J ; sin θ J À Á represents the polarity vector with θ J being its direction; ε T is the intensity of translational noises and β T J t ð Þ are independent unit-variance Gaussian white noise vectors; n J refers to the number of cells that contact cell J, and P J2C i computes an summation over all neighboring cells C i of vertex i. The potential energy U of the system stems from cell area elasticity, cell contractility, and intercellular tension, and can be expressed as 47,57-60 where K a represents cell area stiffness; K c denotes cell contractile modulus; and Λ quantifies the interfacial tension between neighboring cells; A 0 refers to the preferred area and A J is the current area of the Jth cell; L J is the perimeter of the Jth cell; and l ij is the edge length of the cell-cell interface. Considering the LA and the CIL effects among cells, the polarity direction θ J evolves as 47 where μ LA and μ CIL quantify the intensities of LA and CIL, respectively. θ vel ð Þ K is the motion direction of cell K, and α J;K ¼ arg r J À r K À Á denotes the argument of the cell-cell relative direction pointing from cell K to cell J. C J represents the neighbor collection of cell J. ε R denotes the intensity of rotational noise, and η R J t ð Þ are independent unit-variance Gaussian white noises. In addition, T1 topological transition is involved to account for cell neighbor exchange (Supplementary Fig. 7b) Our simulations show that the theoretical model can capture the primary energetic characteristics of 2D mesoscale cell turbulence. First, the dynamic swirling pattern can be reproduced, as reflected by cells' trajectories (Fig. 5a), as well as the velocity field (Fig. 5b). Second, to mimic the jamming transition of cellular motions over time, we decrease the motility of cells (i.e., v sp ) gradually. Consequently, it predicts a linear relation between the kinetic energy E and the enstrophy Ω (Fig. 5c). Furthermore, our model suggests that the characteristic length scale ξ a can be regulated by cell-cell social interactions (Fig. 5c). For example, enhancing LA effect tends to enlarge the characteristic length scale ξ a (Fig. 5c). However, the regulation of cell-cell social interactions leads to a narrow-range variation (~2 cell length) of ξ a , in consistency with our experimental measurements for different cell monolayer systems (Fig. 3a). Our simulations also show that the PDFs of kinetic energy and enstrophy follow the q-Gaussian distribution rather than the Maxwell-Boltzmann distribution, especially when the CIL effect is strong; enhancing the LA effect could break the q-Gaussian statistics (Fig. 5d, e). As well, the active vertex model reproduces the power decaying law of the energy spectrum at large wavenumbers, with the scaling exponent close to that measured in experiments (Fig. 5f). These results suggest that cell-cell social interactions could play key roles in the emerging energetic statistics of collective cell flows.
Remarkably, in the limiting case when the LA effect is too strong, cells will self-organize into a solid-like flocking mode rather than a turbulent motion mode, leading to the energetic landscape strikingly different from our experimental measurements (Supplementary Fig. 8). Therefore, our simulation results point to that the energetic statistics of 2D mesoscale cell turbulence are primarily attributed to active cell motility and cell-cell social interactions, but its experimental validation and underpinning molecular mechanisms await further work.

Discussion
Mesoscale cell turbulence has attracted mounting attention due to its essential significance for many dynamical processes, such as multicellular organization 23 , nutrient mixing 24 , and thrombopoiesis 25 , but its energetic statistics remain incompletely understood. Using living cell imaging on confluent cell monolayer systems consisting of different cell types and various substrates, we probe the energetics of mesoscale cell turbulence emerging in 2D cell sheets. We find that such active turbulence possesses unique energetic features markedly different from classical turbulence in passive fluids. The PDFs of the kinetic energy and enstrophy of cell flows over time are found to obey the q-Gaussian distributions. Complementary to a recent finding that the distribution of cell geometry during epithelial jamming follows a k-gamma distribution 61 , our finding sheds light on cellular selforganization from a perspective of energetics.
Recent study on human bronchial epithelial cell monolayers showed that the enstrophy of collective cell motions is linearly related to the kinetic energy during monolayers' jamming transition 23 . Our results suggest that this linear relation is virtually common to vastly different cell lines, both epithelial and nonepithelial. Intriguingly, this linear constraint is found to conserve in coherent cells crawling on substrates with a wide range of stiffness, further emphasizing the generality of collective cell dynamics. It has been known that non-epithelial cells (C2C12 and NIH-3T3 cells) usually display relatively weak intercellular adhesions and substrate stiffness regulates cell dynamics through cell-substrate adhesion primarily 41,42,62,63 . Therefore, adhesive interactions between cells and between cells and environments may play a trivial role in the emerging common energetics across these epithelial and non-epithelial monolayers. However, our previous study showed that in the existence of free space, the probability distribution of cell velocities can be skewed: at regions close to the free space, cell velocities do not follow the q-Gaussian distribution; whereas away from the free space, the probability distribution of cell velocities recovers to the q-Gaussian 30 . This is the same case for energetic statistics, suggesting that the unique energetic statistical features of collective cell flows arise from the cell-cell interactions. Nevertheless, what's the specific cell-cell interaction that dictates the unique energetic statistical features of collective cell flows needs to be further explored.
In addition, in all cell lines, the kinetic energy decreases over time as maturation takes place. This is indicative that all cell lines become increasingly more jammed as collective cellular dynamics come to a halt. However, the jamming processes for different cell lines appear to take different dynamical "paths" since the kinetic energy-enstrophy relationship is characterized by a different parameter ξ a for different cell lines (Fig. 3a). This could be attributed to the difference of cell-cell social interactions, according to our simulations (Fig. 5c). This suggests that cell jamming 18,29,61,64 cannot be merely characterized by a dynamical slow down, but rather requires at least one additional dynamical parameter to characterize the degree of cooperativity. These results are in agreement with recent work showing that human airway epithelial cells can also undergo solid-to-fluid transitions via different dynamical and morphological paths 65 .
Our theoretical model and numerical simulations reveal that polarization-based cell-cell social interactions, including CIL and LA effects could contribute crucially to the energetic statistics of collective cell flows. However, the molecular mechanisms of LA and CIL have not been fully understood so far 13,52 . To grasp a whole picture of mesoscale cell turbulence, a multiscale model that involves molecular events and cellular behaviors simultaneously is highly deserved. Besides, our study never rules out the possibility that other intercellular interactions, e.g., contact following of locomotion 15 and contact enhancement of locomotion 66 , and even individual cellular properties, such as cell memory 49,67,68 , may regulate the energetic hallmarks of collective cell motions, which merits future investigation.
Taken together, this work could deepen our understanding of multicellular dynamics in diverse physiological and pathological processes, such as tissue morphogenesis, embryo development, and cancer metastasis from the viewpoint of energetics. The combined experimental, theoretical, and numerical results presented here may provide both qualitative and quantitative guidance for future efforts toward identifying the basic principles of dynamic evolution in living cell populations.
We culture cells on diverse substrates of different stiffness to explore the role of substrate stiffness in the energetic statistics of 2D mesoscale cell turbulence. Substrates including PA gels, PDMS, plastics, and glasses are used. In details, PA gel substrates are prepared according to the standard protocols reported previously 69,70 . We manufacture PA gels of stiffness~10 and 40 kPa, using the relative acrylamide (Sigma, 40% w/v) and bis-acrylamide (Sigma, 2% w/v) concentrations, as suggested in previous protocols 69,70 . The PA gel substrates are functionalized with collagen I (500 μg ml −1 ) before seeding cells for cell attachment. PDMS substrates (Sylgard 184; Dow Corning) are prepared by mixing the base and curing agent in a ratio of 1:10 (wt/wt), which results in a stiffness of 2 MPa (ref. 71 ), followed by manually stirring and degassing in a vacuum oven to remove air bubbles. A thin layer of PDMS is poured spin coated over a 35 mm Petri dish, degassed again, and cured at 80°C for 2 h. The Petri dish coated with a PDMS layer is exposed to UV for 30 min and then functionalized with fibronectins (Corning, 50 μg ml −1 ) before seeding cells for cell attachment. For the plastic substrate, cells are cultured on a 35 mm Petri dish (Corning) of stiffness~1 GPa. For the glass substrate, cells are cultured on a 35 mm confocal Petri dish (Corning) of stiffness~70 GPa.
Cells are seeded on a Petri dish (Falcon), or PA gel substrates functionalized with collagen I, or PDMS substrates functionalized with fibronectin, or confocal dish with glass bottom (Nest) for enough time (typically 1 day) to grow into a confluent monolayer of cell density~1:4 10 5 cm À2 . Subsequently, phase contrast images are acquired with a ×10 objective, using an Olympus IX83 inverted fluorescence microscope. Two successive images of the same field are taken at a time interval of 1 min. The FOV is of size 1.33 mm × 1.33 mm, much larger than the spatial correlation length (200-300 µm) of mesoscale cell turbulence 16,18,19 , thus large enough to access their statistical properties.
Particle image velocimetry analysis. In the PIV analysis, a confluent cell monolayer system is treated as a continuum fluid, and instantaneous cell velocities (average displacement within 1 min) at each grid point are given by cross-correlation calculation. In details, raw phase contrast images are firstly preprocessed by contrast enhancement followed by high-pass Gaussian filtering to conserve the high-frequency component and remove the background noise, with the kernal size of~1 cell (≈20 pixels in length). Velocity field is then computed from these preprocessed images using PIV. In the cross-correlation calculation of PIV, the interrogation window size is taken as 64 × 64 pixels, which roughly corresponded to a region of 3 × 3 cell length and is sufficiently large to contaiñ 10 cells, but still small enough to resolve spatial flow field structures on the order of ten cell length 40 . To obtain a more continuous velocity field for calculating the energy spectra and the velocity spatial correlation function, we set a small spacing (four pixels) between two consecutive interrogation windows. Outliers of the obtained velocity vectors are abolished and replaced by fitting values based on the neighboring velocity vectors. Further, we also correct the velocity field for small systematic drift effects by subtracting the mean velocity, as suggested by previous studies 16,18,40 . We have identified that such subtraction to eliminate the driftrelated biases does not affect the statistics of velocity field. Custom PIV software is written in MATLAB. We have verified that our PIV code possesses good accuracy (see Supplementary Note 1 and Supplementary Fig. 1).
Calculation of vorticity field. Based on the velocity field associated at regular grid points, the vorticity field, ω x ð Þ ¼ ∂v y =∂x À ∂v x =∂y, is numerically calculated using the standard central differentiation scheme. Specifically, the vorticity at the grid point (i, j) is calculated as where v i;j x (v i;j y ) and ω i,j are the velocity and vorticity at the grid point (i, j), respectively; Δx and Δy are the spacings between the neighboring velocity vectors in the x direction and y direction, respectively. For the simulation results, the velocity field associated at irregular grid points is firstly reconstructed to one associated at regular grid points (see Supplementary Note 4 and Supplementary  Fig. 9).
Calculation of energy spectrum. The energy spectrum E(k) is calculated based on Eq. (5), by virtue of the Wiener-Khinchine theorem. In details, we firstly calculate the spatial Fourier transformv k; t ð Þ of velocity field v r; t ð Þ; then the Fourier transform of the two-point velocity correlation function is calculated as S k ð Þ 1 v * k; t ð ÞÁv k; t ð Þ D E t , withv * k; t ð Þbeing the complex conjugate ofv k; t ð Þ; finally, we calculate the energy spectrum as E k ð Þ ¼ kS k ð Þ= 4πL 2 ð Þwith L 2 the area of FOV. For the simulation results, the velocity field associated at irregular grid points is firstly reconstructed to one associated at regular grid points (see Supplementary Note 4 and Supplementary Fig. 9).
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
All experimental data and any related experimental background information not mentioned in the text are available from the authors on reasonable request.