Pulsations and flows in tissues as two collective dynamics with simple cellular rules.

Collective motions of epithelial cells are essential for morphogenesis. Tissues elongate, contract, flow, and oscillate, thus sculpting embryos. These tissue level dynamics are known, but the physical mechanisms at the cellular level are unclear. Here, we demonstrate that a single epithelial monolayer of MDCK cells can exhibit two types of local tissue kinematics, pulsations and long range coherent flows, characterized by using quantitative live imaging. We report that these motions can be controlled with internal and external cues such as specific inhibitors and substrate friction modulation. We demonstrate the associated mechanisms with a unified vertex model. When cell velocity alignment and random diffusion of cell polarization are comparable, a pulsatile flow emerges whereas tissue undergoes long-range flows when velocity alignment dominates which is consistent with cytoskeletal dynamics measurements. We propose that environmental friction, acto-myosin distributions, and cell polarization kinetics are important in regulating dynamics of tissue morphogenesis.

Introduction

During development, tissue dynamics shapes organs with remarkable choreography in space and time (Rauzi et al., 2010; Aigouy et al., 2010; Roh-Johnson et al., 2012; Donà et al., 2013; He et al., 2014). Groups of hundreds/thousands of cells flow like liquids in developing embryos in Drosophila (Etournay et al., 2015), zebrafish (Behrndt et al., 2012), and insects (Münster et al., 2019). Solid-like pulsatile movements are also reported in Caenorhabditis elegans (Michaux et al., 2018); in the developing heart of zebrafish, where collective oscillatory motion is required for its proper development (Boselli et al., 2017); and, more generally, in muscle layers spanning epithelial layers (Zhang and Labouesse, 2012). These millimeter-scale changes in shapes involve the acto-myosin cytoskeleton, which drives tissue movement and generates stress at the organism level (Martin et al., 2009; Roh-Johnson et al., 2012; Rauzi et al., 2010). These motions and deformations drive the shaping of tissues, eventually establishing organ functions through cell differentiation. Altogether, flows and pulsations are general phenomena in morphogenesis.

Strikingly, these changes in shapes are reiterated in evolution with certain ‘families’ of tissue transformation, such as elongation and gastrulation, among many morphogenetic events. Even if their genetic backgrounds are different, the basic cellular mechanisms are expected to be shared among different biological model systems (Trepat and Sahai, 2018). Indeed, simple physical rules to achieve morphogenesis are known: (1) tissues are composed of cells that change their shapes over time, (2) cells in tissues interact with their neighbors by adhering to each other via adherens junctions and to the surrounding extracellular matrix through focal contacts (Geiger et al., 1995). It is puzzling, therefore, that the same monolayer of cells connected to each other and the matrix can undergo distinct dynamics. Regardless, the subtle complexity of cellular interactions that remain to be elucidated result in self-organization at the tissue level and minor perturbations in the program can lead to very different tissue morphogenesis reported in vivo (Graner and Riveline, 2017). Indeed for decades the existence of cellular flows and pulsations with groups of associated cells has been recognized in developmental biology across model systems.

Various biological systems share a number of dynamical morphogenetic events in vivo, but they are rarely investigated with a unified physical model to observe their diversity. For instance, during blood vessel formation, seven events were reported to contribute to angiogenesis, i.e., cell migration, anastomosis, apical membrane fusion, cell elongation and rearrangements, cell splitting, cell self-fusion, and cell division (Geudens and Gerhardt, 2011; Betz et al., 2016). Each of them may contribute to the formation of the vessel and the associated cellular flows. More generally, this so-called mechanism of ‘branch formation’ is essential to shape organs, and goes beyond the formation of the vascular systems (Affolter et al., 2009; Ochoa-Espinosa and Affolter, 2012). In fact, during branch formation, the lungs, kidneys, and mammary glands undergo these different dynamic events at the single-cell level, leading to global cellular reorganizations. Cells at the ‘tip’ of the branches are mesenchymal (behaving as migrating single cells) and they are connected to the ‘stalk’ made of epithelial cells (behaving as a connected tissue), which may undergo passive motion and/or rearrangements through neighbor exchanges and cell proliferation for tracheal branching and vertebrate angiogenesis (Affolter et al., 2009; Ochoa-Espinosa and Affolter, 2012). However, the identification of conserved principles still lacks a proper understanding of physical rules for cell level interactions that lead to global movements.

Several types of single-cell dynamics in morphogenesis, such as the growth of protrusions through filopodia and lamellipodia, binding to the surrounding extracellular matrix with focal contacts, and connecting with the neighboring cells through adherens junctions, have been reported (Baum and Georgiou, 2011). Cells can also have motility either as independent individuals or collectively in a group. In the latter case, it is often difficult to disentangle the motion of cells as arising from, for example, either active single-cell motion or from cell proliferation because of an increase in cell density and the associated reorganization of tissues with gradients in cell density stress (Tlili et al., 2018). Thus, in the search for quantitative methods to characterize collective motions beyond experiments with live microscopy, new approaches from theoretical physics for active matter are needed (Shaebani et al., 2020). Indeed, tight and quantitative comparisons with numerical simulations are known to be useful in testing these collective dynamics of groups of cells emerging out of interactiion rules at the cellular levels.

Another level of complexity is related to the cell state, which may change during development, for example, through differentiation and mesenchymal-to-epithelial transitions (Irianto et al., 2016). Changes in cell states can be either chemical or mechanical or both during the subsequent steps following exit from pluripotency. Biological changes in cells are also accompanied by mechanical modifications and alterations in contractility and intercellular adhesion because of changes in protein expression (Irianto et al., 2016; Kalluri and Weinberg, 2009). However, if the chemical state alters the mechanical state of the cells at any time, physical principles governing these cellular dynamics are robust and still at play irrespective of the cell states. Thus, the search for physical principles is essential to decipher the rules of tissue morphogenesis during development. Hence, unifying experiments tightly coupled to models are expected to unravel new means by which cells organize themselves in tissue.

Toward this aim, we use the dynamics of epithelial Madin-Darby Canine Kidney (MDCK) monolayer in vitro on flat 2D surfaces as a paradigm for morphogenesis in general (du Roure et al., 2005; Angelini et al., 2010; Poujade et al., 2007; Tlili et al., 2020). Collective effects in epithelial monolayers in vitro have been reported with new insights into the physics of active matter and developmental biology (Jain et al., 2021). The concepts of topological defects and their importance in setting new rules for morphogenesis (Fardin and Ladoux, 2021) and the influence of cellular flow patterns on stress distribution and cell fates (Chan et al., 2017) are just two of the many examples which show the impact of these new ideas. In addition, these concepts have also served as paradigms for new phenomena in vivo, for example, in identifying the key function of topological defects in hydra development (Maroudas-Sacks et al., 2021), or determining the principles behind the role of cellular flows and waves in vivo in major morphogenetic events in Drosophila (Collinet and Lecuit, 2021).

In this study, we focused on two dynamics of collective cell movements, pulsations and flows, which co-exist in standard 2D cell culture in vitro. We reasoned that the experimental control of each type of motion would guide us in understanding their cellular origins and successfully tested this assumption. First, we generated pulsations by modifying the friction between the tissue with surfaces patterned with fibronectin at the spontaneous coherence length scale of the monolayer. We then induced cellular flows by exposing the monolayer to cytoskeletal drugs. We also generated pulsations in conditions of flows by using laser ablation. To shed light on the underlying cellular mechanisms, we performed numerical simulations based on a unified vertex model with cell motility and successfully reproduced both types of motion by modulating just two factors.

