Wave Propagation of Junctional Remodeling in Collective Cell Movement of Epithelial Tissue: Numerical Simulation Study

Introduction

The morphogenesis of embryonic development involves epithelial tissue deformations (Tomer et al., 2012). Epithelial tissue consists of epithelial cells that adhere to each other through cell-cell junctions, such as adherence junctions and tight or septate junctions. These intercellular junctions are connected to the intracellular actin cables of individual cells (Harris and Tepass, 2010; Takeichi, 2014), which generate a mechanical tension that is exerted on the cell-cell junction. The mechanical forces of an individual cell are further transmitted to its neighboring cells through the junctions. Consequently, the cell monolayer achieves mechanical integrity as an epithelial tissue. While maintaining the balance of mechanical forces in the tissue, an epithelial sheet can dynamically deform to elongate and fold during development. Such epithelial deformations are often achieved through the positional rearrangements of epithelial cells (Bertet et al., 2004; Nishimura et al., 2012; Tomer et al., 2012). Such cell rearrangements in a cell sheet occur through the extension and contraction of cell-cell junctions in the apical plane, which lead to junction exchange between neighboring cells, which is called a T1 transition, while maintaining the monolayer integrity (Bertet et al., 2004; Blankenship et al., 2006; Rauzi et al., 2008). Therefore, the remodeling of the network of cell-cell junctions in the apical plane is a key mechanism by which epithelial tissue achieves a large deformation.

A typical example of a deformation driven by the junctional remodeling is convergent extension, by which a tissue elongates in one direction while shrinking in the perpendicular direction (Tada and Heisenberg, 2012; Collinet et al., 2015). During this process, the cell junctions perpendicular to the extension direction shrink and cells intercalate with each other along the direction of tissue convergence (Munro and Odell, 2002a,b). Consequently, the tissue achieves a large deformation without changing the cell size or shape and thus the internal stress. A well-studied example is the germ band extension (GBE) occurring during Drosophila early embryogenesis. During GBE, epithelial cells undergo intercalation directed along the dorso-ventral (DV) axis. A high accumulation of non-muscle myosin II (Myo-II) at the anterior and posterior cell boundaries increases the strength of the junctional contraction, causing a decrease in junctional length. These DV-oriented cell intercalations cause the tissue to narrow along the DV axis and lengthen along the anteroposterior (AP) axis, resulting in GBE. In other examples of large scale morphogenesis driven by the junctional remodeling, mediolaterally oriented cell intercalation contributes to kidney tubule elongation in Xenopus (Lienkamp et al., 2012), and polarized apical cell constriction drives neural tube invagination in the chick (Nishimura et al., 2012).

Such cell intercalation accompanied by the junctional remodeling is also a driving mechanism for the collective cell movement of epithelial tissue. During the morphogenesis of Drosophila male terminalia, the genitalia undergo a 360° clockwise rotation, which induces dextral spermiduct looping (Suzanne et al., 2010; Kuranaga et al., 2011). The Drosophila genitalia rotation (DGR) is achieved by the collective clockwise movement of surrounding epithelial cells. We previously reported that this collective cell movement is driven by polarized cell intercalation at the right oblique cell boundaries in the surrounding epithelial tissue (Sato et al., 2015a). The moving cells intercalate while remaining attached to their neighboring cells. Most of the remodeled cell boundaries form right oblique angles with the AP axis and show Myo-II accumulation. In addition, numerical simulations revealed that such diagonally polarized cell intercalation is sufficient to induce unidirectional cellular movement (Sato et al., 2015a). We also revealed that such left-right asymmetry of the cell boundary motion accompanied by AP asymmetry of the tissue is indispensable for the unidirectional movement (Sato et al., 2015b). Since epithelial cells also have the asymmetry of apico-basal polarity, the left-right asymmetry of Myo-II accumulation and resultant cell boundary motions in the planer plane can be referred to as the chirality, or the handedness, of the cells (Sato et al., 2015b).

Both GBE and DGR are induced by the cell intercalation. However, the dynamicity of cells in the tissues shows a strong contrast between the two situations. GBE is induced by a tissue deformation involving different aspect ratios of singly rearranged cells. In contrast, DGR involves the movement of a cell collective. Specifically, small cells of about 5-μm diameter move a distance of 300 μm or more in over 12 hours (Kuranaga, 2012; Sato et al., 2015a). What makes this difference in the dynamicity of the two systems? One obvious difference is the direction of cell intercalation; it is perpendicular to the direction of cell movement in GBE (Collinet et al., 2015), while diagonal in DGR (see Figure 4 in Sato et al., 2015a). Consequently, the cell intercalation events take place transiently in GBE until the tissue deformation finishes, while in DGR the cell intercalation events can continue to occur. However, this diagonal direction of cell intercalation is not sufficient to explain the collective cell movement that lasts for more than 12 hours. One factor underlying the difference in these processes is predicted to be the spatiotemporal dynamics (time order and distribution) of the cell intercalation. What determines the spatiotemporal dynamics of the cell intercalations?

A molecular signaling activity can regulate the spatiotemporal dynamics of junctional remodeling and cell rearrangement. But, it can also be organized spontaneously through mechanical processes. Here, we referred to such a spontaneous spatiotemporal organization of cell intercalation without relying on molecular signaling activity as self-organization. After a certain period of time, a series of T1 transitions gives a viscous property or plasticity to a tissue that enables large-scale deformation. In contrast, in the shorter term, elastic behavior can appear in response to the T1 transition of an individual cell that is induced by intracellular molecular signaling. When individual T1 transition is induced, mechanical forces are exerted on the surrounding cells. Such forces can be transmitted through the network of cell-cell junctions. The transmission of such a force usually ceases within the distance of a few cells (Farhadifar et al., 2007). If it further triggers another T1 transition, however, the transmission can occur over a longer distance. Consequently, the junctional remodeling and cell rearrangement could spread throughout an entire tissue in a spatiotemporally correlated way. It is therefore possible that the large scale reorganization of epithelial tissue is not only instructed by molecular signaling activity of axis information, but also involves self-organization through mechanical coupling among cells. However, such self-organized spatiotemporal dynamics of cell intercalations have not been described in our previous study. Therefore, the question that we address in this paper is whether the transmission of mechanical forces due to the cell intercalation can occur without spatiotemporal signal instruction, and how the cell rearrangements by junction remodeling are spatiotemporally organized during the collective epithelial cell movement in DGR.

In this paper, we investigated theoretically the spatiotemporal dynamics of cell movement in the ring-like epithelial tissue that surrounds the male genitalia, using a mathematical model that we described previously (Sato et al., 2015a,b). We show for the first time that the T1 transitions and cell rearrangement propagate in space in the same direction of the collective cell movement, leading to the propagation of cell velocity change in space with speed faster than the collective cell movement.

Model and Methods

Numerical Model of the Epithelial Cell Monolayer with Left–Right Cell Asymmetry

