We consider situations where cells are on a substrate; some cells exist individually on the substrate, while other cells are attached to each other to form a cluster. We focus on the dynamics of the peripheries of the cells, which movements are determined by force balance on the cell boundaries. Each cell has a polarity, and depending on the direction of the polarity, each cell changes the strength of contraction force on the cell boundary. We are concerned with two setups of cell polarity. One is that there is a chemoattractant gradient in a definite direction on the substrate, say the x-direction, and all cells have polarity in the x direction (Sections 3.2–3.6). The other is that cells that form a cluster have different directions of polarity. To be specific, the direction of polarity of each cell tilts with respect to the center of the cluster (Section 3.7). This tilted polarity may be achieved by chance or by using chirality of the cells (Taniguchi et al., 2011; Tee et al., 2015; Wan et al., 2016). A characteristic point of our model is that force balance holds at any parts of the cell boundaries at any time.

In our model, cells on a substrate are represented by polygons; specifically, the α -th cell is represented by a polygon that consists of N α segments and N α nodes (Figure 1A). The number N α can differ from cell to cell and changes with time, based on the following rules: if the segment under consideration becomes longer than some critical length l l o n g * , it is split into two segments by creating a new node at the center of the segment, and if the segment under consideration becomes shorter than some critical length l s h o r t * , it is replaced by a single node whose position is the center of the segment before the replacement (Figure 1B). Using appropriate l l o n g * and l s h o r t * values, the surface of the cell becomes smooth enough, and the dynamics of this discretized model are consistent with those observed in continuous models, as shown in Figures 1C, 2B, C. When a cell adheres to an adjacent cell, the segments and nodes on the cell boundary are shared between the two cells; that is, they have the same segments and nodes at their boundary (Figure 1A). Quantities assigned to the cell boundary, such as surface tension and rigidity of the cell boundary, are obtained by considering each quantity on the surfaces of the two cells, as in the existing vertex models (Fletcher et al., 2014).

Cell boundaries experience two types of forces: one is the frictional forces that are expressed by the dissipation function given in Eq. 1, and the other is the mechanical forces that are expressed by a potential function U in Eq. 2. Frictional forces arise from both external and internal factors. When the segments comprising the cells move with respect to the substrate, each segment bears a frictional force from the substrate, whose strength is proportional to the length and velocity of the segment. The frictional coefficient can depend on the direction of movement relative to the direction of the segment; that is, when the segment moves parallel to the direction of the segment, the frictional coefficient is η ∥ , whereas when the segment moves perpendicular to the direction of the segment, the frictional coefficient is η ⊥ . In addition, it is assumed that each segment has an internal friction; when a segment shrinks or expands, a resistance force arises within the segment. The strength of the internal friction is proportional to the strain rate of the segment, with the friction coefficient μ . Elongation or shrinkage of the segment is assumed to be an affine transformation (i.e., the segment is homogeneously elongated or shrunken). Under these assumptions, we can calculate the dissipation function W , which gives the frictional forces on segments in terms of only the positions and velocities of nodes on the cell surfaces (see Supplementary Appendix A for details), as W = ∑ i j η ∥ i j l i j 2 X ˙ i j ⁡ cos ⁡ θ i j + Y ˙ i j ⁡ sin ⁡ θ i j 2 + l ˙ i j 2 12 + η ⊥ i j l i j 2 [ X ˙ i j 2 + Y ˙ i j 2 + l ˙ i j 2 12 + θ ˙ i j l i j 2 12 − X ˙ i j ⁡ cos ⁡ θ i j + Y ˙ i j ⁡ sin ⁡ θ i j 2 + l ˙ i j 2 12 ] + μ i j 2 l ˙ i j 2 l i j } . (1) Here, l i j is the length of the segment that connects the ith and jth nodes; that is, l i j = r i − r j , where r i = x i , y i is the position of the ith node. R i j = X i j , Y i j is the position of the center of the segment i j , given by R i j = r i + r j / 2 , and θ i j is the angle between the x-axis and the vector from the i -th node to the j -th node. The dot over a quantity represents its time derivative, and the symbol i j under the summation symbol means that the sum is taken over all the segments in the system. Note that the sets of R i j , θ i j and l i j are expressed by the set of positions of nodes, r i . Thus, as mentioned above, W is a function of r i and r ˙ i .

