Raghavan Thiagarajan

Alka Bhat

Guillaume Salbreux

Mandar M Inamdar
email: minamdar@iitb.ac.in

Daniel Riveline
email: riveline@unistra.fr

Control and mechanisms of pulsatile flows in epithelial monolayers

Epithelial cells flows are observed both in vivo and in vitro and are essential for morphogenesis. Here, we show that pulsatile flows involving local contraction and expansion of a tissue can arise in vitro in an epithelial monolayer of Madine Darby Canine Kidney (MDCK) cells. The strength of pulsation can be modulated through friction heterogeneity by observing the monolayer dynamics on micro-contact printed fibronectin grids with dimensions matching the length-scale of spontaneous oscillations. We also report pulsations by inducing wound closure in domains of similar size with micro-fabricated pillars. In contrast, strongly coherent flows can be induced by adding and washing out acto-myosin cytoskeleton inhibitors. To gain insight into the associated cellular mechanisms, we fluorescently label actin and myosin. We find that lamellipodia align with the direction of the flow, and tissue-scale myosin gradients arise during pulsations in wound-healing experiments. Pulsations and flows are recapitulated in silico by a vertex model with cell motility and polarisation dynamics. The nature of collective movements depends on the interplay between velocity alignment and random diffusion of cell polarisation. When they are comparable, a significant pulsatile flow emerges, whereas the tissue undergoes long-range flows when alignment dominates. We conjecture that the interplay between lamellipodial motile activity and cell polarization, with a possible additional role for tissuescale myosin gradients, is at the origin of the pulsatile nature of the collective flow. Altogether, our study reveals that monolayer dynamics is dictated by simple rules of interaction at cellular levels which could be involved in morphogenesis.

I. Introduction

Tissue dynamics during development shape organs and show remarkable correlations in space and in time [START_REF] Adam C Martin | Pulsed contractions of an actin-myosin network drive apical constriction[END_REF][START_REF] Aigouy | Cell flow reorients the axis of planar polarity in the wing epithelium of drosophila[END_REF][3][START_REF] Donà | Directional tissue migration through a self-generated chemokine gradient[END_REF][START_REF] He | Apical constriction drives tissue-scale hydrodynamic flow to mediate cell elongation[END_REF]. These phenomena involve the active cytoskeleton which drives tissue movement and produce stress over large scales [START_REF] Adam C Martin | Pulsed contractions of an actin-myosin network drive apical constriction[END_REF]3,[START_REF] Rauzi | Planar polarized actomyosin contractile flows control epithelial junction remodelling[END_REF]. In addition, groups of hundred cells undergo flows in developing embryos, in Drosophila, zebrafish, and C. elegans. These dynamics have clear implications in their shapes, in their functions, and in their differentiations. In some cases, pulsatile flows have been reported, e.g., in the developing heart in zebrafish, where collective oscillatory motion is required for proper development [START_REF] Duchemin | Mechanotransduction in cardiovascular morphogenesis and tissue engineering[END_REF]. It is puzzling that cells connected to each other and to the extracellular matrix can lead to distinct types of dynamics. External cues breaking the symmetry through chemical gradients and/or local anisotropy in friction may contribute to determine different phases with the very same cellular material. Altogether, flows and pulsations reflect long range interactions across tens to hundreds of cells.

The coupling between neighbouring cells is expected to transduce the local cellular forces into a collective response of tissue. If a cell moves inside a confluent monolayer its velocity and its polarity will affect and will be affected by its neighboring cells. Transposed to a collection of cells it is expected that these local interactions should yield the global response. Some cellular parameters are already known to contribute to the long-range correlation such as persistent cell motion, polarity cues associated with lamellipodia, and anisotropy in myosin distribution. However, the explicit links between cell motility/polarity and the global response together with their molecular origins are still poorly understood.

An important challenge is to identify experimental setups and the relevant theoretical framework to test cellular mechanisms. Recently, many theoretical approaches have been designed to predict these dynamics [START_REF] Shaebani | Computational models for active matter[END_REF].

For example, the theory of active gels which considers living matter as an active continuum has successfully been implemented for understanding early stage of development in C. elegans [START_REF] Mayer | Anisotropies in cortical tension reveal the physical basis of polarizing cortical flows[END_REF] and large scales flows in tissues in Drosophila [START_REF] He | Apical constriction drives tissue-scale hydrodynamic flow to mediate cell elongation[END_REF]. Alternative models have also contributed to explicitly match cellular behavior and in silico cell dynamics in adapted versions of the so-called Vicsek model [START_REF] Sepúlveda | Collective cell motion in an epithelial sheet can be quantitatively described by a stochastic interacting particle model[END_REF]. Both types of approaches as well as alternatives complement each other [START_REF] Shaebani | Computational models for active matter[END_REF]. The current study aims at comparing dynamics of cell in vitro to control experimental conditions and numerical simulations to reveal mechanisms underlying collective effects in tissues.

In this study, we followed the dynamics of epithelial Madine Darby Canine Kidney (MDCK) monolayer in vitro on flat 2-D surfaces. We observe two characteristic features of the flow: an oscillatory component of the flow divergence, which corresponds to local pulses of contraction and expansion and a large correlation length of the velocity field. We aim to experimentally modulate these two features by modifying the friction between the tissue and the substrate in a patterned manner by printing fibronectin on the surfaces, or by exposing the monolayer to cytoskeletal drugs. To shed light on the underlying cell dynamics for this type of collective flow, we perform numerical simulations based on a vertex model to see how the features of the flow field can be recapitulated.

II. Characterizing the spontaneous flow in epithelial monolayers

As noted previously [START_REF] Thomas E Angelini | Cell migration driven by cooperative substrate deformation patterns[END_REF][START_REF] Mathieu Poujade | Collective migration of an epithelial monolayer in response to a model wound[END_REF], MDCK confluent monolayers undergo complex flow patterns changing rapidly in time (Fig. 1a,1c). Plotting the divergence and the velocity fields showed that the flow field is correlated on a characteristic length scale of ≈ 200µm (Fig. 1a). Focusing on regions of this size (Figs. 1b,1c, 1e, 1f), we noticed regions with visible pulses of contracting and expanding tissue (Fig. 1b,1c and Movie 1, where pulses appear to have a period of about 5 hours) or regions where the flow appeared to be more uniform and correlated over a large length scale (Figs. 1e, 1f and Movie 2). In pulsatile regions we also observed +1 topological defects in the flow field (Fig. 1b) and we typically observed 2 -4 pronounced pulses in the flow divergence (see Fig. 1c and the autocorrelation plot in the inset).

