Interplay between cell height variations and planar pulsations in epithelial monolayers

Biological tissues change their shapes through collective interactions of cells. This coordination sets length and time scales for dynamics where precision is essential, in particular during morphogenetic events. However, how these scales emerge remains unclear. Here, we address this question using the pulsatile domains observed in confluent epithelial MDCK monolayers where cells exhibit synchronous contraction and extension cycles of ≈5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx 5$$\end{document} h duration and ≈200μm\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx 200~ \upmu \mathrm{m}$$\end{document} length scale. We report that the monolayer thickness changes gradually in space and time by more than twofold in order to counterbalance the contraction and extension of the incompressible cytoplasm. We recapitulate these pulsatile dynamics using a continuum model and show that incorporation of cell stiffness dependent height variations is critical both for generating temporal pulsations and establishing the domain size. We propose that this feedback between height and mechanics could be important in coordinating the length scales of tissue dynamics.


Introduction
Cells as individual units display a variety of morphological changes during sensing [1], migration [2] and division [3] phenomena. When several individual cells are placed in close proximity to each other, each cell exhibits a complex behavior that is often different from isolated conditions [4]. This complexity further increases when hundreds to thousands of cells are organized in a continuous tissue. The collective dynamics exhibited by tissue is distinct from specific characteristics of individual cells [5]. Understanding the evolution of characteristics while transitioning from one-to-many cells is a challenge. In order to characterize these emerging properties, a combination of experimental and theoretical approaches is important to show generic principles on simple systems. In this context, spontaneous oscillations are attractive phenomena because of their clear spatiotemporal characteristics [6].
Oscillations are critical for dynamic rearrangements during the development of several organisms [7][8][9][10]. Length scales of oscillations span from molecular assemblies of proteins to group of cells in tissues [11][12][13][14][15]. Similarly, the periods of oscillations scale from seconds to hours [16,17]. Several studies using experiments and theory have explored the biological pathways and the a e-mail: minamdar@iitb.ac.in (corresponding author) b e-mail: riveline@unistra.fr (corresponding author) underlying molecular actors [18][19][20]. Although different works propose various mechanisms that set the length and time scales associated with these pulsations, a thorough understanding of their origins is still missing.
Spontaneous planar pulsations of cells are shown to occur in Madin-Darby Canine Kidney (MDCK) monolayers [21][22][23][24][25][26]. We had previously shown [17] that pulsations appear randomly on regular culture without any coating. We also reported that these pulsations can be controlled by plating them on micro-patterns of controlled dimensions by matching pulsations sizes with grid size. Using this setup, now we show that MDCK cells also undergo height variations that inversely correlate with changes in area. We quantify the nature of these pulsations by using the power spectrum of the velocity divergence field and observe the appearance of characteristic length and time scales of ≈ 5 h and ≈ 200 µm.
Several theoretical models, both discrete and continuum, have been developed to get insights into different aspects of tissue pulsations and wave propagation [27][28][29]. The discrete descriptions include active vertex [17,26] and cellular phase field models [24], whereas the active continuum frameworks use underlying viscous or elastic rheology of the tissue to model the phenomena [23,28]. The non-equilibrium activity is implemented in these models using cell polarity, motility and active stress [30]. These active components are generally coupled to tissue kinematics and sometimes also with an additional chemical field such as myosin density [29]. In this paper, we developed a simplest version of a 2−D continuum model in which the tissue deformation kinematics depend on cell area and an active stress term. In such planar models, the tissue could be modeled to either be compressible or incompressible in terms of area or volume. In the current work, we model the tissue to be area compressible in 2D but incompressible in 3D along our observations of cell deformation. On the basis of this assumption of 3D volume conservation with our estimate, we also include a component to the planar tissue stress that we model to arise from spatial gradients in the thickness of the epithelial monolayer. We find that this contribution to the stress provides a part of the restoring force for pulsations and could play a critical role in setting the length scale for the pulsations.
The paper is organized as follows. In Sect. 2, Materials and methods, we detail the experimental and analysis techniques used in the paper. In Sect. 3, Results, we present the main experimental findings of this work regarding planar pulsations and height variations in the MDCK monolayer. Here, we also develop a continuum model for tissue pulsations and compare its findings with the experimental outcomes. Finally, in Sect. 4, Discussion, we provide an overview of our findings, more specifically the potential role of height variations in monolayer pulsations. We also briefly compare our approach with other related works on monolayer pulsations. Detailed derivations of the equations associated with the continuum model are provided in "Appendix."

