Dynamic self-organization of migrating cells under constraints by spatial confinement and epithelial integrity

Abstract Understanding how migrating cells can establish both dynamic structures and coherent dynamics may provide mechanistic insights to study how living systems acquire complex structures and functions. Recent studies revealed that intercellular contact communication plays a crucial role for establishing cellular dynamic self-organization (DSO) and provided a theoretical model of DSO for migrating solitary cells in a free space. However, to apply those understanding to situations in living organisms, we need to know the role of cell–cell communication for tissue dynamics under spatial confinements and epithelial integrity. Here, we expand the previous numerical studies on DSO to migrating cells subjected spatial confinement and/or epithelial integrity. An epithelial monolayer is simulated by combining the model of cellular DSO and the cellular vertex model in two dimensions for apical integrity. Under confinement to a small space, theoretical models of both solitary and epithelial cells exhibit characteristic coherent dynamics, including apparent swirling. We also find that such coherent dynamics can allow the cells to overcome the strong constraint due to spatial confinement and epithelial integrity. Furthermore, we demonstrate how epithelial cell clusters behave without spatial confinement and find various cluster dynamics, including spinning, migration and elongation. Graphical abstract


Introduction
The emergent mechanisms by which living systems acquire complex structures and functions have not been completely understood yet. Recently, the importance of self-organization processes in biological systems has been recognized not only in the physics community, but also in biology community [1][2][3]. Living systems need to demonstrate emergence of dynamic properties such as coherent dynamics and formation of dynamic patterns, as well as static properties. For example, cells often move cooperatively and spontaneously organize their coherent dynamics [4][5][6][7][8][9][10][11][12]. Accumulated experimental knowledge suggested that intercellular contact communication plays crucial roles for such dynamic selforganization of migrating cells [11,12], and the theoretical model which can recapture those dynamic phenomena in solitary cells has been provided [13,14].
One of the next challenges will be to investigate how dynamic self-organization occurs in living organisms. In multicellular living organisms, establishment of coherent dynamics can occur not only by freely-migrating solitary cells, but also by cells under various constraints. An important class of such cells is epithelial cell monolayers, where cells are tightly packed in a sheet-like tissue without gaps and adheres each other mechanically. A classic example of an epithelial cell line which a e-mail: mbithi@nus.edu.sg (corresponding author) exhibits this organization is Madin-Darby Canine Kidney (MDCK) cells [15,16]. It is indeed known that such epithelial cells exhibit coherent dynamics [15,16], which can be relevant to in situ biological processes like gastrulation [17]. Uncovering the mechanisms by which coherent migratory dynamics occur in epithelial cells may provide us new ways to understand mechanistic aspects in such biological processes. In addition to this, as living organisms have a finite size, it is important to consider spatial confinement as a constraint for establishing coherent dynamics. However, still there is a lack of knowledge regarding how cell-cell contact communication can contribute to dynamics of migrating cells under these constraints.
Here, we aim to expand on our earlier work demonstrating theoretical understanding of dynamic selforganization of migrating cells through cell-cell contact communication [14] to cases where migrating cells are subject to constraints of spatial confinement and/or epithelial integrity, i.e., mechanical constraint due to the integrity of an epithelial tissue. After investigating coherent dynamics in cells under spatial confinement but without epithelial integrity, we extended the study to investigate coherent dynamics under epithelial integrity. We assume that the cells interact with each other through mechanical mutual exclusion and two ubiquitous types of cell-cell contact communication, namely contact following and contact inhibition of locomotion [13,14]. For epithelial cell's case, numerical 16 Page 2 of 16 Eur. Phys. J. E (2022) 45 :16 simulations are performed based on the mathematical model of cell migration in an epithelial tissue using the following assumptions: (i) Each cell consists of apical and basal sides. (ii) Epithelial integrity at the apical side is modeled based on the cellular vertex model in two dimensions [18], and for brevity the constraints due to this vertex model at the apical side are termed apical constraint. (iii) The basal side of each cell is modeled by a single disk which spontaneously move on the twodimensional substrate according to its intrinsic polarity [19]. (iv) The apical and basal sides of each cell are connected to each other through a linear elastic spring. (v) Cell-cell contact communication as mentioned above takes place through the basal sides of cells.
Before closing this section, we make comments regarding with a brief review theoretical modeling. Epithelial tissues play an essential biological roles in shaping the multicellular body. Based on this important role, a model for migrating epithelial cells which implements the polygonal shape of cells by Voronoi tessellation (using the points which spontaneously move around), later termed the self-propelled Voronoi model, has been proposed in Ref. [20], and this model was further investigated in details from the perspective of statistical physics [21]. By assuming the tendency that cell polarity direction aligns to the actual direction of migration, flocking behavior was also investigated using this model [22]. Variations on this self-propelled Voronoi model have also been proposed. One example is based on Delaunay triangulation instead of Voronoi tessellation [23]. Furthermore, models combining spontaneous cell motility with purely the vertex description, which does not assume the cell shapes described by Voronoi volumes, have been introduced [24,25]. These concepts, i.e., a combination of spontaneous motility and polygonal shape description of cells, may also be a powerful tool to study cellular dynamic self-organization in an epithelial tissue, where the cells interact through the apical connections. However, for this aim, the motility of the basal side should be appropriately described; it is not the subordinate of the apical configuration, but can show independent migratory activity [26]. Therefore, the model used in this article takes into account such mutual independence of apical and basal dynamics by assuming both of them to be independent variables.
This article is organized as follows: The theoretical model which we investigate in this article is described in Sect. 2, with numerical results in Sect. 3. In particular, we investigate the cases of freely migrating cells under confinement, migrating epithelial cells (i.e., cells with apical constraints) under confinement, and migrating epithelial cells in a free space in Subsects. 3.1, 3.2, and 3.3, respectively. Section 4 is devoted to a brief summary and discussion. "Appendix A" provides the numerical results for freely migrating cells in a periodic system boundary, similar to Ref. [14] but with a smaller system size and a different set of parameters. "Appendix B" provides the numerical results of the case with apical junctional tension fluctuation. "Appendix C" provides the comprehensive discussion on how we determine the state boundaries in this article.