Results

Characterizing the spontaneous flow in epithelial monolayers

We seeded MDCK cells on glass substrate and acquired time-lapse images after the monolayer reached confluency. We start with an initial seeding concentration of over and allow the monolayers to become confluent. As noted previously (Angelini et al., 2010; Poujade et al., 2007), confluent monolayers of MDCK underwent complex flow patterns that changed rapidly over time and eventually slowed down because of cell proliferation. To quantify the spatiotemporal patterns of cellular movements, we acquired the time-dependent velocity field using particle image velocimetry (PIV) on the time-lapse images (Figures 1A, 1D and S3A). The MDCK monolayer seemed to exhibit temporally periodic patterns of spatially contractile and expanding cellular movements. To quantify local area changes in the tissue because of these pulsatory patterns, we used the velocity field to calculate the flow divergence (see velocity field and correlation functions section in STAR Methods for details). The divergence and velocity fields showed that the flow field was correlated in space on a characteristic length scale of m (Figures 1C and 1F). Focusing on the domains of this size (Figures 1A and 1D), we noticed regions with visible pulses of contracting and expanding tissues (Figures 1A, 1C, and S3B and Video S1) for approximately 5 h, or areas where the flow appeared to be more uniform and correlated over a larger length scale (Figures 1D and S3C and Video S2). Although the pulsations were changing locations randomly, there were many stable pulsations that did not change their location over the observed dominant period of h, i.e., the divergence maxima was observed in the same location.

Figure 1

Spontaneous pulsations and flows in an epithelial monolayer

(A), (B), (D), and (E) represent the local characteristics of the monolayer. (A) shows snapshots of the contraction and expansion phases of one period of a pulsating domain, and (D) shows snapshots of flows occurring in the same domain at a different time point. In (A) and (D): the first row corresponds to phase-contrast images, and tracking highlights the back-and-forth movement of cells during pulsation in (A) and flow trajectory in (D); the second row shows the winding number analysis for the snapshots in the first row, and the third row shows the corresponding divergence and velocity magnitude plots, respectively, for (A) and (D). Scale bar . Time is in hh:mm. In those panels indicating the nature of the divergence field, blue indicates contraction, and red indicates expansion; similarly, in those panels indicating velocity magnitude, red and blue indicate low and high velocities, respectively. The numerical values are specified in the corresponding color bars. (B) The mean divergence plot for the pulsatile domain over 48 h represented in (A) and in Video S1. (E) The mean velocity magnitude plot of the flow domain shown in (D) and in Video S2. Arrows in (B) indicate the data points corresponding to snapshots in (A). (C) and (F) represent the global characteristics of the monolayer and show the divergence and velocity correlations, respectively, for distance and time for the same experiment over 48 h. See also Figures S1 and S3, and Videos S1 and S2.

Video S1. Pulsatile domain extracted from the MDCK-E-cadherin-GFP monolayer, related to Figure 1

Cells spontaneously exhibit collective contraction and extension phases on reaching confluency. Time in hh:mm.

Video S2. Domain of cells exhibiting flows, extracted from MDCK-E-cadherin-GFP monolayer, related to Figure 1

Cells display mainly flow behavior. Time in hh:mm.

We use the structure of velocity field vis-a-vis streamlines and the presence of topological defects (see second panel of Figures 1A, 1D, and S1, and winding number section in the STAR Methods) to distinguish between the two types of tissue kinematics, pulsations and flow, presented in Figures 1A and 1D (first panel), respectively. In the monolayer region with pulsations (Figure 1A), we typically saw 2–4 pronounced pulses in the velocity divergence (see Figure 1B and the autocorrelation plot in Figure S3D). Moreover, the structure of the velocity field is complex, and consequently we observe a number of +1 and −1 defects (Figure 1A, second row). On the other hand, in the monolayer domain with flow (Figure 1D, second panel), the velocity streamlines have a simpler structure resulting in far fewer topological defects. However, we clarify that even in this case, there is an underlying pulsatory velocity pattern (Figure S3C), but it is overwhelmed by the overall flow of the cells. Nonetheless, the coherent flow within the given domain does not imply large-scale directionality because the length scale associated with the velocity is around (Figures 1C and 1F), which is chosen to be the approximate size of the domain in Figures 1A and 1D (first panel), and is much smaller than the total field of view of the monolayer (Figure S3A).

The classification of monolayer velocity patterns within the given domains as pulsations or flows can provide a basis for further perturbation experiments. However, to quantify the tissue kinematic behavior over wider length scales, we calculated correlation functions in space and time for the velocity field and its divergence for the entire region under experimental consideration (Figures 1F and 1C, respectively, and Figure S3A). We determined the correlation function of the divergence at the same spatial but different time points, , as a function of time-interval τ. It exhibited oscillatory behavior with a period of about 5 h, indicative of periodic pulsatile movements occurring in the tissue. By calculating the spatial correlation function, , as a function of , we found that the divergence was also correlated in space on a length scale of (Figure 1C), which corresponds to approximately 15 cell lengths (see velocity field and correlation functions section in STAR Methods). Similarly, we noted a spatial decay in the velocity correlation function of that is consistent with the inherent correlation length for divergence (Figure 1F). The velocity field exhibited a temporal decay of a few hours, and is larger than the correlation time of the divergence. This finding is consistent with our observation of regions moving coherently with a relatively uniform velocity on long time scales (Figures 1D, 1E, and S3E). Taken together, divergence and velocity correlation plots highlight different collective features of cell motion in the monolayer. Meanwhile, we also found that the underlying rotational strength of the flow field was comparable to divergence both in terms of correlation length and time (Figure S3F).

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. To this end, by using micro-contact printing (see micro-contact printing section in STAR Methods), we deposited the extracellular matrix protein fibronectin as square motifs of controlled dimensions, called m, m, m (see Figures S4A–S4C, 2A, S5A, and S5D and Videos S3, S4, and S5). We kept the width of the fibronectin regions to be the same for all three cases while varying the size G of the central region. We reasoned that fibronectin coating would locally increase the friction between the tissue and its substrate and separate the entire tissue region into areas of high and low friction (Figures S4A–S4C). We chose the dimension of these patches to be around the size of the domains of spontaneous cell pulsations. We then quantified the temporal and spatial correlation of the flow divergence for the entire region. As in the control case, here also we observed the presence of oscillations with period of around 5 h (Figures 2C, S5C, S5F, and S4E). But the correlation lengths were less variable across experiments with Fibronection (FN) grids, compared to control (Figure S4D), implying that the difference in friction because of FN was tuning the size of pulsatile domains. 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 m case (Figure 2 and Video S4). In contrast, when the grid-size was smaller (m, Figures S5A–S5C and Video S3) or larger (m, Figures S5D–S5F and Video S5), pulsations were less pronounced. It is worth noting that the magnitude of divergence of tissue flow was increased in the larger grid case () (Figures 2B,S5E, S5G and S5I) compared to the smaller grids () (Figures S5B–S5H). This observation is consistent with the notion of lower overall friction between the tissue and its substrate when the fraction of substrate area covered by fibronectin is lower.