We describe the model for the rotational motion of Drosophila genital disc that we have introduced in the previous work (Sato et al., 2015a). We use a vertex cell model in two-dimensional space to simulate the in-plane motions of cells in an epithelial monolayer tissue (Nagai and Honda, 2001; Farhadifar et al., 2007). When we look at the epithelial cell monolayer from the apical side (Figure 1A), individual cells exhibit a polygonal shape. Based on this observation, the cell shapes are described as polygons with vertices and edges (Figure 1B). A polygon in the model is specified by the location of the i th vertex, r_i = (x_i, y_i), and the bond 〈kl〉 that connects vertices k and l. By considering the force balance between the potential force −∂E({r_i}_i, {γ_kl}_〈kl〉)/∂r_i and the friction force of the simplest form −η_ir_i/dt, the time evolution equation of the vertex model is given by

$\begin{array}{ll}\frac{d{\text{r}}_i}{dt}=-\frac{1}{{\eta }_i}{\frac{\partial E({\{{\text{r}}_i\}}_i,{\{{\gamma }_{kl}\}}_{\langle kl\rangle })}{\partial {\text{r}}_i}|}_{{\gamma }_{kl}=\widehat{\gamma }({\theta }_{kl})}& (1)\end{array}$

Here, the potential function E is given by

$\begin{array}{ll}E({\{r_i\}}_i,{\{{\gamma }_{kl}\}}_{\langle kl\rangle })=E_{\text{p}}+E_{\text{t}}+E_{\text{b}}.& (2)\end{array}$

FIGURE 1

Figure 1. Illustration of the mathematical model. (A) Actual cell–cell junctions visualized with E-cadherin. The epithelial cell monolayer resembles polygons whose shared sides correspond to cell–cell junctions. (B) The cell vertex model in two dimensions, modeling the polygon-like structure of the network of cell–cell junctions. The vertex index is indicated by i and j, and the bond index is given by 〈i, j〉. (C) Bond-dependent tension strength. When the bond is tilted in the clockwise direction by ${\theta }_0={\mathrm{4.html5}}^{°}$ from the AP axis, the line tension of the bond is maximum. (D) The heterogeneity in line tension strength in the model. (E) Process of junctional remodeling, called a T1 transition. Lines indicate the edges of cells. (F) T1 transition process in the numerical simulation.

The first term on the right hand side of Equation (2), E_p, is the contribution from the hydrostatic pressure within a single cell given by

$\begin{array}{ll}{E}_{\mathrm{\text{p}}}=\frac{K_a}{2}{\displaystyle \sum _{\alpha =1}^N}{(A_{\alpha }/A_0-1)}^2,& (3)\end{array}$

with the actual area A_α of the α-th cell, the preferred area A₀, and the coefficient K_a controlling the strength of the pressure. The integer N indicates the total number of cells. The second term E_t is the potential function for the total tension on cell-cell adhesions, given by

$\begin{array}{ll}{E}_{\mathrm{\text{t}}}=\frac{K_p}{2}{\displaystyle \sum _{\alpha =1}^N}{(L_{\alpha }/L_0-1)}^2+{\displaystyle \sum _{\langle ij\rangle }}{\gamma }_{ij}(t){\ell }_{ij}.& (4)\end{array}$

The first term is the contribution due to the mismatch between the actual perimeter L_α of the α-th cell and the preferred perimeter L₀, with the coefficient K_p controlling the tension strength. The second term is another contribution due to the cell's spontaneous and bond-specific transport of molecules, such as myosin II and cadherin, which exert tension with the bond-specific tension strength γ_kl (Figure 1C), which can be time dependent, as described in detail below.

The last term E_b is the potential function describing the boundaries of the system as explained below. The epithelial monolayer forms a ring-like tissue. It is flanked by two circular walls with different radii (R_out and R_in for the outer and inner walls, respectively), both of which are centered at the origin (0, 0). These walls mimic the concentric ring-like tissue structure that surrounds the Drosophila male genitalia at the larval stage (Sato et al., 2015a). We assume that additional forces are exerted on the vertices i located along the outer wall ($i\in {O}$) and the inner wall (i ∈ J), described by the potential function as

$\begin{array}{ll}{E}_{\mathrm{\text{b}}}=\frac{k_{\mathrm{\text{out}}}}{2}{\displaystyle \sum _{i\in {O}}}{(r_i-R_{\mathrm{\text{out}}})}^2+\frac{k_{\mathrm{\text{in}}}}{2}{\displaystyle \sum _{i\in J}}{(r_i-R_{\mathrm{\text{in}}})}^2,& (5)\end{array}$

where $r_i=\sqrt{x_i^2+y_i^2}$ is the distance of the i th vertex from the origin (0,0), and k_out and k_in are coefficients. We also assume that the friction coefficients η_i are different between vertices along the outer and inner walls.

The bond-specific tension γ_kl includes regulatory processes that are cell-autonomous and not described by potential functions. We thus consider that the force in Equation (1) is given as a derivative with respect to the position r_i under the condition that γ_kl is independent of r_i. After deriving the force, we consider that the tension strength γ_kl depends on the direction θ_kl of the bond 〈kl〉 around the AP axis of the tissue, as ${\gamma }_{kl}(t)=\widehat{\gamma }\left[{\theta }_{kl}({\text{r}}_k,{\text{r}}_l)\right]$. This consideration breaks the conservation of the potential E in Equation (2), and keeps the system from relaxing to the equilibrium state (Sato et al., 2015b). The specific form of γ_kl(t) in this article is given by

$\begin{array}{ll}{\gamma }_{kl}(t)={\gamma }_c\mathrm{\text{co}}{\mathrm{\text{s}}}^2({\theta }_{kl}(t)-{\theta }_0),& (6)\end{array}$

where ${\theta }_0=+45°$ and γ_c is a constant giving the maximum tension. We define the sign of the angle θ_kl in the clockwise direction around the AP axis. We also consider the case in which temporal fluctuations are present in the tension γ_kl(t), given by

$\begin{array}{ll}{\gamma }_{kl}(t)={\gamma }_c\mathrm{\text{co}}{\mathrm{\text{s}}}^2(2\pi f_{kl}t+{\delta }_{kl})\mathrm{\text{co}}{\mathrm{\text{s}}}^2({\theta }_{kl}(t)-{\theta }_0),& (7)\end{array}$

where f_kl is the frequency and δ_kl is the initial phase of the bond 〈kl〉. The frequency and initial phase are given by uniformly distributed random numbers in f_kl ∈ [0, 1] and δ_kl ∈ [0, 2π] in the presence of temporal fluctuations (Figure 1D).

When the length of a bond becomes shorter than a threshold ϵ during the time evolution according to Equation (1), the bond rotates 90 degrees around its midpoint, and the five bonds of the two vertices at the rotated bond reconnect to achieve a T1 transition (Figure 1E).

