Collision-induced torque mediates transition of chiral dynamic patterns formed by active particles

It is still challenging to control dynamic self-organization patterns of self-propelled particles. Although varieties of patterns associated with chirality have been observed, essential control factors determining patterns remain unclear. Here, we explore numerically how torque upon particle collision affects dynamic self-organization. Based on the particle-based model with both collision-induced torque and torque in self-propulsion, we find that introducing collision-induced torque turns homogeneous bi-polar orientation templated by bi-directional alignment into rotating mono-polar flocks.


I. INTRODUCTION
Controlling the patterns that emerge through dynamic self-organization of motile objects, as observed in the collective motion of living organisms [1][2][3], has been a long-standing great challenge. Recently, the dynamic patterns observed in nature have been contemplated in artificial systems by using the cytoskeletal filament microtubules (MTs) and F-actins. Upon propulsion by their associated biomolecular motors, the cytoskeletal filaments exhibited collective motion mediated by attractive interaction [4][5][6][7][8][9] or crowded conditions [10][11][12].
Local alignment interactions, i.e., collision-induced alignment of two self-propelled filaments, play important roles in the emergence of their collective motion [10][11][12]. The resulting patterns of circular mesoscopic structures, streams, and vortices exhibited local or global rotational motion along the clockwise (CW) or counterclockwise (CCW) direction [6-9, 12, 13] that resembled the coherent motion exhibited by the living organisms in nature.
In nature, a rotational force or torque of cytoskeletal filaments is known to play an important role in chiral morphogenesis of cells, tissues, and organisms [14][15][16]. Similarly, in the artificial systems made of cytoskeletal proteins, chiral collective behaviors have often been observed [6,7,9,12,13,17]. Such observations and demonstrations motivated to explore the collective dynamics of chiral self-propelled objects theoretically by employing analytical and numerical approaches. Among these efforts, the most popular strategy was the introduction of chirality, i.e., the left-right asymmetry, to the motion of the objects translocating in two dimensions [18][19][20][21][22][23][24]. Chirality was also introduced in relation to the mutual interactions of the objects [25]. A more straightforward strategy has been to assume the objects with the finite-size and shape explicitly, instead of assuming the point-objects, and put chirality or LR asymmetry in the object shape [26][27][28]. Indeed, this approach has been useful in studying artificial systems in which the objects are of well-defined shape in a mechanistic sense; e.g., gliding assay of cytoskeletal filaments [28]. Such self-propelled objects with a chiral shape can exhibit varieties of dynamic self-organization.
In view of these two approaches to introduce chirality, the shape chirality can influence both the spontaneous motility and interaction of the objects. However, which chirality plays the crucial role in facilitating the emergence of dynamic self-organization patterns has not been figured out yet. For example, recent study demonstrates well-defined chiral mono-polar flocking of MTs, or dense MT cluster in which motility directions of MTs are aligned uni-directionally, in a gliding assay on kinesins [29], but it remains unclear how the chirality of MTs contributed to such mono-polar flocking. Indeed, the difference of chirality in spontaneous motility and interaction of the objects can be a hint of this as follows. Experimental observations in Ref. [29] suggested mono-polar flocking is attributed to chirality in interaction. In contrast, Refs. [13,17] report homogeneous bi-directional orientation of MTs with chiral rotational motion, in a gliding assay on kinesins, and such rotating bidirectional orientation can be explained by chirality in spontaneous motion [13]. Like this, to obtain a comprehensive understanding of the factors that determine such differences in dynamic self-organization patterns with chirality, it appears inevitable to dissect the effect of chirality on the motility and interactions of self-propelled objects.
Here, we have demonstrated a systematic in-silico study on collective motion of selfpropelled particles (SPPs) each of which has an intrinsic chirality of the both types as mentioned above; namely, chirality in self-propulsion and interaction. We consider SPPs with intrinsic polarity where the SPPs move on a two-dimensional substrate. The SPPs interact with each other through isotropic core repulsion and bi-directional alignment. We chose such bi-directional alignment to avoid the emergence of mono-polar phase purely by the alignment interaction, which is the case for MTs in a gliding assay. Two types of torque, self-propelled torque (ST) and collision-induced torque (CT) (See Fig. 1 and below for more details), are applied to the particles as left-right (LR) asymmetric motility due to twodimensionality. We have investigated the emergence of patterns by tuning the strength of ST and CT without manipulating their alignment interactions. We found that when the CT is introduced, transition from bidirectional orientation to mono-polar flocking takes place although the alignment interaction is bi-directional. The emergence of mono-polar flocking mediated by chirality was reported in Ref. [28], where the authors investigated a mixture of two types of filaments having opposite chirality. On the contrary, the results presented in this article predict another mechanism to account for the chirality-induced mono-polar flocking, according to which such a combination of two types of filaments is unnecessary.
Before moving on, we recall the robust features experimentally observed in the collective motion of chiral MTs driven by kinesins, presented in Ref. [29]. Chirality in the MT structure was introduced by polymerizing tubulins with a certain nucleotide, GMPCPP.
The GMPCPP-MTs are found to align upon collision and eventually form mono-polar flocks.
Notably, these flocks rotate dominantly in the CCW direction [29]. In addition, when the imply a correlation between the collision-induced torque and mono-polar flocking, which has motivated us to build up the simulation setting in this study.