The frictional coefficients η ∥ i j , η ⊥ i j , and μ i j in Eq. 1 can vary from segment to segment. However, our main objective is to show some simple applications of this model. We, therefore, kept these coefficients the same for every segment, such that η ∥ i j = η ⊥ i j = η and μ i j = μ throughout this paper. Differentiating W in Eq. 1 with respect to the velocity the i -th node gives the frictional force F f r i c t i o n i on the i -th node as F f r i c t i o n i = − ∂ W / ∂ r ˙ i (Landau and Lifshitz, 1976). Note that the dissipation function W in Eq. 1 considers the length dependence of the frictional force on the segment. Thus, even if a segment i j is divided into two segments, i k and k j , by adding a new node k , the total frictional forces on the segments i k and k j are the same as those on the segment i j before the division, provided movement of segments i k and k j is the same as that of segment i j . In this sense, the frictional force on the cell surface does not depend on the number of nodes on the cell surface and this model satisfies some continuous limits on the frictional forces.

The other mechanical forces on the cell surface, such as contraction forces coming from the actomyosin network beneath the plasma membrane and hydrostatic pressure from the cytosol, are represented by the following effective potential U, as follows: U = K 2 ∑ α A α − A α 0 2 + K p 2 ∑ α L α − L α 0 2 + 1 2 ∑ i j k κ i j k r j − r i l i j − r k − r j l j k 2 2 l i j + l j k + ∑ i j γ i j t l i j , (2) where K , K p , and κ i j k are non-negative constants. This form of U is quite similar to that given in existing vertex models (Farhadifar et al., 2007; Fletcher et al., 2014; Sato et al., 2015b), although it differs in the introduction of the third term, which represents the bending energy of cell membranes, and the form of γ i j t , which expresses the strengh of surface tension on the cell membrane and is related to cell polarity.

The first term in Eq. 2 represents a pressure acting on the segments arising from the area difference between the current area of the cell, A α , and its preferred constant value, A α 0 . This pressure can be interpreted as the hydrostatic pressure in the cytosol on the membrane of the cell. The α under the summation symbol indicates that the sum is taken over all cells in the system. The second term in Eq. 2 expresses the property that the cell tends to keep the perimeter length L α at some preferable length L α 0 . This term represents the tendency to conserve the amount of cell membrane (Fletcher et al., 2014; Sato, 2017). The third term in Eq. 2 expresses the bending energy of the membrane. The strength of the membrane against bending is characterized by the coefficient κ i j k , such that, if κ i j k is large, the part of the membrane under consideration, expressed by segments i j and j k , is difficult to bend, and vice versa. The symbol i j k under the summation symbol indicates that the sum is taken over all pairs of segments that are connected and adjacent to each other. Lastly, the final term in Eq. 2 expresses the surface tension acting on the membrane, which is the result of both contraction forces from the actomyosin networks beneath the membrane and adhesion of the cell to outside objects, such as adjacent cells or substrate. Surface tension strength, expressed by γ i j , is controlled by the cell via altered expression of proteins, such as myosin, actin, cadherin, and integrin. When contraction force is strong or adhesion is weak at segment i j , γ i j becomes large, whereas when adhesion at segment i j is strong or contraction force is weak at segment i j , γ i j becomes small.

Using this value of γ i j , we can express the polarity of the cell. For example, consider the case where a cell has a polarity in the x-direction, and assume that the cell has a weak contraction force at the leading edge ( x > x 0 , where x 0 is the x-component of the position of the cell center) and a strong contraction force at the rear ( x < x 0 ). This situation is expressed by letting the value of γ i j depend on the relative position of the cell membrane with respect to the center of the cell. That is, a small value is assigned to γ i j at the front of the cell, and a large value is assigned to γ k m at the rear. In this paper, we consider several cases with different assigned γ i j values. If we accept the form of U given in Eq. 2, the total mechanical force F U i from U on the i -th node is given by F U i = − ∂ U / ∂ r i .