The numerical simulation is performed based on Equations (1)–(5) by the simple Euler method with time discretization of the increment dt = 0.001. The T1 transition occurs when the length ℓ_kl of the bond 〈kl〉 becomes shorter than a given length ϵ = 0.005 (see also Figure 1F). The friction constant is set to be η_i = 1 except for the outer and inner walls. We choose the units of length and time so that the radius of outer wall and the coefficient for tension are given by R_out = 1 and γ_c = 1 (η_iR_out/γ_c = 1), respectively. The friction at the outer wall is given by η_i = 100, and at the inner wall as η_i = 10. The total number of cells is given by N = 450. Additional parameter values are given by R_in = 0.5, k_out = 100, k_in = 100, K_a = 10, $A_0=(\pi R_{\mathrm{\text{out}}}^2-\pi R_{\mathrm{\text{in}}}^2)/N~0.0052$, and $L_0=\sqrt{4\pi A_0}~0.26$. The strength of the tension K_p is shown in the figures. This set of parameter values are basically comparable with our previous work (Sato et al., 2015a) except for the total number N of cells, which was 168 in the previous work. By increase of N, we are able to investigate the mechanical processes during the collective cell movement by measuring spatiotemporal distributions and correlations of the velocity fluctuations of cell motility and the T1 transition frequency.

Note that the collective cell migration of genital disc of Drosophila is considered to be driven mainly by the activity at the apical side (Sato et al., 2015a). Therefore, in our model active processes were considered for the apical processes, such as the remodeling of cell-cell junctions. For basal processes, only a passive resistance force was considered between cells and the basal extracellular matrix and cells. A collective migration of epithelial cells may also be induced by protrusive activities along the basement membrane. Motivated by such basal processes, directional driving forces of individual cells have been introduced to a cell vertex model in self-propelled Voronoi model (Li and Sun, 2014; Garcia et al., 2015; Bi et al., 2016).

Analysis of Simulation Results

In the following sections, from numerical simulations, we determine the cell angular velocity v, the elliptical cell shape anisotropy s, the frequency of T1 transition, n_T1, and the bond chirality c. The angular velocity v_α of the α-th cell is given by v_α = dψ_α(t)/dt around the origin (0,0) (genital disc center), where ψ_α is the angular position of the α-th cell around the origin, and the cell position is given by the centroid of the polygonal cell. The angular velocity v_α is positive when the cell moves in clock-wise direction. Then, v is given as the average of v_α among the cells. For the elliptical cell anisotropy s_α of the α-th cell, we first quantify the moment matrix of inertia M of the polygonal cell, and s_α is obtained by the half difference of the two eigenvalues of M. Thus, s_α measures the deviation of the cell shape from that of a circle. The variable s is given as the average of s_α among the cells. For the frequency of T1 transition, n_T1, we count the number of T1 transition per unit time. For the bond chirality c, we count the number of tilted cell-cell junctions as described previously (Sato et al., 2015a). In this article, it is given by c ≡ N_cl/N_ccl − N_ccl/N_cl, where N_cl and N_ccl are the numbers of cell-cell junctions tilted clockwise and counterclockwise, respectively, from the AP axis. When the distribution of the bond directions is uniform with respect to the AP axis, c is zero, while c is not zero when the bond directions are biased. The mathematical details of these quantities are given in the Appendix.

Results

Oscillation in Collective Cell Movement

As reported previously, the entire ring-like epithelial tissue that surrounds the genitalia rotates in the clockwise direction (Figure 2A, Supplementary Movie 1). To see the dynamics of collective cell movement quantitatively, we first measured the global angular velocity for the case without temporal fluctuation in the tension, as plotted over time in Figure 2B. Here, the global angular velocity is the average of the angular velocity of all the cells. The global angular velocity exhibited a dependence on K_p as reported previously (Sato et al., 2015b). In Figure 2C (red symbols), the global angular velocity averaged over time was plotted against K_p. The global angular velocity showed a maximum value at around K_p = 5. As K_p increased further, the state is discontinuously changed from the rotating state to a state without rotation at around K_p ≈ 15. Such an abrupt change is called subcritical transition. This subcritical characteristic of the transition was also seen in the hysteresis curves (blue and green lines).

FIGURE 2

Figure 2. Rotational motion of the genital disc in the absence of line tension fluctuation. (A) Time-lapse images of genital disc rotation. (B) Time evolution of global angular velocity for various values of K_p. (C) K_p-dependence of the global angular velocity. The angular velocity was averaged over the time window from t = 60 to 100 (red). Error bars indicate the standard deviation in the time course. Hysteresis curves for the global angular velocity were plotted for increasing and decreasing parameter values (blue and green lines, respectively).

How are the cell motions organized during the collective movement? One possibility would be that the cells would move at the same speed. Another possibility would be that they move in the clockwise direction randomly without a cell-cell correlation. Interestingly, the global angular velocities exhibited temporal oscillatory behaviors as shown in Figure 2B, indicating that some temporal organization is present in the cell movements. To examine this oscillatory behavior more closely, in addition to the global angular velocity, we plotted the time series of the T1 transition frequency, n_T1(t), the cell shape distortion s(t), and the bond chirality c(t) in Figure 3A (see Appendix for the detailed explanation of these quantities). To see the oscillatory behavior more statistically, the temporal auto-correlation functions were plotted in Figure 3B for the four quantities. The auto-correlation function quantifies the correlation of a quantity at different time points with lag time Δt. This analysis indicated that these quantities exhibit oscillation with the same oscillation period of t ≈ 2. To determine the temporal order of the events, we next studied the temporal cross-correlation functions among the four quantities, shown in Figures 3C,D. The cross-correlation function quantifies the correlation of two quantities at different time points with lag time Δt. This analysis revealed that as the bond chirality c increases, cell distortion is accumulated, as evidenced by the increase in s (green line in Figure 3D), accompanied by a decrease in the angular velocity, as indicated by the negative peak with a positive delay time in the blue line in Figure 3C. Such cell distortion can be released by T1 transition. Figure 3D (blue line) shows that the correlation between n_T1 and s reaches its maximum positive value when the delay time is zero, indicating that an accumulation of distortion accompanies the increase in the number of T1 transition. Then, the T1 transition enhances the change in cell positions, causing an increase in the angular velocity v(t). The correlation between v(t) and n_T1(t) reaches its maximum value at a small positive delay time (red line in Figure 3C).

FIGURE 3

Figure 3. Dynamics and correlation functions for K_p = 7.5 in the absence of line tension oscillation, showing oscillations of the rotation velocity, T1 transition frequency, shape deformation and chirality in bond angles and correlations between them with time delays. (A) Time evolution of the global angular velocity v, (red), the number of T1 transitions per unit time, n_T1 (green), elliptical anisotropy s (blue), and bond chirality c (purple). For these quantities of x = v, n_T1, s, or c, we plotted the scaled values as (x − 〈x〉)/σ_x with the average 〈x〉 and standard deviation σ_x. For visibility, the plots were offset by 2 (n_T1), 4 (s), and 6 (c). (B) Auto-correlations of v, n_T1, s, and c. (C) Cross-correlation functions between v and the other quantities, given by C_xv(Δt) = 〈δx(t)δv(t + Δt)〉/σ_xσ_v, where x is n_T1, s, or c (red, green, and blue, respectively). Here, the positive time difference Δt means that x precedes v. δx = x − 〈x〉. (D) Cross-correlation functions given by C_snT1(Δt) = 〈δs(t)δn_T1(t + Δt)〉/σ_sσ_nT1, C_cs(Δt) = 〈δc(t)δs(t + Δt)〉/σ_cσ_s, and C_cnT1(Δt) = 〈δc(t)δn_T1(t + Δt)〉/σ_cσ_nT1.