To quantify this behaviour, we calculated correlation functions in space and in time for the velocity field and its divergence (Fig. 1d,1g). We found that the correlation function of the divergence at the same spatial point but different time points, div(x, t)div(x, t + τ ) x,t as a function of τ exhibited an oscillatory behaviour with a period of about 5 h, indicative of periodic pulsatile flows occuring in the tissue. Calculating the spatial correlation function div(x + ∆x, t)div(x, t) x,t as function of ∆x, we found that the divergence is correlated on length scale of about 200 µm (Fig. 1d) which corresponds to around 15 cell lengths (see Section VIII Materials and Methods). Plotting the correlation function for the velocity field, we noted a spatial decay of ≈ 200 µm, consistent with the characteristic correlation length of the divergence. We also noted that the velocity field exhibited a temporal decay of a few hours, larger than the correlation time of the divergence. This is consistent with our observation of regions moving coherently with a relatively uniform velocity on long time scales (Fig. 1e,1f). Altogether divergence and velocity correlation plots allowed to highlight different collective features of cell motion in the monolayer.

III. Pulsatile flows can be modulated by spatially periodic tissue-substrate interaction

We next sought to test whether the period and magnitude of pulsatile flows could be influenced experimentally. Considering the cell monolayer as a layer of active material interacting with the underlying substrate, we reasoned that modulating the tissue-substrate interaction in a spatially periodic fashion would allow to influence the tissue-scale collective flow. For this, by using micro-contact printing (see Section VIII Materials and Methods), we deposited the extracellular matrix protein fibronectin on square motifs of controlled dimensions, called G -75 µm, G -150 µm, G -300 µm (see Figs. 2a-c, and Movies 3-5). We reasoned that fibronectin coating would locally increase the friction between the tissue and its substrate. We chose these values for the size of lower adhesion patches so that they were adjusted to the dimension of domains of spontaneous cell pulsations.

We then quantified the temporal and spatial correlation of the flow divergence, as we had done previously with uniform substrate coating (Figs. 2f, 2i, 2l). Remarkably, correlation plots of divergence indicated that the oscillatory component of the flow, as quantified by the dampening of the oscillation in the temporal correlation of the divergence, was maximal when the size of the motifs matched the spontaneous extension of pulsations, in the G -150µm case (Fig. 2g-i, Movie 4) In contrast, when the grid-size was smaller (G -75µm, Figs. 2d-f, Movie 3) or larger (G -300µm, Figs. 2j-l, Movie 5), pulsations were less pronounced. It is worth 2 noting that the magnitude of divergence of tissue flow is increased in the larger grid case (Figs. 2h,2k) compared to smaller grids (Fig. 2e), consistent with the notion of a lower overall friction between the tissue and its substrate when the fraction of substrate area covered by fibronectin is lower.

In these different conditions, more than ≈ 50% of the grids exhibited sustained oscillations. Oscillations were lasting for ≈ 48 hours, corresponding to a couple of cycles in the G -75 µm case, 6 cycles in the G -150 µm case, and a couple of cycles in the G -300 µm case Possibly, decay in oscillation amplitude was associated to cell proliferation, which could lead to reduced tissue motion. We also noticed that the number of ±1 defects was larger in the G -75 µm and G -300 µm cases than in the G -150 µm case (Figs. 2d, 2g, and 2j). Finally, in the 'resonant' case of G -150 µm (Movie 6), we aimed at decreasing friction in the central, non-adhesive region of the pattern by incubation with PLL-g-PEG, which passivates the surface. This process led to 'agitated' pulsations, involving localised rapid motions of cells within the passivated square, and this further supports the role of friction in the process. These results support the notion that friction can influence the nature of cellular movements. Altogether these results show that the magnitude of pulsatile flow can be controlled by spatially modulating tissue-substrate interaction.

IV. Generation of long-range flows in epithelial monolayers

We then wondered if the correlation length of the flow could be experimentally modified. Interestingly, we found that incubation for 5 hours with cytoskeleton inhibitors, followed by washing, could strongly modulate the correlation length of tissue-scale velocity (see Figs. 3a-g and Movies 7 and 8). Specifically, upon addition of blebbistatin (see Section VIII Materials and Methods), tissue motion was significantly reduced (Fig. 3a). However, subsequent wash-out of blebbistatin led to the appearance of flows with a longer spatial correlation and more persistence in time (Figs. 3a-d and Movie 7), as quantified by the velocity correlation functions.include A similar effect was also observed following addition and wash-out of latrunculin-A (Fig. 3f, 3g, and Movie 8).

Assuming that empty space could lead to pulsatile motion of cells, we generated a wound using laser ablation after wash-out of cytoskeletal inhibitor when the tissue exhibit a strongly coherent flow (see Fig. 3h and Section VIII Materials and Methods). Interestingly, cells then underwent back and forth motions with contraction and relaxation after ablation, in a way reminiscent of pulsatile flows observed without treatment (Fig. 3i). This was further confirmed by the measure of divergence correlation (Fig. 3j).

V. Cellular mechanisms behind pulsatile flows

We then sought to generate controlled wounds in the tissue to test their effect on the flow properties. Cells were plated on substrate with microfabricated pillars which were then removed, generating regularly placed empty spaces in the tissue. We then noticed that cells around the site of wound closure underwent pulsatile flows (Figs. 4a-b, and Movie 10). We took advantage of this set-up to visualize the cellular myosin distribution in cells participating to pulsatile flows (Fig. 4c). We noticed that contracting zones appeared to exhibit an increased density of myosin. We also saw the transient appearance of radial gradient of myosin fluorescence intensity around the wound, with larger intensity values towards the center of the wound in contracting phase. This gradient was more homogeneous in the expanding phase (Figs. 4c-d and Movie 10). These results suggest that myosin gradients could play a role in driving pulsatile flows in the tissue.

We also tested whether cell lamellipodia orientation was correlated with the tissue flow. To follow lamellipodia dynamics, we used mosaic experiments with cells transiently transfected with Lifeact-mCherry (Figs. 4e-f and Movie 11, also see Section VIII Materials and Methods). Using blebbistatin wash-out to induce largescale correlated flows, we found that lamellipodia direction indeed correlated with movement of cells (Fig. 4e). We chose to quantify the lamellipodium orientation by tracking and segmenting the outline of cells labelled with Lifeact-mCherry and calculate a nematic order parameter from the cell contour. We found that the cell shape nematic axis was then typically aligned with the direction of the flow (Fig. 4f-h). These observations support that lamellipodia extend along the direction of flows (Fig. 4f-h).

In summary, these observations support the idea that pulsatile contraction and expansion are associated with variations in cell myosin intensity and that cells are polarized along the direction of the flow.

VI. Numerical simulations to model tissue pulsatile flows