In this model, we assume that the inertial force of the cell surface is negligible compared with the mechanical forces in question, and thus, all forces acting on the i -th node coming from U and W must be balanced at all times. That is, − ∂ W ∂ r ˙ i − ∂ U ∂ r i γ j k = γ ^ j k = 0 (3) holds for all i’s at any given time. Here, the symbol γ j k = γ ^ j k in the second term indicates that, after the derivative with respect to r i , each γ j k in the second term is replaced by an explicit value γ ^ j k , which is expressed as a function of r i (Fletcher et al., 2014; Sato et al., 2015b; Sato, 2017). This operation means that γ j k obeys some other dynamics that are much faster than those of r i , such that γ j k immediately reaches a value determined by r i (Sato, 2017). If we accept that the dynamics of localization of molecules related to contraction, such as myosin, on the cell surface (tens of seconds) is much faster than that of cell membrane movement (a few minutes), this operation may be allowed. Eq. 3 gives the time evolution equations for r i .

To numerically solve Eq. 3, we first nondimensionalize the variables and parameters appearing in the equation, using the following units: length, A 0 / π ; time, η / K 0 A 0 / π ; and energy, K 0 A 0 / π 2 ; where K 0 is a typical value of K solving for U in Eq. 2. Hereafter, we imply that the variables and parameters appearing in our model are nondimensionalized with the above units. For example, A 0 = π and η = 1 . We then numerically solve the nondimensionalized Eq. 3 for r i by the Euler method, with the step size dt = 1/5000 or 1/10000. For each step, we further determine whether the length of each segment satisfies the replacement criteria l l o n g * < l i j or l s h o r t * > l i j . Any segments that satisfy either of these two conditions are replaced by the rules shown in Figure 1B.

To test our model, we first determine whether the dissipation function W given in Eq. 1 appropriately expresses the frictional forces on the cell surface. For this we consider a case where a circular cell has a constant surface tension γ = 1 and no constraints on its area A and perimeter L ; that is, we set K = 0 , K p = 0 , and κ i j k = 0 in Eq. 2. In this case, the circular cell shrinks with some speed, while maintaining its circular shape. The time series of the radius of the circular cell, r t , is analytically obtained (Supplementary Appendix B), resulting in r t = μ η W 0 η μ e − 2 μ γ 0 t − η 2 r 0 2 + μ ⁡ log ⁡ r 0 , where W 0 x is the Lambert W function (the product logarithm) that satisfies x = W 0 x   e W 0 x for x > − 1 / e . Comparing this analytical solution with the results of numerical simulation for Eq. 3 reveals good agreement (Figure 1C). Thus, we conclude that the friction force expressed by Eq. 1 appropriately expresses the frictional force experienced by the cell membrane, at least for this simple test case.

Next, we investigate a more realistic case wherein the cell under consideration is circular with a constant area A = A 0 and has a polarity in the x-direction. Constant area of the cell is achieved by using the parameters K ≫ 1 , K p = 0 , and κ i j k = 0 in Eq. 2, and cell polarity is expressed by the direction-dependent surface tension γ ^ i j in Eq. 3 as γ ^ i j = γ 0 − a ⁡ cos ⁡ θ i j α , (4) where γ 0 is a positive constant, a is a non-negative constant, and θ i j α is the angle between the x-axis and the vector connecting the center of cell α that contains segment i j and the center of segment i j (see Figure 1A). The value of a represents the degree of the polarity, such that, when a = 0 , the cell is isotropic and has no polarity, whereas when a > 0 , the cell boundary at a relatively large x has a weak surface tension, and the cell boundary at a relatively small x has a strong surface tension. Hereafter, we refer to the region of cell boundary at a relatively larger x as the “front” of the cell, and the region of the cell boundary at a relatively smaller x as the “rear” of the cell. If a = 0 , the shape of the cell is circular due to the constant area and constant surface tension, γ 0 , on the cell surface. In this simulation, we set a as a small positive value ( a ≪ γ 0 ), such that the circular shape of the cell is still retained.

