On fundamental cellular processes for emergence of collective epithelial movement

Collective motion is an emergent property, in which a multitude of individuals interacting with one another synchronize their movement in a group. It is observed in many living organisms over a wide range of length scales: birds in flocks, sperm cells of sea urchins, bacterial colonies, and cytoskeletal filaments (Parrish and Edelstein-Keshet, 1999; Riedel et al., 2005; Schaller et al., 2010; Zhang et al., 2010; Sumino et al., 2012). A characteristic common among these phenomena is that adjacent local interactions produce well-ordered behaviors in a long range, such as homogeneous clustering movements. For the emergence of such collective motions, a simple theoretical framework was first pioneered by Vicsek and colleagues (Vicsek et al., 1995). They considered particle dynamics under simplified assumptions motivated by biological systems; they showed a kinetic phase transition from disorder to order phases controlled by the increase in density. The versatility of a density-dependent mechanism for collective motion has been demonstrated experimentally in diverse biological group systems. However, situations involving multicellular motion of epithelial tissue are different.

Unlike those self-propelled particles with a high degree of freedom, epithelial cells form clusters through cell–cell adhesion machineries to maintain tissue integrity. This means that cell density in epithelial tissues is of little variation. Despite such unity, epithelial tissue movements include a variety of degrees of order; remarkably, dynamic movement of epithelial cells relative to one another can be observed in the development of the submandibular salivary gland (Larsen et al., 2006; Rørth, 2009). Then, what controls such different multicellular movements of epithelial tissue?

In previous studies, growth factors and adhesion-related molecules involved in multicellular dynamics have been identified by gene knockin/knockdown and pharmacological experiments (Vitorino and Meyer, 2008; Du et al., 2010; Theveneau and Mayor, 2012). MAP kinase, for example, propagates through epithelial sheets as a wounding cue, and the transmitted chemical signal gives the cells directions to move (Matsubayashi et al., 2004; Nikolić et al., 2006). In recent years, mechanical aspects have been focused on, and direct measurements of mechanical forces revealed that cells generate a long-range gradient of intercellular tension and actively remodel to minimize the free energy attributed to intercellular stress, called plithotaxis (Trepat et al., 2009; Tambe et al., 2011; Trepat and Fredberg, 2011). Although their contributions are significant for clarifying molecular evidence or innate cellular properties, a comprehensive understanding for the emergent transition of multicellular collectivity between disorder and order motion is still far from complete.

In this study, we employed integrated approaches combined with accurate quantification techniques and a mathematical model that implemented the quantitative data to reveal fundamental cellular processes underlying the emergence of multicellular collectivity. First, we showed that collective behaviors of cultured human skin cells in a monolayer sheet are dependent on Ca²⁺ concentrations. Next, we quantitatively measured cellular responses to the Ca²⁺ treatments, such as cell growth, cell division, and the strength of intercellular adhesion. In this part, we applied for the first time an original estimation method using femtosecond laser-induced impulse to measure the strength of intercellular adhesion. Incorporating the measured data into a minimal mathematical model, we demonstrated that strengthening intercellular adhesion is sufficient for generating collective cell movement. Finally, the model suggested that not only cell adhesion but also cell migration should be considered for the emergence of collective cell movement. On the basis of these analyses, fundamental cellular processes for the emergence of collective cell movement were summarized.

