Summary

In all animals, collective cell movement is an essential process in many events, including wound healing and embryonic development. However, our understanding of what characterizes the emergence of multicellular collective behavior is still far from complete. In this article we showed the fundamental cellular processes that drive collective cell movement by means of integrated approaches, including precise quantification measurements and mathematical modeling of measured data. First, we observed the dependence of the collective behaviors of cultured human skin cells on Ca2+ concentrations. When the culturing area confined by a PDMS sheet was suddenly expanded by removing the sheet, the group of cells moved to the expanded area with higher collectivity at higher Ca2+ concentrations. Next, we quantitatively measured cellular responses to the Ca2+ treatments, such as cell growth, cell division, and the strength of intercellular adhesion. Using a femtosecond-laser-based assay, an original method for estimating intercellular adhesion, we found that the strength of intercellular adhesion has an approximately 13-fold range in our treatments. Incorporating the quantitative data into a mathematical model, we then confirmed that the model well reproduced the multicellular behaviors we observed, demonstrating that the strength of intercellular adhesion sufficiently determines the generation of collective cell movement. Finally, we performed extensive numerical experiments, and the results suggested that the emergence of collective cell movement is derived by an optimal balance between the strength of intercellular adhesion and the intensity of cell migration.

Introduction

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 Ca2+ concentrations. Next, we quantitatively measured cellular responses to the Ca2+ 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.

Results

Measurement of multicellular movement

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 Ca2+ concentration increases around wounded skin, we examined cellular behavioral responses to different Ca2+ concentrations included in culture medium, and the concentrations of Ca2+ 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 Ca2+ concentrations became proportional to the Ca2+ solutions in the medium, suggesting that the cells for each treatment entered into different states (Fig. 1B, regression line: y = 0.2337x+0.7245, R2 = 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.

Epithelial behavior in response to Ca2+ treatments.

Fig. 1.
Epithelial behavior in response to Ca2+ treatments.

(A) Experimental settings. NHEK cells were cultured in a space confined by PDMS sheets. After Ca2+ induction for 18 hours, the sheets were removed. (B) Intracellular Ca2+ concentrations for the Ca2+ included in the medium. Panels in the graph show fluorescent images of the Ca2+ indicator. Error bar: s.d. (C) Time-series images of epithelial sheet for different Ca2+ treatments. Scale bar: 200 µm. (D) Dose response of the moving distance of the forward edge for Ca2+ treatments.

Fig. 1.
Epithelial behavior in response to Ca2+ treatments.

(A) Experimental settings. NHEK cells were cultured in a space confined by PDMS sheets. After Ca2+ induction for 18 hours, the sheets were removed. (B) Intracellular Ca2+ concentrations for the Ca2+ included in the medium. Panels in the graph show fluorescent images of the Ca2+ indicator. Error bar: s.d. (C) Time-series images of epithelial sheet for different Ca2+ treatments. Scale bar: 200 µm. (D) Dose response of the moving distance of the forward edge for Ca2+ treatments.

Since the shape of the forward edge took on a wavelike appearance in lower Ca2+ treatments although it remained relatively straight in higher Ca2+ 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 Ca2+ 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 Ca2+ concentrations, the trace lines have biased directions and there are fewer intersections at the high Ca2+ 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).

Measurement of cell motion.

Fig. 2.
Measurement of cell motion.

(A) Image of stained nuclei by Hoechst dye. (B,C) Direction and speed of cell motion obtained by the PIV analysis. (D–G) Tracking of movement of individual cells from the time of PDMS sheet removal (upper) to 18 hours later (bottom). Ca2+ concentrations: (D) 0.075 mM, (E) 0.15 mM, (F) 0.75 mM, and (G) 1.5 mM. (H) Degree of forward movement. (I) Spatial correlation length. (J) Speed of cells. Scale bars: 100 µm.

Fig. 2.
Measurement of cell motion.

(A) Image of stained nuclei by Hoechst dye. (B,C) Direction and speed of cell motion obtained by the PIV analysis. (D–G) Tracking of movement of individual cells from the time of PDMS sheet removal (upper) to 18 hours later (bottom). Ca2+ concentrations: (D) 0.075 mM, (E) 0.15 mM, (F) 0.75 mM, and (G) 1.5 mM. (H) Degree of forward movement. (I) Spatial correlation length. (J) Speed of cells. Scale bars: 100 µm.

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 Ca2+ 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 Ca2+ 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 Ca2+ increased. To summarize the relationship between the Ca2+ 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 Ca2+ increased.