In order to identify the origins of the pulsatile collective flow, we ran numerical simulations of the tissue flow based on a vertex model. In our model, cells are assigned a polarity vector which changes over time according to a rule of local alignment and a polarity diffusion, and cells are subjected to a motile force along the polarity axis (Fig. 5a). Our model has seven parameters: preferred cell area (A 0 ), cell area modulus (K 0 ), line tension of cell interfaces (Λ), cell motility (v 0 ), polarity alignment rate (ξ), rotational diffusion of polarity (D r ), and substrate friction coefficient (η). We note that from these parameters a set of time scales can be defined (see Section VIIIG Materials and Methods).

We first tested the simulations with periodic boundary conditions and no rotational noise on the cell polarity vector, and we failed to generate pulsations. However, flows of cells were observed in many of the versions [START_REF] Marchetti | Hydrodynamics of soft active matter[END_REF]. This finding suggests that the generation of pulsations required another element, provided here by the rotational noise on cell polarity. We found that the absolute and the relative values of ξ and D r are important in setting the nature of collective migration modes in the tissue (see Section VIII Materials and Methods). High values of ξ promotes polarity alignment in some direction and hence long range flows. On the other hand, high values of D r promotes fast and short range oscillations. When the relative values of D r and ξ are well-balanced, we find temporally periodic and spatially oscillatory migration patterns of cells similar to experimental dynamics.

We next sought the reproduce three experimental conditions (Fig. 5 and Movies 12-15). We numerically mimicked standard monolayer, monolayer on grid and flows after blebbistatin wash-out. To perform these numerical experiments, we tested the parameters which could allow to reproduce experimental data. To be quantitative, we also plotted the divergence correlation and the velocity correlation as a function of time and space. The three situations are represented on Figs. 5b-j, and they all show good agreement with our data. We needed to have ξ Dr ∼ 1, grid dependent friction, and ξ Dr 1 to produce, respectively, Fig. 5b-d, Fig. 5e-g, and Fig. 5h-j.

Altogether these results reproduce the main observations of our experiments. We conclude that a simple polarity dynamics could account for the pulsatile nature of the flow. More specifically, modulating substrate friction, polarity diffusion and the strength of polarity-velocity alignment led to a variation in the apparent pulsatile behaviour of the tissue.

VII. Discussion

Here we report pulsatile flows occuring in epithelial monolayers, which can be modulated by acting on substrate friction, inducing wounds in the tissue healing, and addition and wash-out of cytoskeletal inhibitors. These perturbations affect both the external environment of the tissue and internal cellular properties. We also found that these dynamics are quantitatively reproduced with a minimal vertex model which encoded diffusion in polarity and polarity alignment as key ingredients.

These results go along other studies performed on MDCK epithelial monolayers with a variety of theoretical approaches [START_REF] Shaebani | Computational models for active matter[END_REF]. For example, various studies reported oscillations in vitro and they characterized them with various metrics yielding about the same length and time scales as ours [START_REF] Steven M Zehnder | A langevin model of physical forces in cell volume fluctuations[END_REF][START_REF] Steven M Zehnder | Multicellular density fluctuations in epithelial monolayers[END_REF][START_REF] Steven M Zehnder | Cell volume fluctuations in mdck monolayers[END_REF][START_REF] Peyret | Sustained oscillations of epithelial cell sheets[END_REF][START_REF] Petrolli | Confinement-induced transition between wavelike collective cell migration modes[END_REF][START_REF] Notbohm | Cellular contraction and polarization drive collective cellular motion[END_REF][START_REF] Rodríguez-Franco | Long-lived force patterns and deformation waves at repulsive epithelial boundaries[END_REF]. Their theoretical approaches were based on Langevin harmonic oscillator, continuum active elastic tissue, phase field active Brownian particles, and Voronoi active Brownian model. Here we have experimentally controlled and reproduced numerically using an active vertex model how external and internal tissue perturbation can lead to drastic changes in the properties of the collective tissue flow, and particularly on pulses of contraction and expansion.

Collective cellular motions are common in developmental biology [START_REF] Scarpa | Collective cell migration in development[END_REF]. However, the focus is generally on the genetic and signalling pathways associated with their onset. Still such dynamics are recognized to be essential in key morphogenetic events, for example, Drosophila trachea development and heart beating of zebrafish [START_REF] Duchemin | Mechanotransduction in cardiovascular morphogenesis and tissue engineering[END_REF][START_REF] Scarpa | Collective cell migration in development[END_REF]. Our in vitro results suggest that our readouts could be used in vivo to test mechanisms for these collective modes during morphogenesis. For example, pulsation could lead to local stresses associated with out-of-plane movements and with delamination, like in amnioserosa [START_REF] Moreno | Competition for space induces cell elimination through compaction-driven erk downregulation[END_REF][START_REF] Meghana | Integrin adhesion drives the emergent polarization of active cytoskeletal stresses to pattern cell delamination[END_REF]. Considering the large number of model systems available, ranging from Drosophila, C. elegans, to zebrafish and mouse, these self-organization rules could be explored experimentally and compared across organs and species.

This approach goes beyond the cellular scale and set the multi-cellular scale as the relevant length to consider to capture phenomena. Our combined study with control of tissue mode experimentally with close reproduction with numerical simulations led to these new readouts. They could not be anticipated prior cellular and numerical simulations and this shows that collective effects are satisfactorily addressed by integrated study with experiments and models.

These results go along the recent interest in the reported collective effects of epithelial monolayers using either topological defects or velocity correlations as readouts [START_REF] Beng | Topological defects in epithelia govern cell death and extrusion[END_REF][START_REF] Kawaguchi | Topological defects control collective dynamics in neural progenitor cell cultures[END_REF]. We add to these results the external and internal controls which could open the design of new experiments in vivo to determine the shape of specific organs in new directions and organisations. In turn, this may lead to new rules for setting cellular states, along the hypothesis that self-organisation would result from the interplay between cellular differentiation and mechanical processes [START_REF] Shahbazi | Self-organization of stem cells into embryos: a window on early mammalian development[END_REF].

VIII. Materials and Methods

A. Cell culture and cytoskeletal drugs MDCK labelled with Green Fluorescent Protein (GFP) for E-Cadherin (gift from Nelson lab, Stanford University) and MDCK labelled with GFP for myosin (gift from Shigenobu Yonemura lab, RIKEN) were cultured at 37 • C in low-glucose Dulbecco's modified Eagle's (DMEM) medium (31885-049, Inivtrogen) supplemented with 10% Fetal Bovine Serum (FBS) (10309433, HyClone) and 1% antibiotics (Penicillin-Streptomycin; 11548876, Invitrogen). Cells were maintained at a minimal confluency of 40% -50% and were trypsinized at 70% -80% confluency for experiments. Live experiments were performed in Leibovitz L-15 medium (11540556, Invitrogen) and unless mentioned, the initial seeding density was maintained at 10 6 cells 175 µm 2 for all experiments. In order to inhibit cytoskeletal proteins, these inhibitors were used with the following concentrations: blebbistatin (100 µM) (B0560, Sigma Aldrich), latrunculin A (1 µM) (L5163, Sigma Aldrich).