To precisely quantify cellular dynamics, we employed a simple experimental setting, in which multicellular movement occurred within a monolayer sheet on a dish. Motivated by an earlier work (Poujade et al., 2007), we adopted the free-injury experimental system, in which removing the polydimethylsiloxane (PDMS) sheet allowed confluent cultured cells to move to a free space (Fig. 1A). In our experiments, normal human epidermal keratinocytes (NHEK) derived from human skin were used. On the basis of the fact that Ca²⁺ concentration increases around wounded skin, we examined cellular behavioral responses to different Ca²⁺ concentrations included in culture medium, and the concentrations of Ca²⁺ were determined as 0.075, 0.15, 0.75, and 1.5 mM according to relevant references (Grzesiak and Pierschbacher, 1995; Lansdown et al., 1999). After the incubation for 18 hours, the intracellular Ca²⁺ concentrations became proportional to the Ca²⁺ solutions in the medium, suggesting that the cells for each treatment entered into different states (Fig. 1B, regression line: y = 0.2337x+0.7245, R² = 0.9784; n = 6 for each treatment). As shown in Fig. 1C, peeling away the PDMS sheets allowed the cells to begin to move (supplementary material Movie 1). This movement included composite cellular actions, such as cell division, cell migration, and relaxation from the compressed state due to the potential volume constraint.

Since the shape of the forward edge took on a wavelike appearance in lower Ca²⁺ treatments although it remained relatively straight in higher Ca²⁺ treatments, we took the average distance from the original position to the expanded edge in order to define the moving distance of the forward edge. The analysis revealed that the forward edge in the group of cells moved farther for higher Ca²⁺ concentrations (Fig. 1D, Jonckheere–Terpstra test, P<<0.001 for each time; n = 3 for each treatment).

To quantify the dynamics of individual cells, we conducted particle image velocimetry (PIV) analysis for the time-lapse images of NHEK movement captured every 10 minutes (Fig. 2A; supplementary material Movie 2). Using the results of the PIV analysis, we computed the average velocity of particles within a small fixed grid between two temporal sequential images (Fig. 2B,C; supplementary material Movie 3). The results matched well with our manual tracing data of the movement of several cells. While the trace lines tend to spread in all directions and there are a relatively large number of trace-line intersections at the low Ca²⁺ concentrations, the trace lines have biased directions and there are fewer intersections at the high Ca²⁺ concentrations (Fig. 2D–G). This clearly reflects the time-series data of the angle distribution of cell movement obtained by the PIV analysis (supplementary material Movie 4).

Then, we introduced two order parameters to measure the directional and cohesive cell movement, which characterize the collectivity of multicellular movement: (1) the degree of forward movement and (2) the spatial correlation length. The more directional/cohesive the movement becomes, the greater the degree of forward movement/spatial correlation length. See Materials and Methods for the definitions of these quantities. In this paper we regard these two quantities and the moving distance of the forward edge as the measures of collective cell movement. Fig. 2H,I show that more directional and cohesive movement was generated in the higher Ca²⁺ treatments (Jonckheere–Terpstra test, Fig. 2H: P = 0, Fig. 2I: P = 0.001; n = 4 for each treatment). Note that there were no significant differences in the speed of the movement among the Ca²⁺ treatments (Fig. 2J; Kruskal–Wallis test, P = 0.497; n = 4 for each treatment); nevertheless, greater forward movement of cell mass and more collective movement were achieved as the concentration of Ca²⁺ increased. To summarize the relationship between the Ca²⁺ stimulus and the multicellular movement of NHEK, the individual cells within the monolayer sheet moved to the open space with unity, and the cluster of cells moved farther as Ca²⁺ increased.

To reveal what were altered by the Ca²⁺ stimuli, we examined cellular processes directly related to the movement in the population level: cell–cell adhesion, cell division, and cell growth. First, we estimated the strength of intercellular adhesion quantitatively by a femtosecond laser-induced impulse (Hagiyama et al., 2011). When an intense femtosecond laser pulse was focused on culture medium under a microscope, a local explosion of the culture medium was induced at the laser focal point, and a resultant stress wave, with a magnitude around micro-Newton, propagated at the vicinity of the laser focal point (Hosokawa et al., 2011). When the laser was focused at a few tens of µm from the intercellular junction, the junction was broken up by the stress wave as shown in Fig. 3A–D (supplementary material Movie 5). When we assumed that the stress wave propagates spherically as a short wave packet, the energy to break the intercellular adhesion was inversely proportional to the square of distance between the focal point and the intercellular junction d. We here estimated the maximum distance d_max at which the intercellular adhesion was broken. Since the energy to break the intercellular adhesion can be regarded as the energy to maintain the adhesion, we can acquire the strength of intercellular adhesion.