Spatial Propagation of Cell Rearrangement

So far, we have studied the global quantities averaged over space. To see whether the cell movement is organized spatially, we next examined the cell movement and rearrangements in a local area. In Figure 4, the local averages of the angular velocities and the T1 transition frequency, v(θ, t) and n_T1(θ, t) are shown in Figures 4A,B, respectively (the last 10 units of time are displayed). We found small patches indicating propagations of increase in both angular velocity and T1 transition frequency generated at different positions. We then scrutinized these propagation patterns in more quantitative ways as follows. We first considered the temporal autocorrelations of the local angular velocity v (Figure 4C, red line), which also exhibited an oscillatory behavior with period of about 2 units of time, consistent with the case of global velocity (Figure 3B). Then, we calculated the spatiotemporal autocorrelation function of the local angular velocity as shown in Figure 4E (a schematic explanation of spatiotemporal correlation function is shown in Supplementary Figure 1). The peak in the correlation at t = 0 and Δθ = 0 exhibited a propagation in the clockwise direction Δθ > 0 with a velocity of about 2 (rad) per unit time. Interestingly, this propagation speed was much faster than the angular velocity of cell movement, which was at most about 0.05 (rad) per unit time (Figure 2C, K_p = 5). Considering that in the present model the potential force is balanced with the frictional force that is proportional to the velocity, the propagation of velocity change indicates that the force is transmitted through the tissue as in Introduction.

FIGURE 4

Figure 4. Spatial propagations of cellular motility and T1 transition frequency. (A,B) Spatiotemporal distributions of local angular velocities (A) and frequency of T1 transitions (B) for K_p = 7.5. (C–G) Spatiotemporal correlation functions for K_p = 7.5 in the absence of oscillation of chiral line tension. A schematic explanation of spatiotemporal correlation function is shown in Supplementary Figure 1. (C) Temporal auto-correlation of the local angular velocity (red line) and the T1 transition (green line). (D) Temporal cross-correlation function. (E) Spatiotemporal auto-correlation of the local angular velocity. (F) Spatiotemporal auto-correlation of the local frequency of T1 transition. The max/min color was plotted when the color value exceeded the max/min of the color bar. (G) Spatiotemporal cross-correlation C_nT1v(Δθ, Δt) = 〈δn_T1(θ, t) δv(θ+Δθ, t+Δt)〉/σ_nT1σ_v between fluctuations of the local angular velocity and the frequency of T1 transition. θ and Δθ here are defined in the clockwise direction. (F,G) The auto-correlations (F) and the cross-correlations (G) along the time axis of Δθ = 0.

To see the propagation of the T1 transition statistically, we next plotted the temporal and spatiotemporal auto-correlation functions in Figure 4C (green line) and Figure 4F, respectively. When compared with the local velocity [Figure 4C (red line) and Figure 4E], the oscillation and propagation of the T1 transitions were not so evident. Then, the cross-correlation function clearly indicated a correlation between the local quantities (Figures 4D,G), indicating the oscillation and propagation behaviors of the T1 transition. The temporal cross-correlation function reached a maximum value with a small positive lag time (Figure 4D), and the correlation propagated spatially in the clockwise direction (Figure 4G). This analysis of cross-correlation functions indicated that after the propagation of the T1 transition passes through, the propagation of the increase in angular velocity occurs.

To see the mechanism of the global oscillation in Figure 3, we next considered the relationship between this spatial propagation of cell rearrangement and the global oscillation. For both the global and local angular velocities, the auto-correlation functions without normalization by their variances were plotted in Figures 5A,B for K_p = 2.5 and 7.5, respectively. We found that the peak correlations at Δt ~ 1.9 and 3.7 as well as the variances (Δt = 0) of the local velocity (blue lines) were about 10 times larger than those of the global velocity (red line). We then examined the K_p-dependence of the global and local velocities. In Figure 5C, for various K_p values, we plotted the coefficient of variation σ_v/〈v〉 (CV) for the global and local velocities as red crosses and blue circles, respectively. Here, σ and 〈v〉 were the standard deviation and the temporal average of velocity, respectively. As shown in Figure 2C, the global angular velocity decreased as K_p decreased below K_p = 5. The CV for the local velocity was about ten times larger than that for the global velocity. These differences in both the auto-correlation functions and CV could stem from the propagation of cell rearrangement. Changes in the local velocity showed temporal oscillatory fluctuations, which propagated spatially. Because of this propagation, the phase of oscillation in the velocity was distributed through the entire system. Therefore, the amplitude of oscillation in the global velocity was averaged out by the summation of local velocities in the entire system. The proportion of the different phases of oscillation in the entire system was time-dependent. Such time-dependence might arise due to the fact that the system size is incommensurate with the wavelength. For a sufficiently large system, the oscillation in the global velocity may become negligible. Such inconsistency was clearly seen in the spatiotemporal profile of the local velocity shown in Figure 4A, where continuous propagations were interrupted by a sudden change in the phase of oscillation.

FIGURE 5

Figure 5. Comparison of global and local behaviors, suggesting the relationship between the spatial propagation and the global oscillation of cell rearrangements. (A,B) Auto-correlation functions $C_{vv}^{*}(\Delta t)=\langle \delta v(t)\delta v(t+\Delta t)\rangle$ of the angular velocities for K_p = 7.5 (A) and K_p = 2.5 (B). (C) Standard deviation scaled by the average time courses of the angular velocity. Red symbols, for global velocity; blue symbols, local velocity. Inset shows their standard deviation.

Propagation of Cell Rearrangement Depends on the Chirality

The direction of the collective cell movement depends on the chiral property θ₀ in tension in Equation (6). Figure 6A shows that the angular velocity v and the bond chirality c exhibited their maximum and minimum values at ${\theta }_0={45}^{°}$ and ${\theta }_0={135}^{°}$, respectively, and that the T1 transitions occurred most frequently at both these angles, as reported previously (Sato et al., 2015b). The bond chirality c disappearred at ${\theta }_0=0,90,{180}^{°}$, as did the angular velocity and T1 transition. In Figures 6B–D, the spatiotemporal correlations of the local angular velocities were plotted. We found that the propagation speed of the spatiotemporal correlation in local velocity fluctuations depended on the chiral property θ₀ in tension.

FIGURE 6

Figure 6. Dependence of collective cell movement on the chirality in the absence of tension fluctuation. (A) Dependence of average quantities on the angle θ₀ in Equation (6): v (red), n_T1 (green), s (blue), and c (purple). To calculate the K_p-dependence of the angular velocity, the time window from t = 60 to 100 was used for this average. Vertical bars indicate standard deviation in the time course. (B–D) Dependency of spatiotemporal correlation functions of fluctuations in locally averaged cell angular velocity $C_{vv}(\Delta \theta ,\Delta t)=\langle \delta v(\theta ,t)\delta v(\theta +\Delta \theta ,t+\Delta t)\rangle /{\sigma }_v^2$ for K_p = 7.5 in the absence of oscillation of chiral line tension; (B) for ${\theta }_0={22.5}^{°}$, (C) for ${\theta }_0={45}^{°}$, and (D) for ${\theta }_0={67.5}^{°}$.

The Cell Rearrangement Propagates through the System with Tension Fluctuations