We then examined the steady state of the cell with these model parameters. Numerical simulations reveal that the cell moves in the positive x-direction, with a constant speed, while keeping a circular shape in the steady state (Figure 2A). The driving force for this movement is the direction-dependent surface tension in Eq. 4 and the resulting cell surface flow from the front to the rear (Supplementary Movie S1). That is, in the front region of the cell, surface tension is relatively weak due to the parameters of Eq. 4, such that the surface in the front region tends to expand. In contrast, surface tension in the rear region of the cell is relatively strong, and hence, the surface in this region tends to shrink. This surface tension-dependent tendency of the cell surface to expand or shrink produces a flow of cell membrane, in which the front region is always expanding (in some sense, blebbing is continuously occurring in front of the cell), whereas the rear region is always shrinking (Supplementary Movie S1). These behaviors act as a source and sink of cell surface and yield a flow of cell surface from front to rear. In addition, within this system, there is a frictional force between the cell surface and the substrate. Thus, if the cell surface moves in some direction, the whole cell experiences forces that move it in the opposite direction, based on the action–reaction principle.

The velocity V of this cell movement resulting from the direction-dependent surface tension can be analytically calculated, by treating the cell as a continuous circular object (see Supplementary Appendix C for details). The result ist V = a / η R   1 + 2 μ / η R 2 , 0 , (5) where R is the radius of the circular cell. We then compared this analytical solution with the cell speed obtained by numerically solving Eq. 3 and found that these values are in good agreement with each other (Figures 2B, C). This indicates that our discrete model described by Eqs. 1–3 can describe continuous cell surface dynamics if we choose appropriate values for parameters, such as l l o n g * and l s h o r t * .

We further find that single-cell migration induced by direction-dependent surface tension, described in Eq. 4, occurs even when the cell shape is elliptical (Figure 2D). We can model an elliptical-shaped cell using finite values for K p and κ i j k in Eq. 2 and by ensuring the preferred cell perimeter L 0 is longer than that of a circle with the area A 0 , 2 π A 0 . Overall, the mechanism for cell movement is basically the same as for circular cells, that is, the surface flows from front to rear due to direction-dependent surface tension, described in Eq. 4. In this simulation, we set the initial configuration of the cell, such that the long axis is aligned in the y-direction (Figure 2D; t = 0). Initially, up until t = 8, the elliptical cell moves in the x-direction, while keeping the short axis of the cell aligned with the x-axis. However, at around t = 25, the long body axis begins to incline toward the x-direction, and the direction of cell movement also begins to incline; that is, the cell has the y-component of velocity (Supplementary Movie S2). In the final stage, the long body axis is completely oriented to the x-axis, and the cell moves in the x-axis again. The final speed of the cell is faster than in the earlier stage, where the long axis of the cell is perpendicular to the x-axis. This simulation indicates that if the cell shape is elliptical and the direction of cell polarity is fixed in the x-axis, cell movement in which the short axis is oriented to the x-axis is slow and unstable, whereas cell movement in which the long body axis is oriented to the x-axis is fast and stable. In general, the parameters of Eq. 4 make cell movement fast if the long axis of the cell is directed in the direction of cell polarity. This property clearly appears in the next case, where the cell is sandwiched by two parallel walls.