Cell culture
MDCK-E-Cadherin-GFP cells (Nelson lab., Stanford) were cultured in low-glucose Dulbecco's Modified Eagle's Medium (DMEM) (31885-049, Invitrogen) supplemented with 10% FBS (10309433, HyClone) and 1% antibiotics (Penicillin-Streptomycin; 11548876, Invitrogen) at 5% CO 2 in 37 • C. The culture was maintained by replating cells before they reached confluence. For imaging, 1 million cells were added to 175 mm 2 area of coverslips (CS) mounted on imaging chambers and allowed to settle. After 1 h, the non-attached cells were washed and the media was changed to L-15 (11540556; Invitrogen) supplemented with 1% FBS. At this stage, the mean cell density was about 1 cell per 1000 µm 2 . In the next 12-16 h following seeding, cell proliferation leads to the formation of confluent monolayer without leaving any empty space. The pulsation phase appears as soon as the monolayer is confluent and lasts about 10 h. Later on as a consequence of cell proliferation and when the cell density is high enough to prevent largescale movements (corresponding to a mean cell density of about 6.5 cells per 1000 µm 2 ), the monolayer enters a jammed state consistent with [28,31].

Micro-patterning
The grid with the required dimension was fabricated as SU-8 photoresist (MicroChem) molds on silicon wafer using the standard photolithography procedure. Then, stamps of polydimethyl siloxane (PDMS) (Sylgard 184; Dow Corning) made with a pre-polymer to cross-linker ratio of 1 : 9 (V/V) were obtained by replica molding from the silicon wafer. The micropatterning was performed using standard soft lithography procedure as described in [17]. Briefly, coverslips were treated with "Piranha" (3:7 parts of H 2 O 2 (516813; Sigma) and 7% H 2 SO 4 (258105; Sigma)) for 10 min, followed by careful washing and sonication in double distilled water for 5 min. The coverslips were dried using N 2 blower before functionalizing with (3-mercaptopropyl)trimethoxysilane (S10475, Fluorochem) and stored in a dry and clean glass petri dish at 65 • C for 2 h. Meanwhile, a drop of 10 µg/ml of TRITC-labeled Fibronectin (FNR01; Cytoskeleton) was added to the plasma-treated PDMS stamps to allow the FN to settle on the stamp. After 1 h, the remaining liquid was removed and the stamp was carefully dried using N 2 blower. The stamp was then brought in contact with the activated coverslip and left untouched for 30 minutes in order to transfer the FN to the coverslip. Finally, the coverslip was incubated with 100 µg/ml solution of poly-L-Lysine-grafted-polyethylene glycol (pLL-g-PEG) (SuSoS) in 10 mM HEPES. After 20 minutes, the coverslip was rinsed three times with phosphate-buffered saline (PBS) (11530486; Invitrogen) and stored immersed in PBS in 4 • C for up to 5 h before seeding the cells.