From this analysis, we obtained the relative energy necessary to maintain intercellular adhesion for each Ca²⁺ treatment. Indeed, the strength of intercellular adhesion increased in response to the applied Ca²⁺ concentration (Fig. 3E; n = 22, 24, 22, and 26 for each 0.075, 00.15, 0.75, and 1.5 mM Ca²⁺ concentrations). Normalization of those values revealed that the strength of the intercellular adhesion was approximately 13-fold between the lowest and highest Ca²⁺ treatments (Fig. 3F). The relative intercellular adhesion strength for each treatment W, i.e. W^X for X = 0.075, 0.15, 0.75, and 1.5, was quantified and applied to mathematical modeling using the parameters mentioned below.

Furthermore, by immunofluorescence labeling, we confirmed that the activity of some intercellular adhesion molecules increased as the Ca²⁺ concentration increased (Fig. 3G). For example, E-cadherin/occludin, representative adhesion molecules forming the adherens junctions/tight junctions, were more activated in the higher Ca²⁺ treatments. Also activated was vinculin, a main actin-binding protein of apical junctions that is suggested to be a marker of mechanical cues from adjacent cells (Yonemura et al., 2010). Thus, some of the molecules involved in intercellular adhesion certainly contributed to the establishment of mechanical strength of intercellular adhesion.

Next, we examined the number of dividing cells using the BrdU labeling method. We counted BrdU-positive cells soon after the removal of the PDMS sheet or 24 hours after the removal. We found uniform distributions of the number of BrdU-positive cells in both space and time (Fig. 4A,B). In addition, there were no significant differences among Ca²⁺ treatments, indicating that cell division was not a key factor in the generation of collective cell movement (Fig. 4C; Kruskal–Wallis test, 2 hr: P = 0.52, 24 hr: P = 0.773; n = 5). Also, we found that there were no significant differences in single-cell volume or shape among the treatments. This suggests that changes in cell growth and shape did not contribute to the emergence of collective cell movement either (Fig. 4D–H; Kruskal–Wallis test, P = 0.069; n = 30 for each treatment).

Although there appears to be a positive correlation between the collectivity of multicellular movement and the strength of intercellular adhesion, it is unclear whether or not intercellular adhesion is a sufficient factor for generating collective cell movement. One appropriate way to demonstrate sufficiency is a constructive approach. Therefore, we used mathematical modeling to implement and utilize our quantitative data.

We employed the cellular Potts model, which represents each cell morphology as a cluster of connected lattice sites; in this model, the dynamics of a multicellular system proceed in the form of a generalized energy H to be minimized (Glazier and Graner, 1993; Merks and Glazier, 2005; Hirashima et al., 2009). The energy in our model is composed of the minimal factors necessary to capture the multicellular dynamics, such as interfacial energy, cell volume constraint, cell division, and cell migration (Fig. 5A). The interfacial energy per unit length J determines the adhesion strength of the interacting materials, such as cell–cell or cell–medium, and λ_v represents resistance to cell compression/expansion. The system transition occurs stochastically by a lattice-based Monte Carlo method; that is, the labeled value of a randomly chosen lattice site is attempted to be replaced by a different labeled value of its neighboring lattice site randomly chosen ⁠. The transition is realized by evaluating the change in energy ΔH associated with its replacement. Even if ΔH is positive, the replacement occurs at a given probability. This stochastic state transition reflects biological noise such as membrane fluctuations and is determined using a parameter for the magnitude of biological fluctuation β. See Materials and Methods for more explanation.