II. MODEL AND METHOD
We consider N SPPs, which are located at x j = (x j , y j ) (j = 1, 2, · · · , N) and have intrinsic polarity q j = (cos θ j , sin θ j ) [30], in a regular rectangle space with periodic boundaries in two dimensions [13,[31][32][33]. We assume that locations and polarity directions of SPPs (j = 1, 2, · · · , N) obey and respectively, with q ⊥ j = (− sin θ j , cos θ j ). Equation (1) assumes the over-damped dynamics, and each object moves along the polarity q j with a constant velocity v 0 [30] in the absence of volume exclusion interactions. The volume exclusion, by which two SPPs mechanically interact with each other, is implemented by The summation j ′ (n.j) runs for all the neighbors j ′ of j-th SPP, defined by |∆x j,j ′ | < r with ∆x j,j ′ = x j − x j ′ . Equation (1) also assumes that each SPP hardly moves along the direction perpendicular to the polarity direction, using the rescaled anisotropic friction with the ratio R ζ = ζ /ζ ⊥ of friction coefficients in parallel ζ and perpendicular directions ζ ⊥ [13]. (⊗ is the tensor product, and I is the identity matrix.) This anisotropy in friction has been implemented to phenomenologically reflect the observation that, when a gliding MT collides with another MT from its side, the colliding MT either stops to move or gets over the other, and the collided MT does not move into the direction perpendicular to its polarity [13]. This anisotropy is not essential for the phenomenology focused in this paper, as shown below in Appendix ("The case with isotropic mobility"). Equation (2) assumes that polarity is spontaneously established for each SPP with a fixed amplitude, |q i | = 1, and only its direction can evolve over time [30].
The first term indicates bidirectional alignment interaction [13], See Appendix ("Details of the theoretical model") for details. The coefficient α AL indicates the strength of alignment. The second term ξ j (t) represents Gaussian white noise with In this study, we also apply two different types of torque, self-propelled torque (ST) and collision-induced torque (CT), to the particles, which are represented by the last two terms in Eq. (2) [ Fig. 1]. Note that, since SPPs are gliding on the substrate, the special directionality exists in the z-axis and here we are focusing on only the other two dimensions; reflecting this fact, we assumed torque as representation of LR asymmetric motility. In ST, we assume that intrinsic polarity q j of each SPP rotates with a given speed ω ST in either CCW and CW direction. ω ST denotes the strength of ST [ Fig. 1 left], which we assume is a given constant. It is to be noted that ST has been observed in the gliding assay of MTs [6]. We further assume that, when SPPs collide, another torque is exerted on their intrinsic polarities, which we name CT. Ω  [29], an increase in the mean curvature of the trajectory was observed upon increasing the MT-density, which suggests that our assumption for CT is not artificial. In the absence of these LR asymmetries, this model is essentially the same as that given in Ref. [13].
We numerically calculate Eqs.

III. RESULTS
First, we examine the possible patterns which may emerge based on our theoretical model system. Typical snapshots are shown in Fig. 2. For a two-dimensional hard-disk system, the formation of an ordered phase requires a packing fraction higher than 0.7 (ρ = 0.89) [34]. Despite a low global particle density, ρ of 0.2, mono-polar flocks are observed in silico, as shown in Fig. 2 When the ST and CT were opposite to each other, we found an island region in the lower right of Fig. 3(b) and (d) for the flocks rotating in the CW direction. The doughnut-like shapes in the island are the same as those in the mono-polar area on the upper left region in Fig. 3(d). However, the SPP number fraction in the polar order region of the island area is larger than that of the upper left region. Thus, the collision probability in the flock increases, and the density increases by the collisions as ω ST becomes large. Therefore, the increase of the collision frequency caused by the ST leads to another mono-polar flocking phase in combination with the CT of the opposite direction.In contrast, when the CT is the same direction as the ST, as shown in Fig. 3(c), this island area is not observed. This is consistent with the above statement because the density seems to always decrease for increasing ω ST in this case.
As mentioned above, the rotating flocks are observed for high values of ω CT [ Fig. 3(a,b)].
As are canceled by each other. Furthermore, the slope of the line was roughly 0.3. We also evaluated the mean contact number around each particle in mono-polar flocks m i i∈MPF in the numerical calculation, for the parameter window exhibiting the polar order, and it was around 2 to 5. Hence, the above slope agrees with the condition for vanishing torque, It is to note that, since the global particle density, ρ has been set at 0.2, this means the contact number and the slope are the consequence of flocking.
Finally, we discuss the mechanism behind the violation of bidirectional orientation and the formation of mono-polar flocks mediated by only the CT based on the three-SPP simulation (See Appendix "Three-particle simulation" for the details). The increase of CT changes a stable state from bidirectional orientation to the rotating mono-polar flock [ Fig. 2(a)].
Flocking is mediated by the alignment interaction and CT. Although both the mono-polar and bidirectional orientation of motion can be stabilized by alignment interaction, simulations for a few SPPs in Fig. 4(a) and Appendix ("Three-particle simulation") reveal that, in the absence of CT, bi-polar orientation is stable and mono-polar flocks are rarely formed.
On the contrary, CT rotates the direction of movement of the SPPs moving in the same direction as a cluster, which breaks the bi-polar orientation and provides more chances of mono-polar flocking [ Fig. 2(b)]. For particles moving in the same direction, alignment interaction worked as an effective attractive interaction and maintained the mono-polar flock once it is formed. Therefore, when alignment interaction is strong enough and the density of SPPs is large enough, this effective alignment among the SPPs moving in the same direction may result in phase separation, and allows the emergence of the mono-polar flocks with high local density and large local polar order. In fact, such dependency on the associated parameters is seen in Appendix ("A few notes regarding density in flocks"). Note that we could not find the significant difference between the cases with and without ST [compare red circles and blue triangles in Fig. 4(b)], which is a stark difference from the results shown in Fig. 3(c,d), where we simulated the many-SPP case (N = 20, 000) and observed the ω ST -dependence of the SPP number fraction within mono-polar flocks. This difference implies that the many-body effect is indispensable to recapitulate the ω ST -dependence of the threshold ω CT seen in Fig. 3.

IV. CONCLUSIONS
In conclusion, through an in silico study, we have clarified how the types of torque of SPPs, i.e. torque due to the collision and torque associated with their self-propulsion, can affect their coherent dynamics. By varying the magnitudes of these two types of torque, we have discovered that there is a transition between different forms of coherent dynamics that are manifested by homogeneous bidirectional orientation and mono-polar flocking, which have been observed in microtubule-kinesin gliding assays in Refs. [13,17] and Ref. [29], respectively. When the self-propelled torque is dominant, SPPs maintain their homogeneous bidirectional orientation, although their direction rotates. We discovered that an increase in collision-induced torque breaks the homogeneous bidirectional order and the stabilized mono-polar flocks. Our results clarify the role of collision-induced torque in the emergence of their coherent dynamics, and resolve the discrepancy in the observations mentioned above.
The findings in this work point out the importance of the type of torque as a control factor in the dynamic self-organization patterns of SPPs.
V. APPENDIX

Details of the theoretical model
We consider N particles in a square box with periodic boundaries in two dimensions, and assume that each particle is a self-propelled particle (SPP) with an intrinsic polarity along which the domain tries to move [13]. The location and intrinsic polarity of the j-th particle are described by x j = (x j , y j ) and q j , respectively (j = 1, 2, · · · , N). We assume that the velocity v j , which determines the x j 's time evolution as and the polarity q j of the j-th particle obeys and under the constraint |q j | = 1, respectively, for every j. Equation (6) assumes the overdamped dynamics, and each particle moves with a constant velocity v 0 in the absence of volume exclusion interactions. Equation (6) also assumes that each particle hardly moves along the direction perpendicular to the direction of q j , which has been implemented by the (rescaled) anisotropic friction tensor Θ(q j ) = q j ⊗ q j + R −1 ζ (I − q j ⊗ q j ) with the ratio R ζ = ζ /ζ ⊥ of friction coefficients in parallel ζ and perpendicular directions ζ ⊥ . Here, ⊗ is the tensor product, and I is the identity matrix. (The reason we introduced such anisotropy in friction is as follows: Ref. [13] observed that, when the gliding microtubule collides another microtubule from its side, the colliding microtubule either stops moving or gets over the other. In other words, phenomenologically, the microtubule collided from its side does not move into the direction perpendicular to its polarity. Since this work is motivated by the observations in microtubule gliding assay, as we detailed in Introduction, we reflected this fact phenomenologically by using the anisotropic mobility, or inverse friction, and below setting the friction perpendicular to the particle's intrinsic polarity to be much higher than the parallel counterpart, with which the particle hardly moves to the perpendicular direction indeed. However, it is to be noted that this anisotropic mobility is not essential in our main result.) The term ξ q j (t) represents the noise, for which, for simplicity, we assume a white Gaussian noise with ξ q j = 0 and where the subscripts k and l specify the directions, k, l = x, y, with the statistical average · . The coefficient D indicates the noise strength. The particle-particle interactions are given by J VE j and J AL j , which represent the volume exclusion and bidirectional alignment interaction with the interaction ranges r, respectively. The volume exclusion is given by The summation j ′ (n.j) runs for all the neighbors j ′ of j-th particles, defined by |∆x j,j ′ | < r with ∆x j,j ′ = x j − x j ′ . Here, we have assumed not rigorous but soft volume exclusion, and the softness is controlled by the factor β. Bidirectional alignment interaction is given by [13] The coefficient α AL indicates the strength of alignment. We can scrutinize the meaning of this term by rewriting it in the potential form; with The summation j,j ′ ∈n.p. runs for all the neighboring j-th and j ′ -th particle pairs, defined again by |∆x j,j ′ | < r. As Eqs. (11) and (12) indicate, this interaction term align the polarities of neighboring pair of particles in the bidirectional way, i.e. toward either θ j −θ j ′ = 0 or θ j − θ j ′ = π. In this study, we naively assume given constants α AL and β for the interaction, but there are other possible choices. For example, in some literature on collective motion of SPPs, the interaction is defined in the way that its strength depends on the SPP speed [35]. This may be the point which one has to be careful when investigating the phase diagram over e.g. v 0 . In this paper, we focus on only the torque strengths, mentioned below, so that this choice may not affect the conclusion.
The last two terms of Eq. (7) are the key terms of this study, which give rise to chirality.
As mentioned in Introduction, in this article, we assume that chirality affects the system's dynamics through the two different ways: one is the self-propelled torque (ST) ω ST ≡ ω ST q ⊥ j , and the other is the collision-induced torque (CT) Ω CT j ≡ Ω CT j q ⊥ j . Here, q ⊥ j is one of the unit vectors perpendicular to the polarity, given by q ⊥ j ≡ (− sin θ j , cos θ j ). The constant ω ST is the strength of ST. On the other hand, Ω CT j is given by Ω CT j = ω CT m j with the number m j of particles within the range r from the focused particle (j), and the constant ω CT controls the strength of CT. (See Fig. 1 for the schematics.) It is worth noted that, in two dimensions, there are two unit vectors perpendicular to a certain reference vector, corresponding to left or right. The above definition of q ⊥ j selected out one of them. Thus, the existence of q ⊥ j in the definitions of ω ST and Ω CT is expressing the origin of chirality. Equation (7) can be rewritten as the time evolution of the angles θ j of the polarity directions, using q j = (cos θ j , sin θ j ). For this purpose, we may take the inner products of q ⊥ j and the both sides of Eq. (7). As a result, Eq. 2 is obtained, with ξ j ≡ ξ j · q ⊥ j . The new noise term ξ j (t) represents the angular noise, and since ξ j = −ξ q x,j sin θ j + ξ q y,j cos θ j , it is again a Gaussian white noise with ξ j = 0 and Note that D −1 now characterizes the persistence time of polarity direction, or characteristic time of angular fluctuation.
Equation (6) and Eq. 2, i.e., the angular description of Eq. (7), can be rewritten by the dimensionless forms as and respectively, wherex j = x j /X,t = t/T ,∆x j,j ′ = ∆x j,j ′ /X,β = βT /X,α AL = α AL T , with characteristic length X ≡ r and time T ≡ r/v 0 . The noise term is also rescaled into the new notationξ j (t), which is a Gaussian white noise satisfying ξ j = 0 and withD = DT . In the main text, we applied the same nondimensionalization by putting r = 1 and v 0 = 1. The propulsion strength is often quantified by using Péclet number [18,36]. The Péclet number, or specifically the rotational Péclet number [18], is defined by Pe ≡ v 0 τ p /r, with the migration persistence time τ p [36]. Pe is given by the inverse of the dimensionless noise, Pe =D −1 , because τ p = D −1 in our case. The parameter values which we used in our simulations here correspond to Pe = 100.