B. Micro-contact printing

Micropatterning of fibronectin on glass coverslips was performed using the micro-contact printing procedure [START_REF] Caballero | Protrusion fluctuations direct cell motion[END_REF]. First, motifs were microfabricated as SU-8 (MicroChem) molds on silicon wafer using standard photolithography technique. Then they were replica molded on PolyDiMethylSiloxane (PDMS) (Sylgard 184, Dow Corning; pre-polymer:crosslinker ratio -9 : 1 (V/V)) for micropatterning. During the micro-contact printing procedure, coverslips were hydrophilized with "Piranha" activation (Piranha: 3:7 parts of H 2 O 2 (516813, Sigma) and H 2 SO 4 (258105, Sigma)) followed by functionalization with (3-Mercaptopropyl)trimethoxysilane (S10475, Fluorochem) for 1 h by vapor phase deposition. Coverslips were then dried at 65 • C for ≈ 2 h. Meanwhile, PDMS stamps with microstructures were hydrophilized by oxygen plasma activation (Diener electronic) and incubated with TRITC or HiLyte 488-labelled Fibronectin (FNR01/FNR02, Cytoskeleton) (10 µg/ml) for 1 h. Then the dried stamp was brought in contact with the activated coverslip and a gentle pressure was applied to transfer the fibronectin patterns from PDMS stamps to the coverslip. The coverslip was then stored in Phosphate Buffered Saline (PBS) (11530486, Invitrogen) solution at 4 • C until the time of cell seeding. If the non printed zones were to be passivated, the micropatterned coverslip was brought in contact with pLL-g-PEG (poly-L-Lysine-grafted-PolyEthylene Glycol) (100 µg/ml) (SuSoS) in 10mM HEPES (pH 7.4) solution for 20 min. Then the coverslips were rinsed with PBS before proceeding for cell seeding.

C. Wound healing experiments

Wound healing experiments were performed using microfabricated pillars [START_REF] Anon | Cell crawling mediates collective cell migration to close undamaged epithelial gaps[END_REF]. The PDMS stamp with micropatterns was incubated with 0.2% pluronic acid (P2443, Sigma Aldrich) for 1 h and then rinsed with PBS. The stamp was then activated using an O 2 plasma cleaner and the side with micropillars was brought in contact with an activated glass coverslip for bonding. The bonded stamp and coverslip were placed in 65 • C for 30 minutes in order to strengthen the bonding. After UV sterilization of stamp+coverslip, cells at a concentration of approximately 10 6 /20 µl were added at the interface between stamp and coverslip. Cells filled the space between micropillars and they were allowed to settle for 45 min. Once the cells have settled, fresh DMEM media was added to the sample which was then stored at 37 • C for ≈ 12 h to allow cell spreading between micropillars. After ≈ 12 h, the PDMS stamp was carefully peeled off from the coverslip and the sample was finally rinsed and filled with L-15 media before starting image acquisition.

D. Optical setups and imaging conditions

Acquisition by phase contrast and fluorescence microscopy was performed using inverted Olympus CKX41 microscopes. The setups were equipped with a manual stage or a Marzhauser Wetzler stage with a stepper motor (MW Tango) enabling multi-point acquisitions, and an Xcite Metal-halide lamp or OSRAM mercury arc lamp for fluorescence acquisition. The phototoxicity during acquisition was prevented by using 5 two shutters (Uniblitz VCM-D1 / VMM-D1 or ThorLabs SC10) separately for white light and fluorescence lamp synchronized with a cooled charge-coupled device (CCD) from Hamamatsu C4742-95 / C8484-03G02 or Scion corporation CFW1612M. The devices were controlled by custom made scripts using either Hamamatsu Wasabi or Micromanager interfaces. In order to obtain a large field of view, we took objectives with magnifications of 4x (0.13 NA Phase, Olympus) and 10x (0.25 NA Phase, Olympus). For phase contrast imaging, we also used an incubator microscope -SANYO MCOK-5M, integrated with a 10x phase contrast objective, and a Sentech XGA color CCD to scan multiple samples illuminated by LED. The unit is controlled by an MTR-4000 software. For laser ablation experiments, we used a Leica TCS SP5 confocal system equipped with Z-galvo stage and multiphoton infrared femtosecond pulsed laser (Coherent Chameleon Ultra II). The SP5 setup is mounted with an inverted Leica DMI6000 microscope, PMT and HyD detectors that are controlled by LAS AF acquisition interface. When the monolayer starts to exhibit flows after blebbistatin wash-out, cells within an area of ≈ 150 µm × 150 µm were removed by ablation. The ablation was performed with smaller Region of Interests (ROI) within the ≈ 150 µm × 150 µm area to avoid damaging the neighboring cells. The ablation was performed with Leica 63x (1.4 NA) oil objective at 800 nm wavelength and ≈ 4300 mW laser power. All experiments were performed at 37 • C and images were acquired at 5 min, 10 min or 20 min interval between frames. In order to prevent evaporation of media due to long-term acquisitions (48 h -60 h), samples were covered with a glass Petri dish or by adding 2 ml of mineral oil (M8410, Sigma Aldrich).

E. Velocity Field and Correlation Functions

For all the analysis, the time point at which the monolayer reached complete confluency was assigned as time zero (t 0 ). In order to obtain velocity field v(x, y, t) in the monolayer, Particle Image Velocimetry (PIV) was performed using the open source software PIVlab [START_REF] Thielicke | PIVlab-towards user-friendly, affordable and accurate digital particle image velocimetry in MATLAB[END_REF]. Specifically, the velocity components (v x , v y ) at particular time t were obtained by using Direct Cross Correlation (DCC) algorithm on the monolayer images at times t and t + ∆t for interrogation window-size and step-size of 64 pixels and 32 pixels, respectively. The grid positions where no appropriate velocities could be obtained from the DCC algorithm were provided with velocity values via interpolation from the neigbouring grid points. Divergence (div = ∇ • v = ∂vx ∂x + ∂vy ∂y ) for the velocity field was calculated numerically using finite difference derivatives on the square grid plots. Thus for each of v x , v y and div, we obtain a N y × N x × N t size matrix, where N x and N y are the number of grid points in the x and y directions, respectively, and N t is the number of time frames. The spatio-temporal velocities and divergence thus obtained were further processed to calculate correlation functions in order to identify the underlying length and time scales in the following manner.

The velocity correlation function was obtained as:

C vv (x, y, t) = v(x , y , t ) • v(x + x, y + y, t + t) (x ,y ,t ) , (1)
where the averaging was done over the space-time grid points x , y , and t . The full correlation function was then azimuthally averaged to obtained C vv (r, t), where r = x 2 + y 2 . The spatial correlation function C vv (r) = C vv (r, 0), whereas the temporal correlation function C vv (t) = C vv (0, t). Please note that for convenience we use the same expression C vv for the different variants of the correlation function. Correlation functions for divergence were similarly obtained starting from the full expression

