Collective cell migration of epithelial cells driven by chiral torque generation

Various multicellular tissues show chiral morphology. Experimental studies have shown this can originate from cell chirality. However, no theory has been proposed to connect the cellular chiral torque and multicellular chiral morphogenesis. We propose a model of confluent tissue dynamics with cellular chiral torque. We found that cells migrate unidirectionally under a gradient of cellular chiral torque. While the migration speed varies depending on the tissue's mechanical parameters, it is scaled solely by a structural order parameter for liquid-to-solid transition in confluent tissues.

The establishment of left-right (LR) asymmetry in tissues and organs is an intriguing event in development, which involves the coordinated activity of cells and molecules [1][2][3][4]. Since proteins as the basic components of biological systems are chiral molecules, the breaking of LR symmetry can be a collective property orchestrated by these chiral molecular components [5,6].
For the LR asymmetry in tissue morphogenesis, it has been shown that in Drosophila, the embryonic hindgut twists unidirectionally [7], and the male genitalia undergoes a clockwise rotation when looked from the apical side [8,9]. These are epithelial tissues, which consist of epithelial cells that adhere to each other through cell-cell junctions [10]. During the hindgut twisting, the hindgut epithelial cells exhibit chiral properties in their shape and other properties [7]. In the genitalia rotation, the surrounding epithelial cells move collectively in the clockwise direction driven by LR asymmetric positional rearrangements of cells, which rotates the genital disk [9,11]. Even more generally, such chiral cell rearrangements can induce unidirectional cell migration within an epithelial tissue [12]. The above examples suggest that a collective behavior of cellular chirality gives rise to the tissue LR asymmetry.
At the single cell level, several types of cells have been reported to exhibit chiral dynamics. Examples include the extension of neurites in cultured nerve cell [13], the chiral movement of C. elegans cell cortex [14], nuclear rotation in melanophores of zebrafish [15], and actin cytoskeleton swirling in human foreskin fibroblasts (HFFs) cultured on a micro pattern [16]. Remarkably, in zebrafish melanophores and the HFFs on the micro pattern, the single cells can cell-autonomously generate chiral torques in such a way that the apical side (top side) of cells exhibits rotations with respect to the basal side (bottom side) which adheres to substrates. Although a variety of specific mechanisms generate chiral torque at the cellular level, many of these chiral properties are governed by active chiral processes of actomyosin cytoskeleton [17][18][19].
Since actomyosin is a ubiquitous component of eukary-otic cells, it is natural to consider that chiral torque generation is not a feature restricted to the specific cells mentioned above. In particular, we consider a situation where chiral torque is generated by individual single epithelial cells in a tissue, motivated by the LR symmetry breaking of the epithelial tissues [7][8][9], and by the fact that single cells can generate chiral torque [13][14][15][16]. Based on the cell vertex model (CVM) [20], we propose a theoretical model of dynamics of an epithelial tissue with a chiral torque generated by individual single cells. Then we investigate the tissue scale dynamical emergent properties.
To model the dynamics of an epithelial tissue with chiral torques generated by individual single cells, we use a two dimensional (2D) CVM. In the CVM, cells are described by polygons with vertices and edges. The position of the ith vertex is represented by r i . The force acting on each vertex is given by the derivative −∂E({ r i })/∂ r i of a potential function E({ r i }), given by Here, the first term on the right hand side describes the area elasticity with the area A α of cell α, the elastic modulus K, and the preferred area A 0 . The second term defines the perimeter elasticity with the perimeter length P α of cell α, the elastic modulus K p and the preferred perimeter length P 0 . We assume for simplicity that these cellular mechanical properties K, K p , A 0 , P 0 are spatially homogeneous. N is the total number of cells. The third term is introduced to explicitly describe the fluctuation in the line tension ∆Λ ij (t). Here, ℓ ij is the length of a cell edge between ith and jth vertices. We introduce a fluctuating tension ∆Λ ij (t) as a colored Gaussian noise with ∆Λ ij (t) = 0 and ∆Λ ij (t 1 )∆Λ kl (t 2 ) = δ ik δ jl σ 2 e −|t1−t2|/τ . Such timecorrelated fluctuation of line tension is reported in [21]. We set τ = 1 in this letter for simplicity. Hereafter, we choose the units of length and forces such that the elastic modulus K p and the preferred cell area A 0 are unity, i.e., K p = 1 and A 0 = 1. With these units, P 0 is a so-called target shape index, which is a ratio between perimeter and square root of area [22,23]. If a single cell tends to be circular, P 0 = 2 √ π ∼ 3.54. We set P 0 = 3.54 and K = 10 to avoid large cellular shape deformation unless otherwise noted.
In the case of torque generation by isolated single cells, as reported in [15,16], the apical surface rotates unidirectionally with respect to the basal substrate, which is driven by the torque generated by actomyosin network. In the epithelial tissue, we also consider that a torque force is generated at the apical side in a specific direction with respect to the basal side, which adheres to the extracellular matrix (ECM). To represent such a torque in the framework of the 2D CVM, we introduce the torque force around the cell center as shown in Fig. 1(a). The simplest form may be given by where ν α is the coefficient of the torque force generated by the cell α, r α g is the area centroid of cell α, and n is a unit normal vector from the basal to apical sides.
The time evolution equation for r i is obtained by considering the force balance between the frictional force ηd r i /dt with friction constant η, the potential force −∂E({ r i })/∂ r i derived from Eq.(1), and the torque force T i given by Eq.(2) as follows; When the length of a cell edge falls below a threshold l th = 0.03 during a time-evolution according to Eq.(3), a T1 transition is performed by flipping the edge by 90 • [20]. We consider a simple configuration of a rectangular epithelial sheet with size L x and L y in the x-and ydirection, respectively, as shown in Fig. 1(b). For the x-direction, we apply the periodic boundary condition at x = 0 and x = L x . The rectangular area size is equal to the number of cells so that the average cell area is set to beĀ α = 1. This configuration corresponds to the tube geometry, a biologically ubiquitous structure of organs, which exhibits chiral morphogenesis as twisting of a heart tube and hindgut. We simply consider that the cells are attached to the boundary. We hence fix the ycoordinates of vertices on the bottom and top boundaries to y = 0 and y = L y , respectively. η is set to 1 unless otherwise noted. We calculate the time-evolution Eq. (3) with the time step ∆t = 0.01. We prepare initial cellular configurations packed with regular hexagonal cells to calculate the dynamics with σ = 0 and ν α = 0 to relax the system, and then set σ and ν α to the target values. The numbers of cells in the initial hexagonal configuration in the x-and y-directions are respectively denoted as N x and N y , and we set N x = 10 for all the simulation.
We first investigated the role of chiral torque generation on the dynamics of cells when the strength ν α of chiral torque generation is spatially homogeneous. As shown in Fig. 2(a), we found that the cells migrate bidirectionally when the fluctuation in line tension is present (ν α = 0.2, σ = 0.3, N y = 6). With negative ν α , the cellular behavior was reversed, confirming that the chiral torque generation is the driving force of the cellular migration. For sufficiently large system (N y = 40), the flow profile rapidly decays near the boundary, and the bulk velocity vanishes (Fig. 2(b)). This indicates that the torque force and the potential force are balanced in the bulk, while not balanced near the boundary. On the top and bottom boundaries, the torque forces on the vertices are exerted in the rightward and leftward direction, respectively, as depicted in Fig. 2(c), leading to the bidirectional cellular flow at the boundaries.
In the bulk, to understand how the torque force and the potential force are balanced, we consider a meanfield model where the torque force is exerted on a vertex O surrounded by three regular hexagonal cells as shown in Fig. 1(a), without the line tension fluctuation. When ν α is homogeneous, we readily find that the torque force exerted on the vertex O vanishes due to the 3-fold symmetry. Therefore, even though the torque force is present in the bulk due to cellular deformation, it should be so small that the deformation of cells can readily balance it.
We note that without the line tension fluctuation (σ = 0), all the cells are just deformed in the asymmetric fashion without continuous cellular flow (See Supplemental Material [24]). Hence, the cell rearrangements induced by the stochastic fluctuation of the line tension promote the continuous cellular flow.
We next consider a condition in which collective cell migration appears in a defined direction, induced by chiral torque generation. According to a symmetry argument [12], the symmetry along the y-axis has to be broken to achieve a unidirectional cellular movement along the x-axis. In this letter, we consider a situation where chiral torque generation depends on the position of the cells along the y-axis. We simply assume a linear form given by ν α = −λ(y α g − L y ), where y α g is the y-coordinate of the center of cell α. Here, we impose a boundary condition in which the friction coefficients on the top and bottom boundaries are considerably higher (η = 10000) to avoid the effect of the boundary torque force as we discussed above (Fig. 2(c)). With the torque gradient, we found that cells migrate unidirectionally in the direction perpendicular to that of the torque gradient ( Fig. 3(a)). A time-averaged flow profile is shown in Fig. 3(b) (λ = 0.01, σ = 0.3, N y = 40). With negative λ, the direction of the cellular migration is reversed. As λ increases, the steady-state cellular velocity V ave in the x-axis increases almost linearly (Fig. 3(c)). These results confirm that the chiral torque gradient is the driving force for this unidirectional collective cell migration. We also found that without the line tension fluctuation, the cells only deform without continuous flow (See Supplemental Material [24]). Hence, the cell rearrangements induced by the stochastic fluctuation in the line tension are necessary for the continuous cellular flow.
In contrast to the case in which the strength of torque strength is homogeneous, the flow profile shown in Fig. 3(b) indicates that the cells can migrate in the bulk region under the torque gradient. To see the mechanism of the unidirectional cellular migration, we consider the mean field model of three cells (Fig. 4(a)). We set the position of the target vertex O as the origin of the coordinate axes, and impose a linear torque gradient λ along the y-direction. Without loss of generality, we define the center of each surrounding cell as r α g = l(cos (2πα/3 + β), sin (2πα/3 + β)), where l = 2 √ 3/3 is the edge length for the hexagonal regular cell with area 1, β is an arbitrary constant and the cell number α is indexed from 1 to 3. In this configuration, since the force derived from the potential function E vanishes in total owing to the 3-fold symmetry, only torque forces from the 3 cells contribute to the force exerted on the vertex O. After a straightforward calculation, we obtain the force f = (3λl 2 /2, 0). Consequently, the mean field model predicts that, when the cells are packed in a regular hexagonal pattern without any boundary, the cells migrate with the speed v theory = 3λl 2 /2η.
We tested this prediction numerically using the model without fluctuation in line tension. We apply a small value of P 0 = 2.90 to make the cell shape close to regular hexagon. To suppress the boundary effect, we set the friction coefficients of the vertices on both of the bottom and top boundaries equal to those in the bulk, and set the torque force on the vertices on the boundaries to zero. As shown in Fig. 4(b), we found that, by increasing the system size N y , the average cell speed V ave approaches to v theory for each λ, confirming that the spatial gradient of torque generated by the cells drives the cellular migration. The deviation from the theoretical value even for the large system size is owing to the inevitable cellular deformation from the regular hexagon assumed in the mean field model. In the inset of Fig. 4(b), we confirmed the cellular deformation by calculating the av- erage q α of the cell shape q α = P α / √ A α , which equals to q hex = 2 √ 2 4 √ 3 ∼ 3.72, when the cells are regular hexagon.
From the above analysis, we found that the cell shape affects the cellular migration velocity. The cell shape can be modulated by altering the target shape index P 0 or the noise strength σ. Therefore, we investigated how the migration speed V ave depends on P 0 and σ (inset of Fig. 5). Here, we set N y = 40, λ = 0.01 and applied the same boundary conditions that were used in Fig. 3. We found that an increase of either P 0 or σ enhances the migration speed (inset of Fig. 5). It has been reported that an increase of P 0 turns the energy barrier for T1 transition to be lowered to almost zero above a critical value P * 0 ∼ 3.81, which is liquid-to-solid transition in confluent tissues [23]. Thus, as P 0 increases, cell rearrangements occur more frequently. The energy barrier is also overcome by increasing the line tension fluctuation σ, leading to an increase of the cell rearrangements frequency. Hence, the migration velocity increases as either P 0 or σ increases (inset of Fig. 5).
To see the influence of P 0 and σ in a unified way, we pay attention to the cell shape q α obtained from the numerical simulations. In Fig. 5, we plotted the migration velocity V bulk in the bulk layers against the cell shape q α BL averaged in the boundary layers, and found that the velocity curves in the inset of Fig. 5 are surprisingly collapsed in a single line. Hence, q α BL is an excellent indicator of the migration velocity in our model. Here, we define the boundary and bulk layers from the flow profile, and exclude the cells on the bottom and top boundaries for calculation of q α BL , since the cells are strongly deformed by the flat boundaries (See Supplemental Material [24]).
We discuss how the bulk velocity V bulk is uniquely determined by the cell shape q α BL in Fig. 5. The cellular velocity should be determined by the balance between the bulk torque force and how easily cells rearrange over the energy barrier of T1 transition in the boundary layers. In the present model, since the torque force depends on the cell shape, the bulk torque force should be determined by the cell shape q α BU averaged in the bulk. The cell shape q α has been also suggested as an indicator of how easily cells rearrange over the energy barrier of T1 transition [25,26]. Hence, the energy barrier of T1 transition in the boundary layers depends on q α BL . Consequently, considering that q α BL ≈ q α BU is satisfied (See Supplemental Material [24]), and hence both bulk torque force and energy barrier of T1 transition depend on q α BL , q α BL uniquely determines the migration velocity.
In this letter, we have reported that a dynamical chiral property at the single-cell level, which is cellular chiral torque generation, can induce the collective cellular migration. Our model predicts that, when the strength of torque is spatially homogeneous, the torque forces in the bulk almost balanced, and those generated by the cells at the boundaries drive the bidirectional cellular migration. Another prediction is that, under the gradient of chiral torque strength, the cells can migrate unidirectionally and perpendicularly to the gradient driven by the torque force generated in the bulk. Although the mechanism of the torque generation at the single-cell level is not fully revealed experimentally, previous studies reported that a combination of cytoskeleton and motor proteins generates the cellular torque. Activity of such cytoskeleton and motor proteins is regulated by various biochemical pathways such as Rho signaling pathway [27]. Hence, the spatial gradient of regulatory molecules should produce a gradient of the strength of the cellular torque. In in vivo systems, an epithelial tissue is attached to other different tissues, and hence such attached tissues probably emit biochemical signals to generate the concentration gradient of molecules regulating the strength of the cellular torque. Since both the cellular chiral torque generation and the gradient of the cellular torque are considered to be ubiquitous in biological systems, we expect that the mechanism of the LR symmetry breaking in tissue dynamics we have proposed plays an essential role in the chiral collective cellular movements, such as a twist of epithelial tube [7] and a unidirectional epithelial cellular flow [9], during development.