At every Monte Carlo step (MCS), assumed to be one minute, whether or not cell division occurs was trialed under constant probability μ. In our model, once a cell divides, it enters a suspension period for 24 hours, corresponding to the NHEK cell cycle (Dover and Potten, 1983; Dover and Potten, 1988).

In comparison with two extremes of the relative strength of intercellular adhesion W, we found that the model captured the morphological and dynamical features of NHEK movement. For the case of weakest cell–cell bonding W^0.075, the forward edge showed a wavelike shape and the direction of cell movement became scattered, resulting in irregular connections of neighbors. On the other hand, for W^1.5, the forward edge was relatively straight and the moving direction became highly biased, with a regular geometry of neighbor connections (Fig. 5C,D; supplementary material Movie 7). This is because the lattice replacement in the case of low-energy intercellular bonding is more susceptible to stochasticity than that in the case of high-energy intercellular bonding. Although the front–rear axis of each cell tends to be aligned sequentially from the front through the process in which cells move to the open space, the sequential axis alignments cannot be transmitted correctly in the case of weak intercellular adhesion because of the agitation by biological fluctuations.

Extensive numerical simulations over a wide parameter range showed that the model well explained the dependence of the measures for multicellular dynamics on the strength of intercellular adhesion (n = 2000 for each value of relative adhesion strength, Fig. 5E–H). Indeed, in response to the increase in cell–cell adhesion strength in the model, the movement distance of the forward edge became long, and the two order quantities for the collectivity of multicellular movement also increased as the experimental data show (Fig. 5E–G). The decreasing speed of individual cells as the Ca²⁺ concentration increased was not the same as in our observed data (Fig. 5H). However, the fact that the forward edge of collective cells moved faster for the higher Ca²⁺ treatments means that the speed of individual cells is not a dominant factor in collective cell movement but the directional and cohesive movement can be a dominant factor. We conclude that the strength of intercellular adhesion can be a sufficient factor for the emergence of collective cell movement.

For the further investigation of parameter dependence, we mapped each measure for collective cell movement onto two parameters, γ/λ_m and λ_v/λ_m, both of which represent cell immobility. The former is determined by intercellular adhesion and the latter by cell rigidity (Fig. 6A–D).

We found that there were optimal balances between the strength of intercellular adhesion and the intensity of cell migration for the generation of collective cell movement in each measure (Fig. 6A–D). Note that the gradient of each measure along γ/λ_m is greater than that along λ_v/λ_m, and that γ/λ_m can relatively scale each measure (Fig. 6E–H). Furthermore, even if the cell division rate varies from μ = 0 to twice as much as the observed data, optimal balances still exist, indicating that optimality should be general (Fig. 6E–H).

In summary, the mathematical model suggests that the balance between the strength of intercellular adhesion and the intensity of cell migration can determine the emergence of collective cell movement (Fig. 6I). For the case of tight cell–cell adhesion and weak cell migration, a flock of cells gets stuck, preventing the movement of a cell cluster. Oppositely, for the case of loose cell–cell adhesion and strong cell migration, each cell tends to move in a disperse direction without cooperation; therefore, directional and cohesive movement cannot be produced.

In this paper, we have shown that a balance between intercellular adhesion and cell migration is required for the emergence of collective cell movement. The importance of intercellular adhesion has been corroborated by a systematic approach based on accurate data acquisition of cellular processes and on mathematical modeling; however, that of cell migration remains to be validated experimentally. Because the NHEK in a single cell state became clearly less active than in a multicellular state, we could not measure the single cell movement. We consider that the detailed analysis of molecular and mechanical factors about cell migration leads to a better understanding of multicellular movement.

We emphasize that the intensity of cell migration λ_m introduced in the model contains multiple interpretations. λ_m can be regarded as the accuracy of front–rear polarity based on the localization of polarity-related signals and subcellular organelles. Also, it represents the magnitude of mechanical force required for cell migration driven by actin and microtubule accumulation as well. Although we presumed a constant intensity of cell migration λ_m without specific interpretations in this study, further analysis of the spatio-temporal data of λ_m enables its more detailed descriptions for multicellular dynamics.