Motivated by results from published experiments (Bergert et al., 2015; Liu et al., 2015; Sakamoto et al., 2022), we next examined the scenario in which a cell is sandwiched by two parallel walls. In this case, the wall is expressed using the potential U w a l l y = y − d Θ y − d − y + d Θ − y − d , where d is one-half the distance between the walls, and Θ x is the Heaviside step function, defined as Θ x = 1 , for x > 0 , and Θ x = 0 , for x ≤ 0 . By adding the terms ∑ a l l   n o d e s   i ′ s U w a l l y i to Eq. 2, the cell surface close to the walls experiences a repulsive force from the wall potential. In this setup, the area of the cell is held constant as A = A 0 , which is achieved with the same parameters as used for the simulation with a circular cell. Under these conditions, the cell also migrates in the x-direction in the steady state (Figure 3A), where the speed of the cell increases as d decreases (Figure 3B), consistent with published experimental results (Liu et al., 2015; Sakamoto et al., 2022). To demonstrate how and why cell speed increases with decreased d , we derive the analytical expression for the speed of the cell (see Supplementary Appendix C), which is given by V x = 2 π η d ∫ 0 π / 2 ⁡ sin ⁡ ξ 1 d γ θ f o r w a r d ξ 1 d ξ 1   d ξ 1 + ∫ π / 2 π ⁡ sin ⁡ ξ 2 d γ θ b a c k w o r d ξ 2 d ξ 2   d ξ 2 + 2 Δ γ π η d , (6) where θ f o r w a r d ξ 1 = arctan d ⁡ sin ⁡ ξ 1 d ⁡ cos ⁡ ξ 1 + L x , θ b a c k w o r d ξ 2 = arctan d ⁡ sin ⁡ ξ 2 d ⁡ cos ⁡ ξ 2 − L x , and Δ γ = γ π − θ 1 − γ θ 1 = 2 a L x / L x 2 + d 2 . The angles θ f o r w a r d , θ b a c k w a r d , ξ 1 , and ξ 2 used here are defined as shown in Figure 3C. The values L x and d are the half-length of the straight part of the cell and the radius of the circular part of the cell, respectively. Because the total area of the cell is kept at A = A 0 , L x and d are related as 4 L x d + π d 2 = A 0 .

Results from analytical expression Eq. 6 are in good agreement with the results of numerical simulations (Figure 3B). Thus, we use Eq. 6 to interpret the cell speed in this scenario. The first term in Eq. 6 is the contribution from surface flow on the two semicircles of the cell, and the second term is the contribution from surface flow on the straight parts of the cell. When d becomes small, the semicircles also become small, and the flow speed of cell surface along the semicircles becomes slow due to the surface tension gradient along the small semicircles. Thus, the first term in Eq. 6 becomes small when d decreases (Figure 3B, yellow dashed curve). In contrast, Δ γ in the second term in Eq. 6, which is the surface tension difference between the two edges of the straight parts of the cell (Figure 3C), increases with decreased d (see Eq. 4). Furthermore, d in the denominator of the second term in Eq. 6, which represents the resistance of the semicircle parts to the cell movement, also enlarges this term as the value of d decreases. Thus, the whole second term in Eq. 6 drastically increases when d decreases (Figure 3B, green dashed curve). The tendency that the migration speed v of the sandwiched cell increases with decrease in d appears even when we take a different setup of γ ^ i j . For example, we consider the case where the surface tension is constantly increased along the cell surface form front to rear, in which the minimum and maximum surface tensions are kept to be constant, i.e., γ ^ i j = 0.8 at the front and γ ^ i j = 1.2 at the rear. These values are the same as those in the simulations in Figure 3B. The results of numerical simulations with this setup is given in Supplementary Figure S2, where as in Figure 3B v increases with decrease in d , while the d dependence of v is somewhat moderate compared to Figure 3B. This same tendency of v implies that the increase in v of the sandwiched cell with decrease in d may be a robust phenomena, as frequently observed in experiments (Bergert et al., 2015; Sakamoto et al., 2022).