I. DEFINITION OF THE AREA CENTROID OF A CELL
We define the cell center of cell α using the area centroid as r α g = 1 6Aα Nα µ=1 ( r α µ + r α µ+1 )( r α µ × r α µ+1 ) · n. Here, µ is the index of the N α vertices around the cell α, indexed in the counterclockwise direction. The position vector of each vertex is represented by r α µ and r α Nα+1 = r α 1 .

II. PREPARATION OF THE INITIAL CELLULAR CONFIGURATION
We prepared the initial condition as shown in Fig. S1(a). In the initial condition, the cells in the bulk are set to regular hexagonal shape, while the cells at the boundaries, which are labeled with red dots in Fig. S1(a), are set to half of the regular hexagonal cells. The average cell areaĀ is set to 1. Then, we calculate the dynamics with σ = 0 and ν α = 0 to relax the system for 50 time units, and then we set σ and ν α to the target values. Fig. S1(b) is the configuration after the relaxation (P 0 = 3.54).

IV. DEFINITION OF THE BOUNDARY LAYERS
We systematically defined the boundary layers near either the bottom or the top boundaries. We discretized the area into N y layers from the 0th layer to the (N y − 1)th layer. The ith layer is defined as the layer ranging from y = iL y /N y to y = (i + 1)L y /N y (i = 0 ∼ (N y − 1)). The boundary layers are defined as the 0th∼N b th and N t th∼(N y − 1)th layers, respectively for those near the bottom and top boundaries.
We calculated the average velocity of each layer by averaging the velocity of cells in the layer in time, and drew the velocity profiles as shown in Fig. S3(a) (N y = 40, λ = 0.01, P 0 = 3.33). Using the velocity profile, we defined the boundary layers via the following procedure.
1. We calculated the approximate bulk velocity V ′ bulk and the approximate standard deviation σ ′ bulk of the bulk velocity by averaging the velocity of the 15th ∼ 24th layers.
2. We smoothed the velocity profile to obtain v = v sm (i) by a simple moving average, where the mean is taken from two data on either side of a central value.
3. We determined N b as the maximum integer i < 15 which satisfies either v sm (i) < V ′ bulk − 3σ ′ bulk or v sm (i) > V ′ bulk +3σ ′ bulk . Also, N t is determined as the minimum integer i > 24 which satisfies either v sm (i) < V ′ bulk −3σ ′ bulk or v sm (i) > V ′ bulk + 3σ ′ bulk .
In Fig. S3(a), we show the N b th and N t th layers, determined by the above procedure, with the square markers on examples of velocity profiles (N y = 40, λ = 0.01, P 0 = 3.33). V. RELATIONSHIP BETWEEN qα BL AND qα BU In Fig. S3(b), we show examples of spatial profiles of the shape index q α . We find that the shape index for the top and bottom layers are outlier due to the flat boundaries which induce large cell deformation. To avoid the outliers, we eliminated the top and bottom layers when we calculated q α BL .
In Fig. S4(a), we show the dependence of q α BU on the target shape index p 0 for the data in Fig. 5 in the main text (N y = 40, λ = 0.01, P 0 = 3.33). We found that larger σ and p 0 provide larger q α BU .
In Fig. S4(b), we show the relationship between q α BL and q α BU for the same data set. We find q α BU ≈ q α BL . As shown in Fig. S3(b), the spatial profiles of the shape index are nearly homogeneous except on the top and bottom layers. Hence, we obtain q α BU ≈ q α BL .