C dd = div(x , y , t ) × div(x + x, y + y, t + t) (x ,y ,t ) . (2)
Length and time scales for pulsations and flows were extracted with the following rule. If plots exhibited clear minima, we took this value. Alternatively, the intercept with the x-axis was extracted. We checked that the length and time values were consistent with actual motions measured on movies.

F. Winding Number

From visual observations of the monolayer time-lapse movies, we noticed that pulsations in the tissue had focal points that corresponded to the locations where the velocities were lower than the surrounding regions. We formally quantified the location of these singular points by calculating the winding number for the velocity field on the smallest closed cell of the underlying rectangular grid as follows. As shown in Supplementary Figure 1, there are four nodes i ∈ {1, 2, 3, 4} with velocity v i at the respective node. The change in the angle between the velocity vector at a given node and the next was estimated as

∆φ i = sgn(v i × v i+1 ) cos -1 (v i • v i+1 ) (3)
The expression sgn(v i × v i+1 ) determines if the sense of rotation from i to i + 1 is clockwise or anti-clockwise.

The second expression cos -1 (v i •v i+1 ) gives the smallest angle between i and i+1. The total winding number k over the contour is

k = 1 2π 4 i=1 ∆φ i (4)
where the loop is traversed in an anti-clockwise sense. If a singular point is present inside this loop then k takes an integer value -else it is zero. The singular points and their nature were further confirmed by plotting the streamlines of velocity fields. They broadly correspond to two main scenarios: (i) vorticity or source/sink (k = +1) and (ii) saddle point (k = -1). We found that the focal points of pulsation typically had a strength of k = +1 and corresponded either to a source ∇ • v > 0 or a sink ∇ • v < 0.

G. Computational Model

We use a cell-based vertex model to understand the various experimental observations. The model is as follows.

Basic Model

The basic vertex model for the pulsating colony consists of epithelial cells represented as polygons (Fig. 5a). Each cell polygon is constituted of vertices. Every pair of contiguous vertex is shared by an edge which in turn has a cell on each of its side. The position of the vertices and the connectivity of the cells define the geometry and topology of the vertex model. The mechanical properties are assigned to the minimal vertex model via a work function W :

W = N α=1 K α 2 (A α -A 0 α ) 2 + αβ Λ αβ l αβ . (5)
The first term of the work function models the resistance from a cell α with actual area A α to deviation from its preferred area A 0 α and K α is the corresponding area modulus of the cell. The second term represents the effective interaction energy between two connected cells α and β that results from a combination of cell-cell adhesion energy and acto-myosin contractility along the shared edge αβ. The force acting on a particular vertex i with position r i is given as

F basic = - ∂W ∂r i . (6)

Cell motility

It is known that cells within the colony get polarised and have a tendency to migrate. The motility behavior of an isolated polarised cell can be modeled using a simple description in which the cell moves with a speed v 0 along its polarity direction p. In a vertex model this self-propulsion tendency of cells is converted to an effective motile force on its constituent vertex i as

F motile = ηv 0 1 N i cellβ pβ , (7)
where summation is over the polarities p β , respectively, of the N i cells that contain the vertex i, and η is the effective friction coefficient between the cell vertex and the substrate [START_REF] Daniel M Sussman | cellgpu: Massively parallel simulations of dynamic vertex models[END_REF][START_REF] Comelles | Epithelial colonies in vitro elongate through collective effects[END_REF].

Dynamical Evolution

Initially, N = 3060 cells are created in a periodic box of size L x × L y . Each cell is provided a polarity p = cos φe x + cos φe y where the angle φ is chosen randomly from [0, 2π]. The time evolution of vertex positions r i is based on the equation

η dr i dt = F basic + F motile . (8)
To temporally evolve the polarity of the cells we make use of the experimental observations depicted in Supplementary Figure 2. Specifically, we assume that the polarity vector pα for a cell α has two tendencies:

7

(i) to align with the direction vα of its instantaneous velocity and (ii) to undergo rotational diffusion, which mathematically translates to the following differential equation [START_REF] Henkes | Active jamming: Self-propelled soft particles at high density[END_REF][START_REF] Ss Soumya | Coherent motion of monolayer sheets under confinement and its pathological implications[END_REF]

dp α dt = (ξm α align + ξ α rand )p ⊥ α , (9)
m α align = pα × vα • êz and êz = pα × p⊥ α . (10)
Here, ξ is the alignment strength of polarity and ξ rand is the diffusive, rotational, Gaussian noise with strength D r . The differential equations are solved using an explicit forward scheme with a time-step ∆t. Both the absolute and the relative values of ξ and D r are crucial in dictating the nature of collective migration modes in the tissue. High values of ξ promote polarity alignment in a particular direction and hence long range flows. On the other hand, high values of D r promote fast and short range oscillations. When the relative values of D r and ξ are well-balanced, we find temporally periodic and spatially oscillatory migration patterns of the cells that mimic the experimental dynamics. In the friction grid experiments, we mimicked the gridlike patterns that were created in the experiments with regions of high and low substrate friction. This was achieved in the simulations by increasing the value of η in the fibronectin patch relative to the value outside of the patch. While modifying the value of η, we ensured that the magnitude of the motile force ηv 0 remains the same both inside and outside of the fibronectin patch. Tissue dynamics in typical simulations is reported along with experimental movies (Movies 12-15.) The time-lapse images from the simulations have dimension of 1000 × 1000 px 2 and are collected approximately after every 50 simulation time-steps. As in the experiments, the images from simulations are analyzed using PIV to get the velocity field [START_REF] Thielicke | PIVlab-towards user-friendly, affordable and accurate digital particle image velocimetry in MATLAB[END_REF][START_REF] Moisy | PIVMat toolbox for Matlab[END_REF]. This velocity field is smoothened using Gaussian filter of size 2px using pivmat to reduce the noise arising from short-range spatial noise in the movements of individual cells [START_REF] Moisy | PIVMat toolbox for Matlab[END_REF]. The divergence of this smoothed velocity field is obtained numerically using finite difference method [START_REF] Moisy | PIVMat toolbox for Matlab[END_REF].

Parameter values for the vertex model

A detailed list of the parameters used in the system and their corresponding non-dimensionalisation with respect to some length (l 0 ), time (t 0 ) and force (f 0 ) scales is shown in the accompanying table. The force scale f 0 is arbitrary since the force parameter appears on both left and right hand sides of the dynamical equation and hence need not be chosen. For the simulations we choose the non-dimensionalised simulation time-step ∆t = 1. The most obvious choice of the length scale l 0 is related to the cell-dimension √ A 0 . Based on these parameters we can also extract additional length (L) and time (τ ) scales. 8