Image acquisition
The fluorescent images for the height measurements were acquired using a Leica SP8 X line scanning confocal system equipped with Z galvanometric stage, controlled by LAS X interface and equipped with Oko lab stage top incubator. Using a 63x (1.4 NA, oil) or 40x (1.3 NA, oil) objective, a tile scan covering the size of a FN grid was performed with z steps of 0.5 µm. To get a larger field of view for computing the power spectrum, the phase contrast acquisitions were obtained using one of the microscopes with the corresponding objectives: inverted Olympus CKX41 (4x, 0.13 NA Phase, Olympus); SANYO MCOK-5M incubator microscope (10x, 0.25 NA Phase, Olympus). These microscopes are equipped with a CCD camera (Hamamatsu C4742-95/C8484-03G02; Sentech XGA) and operated by Hamamatsu Wasabi/Micromanager and MTR-4000 software, respectively. The interval between acquisitions were either 5 min, 10 min or 20 min, and all the imaging was performed at 37 • C for 48 h. To prevent evaporation of media during imaging, either mineral oil (M8410, Sigma-Aldrich) was added or the imaging chamber was covered with a glass Petri dish.

Image processing
The images were analyzed using FIJI application [32], and the power spectra were computed using MATLAB software. The schemes were drawn using Inkscape, and the plot showing the difference in height was obtained using GraphPad Prism 7. The height variations in the monolayer were obtained using a custom written FIJI plugin. The plugin allows visualization of changes in height as variations in intensities by encoding the height as intensity value to the corresponding pixel. In order to get the power spectrum of the pulsations, the velocity field was obtained by particle image velocimetry (PIV) using the MATLAB-based PIVlab application [33]. All experiments were repeated a minimum of three times, and the typical results are reported.

Results
In this section, we present the main findings of the current work on (i) localization of planar pulsations on patterned substrates, (ii) the correlation of height variation in the monolayer with planar contraction/expansion, (iii) quantification of the length and time scales associated with the monolayer pulsations and (iv) a continuum model that recapitulates experimental results.

Planar pulsations in MDCK monolayers are localized on a substrate patterned with fibronectin and PEG
We prepared MDCK monolayers with controlled density on coverslips with patterns of fibronectin (FN) (see Sect. 2 and Fig. 1). These adhesive patterns were grids of 120 µm thickness surrounding a square gap of 150 µm, passivated by pLL-g-PEG (PEG) that reduces adhesion of cells to the coverslips (Fig. 1). This setup generates a region of low friction connected to zones of high friction. As we previously reported [17], the FN grid led to spontaneous pulsations of MDCK cells spanning 10-15 cells in diameter. To gain insight into cell shape changes in 3D, we imaged the pulsating domain over 48 h using confocal microscopy (Fig. 2a). This experiment allowed us to gain insight into the potential correlation between the change in cell height and pulsations.

Pulsations in MDCK monolayers correlate with variations in cell height
During the pulsations, we found that the height of cells in the monolayer underwent fluctuations. This was confirmed by checking the appearance of apical side of cells at different axial planes over time (Fig. 2a). During the contraction phase, cells were visible already at z-plane 8 µm, while in the extended phase cells were visible only at the z-plane 2.5 µm (Fig. 2c). The maximum height reached by cells in the monolayer was about 10 µm (Fig. 2a). In addition, to quantify the height profiles, we encoded the cell height as intensity values ( Fig. 2b and 2c and (movie 1)). The results confirmed our visual observation, and the mean height profile of cells measured over a region of interest across time showed variations in mean cell height (Fig. 2b, c). Specifically, the height profiles between the contraction and extended phases were more than twofold different (Fig. 2d). This confirmed that pulsations are correlated with cell flattening with reduced height during the extension phase followed by cell squeezing with increase in height during the compression phase. The cell volume was estimated to be roughly constant, at around 400 µm 3 .