Figure 2

Substrate friction determines the localisation of pulsations

(A and B) represent the local characteristics, and (C) illustrates the global characteristics. (A) shows a typical sub-unit of the fibronectin (FN) grid shown in Figure S3B. Here, the first row shows snapshots of one contraction expansion cycle of a pulsatile domain in a single grid unit where tracking highlights the cell trajectories during pulsations. Similarly, the second and third rows show the topological defects and divergence, respectively, for the snapshots in the first row. Color bar indicates the scale for divergence maps where blue and red indicate contraction and expansion, respectively. Scale bar, and time in hh:mm. (B) shows the mean divergence plot that highlights the variations in divergence over 48 h for the FN grid domain in (A) and in Video S4. Arrows indicate data points corresponding to the snapshots in (A). (C) The time correlation function for the corresponding divergence field exhibits striking oscillatory behavior with a period of approximately 5 h. The divergence correlation function in space also shows a minimum at approximately , which reflects the periodic pulsatory pattern in space. See also Figures S4, S5, and S9 and Videos S3, S4, S5, and S6.

Video S3. FN grid (G-75; W-120) condition, related to Figure 2

Pulsatile domain extracted from MDCK-E-cadherin-GFP monolayer plated on fibronectin grids of width (W) and a gap (G) (non-fibronectin area) of . Fibronectin grid is shown in red and the gap is shown in gray. Time in hh:mm.

Video S4. FN grid (G-150; W-120) condition, related to Figure 2

Pulsatile domain extracted from MDCK-E-cadherin-GFP monolayer plated on fibronectin grids of width (W) and a gap (G) (non-fibronectin area) of . Fibronectin grid is shown in red and the gap is shown in gray. Time in hh:mm.

Video S5. FN grid (G-300; W-120) condition, related to Figure 2

Pulsatile domain extracted from MDCK-E-cadherin-GFP monolayer plated on fibronectin grids of width (W) and a gap (G) (non-fibronectin area) of . Time in hh:mm.

In these different conditions, more than of the grids exhibited sustained oscillations. Oscillations occurred within h, corresponding to a couple of cycles in the m case, six cycles in the m case, and a couple of cycles in the m case. Decay in oscillation amplitude was possibly associated with cell proliferation leading to reduced tissue motion over time. We also noticed that the number of defects was larger in the m and m cases than in the m case (Figures 2A, S5A, and S5D). Finally, in the ‘resonant’ case of m, we aimed to decrease friction in the central, non-adhesive region of the pattern through incubation with poly(l-lysine) grafted with poly(ethylene glycol) (PLL-g-PEG), which passivates the surface. This process led to ‘agitated’ pulsations, involving localized rapid motion of cells within the passivated square (Video S6), and this further supports the role of friction in the process. These results strongly suggest that friction can influence the nature of cellular movements. Taken together, these results demonstrate that the magnitude of pulsatile flow can be controlled by spatially modulating tissue-substrate interaction.

Video S6. FN grid (G-150; W-120) condition with passivation, related to Figure 2

Pulsatile domain extracted from MDCK-E-cadherin-GFP monolayer plated on fibronectin grids of width (W) and passivated regions (G) (pLL-g-PEG coating) of . White box shows the passivated area. Time in hh:mm.

Generation of long-range flows in epithelial monolayers

We then examined whether the correlation length of the flow could be experimentally modified. Interestingly, we found that incubation for 12 h with cytoskeleton inhibitors, followed by washing, strongly modulated the correlation length of the tissue-level velocity (see Figures 3A, 3B, and S6A–S6C, Videos S7 and S8). Specifically, on the addition of blebbistatin (see cytoskeletal inhibitors section in STAR Methods), tissue motion was significantly reduced (Figure S6A). However, about 2 h after subsequent washout of blebbistatin led to the appearance of flows with a longer spatial correlation but without change in persistence over time (Figures 3A, 3B, S6A, S4H, and S4I and Video S7), as quantified by the velocity correlation functions. The myosin retrieves its activity within minutes after washout but the cell and the monolayer have their own time to integrate activity and cause flows on length scale. A similar effect was also observed following the addition and washout of latrunculin-A (Figures S6B, S6C, S4H, and S4I and Video S8). We propose that blebbistatin washout promotes de novo monolayer collective order which is distinct from cellular arrangement obtained with standard plating protocol in which cells are randomly distributed at the early time points. We suggest that the formation of new cell-cell contacts after washout resets individual cell polarities to generate long range flow patterns.

Figure 3

The transition from pulsations to flows and back to pulsations

(A and B) Resetting myosin activity (incubation with blebbistatin for 12 h and washout) leads to the transition from pulsations to flows. (A) and (B) highlight the local and global flow characteristics, respectively. The first row of (A) shows snapshots of phase-contrast images where the tracking highlights flow behavior. The second-row shows streamlines without the presence of topological defects. Scale bar, . See also the associated Video S7. (B) shows the spatial and temporal correlation plots for velocity over 36 h. (C–E) Ablation of tissue performed after flows that resulted from resetting myosin activity leads to pulsatile behavior. (C) shows snapshot of phase contrast image with the ablated spot highlighted by dotted lines. Scale bar, . Time is in hh:mm. See the associated Video S9 and Figure S6D corresponds to the ablated region in (C). The first row of (D) shows snapshots of phase-contrast images where the tracking highlights back-and-forth pulsatile movements. The second and third rows show winding number and divergence respectively. Scale bar, . Time is in hh:mm; (E) shows the mean divergence plot over 24 h for the pulsatile domain shown in (D). Arrows indicate the data points corresponding to snapshots in (D). See also Figures S6 and S9 and Videos S7, S8, and S9.

Video S7. Blebbistatin washout condition, related to Figure 3

The video initially shows a pulsating domain in control condition from MDCK-E-cadherin-GFP monolayer. After addition of blebbistatin, pulsations are arrested. When blebbistatin is washed out after 12 h incubation, the same domain starts to exhibit long range flows. Time in hh:mm.

Video S8. Latrunculin-A washout condition, related to Figure 3 and S6

The video initially shows a pulsating domain in control condition from MDCK-E-cadherin-GFP monolayer. After addition of latrunculin-A, pulsations are arrested. When latrunculin-A is washed out after 12 h incubation, the same domain starts to exhibit long range flows. Time in hh:mm.