L ≡ Λ K 1/3 ; ηv 0 K 1/3 ; v 0 ξ ; v 0 D r (11) τ ≡ 1 ξ ; 1 D r ; v 0 √ A (12)

Comparison of simulation numbers with experimental values

There is an excellent qualitative correspondence between the experimental findings (Fig. 1, 2, 3) and the simulation results (Fig. 5). Below, we make quantitative connections between the experimental and simulation results by first linking the length and time-scales for plots of Fig. 5 with the experimental units.

The simulation velocity field of cells was quantified using PIV on images of 1000 × 1000 px 2 , sampled after ∆f ≈ 50∆t (simulation time-step). The distance R in the correlation function of Fig. 5d is reported in px units and time t in terms of the duration between two frames (∆t f ). Since there are N = 3060 cells in the 10 6 px 2 image, each cell length l cell corresponds to

l cell ≈ 10 6 3060 px = 18 px. ( 13
)
The minima of the correlation function for Fig. 5c is between 150 -200 px. Hence, the collective pulsatile movement occurs approximately on the scale of 8 -10 cell lengths. This number is similar to the experimentally observed collective length-scale of pulsation. Moreover, by noting the l cell ≈ 20 µm, we also estimate that

1 px = 20 18 µm ≈ 1 µm, (14)
thus giving correlation length of approximately 150 -200 µm which is consistent with the experimental findings. This also provides the experimental measure for distance R in Fig. 5. Going along these lines, we now make an estimate of the length scale (l 0 ) that is implicit in the simulations. From Table 1, the preferred area of a single cell A 0 = 1000l 2 0 . Taking A 0 ≈ 20 2 µm 2 , we obtain

l 0 ≈ 0.6 µm. (15)
The minima for the correlation function in Fig. 5c for the control case is at around 20∆t f (frame sampling). This corresponds to a period of ≈ 40∆t f ≈ 40 × 50∆t = 2000∆t, as the duration between two frames is around 50∆t. Since the experimentally observed pulsation time-period is ≈ 5 h, we estimate

∆t ≈ 5 2000 h = 1 400 h. (16)
This is also the same as the time-scale t 0 implicit in the simulations (Table 1). Also, the time difference between two frames of simulation time-lapse images is

∆t f ≈ 50∆t ≈ 0.125 h. (17)
This is the experimental measure for the time in Fig. 5. We check the consistency of this number as follows.

In simulation units, the value of cell motility v 0 = 0.1. From Table 1 we can note that

v 0 = v 0 t 0 l 0 , (18)
⇒ v 0 = 0.1 × 0.6 1/400 µm • h -1 ≈ 25 µm • h -1 , (19)
where we have used the estimates of l 0 and t 0 made above. This number is consistent with the experimental findings for cell speed.

H. Velocity-polarisation correlation for migrating cells

The quantification for the alignment of lamellipodia direction with the local flows is shown in Figure 4g and4h. This quantification is obtained as follows. We obtain the velocity field of the flow domain and cell contour separately and plot them together. Velocity field is obtained from the PIV of the domain using phase contrast images. From this, we compute the mean flow direction for the domain. For the lamellipodia direction, since cell direction is set by the leading front, we obtain the effective cell direction as a representative of the lamellipodia direction. We obtain the PIV of the cell contour which gives the vectors along the cell contour during cell movement. From this, we get the mean direction/angle of cell movement. This angle is used to represent the lamellipodia direction. We visually verified that the lamellipodia is always at the cell front in the direction of cell movement.

We also used another method to obtain the polarisation of a cell as defined in Supplementary Figure 2. Individual cells are tracked after blebbistatin wash-out over a duration of approximately 10 hours. First cells are outlined. We then obtain the centroid of the cell for every time-step t as follows:

r c = 1 N b r b , ( 20
)
where b is the index, respectively, corresponding to the N boundary pixels. The polarity of the cell is then defined in terms of its nematic tensor :

q xx = 2π 0 cos 2φdφ, (21)
q xy = 2π 0 sin 2φdφ, (22)
where the angle φ is defined between x axis and the line connecting the centroid and the boundary pixel as shown in Supplementary Figure 2a. The polarity can then be defined as p = cos θê x + sin θê y , where θ = 0.5 tan -1 (q xy /q xx ). The direction of cell velocity v is obtained from the instantanous displacement of the centroid between subsequent time steps. The correlation between the polarity and velocity direction thus defined is calculated by getting the frequency distribution over all time steps of the quantity

cos 2 ∆θ = (v • p) 2 . ( 23
)
The combined histogram for this quantity from multiple cases is obtained in Supplementary Figure 2. The distribution is skewed towards cos 2 ∆θ = 1, i.e., ∆θ, the angle between the velocity and polarity is systematically biased towards smaller angles (Supplementary Figure 2b. g) and (j) highlight the changes in nature of divergence field on changing the FN grid condition. The corresponding FN grid conditions are mentioned on the sides. In each of these panels, first row shows snapshots of one contraction -extension cycle of a pulsatile domain in a single grid unit where tracking highlights the cell trajectories during pulsations. Similarly, second and third rows show the winding numbers and divergence maps respectively for the snapshots in first row. Colour bar at the bottom provides the scale for divergence maps where blue and red indicate contraction and extension respectively. Scale bar 100 µm and time in hh:mm. See also the associated Movies (3), ( 4) and [START_REF] He | Apical constriction drives tissue-scale hydrodynamic flow to mediate cell elongation[END_REF]. (e), (h) and (k) show mean divergence plots that highlight the variations in divergence period over 48 h for different FN grid conditions. Arrows indicate data points corresponding to the domains shown in (d), (g) and (j) respectively. Insets show the autocorrelation. (f ), (i) and (l) represent the global characteristics. (f ), (i) and (l) show the spatial and temporal correlation plots for divergence to highlight the variation between different FN grid conditions. Each of these (f ), (i) and (l) plots correspond to the experiment from which the exemplary pulsatile domains (d), (g) and (j) were isolated. In each of these panels, first rows show snapshots of phase contrast images where the tracking highlights flow behaviour. Second rows show streamlines without the presence of topological defects. Scale bar 100 µm. Time is in hh:mm. See also the associated Movies [START_REF] Duchemin | Mechanotransduction in cardiovascular morphogenesis and tissue engineering[END_REF] and [START_REF] Shaebani | Computational models for active matter[END_REF]. (e)-(g) represent the global characteristics of monolayer and show the spatial and temporal correlation plots for velocity magnitude over 36 h. (h)-(j) Wound healing and the transformation of flows to pulsations; these panels describe the behaviour of monolayer after the closure of wound-induced by laser ablation, 12 h after blebbistatin wash-out. (h), (i) and (j) represent the local characteristics of monolayer over 24 h; (h) shows snapshot of phase contrast images of the ablated spot and the surrounding area. Wound due to ablation is highlighted by a dotted line in the first time point which eventually heals (second time point). In the consecutive time points, a smaller box is drawn to highlight a region exhibiting pulsatile behaviour. Scale bar 100 µm. Time is in hh:mm. See the associated Movie [START_REF] Mayer | Anisotropies in cortical tension reveal the physical basis of polarizing cortical flows[END_REF]. (i) describes pulsation in the region highlighted by a smaller box in (h). First row shows snapshots of phase contrast images where tracking highlights back and forth cell movement. Second and third rows show streamlines and winding number and divergence map respectively. Scale bar 50 µm. Time is in hh:mm; (j) shows the mean divergence plot over 24 h for the pulsatile domain shown in (i). Arrows represent the data points corresponding to snapshots in (i). Inset shows autocorrelation. [START_REF] Thomas E Angelini | Cell migration driven by cooperative substrate deformation patterns[END_REF]. (g) Plot showing the alignment of lamellipodial orientation (obtained from the effective cell direction -see Materials and Methods-H) and mean direction of flow field, corresponding to snapshots in (f ). Direction is given by the distribution of angles corresponding to the 2D space of images shown in (f ). Velocity magnitude is given by the concentric circles. Green arrow corresponds to the mean direction of flow field (shown in green in (f )) and red arrow corresponds to the mean orientation angle of lamellipodia (shown by the red cell in (f )). (h) Plot with all the time points during the time of flow lifetime (12 h), showing the alignment of lamellipodia direction (obtained from the effective cell direction -see Materials and Methods-H) in x axis with the mean flow direction in the domain (y axis), for the cell shown in 4e-f. The linear tendency confirms the alignment of lamellipodia with the flows. Scale bar 50 µm and time is in hh:mm.