Measurement of cellular processes: intercellular adhesion, cell division, and cell growth

To reveal what were altered by the Ca2+ 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 dmax 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.

Measurement of strength of intercellular adhesion.

Fig. 3.
Measurement of strength of intercellular adhesion.

(A) Schematic represents the propagation of the impulse generated by the femtosecond laser pulse into the medium. The distance between the laser focal point and the cellular interface is expressed as d. (B–D) Microphotographs (top) and their graphics (bottom) before and after the impulse loading. The black arrows point at the laser focal points. Note that there appeared to be a cleaved interface at the cellular junction (red arrow) at 33 ms after the impulse loading (C). (E) The threshold energy to maintain the intercellular adhesion, which is inversely proportional to the square of dmax, as dependent on the Ca2+ treatments. (F) Normalized strength of intercellular adhesion for the Ca2+ treatments. (G) Immunofluorescent images of active expression of adhesion-related molecules. Blue: DAPI. Scale bars: 20 µm for B–D, 50 µm for G.

Fig. 3.
Measurement of strength of intercellular adhesion.

(A) Schematic represents the propagation of the impulse generated by the femtosecond laser pulse into the medium. The distance between the laser focal point and the cellular interface is expressed as d. (B–D) Microphotographs (top) and their graphics (bottom) before and after the impulse loading. The black arrows point at the laser focal points. Note that there appeared to be a cleaved interface at the cellular junction (red arrow) at 33 ms after the impulse loading (C). (E) The threshold energy to maintain the intercellular adhesion, which is inversely proportional to the square of dmax, as dependent on the Ca2+ treatments. (F) Normalized strength of intercellular adhesion for the Ca2+ treatments. (G) Immunofluorescent images of active expression of adhesion-related molecules. Blue: DAPI. Scale bars: 20 µm for B–D, 50 µm for G.

From this analysis, we obtained the relative energy necessary to maintain intercellular adhesion for each Ca2+ treatment. Indeed, the strength of intercellular adhesion increased in response to the applied Ca2+ concentration (Fig. 3E; n = 22, 24, 22, and 26 for each 0.075, 00.15, 0.75, and 1.5 mM Ca2+ concentrations). Normalization of those values revealed that the strength of the intercellular adhesion was approximately 13-fold between the lowest and highest Ca2+ treatments (Fig. 3F). The relative intercellular adhesion strength for each treatment W, i.e. WX 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 Ca2+ concentration increased (Fig. 3G). For example, E-cadherin/occludin, representative adhesion molecules forming the adherens junctions/tight junctions, were more activated in the higher Ca2+ 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 Ca2+ 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).

Measurement of cell division and cell growth.

Fig. 4.
Measurement of cell division and cell growth.

(A,B) BrdU assay. Red: BrdU-positive cells, Blue: nuclei. The samples were collected from the area several millimeters behind the forward edge (A) and from the area of the edge (B). (C) Fraction of BrdU-positive cells for the Ca2+ treatments at two time points: 2 hours and 24 hours after the removal of the PDMS sheet. Error bars: s.d. (D–G) Images of a single cell in the Ca2+ treatments. Ca2+ concentrations: (D) 0.075 mM, (E) 0.15 mM, (F) 0.75 mM, and (G) 1.5 mM. (H) Cell size for the Ca2+ treatments. Error bars: s.d. Scale bars: 100 µm for A,B and 20 µm for D–G.

Fig. 4.
Measurement of cell division and cell growth.

(A,B) BrdU assay. Red: BrdU-positive cells, Blue: nuclei. The samples were collected from the area several millimeters behind the forward edge (A) and from the area of the edge (B). (C) Fraction of BrdU-positive cells for the Ca2+ treatments at two time points: 2 hours and 24 hours after the removal of the PDMS sheet. Error bars: s.d. (D–G) Images of a single cell in the Ca2+ treatments. Ca2+ concentrations: (D) 0.075 mM, (E) 0.15 mM, (F) 0.75 mM, and (G) 1.5 mM. (H) Cell size for the Ca2+ treatments. Error bars: s.d. Scale bars: 100 µm for A,B and 20 µm for D–G.