To better understand why and how cells move due to surface flow while keeping the forces balanced at all times, we surveyed the properties of cell movement observed in this model. Although these properties derived below are provided in a continuous form, it is not difficult to translate them into a discrete form. As noted above, our model neglects all inertia of the cell surface, so that the forces on any element of the cell surface are balanced at any time, and more importantly, the total force on the cell surface from the substrate must vanish at any time. Because the force on the cell surface from the external object—the substrate—is friction only, we have the following equality: ∫ Ω η v ξ , t ∂ s ξ , t ∂ ξ d ξ = 0 , (7) at any time t . Here, ξ is the material coordinate that is assigned to the element of cell surface under consideration and s ξ , t is the counter length of the arc of cell surface between some reference point on the cell surface and the cell surface element specified by ξ at time t . Additionally, v ξ , t is the velocity of the surface element at ξ and t ; Ω under the integral symbol indicates that the range of the integral is the whole cell surface. We then consider the situation where cell movement reaches a steady state, with cell velocity and cell shape constant in time, focusing on some time t 0 in this steady state. At t = t 0 , we reassign the material coordinate ξ on the cell surface, such that ξ coincides with the counter length s ; that is, s ξ , t 0 = ξ . (8)

With this situation, let us consider the time evolution of the cell surface density, ρ . In general, ρ obeys the following time evolution equation ∂ ρ ξ , t / ∂ t = − ∂ r / ∂ ξ ⋅ ∂ v / ∂ ξ / ∂ s / ∂ ξ 2 ρ + J ξ , t , (9) Where r ξ , t is the position vector of the surface element at ξ and t , and J is the flux of cell surface component from the cell inside to the cell surface per unit length. The derivation of Eq. 9 is shown in Supplementary Appendix D. In the steady state, ρ is constant in time, say ρ = ρ 0 , and the time derivative ∂ ρ / ∂ t is zero. Thus, from Eqs. 8, 9 we obtain J ξ , t 0 = ρ 0   ∂ r / ∂ ξ ⋅ ∂ v / ∂ ξ , (9.1) which describes the relationship between J , r , and v in the steady state. Because of the constant cell shape in the steady state, the total flux of the cell surface components must vanish at any time. Thus, ∫ Ω J ξ , t 0 d ξ = 0 . (10)

Interestingly, J in Eq. 9.1 is related to the cell migration velocity V in the steady state for the following reasons. In the steady state, the cell surface center of mass also moves with the same velocity V . However, as indicated in Eqs. 7, 8, there is no net velocity of components on the cell surface, i.e., ∫ Ω v ξ , t 0 d ξ = 0 . This implies that the shift in cell surface center of mass is, in fact, achieved by the inflow and outflow of cell surface components from inside the cell. Thus, the integral of the product of J ξ , t 0 and r ξ , t 0 over the whole cell surface gives the transport rate of cell surface components through the inside of the cell. Dividing this quantity by the total amount of components on the cell surface, ρ 0 L c e l l , yields the velocity of the cell surface center of mass. Therefore, V = ∫ Ω r ξ , t 0 J ξ , t 0 d ξ / ρ 0 L c e l l , (11) where L c e l l is the cell perimeter in the steady state. Eq. 11 is interesting because this relation connects apparently different quantities, the cell migration velocity V and the turnover rate J of the cell surface. Indeed, we can confirm that Eqs. 7, 10, and 11 exactly hold for the circular cell migration case with μ = 0 as shown below. From the analytical results given in Equation (C.16) in Supplementary Appendix C, r and v of the circular cell in the steady state are given as r ξ = a η R t + R ⁡ cos ξ R , R ⁡ sin ξ R v ξ = a η R − 2 a η R sin ξ R 2 , 2 a η R sin ξ R cos ξ R . (12)

From Eqs 9.1–12, we have J = ρ 0   2 a η R 2 cos ξ R . (12.1)

Inserting Eqs. 12–12.1 into Eq 11 gives V = a η R , 0 , (12.2) where we have used L c e l l = 2 π R and the fact that the range of ξ is 0,2 π R . Eq. 12.2 coincides with Eq. 5 for the case of μ = 0 .