Planar pulsations in MDCK monolayers have characteristic time and length scales
As described earlier, the pulsations observed in the epithelial monolayer corresponded to the periodic formation of contracting and expanding domains within the monolayer. To quantify the pulsations, we first used particle image velocimetry (PIV) to obtain velocity field v(x, y, t) of the cells at grid points x i , y j at time frame t k . We then numerically obtained the planar divergence d(x, y, t) = ∇ · v since it is the relevant measure of contraction/expansion in the monolayer due to the velocity field. Representative divergence fields in space at a particular time and at a given space point as a function of time are shown in Fig. 3a, b, respectively. In order to extract the spatiotemporal characteristics of divergence, we numerically obtained its power spectrum S d (q x , q y , ω) through discrete Fourier transform, where q x and q y are the wavenumbers, respectively, along x and y, and ω is the frequency component of the power spectrum. To quantify the temporal behavior of d, we then calculated S d (ω) by summing S d (q x , q y , ω) over all q x and q y . The spatial characteristics of d are quantified via S d (q) by summing S d (q x , q y , ω) over all ω and also over q x and q y such that q 2 x + q 2 y = q 2 . Note that for convenience, we keep the same notation S d , irrespective of the associated argument. The power spectra S d (ω) and S d (q) are plotted in Fig. 3c, e, respectively. The blue line and the surrounding gray region correspond, respectively, to the average and the standard deviation of S d (ω) and S d (q) over three separate experimental repeats. It can be seen that S d (ω) has a peak at ω/2π ≈ 0.2 h −1 that corresponds to a pulsation period of ≈ 5 h. The peak corresponding to S d (q) is not as sharp and also shows some variation with respect to the experimental repeats. However, S d (q) peaks at q/2π ≈ 0.005 µm −1 which corresponds to a periodic pattern in divergence with a length scale of ≈ 200 µm.
Thus, we observe from the experiments that MDCK monolayers undergo planar pulsations with characteristic time and length scales. We also found that the planar pulsations were accompanied with modifications in the monolayer thickness that was both space and time dependent. Moreover, it was also shown in our earlier work that gradients in myosin concentration were involved in monolayer pulsations [17]. Hence, in order to get a better understanding of how planar deformations, thickness modifications and contractile mechanisms together contribute to the observed spatiotemporal pulsatory patterns in epithelial monolayers, we developed a continuum model based on a previous work [34].