Mathematical modeling and model confirmation

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.

Mathematical modeling and model confirmation.

Fig. 5.
Mathematical modeling and model confirmation.

(A) Schematics of the cellular Potts model. (B) Migration energy and update rule of front–rear axis in the model. (C) Time-series images of cellular sheet to W0.075 and W1.5 in the simulation. Scale bar: 100 µm. (D) Direction of cell movements and the connections among neighboring cells to W0.075 and W1.5, taken from the black boxes in panel C. (E) Changes in the movement distance of the forward edge to different strengths of intercellular adhesion. (F) Degree of forward movement, (G) spatial correlation length, and (H) speed of cells at different strengths of intercellular adhesion. Samples are randomly collected from the following parameter ranges: λv = 5.0–10.0, λm = 1.0–5.0, T = 1.0–2.0, μ = 10−4. Error bars: s.d.

Fig. 5.
Mathematical modeling and model confirmation.

(A) Schematics of the cellular Potts model. (B) Migration energy and update rule of front–rear axis in the model. (C) Time-series images of cellular sheet to W0.075 and W1.5 in the simulation. Scale bar: 100 µm. (D) Direction of cell movements and the connections among neighboring cells to W0.075 and W1.5, taken from the black boxes in panel C. (E) Changes in the movement distance of the forward edge to different strengths of intercellular adhesion. (F) Degree of forward movement, (G) spatial correlation length, and (H) speed of cells at different strengths of intercellular adhesion. Samples are randomly collected from the following parameter ranges: λv = 5.0–10.0, λm = 1.0–5.0, T = 1.0–2.0, μ = 10−4. Error bars: s.d.

The migration energy Em is defined as:
formula
(1)
where represents the vector of the front–rear axis of cell migration, represents the vector connecting from the position of cell to a candidate position of the cell shifted by the replacement, and λm is the intensity of cell migration. The tilde over each letter indicates the unit vector (Fig. 5B). To update the front–rear axis accompanied by the state transition, we propose the following:
formula
(2)
where is the front–rear axis and the suffix i represents the cells that change position by the lattice replacement. We determined Eqn 2 according to an earlier study (Szabó and Czirók, 2010), and we also confirmed that Eqn 2 well reproduced characteristic dynamics of a few NHEK cells (supplementary material Movie 6).

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 W0.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 W1.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 Ca2+ 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 Ca2+ 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.

Model predictions

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

Model prediction.

Fig. 6.
Model prediction.

(A–D) Parameter dependence of measures for collective cell movement on γ/λm and λvm: (A) movement distance of forward edge, (B) degree of forward movement, (C) spatial correlation length, and (D) speed of cells. (E–H) Parameter dependence on γ/λm for various cell division rates μ, i.e. μ = 0, 10−5, 10−4, and 2×10−4. (I) Schematic summary of the model prediction. For the emergence of collective cell movement, an optical balance between the strength of intercellular adhesion and the intensity of cell migration is required.

Fig. 6.
Model prediction.

(A–D) Parameter dependence of measures for collective cell movement on γ/λm and λvm: (A) movement distance of forward edge, (B) degree of forward movement, (C) spatial correlation length, and (D) speed of cells. (E–H) Parameter dependence on γ/λm for various cell division rates μ, i.e. μ = 0, 10−5, 10−4, and 2×10−4. (I) Schematic summary of the model prediction. For the emergence of collective cell movement, an optical balance between the strength of intercellular adhesion and the intensity of cell migration is required.

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 λvm, 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.

Discussion

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.

Materials and Methods

Cell culture

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% CO2. KBM-Gold basal medium without Ca2+ (Clonetics) was prepared for media that included different Ca2+ concentrations.

Intracellular Ca2+ measurement

Fluo-4 NW (Invitrogen) was used as a fluorescent indicator for the intracellular free Ca2+. 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).

Particle Image Velocimetry (PIV) analysis

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 pixel2) was broken into 4096 regions of 16×16 pixel2 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).

Quantities for multicellular movement

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.