IX. Figures and movie captions

Following sequence is followed for presentation below: 1 .

 1 Main figure captions 2. Supplementary figure captions 3. Movie captions 4. Main figures 5. Supplementary figures 13

Figure 1 .

 1 Figure 1. Spontaneous pulsations and flows in epithelial monolayer. (a) Snapshot of a confluent MDCK monolayer in phase contrast and the corresponding divergence and velocity magnitude colour maps. Scale bar, 200 µm. (b), (c), (e) and (f ) represent the local characteristics of monolayer. (b) shows snapshots of contraction and extension phases of one period of a pulsating domain and (e) shows snapshots of a flow domain. In and (e): First row corresponds to phase contrast images and tracking highlights the back and forth movement of cells during pulsation in (b) and flow trajectory in (e); Second row shows the streamlines obtained while computing winding number for each snapshot corresponding to first row. In (b), the position of a topological defect with winding number +1 indicates the position of pulsation centre and in (e), the absence of winding number confirms the flow behaviour; Third and fourth rows show the corresponding divergence and velocity magnitude plots. Scale bar 50 µm. Time is in hh:mm. In all panels where colour bars indicate the nature of divergence field, blue indicates contraction and red indicates extension; and similarly, where colour bar indicates velocity magnitude, red and blue indicate low and high velocities respectively. (c) shows the plot of total divergence for the pulsatile domain over 48 h represented in (b) and in Movie (1) and (f ) shows the velocity magnitude plot of flow domain shown in (e) and in Movie (2). Arrows in (c) show the data points corresponding to snapshots in (b). Insets show the autocorrelation of divergence in (c) and the autocorrelation of velocity magnitude in (f ). (d) and (g) represent the global characteristics of monolayer and respectively show the divergence and velocity correlations for distance and time for the same experiment over 48 h.

Figure 2 .

 2 Figure 2. Substrate friction determines the intrinsic length scale of pulsations. (a), (b) and (c) show different fibronectin (FN) grid conditions. Each grid unit -highlighted in dotted lines in (a), (b) and (c) has a non-FN area or a gap (G) shown in black, and a surrounding FN region shown in red. In all conditions, the FN width (W) is maintained constant at 120 µm while the gap (G) is varied: 75 µm in (a); 150 µm in (b); and 300 µm in (c). Scale bar, 200 µm. (d)-(e), (g)-(h) and (j)-(k) represent the local characteristics. (d), (g) and (j) highlight the changes in nature of divergence field on changing the FN grid condition. The corresponding FN grid conditions are mentioned on the sides. In each of these panels, first row shows snapshots of one contraction -extension cycle of a pulsatile domain in a single grid unit where tracking highlights the cell trajectories during pulsations. Similarly, second and third rows show the winding numbers and divergence maps respectively for the snapshots in first row. Colour bar at the bottom provides the scale for divergence maps where blue and red indicate contraction and extension respectively. Scale bar 100 µm and time in hh:mm. See also the associated Movies (3), (4) and (5). (e), (h) and (k) show mean divergence plots that highlight the variations in divergence period over 48 h for different FN grid conditions. Arrows indicate data points corresponding to the domains shown in (d), (g) and (j) respectively. Insets show the autocorrelation. (f ), (i) and (l) represent the global characteristics. (f ), (i) and (l) show the spatial and temporal correlation plots for divergence to highlight the variation between different FN grid conditions. Each of these (f ), (i) and (l) plots correspond to the experiment from which the exemplary pulsatile domains (d), (g) and (j) were isolated.

14

 14

Figure 3 .

 3 Figure 3. Transition from pulsations to flows and back to pulsations. (a)-(g) Resetting myosin activity leads to transition from pulsations to flows. (a), (b) and (c) show snapshots of phase contrast images overlapped with velocity vectors and below are the corresponding velocity magnitude maps. The orientation of velocity vectors indicates the direction of cell movement and the length indicates magnitude. (a) shows a snapshot before addition of where velocity vectors exhibit a pulsatile pattern. (b) shows a snapshot in the presence of blebbistatin where vectors are very small, (c) shows a snapshot after the wash-out of blebbistatin where the velocity vectors show similar orientation over large lengths. Scale bar, 200 µm. (d)-(e) and (f )-(g) report the flow behaviour after washing out blebbistatin and latrunculin-A respectively, in two different experiments. (d) and (f ) represent the local characteristics of monolayer.In each of these panels, first rows show snapshots of phase contrast images where the tracking highlights flow behaviour. Second rows show streamlines without the presence of topological defects. Scale bar 100 µm. Time is in hh:mm. See also the associated Movies[START_REF] Duchemin | Mechanotransduction in cardiovascular morphogenesis and tissue engineering[END_REF] and[START_REF] Shaebani | Computational models for active matter[END_REF]. (e)-(g) represent the global characteristics of monolayer and show the spatial and temporal correlation plots for velocity magnitude over 36 h. (h)-(j) Wound healing and the transformation of flows to pulsations; these panels describe the behaviour of monolayer after the closure of wound-induced by laser ablation, 12 h after blebbistatin wash-out. (h), (i) and (j) represent the local characteristics of monolayer over 24 h; (h) shows snapshot of phase contrast images of the ablated spot and the surrounding area. Wound due to ablation is highlighted by a dotted line in the first time point which eventually heals (second time point). In the consecutive time points, a smaller box is drawn to highlight a region exhibiting pulsatile behaviour. Scale bar 100 µm. Time is in hh:mm. See the associated Movie[START_REF] Mayer | Anisotropies in cortical tension reveal the physical basis of polarizing cortical flows[END_REF]. (i) describes pulsation in the region highlighted by a smaller box in (h). First row shows snapshots of phase contrast images where tracking highlights back and forth cell movement. Second and third rows show streamlines and winding number and divergence map respectively. Scale bar 50 µm. Time is in hh:mm; (j) shows the mean divergence plot over 24 h for the pulsatile domain shown in (i). Arrows represent the data points corresponding to snapshots in (i). Inset shows autocorrelation.