We should keep in mind that individual epithelial movement in a crowded situation is determined by relative movement of the surrounding cells because their movement is constrained through intercellular adhesion that triggers mechanical mutual interaction between neighboring cells rather than interfacial contact alone. Thus, multicellular dynamics cannot be expressed by only an aggregation of single cell movement. Also, we should note that the direction of front–rear polarity each cell has is altered through the relative movement. Along with the analysis for the intensity of cell migration, we need further studies about the dependence of the dynamics of front–rear axis in multicellular movement on the strength of intercellular adhesion, although the novel update rule of front–rear axis proposed in our model can successfully reproduce the multicellular dynamics.

For those analyses, various approaches will be more necessary in order to deepen our understanding of collective cell movement. In addition to prevailing techniques of genetic engineering, measurement and controllable device are need for quantification of biological phenomena, and mathematical model with the obtained data certainly become a powerful tool for the understanding and the prediction (Szabó et al., 2006; Bindschadler and McGrath, 2007; Lee and Wolgemuth, 2011; Vitorino et al., 2011; Puliafito et al., 2012; Vedula et al., 2012; Vignaud et al., 2012). We believe that this study serves the groundwork in terms of integrated approach between experiments and mathematical study.

NHEK cells (Clonetics) were cultured in the manufacturer-recommended medium, KGM-Gold (Clonetics) with penicillin (10 units/ml)–streptomycin (10 µg/ml) mixed solution (Nacalai Tesque) at 37°C in 5% CO₂. KBM-Gold basal medium without Ca²⁺ (Clonetics) was prepared for media that included different Ca²⁺ concentrations.

Fluo-4 NW (Invitrogen) was used as a fluorescent indicator for the intracellular free Ca²⁺. We applied the Fluo-4 NW into confluent cells and incubated the cells at 37°C for 30 minutes to obtain optimal results, as instructed in the product manual. The Fluo-4 NW was excited by a 488 nm line of argon-ion laser (CVI Melles Griot).

The PIV and statistical analysis code were written in MATLAB (MathWorks) based on MatPIV1.6.1, freeware distributed under the terms of the GNU general public license. In the PIV analysis, each image (512×512 pixel²) was broken into 4096 regions of 16×16 pixel² with 50% overlap for each region. This was determined to optimize both the number of cells in the observation window, approximately 60–70 cells in one dimension, and the accuracy dependent on the size of the evaluation window. We used the minimum quadratic difference method to maximize the accuracy (Gui and Merzkirch, 2000).

Three quantities for the quantification of multicellular movement were introduced. We calculated the moving distance of the forward edge as follows. The net moving distance of the tip of the cell mass along the x-axis was averaged over the y-axis. We manually determined the positions of the tips and took at least three positions from each sample in the measurement. In the computational simulations, we regarded cells that had fewer than four neighbors as ones that were away from the cluster, and excluded them from the calculation.

Statistical analyses were performed with a free software R using nonparametric alternatives to one-way analysis of variance. We adopted the Jonckheere–Terpstra/Kruskal–Wallis test when we emphasized difference/no difference among treatments.

The cultured NHEK cells were fixed and permeabilized by BD Cytofix/Cytoperm (BD Biosciences), and were treated with rabbit monoclonal anti-E-cadherin antibody (24E10, Cell Signaling Technology), mouse monoclonal anti-occludin (OC-3F10, Invitrogen), and mouse monoclonal anti-vinculin (7F9, Millipore). The first antibodies were reacted with AlexaFluor488 goat anti-rabbit IgG (Invitrogen) and AlexaFluor568 goat anti-mouse IgG (Invitrogen). DAPI (Millipore) and Hoechst33342 (Invitrogen) were used for nuclear staining.