The degree of forward movement is defined as follows:
formula
(3)
where i is the index of cells, is the x component of the unit vector of cell displacement per 10 minutes, and N is the total number of cells in the observation window. A tilde and indicate the unit vector and time average, respectively. We assumed that single evaluation windows in the PIV analysis each correspond to single cells.
For the spatial correlation length, we first calculated the spatial correlation function φ:
formula
(4)
where is a deviation from its spatial mean . The spatial index k was converted to the distance in radial directions R. We then define the spatial correlation length as follows:
formula
(5)

Statistical analyses

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.

Immunofluorescence of cell adhesion molecules

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.

Cell–cell adhesion measurement

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 dmax 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.

Cell proliferation measurement

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.

Cell size measurement

NHEK cells were cultured at relatively low density and were induced by the Ca2+ 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).

Cellular Potts model

In the cellular Potts model, a state transition resulting from lattice replacement is preferred for decreasing the energy of the system H. H is defined in our model as:
formula
(6)
where J is the interfacial energy between cell–cell or cell–medium, λv is the magnitude of cell volume constraint, Vσ is the current cell volume, V0 is the ideal cell volume, and δ represents the Kronecker delta. The change in energy ΔH resulting from the state transition contains the energy for the cell migration:
formula
(7)
where Em is defined in Eqn 1. To realize the minimum system energy at equilibrium, the state transition is performed stochastically with biological noise such as cytoskeletal fluctuations in the following:
formula
(8)
where β represents the magnitude of the biological noise. See Glazier and Graner, Merks and Glazier, and Hirashima and colleagues for more detail about the cellular Potts model (Glazier and Graner, 1993; Merks and Glazier, 2005; Hirashima et al., 2009).

Parameters in computer simulation

We set one pixel of the simulation space as 2 µm, and the computer simulation was conducted in a space 600×600 pixel2, corresponding to the size of the observation window (1.2×1.2 mm2). The ideal cell volume was determined by images, i.e. V0 = 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. V0 = 95. The initial variation of V0 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 γ = Jcell–mediumJcell–cell/2, where Jcell–cell and Jcell–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 Ca2+ treatments, i.e. WX = γX0.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 Ca2+ concentrations (Fig. 3G). We set Jcell–medium = 1.5Jcell–cell, leading to γX = JXcell–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.

In our model, a cell divides with a constant probability μ at each MCS. For our simulations, the cell division probability was given as μ = 10−4 unless otherwise noted. The value was derived as follows. We here consider the expected increased fraction of the size of the cell mass from the initial state at a certain MCS. Because the fraction can be estimated from the measured data as at most 0.1 by 1440 MCS, we consider the following equation:
formula
(9)
The left-hand side of the equation indicates the probability that a single cell experiences at least one division by 1440 MCS. Once cell division occurs, the next cell division can be trialed after 24 hours in our simulations, because cell-cycle times in the confluency are on the order of 24 hours (Dover and Potten, 1983; Dover and Potten, 1988).

Acknowledgements

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.

References