Hysteresis Analysis
In Fig. 5(a), the SPP number fraction in mono-polar flocks, n MPF , is plotted with changing ω CT up and down, which shows the clear hysteresis. This result suggests that the emergence mono-polar flocking from the homogeneous bi-directional orientation (or vice versa) has the first order transition-like nature. This is consistent with the appearance of thick transition region in Fig. 3.

A few notes regarding density in flocks
In the case that ST and CT are the same directions with each other (filled circles) in The density increases as the absolute value of the mean rotation speed increases. The collision becomes more frequent and stronger as the torque becomes greater. We can regard that the effective attraction between the particles is strong under such conditions. In the case that ST and CT are the opposite directions with each other (filled diamonds in Fig.   5(b)), this feature is remarkable. Therefore, we find the density 2.5 or more. The large value suggests, the particles overlay each other due to the strong collision. Here, the red diamonds for the high-density flocks belong to the island region of monopolar flocks in Fig. 3(b) and (d) because the red diamonds mean ω ST > 1.5 × 10 −2 . The strong collisions in the cluster caused by the CW-ST and the CCW-CT are confirmed in Fig. 5(b). These high-density flocks are clearly distinguished from the low-density homogeneous phase.
Three-particle simulation creasing ω CT on average. In Fig. 6(d), time evolution of the polar order R p (red solid curve), nematic order R n (green dotted curve) and contact numbers (blue broken curve) corresponding to the sample dynamics shown in Fig. 6(b), or Fig. 4(a). The polar-and nematic-order parameters, R p and R n , are defined as and respectively. The summation i=1,2,3 runs over three SPPs labelled by i = 1, 2, 3. We assume that two SPPs are in contact with each other if the distance of those two SPPs is smaller than the interaction range r, and contact number is defined as the number of such pairs in contact. Figure 6(d) also shows the mono-polar flock and bidirectional orientation regimes. Here, the mono-polar flock regime has been defined as the time window during which R p > R Th and the contact number is 2 or 3. The bidirectional orientation regime has been defined as the time window during which R n > R 2 Th , R p < R Th and the contact number is equal to or higher than 1. (We here set R Th = 0.9 again.) Although there is another short bidirectional orientation regime around t = 100 in Fig. 6(d), we skipped plotting it for better visibility. In Fig. 4(b), we plotted the probability P (t<128) by which the three SPPs form the mono-polar flock at least a time by t = 128 for various ω CT and ω ST . To define P (t<128) for each parameter set, we have counted the number of the samples which has mono-polar flock regimes at least a time until t = 128 and divided it by the total number of samples (128 samples).

The case with isotropic mobility
As expected by the explained mechanism, the formation of mono-polar flocks by CT is not relying on the anisotropic setting of friction in Eq. (1) or (6). In Fig. 7, we performed the simulations with the isotropic setting, R ζ = 1.0, and indeed obtained the similar snapshots.      reproduced while typical morphology of each mono-polar flock seems different from that seen in the anisotropic-friction case, where the flock seems to be more elongated (Fig. 2).