The above simulations focus on migration of individual cells; however, cells often move together with other cells. Therefore, we next examined whether cell migration due to direction-dependent surface tension occurs even when multiple cells are attached to each other, forming a cluster. Specifically, we investigated the case where the number N c e l l of cells on the substrate is N c e l l = 2,3,4,10 with some initial configuration in which the cells are attached to each other (Figure 4C). In these simulations, we have set the surface tension on the boundary between the cells to a constant value (i.e., γ = γ b ), by assuming that the same type of cells is adhering to each other with some characteristic strength (Figure 4A). The other cell boundaries, which do not contact other cells, have the direction-dependent surface tension specified by Eq. 4. As in the case of single-cell migration investigated in Sections 3.2–3.4, each cell in this system keeps its area constant and has no constraint on its perimeter; that is, K = 100 , K p = 0 , and κ i j k = 0 in Eq. 2. Under these conditions, multiple cells move in the x-direction, while maintaining attachment between cells in the steady state (Figure 4C, Supplementary Movies S3–S10). Moreover, the mechanism of multiple-cell migration is basically the same as that of single-cell migration. That is, due to the direction-dependent surface tension of each cell, the surface of the cell cluster flows from front to rear, and this flow drives movement of the whole cluster.

The shape and velocity of the cluster during this movement depend on the value of γ b . When γ b decreases, the shape becomes round, and speed becomes slow (Figures 4B, C). The roundness of the cell cluster at small γ b results from the constant area and lack of constraint on the perimeter of each cell. Due to these parameters, when γ b becomes small, boundaries between cells tend to extend, and the whole cell cluster behaves like one object. Further, because the cell cluster surface has a constant surface tension, γ 0 , from Eq. 4, which reduces the cluster perimeter as much as possible under the constant area, A c l u s t e r = N c e l l A 0 , the whole cluster is nearly a circle. In contrast, when γ b becomes large, boundaries between cells tend to shrink, and outer cells within the cluster tend to be round, with the surface tension γ 0 (Figure 4C).

The slow cluster movement speed at small γ b results from the roundness of the cell cluster. As shown in Sections 3.2–3.4, surface flow-induced cell migration, in general, becomes fast when the cell is elongated in the x-axis. This tendency is most evident in the case of a cell sandwiched by two walls, described in Section 3.4, where Δ γ in Eq. 6 acts as the driving force for movement. Similarly, in the case of cluster migration, the surface tension gradient along the cell surface (and the resultant flow of cell surface) is also the driving force for movement. Thus, if cells within the cluster are relatively elongated in the x-direction, the speed of cluster movement is relatively fast. However, when the cluster is round, such as for N c e l l = 2 , 4 and γ b = 0.1 (Figure 4C), each cell in the cluster is relatively elongated in the y-axis, and thus, the driving forces are small, and the movement becomes slow.

To this point, we have examined only scenarios in which all cells on a substrate have the same direction of polarity (e.g., the x-direction). However, because clustered cells may also show distinct polarities, we next investigated the case wherein direction of polarity for individual cells within a cluster differs from cell to cell. In particular, we considered four cells comprising a cluster, in which each individual cell has a polarity that is directed perpendicular to the direction from the center of the cell to that of the cell cluster (Figure 5A). Surface tension for cells within this situation is expressed as γ ^ i j = γ 0 − a ⁡ cos θ i j α , c − 3 π / 2 , (13) where θ i j α , c is the angle between the vector from the center of cell α that contains segment i j to the center of the cluster and the vector from the center of cell α to the center of segment i j (Figure 5A). The other parameters for this simulation are the same as those for cluster migration investigated in Section 3.6; that is, boundaries between cells within the cluster have a constant surface tension, γ i j = γ b , and each cell in the cluster keeps its area constant and has no constraint on its perimeter, with K = 100 , K p = 0 , and κ i j k = 0 in Eq. 2.