We had observed that pulsations in the MDCK monolayer started around the time when the cells were about to become confluent by the closure of gaps or ‘wounds’ in the tissue. Assuming that an empty space can lead to pulsatile motion of cells, we generated a wound using laser ablation after washout of the cytoskeletal inhibitor when the tissue exhibited a strongly coherent flow (see Figure 3C and S6D, and optical setups and imaging conditions section in STAR 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 (Figures 3D and Video S9). This observation was further confirmed by the measurement of divergence and its correlation (Figures S4F, S4G, 3E, and S6E). These results indicate that a combination of cytoskeletal remodeling and physical perturbations can toggle collective cell movements between pulsatile and long-range flows.

Video S9. Transition from flows to pulsations, related to Figure 3

After 12 h following blebbistatin washout, a large wound () is induced by laser ablation in MDCK-E-cadherin-GFP monolayer. The dotted box in the beginning of video shows the ablated spot. During the wound healing process that followed laser ablation, cells exhibited oscillatory phase highlighted by smaller continuous rectangle in the top right. Time in hh:mm.

Cellular mechanisms behind pulsatile flows

Following the same idea as laser ablation, where wound closure led to pulsations, we then sought to generate controlled wounds in the tissue to test their effect on the flow properties. Cells were plated on the substrate with microfabricated pillars (Anon et al., 2012), 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 (Figures S7A and S7B and Video S10). We took advantage of this set-up to visualize the cellular myosin distribution in cells participating in pulsatile flows (Figure 4A). We noticed that contracting zones appeared to exhibit an increased density of myosin. We also observed the transient appearance of a radial gradient of myosin fluorescence intensity around the wound, with larger intensity values toward the center of the wound in the contracting phase. This gradient was more homogeneous in the expanding phase (Figures 4B, 4C, S7C, and S7D and Video S10). These results suggest that myosin gradients participate in driving pulsatile flows in the tissue.

Figure 4

Molecular actors associated with pulsations and flows

(A–C) Fluctuation in myosin density during the contraction and expansion phases of pulsation. Micropillar assay was used to create wounds in the monolayer. Following wound closure, pulsations were observed in these sites. (A) Snapshots showing the fluctuation in myosin-GFP intensity. The first and the second time points show the wound and its eventual closure, respectively. Subsequent time points show the densification and diffusion of myosin. Scale bar, and time is in hh:mm. (B) shows myosin intensity profiles obtained along the yellow stripe in (A). Myosin intensity is initially zero at the center () of the wound (indicated by blue arrows) and increases to the half-way mark during wound closure. However, as the pulsation begins with contraction (at 12:40 in (A)), myosin intensity reaches the maximum and diffuses back to the half-way mark as the domain extends. See the associated Video S10, Figures S7A and S7B. (C) Divergence and intensity plots for the pulsatile domain shown in (A). The green region represents the wound closure phase, and the pink region represents the pulsation phase. (D and E) Orientation of lamellipodia aligning with the flow direction (on resetting myosin activity through incubation with blebbistatin for 12 h and washout). The first two snapshots in (D) represent the timepoint before flow initiation; the transition to flows is shown by the third time point and the subsequent time points show the flow phase. (D) The lamellipodial orientation of a cell transiently transfected with LifeAct-mCherry, shown in red, is superimposed with the direction of flow field vectors (green). Scale bar, , and time is in hh:mm. (E) Plot showing the linear tendency between lamellipodia orientation and the mean flow direction of the domain. The lamellipodia direction (obtained from the effective cell direction - see velocity-polarisation correlation in experiments section in the STAR methods) in the x-axis, is obtained for the cell shown in (D). See also Figures S7 and Video S10.

Video S10. Myosin density fluctuation during pulsation, related to Figure 4

Pulsatile domain extracted from MDCK-myosin-GFP monolayer. Monolayer undergoes wound closure with removal of micropillars followed by pulsatile activity. Phase contrast images of the cells are shown (left) and GFP-tagged myosin (right). During the pulsation that follows wound closure, accumulation of myosin can be seen with contraction. The accumulated myosin diffuses back to original state during extension phase. Time in hh:mm.

We also tested whether cell lamellipodia orientation was correlated with tissue flow. To follow lamellipodia dynamics, we used mosaic experiments with cells transiently transfected with Lifeact-mCherry (Figures 4D and S7E and Video S11; also see cell culture and sample preparation section in STAR Methods). Using blebbistatin washout to induce large-scale correlated flows, we found that lamellipodia direction indeed correlated with the movement of cells (Figures S7E, S7F, 4E, S7G, and S7H, see velocity-polarisation correlation in experiments section in STAR Methods). We chose to quantify the lamellipodium orientation by tracking and segmenting the outline of cells labeled with Lifeact-mCherry and calculated a nematic order parameter from the cell contour. We found that the nematic axis associated with cell elongation was then typically aligned with the direction of the flow (see velocity-polarisation correlation in experiments section in STAR Methods, Figure S2). These observations support that lamellipodia extend along the direction of flows.

Video S11. Orientation of lamellipodia with the flow field, related to Figure 4

Domain exhibiting flow behavior extracted from MDCK-E-cadherin-GFP monolayer transiently transfected with LifeAct-mCherry. The cell shown in red is transfected with LifeAct-mCherry. The flow field vectors are shown in green. Time in hh:mm.

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.

Numerical simulations to model tissue pulsatile flows

To identify the origins of the pulsatile collective flow, we ran numerical simulations of the tissue flow based on a vertex model (Comelles et al., 2021). In our model, cells were assigned a polarity vector that evolved over time according to a local alignment rule and polarity diffusion. The cells were subjected to a motile force along the polarity axis (Figure 5A).

Figure 5

Simulations of pulsations in friction grid and flows, using active vertex model

(A) Basic schematic of the vertex model and the source of forces on the vertices. The outcome of simulations for friction grid (B–D) and flows (E–G) correspond to their experimental counterparts in Figures 2 and 3, respectively. In (B), the first row shows: periodic contraction and expansion of a local domain with the back-and-forth movement of cells. The second and third rows show the topological defects from winding number analysis and divergence, respectively. Color coding: blue indicates contraction and red indicates expansion in the panels showing divergence field. See the associated Video S13.

(C) Total divergence for the local domain shown in (B) as a function of time. (D) The time and space correlations for divergence over the entire simulation. As observed in the experimental findings, length and time scales emerged from the divergence correlation plots. In (E), the first row shows the presence of long-range flows in the selected region. The second and third rows show the absence of topological defects and velocity magnitude, respectively. Color coding: red and blue indicate low and high velocities in those panels representing velocity magnitude. See the associated Video S15. (F) The corresponding velocity magnitude within the domain in (E) is nearly constant in time. (G) Velocity correlation function over the entire simulation shows persistence (long correlation) in both space and time. For a detailed discussion on experimental units of time and length, see additional details for the computational model section in the STAR Methods. See also Figure S8 and Videos S12, S13, S14, and S15.

In the basic vertex model of the confluent monolayer, the epithelial cells were represented as polygons that were created with vertices (Figure 5A). Every pair of adjacent vertices was shared by an edge which, in turn, had a cell on each of its side. The position of the vertices and the connectivity of the cells defined the geometry and topology, respectively, of the tissue. The mechanical properties were assigned to the minimal vertex model via a work function W:

The first term of the work function modeled the isotropic mechanical resistance with stiffness from a cell α having actual area to any deviations from its preferred area . The second term represented the effective interaction energy between two connected cells α and β that resulted from a combination of cell-cell adhesion energy () and acto-myosin contractility along the shared edge length . The force acting on a particular vertex i with position was

The cells within the colony were polarized and had a tendency to migrate (Figure 5A). The motility behavior of an isolated polarized cell was modeled using a simple description in which the cell moved with speed along its polarity direction . For the vertex model, this self-propulsion tendency of cells was converted to an effective motile force on its constituent vertex i aswhere summation was over the polarities , respectively, of the cells that contain the vertex i, and η is the effective friction coefficient between the cell vertex and the substrate (Sussman, 2017; Comelles et al., 2021).

Initially, cells were created in a periodic box of size . Each cell is provided a polarity , where the angle ϕ was chosen randomly from . The time evolution of vertex positions was based on the equation

To temporally evolve the polarity of the cells we utilized the experimental observations depicted in Figure S2. Specifically, we assumed that the polarity vector for a cell α has two tendencies – (1) to align with the direction of its instantaneous velocity and (2) to undergo rotational diffusion – which mathematically translated to the following differential equation (Henkes et al., 2011; Soumya et al., 2015)

Here, ξ was the alignment strength of polarity and was the diffusive, rotational, Gaussian noise with strength . The differential equations were solved using an explicit forward scheme with a time-step . See STAR Methods for additional details. Our model had seven parameters: preferred cell area (), cell area modulus (), line tension of cell interfaces (), cell motility (), polarity alignment rate (ξ), rotational diffusion of polarity (), and substrate friction coefficient (η). From these parameters, a set of time scales can be defined (see STAR Methods in the additional details for the computational model section and Table 1).

Table 1

Parameters for the vertex model

Parameter	Symbol	Non-Dimensionalised form	Numerical Value
preferred cell area	A_0	A0′=A_0l02	1000
cell area modulus	K	K^′=Kl03f_0	0.003
edge tension	Λ	Λ^′=Λf_0	100
substrate friction	η	η^′=ηl_0f_0t_0	500−2000
motile speed	v_0	v0′=v_0t_0l_0	0.025−0.1
polarity-velocity alignment	ξ	ξ^′=ξt_0	0.05
polarity rotational noise	D_r	Dr′=D_rt_0	0.003−0.02
simulation time-step	Δt	Δt^′=Δtt_0	1
small edge length during T1	ϵ_close	ϵclose′=ϵ_closel_0	2
open edge length during T1	ϵ_open	ϵopen′=ϵ_openl_0	3

We first tested the simulations with periodic boundary conditions and no rotational noise on the cell polarity vector and failed to generate pulsations. However, flows of cells were observed in many of the versions (Marchetti et al., 2013). 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 were important in setting the nature of collective migration modes in the tissue (see additional details for the computational model section in STAR Methods). High values of ξ promote polarity alignment in some direction, and hence long-range flows. On the contrary, high values of promoted fast and short-range oscillations. When the relative values of and ξ were well-balanced, we found temporally periodic and spatially oscillatory migration patterns of cells similar to that experimentally observed.

We next sought to reproduce three experimental conditions (Figures 5, S8, and Videos S12, S13, S14, and S15). We numerically mimicked standard monolayer (Figures S8A–S8C and Video S12), monolayer on grid (Figures 5B–5D and Video S13), and flows after blebbistatin washout (Figures 5E–5G and Video S15). To perform these numerical experiments, we tested the parameters that could allow reproduction of the 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 in Figures 5D, 5G, and S8C, all of which show good agreement with our data. We needed to have , grid-dependent friction, and to produce, respectively, Figures S8B, S8F, S8C, 5C, S8D, 5D, 5F, S8E, and 5G. Interestingly, we could also mimic the agitated movement of cells observed in a fibronectin grid with a PEG coating at the center, which is expected to reduce cell-substrate friction in that region (Video S14). We note that it is the differential friction between the central region and the adhesive patch of the simulated grid that leads to the modification in the pattern of cellular movements. Interestingly, the velocity patterns were more pronounced when the difference in the friction between the two regions was greater. We use periodic boundary conditions for all the simulations. In Videos S12, S13, S14, and S15 in which simulations are compared with their experimental counterparts, a smaller spatial domain is extracted from the simulation region of the complete tissue.

Video S12. Isolated pulsatile domains from the simulation (right) and experiment (left), related to Figures 5 and 1

The parameters for the simulation are: , , , , , , , .

Video S13: Isolated domains from simulation (right) and experiment (left) for the FN grid condition, related to Figures 5 and 2

The FN grid condition is simulated by increasing the motility () and mobility () of the cells in the center relative to those at the periphery (FN coated). The parameters for the simulation are: , , , , , , , , , . The dimension of the central square and the total thickness (left + right) of the fibronectin grid is kept to be the same.

Video S14. Isolated domains from simulation (right) and experiment (left) for the FN grid condition with PEG, related to Figures 5 and 2

The FN grid condition is simulated by increasing the motility () and mobility () of the cells in the center relative to those at the periphery (FN coated). The contrast in these properties between the central and the peripheral region is enhanced when compared to the regular grid in Video 13. The parameters for the simulation are: , , , , , , , , , . The dimension of the central square and the total thickness (left + right) of the FN grid is kept to be the same.

Video S15. Isolated domains from simulation (right) and experiment (left) for the flows, related to Figures 5 and 3

The parameters for the simulation are: , , , , , .

Our model assumptions and predictions are consistent with experimental observations. When the velocity direction and polarity at any position in monolayer were weakly correlated, pulsation was the main mode of cell movements (Figures S9A and S9B and velocity-polarisation correlation in experiments section in the STAR Methods). In contrast, when the velocity-polarity correlation was strong, cell flows were dominant (Figures S9C and S9D). In our model, the velocity alignment term ξ tends to align cell polarity with its velocity direction, whereas the rotational noise perturbs this alignment. Taken together, these observations support our modeling results that when velocity alignment strength is dominant , we obtain long range cellular flows. However, when the strength of velocity alignment and random rotation are comparable , we generate pulsatory cell movements.

Taken together, these results reproduce the main observations of our experiments. We conclude that a simple combination of cell motility and 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 variations in tissue kinematics.

Discussion

Pulsations and flows in tissues are striking examples of collective motions. Although they have independently been studied earlier in MDCK monolayers (Petitjean et al., 2010; Poujade et al., 2007; Angelini et al., 2010; Zehnder et al., 2016, 2015; Peyret et al., 2019; Notbohm et al., 2016; Rodríguez-Franco et al., 2017; Petrolli et al., 2019), in this article, we generated both types of multi-cellular dynamics in the same epithelial monolayers. We controlled each of them by spatially modulating substrate friction and by introducing in and then washing-out cytoskeletal inhibitors. We revealed that these perturbations affected the external environment and internal cellular properties of the monolayers, leading to these distinct robust collective motions. We also successfully reproduced these dynamics in a minimal vertex model with cell motility that took rotational diffusion of polarity and polarity-velocity alignment as essential components. In this model, a simple change in respective weights of these two polarity mechanisms resulted in a switch between pulsations and long-range flows, showing a way to control collective cell motions in monolayers. Thus, our study shows that the very same epithelial layer can switch its dynamics using simple rules, and pulsations and flows appear as emergent properties of epithelial monolayers. Analogously, we could expect that tissues in vivo could also switch their dynamics in time and space with the mechanisms we report by tuning local friction or modulating polarity dynamics in the same organ. Future studies can test these hypotheses in vivo by tuning the local interaction between the neighboring tissues by locally overexpressing ECM proteins or by locally incubating organs with myosin inhibitors and visualizing their recovery live. Thus, our findings open up a new perspective on using these readouts and numerical simulations to test morphogenetic events in general.

It is interesting to note that the domains for pulsations and flows are about in length scales for time scales of about few hours. These similar values can be obtained by scaling arguments. The Maxwell time τ is the ratio of viscosity overYoung’s modulus. Taking viscosity of about and Young modulus of about (Guevorkian et al., 2010; Stirbat et al., 2013), we obtain τ of the order of hours. Cell motility speeds are about . The product gives . Altogether, these estimates support a simple viscoelastic behavior.

Collective cell migrations in vivo occur in an integrated manner, involving signaling and mechanics pathways. For example, in zebrafish, that the so-called lateral line undergoes collective migration through chemical gradients generated by the tissue itself (Lecaudey et al., 2008; Donà et al., 2013). Similarly, in Xenopus laevis neural crest cells, stiffening is involved in setting the mesenchymal to epithelial transition through mechanisms involving an interplay between tissue mechanics and local signaling (Barriga et al., 2018; Scarpa and Mayor, 2016). Still, mechanical events alone were sufficient to capture some morphogenetic events. Hence, the search for physical rules is essential to capture the deformation kinematics of tissues. For instance, cell motility per se was associated with the elongation of the amniote axis. Similarly, differential random motion of cells was reported to control somitogenesis in chicken embryos (Bénazéraf et al., 2010; Xiong et al., 2020; Bénazéraf et al., 2017). Even if gradients in fibroblast growth factor (FGF) signaling were also present as well in this case, the cellular dynamics were sufficient as the main driving force for elongation, supporting the need to look for physical principles piloting morphogenesis. Germ cells were also shown in many organisms to move toward the gonads through active and passive processes associated with the dynamics of the surrounding tissues (Paksa and Raz, 2015; Grimaldi and Raz, 2020). Our generic mechanisms of tissue dynamics for identical cell lines could serve as a framework to revisit these examples guided by the quantification procedure and theoretical methods developed in this study.

In this context, we herein report two types of collective cell motion and provide a generic framework to explain them. Pulsatile movements were observed in various conditions. For example, in Drosophila, pulsed contractions guided by acto-myosin are involved in apical constriction and abdominal morphogenesis (Martin et al., 2009; Solon et al., 2009; Pulido Companys et al., 2020). Moreover, oscillation may vary across the tissues, as shown in Drosophila abdominal morphogenesis (Pulido Companys et al., 2020), and the accompanying myosin gradients are reminiscent of the ones we report in the current study. Elsewhere, oscillations within filopodia during the heart formation of Drosophila are essential to optimize morphogenesis in a proof-reading mechanism (Zhang et al., 2018, 2020). Cardiomyocytes also undergo these oscillatory motions both in vivo and in vitro (Nitsan et al., 2016). This expanding list of pulsation incidences suggests that a self-organized ‘pacemaker’ associated with a simple pulsation mechanism of physical origin, which span various tissues, model organisms, timescales, and length scales, could be involved in these situations.

We propose the same extension from our study for describing flows in vivo(Ishii et al., 2020). A previous study has shown how forces drive epithelial spreading in zebrafish gastrulation through a flow-friction mechanism (Behrndt et al., 2012). During epiboly in the insect Tribolium castaneum, part of the tissue becomes solid and pulls the associated tissues undergoing rearrangement and direct flow (Jain et al., 2020; Münster et al., 2019). A similar mechanism was suggested during wing formation in Drosophila with the hinge and the connected blade (Etournay et al., 2015). Force balance associated with adhesion between the blastoderm and the vitelline envelope led to this long-range flow. In another context, acto-myosin cables were shown to lead to collective movement of cells in a manner reminiscent to a liquid flow (Saadaoui et al., 2020). These situations show similar cell flows to that in the present study. Our study suggests that a common mechanism across scales, tissues, and species may be related to the greater contribution of velocity alignment than rotational diffusion in polarity. Lamellipodia in cells and its actin dynamics could be the main motor in these long-range flows.

This search for simplicity in the multi-cellular mode aims at extracting generic principles of morphogenesis. The MDCK epithelial monolayer in this study is an excellent system to mimic any epithelial dynamics in vivo. As mentioned above, collective cellular motions and deformations are common in developmental biology. However, the focus is generally on the genetic and signaling pathways associated with their onsets. Our in vitro results suggest that our readouts of diffusion in polarity and cell alignment can be used in vivo to test mechanisms in complex systems during morphogenesis. Because these parameters are now identified, collective cell dynamics can be tested with the same framework across species. Considering the large number of model systems available, including Drosophila, C. elegans, to zebrafish, Xenopus, and mice, these self-organization rules could be explored experimentally and compared across organs and species.

Our current study focused mainly on the physical parameters (substrate friction, cell motility, diffusion in polarity, and polarity-velocity alignment) controlling these types of motion and their cellular readouts through lamellipodia distributions and myosin gradients. In the future, it would be interesting to study the coupling between mechanical processes and cellular signaling, with experimental measurements of both types of events and their inclusions in the associated computational framework (Hino et al., 2020). Indeed, recent studies have highlighted the connection of wave propagation in tissues with signaling. Mechanical events alter signaling activity and these feedbacks are mediated by mechanochemical interactions between cells through extracellular-signal-regulated kinase (ERK) activation (Hino et al., 2020). Investigations on such interplay between mechanical and chemical contractions and waves could generate new paradigms in morphogenesis and allow the identification of conserved principles of mechanics and signaling networks in living matter. Cell density is another important parameter that regulates the characteristics of monolayer dynamics in a time dependent manner Tlili et al., (2018). In our experiments, beyond 48 h, cell proliferation leads to jamming like effect where the large scale dynamics are arrested and only local movements are observed. Along this line, the role of cell density in transitioning between pulsations and flows is an interesting avenue to be explored both in terms of dynamics (Tlili et al., 2018) and the underlying molecular mechanisms (Kocgozlu et al., 2016).

Limitations of the study

In the current study, we generate cellular flows by washing out the myosin inhibitor blebbistatin of the confluent monolayer but do not perfom experiments to directly test the associated mechanisms. Future experiments could modulate cell-cell interactions by incubating DECMA-1 antibody to perturb cell-cell interactions after blebbistatin washout. We anticipate that at high concentrations of DECMA-1, cell-cell junctions would be separated whereby no collective motion would be observed. However, at intermediate detachment levels, we expect that the spatial length and persistent time of the cellular flows would be decreased. In contrast, other MDCK cell lines overexpressing cadherin could be generated in which case we expect that both the flow correlation length and persistent time of the flows would be larger. These experiments would help substantiate the key role of cell-cell junction lifetime in modulating the velocity alignment of cells.

Our experimental findings in this article are based on MDCK epithelial monolayers. However, for a more general understanding, other epithelial cells such as HaCat and Huvec cells could be used. Nevertheless, based on our understanding from MDCK experiments and simulations, we predict that the interplay between the levels of adhesion (cadherin) and force (myosin) would determine the types of dynamics – pulsations or waves. We anticipate the cells with higher levels of adhesion to facilitate polarity alignment and cellular flows, whereas cell lines with the larger forces may facilitate pulsations. At that stage, numerical simulations would allow to quantitatively test these predictions by associating appropriate model parameters.

In the current study, we do not systematically explore the role of myosin levels on tissue kinematics. For that, clones with controlled levels of myosin could be prepared in future experiments. The role of myosin levels on the length and time-scales of cellular flows after blebbistatin washout experiment could also be explored from these experiments. We propose that the over-expression of myosin would decrease the persistence length and time by reducing velocity alignment whereas downregulation of myosin would trigger the opposite results. Future experiments would also help disentangle the potential interplay of myosin roles in velocity alignment and polarity diffusion as is proposed theoretically.

In this work, we do not perform a detailed comparative study of the combined quantitative effect of myosin and cadherin on collective cell movements (pulsations or flows) observed in the monolayer. Mosaic experiments using the clones discussed above would allow to further probe the respective mechanisms that are involved. By using fluorescent reporters for cadherin and myosin with distinct colors, levels of adhesion and force could be measured with respect to the local patterns of monolayer kinematics. From such measurements, a more precise mapping between the pattern of collective cell migration and spatiotemporal levels of myosin and cadherin can be obtained.

In the current work, we do not explore the role of substrate stiffness in quantitatively dictating the patterns of collective cell migration. However, we propose predictions about how the cell migration patterns would be affected by the modulation of substrate stiffnesses. With soft surfaces, we expect pulsations to be facilitated based on the fact that deformations of substrates can accompany the contractile activity of cells. In this case, the oscillation period and the spatial length could both be larger. Analogously, we expect a similar increase in the length and timescale for flow conditions. With stiffer substrate, however, we do not expect major changes because our substrates are already much above the typical Young modulus of our epithelial monolayers.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Daniel Riveline (riveline@unistra.fr).

Materials availability

This study did not generate new unique reagents.

Experimental models and subject details

Cell culture and sample preparation

MDCK cells labeled with Green Fluorescent Protein (GFP) for E-Cadherin (gift from W. James Nelson lab, Stanford University), MDCK cells labeled with actin-mcherry/myosin-GFP (gift from Roland Soldner lab), and MDCK cells labeled with GFP for myosin (gift from Shigenobu Yonemura lab, RIKEN) were cultured at C in low-glucose Dulbecco’s modified Eagle’s medium (DMEM) supplemented with Fetal Bovine Serum (FBS) and Penicillin-Streptomycin antibiotics. Cells were maintained at a minimal confluency of and then trypsinized at confluency for experiments. Live experiments were performed in Leibovitz L-15 medium and the initial seeding density was maintained at for all experiments. For the mosaic experiment, MDCK-Ecadherin-GFP cells transiently transfected with LifeAct-mCherry were mixed with non-transfected cells. All experiments were repeated a minimum of three times.

Method details

Cytoskeletal inhibitors

To inhibit cytoskeletal proteins, these inhibitors were used at the indicated concentrations: blebbistatin () and latrunculin A (). For incubation with inhibitors, media from the dish was replaced with media with the inhibitor. During washout, inhibitor media was carefully removed and the sample was rinsed at least twice with fresh L-15 before finally adding fresh L-15 for further acquisition. In both blebbistatin and latrunculin A conditions, the inhibitor was added to the monolayer after 14–16 h of reaching confluency. The monolayer was then incubated with the inhibitors for 12 h and then washed out. Flows were observed in the following h starting h after washout.

Micro-contact printing

Micropatterning of fibronectin on glass coverslips was performed using the micro-contact printing procedure (Caballero et al., 2014). First, motifs were microfabricated as SU-8 molds on silicon wafers using standard photolithography technique. Next, the replica were molded on polydimethylsiloxane (PDMS) at a pre-polymer, crosslinker ratio – (V/V)) for micropatterning. During the micro-contact printing procedure, coverslips were hydrophilized with “Piranha” activation (Piranha: 3:7 parts of and followed by functionalization with (3-mercaptopropyl)trimethoxysilane for 1 h by vapor phase deposition. The coverslips were then dried at C for 2 h. Next, PDMS stamps with microstructures were hydrophilized by oxygen plasma activation (Diener Electronic) and incubated with TRITC or HiLyte 488-labelled Fibronectin () for 1 h. The dried stamp was then brought in contact with the activated coverslip and a gentle pressure was applied to transfer the fibronectin patterns from the PDMS stamps to the coverslip. The coverslip was then stored in phosphate buffered saline (PBS) solution at C until the time of cell seeding. If the non-printed zones were to be passivated, the micropatterned coverslip was incubated with PLL-g-PEG (poly-L-Lysine-grafted-PolyEthylene Glycol) () in 10mM HEPES (pH 7.4) solution for 20 min. Then the coverslips were subsequently rinsed with PBS before proceeding for cell seeding.

Wound healing experiments

Wound healing experiments were performed using microfabricated pillars (Anon et al., 2012). A PDMS stamp with micropatterns was incubated with pluronic acid for 1 h and then rinsed with PBS. The stamp was then activated using an plasma cleaner and the side with micropillars was brought in contact with an activated glass coverslip for bonding. The bonded stamp and coverslip were stored at C for 30 min to strengthen the bonding. After UV sterilization of the stamp and coverslip, cells at a concentration of approximately were added at the interface between the stamp and the coverslip. The cells filled the space between the micropillars and they were allowed to settle for 45 min. Once the cells settled, fresh DMEM media was added to the sample which was then stored at C for h to allow cell spreading between the micropillars. After h, the PDMS stamp was carefully peeled off from the coverslip and the sample was finally rinsed and filled with L-15 media before imaging.

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 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 used 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 washout, cells within an area of were removed by ablation. The ablation was performed with smaller Region of Interests (ROI) within the area to avoid damaging the neighboring cells. The ablation was performed with Leica 63x (1.4 NA) oil objective at 800 nm wavelength and mW laser power. All experiments were performed at C and images were acquired at 5, 10 or 20 min interval between frames. In order to prevent evaporation of media due to long-term acquisitions (48–60 h), samples were covered with a glass Petri dish or by adding 2 ml of mineral oil.

Velocity field and correlation functions

For all the analysis, the time point at which the monolayer reached complete confluency was assigned as time zero (). In order to obtain velocity field in the monolayer, Particle Image Velocimetry (PIV) was performed using the open source software PIVlab (Thielicke and Stamhuis, 2014). Specifically, the velocity components (, ) at particular time t were obtained by using Direct Cross Correlation (DCC) algorithm on the monolayer images at times t and 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 () for the velocity field was calculated numerically using finite difference derivatives on the square grid plots. Likewise, the curl of the velocity field, , where ω is the vorticity, is also numerically calculated. Thus for each of , and div, we obtain a size matrix, where and are the number of grid points in the x and y directions, respectively, and is the number of time frames. The spatiotemporal 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:where the averaging was done over the space-time grid points , and . The full correlation function was then azimuthally averaged to obtained , where . The spatial correlation function , whereas the temporal correlation function . Please note that for convenience we use the same expression for the different variants of the correlation function. Correlation functions for divergence were similarly obtained starting from the full expression

Length and time scales for pulsations and flows were extracted from the minima of the plots. We checked that the length and time values were consistent with actual motions measured on videos.

Winding number

From visual observations of the monolayer time-lapse videos, 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 Figure S1A, there are four nodes with velocity at the respective node. The change in the angle between the velocity vector at a given node and the next was estimated as

The expression determines if the sense of rotation from i to is clockwise or anti-clockwise. The second expression gives the smallest angle between i and . The total winding number k over the contour iswhere the loop is traversed in an anti-clockwise sense (Figure S1B). 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 ( and (ii) saddle point (). We found that the focal points of pulsation typically had a strength of and corresponded either to a source or a sink .

Additional details for the computational model

The active vertex model with cell motility is defined using (Equation 1), (Equation 2), (Equation 3), (Equation 4), (Equation 5), (Equation 6) of the main paper. In this model, both the absolute and the relative values of ξ (polarity velocity alignment) and (polarity rotational diffusion) 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 promote fast and short range oscillations. When the relative values of 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 grid-like 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 approximately remained the same both inside and outside of the fibronectin patch. Tissue dynamics in typical simulations is reported along with experimental videos (Video S12. Isolated pulsatile domains from the simulation (right) and experiment (left), related to Figures 5 and 1, Video S13: Isolated domains from simulation (right) and experiment (left) for the FN grid condition, related to Figures 5 and 2, Video S14. Isolated domains from simulation (right) and experiment (left) for the FN grid condition with PEG, related to Figures 5 and 2, Video S15. Isolated domains from simulation (right) and experiment (left) for the flows, related to Figures 5 and 3).

The time-lapse images from the simulations have dimension of 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 (Moisy, 2017; Thielicke and Stamhuis, 2014). This velocity field is smoothened using Gaussian filter of size using pivmat to reduce the noise arising from short-range spatial noise in the movements of individual cells (Moisy, 2017). The divergence of this smoothed velocity field is obtained numerically using finite difference method (Moisy, 2017).

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 (), time () and force () scales is shown in the accompanying table.The force scale 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 . The most obvious choice of the length scale is related to the cell-dimension . Based on these parameters, we can also extract additional length (L) and time (τ) scales.

Comparing simulation and experimental values

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

The simulation velocity field of cells was quantified using PIV on images of , sampled after , where is the simulation time-step and denotes non-dimensionalised value. The distance R in the correlation function of Figure 5D is reported in px units and time t in terms of the duration between two frames (). Since there are cells in the image, each cell length corresponds to

The minima of the correlation function for Figure 5D is between . Hence, the collective pulsatile movement occurs approximately on the scale of cell lengths. This number is similar to the experimentally observed collective length-scale of pulsation. Moreover, by noting the , we also estimate that

thus giving a correlation length of approximately , which is consistent with the experimental findings. This also provides the experimental measure for distance R in Figure 5. Going along these lines, we now make an estimate of the length scale () that is implicit in the simulations. From Table 1, the preferred area of a single cell . Taking , we obtain

The minima for the correlation function in Figure S8C for the control case is at around (frame sampling). This corresponds to a period of , as the duration between two frames is around . Since the experimentally observed pulsation time-period is h, we estimate

Since, we used , the value of is also the same as the timescale that is implicit in the simulations (Table 1). Also, the time difference between two frames of simulation time-lapse images is

This is the experimental measure for the time in Figure 5. We check the consistency of this number as follows. In simulation units, the value of cell motility . From Table 1 we can note thatwhere we have used the estimates of and made above. This number is consistent with the experimental findings for cell speed.

Velocity-polarisation correlation in experiments

The quantification for the alignment of lamellipodia direction with the local flows is shown in Figures 4E, S7F, S7G, and S7H. This quantification is obtained as follows. We obtain the velocity field of the flow domain and cell contour separately and plot them together. The 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 movement direction as a representative of the lamellipodia direction. Individual cells are tracked after blebbistatin washout over a duration of approximately 12 h, and the cells are outlined.

We then 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 are 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 Figure S2. As described in the previous paragraph, we already have tracked cells with outlines. We then obtain the centroid of the cell for every time-step t as follows: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:where the angle ϕ is defined between x axis and the line connecting the centroid and the boundary pixel as shown in Figure S2A. The polarity can then be defined as , where . The direction of cell velocity is obtained from the instantaneous 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 overall time steps of the quantity

The combined histogram for this quantity from multiple cases is obtained in Figure S2. The distribution is skewed toward , i.e., , the angle between the velocity and polarity is systematically biased toward smaller angles (Figure S2B).

We also obtained the connection between cell polarity, defined in terms of collective shape anisotropy, and cell velocity across different experimental conditions (Figure S9). We first acquired polarity field of a group of cells through cell shape anisotropy as calculated from the structure tensor of the time-lapse images for fibronectin grid and blebbistatin washout conditions. The actual analysis was performed using the OrientionJ plugin of FIJI (Püspöki et al., 2016; Schindelin et al., 2012). The velocity field v of the cells was already obtained for these time-lapse images using PIVlab as discussed in the main paper. The polarity and velocity were both sampled on a grid with indices across the tissue domain examined for each of the experimental condition. To make an esimate of the degree of alignment between polarity and velocity orientation , we calculated the cross-correlationwhere the average was done over the the spatial grid points for all time-frames t. Note that, as in the earlier paragraph, was defined to satisfy equivalence for cell shape anisotropy. We found that the value of was approximately and 0.55, respectively, for grid and blebbistatin experiments.

Quantification and statistical analysis

Graph plotting and statistical analysis were performed using GraphPad Prism 7 software. Kruskal-Wallis test was performed to test the significance between control and different conditions in Figures S4D–S4I. The statistical significance is indicated as ∗p < 0.05. The number of experimental repeats is mentioned in the figure caption and the data is expressed as Mean value SD.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Chemicals, peptides, and recombinant proteins
Dulbecco’s Modified Eagle’s Medium (DMEM)	Invitrogen	31885-049
Fetal Bovine Serum (FBS)	HyClone	10309433
Penicillin-Streptomycin	Invitrogen	11548876
Leibovitz L-15 medium	Invitrogen	11540556
Blebbistatin	Sigma Aldrich	B0560
Latrunculin A	Sigma Aldrich	L5163
PolyDimethyl Siloxane (PDMS)	Dow Corning	Sylgard 184
SU-8	MicroChem	N/A
Hydrogen Peroxide	Sigma Aldrich	516813
Sulphuric Acid	Sigma Aldrich	258105
(3-mercaptopropyl)trimethoxysilane	Fluorochem	S10475
TRITC labelled Fibronectin	Cytoskeleton	FNR01
HiLyte 488 labelled Fibronectin	Cytoskeleton	FNR02
Phosphate Buffered Saline (PBS)	Invitrogen	11530486
PLL-g-PEG (poly-L-Lysine-grafted-PolyEthylene Glycol)	SuSoS	N/A
Pluronic acid	Sigma Aldrich	P2443
Mineral oil	Sigma Aldrich	M8410
Experimental models: Cell lines
MDCK-GFP-E-Cadherin	James W. Nelson lab., Stanford University	N/A
MDCK-mCherry-Actiin-GFP-Myosin	Roland Wedlich-Söldner lab., University of Münster	N/A
MDCK-GFP-Myosin	Shigenobu Yonemura lab., RIKEN	N/A
Software and algorithms
Fiji	Schindelin et al. (2012)	https://imagej.net/software/fiji/
OrientationJ plugin, Fiji	Püspöki et al. (2016)	http://bigwww.epfl.ch/demo/orientation/
MATLAB	N/A	https://se.mathworks.com/products/matlab.html
PIVlab, MATLAB	Thielicke and Stamhuis (2014)	https://se.mathworks.com/matlabcentral/fileexchange/27659-pivlab-particle-image-velocimetry-piv-tool-with-gui
Pivmat	Frédérick Moisy., 2017	http://www.fast.u-psud.fr/pivmat/
Other
Code for the correlation analysis	This paper	https://github.com/minamdar/Pulsations-Flows