During the rotation of Drosophila genital disc, the contraction of the cell-cell junctions has been shown to be accompanied by temporal fluctuations (Sato et al., 2015a). The myosin II intensity at the junctions has also exhibited temporal fluctuations in an inversely correlated way. Therefore, we next asked whether the same propagation of velocity change and T1 transitions occur in the collective cell movement of epithelial cells under this condition. To consider the fluctuating contraction, we introduced a temporal fluctuation in the tension γ_kl as shown in Equation (7). The numerical simulation of the case with tension fluctuation is shown in Supplementary Movie 2. We found that an oscillation of the global angular velocities also occurred in this condition with tension fluctuation (Figure 7A), although the oscillation was more irregular than in the case without tension fluctuation (Figure 2B). As K_p increased, the global angular velocity continuously decreased without an apparent transition (Figure 7B, in contrast to the previous case without fluctuations in Figure 2C).

FIGURE 7

Figure 7. Global angular velocity in the presence of line tension fluctuation. (A) Time evolution of the global angular velocity of an entire ring-like tissue for various values of K_p. (B) K_p-dependence of the global angular velocity. To calculate the K_p-dependence, the angular velocity was averaged over the time window from t = 20 to 60. Error bars indicate standard deviation in the time course.

In Figure 8A (red line), the temporal evolution of the global angular velocity was plotted. As evident from the auto-correlation functions in Figure 8B, the rotating motion was more irregular compared to the case without tension fluctuation (Figure 3). The temporal cross-correlation between v and n_T1 (red line in Figure 8C) indicated that the angular velocity increased or decreased suddenly after the frequency of T1 transition increased or decreased, respectively. The cross-correlation between s and n_T1 reached its positive maximum value with a negative delay time (blue line in Figure 8D), suggesting that the cell distortion increased or decreased after the frequency of T1 transition increased or decreased, respectively.

FIGURE 8

Figure 8. Dynamics and correlation functions for K_p = 7.5 in the presence of line tension fluctuation, showing a weak oscillation in cell motility and the correlations with other quantities. (A) Time evolutions of the average angular velocity (red curve), the number of T1 transitions per unit time (green curve), the average s (blue curve), and the bond chirality c (purple curve) in the time window from t = 50 to 60. The quantities x = v, n_T1, s, c plotted in this figure were rescaled as (x − 〈x〉)/σ_x, respectively, with the average 〈x〉 and standard deviation σ_x. The plots are offset by 2 (n_T1), 4 (s), and 6 (c), for visibility. (B) Auto-correlations of v, n_T1, s, and c. (C) Cross-correlation functions between v and the other quantities, given by C_xv(Δt) = 〈δx(t)δv(t + Δt)〉/σ_xσ_v, where x is n_T1, s, or c (red, green, and blue, respectively). Here, the positive time difference Δt means that x precedes v. δx = x − 〈x〉. (D) Cross-correlation functions given by C_snT1(Δt) = 〈δs(t)δn_T1(t + Δt)〉/σ_sσ_nT1, C_cs(Δt) = 〈δc(t)δs(t + Δt)〉/σ_cσ_s, and C_cnT1(Δt) = 〈δc(t)δn_T1(t + Δt)〉/σ_cσ_nT1.

To examine the propagation of velocity fluctuation, the temporal auto-correlation functions of the local angular velocity v was depicted in Figure 9A; it quickly decayed with oscillation with a period of about 1.2, which agreed with the small peak in the auto-correlation function of the global angular velocity in Figure 9B (red line). In Figure 9C, although the peak at the origin also decayed quickly in space, it exhibited a propagation in both clockwise and counter-clockwise directions, in contrast to the previous case without tension fluctuation. The propagation in the counter-clockwise direction ceased in about one unit of time and 1 rad. In contrast, the propagation in the clockwise direction, seen as the yellow region, continued for more than five units of time (Δt = 5). The propagation speed was much slower than in the previous case without tension fluctuation. The spatiotemporal auto-correlation function of the local frequency of T1 transition did not indicate a clear propagation behavior (Figure 9D). However, the integral of the correlation function with respect to time showed that the correlated behavior was seen more in the clockwise direction than the other (Figure 9F). The cross-correlation between the angular velocity and the T1 transition was plotted (Figure 9E), and showed that the correlated behavior indicated by yellow was more evident in the clockwise direction. Although the propagation speed was much slower than in the case without tension fluctuation, it was still much faster than the angular velocity of cell movement (about 0.025 rad per unit time for K_p = 7.5). We conclude that the propagation of the cell rearrangement is a robust property of the epithelial cells' collective cell movement.

FIGURE 9

Figure 9. Spatiotemporal correlation functions for K_p = 7.5 in the presence of chiral line tension fluctuation, showing the spatial propagation of cell rearrangements. (A) Temporal auto-correlation function of the local angular velocity and the local frequency of T1 transitions (red and green lines, respectively). (B) Temporal cross-correlation function between the two quantities. (C,D) Spatiotemporal correlation function of the local angular velocity (C) and the local frequency of T1 transition (D). (E) Spatiotemporal correlation function between the two quantities. (F) Integrated spatiotemporal correlation $\overline{C_{n_{\mathrm{\text{T}}1}{n}_{\mathrm{\text{T}}1}}}(\Delta \theta )$ of (D) along the time difference between Δt = 1 and Δt = 5, given by $\overline{C_{n_{\mathrm{\text{T}}1}{n}_{\mathrm{\text{T}}1}}}(\Delta \theta )=\underset{1}{\overset{5}{\int }}{C}_{n_{\mathrm{\text{T}}1}{n}_{\mathrm{\text{T}}1}}(\Delta \theta ,\Delta t)d(\Delta t)$, where C is the spatiotemporal correlation function. To obtained the plot, 132 simulations were carried out with different values of the random parameters, the frequency f_kl and the initial phase δ_kl of tension fluctuation. The integration $\overline{C_{n_{\mathrm{\text{T}}1}{n}_{\mathrm{\text{T}}1}}}(\Delta \theta )$ was first performed for each simulation. Then, the means over the 132 simulations were plotted (marks) against the angle difference Δθ. Vertical bars indicate 95 percent confidence intervals of the means among the 132 simulations.

Finally, we examined the dependence of the collective cell movement on the chiral property θ₀ in tension in Equation (6) in the presence of tension fluctuation. As shown in Figure 10A, the angular velocities v, the bond chirality c, the T1 transition n_T1, and the cell shape deformation s exhibited smooth dependences on the chiral direction θ₀ with the same tendency as in the case without tension fluctuation shown in Figure 6A. In Figures 10B–D, the spatiotemporal correlation of the local angular velocities was plotted for different values of θ₀, and showed that the propagation of the fluctuation in the local velocity fluctuations occurred in the clockwise direction for these cases.

FIGURE 10

Figure 10. Dependence of collective cell movement on chirality in the presence of tension fluctuation. (A) Dependencies of average quantities on the angle θ₀ in Equation (6). (B–D) Dependence of spatiotemporal correlation functions of fluctuations in the local angular velocity $C_{vv}(\Delta \theta ,\Delta t)=\langle \delta v(\theta ,t)\delta v(\theta +\Delta \theta ,t+\Delta t)\rangle /{\sigma }_v^2$ for K_p = 7.5, (B) for ${\theta }_0={22.5}^{°}$, (C) for ${\theta }_0={45}^{°}$, and (D) for ${\theta }_0={67.5}^{°}$. See Figure 6 for details.