Continuum model for spontaneous pulsations in MDCK monolayer
The experimentally observed pulsations in epithelial monolayers could be quantified in terms of the velocity divergence of the tissue flow that captures the dynamics of area expansion and contraction of tissue domains. It was also observed earlier [17] that myosin contractility was associated with these collective deformation modes of the monolayer. Hence, we now develop a simple continuum model with cell area a, contractile matter ζ and velocity field v as three variables that are functions of x, y and t. The mass conservation equation can be written as by neglecting the contributions from cell division and death. The stress σ ij within the tissue should satisfy the momentum balance equation that is given by Here, we assume that the imbalance in the internal tissue stress is countered by the passive fluid friction between the tissue and the substrate. Cell-substrate interaction could also produce active substrate traction, possibly arising from cell motility. However, for simplicity, in our model we do not include this potential contribution. In the above equation as well as the ones below, we use indicial notation for vectors and tensors with the indices i, j representing x, y coordinates. Since we are concerned only with the contractile and expansive movements of the tissue, the stress is modeled as an isotropic quantity Here, a 0 is the homeostatic area of the cells, K is the monolayer area modulus, K h is the gradient elasticity modulus that energetically penalizes the height variations in the monolayer (also see "Appendix C"), and ζ is the myosin dependent contractile active stress in the tissue. Finally, we discuss the dynamics of the contractile material c(x, y, t). For simplicity, we model the active stress ζ = Bc, where B is a proportionality constant. Consequently, to prevent the introduction of a new field variable c(x, y, t), we simply replace c with ζ without losing generality. Hence, the dynamics of the active stress ζ in the tissue is given by Note that the use of ζ instead of c would involve a simple linear scaling of the parameters β and β c with the factor B. Here, τ is the time scale for ζ to reach its homeostatic value ζ 0 , β > 0 is the coupling term that modulates ζ depending on deviation of cell area from a 0 , and ξ c is the chemical noise with strength β c . We then combine the above equations and linearize them for small perturbations δa = a − a 0 and δζ = ζ − ζ 0 . We then replace β βc δa a0 → a and β βc δζ ζ0 → ζ and non-dimensionalize length and time with = Kζ0 αβ andτ = ζ 0 /β, respectively. Finally, we get the following non-dimensionalized equations (see "Appendix A" for details) Here, the non-dimensionalization resulted from the following substitutions of parameters in terms of the original parameters.
Without the noise term ξ c , the above set of equations are satisfied for a = 0 and ζ = 0 and result in v i = 0. However, for nonzero noise, the system is perturbed and the feedback between a and ζ results in persistent fluctuations having dominant length and time scales that, as we show below, depend on the tissue parameters. Thus in our model, the noise plays an essential role in sustaining the pulsatory movements. By modeling ξ c to be correlated in time and uncorrelated in space (see "Appendix A"), and doing Fourier analysis of these set of equations, the power spectrum S d (q x , q y , ω) of velocity divergence can be obtained exactly. Further, and the denominator The noise term ξ c can have its origin, for example, in the modulation in myosin levels at cell-cell contacts that are known to cause junctional fluctuations [35,36]. The theoretical estimates of S d (ω) and S d (q) are presented in Fig. 3d, f, respectively, and show qualitative similarity with their experimental counterparts that include the presence of peak in S d (ω) and S d (q), thus providing characteristic time and length scales for the pulsatory cellular movements. Moreover, as is experimentally observed, both S d (q) and S d (ω) decay to zero for larger values of q and ω, respectively. However, unlike for the experimental data in which S d (ω = 0) = 0 (Fig. 3c), the theoretically estimated value is zero (Fig. 3d). By performing simple asymptotic analysis on S d , the dominant length (l 0 ), in terms of the original parameters, can be calculated as giving the approximate location of the peak for S d (q) in Fig. 3f (see Eqs. 7 and A.26). We obtain l 0 of about 100 µm, by taking typical values reported or estimated for the four parameters (see "Appendix D"). It can be seen that for a pulsating length scale to emerge in the monolayer, K h , the gradient modulus that is associated with spatial height variation is necessary as per our model. Similarly, the other important quantities for the emergence of pulsation length scale are the turnover time scale τ for myosin and the time scale τ c for the chemical noise. It can be seen that if any of these quantities tend to zero, then the peak of S d (q) is pushed toward higher values of q and eventually disappears. Moreover, when τ c → 0, i.e., the chemical noise is uncorrelated or white, the decay in S d (ω) and S d (q) is much slower than is experimentally observed. Thus, the qualitative behavior of the spatiotemporal patterns exhibited by the tissue in our model depends on a combination of tissue mechanical properties and characteristics of the chemical noise. We point out that the amplitude of the tissue response would be proportional to the noise amplitude β c . However, since we plot only the normalized values of power spectra, the noise amplitude does not explicitly show up in Fig. 3d, f. The details of this calculation and other estimates for the power spectrum including the experimental comparison are provided in "Appendix A." Based on these estimates and additional numerical exploration, we find that the qualitative behavior of S d (q) and S d (ω) (Fig. 3d, f) is similar to their experimental counterparts (Fig. 3c, e) for a wide range of parameters. We chose the typical parameters that showed this trend (Fig. 3). We note that although the experimental data (Fig. 3c, e) are obtained for the monolayer patterned substrate as shown in Fig. 1, the theoretical model does not account for the variations in substrate frictions that are expected to be inherent in the experiments. Hence, we also calculated S d (ω) and S d (q) from our previous data on non-patterned substrate [17]. We find that the differences between S d (ω) and S d (q), especially with respect to the location of maxima in q and ω, are not significant (see "Appendix B").