Model
As described in Sect. 1  Schematic illustration of the case with epithelial integrity. d Sample snapshots of numerical simulations based on this model. Circles indicate the cell periphery (i.e., half interaction range) at the basal side, and arrows represent the intrinsic polarity q of each cell mutual exclusions (Fig. 1a) and intercellular communication (Fig. 1b) on a two-dimensional flat substrate. Later, in part, we consider an epithelial monolayer in which each cell is spontaneously migrating (Fig. 1c). The basal side of each cell is modeled by a so-called selfpropelled particle, meaning that the particle is spontaneously moving on the substrate according to its intrinsic polarity [14,19]. This mimics the spontaneous migration ability of each cell on the basal substrate. The apical side of the cell is modeled based on the vertex dynamics framework in two dimensions, which can simulate in-plane motions of cells in an epithelial monolayer tissue [18]. In this framework, since the individual cells in an epithelial monolayer often have straight cellcell junctions at the apical face, the shape of each cell is described by a polygonal shape. The apical and basal sides of a cell are connected with each other through a linear elastic spring (Fig. 1c). The system is assumed to have a boundary over which the cells cannot travel (Fig. 1d), except when considering migrating cells in a free space (Sect. 3.3). We adopt the rectangular boundary with the aspect ratio of 20/7 for this, except for Sect. 3.3. This article aims at investigating the impact of contact intercellular communication for DSO of migrating cells under constraints. In particular, we focus on the following manners of intercellular communication, called contact inhibition/attraction of locomotion (CIL/CAL) and contact following (CF), as with our previous work [14]. CIL: When two cells come in contact, they modulate their polarities to avoid overlap and migrate away as if they are scattered (Fig. 1b left) [5,6,[27][28][29]. This is called CIL. It is known that, for example, CIL takes place in Xenopus neural crest cells [5] and is involved in coordinated migration and directional migration in response to external cues [6,27]. CAL: When two cells come in contact, they modulate their polarity such that they are mutually attracted (Fig. 1b center) [8]. This is called CAL and may result in their behaviors in population. CF: When a cell contacts another cell, the cell at the back chases the cell at the front but not vice versa (Fig. 1b right). This is called CF and has been observed in several cells [11,12,30]. CF is involved in their DSO, such as formations of rotating mount [11] and traveling density band [12].
We assume that, except for the self-propulsion mechanism at the basal side, the system obeys variational dynamics with the potential function E composed of three parts Each contribution in the potential, E vertex for the apical side, E VE for the basal side and E A−B for the apicobasal connection, is explained in detail in the following subsections. Here, r i (i = 1, 2, . . . , M with the total number of vertices M ) and X α (α = 1, 2, . . . , N with the number of cells N ) denote the variables specifying the configurations at the apical and basal sides of cells, respectively, the details of which are also explained in the following subsections. This manner of implementation of epithelial integrity, for which we introduce two independent layers of variables, reflects the experimental observations for dynamics of apical and baso-lateral sides of a cell, which can be independent with each other, as described in Introduction. This implementation may inherently be able to cause characteristic dynamics replying on the relative motion of apical and basal sides of each cell. However, we have not observed such characteristic dynamics requiring two layers of variables so far. Moreover, in this article, we implement epithelial integrity by the vertex model, following the convention for the simulations of epithelial tissue dynamics. However, the integrity can be also implemented using the disk-disk adhesion [12]. It might be another interesting question how this difference in implementation of cell-cell adhesion can affect the results.

Model for basal motility: Self-propelled particle model
For this part, we consider N particles in the twodimensional space surrounded by the system boundary. The particle with the indices α (α = 1, 2, . . . , N) represents the basal body of the α-th cell, and the location of it is described by X α = (X α , Y α ). The velocity v α ≡ dX α /dt and the intrinsic polarity q α of the α-th particle obey and respectively [12][13][14]31,32].
Equation (2) is the equation of motion, where Γ is the mobility coefficient and e is the coefficient representing effective driving force for migration, and the last term is the mechanical interaction between the particles with the potential E in Eq. (1). In the potential E, the parts relevant for the basal side dynamics are those for the volume exclusion E VE (Fig. 1a) and for the apico-basal connection E A−B (Fig. 1c). (For E A−B , see Subsect. 2.3.) For simplicity, we idealize each cell as a soft-core repulsive disk with a fixed diameter r, thus assume where r is the cell-cell interaction range and roughly means the cell diameter. The summation <α,α > runs for all the pairs of neighboring cells (α-th and α -th) satisfying |ΔX α α | < r. The second term, cell-wall volume exclusion, is defined similarly but with |ΔX wall,α | < r/2 for the interaction range.
Note that Eq. (2) can be rewritten as with the interaction term J v α = −Γ ∂E/∂X α and the coefficient v 0 ≡ Γ e giving the fixed speed of migration for the J v α = 0 case, or a migrating cell without any interaction. From Eq. (4), J v α is given by where ΔX α α = ΔX α α /|ΔX α α | and β ≡ Γ β . The summation α (n.α) runs for all the neighbors of αth cell satisfying |ΔX α α | < r. This definition is also applied to wall (n.α) but with the range r/2.
Equation (3) describes the time evolution rule for the intrinsic polarity q α of the α-th particle, or at the basal side of α-th cell. The first term of the right-hand side in Eq. (3) controls the spontaneous formation of the polarity of each particle, with the coefficient I q specifying the strength of polarity formation and the optimum magnitude Q. The second term J q α represents interaction between the particles through the intrinsic polarity, which is further specified below. The last term ξ α (t) specifies the noise for the polarity; this noise is assumed to be a white Gaussian noise with ξ α = 0 and where the subscripts k and l specify the directions, k, l = x, y. The coefficient D indicates the noise strength and has the unit of inverse time.
In this article, for polarity interaction J q α we adopt the following combination of contact communication through the basal sides of cells. Following Ref. [14], we assume that this term J q α is given by the additive combination of CF and CIL/CAL terms, i.e. J q α = J CF α + J CIL α , and these two terms are defined as follows; J CF α represents the CF contribution and J CIL α represents the CIL/CAL contribution (for α CIL > 0/< 0, resp.) Again, the summation α (n.α) runs for all neighbors of α-th cell satisfying |ΔX α α | < r. Here, α CF is the strength of contact follow (CF) and should be α CF ≥ 0, and α CIL is the strength of contact inhibition/attraction of locomotion (CIL/CAL). For CIL, α CIL > 0, whereas for CAL, α CIL < 0.
Here we briefly review how each of these terms affects the population dynamics of freely migrating cells, or self-propelled particles, described in our earlier work [14]. The CF term gives rise to a nonreciprocal local interaction with forward-backward asymmetry (Fig. 1b  right). In terms of symmetry, this is similar to the chasing interaction assumed in the escape-and-pursuit and cognitive flocking models [33][34][35][36], and can induce, e.g., chain migration. Indeed, our previous work [14] reported the rotating rings and spirals and dynamic assemblies with snake-like stripe shapes due to this CF term. The CIL term gives rise to an effective shortrange repulsion acting on the intrinsic polarity q (Fig.  1b left) when the coefficient α CIL is positive. This is similar to local cognitive repulsion in the animal group model [34]. CIL combined with the volume exclusion effectively causes alignment due to the pair behavior analogous to inelastic collision [13,[37][38][39]. When the coefficient α CIL is negative, this CIL term works as CAL, which is here simply assumed as the inverse direction of CIL ( Fig. 1b center), and gives rise to an effective short range attraction on q. Aggregation patterns with various dynamics were obtained due to CAL [14]. Furthermore, our model for the freely migrating case [14] exhibits polar traveling bands and homogeneous polar flocking for α CIL > 0 (repulsion), whereas aggregation for α CIL < 0 (attraction), which is similar to the collective motion due to short-range selective attractionrepulsion [40].
CIL is recognized as an important cell-cell type of communication for coherent cell migration and has been introduced in several earlier theoretical works [24,41]. The definition of CIL given by Eq. (9) is similar to that of Ref. [24] except for the following details. Our definition has cell-cell distance dependency and is additive among the cells in contact, which are different from Ref. [24]. The definition of the cells in contact is also different; our definition is based on the distance between the centroids of the two cells at basal sides, whereas the definition in Ref. [24] is based on the shared vertex at apical sides. The model in Ref. [41] uses the theoretical framework which explicitly describes the deformable cell shape, and CIL is implemented by explicitly introducing the density distribution of polarity molecules and its corresponding dynamics.
In what follows, we assume I q is infinitely large, and hence the polarity amplitude |q| is fixed and only its direction can vary (i.e., I q → ∞ and |q α | = Q).
This allows us to replace Eq. (3) by the time evolution of the angles θ j of the polarity directions, q α = (Q cos θ α , Q sin θ α ). This replacement makes the numerical calculation stable since we can avoid the large coefficient I q . Note that, since Eq. (5) does not depend on the magnitude of the polarity vector Q, we can rescale Eq. (3) and set Q = 1 without loss of generality, just by redefiningq α = q α /Q,Ĩ q = I q Q 2 andξ q α = ξ q α /Q and the coefficients in J q α accordingly. Together, Eq.
(For the brevity of notations, the tilde· on each variable has been omitted.) Here, ξ α is the angular noise defined by ξ α = −ξ q x,α sin θ α + ξ q y,α cos θ α , which satisfies ξ α (t) = 0 and Note that D −1 now characterizes the persistence time of polarity direction.

Model for apical constraint: Vertex model
For this part, we apply the cellular vertex model in two dimensions (x, y) [18,[42][43][44][45][46][47][48]. We describe the cell shapes by the 2D polygons as mentioned above, which is specified by the locations of vertices r i = (r x i , r y i ) (for i-th vertex; i = 1, 2, . . . , M) and their connectivity kl representing the edge which connects vertexes k and l. The time evolution of each vertex is given by with the mobility coefficient κ.
In the potential E, the parts relevant for the apical side dynamics are those defined within the apical side E vertex and for the apico-basal connection E A−B (Fig. 1c). (For E A−B , see Sect. 2.3.) The pure apical contribution E vertex is given by four different sub-parts as and Eq. (12) is now written as The first term of Eq. (13), E pres , represents a weak constraint giving optimum in-plane area of each cell. This is given by with the actual area A α of the αth cell, the preferred area A 0 , and the coefficient K a controlling the strength of the pressure. The area A α is calculated by the for- where the summation j∈α−cell runs for all the vertices i (numbered in the counter-clockwise order) belonging to αth cell. The second term E tens represents the potential constraining the total perimeter on cell-cell junction, given by with the actual and optimum perimeters of the α-th cell L α and L 0 , respectively, and the coefficient K p controlling the tension strength. The third term E fluc is introduced to implement the fluctuation in line tension in each cell-cell junction, for which we assume Here, ij is the length of the bond connecting vertices i and j, ij ≡ |r i − r j |. The time-evolution of junctionspecific tension γ ij is assumed to obey where η ij (t) is a white Gaussian noise satisfying η ij (t) = 0 and η ij ( The parameters τ and σ represent the relaxation time and strength of noise, respectively. The last term E wall is the potential to implement the boundaries of the system, for which we assume Here, X + and X − are the locations of the walls along x-axis, whereas Y + and Y − are the locations of the walls along y-axis. k wall,1 and k wall,2 are the coefficients specifying the strength of the confinement by the wall. (In general, the coefficients k wall,1 and k wall,2 can be different with each other, but later we assume the same value for both. The choices of these values may not cause anything significant as long as the values are enough large.) The summations i∈X + and i∈X − run over all vertices which locate at r x i > X + − δ and r x i < X − + δ, respectively, and belong to only two cells. The summations i∈Y+/− are defined similarly. The parameter δ is introduced to specify the cells located near the system boundary. In this model, with positive δ, the vertices at the boundary adhere to the boundary wall.
Cell rearrangement is implemented by junctional remodeling (JR) in the standard way. When a bond becomes shorter than a threshold length during the time evolution of vertex locations, the bond is rotated by 90 degrees around its midpoint. Simultaneously, the five edges connected to either of two vertices at the JR site are reconnected such that the T1 transition is achieved; see e.g. Refs. [18,46,49].

Connection between apical and basal sides
We naively assume that the basal side (self-propelled particle) of the α-th cell is connected to the vertices at the apical side of the same cell by a linear elastic spring. The spring is represented by the following potential function where the coefficient K means a spring constant, R α represents the centroid of the vertices of the α-th cell , and the summation i∈α−cell runs for all the vertices i (numbered in the counter-clockwise order) belonging to α-th cell.

Results
In this section, we present the results obtained by numerical simulations of the model presented in Sect. 2 [Eqs. (5), (10) and (14) with Eqs. (6), (8), (9), (11), (15), (16), (17), (18), (19) and (20)]. The parameter values which we use for the simulations in this article are as follows. Cell-cell communication strengths α CF and α CIL are the control parameters in this article and are varied over 0.0 ≤ α CF ≤ 1.0 and −2.0 ≤ α CIL ≤ 2.0. The spring coefficient of the apico-basal connection K is another control parameter, which is varied within the range from 1.0×K 0 to 6.0×K 0 with a certain K 0 defined right below to set the force unit. The shape factor P ≡ L 0 / √ A 0 is the other control parameter and is varied around the so-called solid-to-fluid transition threshold (P ∼ 3.812) [52]; 3.700 ≤ P ≤ 3.880. We set the units of length, time and force by A 0 , √ A 0 /v 0 and K 0 √ A 0 , respectively, and hence, in the equations, always put A 0 = 1, v 0 = 1 and K 0 = 1 (i.e., when K = 1.0, the migration driving force of a single cell has ability to displace the centroid of the basal side from that of the apical side over the unit length, which is roughly the diameter of a cell). The total number of cells is N = 140. The system size is 2X NA 0 along x and y, respectively. (Only in "Appendix A", the system size is √ NA 0 along both x and y.) The energy scales of the vertex model are assumed enough large, as K a = 150.0 and K p = 150.0, which sets the relaxation time scale of the apical side dynamics to be much smaller than the unit time. The other parameters are given as follows unless otherwise mentioned: β = 5.0, β w = 100.0, r = (1 + v 0 /β)r 0 with r 0 = 2.0 A 0 /π (r corresponds to the diameter of a cell at the basal side, whereas r 0 roughly corresponds to a typical cell diameter at the apical side, respectively), For the numerical integration, time is discretized with the increment dt = 0.0001 and dt = 0.001 for Fig. 7 and the other results, respectively. The maximum number of steps for each simulation run is 640, 000 (in terms of time, t = 640), 1, 600, 000 (t = 160) and 160, 000 (t = 160) for Figs. 2f, 7 and the other results, respectively.

Freely migrating cells under spatial confinement
We start with studying the case when the apical constraint is absent so that the cells can freely migrate. In other words, this is the limit when the apical side dynamics is perfectly dominated by the basal side dynamics and the apical restriction does not affect the basal side dynamics at all.  of the particles which do not induce motility-induced phase separation [50,51]. For (α CIL , α CF ) = (2.0, 0.0), the cells form the dense shell-like structure near the wall, as indicated by the broken line, which provides the inner cells with more spaces and hence the inner cells can become more motile (Fig. 2b). In contrast, in the presence of enough strong contact follow, e.g. for (α CIL , α CF ) = (1.0, 1.0), the cells accumulate in the half side of the system, which suppress the cells' motility (Fig. 2c). Even with such large α CF , when the α CIL is enough close to zero or negative, e.g., (α CIL , α CF ) = (0.0, 1.0), the coherent swirling motion is observed in the cell accumulation (Fig. 2d). In this swirling regime, cells have much higher motility than the three regimes mentioned above. (Therefore, in what follows, we call this regime/state highly-motile regime/state.) Furthermore, when we decrease α CF from this regime, we see the situation where swirling is not obvious but still stream line is visible (Fig. 2e).
In the highly-motile regime, the side in which cells are We now focus on the motility of cells depending on the cell-cell communication parameters. Due to the spatial confinement and cell-cell volume exclusion, cells in the current situation are basically much less motile than a single freely migrating cell without any constraint. We will see how and when the motility is improved back associated with the aforementioned coherent dynamics induced by cell-cell communication. We measure the mean square displacement of each cell (α-th cell) which is the rescaled MSD 10%,t=40 (αCIL,αCF) such that it gives 0 and 1 (−1) for MSD 10%,t=40 (αCIL,αCF) /MSD 10%,t=40 (0,0) = 1 and 1 ( 1), respectively. Here, the solid curve near α CIL = 0.2 in Fig. 3b, d shows the state boundaries of highlymotile and stuck regimes. The other state boundary appears at the upper edge of the region highlighted by the gray hatching.
These state boundaries surrounding the highly motile regime are determined as follows. The boundary near α CIL ∼ 0 is identified by the positive to negative jump of the Z value for increasing α CIL , as shown in Fig. 3e. To determine the other boundary, there is still an ambiguity. It may be determined as the line where the Z value suddenly drops from 0.4 to 0.2 in Fig. 3e for decreasing α CIL . However, it is difficult to identify whether the regime right below this boundary, including the case shown in Fig. 2e, is genuinely different from highly motile or disorder regime, and this is why we indicate such regime by the gray hatching in Fig. 3d. Note that the "streaming state" boundaries for the periodic boundary case given in "Appendix A," determined by the analysis of the apparent global polar order, match these boundaries. This fact supports our definitions of the state boundaries.
The other state boundary of the stuck regime is identified by the degree of asymmetric cell distribution Tini dt|c(t)| shown in Fig. 3f. See also Fig. 3e inset.
(More comprehensive discussions about the way the state boundaries are determined are given in "Appendix C".)

Epithelial cells under spatial confinement
We here investigate the case with apical constraint. We set K = 4.0, which is larger than the unit 1.0 by several factors so that the discrepancy between the centroid positions of the basal and apical sides of each cell can be ignored. To check the consistency with literature, we firstly analyze the MSD in the absence of any cell-cell communication. The result is shown in Fig. 4a for various shape index P . At t 10, the plateau is observed for P ≤ 3.812, whereas the diffusive behavior appears for P > 3.812. This result is consistent with the solid-to-fluid transition in early literature [52] (for small driving force limit in their result). Note that this transition shape index P = 3.812 matches those for the self-propelled Voronoi model at the small driving force limit [21]. This result may reflect the fact that the energy scales for apical constraints used in this article are so large (K a = K p = 150) that the spontaneous motility at the basal side does not affect the dynamics at the apical side well without coherent dynamics among cells.
Next, we study the case with cell-cell communication at the basal side. We firstly focus on the case with only contact follow α CF . Figure 4b shows the snapshots when the contact follow is incorporated, as (α CIL , α CF ) = (0.0, 1.0). As same as the free cell's case (Fig. 2d), coherent swirling dynamics is observed. Figure 4c plots the MSD curves for α CF = 1.0 for various shape index P . In this case, even when the shape factor P is smaller than the transition threshold 3.812, the MSD curve exhibits the diffusion-like tendency at t 10. These tendencies are also captured in the typical trajectories of cells shown in Fig. 4d. For α CF = 1.0, cells can still swirl under the apical constraint. Even in the solid regime, trajectories of some cells show large displacements in the presence of contact follow (e.g., P = 3.790 and α CF = 1.0; bottom-left sub-figure in Fig. 4d). This result suggests that dynamic self-organization by contact follow communication can provide cells with collective driving force enough strong to overcome the apical constraints. Remember that we here set the vertex energy scales so large, or saying in another way, the self-propelled driving force is so weak, that migration ability of cells in the absence of intercellular communication does not affect the solid-to-fluid transition threshold P .
We next perform numerical simulations for a various set of α CF and α CIL . The dependence of MSD at t = 40.0 on both α CF and α CIL is shown in Fig. 5a for a fluid phase, P = 3.880. Although the tendency of the color map is similar to that in the free cell's    (Fig. 3b), the highly motile regime is extended to the positive-α CIL side, e.g., up to α CIL ∼ 1.0 for α CF = 1.0. We also plot Z defined by Eq. 22 in Fig. 5b (for P = 3.880). The parameter windows where cells are highly motile follows the similar trends with the free case, but it is extended to higher α CIL region as mentioned above. Compare the location of the highly motile regime boundary for this case (α CIL > 0.0; see Fig. 5a, b) with the corresponding one for the free case (α CIL ∼ 0.0; see Fig. 3b, d, f). This may be because the local cell density cannot be too high due to the apical area constraint. Note that this reasoning came from the following arguments. In our results, CIL prevents the swirling dynamics as we saw for freely migrating cells (Fig. 3b, d, f). (This may further stem from the CILinduced effective alignment and the resultant emergence of the flocking state [13,14,[37][38][39]. See "Appendix A" and Ref. [14] for more details about this.) This physical fact, together with the local cell density dependence of effective CIL strength, may be the reason of this enhancement of highly motile regime due to the apical area constraint. Although such density-dependence of effective CIL strength may be caused by the distancedependent CIL coefficient in Eq. (9) in our case, the same result is expected as long as CIL strength effectively has the local density dependence for some mechanisms. It may be an interesting future direction to study how robust is this result for various assumptions for CIL. Figure 5c, d shows the maps of MSD(t = 40) and Zvalues for various P . As shown in Fig. 5c top, decreasing P just reduces the MSD value almost homogeneously over (α CIL , α CF ). The color maps of Z-value for various P given in Fig. 5d shows that the optimum P exists. It is around the solid-fluid transition point, P = 3.812. This optimality may appear because even the coherent dynamics is hindered for decreasing P when P < 3.812 (the solid regime), whereas the base motility in the absence of cell-cell communication becomes larger for increasing P when P > 3.812 (the fluid regime).
We performed the similar analysis for the case with the junctional tension fluctuation σ = 0.5, as shown in "Appendix B". In this case, the transition P becomes smaller (around P ∼ 3.790), and indeed the optimum P for the Z value becomes smaller (P ∼ 3.790) as well. "Appendix B" provides more results for the case with junctional tension fluctuation. Note also that we could not find any qualitative change of the results by this introduction of junctional tension fluctuation, besides the shift of transition P threshold.
Furthermore, in this mathematical model, we can tune the strength K of the elastic connection between apical and basal sides of each cell. Such K-dependency is briefly investigated in Fig. 6. Figure 6a, b compares MSD 10%,t=40 for K = 2.0 and 4.0. In the highly motile regime, MSD 10%,t=40 for K = 4.0 is higher than that for K = 2.0. This is a striking contrast to the tendency in the disorder regime with lower α CF , where the MSD 10%,t=40 for K = 4.0 is lower than that for K = 2.0, as shown in Fig. 6a, b, which is indeed natural because the basal dynamics is less restricted. Such characteristic K-dependency in the highly-motile regime is explicitly plotted in Fig. 6c. For K > 2.0, MSD 10%,t=40 rises for increasing K.
For further smaller K, i.e., K < 2.0, MSD 10%,t=40 seems to slightly increase for decreasing K, but this may be just because the basal side of a cell can move more easily around the centroid of its apical side. Indeed, the average distance S between the centroids of the apical and basal sides of each cell is plotted in Fig. 6c inset, showing the monotonic increase for decreasing elasticity K.   Here is a comment on the magnitude of S in this simulation. Within the range of K used in this simulation, S is near or less than 0.5 on average. This guarantees the validity of the model assuming a "phantom" spring between apical and basal sides.

Epithelial cells in a free space
In this subsection, we briefly study epithelial cell dynamics without the system boundary and demonstrate the cluster dynamics. When α CF = 0.0, characteristic dynamics of a cluster is not observed for any α CIL ; the clusters just waver around (Fig. 7a-c). When α CF is enough large, as α CF > 0.2, the cluster exhibits characteristic dynamics. For negative α CIL , the cluster tends to spin (Fig. 7d). For positive α CIL , the cluster tends to migrate or elongate, i.e. the two sides of a cluster travel to the different directions (Fig. 7e or f, respectively). These tendencies are roughly consistent with what we investigated for the case of freely migrating cells in Sect. 3.1. The coherent dynamics with apparent swirling for the freely migrating cells may correspond to the spinning in this case, whereas the stuck, or global polarly ordered flocking under a periodic boundary condition, may correspond to the cluster migration and elongation. It is apparently stochastic whether cluster migration or elongation takes place.

Summary and discussion
In this article, we extended our previous theoretical investigation on dynamic self-organization of migrating cells through cell-cell contact communication [14] to cases with spatial confinement and/or epithelial integrity, and studied what coherent dynamics are caused and how this contributes to tissue dynamics. Epithelial monolayer dynamics were numerically simulated based on a mathematical model that combined self-propelled particle models with cell-cell contact communication for basal-side dynamics and the 2D cellular vertex model for apical-side dynamics.
Our main focus was whether and how the coherent dynamics due to cell-cell contact communication can help the cells overcome the suppression of motility by the constraints mentioned above. Both with and without epithelial integrity, we identified parameter regimes in which cells exhibited highly motile coherent dynamics, including swirling in the presence of contact following. We found that such coherent dynamics can help the cells overcome constraints due to spatial confinement and epithelial integrity and enhance the cell motility. This happened even when the parameter values of the  vertex model (shape index) at the apical side were given to set the tissue in the solid regime.
By taking advantage of our model, we also investigated how the obtained coherent dynamics depend on the strength of the elastic connection between apical and basal sides K and the junctional tension fluctuation. Junctional tension fluctuation did not change the results significantly. Cell motility in the highly-motile regime, as measured by mean square displacement, depends on the strength K of the apico-basal connection. Furthermore, we demonstrated the potential behaviors of epithelial cell clusters without spatial confinement according to dynamic self-organization. Our simulation suggested that cell clusters adopt either spinning, migration, and elongation behavior.
Coherent rotational motions of migrating epithelial cells were observed in Ref. [20] by numerical simulations, under the assumption that each cell has a memory of past motile directions (with a finite memory time) and drives itself along the average of such motile directions in memory. The swirling behavior due to contact following which we observed here is distinct from such coherent rotation. Coherent rotation in Ref. [20] is a global rotation by which the migration direction of the cells follows the direction of the system boundary wall, whereas our swirling is a local circulation. It may be interesting for a future study to see what cell behaviors arise by a combination of these two different mechanisms of coherent dynamics. We expect that such study clarifies the relation between the apparent swirling dynamics observed in this article, which essentially requires contact following, and the local swirling dynamics due to CIL combined with the alignment of the polarity-to motility-directions [24].
To make the theoretical results comparable with various experiments using an epithelial monolayer, future studies will need the following technical improvements. In this article, we used much larger energy scale for the constraints of apical side area and perimeter than the basal driving force for migration such that, even in the presence of cell-cell communication, the cell area was not changed largely. However, in reality, each cell area in a monolayer can fluctuate largely and this may be important for phenomena such as collective oscillation [16]. In addition, our theoretical system comprises 140 cells, but this may be too small to compare with experimental systems required to study collective behavior. These settings were needed in this study because the illegal process of cellular vertex dynamics mentioned in the preface of Sect. 3 occurs more frequently for smaller apical energy scale and larger system size. By improving the technical details of the model and making the numerical simulation more stable, we will be able to reliably investigate wider varieties of epithelial cells dynamics.
In this article, we considered only the contact cellcell communication through the basal side. However, in real tissues, other channels of cell-cell communication may be also important. Since cells are connected by the apical junctions, they can communicate through both chemical and mechanical signaling through apical side. Moreover, the cells can secret the chemicals to the extracellular environment which might facilitate communication between cells in a distance. Although the type of cell-cell communication considered in this article is limited, our results imply that dynamic selforganization due to cell-cell communication can overcome the suppression of motility by various constraints in a tissue, which could be extended to more generic cases.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

A Freely migrating cells in a periodic boundary
In the previous work [14], we studied cellular dynamic selforganization in a periodic boundary and found that apparent global orientation of cell polarities is a useful measure to classify several states. In this appendix, we study the global orientational order of the similar system but with the small system size and parameter-value set used in the main text of this article (N = 140 cells; cf., in Ref. [14], we used N = 10, 000 − 80, 000).
The results of the numerical simulations are given in Fig. 8. Figure 8a shows the polar order parameter against αCIL and αCF with Tini = 80 and T end = 160. · means the ensemble average over 16 independent realizations. We find that as αCF is increased, the polar order arises. Time evolution of instantaneous order parameter of three representative cases (corresponding white points in Fig. 8a) are shown in Fig. 8b. For αCIL > 0.2, the stable polar order is established whereas, for αCIL < 0.2, the instantaneous order parameters fluctuates. Corresponding snapshots for αCIL = −2.0 and αCIL = 2.0 are displayed as the insets, which again shows the clear difference; the case for αCIL < 0.2 lacks the clear order which is seen for αCIL > 0.2 but exhibits converging streams. Therefore, it is expected that there is the boundary between the genuinely ordered regime and dynamic streaming regime around αCIL ∼ 0.2. In Ref. [14], by carefully studying the systems size dependency, we realized that the dependence of apparent global order parameterRp on the parameter such as αCIL is a good measure to identify such state boundary. Hence, we plotRp against αCIL in Fig. 8c and find that αCIL-dependence ofRp changes around αCIL ∼ 0.2. To see this fact explicitly, we plot ∂Rp/∂αCIL in Fig. 8d. Indeed, for αCF > 0.2, ∂Rp/∂αCIL starts to rapidly increase near αCIL ∼ 0.2 as αCIL is decreased. Collectively, we determine the state boundaries among disorder, order and streaming states as shown by the solid curves in Fig. 8a, d. B Results for the case with line-tension fluctuation at the apical cell-cell junctions   The difference is observed only in the quantitative features. In particular, the shape index P at the solid-fluid phase transition seems to be shifted to the lower P , around P ∼ 3.790, as shown in Fig. 9a. (cf. In the absence of the tension fluctuation, transition occurs at P = 3.812.) The optimum P giving the largest Z value is also shifted similarly, as shown in Fig. 10d bottom.

C How the state boundaries are determined
This appendix, together with Fig. 11, summarizes how we have divided the states in the main text of this article both for freely migrating cells (Fig. 11a-d) and the cells under apical constraints (Fig. 11e-h). In this article, we focused on the motility of the cells quantified by mean square displacements (MSD). As shown in Fig. 11a, e, g and high-  (g, h). a, e, g Mean square displacements at t = 40. b, f, h Z values. c Degree of asymmetric cell distribution c. In each, left, center and right rows show the colormap over (αCIL, αCF), dependency on αCF for various given αCIL and αCIL for various given αCF, respectively. The circles and inverse triangles highlight the state boundaries. d Sample snapshot of the regime where the outermost cells form shell-like layer lighted by the circles, we find that, when we increase the αCIL from −2.0 to +2.0, the mean square displacement at the fixed time, MSD(t = 40.0), suddenly drops at a certain αCIL for αCF ≥ 0.3; for example, at αCIL ∼ 0.2 for αCF = 1.0 (purple curve) in Fig. 11a right. These sudden drops of the motility corresponds to the shift from the state with high motility (Fig. 2d, e for freely migrating cell case) to that with low motility (Fig. 2a-c for freely migrating cell case). In this article, the former state is called the highlymotile state whereas the later state(s) called the disorder or stuck state. To further explicitly quantify these states, we used the quantity Z in Eq. (22), which is similar to measuring MSD itself at a fixed time but takes care of the following issues: -Depending on the constraint situations (freely migrating or under apical constraints) and parameter values, the MSD value at the fixed time in the absence of any intercellular communication, i.e., (αCF, αCIL) = (0, 0), varies. To avoid this, MSD value in the absence of any intercellular communication was subtracted. -The absolute value of MSD varies by orders of magnitudes depending on the situations. To avoid this, in the definition of Z, the value was rescaled such that it ranges within −1.0 and +1.0. -As shown in Fig. 3c, most of the cells are apparently less motile. In the definition of Z, we only used the cells with the top 10 percent of the MSD over all the samples.
Empirically, this definition of Z works well to classify the high-and low-motility states; As shown in Fig. 11b, f, h, around the state boundary, the sign of this Z's value flips. According to this empirical characteristics, we decided to use the sign of Z to define the boundary between the highand low-motility regimes.
Be cautious that we have another state boundary for the cases with freely migrating cells as follows. Plotting the degree of asymmetric cell distribution c exhibits the sudden changes of the slopes over αCF, indicated by the inverse triangles in Fig. 11c right. In the parameter region where αCF is lager than this boundary, we find the low MSD or negative Z as well as high value of c. This means that, in this regime, the cells are not only asymmetrically distributed but also less motile. Such qualitative difference around this boundary is consistent with the difference between the states shown in Fig. 2b, where no asymmetric distribution is observed, and Fig. 2c, where the cells are asymmetrically distributed and stuck at a system edge. In the state shown in Fig. 2b, apparently the outermost cells are forming high-density shell-like layer. When the cells form such high-density shell-like layer, the motility becomes higher due to the lower density in the bulk regions as illustrated in Fig. 11d. This is why the positive Z takes place at the bottom right region of Fig. 11b left. We excluded this positive-Z region from the above arguments to dissect the highly motile and the other states.
There are other ways to define the highly-motile states. In Ref. [24], so-called cage relative mean square displacement and the vortex density are used to define it. This definition can further dissect the global and local coherent motions. Although this article did not discuss this aspect, it is a possible future direction to investigate how the pair of contact following and contact inhibition/attraction of locomotion, which are the types of intercellular communication used in this article, can affect the locality of coherent motility of cells.