Discussion

In this study, we investigated the spatiotemporal self-organization of junctional remodeling and cell rearrangement in the collective cell movement of epithelial cells based on a numerical cellular vertex model. For this analysis, we focused on the rotating motion of a ring-like tissue composed of epithelial cells that exhibit a chiral contraction of cell-cell junctions. In the absence of spatiotemporal patterns in molecular signaling activity, the spatiotemporal organization of T1 transition occurred uninterruptedly in time, so that the tissue rotated in a specific direction, as we reported previously (Sato et al., 2015a). We further found in this study that the sequence of T1 transitions occurred as a traveling wave in the following way. When a cell changed its position owing to contraction and T1 transitions, it increased stress at the surrounding cells, leading to an increase in the frequency of T1 transition at the neighboring cells. Chiral contraction could bias the propagation direction by modulating this frequency in a direction-dependent manner. Consequently, the cell rearrangements propagated in a specific direction as a traveling wave, which drove the rotation of the epithelial ring-like tissue. This traveling wave behavior occurred in both the absence and presence of tension fluctuation.

Interestingly, the speed of the traveling wave of T1 transition and cell rearrangement in this mechanism was found to be much faster than the motile speed of individual cells, as shown in Figure 11 and Supplementary Movie 1. This finding indicates that in some developmental processes, morphogenetic information may be transmitted much faster and farther through the mechanical coupling reported in this article than through cell movement alone.

FIGURE 11

Figure 11. Propagation of T1 transition and cell rearrangements. T1 transition occurs between cells (indicated by colored dots) that undergo cell intercalation. The propagation of T1 transition is much faster than the cell movement (see red-colored cell). This figure was obtained from our numerical simulation in the absence of tension fluctuations (Kp = 7.5).

Our simulation result naturally arises a question as to whether the similar traveling wave behavior can be observed experimentally, in particular, the genitalia rotation in Drosophila male terminalia (Sato et al., 2015a). In the presence of tension fluctuation, the propagation of velocity change is difficult to be seen in a single time course as shown in Supplementary Movie 2, contrasting to the case without tension fluctuation as shown in Supplementary Movie 1. To see the traveling wave behaviors in the spatiotemporal correlation function, we had to take statistics for sufficiently enough number of samples. This is also the case for experimental situations, where fluctuations are inevitably present. To confirm whether the traveling wave behavior occurs in the tissue, a statistical analysis for several samples with image segmentation needs to be carried out, which remains to be a future topic for the study of tissue morphogenesis dynamics.

Another mechanism of collective epithelial cell movement is observed in Drosophila follicular epithelium (Haigo and Bilder, 2011; Cetera et al., 2014; Chen et al., 2016; Barlan et al., 2017). Elongating follicles are known to rotate circumferentially around their long axes driven by the epithelial cell migration (Haigo and Bilder, 2011). For this, the cells use the mechanism taking advantage of basal protrusion activities, which is explained in the Model and Method section. The same mechanism of collective epithelial cell movement is reported for the in-vitro monolayer sheet of MDCK cells (Serra-Picamal et al., 2012). Interestingly, a mechanical wave propagation, such as propagations of cellular velocities, rate of cell deformation and monolayer stresses, has been observed also in such a system (Serra-Picamal et al., 2012). In contrast to these examples, the Drosophila genitalia rotation relies not on the basal protrusive activity but only on the contraction activity of apical junctions (Sato et al., 2015a), as considered in the present paper. Whether a further generic principle for collective cell migrations underlies both apical and basal mechanisms is an intriguing future topic in both studies of biophysics and tissue morphogenesis.

Our model predicts that a non-uniform spatiotemporal patterning of the motile behaviors of cells including locomotion and T1 transition can be observed in the epithelial dynamics. We found that, accompanied by the transmission of an increase or decrease in the frequency of T1 transition (Figures 4F, 9D), the transmission of an increase or decrease in the local locomotion velocity of cells was observed (Figures 4E, 9C). Here we refer to the transmission of an increase or decrease in local locomotion velocity as “motility transmission.” This phenomenon may have analogies in the motion of particles, granules, or other elements in a crowded situation. A typical example is a traffic jam, where the congestion at the rear end spreads in time and space, which has been extensively studied in the field of physics (Sugiyama et al., 2008). Collective cellular motions in crowded situations and their spatiotemporal motility patterns have also been studied for human bronchial epithelial cell (HBEC) monolayers and numerical models of epithelial cells spontaneously migrating on a substrate (Garcia et al., 2015; Bi et al., 2016). A challenging future problem will be to reveal the general rules that determine the directionality and speed of motility transmissions, by comparing the dynamics of transmission and the spreading of motility in various crowded situations. The biological functions and possible applications of this phenomenon are also important future topics.

Author Contributions

TH and TS conceived the project and designed the numerical simulations. TH performed the numerical simulations. TH, EK, and TS wrote the manuscript.

Funding

This work was supported by MEXT KAKENHI Grant Number JP26114003 (EK), and the JSPS KAKENHI, grant numbers JP16K17777 (TH), JP24687027 (EK), JP16H04800 (EK), JP15KT0086 (TS).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fcell.2017.00066/full#supplementary-material

Supplementary Figure 1. Schematic explanation of spatiotemporal correlation function.

Supplementary Movie 1. Numerical simulation of the vertex model without tension fluctuations corresponding to Figure 3.

Supplementary Movie 2. Numerical simulation of the vertex model with tension fluctuations corresponding to Figure 8.

References

Aigouy, B., Farhadifar, R., Staple, D. B., Sagner, A., Röper, J.-C., Jülicher, F., et al. (2010). Cell flow reorients the axis of planar polarity in the wing epithelium of Drosophila. Cell 142, 773–786. doi: 10.1016/j.cell.2010.07.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Barlan, K., Cetera, M., and Horne-Badovinac, S. (2017). Fat2 and lar define a basally localized planar signaling system controlling collective cell migration. Dev. Cell 40, 467–477. doi: 10.1016/j.devcel.2017.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Bertet, C., Sulak, L., and Lecuit, T. (2004). Myosin-dependent junction remodelling controls planar cell intercalation and axis elongation. Nature 429, 667–671. doi: 10.1038/nature02590

PubMed Abstract | CrossRef Full Text | Google Scholar

Bi, D., Yang, X., Marchetti, M. C., and Manning, M. L. (2016). Motility-driven glass and jamming transitions in biological tissues. Phys. Rev. X 6:021011. doi: 10.1103/physrevx.6.021011

CrossRef Full Text | Google Scholar