Discussion
We show that coordinated cell height variations with the planar contraction and extension cycles are potentially important in determining the size of pulsatile domains in monolayer tissues. Shape fluctuations in biological systems range from cell membrane [37] and organelles [38] to single [13] and collection of cells in embryos [15]. During developmental phenomena in which motions of cell collections can be synchronized, such movements underlie crucial morphodynamic events that exhibit outstanding spatial precision [39]. In this context, our experimental findings suggest that variations in cell height are associated with changes in cell area in a nearly incompressible viscoelastic cytoplasm and are consistent with one of our model assumptions. The coupling between the contractile machinery and the cell area, highlighted by β in Eq. 4, indicates that due to the chemical noise β c (Eq. 4), the concentration of the contractile material ζ can change on a cell-by-cell basis leading to varying levels of tension and therefore inhomogeneous height distributions across the pulsatile domains. The variations in height profiles of cells indicate that the height of cell-cell junctions is shortened and extended, respectively, during the extension and contraction phases. Therefore, the contractile stress of the cells combined with cell-cell junction height diminishes away from the pulsation center. In this specific context, the link between junction height and the contractile molecules is an interesting question to be explored, since the accumulation of myosin could be due to intercellular flows [40] or due to transcriptional activity [41]. Finally, it would be interesting to analyze the height variations with varying cell density at plating, since cell density was reported to have direct effects on cell velocities [28] and cell extrusion [31]. In our model, the relaxation of the contractile stress, dictated by the turnover time τ (Eq. 4), correlation time of chemical noise τ c , height gradient stiffness modulus K h (Eq. 3) and substrate friction α are crucial in dictating the inherent size l 0 of the pulsating domain of the monolayer (Eq. 9). Interestingly, although l 0 does not depend on the contractile activity ζ 0 in the tissue, the frequency of pulsation depends on a combination of the active stress and mechanical properties of the tissue (Eq. 10). All these quantities can be experimentally measured for the monolayer under the given experimental conditions. For example, the turnover time can be measured by fluorescence recovery after photobleaching (FRAP) of myosin in cells at different positions on the pulsating domain. The model can be further extended to explicitly include the role of signalling pathways that regulate contractility and deformation dynamics in tissues [36]. In addition, although we ignored the role of cell division and apoptosis on tissue area fluctuations, myosin accumulation that is involved in these process could be important in regulating the dynamics of pulsatile domains in the monolayer [42]. Indeed, not accounting for some of these important quantities and coupling could potentially explain the discrepancy between the experimental values and theoretical predictions ("Appendices A and B," Figs. 3,4,5).
Interestingly, other studies report the origin of pulsations to osmotic pressure from water flux where MDCK cells do not change significantly in height during pulsations [21]. Similarly, in single-cell zebrafish embryos, bleb formation which is primarily due to contractility is also regulated by water flow through aquaporins [43]. In vivo, during dorsal closure by amnioserosa cells in Drosophila, the force required for contraction has been shown to arise from both contractile actomyosin cable and cell volume change regulated by potassium channels [44]. Along this line, our study can be further extended to volume measurements (from apical-basal area and cell height) to determine the respective contributions of contractility and osmotic pressure in pulsations in MDCK monolayers. Such precise measurements can be used to refine the hydrodynamic descriptions for tissue dynamics with different sources of stresses and material flux [45]. Together, this will lead to detailed insights on the role of mechanics in tissue dynamics and organization at different scales.
Acknowledgements The MDCK-GFP-E-Cadherin cells were a kind gift from Nelson lab., Stanford. We are grateful to Guillaume Salbreux for discussions and many key insights into the theoretical model. We thank the Riveline lab. for discussions and comments and acknowledge the help of Basile Gurchenkov, Yves Lutz and Pascal Kessler from the IGBMC imaging facility. We thank the late Marcel Boeglin for writing the Fiji plugin to perform height measurements. R.T. was an IGBMC international

