Symmetry and stability of orientationally ordered collective motions of self-propelled, semiflexible filaments

Ordered, collective motions commonly arise spontaneously in systems of many interacting, active units, ranging from cellular tissues and bacterial colonies to self-propelled colloids and animal flocks. Active phases are especially rich when the active units are sufficiently anisotropic to produce liquid crystalline order and thus active nematic phenomena, with important biophysical examples provided by cytoskeletal filaments including microtubules and actin. Gliding assay experiments have provided a testbed to study the collective motions of these cytoskeletal filaments and unlocked diverse collective active phases, including states with long-range orientational order. However, it is not well-understood how such long-range order emerges from the interplay of passive and active aligning mechanisms. We use Brownian dynamics simulations to study the collective motions of semiflexible filaments that self-propel in quasi-two dimensions, in order to gain insights into the aligning mechanisms at work in these gliding assay systems. We find that, without aligning torques in the microscopic model, long-range orientational order can only be achieved when the filaments are able to overlap. The symmetry (nematic or polar) of the long-range order that first emerges is shown to depend on the energy cost of filament overlap and on filament flexibility. However, our model also predicts that a long-range-ordered active nematic state is merely transient, whereas long-range polar order is the only active dynamical steady state in systems with finite filament rigidity.


INTRODUCTION
Active matter systems are composed of a large number of interacting "agents" that individually convert internal energy into motion.Flocks of birds, herds of animals, and growing or swarming bacterial colonies are some examples of active matter across various length scales [1][2][3][4].The formation of spontaneously ordered collective motion is a common phenomenon exhibited by active systems.Clustering and flocking ordered states seen in many active systems can be modeled using active Brownian particles, in cases where the active agents are spherically symmetric [5][6][7].But active matter also comprises anisotropic systems with agents whose aspect ratios differ strongly from unity, such as swimming rod-shaped bacillus subtilis [8] or E. coli bacterial swarms [9,10], glioma cells in brain tumors [11,12], and cytoskeletal filaments such as actin and microtubules [13][14][15].The anisotropy of these individual active particles is responsible for the orientational order and the distinct behaviors of active nematics [16,17].
Collective phases of anisotropic cytoskeletal biofilaments have been experimentally observed in gliding assay systems, which were originally used to study molecular motors such as kinesin, dynein, and myosin, which walk along cytoskeletal filaments such as actin and microtubules [18][19][20].In a conventional in vitro gliding assay setup, molecular motors are uniformly dispersed on a glass slide where they are held fixed with their tail end, which attaches to molecular cargo in vivo, here attached to the surface of the glass slide.Cytoskeletal filaments glide on these motors, which push the filaments forward * d.a.beller@jhu.edu in a defined direction, making the system microscopically polar.When the densities of motors and filaments are sufficiently high in a gliding assay, diverse active collective phases spontaneously emerge, including polar flocks, nematic and polar lanes, asters, and "spools" of filaments [21][22][23][24].One such active state observed is the active nematic state with long-range orientational order.Understanding how this long-range ordered (LRO) nematic state emerges and evolves from microscopically polar, flexible constituents is the main subject of this paper.
The effects of shape on interacting colloidal particles were first studied in thermal systems of passive, hard, rod-like colloidal particle systems to reveal that higher aspect ratio and higher packing fraction favor long-range nematic order [25].These trends hold also in active systems of self-propelled hard rods; additionally, increasing the self-propulsion strength favors nematic alignment [26][27][28][29].
In gliding assays, however, the hard-rod approximation becomes unsuitable due to large deformations of the active filaments.The effects of polymer flexibility on nematic ordering are notable already in passive, lyotropic liquid crystals, where lower rigidity disfavors order and increases the isotropic-nematic phase transition concentration, compared to hard-rod colloids [30].Rigidity is likewise an important parameter in determining the active phases of self-propelled, semiflexible filaments [24,[31][32][33][34][35].
Past simulation studies of active semiflexible filaments that have reported LRO states [24,31,34] employed an attractive or aligning interaction between nearby active filaments in their microscopic model.These assumptions were justified by the use of depletants, which cause filaments to bundle in parallel or antiparallel orientations, in conventional gliding assay experiments where LRO states have been observed [21,22,24,34].However, recent ex-periments with microtubules on lipid substrates have reported LRO states without the need for depletants [36].This experimental evidence motivates us to study how LRO active nematic states emerge and evolve without imposed aligning or attractive interactions in the microscopic model.
Another important feature of experimental gliding assays is that filaments frequently cross over one another, with one filament partially and temporarily delaminating from the substrate [21,24,[34][35][36].These crossovers are forbidden by models that assume strictly two-dimensional motion of impenetrable filaments.As we show in this work, filament crossovers are essential to LRO states in the absence of aligning or attractive interactions at the filament-pair level.
In this paper, we conduct a detailed study of the active phase behaviors of self-propelled, semiflexible filaments without imposed attractive or aligning interactions between the active agents.We determine the conditions under which these filaments develop collective motions with long-range orientational order, by varying parameters including the filaments' flexibility and area fraction.
We find that a long-range orientationally ordered phase does not arise when the filaments are strictly constrained to two dimensions through volume exclusion.Instead, we observe LRO states only in quasi -two dimensional systems, in which we replace the volume-excluding steric repulsion with a finite energy cost for filament overlap [33,37], in order to account for crossovers.
In our simulations, the symmetry of the long-range orientational order that first emerges may be polar or nematic (apolar) depending on the flexibility of the filaments.Surprisingly, we find the LRO active nematic state to be merely transient, with all LRO systems eventually becoming polar.Our model predicts that the time taken by the system to reach the LRO active polar steady state monotonically increases with filament stiffness.These findings suggest that experimental observations of LRO active nematic states in microtubule gliding assays may be better interpreted as transient behavior, observable due to the finite lifetime of experiments and relatively high rigidity of the filaments.

II. MODEL
We use Brownian dynamics to simulate active (selfpropelled), semiflexible filaments in two dimensions.Each filament is modeled as a bead-spring chain with a Hookean restoring force between adjacent beads and a bending stiffness κ.Each bead self-propels with a constant active force in the chain's local tangent direction, resulting in an active speed v 0 in the absence of other forces.The local tangent direction is calculated at each time step from bead positions and does not have its own relaxation time in this model.
Steric repulsion between beads of different chains is modeled using a modified Weeks-Chandler-Andersen po-tential [32,38]: Here, r ij is the distance between the centers of two beads i and j, and r shift is a parameter introduced to make the U WCA potential finite when r ij = 0; this ensures that chains can interpenetrate.The overlap of chains in our explicitly two-dimensional model mimics the crossover events of a quasi-two-dimensional microtubule gliding assay experimental system.In addition to acting between beads of different chains, U WCA also acts between bead pairs on the same chain separated by more than three bead positions, penalizing self-intersections of the semiflexible chains.
We scale simulation lengths in units of x = 2 1/6 σ, the range of the WCA repulsive potential.In these units, the bead radius is r 0 = 0.5x.The chain length ℓ = 31r 0 is the same for all chains and is kept fixed across all simulation results presented here; polydispersity will be investigated in future work.Time t is scaled with respect to x/v 0 , although in our plots we present time t in units of ℓ/v 0 .The Brownian (over-damped) dynamics updates bead positions through The third term on the right-hand side of Eq. 2 is the stochastic noise term representing random fluctuations in bead positions, with each vector component drawn from a zero-mean Gaussian white noise distribution.We use a second-order stochastic Runge-Kutta (SRK-2) method [39] to update the bead positions, with an adaptive time step.All simulations are conducted in a square box of side length L x with periodic boundary conditions in both the x and y directions.Full details of the bead-springchain model can be found in Appendix A.

A. Dependence of active state on chain penetrability
We begin by examining the strictly two-dimensional scenario in which crossovers are prohibited, with r shift = 0.This choice makes the repulsion potential U WCA (r ij ) diverge when two chains overlap (r ij = 0), thus forbidding chains to interpenetrate when they collide.We vary bending rigidity κ, which in passive systems would favor nematic order as it is increased [30].However, in this active system, we find by varying κ across three orders of magnitude that the system does not exhibit long-range order.Instead, polar flocks are seen for all values of κ, as shown in Fig. 1A and in Supplemental Video 1 [40].
These results indicate that long-range orientational order for mutually impenetrable active semiflexible chains cannot be achieved without an imposed alignment or attractive potential.Active polar flocks similar to those observed here have been reported in various microtubule and actin filament gliding assays both in experiments using depletants and in simulations where the filaments were strongly confined to two dimensions [14,29,33,34,41,42].Similar polar flocks have also been reported in experiments on rod-shaped E. coli bacterial colonies when they were strongly confined to two dimensions [37].
We now introduce chain penetrability by modifying the r ij -dependence of the U WCA interaction potential as shown in Fig. 1B, with r shift fixed at 0.1, so that the potential is finite at r ij = 0. We vary ε to change the potential barrier for chain interpenetration, which affects the frequency of crossovers (Supplemental Video 2 [40]).
At the lowest ε values tested, such as ε = 10 −12 as shown in Fig. 2A, we obtain an isotropic state, as chains simply cross over each other essentially unimpeded.When ε = 10 −10 we observe polar density waves (Fig. 2B) traveling parallel to the chain self-propulsion direction, as previously seen in simulations and experiments involving actin filaments [24] (where the term "polar clusters" is applied to this state).A state with no positional order but long-range orientational order (Fig. 2C), our main interest in this work, is obtained when ε = 10 −8 .We therefore maintain U WCA parameters of ε = 10 −8 , r shift = 0.1 in the following sections of this paper.
At even higher values of ε, we obtain density inhomogeneities in the form of flocks.When ε = 10 −6 , a jammed flocking state appears (Fig. 2D), with characteristics dis-tinct from the polar flocks seen in Fig. 1A.The jammed flocks seen here are formed by counter-moving chains that are dynamically arrested relative to the center of mass of the flock, which then moves as a whole.(In contrast, all the chains in a particular polar flock seen in Fig. 1A are oriented nearly in the forward direction.)The jammed flocking state of (Fig. 2D) can be interpreted as a form of motility-induced phase separation [6], with countermoving filaments impeding each others' motions as they overlap.
The density inhomogeneities seen in Fig. 2E and F are dynamic and they reconfigure into different lanes and flocks multiple times in the course of the simulation.While the polar lane in Fig. 2F spans the system, we observe upon increasing the system size a state like the finite flock of Fig. 2E.We therefore interpret the system-spanning lane of Fig. 2F as a finite-size effect.

B. Long-range-ordered states
With ε fixed at 10 −8 henceforth, we now investigate active states with long-range orientational order (LRO) in greater detail.First, we study the effect of changing bending rigidity κ on the form of the orientational order.In Fig. 3, we see snapshots of the initially formed long-range orientational order (time t * = 200.58[ℓ/v0 ]) for a range of κ values.To quantify the orientational order we compute two different order parameters: nematic order S( t) = exp i2θ t j j and polar order , where θ is the average orientation of all the beads in a chain with respect to the positive xaxis, and subscripts j indicate a system average over all chains.We label a state as orientationally ordered when S( t) ≥ 0.  In Fig. 3, we see that more rigid chains (greater κ) form a LRO nematic state (S(t * ) > 0.7 and P (t * ) < 0.7) whereas more flexible chains achieve a LRO polar state (P (t * ) > 0.7).The same trend is revealed when nematic and polar order parameters, averaged over 20-run ensembles, are plotted against κ (Fig. 4A): whereas ⟨S(t * )⟩ increases modestly with increasing κ, ⟨P (t * )⟩ decreases dramatically.Biopolymers such as microtubules have a higher bending rigidity (persistence length ≈ 1 mm) compared to actin filaments (persistence length ≈ 10-20 µm) [43].While LRO polar states have not been observed in experiment, our results are consistent with the known tendency of actin gliding assays to exhibit short-range polar order whereas nematic order (short-or long-range) is more common in microtubules [14,21,24,34,36,41].
The modest increase in nematic order with increasing rigidity can be attributed to the fact that more rigid chains have greater restoring forces against bending deformations, and are therefore less likely to be significantly deflected from their direction of motion due to collisions or random fluctuations.Box plots in Fig. 4B and C show the median and spread of nematic and polar order parameter values at time t * , when all the simulations have reached long-range orientational order.It is interesting to note that, at larger κ values that gave apparently nematic states, the polar order parameter values (Fig. 4C) are significantly above zero and widely spread.This indicates that LRO nematic order coexists with LRO polar order, and raises the question of whether the LRO nematic state is, in fact, the dynamical steady state.

C. Nature of the active steady state
In order to determine whether LRO nematic order with weak LRO polar order is a dynamically stable active state, we now examine the evolution of orientationally ordered systems over times much longer than t * .Examining the time-evolution represented by the snapshots in Fig. 5A and Supplemental Video 3 [40], we can see that the LRO active nematic state discussed in the previous section is merely transient, and the system reaches an active global polar steady state at longer times.To quantify this evolution of orientational order, we plot in Fig. 5B the polar and nematic order parameters against time for each simulation in a 20-run ensemble, with bending rigidity κ = 70.The black dashed line indicates the time t * when all the simulations have achieved long-range nematic orientational order.However, the nematic (apolar) state is transient, as indicated by P ( t) rising to saturation near 1 at much later times.The active steady state is therefore a LRO active polar state.At intermediate times after t * , the polar order appears to plateau for a significant time interval, though the plateau value of P ( t) exhibits wide variation between simulation runs.
To understand the effect of changing bending rigidity on the steady-state long-range order, we plot in Fig. 6 the ensemble-averaged nematic (Fig. 6A) and polar (Fig. 6B,  C) order parameters against time for multiple values of κ.From Fig. 6A we see that nematic order saturates at around the same time t * for all κ values.Fig. 6B shows that, while polar order eventually saturates for all κ values, the time of saturation depends strongly on κ.Systems with larger κ tend to take much longer times to reach their LRO active polar steady state.To examine the κ-dependence of the P ( t) saturation time, we plot in Fig. 6C the same data as in Fig. 6B but with a linear scale on the time axis.Fitting a saturating exponential function f ( t) = r−ae − t/τ to the ensemble-averaged polar order parameter yields typical time scales τ for saturation of polar order at P ( t) → r.As shown in the inset to Fig. 6C, we find that τ monotonically increases with increasing κ, quantifying the slower polar ordering exhibited by more rigid filaments.The saturation value r varies between simulation runs, and its ensemble average increases slightly with increasing κ, indicating better polar ordering in more rigid filaments.
We next investigate the effect of changing the area fraction ϕ of filaments on the LRO state (Appendix B Fig. 8, 9, 10).We see that the lowest area fraction tested (ϕ = 0.1) takes the longest to reach the long-range orientational order (Fig. 9A).When we simulate the sys- tems for much longer times compared to t * , and fit the ensemble-averaged polar order parameter to a saturating exponential as we did in Fig. 6C, we find that the time scale τ required to reach a LRO polar steady state monotonically decreases with increasing ϕ (Fig. 9C).This effect can be attributed to an increase in the rate of filament interactions for higher ϕ-values, giving the system more opportunities for reorientation and hence reaching the steady state sooner.
The wide range of plateau P -values observed at the time t * of nematic ordering poses a challenge for isolating the effect of ϕ on the timescale for polar ordering.To control for different histories of partial polar ordering, we conducted a set of simulations in which the filaments are initialized in one of two opposing directions, so that the initial global nematic order is perfect, S 0 = 1, and the initial global polar order P 0 is precisely specified (Fig. 10).A polar order saturation time T Pc was defined as T Pc = tPc + τ Pc , where P c is a chosen threshold P -value, tPc is the time at which the increasing P ( t) reaches P c , and τ Pc is the time constant from an exponential fit to P ( t) for all later times.By calculating T Pc for an ensemble of simulations over a range of P 0 and ϕ values, we find that systems with higher area fraction take less time to reach LRO polar order for each P 0 (Fig. 11.) Furthermore, the exponential saturation time constant τ Pc is remarkably consistent across a wide range of P 0 values (excluding cases of P 0 > P c ) at each value of ϕ (Fig. 12).This latter result indicates that the late-time polar ordering dynamics have no dependence on the history of nematic ordering other than through a t-shift that depends on the plateau value P ( t * ).These findings are consistent over a range of choices of cutoff value P c for the exponential fits.
Changing the system size L x does not seem to have any significant effect on the typical time scales needed to reach long-range orientational order (t * ) or to reach the LRO active polar steady state (τ ) (Appendix D Fig. 14).The plateau values of the polar order seen in Fig. 6B are also seen not to depend significantly on ϕ and L x .We additionally note that incorporating anisotropic friction into the Brownian dynamics model [33] does not qualitatively change our results, as seen in Appendix C Fig. 13.

IV. DISCUSSION AND CONCLUSIONS
In this study, our objective was to investigate ordered states in the collective motions of self-propelled, semiflexible filaments.With a microscopic model lacking attractive or aligning mechanisms between filament pairs, we find that crossovers are essential to the emergence of LRO states.In contrast, strictly two-dimensional confinement permits only locally ordered polar flocks [14,29,33,34,41,42].This finding is consistent with a recent experimental study of E. coli bacterial colonies on a surface, which are seen to form LRO nematic states when crossovers are allowed, and polar flocks when crossovers are suppressed [37].
When we allow filaments to interpenetrate in order to model crossovers, the symmetry of the LRO state that initially forms is nematic for more rigid filaments and polar for more flexible ones.However, we find that the LRO nematic state, despite its experimental importance, is only a transient state in our simulations.All LRO states eventually evolve into LRO polar states, which we find to be the unique dynamical steady state.The typical time required for emergence of the polar steady state increases monotonically with filament rigidity.
Our results suggest that experimental realizations of LRO active nematic states in microtubule gliding assays are long-lived transient states rather than dynamical steady states.Due to factors such as ATP consumption, gliding assay experiments necessarily have finite lifetimes.For relatively rigid filaments such as microtubules, our model predicts that the timescale τ for emergence of polar order (and thus, the persistence of an apparently nematic state) may be much longer than the time t ∼ t * at which nematic order emerges from an isotropic initial condition.Comparing the typical length and active velocities of microtubules in gliding assay experiments [36] to our simulation, we estimate that microtubule bending rigidity corresponds to κ ≈ 250 (Appendix E), which is certainly in the high-κ regime according to our results.
Our results also suggest that in the limit of κ → ∞ (active rigid rods) the timescale for polar ordering may diverge, making the active nematic state dynamically stable [26][27][28][29].Long-lived, long-range nematic order in E. coli colonies might be interpreted as an effectively stable state in a system of effectively rigid rods [37].
It remains an open question how strongly our findings depend on the use of periodic boundary conditions and reduced system size compared to typical experiments.Surprisingly, however, we find that changing the system size does not have any significant effect on the typical timescale for emergence of the polar steady state.This suggests that our simulations may be capturing the true physics of the bulk system.We also find that our choice of isotropic friction, at the level of individual beads, is not important to our findings, as introducing anisotropic friction does not qualitatively change our finding that the LRO polar order is the steady state (Appendix C Fig. 13).
A number of interesting questions are raised by the plateaus in polar order P ( t) that characterize the transient LRO nematic state in systems of more rigid filaments, before the late-time emergence of the LRO polar steady state (Fig. 5B and 6B).These plateaus are clear in P ( t) plots from individual simulation runs (Fig. 5B), even though the ensemble-averaged ⟨P ( t)⟩ appears to be monotonically increasing.We have confirmed that the Pvalue of this plateau has no significant dependence on the system size or on the area fraction of filaments (Fig. 9,14).The reason for this plateau in the polar order parameter's evolution, and the dependence of the plateau value on κ, remain open questions.
The variability that we observe in plateau P -values suggests that transient states may have wide variability in their orientational order, in an ensemble of active systems with identical parameters.As nematic order S saturates at earlier times, the polar order P that happens to develop simultaneously is apparently remembered over the long lifetime of the transient LRO nematic state.Increasing the strength of noise in the simulations may help to erase memory of the details of the nematic ordering process and thus reduce differences in order between transient states of different runs, resulting in a narrower distribution of P values.
In this work, we have only focused on the long-range orientationally ordered states with no density inhomogeneities and established that the LRO nematic state is only a transient state.But gliding assays also spontaneously form a variety of ordered states with non-uniform densities [23,24,33,34,36].Our results highlight the potential importance of long-lived transients in active states generally, so it is of interest to investigate whether distinct active behaviors may be found in the time-evolution of other, inhomogeneous active states.
A full understanding of LRO active states, including the role of filament bending rigidity, requires detailed study of binary collisions [24,34,36], which is the subject of ongoing work.This understanding will complement our study in this work of the effect of area fraction ϕ.We find that increasing ϕ reduces the time required for emergence of the LRO polar steady state.This effect can be attributed to an increase in the rate of collisions with increasing ϕ, giving the system more opportunities to reorient filaments.Other realistic complications of interest, including polydispersity, chirality, and growth, make the active states of self-propelled, semiflexible filaments a subject of rich phenomenology even in the simplest ordered phases.
is the angle change between points separated by a small distance ∆s along the contour and κ = 1 2 l p k B T is the bending stiffness where l p is the persistence length of the chain [44].We incorporate this bending energy as a potential acting between every set of three consecutive beads (∆s ≈ 2r 0 ) that penalizes the bending of the chain.
κ = κ ∆s has the units of energy and will be referred to as the bending stiffness hereafter and θ is the bond angle between two adjacent linear springs (Fig. 7).The total force excluding viscous drag is given by where U WCA is given by Eq. 1.The Brownian (overdamped) dynamics [33] updates bead positions through Here v is the velocity of the bead, γ is the drag coefficient, and ξ(t) is the stochastic noise term representing random fluctuations in bead positions, drawn from a zero-mean Gaussian white noise distribution.

Appendix B EFFECT OF AREA FRACTION
We investigated the effect of changing the area fraction of filaments on the LRO order.Fig. 8 shows snapshots of simulations with different area fractions after they have reached orientational order.All the simulations shown here have filaments with κ = 10 and U WCA parameters ϵ = 10 −8 and r shift = 0.1 (Fig. 2C).From these snapshots, we see that the symmetry of the LRO state reached for all the area fractions is polar which is consistent with the polar order seen for κ = 10 in (Fig. 3) where ϕ = 0.3.So increasing the area fraction to ϕ = 0.5 or ϕ = 0.7 does not change the type of orientational order initially reached.
To understand the steady state behavior of the system, we examined the order of the system after a long time compared to t * .The global nematic order parameter in Fig. 9A reveals that all the simulations reach an orientational order but the system with the lowest area fraction (ϕ = 0.1) takes the longest time to reach the global orientational order.This is presumably due to filament reorientation interactions being less frequent in the lower-density system.The LRO polar state is still the steady state of the system as seen from Fig. 9B.The system takes less time to get to the LRO polar state when the area fraction is higher.Fitting a saturating exponential function f ( t) = r − ae − t/τ to the mean polar order parameter, we see that τ decreases with increasing area fraction (Fig. 9C lower inset) whereas saturation parameter r is not seen to follow any trend with change in area fraction (Fig. 9C upper inset).
To better understand the dependence of saturation time (from LRO nematic to LRO polar) on ϕ separately from its dependence on the plateau value of P ( t), we set up systems with varying area fraction ϕ having all filaments oriented initially along one of two opposing directions.These systems have perfect initial nematic order S 0 = 1 and controlled polar order P 0 .Fig. 10 shows the time evolution of nematic and polar order in these systems for two example values of initial polar order, P 0 = 0.1 and P 0 = 0.7.Unlike the ensembleaveraged P ( t) examined in Fig. 6, P ( t) data for individual runs in these controlled simulations are not generally well-fit by a saturating exponential function over the entire simulation time, especially for low P 0 values such as P 0 = 0.1.However, because the late-time behavior is well-described by a saturating exponential, we can still find a characteristic time scale by defining the total time T Pc to saturation of polar order as T Pc = tPc + τ Pc where tPc is the time taken by the rising polar order P ( t) of each simulation to reach a specific cutoff value P c .For time t > tPc we fit P ( t) to a saturating exponential (f ( t) = r − ae − t/τ P ) to obtain the time constant τ Pc .In Fig. 11A, B and C the cutoff P c = 0.5, P c = 0.7 and P c = 0.9 respectively.Over a wide range of choices of cutoff P c , Fig. 11 shows that a higher area fraction leads to lower saturation time for each initial polar order P 0 .Simulation runs initialized with P 0 > P c are excluded from this analysis for each choice of P c .Fig. 12A, B and C show the fit parameter 1/τ Pc for various values of ϕ with cutoffs at P c = 0.5, 0.7, and 0.9, respectively.It can be seen that the system is memoryless in that the saturation time τ Pc is similar for all P 0 .For example, τ 0.7 is very similar for P 0 = 0.7 (fitting over the entire simulation time-range) as for P 0 = 0.1 (in the time-range after P ( t) reaches 0.7).

Appendix C ANISOTROPIC FRICTION
The model described previously assumes isotropic friction on the beads.But as these beads form long filaments with a high aspect ratio, we examine here whether our results are significantly changed by friction anisotropy.For this purpose, we replace Eq. 8 with an anisotropic generalization for a particular bead α, where M α ij is the anisotropic mobility tensor of the αth bead [33]: Here, γ ⊥ and γ ∥ are the perpendicular and parallel drag coefficients.The symbol ⊗ represents the outer product, and t α is the tangent unit vector to bead α, found by where r α is the separation vector between the centers of (α − 1)th and αth connected beads, and 1 < α < N if there are N beads in a chain.For a bead that is located at the ends, t 1 = r 2 and t N = r N .The diffusion coefficient is also anisotropic and for a wormlike chain, it is given by D α ij = k B T M α ij [33].This can be used in Equation 7to compute the positions of beads at the  Using anisotropic friction of γ ⊥ = 2γ ∥ in a 20-run ensemble of simulations with κ = 30 (Fig. 13), we see no qualitative difference compared to simulations with isotropic friction as both Fig. 13A and B reach an LRO polar steady state after being in a transient nematic state at intermediate times.From ⟨S( t)⟩ values of both the cases, we can see that the case with isotropic friction reaches the LRO nematic state at a slightly earlier time compared to the case with anisotropic friction.However, the case with anisotropic friction is seen to reach the LRO polar steady state sooner.We also find that anisotropic friction reduces differences in polar order between transient states of different runs during intermediate times, resulting in a narrower distribution of P -values.There is also a slight increase in ⟨P ( t)⟩ at these intermediate times compared to an ensemble with isotropic friction.We provide here a comparison of our simulation units to microtubule (MT) gliding assay experiments on lipid membranes where a depletion agent was not needed to facilitate microtubule bundling [36].If we equate our simulation filament length ℓ = 15.5x to the average length of MTs in experiments ≈ 10 µm [36], we can obtain an estimate of how the simulation length unit x corresponds to the experimental length units: x ≈ 0.64 µm.Similarly, the time scaling is set by comparing the active speed of the filaments in our simulation (v 0 = x/ t = 1) and the MT velocity in the experiments ≈ 0.5 µm/s.This implies t ≈ 1.28s.To obtain an estimate of mass scaling, we look to the over-damped Brownian dynamics equation 8: [γ] = m/ t.To estimate the drag coef- For the bending stiffness κ, we obtain a relation to the stiffness of MTs in experiments using the persistence length l p ≈ 10 −3 m of MTs.Since bending energy in our model acts between every set of three consecutive beads, we have ∆s = 2r 0 x (refer to the model details in Appendix A).Bending stiffness κ is related to energy units ε by: Thus the bending stiffness of taxol-stabilized MTs in the gliding assay experiment on a lipid membrane in Ref. [36] roughly corresponds to a value of κ ≈ 250 in our simulations.Another common method of stabilizing microtubules is by using guanosine-5 ′ [(α, β)-methyleno] triphosphate sodium salt (GMPCPP), and this is known to increase the rigidity of microtubules compared to taxol [45].The corresponding rigidity in our simulations is therefore κ > 250, and our results for the high-rigidity regime are expected to apply for both MT stabilization methods.

FIG. 1 .
FIG. 1. A. Polar flocks observed when r shift = 0 and ε = 1 in UWCA.The simulation box size Lx = 250, area fraction ϕ = 0.5, and number of chains N = 2016.The color wheel legend, here and in subsequent figures, indicates the the orientation of each filament, averaged over the tangents at each bead.B. Plot of the WCA potential for bead interactions in the cases of steric repulsion (blue) and penetrable filaments (orange).

FIG. 2 .
FIG. 2. Changing the energy cost ε for chains to interpenetrate we obtain the following states: A. Isotropic state with no positional or orientational order.B. Polar wave state where chains move in the same direction as the wave.C. Long-range orientational ordered state.The orange line in Fig. 1B corresponds to the potential used here.D. Jammed flock state.E-F.Dynamic polar flocks and lanes.All the simulations have box size Lx = 100, ϕ = 0.5, ℓ = 15.5, κ = 30 and r shift = 0.1.The color wheel legend indicates the orientation of each filament, averaged over the tangents at each bead.

FIG. 3 .
FIG. 3. States with long-range orientational order obtained for different bending rigidities κ.All simulations have Lx = 100, ϕ = 0.3, ℓ = 15.5, ε = 10 −8 and r shift = 0.1.S(t * ) and P (t * ) are the nematic and polar order parameters, respectively, of the simulation runs pictured here, at time t * roughly corresponding to the saturation time for nematic order.

FIG. 7 .
FIG. 7. Schematic illustration of forces in the active bead-spring chain model.(Left) Illustration of active selfpropulsion and WCA repulsion between pairs of chains.(Right) A zoomed-in view of the region enclosed by the dashed lines at left, to illustrate the stretching and bending spring forces.

FIG. 8 .
FIG.8.Changing area fraction ϕ of filaments.All the snapshots are taken at a time when the system has reached orientational order, which is polar for all cases shown here.The filaments have κ = 10 and the system size Lx = 250.

FIG. 10 .
FIG. 10.Time evolution of a system of κ = 30 filaments initialized with perfect nematic order and with specified values of initial polar order P0, for different area fractions ϕ. A. Nematic order time-evolution for P0 = 0.1.Inset shows a snapshot of the initial condition with filaments oriented either up (purple) or down (green) for an area fraction of ϕ = 0.3.B. Polar order time-evolution for P0 = 0.1.C. Polar order time-evolution for P0 = 0.7.Inset shows the nematic time-evolution and a snapshot of an initial condition with (ϕ = 0.3).The ensemble size is n = 10 for ϕ = 0.1, 0.3 and n = 3 ϕ = 0.5, 0.7, for each value of P0.