Blankenship, J. T., Backovic, S. T., Sanny, J. S. P., Weitz, O., Zallen, J. A., and Zallen, J. A. (2006). Multicellular rosette formation links planar cell polarity to tissue morphogenesis. Dev. Cell 11, 459–470. doi: 10.1016/j.devcel.2006.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Cetera, M., Ramirez-San Juan, G. R., Oakes, P. W., Lewellyn, L., Fairchild, M. J., Tanentzapf, G., et al. (2014). Epithelial rotation promotes the global alignment of contractile actin bundles during Drosophila egg chamber elongation. Nat. Commun. 5:5511. doi: 10.1038/ncomms6511

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, D.-Y., Lipari, K. R., Dehghan, Y., Streichan, S. J., and Bilder, D. (2016). Symmetry breaking in an edgeless epithelium by Fat2-regulated microtubule polarity. Cell Rep. 15, 1125–1133. doi: 10.1016/j.celrep.2016.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Collinet, C., Rauzi, M., Rauzi, M., Lenne, P.-F., and Lecuit, T. (2015). Local and tissue-scale forces drive oriented junction growth during tissue extension. Nat. Cell Biol. 17, 1247–1258. doi: 10.1038/ncb3226

PubMed Abstract | CrossRef Full Text | Google Scholar

Farhadifar, R., Röper, J.-C., Aigouy, B., Eaton, S., and Jülicher, F. (2007). The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing. Curr. Biol. 17, 2095–2104. doi: 10.1016/j.cub.2007.11.049

PubMed Abstract | CrossRef Full Text | Google Scholar

Garcia, S., Hannezo, E., Elgeti, J., Joanny, J.-F., Silberzan, P., and Gov, N. S. (2015). Physics of active jamming during collective cellular motion in a monolayer. Proc. Natl. Acad. Sci. U.S.A. 112, 15314–15319. doi: 10.1073/pnas.1510973112

PubMed Abstract | CrossRef Full Text | Google Scholar

Haigo, S. L., and Bilder, D. (2011). Global tissue revolutions in a morphogenetic movement controlling elongation. Science 331, 1071–1074. doi: 10.1126/science.1199424

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, T. J. C., and Tepass, U. (2010). Adherens junctions: from molecules to morphogenesis. Nat. Rev. Mol. Cell Biol. 11, 502–514. doi: 10.1038/nrm2927

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuranaga, E. (2012). Beyond apoptosis: caspase regulatory mechanisms and functions in vivo. Genes Cells 17, 83–97. doi: 10.1111/j.1365-2443.2011.01579.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuranaga, E., Matsunuma, T., Kanuka, H., Takemoto, K., Koto, A., Kimura, K.-I., et al. (2011). Apoptosis controls the speed of looping morphogenesis in Drosophila male terminalia. Development 138, 1493–1499. doi: 10.1242/dev.058958

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, B., and Sun, S. X. (2014). Coherent motions in confluent cell monolayer sheets. Biophys. J. 107, 1532–1541. doi: 10.1016/j.bpj.2014.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Lienkamp, S. S., Liu, K., Karner, C. M., Carroll, T. J., Ronneberger, O., Wallingford, J. B., et al. (2012). Vertebrate kidney tubules elongate using a planar cell polarity–dependent, rosette-based mechanism of convergent extension. Nat. Genet. 44, 1382–1387. doi: 10.1038/ng.2452

PubMed Abstract | CrossRef Full Text | Google Scholar

Munro, E. M., and Odell, G. (2002a). Morphogenetic pattern formation during ascidian notochord formation is regulative and highly robust. Development 129, 1–12.

PubMed Abstract | Google Scholar

Munro, E. M., and Odell, G. M. (2002b). Polarized basolateral cell motility underlies invagination and convergent extension of the ascidian notochord. Development 129, 13–24.

PubMed Abstract | Google Scholar

Nagai, T., and Honda, H. (2001). A dynamic cell model for the formation of epithelial tissues. Philos. Mag. B 81, 699–719. doi: 10.1080/13642810108205772

CrossRef Full Text | Google Scholar

Nishimura, T., Honda, H., and Takeichi, M. (2012). Planar cell polarity links axes of spatial dynamics in neural-tube closure. Cell 149, 1084–1097. doi: 10.1016/j.cell.2012.04.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Rauzi, M., Rauzi, M., Verant, P., Lecuit, T., and Lenne, P.-F. (2008). Nature and anisotropy of cortical forces orienting Drosophila tissue morphogenesis. Nat. Cell Biol. 10, 1401–1410. doi: 10.1038/ncb1798

PubMed Abstract | CrossRef Full Text | Google Scholar

Sato, K., Hiraiwa, T., Maekawa, E., Isomura, A., Shibata, T., and Kuranaga, E. (2015a). Left-right asymmetric cell intercalation drives directional collective cell movement in epithelial morphogenesis. Nat. Commun. 6:10074. doi: 10.1038/ncomms10074

PubMed Abstract | CrossRef Full Text | Google Scholar

Sato, K., Hiraiwa, T., and Shibata, T. (2015b). Cell chirality induces collective cell migration in epithelial sheets. Phys. Rev. Lett. 115:188102. doi: 10.1103/PhysRevLett.115.188102

PubMed Abstract | CrossRef Full Text | Google Scholar

Serra-Picamal, X., Trepat, X., Conte, V., Anon, E., Vincent, R., Tambe, D. T., et al. (2012). Mechanical waves during tissue expansion. Nat. Phys. 8, 628–634. doi: 10.1038/nphys2355

CrossRef Full Text | Google Scholar

Sugiyama, Y., Fukui, M., Kikuchi, M., Hasebe, K., Nakayama, A., Nishinari, K., et al. (2008). Traffic jams without bottlenecks—experimental evidence for the physical mechanism of the formation of a jam. New J. Phys. 10:033001. doi: 10.1088/1367-2630/10/3/033001

CrossRef Full Text | Google Scholar

Suzanne, M., Petzoldt, A. G., Spéder, P., Coutelis, J.-B., Steller, H., and Noselli, S. (2010). Coupling of apoptosis and L/R patterning controls stepwise organ looping. Curr. Biol. 20, 1773–1778. doi: 10.1016/j.cub.2010.08.056

PubMed Abstract | CrossRef Full Text | Google Scholar

Tada, M., and Heisenberg, C. P. (2012). Convergent extension: using collective cell migration and cell intercalation to shape embryos. Development 139, 3897–3904. doi: 10.1242/dev.073007

PubMed Abstract | CrossRef Full Text | Google Scholar

Takeichi, M. (2014). Dynamic contacts: rearranging adherens junctions to drive epithelial remodelling. Nat. Rev. Mol. Cell Biol. 15, 397–410. doi: 10.1038/nrm3802

PubMed Abstract | CrossRef Full Text | Google Scholar

Tomer, R., Khairy, K., Amat, F., and Keller, P. J. (2012). Quantitative high-speed imaging of entire developing embryos with simultaneous multiview light-sheet microscopy. Nat. Meth. 9, 755–763. doi: 10.1038/nmeth.2062

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendix

We calculated various quantities for our analyses as follows.