Author contribution statement
All authors contributed equally to this paper.

Appendix A: Derivation of equations for the continuum model and analytical estimates
Here, we systematically develop the continuum model described in the main paper and also provide detailed derivations of the relevant equations.
The theory has 2D cell area a, velocity field of the cells vi, internal stress σij and the contractile matter ζ as the field variables [34]. The relevant equations are In the above equations, α is the friction coefficient between the monolayer and the substrate, K is the linear area modulus of the cells, K h is the height modulus of the tissue that is modeled to depend on the area gradient in the tissue (see "Appendix C"). In Eq. A.4, τ is the time required for myosin to reach its homeostatic value ζ0, β is the coupling coefficient that relates cellular area with ζ modulation, and ξc is time-correlated noise satisfying ξc(x, y, t)ξc(x , y , t ) In the limit when the correlation time τc → 0, the noise becomes uncorrelated or white noise. The parameters in the continuum equations and their units in SI are summarized in Table 1-see also in "Appendix D." For small perturbations, δa = a − a0 and δζ = ζ − ζ0, eliminating the velocity term and linearizing in δa and δζ, we get (A.8) Here, we neglect the advection term of the form ζ0∂ k v k for conceptual simplicity in the weak advection limit. Choosing, τ = ζ 0 β as a reference time scale and = Kζ 0 βα as the length scale, we get the following non-dimensionalized dynamical equations. (A.10) We can further rewrite these equations to give Finally, we present the governing equations as where the non-dimensional dynamical variables are written as a ≡ β βc δa a 0 and ζ ≡ β βc δζ ζ 0 . The four, free, nondimensional parameters in these equations are μ ≡ μ in terms of the original dimensional parameters. The correlation time of the colored noise ξc is written as τc ≡ τcβ ζ 0 . Taking Fourier transform of the two coupled equations in space (x → qx and y → qy) and time (t → ω), we get whereã andζ are the space-time Fourier transforms, respectively, of a and ζ, and q 2 = q 2 x + q 2 y . Inverting the C matrix, we get The power spectrum fromã andζ then becomes (A. 18) We note that since q = q 2 x + q 2 y , the power spectrum only in terms of q can now be obtained as Sa(q, ω) = 2πq × Sa(qx, qy, ω). Since the linearized mass balance equation Eq. A.1 in terms of the non-dimensionalized variables and parameters is the power spectrum for velocity divergence becomes, S d (q, ω) = β 2 c β 2 ω 2 Sa(q, ω). Since we have taken ξc to be correlated noise (Eq. A.4), its power spectrum |ξc| 2 = 1/(1 + τ 2 c ω 2 ). Hence, the effective power spectrum for the velocity divergence is given as We can perform a scaling near q → ∞ .