Figure 4 .

 4 Figure 4. Molecular actors associated with cellular movements. (a)-(d) Fluctuation in myosin density during contraction and extension phases of pulsation. (a) Schematic of the micropillar assay. Monolayer was grown in the presence of PDMS micropillars which on removal left free spaces/wounds for the monolayer to close. After wound closure, these sites exhibited pulsations. (b) Wound closure leading to pulsations. First and second rows show snapshots of phase contrast images and the associated divergence maps. In both rows, first and second time points represent the wound and its eventual closure respectively. Last four time points represent the pulsation that follows wound closure. (c) Fluctuation in myosin intensity. First row shows snapshots of myosin-GFP corresponding to first row in (b). Densification and diffusion of myosin can be observed at times 12:40 and 15:10 respectively. This densification and diffusion of myosin correspond to the snapshots showing contraction (blue region) and extension (red region) phases of pulsations in divergence map in (b). Second row shows myosin intensity profiles obtained along the yellow stripe shown in first row. Myosin intensity is initially zero at the centre (0 µm) of the wound (indicated by arrows) and increases to half-way mark during wound closure. But as the pulsation begins with contraction (at 12:40 in first row), myosin intensity reaches maximum and diffuses back to half-way mark as the domain extends. See the associated Movie (10). (d) and intensity plots for the pulsatile domain shown as snapshots in (b) and (c) respectively. Green region represents the wound closure phase and pink region represents the pulsation phase. (b) and (c), scale bar 100 µm and time is in hh:mm. (e)-(g) Orientation of lamellipodia aligns with the direction of flow. First two time points in (e), (f ) and (g) represent the phase before flow initiation; the transition is shown by the third time point; and the last three time points represent the flow phase. (e) First row shows the phase contrast image of MDCK-E-cadherin-GFP cells where a cell transiently transfected with LifeAct-mCherry is shown in red to demonstrate the orientation of lamellipodia before and after initiation of flow. Second row shows the streamlines and winding number. The presence and absence of winding number before (first two snapshots) and after (last three snapshots) initiation of flow, demonstrates the transition to flow behaviour. (f ) Superimposition of lamellipodial orientation in the cell (red) and the direction of flow field vectors (green). The cell shown in red corresponds to the cell shown in first row of (e). See the associated Movie[START_REF] Thomas E Angelini | Cell migration driven by cooperative substrate deformation patterns[END_REF]. (g) Plot showing the alignment of lamellipodial orientation (obtained from the effective cell direction -see Materials and Methods-H) and mean direction of flow field, corresponding to snapshots in (f ). Direction is given by the distribution of angles corresponding to the 2D space of images shown in (f ). Velocity magnitude is given by the concentric circles. Green arrow corresponds to the mean direction of flow field (shown in green in (f )) and red arrow corresponds to the mean orientation angle of lamellipodia (shown by the red cell in (f )). (h) Plot with all the time points during the time of flow lifetime (12 h), showing the alignment of lamellipodia direction (obtained from the effective cell direction -see Materials and Methods-H) in x axis with the mean flow direction in the domain (y axis), for the cell shown in 4e-f. The linear tendency confirms the alignment of lamellipodia with the flows. Scale bar 50 µm and time is in hh:mm.

Figure 5 . 16 Movie 13 .

 51613 Figure 5. Simulation results from active vertex model. (a) Basic schematic of the vertex model and the source of forces on the vertices. The outcome of simulations for (b-d) control (Movie 12), (e-g) friction grid (Movie 13), and (h-j) flows (Movie 15) corresponding to their experimental counterparts in Figures 1, 2, and 3, respectively. Divergence maps are shown within a selected domain from the simulations at different for (b) control and (e) friction grid conditions. Periodic contraction and expansion of a group of cells are visible along with their back and forth movement. Total divergence for local domains in (b) and (e) as a function of time is shown, respectively, in (c) and (f ). The time periodicity is apparent from the corresponding time autocorrelation function shown in the insets of the individual plots. The time and space correlations for divergence over the entire region corresponding to (b)and (e) are shown in (d) and (g), respectively. As was observed in the experimental findings, a length scale emerges both for (d) control and (g) friction grid. The existence of minima in the divergence time correlation functions in (d) and (g) indicates the presence of periodic pulsations. (h) In the simulations corresponding to wash-out experiments, the presence of long-range flows is seen in the selected region. (i) The corresponding velocity magnitude within the domain is nearly constant in time. (j) Velocity correlation function over the entire simulation domain shows persistence (long correlation) in both space and time. For a detailed discussion on experimental units of time and length, see Section VIIIG5.

Acknowledgments

We thank F. Jülicher and S. Soumya for helpful discussions. We thank the Riveline Lab for discussions and help, and we acknowledge the Imaging Platform from IGBMC. We are grateful to W. James Nelson, Shigenobu Yonemura and Sylvie Coscoy labs for sending cell lines. D.R. acknowledges support from CNRS (ATIP), ciFRC Strasbourg, the University of Strasbourg, Labex IGBMC. This study with the reference ANR-10-LABX-0030-INRT has been also supported by a French state fund through the Agence Nationale de la Recherche under the frame programme Investissements d'Avenir labelled ANR-10-IDEX-0002-02. R.T. was an IGBMC International PhD Programme fellow supported by LabEx INRT funds. M.M.I. acknowledges funding from Department of Biotechnology, Government of India (BT/06/IYBA/2012) and Industrial Research and Consultancy Center (IRCC), IIT Bombay. G.S. was supported by the Francis Crick Institute which receives its core funding from Cancer Research UK (FC001317), the UK Medical Research Council (FC001317), and the Wellcome Trust (FC001317). D.R. and M.M.I. also acknowledge hospitality at MPI-PKS, Dresden, where a part of this work was done.

Supplementary Figure 1. Calculation of winding number for the velocity field over the smallest cell of the PIV grid as described in Section VIII F. Supplementary Figure 2.