Average cell angular velocity v(t) = 〈d_ψα/dt〉_global among the entire tissue. Cell angular velocity is defined in this article as the increase in the angle of the α-th cell location in the clock-wise direction in unit time dψ_α/dt around the center of the ring-like tissue. ψ_α is given by the angle between an arbitrary direction and the direction of the center of the cell seen from the center of the disc, i.e., through the relationships $\mathrm{\text{cos}}(-{\psi }_{\alpha })=X/\sqrt{X^2+Y^2}$ and $\mathrm{\text{sin}}(-{\psi }_{\alpha })=Y/\sqrt{X^2+Y^2}$ with the position vector (X, Y) of the center of mass of the cell from the origin (0, 0) at the disc center. We calculate the center of mass of the cell consisting of vertices i = 1, 2, ⋯ , I designated in the counter-clockwise direction by

$\begin{array}{ll}X=\frac{1}{6A}{\displaystyle \sum _{i=1}^I}({\text{r}}_i\times {\text{r}}_{i+1})(x_i+x_{i+1}),& (\text{A}1)\end{array}$

and

$\begin{array}{ll}Y=\frac{1}{6A}{\displaystyle \sum _{i=1}^I}(-{\text{r}}_i\times {\text{r}}_{i+1})(y_i+y_{i+1}),& (\text{A}2)\end{array}$

with the position vector r_i = (x_i, y_i) of the i-th vertex from an arbitrary position and the cell area A. The cell area is calculated by

$\begin{array}{ll}A=\frac{1}{2}{\displaystyle \sum _{i=1}^I}({\text{r}}_i\times {\text{r}}_{i+1}).& (\text{A}3)\end{array}$

The i = (I + 1)-th vertex is identified with the i = 1-th vertex. The symbol × indicates a cross product in two dimensions, i.e., (a₁, a₂) × (b₁, b₂) = a₁b₂ − b₁a₂. The average 〈·〉 is performed among all the cells.

Local average θ(t) = 〈d_ψα/dt〉_local(θ) of the cell angular velocity. We assign the angle θ of an arbitrary position (x, y) into J grids of the width 2π/J, as 2π(j − 1/2)/J ≤ θ < 2π(j + 1/2)/J with j = 1, 2, ⋯ , J. The grid number is set at J = 50. The angle θ is defined by $\mathrm{\text{cos}}(-\theta )=x/\sqrt{x^2+y^2}$ and $\mathrm{\text{sin}}(-\theta )=y/\sqrt{x^2+y^2}$. Furthermore, we take the average 〈·〉_Local of the angular velocity among the cells whose center of mass (X, Y) is located within each grid, 2π(j − 1/2)/J ≤ ψ_α(X, Y) < 2π(j + 1/2)/J, at a later time step to define the velocity. Here, the angle ψ_α of the α-th cell location direction is defined as above.

Total number n_T1(t) of T1 transition events that occur between two adjacent time steps. We define n_T1(t) by counting the total number of T1 transitions that occur between two adjacent time steps and dividing it by the step size dt of time.

Number density n_T1(θ, t) of T1 transitions at time t. We assign the spatial angle into grids as described above and count the number of T1 transitions that occur within each grid, 2π(j − 1/2)/J ≤ ϕ < 2π(j + 1/2)/J in the time interval from t to t + dt. The grid number is set as J = 50 as above. The angle ϕ of the bond location is defined by $\mathrm{\text{cos}}(-\varphi )=X/\sqrt{X^2+Y^2}$ and $\mathrm{\text{sin}}(-\varphi )=Y/\sqrt{X^2+Y^2}$ of the midpoint (X, Y) of the bond being examined. The number density n_T1(θ, t) is given by dividing the number of T1 transitions by the time step size dt.

Elliptical anisotropy s(t) in the shape of each cell. We define the amplitude s of the elliptical deformation in the following way. First, we calculate the moment of inertia M for a uniformly thin plate, the mass density of which per area is scaled to be 1, taking the same shape as the cell being examined (Aigouy et al., 2010).

$\begin{array}{l}M=\left(\begin{array}{cc}{M}_{xx}& M_{xy}\\ M_{yx}& M_{yy}\end{array}\right)\\ =\left(\begin{array}{cc}{\displaystyle \underset{V}{\int }{d}^2\text{r}{(\Delta y)}^2}& -{\displaystyle \underset{V}{\int }{d}^2\text{r}\Delta x\Delta y}\\ -{\displaystyle \underset{V}{\int }{d}^2\text{r}\Delta x\Delta y}& {\displaystyle \underset{V}{\int }{d}^2\text{r}{(\Delta x)}^2}\end{array}\right).& (\text{A}4)\end{array}$

Here, (Δx, Δy) is the position vector of each point from the center of mass of the plate (cell), d² r = d(Δx)d(Δy), and the integral $\underset{V}{\int }$ runs for the entire plate. This tensor M is calculated through the following formulae:

$\begin{array}{ll}{M}_{xx}=\frac{1}{12}{\displaystyle \sum _{i=1}^I}(\Delta {\text{r}}_i\times \Delta {\text{r}}_{i+1})(\Delta y_i^2+\Delta y_i\Delta y_{i+1}+\Delta y_{i+1}^2),& (\text{A}5)\end{array}$

$\begin{array}{l}{M}_{xy}=M_{yx}=-\frac{1}{24}{\displaystyle \sum _{i=1}^I(\Delta {\text{r}}_i\times \Delta {\text{r}}_{i+1})[\Delta x_i(2\Delta y_i+\Delta y_{i+1})}\\ + \Delta x_{i+1}(\Delta y_i+2\Delta y_{i+1})],\end{array}$

and

$\begin{array}{ll}{M}_{yy}=\frac{1}{12}{\displaystyle \sum _{i=1}^I}(\Delta {\text{r}}_i\times \Delta {\text{r}}_{i+1})(\Delta x_i^2+\Delta x_i\Delta x_{i+1}+\Delta x_{i+1}^2),& (\text{A}6)\end{array}$

with the position vector Δr_i = (Δx_i, Δy_i) of each vertex (i-th vertex) from the center of mass of the plate (cell). Here, the cell is assumed to consist of the vertices i = 1, 2, ⋯ , I designated in counter-clockwise order of the position seen from the cell's center, and the i = 1-th vertex is identified with the i = (I+1)-th vertex. The symbol × indicates a cross product in two dimensions, i.e., (a₁, a₂) × (b₁, b₂) = a₁b₂ − b₁a₂. From this analysis, we define the elliptical cell anisotropy s_α of α-th cell as the half difference of two eigenvalues of the tensor M, which is given by

$\begin{array}{ll}{s}_{\alpha }=\sqrt{{(M_{xx}-M_{yy})}^2+4M_{xy}{M}_{yx}}.& (\text{A}7)\end{array}$

This quantity s_α is 0 when the plate is set to be a circle, and when the plate has an elliptically deformed shape from the circle, s_α has non-zero value. Therefore, we use s_α to quantify how the cell shape is elliptically deformed. Finally, the variable s is given as the average of s_α among the cells.

Bond chirality c(t). At each time step, we count the numbers N_cl(t) and N_ccl(t) of cell-cell junctions tilted clockwise and counterclockwise from the AP axis, respectively. We then define the bond chirality c(t) by c(t) ≡ N_cl(t)/N_ccl(t) − N_ccl(t)/N_cl(t). If the system has no chirality, N_cl(t)/N_ccl(t) = N_ccl(t)/N_cl(t) with the statistical average 〈·〉, and hence 〈c(t)〉 = 0.