(A.21)
Upon integrating the ω terms out, we get , (A.22) and near q → 0 we get ; (A.23) upon integrating the ω terms we get (A.24) As derived above, the series expansion for S d (q) is increasing and decreasing in q for small and large values, respectively, of q. Hence, we could expect the location qm of the maxima of S d (q) to approximately be at the intersection of these two curves. Hence, we equate them The wavelength associated with this pattern will be given by l0 = 2π/qm. Following the same procedure as for q, we now check for the asymptotic behavior for S d (q, ω) in terms of ω.
(A. 28) We find that the location of maxima for S d (ω) is obtained quite well when we equate S(q, ω → 0) and S(q, ω → ∞) corresponding q = qm as obtained in Eq. A. 26. From that, we get that the maxima for S d (ω) occur at ωm = 1 τ τc 5 12 [K − τ ](τ τc) Based on Eqs. A.21-A.29, we can get some insights into the role the non-dimensionalized parameters K, K h , μ, τ and τc on the behavior of power spectra S d (ω) and S d (q). From Eq. A.26, we find that, as per our model, the three terms K h , τ and τc are important for the observation of maxima in S d (q) as is observed experimentally (Fig. 3e). Hence, the term corresponding to gradients in area in tissue stress (Eqs. 3 and A.3), the active stress homeostatic term (Eqs. 6 and A.14) and the time-correlation component of the noise (Eq. A.6) play key role in setting the length scale of pulsations in the tissue within our model framework. Additionally, from Eqs. A.22 and A.24, we see that the rate at which S d (q) grows and decays for small and large values, respectively, of q is governed by the parameters τ, K, and K h . From Eq. A.28, it can be seen that although the parameter μ does not play an essential role in governing the qualitative nature of S d (q) and S d (ω), it modulates the decay behavior of S d (ω) for larger values of ω.
In Fig. 3, we plot the normalized power spectrum in terms of the wavenumber q and ω by defining the following.
Note that for notational simplicity, we use the same symbol S d for the different variants of the power spectrum.
In order to quantitatively obtain the low and high q and ω behavior of S d (q) and S d (ω), respectively, we show in Fig. 4 log-log plot of the data presented in Fig. 3a, c. Since, as opposed to the theoretical calculation, S d (ω = 0) = 0 for the experimental data, we fit a power law of the form ω m while neglecting a few data points at the beginning (Fig. 4a). The exponents m ≈ 1.4 and m ≈ −1.7 (Fig. 4a) are not the same as m ≈ 2 (Eq. A.27) and m ≈ −4 (Eq. A.28), respectively, as is predicted by the model. Similarly, the exponents m ≈ 1.5 and m ≈ −4.8 (Fig. 4b) for S d (q) also do not match with the theoretical estimate of m = 5 (Eq. A.24) and m = −3 (Eq. A.22).

Appendix B: Power spectrum for cellular pulsations on non-patterned surface
In our previous work [17], we had studied pulsations on epithelial monolayers. Here, similar to the power spectrum S d (ω) and S d (q) for the divergence of pulsatory cellular movements on patterned surfaces discussed in Sect. 3.3 and plotted in Fig. 3c, e, respectively, we do the equivalent for the non-patterned surfaces in Fig. 5a, c. As shown in Fig. 5, the blue curves and the gray-shaded region correspond to the mean and standard deviation of the quantities over five experimental repeats. The peak for S d (ω) happens at ω/2π ≈ 0.15 h −1 corresponding to oscillation period of approximately 6.5 h. Although the peak is not as well defined for S d (q), its maxima occur at around q/2π ≈ 0.004 µm −1 corresponding to a periodic pattern of approximately 250 µm. The values are not significantly different from that obtained for the patterned surface. We also check the values of the power spectra at larger and smaller values of ω (Fig. 5b) and q (Fig. 5d) and obtain the powerlaw exponents as shown in Fig. 4 for patterned substrate. Although the qualitative rise and decay in q and ω for the power spectra is similar for the patterned and non-patterned surfaces, the power-law exponents are not the same for these two conditions. However, the exponents for S d (q) do exhibit a better match with each other. As for the patterned surface, the exponents also do not match the values predicted by the model.

Fig. 5
Power spectra a S d (ω), c S d (q) and the corresponding log-log plots (b) and (d) for pulsatory cellular movements in MDCK epithelial monolayers on non-patterned surfaces. The nature of the power spectra is similar to that obtained for patterned surfaces in Fig. 3c, e. The blue lines and the gray-shaded region in (a) and (c) correspond to the mean and standard deviation, respectively, of S d (ω) and S d (q) obtained over five separate experimental repeats. The power-law fits for low and high values of (b) ω and (d) q are represented with the blue lines, and the corresponding exponents are indicated as written in Eq. 3. We note that in the absence of volume incompressibility of the monolayer, we cannot rigorously replace the ∇h term with ∇a as described above.
Although, in such a case, we do not expect major qualitative changes in our model except for re-scaling of K h , the use of ∇a term is, strictly speaking, a model approximation.