A femtosecond laser pulse (780 nm, 250 fs) from a regeneratively amplified Ti: Sapphire laser system (IFRIT SP-1, Cyber Laser) was focused on a glass-bottom dish culturing NHEK through a 20× objective lens (N.A. = 0.46, UMPlanFL, Olympus). An impulse was generated by shooting the laser pulse with energy of 280 nJ at the culture medium. We set the focal point at the medium distant from the cell–cell junction and made the position closer to the junction until the junction was broken. Laser impulses were applied at intervals of several minutes with consideration for typical turnover of the adhesion molecules (e.g. 10–15 minutes for cadherins (Lambert et al., 2007)). In this procedure, we determined the maximum length d_max at which the targeted cell–cell junction was broken.

The impulsive force might be loaded not only locally to cell–cell junction but also, more or less, to peripheral region including cell–substrate adhesion, although we could not observe the separation of cell–substrate adhesion or the damage of cytosol. Since the separation of cell–cell junction occurred within milliseconds, the same time order as the propagation of stress wave, it is considered that the separation was not influenced by strain generated in the peripheral region.

We applied 10 µM of bromodeoxyuridine into the cell culture medium when we removed the PDMS sheet or 24 hours after the removal. In both cases, the cells were then incubated for 2 hours. BrdU-positive cells were detected with mouse monoclonal anti-BrdU/IdU (IU-4, Invitrogen) and AlexaFluor568 goat anti-mouse IgG (Invitrogen). Counterstaining was done with Hoechst33342 (Invitrogen). We manually counted the BrdU-positive cells and calculated the fraction.

NHEK cells were cultured at relatively low density and were induced by the Ca²⁺ treatments for 18 hours. Then, single cells within those populations were chosen randomly for the measurement. We manually traced the shapes of single cells and measured the area of each using ImageJ (NIH).

We set one pixel of the simulation space as 2 µm, and the computer simulation was conducted in a space 600×600 pixel², corresponding to the size of the observation window (1.2×1.2 mm²). The ideal cell volume was determined by images, i.e. V₀ = 100. At the initial time of the simulation, we assumed the direction of the front–rear axis had a uniform distribution, and that cell volume was slightly compressed, i.e. V₀ = 95. The initial variation of V₀ did not affect the results qualitatively.

As previously shown (Davies and Rideal, 1963; Glazier and Graner, 1993), the energy to maintain cell–cell adhesion is defined as γ = J_{cell–medium}−J_cell–cell/2, where J_cell–cell and J_{cell–medium} are the interfacial energy between cells and that between the cell and medium, respectively. In our study, the laser-used estimation revealed the relative strength of intercellular adhesion W among Ca²⁺ treatments, i.e. W^X = γ^X/γ^0.075 for X = 0.075, 0.15, 0.75, and 1.5. We assumed that the interfacial energy between cell and medium is proportional to that between cells because the density of active adhesion molecules tends to be related to the applied Ca²⁺ concentrations (Fig. 3G). We set J_{cell–medium} = 1.5J_cell–cell, leading to γ^X = J^X_cell–cell for X = 0.075, 0.15, 0.75, and 1.5; the qualitative results did not change when the constant value changed. For simulations over a wide range of parameter sets, parameter values were given as: γ = 0.1–10.0, λ_v = 0.1–10.0, λ_m = 1.0–10.0, and β = 1–2 unless otherwise noted.

This work was done by the support of Grant-in-Aid by CREST, JST to M.N., that by JSPS Research Fellowship for Young Scientists to T.H., that by JSPS Scientific Research (B) and MEXT Scientific Research on Innovative Areas to Y.H. We thank the following people: M. Denda, S. Denda, and M. Tsutsumi in Shiseido for their very helpful comments, I. Kii and Y. Yamaguchi in Kyoto University for their technical advice, and T. Miura in Kyoto University for partial financial support and the use of a laboratory.