The initial configuration of the cell cluster is the same as in Figure 4C ( N c e l l = 4 ). Numerical simulations with the above parameters further indicate that each cell within the cluster continues to rotate counterclockwise around the center of the cluster while maintaining its attachments (Figure 5B; Supplementary Movies S11, S12). This rotation continues as long as each cell has the polarity established in Eq. 13, and the mechanism of movement is basically the same as observed for cluster migration and single-cell migration. That is, cells move due to the direction-dependent surface tension expressed by Eq. 13; the surface of the cell cluster moves clockwise, and this surface movement generates the driving force for rotational movement of the whole cluster. In this simulation, we set the polarities of the cells in the cluster with Eq. 13. However, in reality it may be more reasonable that the direction of cell polarity changes with time by following some other rules. One such rule is called the velocity alignment mechanism (Camley et al., 2014), in which cells align their polarity to their velocity. We incorporated this rule into our model and performed numerical simulations. The results show that for any initial directions of cell polarity, the cell polarities finally align such that the cluster of cells rotates (Supplementary Movies S13, S14). The direction of rotation depends on the initial distribution of cell polarities in the cluster. That is, the cell cluster rotates clockwise or counter clockwise depending on the initial distribution of cell polarities. Our system does not require cell confinement for this rotational motion, which is different from the results of previous studies (Camley et al., 2014).

As in the case of single-polarity cluster migration examined in Section 3.6, shape and rotational velocity of the cell cluster change with γ b (Figures 5B, C). When γ b is small, the cell cluster becomes round, and the angular velocity for the rotational movement of the cluster becomes large. The reasons for this cluster roundness are the same as those outlined in Section 3.6, that is, the whole cluster behaves like one object due to a lack of constraint on cell perimeter and small γ b . In addition, the constant surface tension, γ 0 , reduces cluster surface as much as possible under the constant whole area, A c l u s t e r = N c e l l A 0 .

The faster rotation of the cell cluster at small γ b originates from the velocity distribution of the cell cluster surface (Figure 5Da), and the reasons for this are as follows. First, cell boundaries within the cluster are classified into one of two types: 1) Cell boundaries that are located at the surface of the cluster and are associated with only one cell, referred to as the “outer cell boundaries” and 2) cell boundaries that are located within the cluster and form boundaries between two adjacent cells, referred to as “inner cell boundaries.” During rotation of the cluster, the driving forces for cluster movement come from the flow of the outer cell boundaries, because only the outer cell boundaries have cell polarity, as indicated by Eq. 13. In the case of the round cluster ( γ b = 0.2 ; Figure 5Da), most velocities of the outer cell boundaries are directed in the tangent of the cluster surface, which comes from the monotonic decrease in surface tension along the surface. These velocities generate a large magnitude of torque compared with those present in the case of γ b = 1.0 (see Figures 5D, E). Here, N z f l o w is the z-component of the torque generated by the surface flow of the outer cell boundaries and is evaluated as N z f l o w = − η ∑ i j o u t e r R i j − r c × R ˙ i j − r ˙ c z l i j , where r c is the position vector of the cluster center, and i j o u t e r indicates that the summation range is over all outer cell boundaries. N z f l o w is balanced with the z-component of torque generated from frictional forces experienced by the inner cell boundaries, denoted by N z f r i c t i o n , and this balance determines the angular velocity, ω , of the cluster rotation. N z f r i c t i o n is roughly expressed as N z f r i c t i o n = − ζ   ω , where ζ is a coefficient given by ζ = η ∑ i j i n n e r R i j − r c 2 l i j . The ζ coefficient is almost constant in time because it is determined by the shape of the network of inner cell boundaries, which is also almost constant in time (see Supplementary Movies S11, S12). The torque balance on the cell surface gives ω = N z f l o w / ζ . Although both N z f l o w and ζ increase with decreased γ b (Supplementary Figure S1), the increase in ζ is more gentle than for N z f l o w , which results from the fact that the size of the network of inner cell boundaries does not change much with γ b (see Figure 5B). Thus, the increase in N z f l o w , which occurs with decreased γ b , increases the angular velocity ω .