Bindschadler
M.
,
McGrath
J. L.
(
2007
).
Sheet migration by wounded monolayers as an emergent property of single-cell dynamics.
J. Cell Sci.
120
,
876
884
.
Davies
J. T.
,
Rideal
E. K.
(
1963
).
Interfacial Phenomena, 2nd edition
.
New York, NY
:
Academic Press
.
Dover
R.
,
Potten
C. S.
(
1983
).
Cell cycle kinetics of cultured human epidermal keratinocytes.
J. Invest. Dermatol.
80
,
423
429
.
Dover
R.
,
Potten
C. S.
(
1988
).
Heterogeneity and cell cycle analyses from time-lapse studies of human keratinocytes in vitro.
J. Cell Sci.
89
,
359
364
.
Du
D.
,
Xu
F. L.
,
Yu
L. H.
,
Zhang
C. Y.
,
Lu
X. F.
,
Yuan
H. X.
,
Huang
Q.
,
Zhang
F.
,
Bao
H. Y.
,
Jia
L. H.
et al.  (
2010
).
The tight junction protein, occludin, regulates the directional migration of epithelial cells.
Dev. Cell
18
,
52
63
.
Glazier
J. A.
,
Graner
F.
(
1993
).
Simulation of the differential adhesion driven rearrangement of biological cells.
Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics
47
,
2128
2154
.
Grzesiak
J. J.
,
Pierschbacher
M. D.
(
1995
).
Shifts in the concentrations of magnesium and calcium in early porcine and rat wound fluids activate the cell migratory response.
J. Clin. Invest.
95
,
227
233
.
Gui
L.
,
Merzkirch
W.
(
2000
).
A comparative study of the MQD method and several correlation-based PIV evaluation algorithms.
Exp. Fluids
28
,
36
44
.
Hagiyama
M.
,
Furuno
T.
,
Hosokawa
Y.
,
Iino
T.
,
Ito
T.
,
Inoue
T.
,
Nakanishi
M.
,
Murakami
Y.
,
Ito
A.
(
2011
).
Enhanced nerve-mast cell interaction by a neuronal short isoform of cell adhesion molecule-1.
J. Immunol.
186
,
5983
5992
.
Hirashima
T.
,
Iwasa
Y.
,
Morishita
Y.
(
2009
).
Dynamic modeling of branching morphogenesis of ureteric bud in early kidney development.
J. Theor. Biol.
259
,
58
66
.
Hosokawa
Y.
,
Hagiyama
M.
,
Iino
T.
,
Murakami
Y.
,
Ito
A.
(
2011
).
Noncontact estimation of intercellular breaking force using a femtosecond laser impulse quantified by atomic force microscopy.
Proc. Natl. Acad. Sci. USA
108
,
1777
1782
.
Lambert
M.
,
Thoumine
O.
,
Brevier
J.
,
Choquet
D.
,
Riveline
D.
,
Mège
R. M.
(
2007
).
Nucleation and growth of cadherin adhesions.
Exp. Cell Res.
313
,
4025
4040
.
Lansdown
A. B. G.
,
Sampson
B.
,
Rowe
A.
(
1999
).
Sequential changes in trace metal, metallothionein and calmodulin concentrations in healing skin wounds.
J. Anat.
195
,
375
386
.
Larsen
M.
,
Wei
C.
,
Yamada
K. M.
(
2006
).
Cell and fibronectin dynamics during branching morphogenesis.
J. Cell Sci.
119
,
3376
3384
.
Lee
P.
,
Wolgemuth
C. W.
(
2011
).
Crawling cells can close wounds without purse strings or signaling.
PLOS Comput. Biol.
7
,
e1002007
.
Matsubayashi
Y.
,
Ebisuya
M.
,
Honjoh
S.
,
Nishida
E.
(
2004
).
ERK activation propagates in epithelial cell sheets and regulates their migration during wound healing.
Curr. Biol.
14
,
731
735
.
Merks
R. M. H.
,
Glazier
J. A.
(
2005
).
A cell-centered approach to developmental biology.
Physica A
352
,
113
130
.
Nikolić
D. L.
,
Boettiger
A. N.
,
Bar–Sagi
D.
,
Carbeck
J. D.
,
Shvartsman
S. Y.
(
2006
).
Role of boundary conditions in an experimental model of epithelial wound healing.
Am. J. Physiol.
291
,
C68
C75
.
Parrish
J. K.
,
Edelstein–Keshet
L.
(
1999
).
Complexity, pattern, and evolutionary trade-offs in animal aggregation.
Science
284
,
99
101
.
Poujade
M.
,
Grasland–Mongrain
E.
,
Hertzog
A.
,
Jouanneau
J.
,
Chavrier
P.
,
Ladoux
B.
,
Buguin
A.
,
Silberzan
P.
(
2007
).
Collective migration of an epithelial monolayer in response to a model wound.
Proc. Natl. Acad. Sci. USA
104
,
15988
15993
.
Puliafito
A.
,
Hufnagel
L.
,
Neveu
P.
,
Streichan
S.
,
Sigal
A.
,
Fygenson
D. K.
,
Shraiman
B. I.
(
2012
).
Collective and single cell behavior in epithelial contact inhibition.
Proc. Natl. Acad. Sci. USA
109
,
739
744
.
Riedel
I. H.
,
Kruse
K.
,
Howard
J.
(
2005
).
A self-organized vortex array of hydrodynamically entrained sperm cells.
Science
309
,
300
303
.
Rørth
P.
(
2009
).
Collective cell migration.
Annu. Rev. Cell Dev. Biol.
25
,
407
429
.
Schaller
V.
,
Weber
C.
,
Semmrich
C.
,
Frey
E.
,
Bausch
A. R.
(
2010
).
Polar patterns of driven filaments.
Nature
467
,
73
77
.
Sumino
Y.
,
Nagai
K. H.
,
Shitaka
Y.
,
Tanaka
D.
,
Yoshikawa
K.
,
Chaté
H.
,
Oiwa
K.
(
2012
).
Large-scale vortex lattice emerging from collectively moving microtubules.
Nature
483
,
448
452
.
Szabó
A.
,
Czirók
A.
(
2010
).
The role of cell–cell adhesion in the formation of multicellular sprouts.
Mathematical Modelling of Natural Phenomena
5
,
106
122
.
Szabó
B.
,
Szöllösi
G. J.
,
Gönci
B.
,
Jurányi
Z.
,
Selmeczi
D.
,
Vicsek
T.
(
2006
).
Phase transition in the collective migration of tissue cells: experiment and model.
Phys. Rev. E Stat. Nonlin. Soft Matter Phys.
74
,
061908
.
Tambe
D. T.
,
Hardin
C. C.
,
Angelini
T. E.
,
Rajendran
K.
,
Park
C. Y.
,
Serra–Picamal
X.
,
Zhou
E. H.
,
Zaman
M. H.
,
Butler
J. P.
,
Weitz
D. A.
et al.  (
2011
).
Collective cell guidance by cooperative intercellular forces.
Nat. Mater.
10
,
469
475
.
Theveneau
E.
,
Mayor
R.
(
2012
).
Cadherins in collective cell migration of mesenchymal cells.
Curr. Opin. Cell Biol.
24
,
677
684
.
Trepat
X.
,
Fredberg
J. J.
(
2011
).
Plithotaxis and emergent dynamics in collective cellular migration.
Trends Cell Biol.
21
,
638
646
.
Trepat
X.
,
Wasserman
M. R.
,
Angelini
T. E.
,
Millet
E.
,
Weitz
D. A.
,
Butler
J. P.
,
Fredberg
J. J.
(
2009
).
Physical forces during collective cell migration.
Nat. Phys.
5
,
426
430
.
Vedula
S. R. K.
,
Leong
M. C.
,
Lai
T. L.
,
Hersen
P.
,
Kabla
A. J.
,
Lim
C. T.
,
Ladoux
B.
(
2012
).
Emerging modes of collective cell migration induced by geometrical constraints.
Proc. Natl. Acad. Sci. USA
109
,
12974
12979
.
Vicsek
T.
,
Czirók
A.
,
Ben–Jacob
E.
,
Cohen
I.
,
Shochet
O.
(
1995
).
Novel type of phase transition in a system of self-driven particles.
Phys. Rev. Lett.
75
,
1226
1229
.
Vignaud
T.
,
Blanchoin
L.
,
Théry
M.
(
2012
).
Directed cytoskeleton self-organization.
Trends Cell Biol.
22
,
671
682
.
Vitorino
P.
,
Meyer
T.
(
2008
).
Modular control of endothelial sheet migration.
Genes Dev.
22
,
3268
3281
.
Vitorino
P.
,
Hammer
M.
,
Kim
J.
,
Meyer
T.
(
2011
).
A steering model of endothelial sheet migration recapitulates monolayer integrity and directed collective migration.
Mol. Cell. Biol.
31
,
342
350
.
Yonemura
S.
,
Wada
Y.
,
Watanabe
T.
,
Nagafuchi
A.
,
Shibata
M.
(
2010
).
alpha-Catenin as a tension transducer that induces adherens junction development.
Nat. Cell Biol.
12
,
533
542
.
Zhang
H. P.
,
Be'er
A.
,
Florin
E. L.
,
Swinney
H. L.
(
2010
).
Collective motion and density fluctuations in bacterial colonies.
Proc. Natl. Acad. Sci. USA
107
,
13626
13630
.

Competing interests

The authors have no competing interests to declare.

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial Share Alike License (http://creativecommons.org/licenses/by-nc-sa/3.0).