## Abstract

Motor proteins drive persistent motion and self-organization of cytoskeletal filaments. However, state-of-the-art microscopy techniques and continuum modeling approaches focus on large length and time scales. Here, we perform component-based computer simulations of polar filaments and molecular motors linking microscopic interactions and activity to self-organization and dynamics from the filament level up to the mesoscopic domain level. Dynamic filament cross-linking and sliding and excluded-volume interactions promote formation of bundles at small densities and of active polar nematics at high densities. A buckling-type instability sets the size of polar domains and the density of topological defects. We predict a universal scaling of the active diffusion coefficient and the domain size with activity, and its dependence on parameters like motor concentration and filament persistence length. Our results provide a microscopic understanding of cytoplasmic streaming in cells and help to develop design strategies for novel engineered active materials.

## INTRODUCTION

Active systems are driven by nonthermal, energy-consuming processes that lead to very rich collective behavior (*1*–*4*). Examples are found from macroscopic length scales (fish and birds) to microscopic length scales (bacteria and algae) to even subcellular structures like the cell’s cytoskeleton. The cytoskeleton is composed of filaments that are dynamically interconnected by passive and active cross-linkers. It provides mechanical stability to biological cells, acts as a force-generating element, and serves as a track network for active intracellular transport (*5*, *6*). Moreover, its dynamics generates internal motion of the cytoplasm that secures nutrient availability and the distribution of organelles (*7*). Fundamental knowledge about the relationship between cytoskeleton structure and dynamics and about the molecular driving forces helps to obtain a deeper understanding of cellular function and dysfunction in vivo and to design active gel materials, e.g., artificial cells (*8*) in vitro. Experimental studies of purified cell extracts containing cytoskeletal filaments and molecular motors show a plethora of dynamical phenomena. The long, thin, and stiff filaments self-organize into various liquid crystalline structures depending on filament concentration, motor type, and presence of passive cross-linkers. In vitro experiments on microtubule-kinesin and actin-myosin mixtures, where polymer-induced depletion interaction drives bundling of filaments and confines them at an oil-water interface, have shown exciting nonequilibrium behaviors, such as persistent spontaneous flows, turbulent-like motion, and the formation of highly mobile +1/2 topological defects, a characteristic signature of active nematics (*9*–*13*).

Various theoretical modeling approaches (*14*–*22*) have been used to reproduce and explain the nonequilibrium features of cytoskeletal filament-motor mixtures. Continuum descriptions based on active gel theory have been amazingly successful to capture essential aspects on the mesoscopic scale, particularly the formation and dynamics of topological defects (*14*, *15*, *17*–*20*). In this approach, activity is incorporated into a model of passive nematics by an active current **J** = ζ∇ · **Q**, where **Q** is the nematic tensor and ζ is the activity strength. This active current originates from the active force dipoles, which generate a particle flux along or against the nematic curvature (*23*, *24*). For coarse-grained descriptions of active nematics at the filament level, nematic symmetry and activity are difficult to reconcile. Models of rods expanding in length and breaking in periodic intervals represent growing filament bundles (*25*) and have been used, as well as models of filaments that intermittently move forward and backward (*26*). At the microscopic level of filaments and motor proteins, the system is not nematic but polar. Studies of concentrated systems modeling polar filaments and molecular motors show the formation of polar bands and persistent motion generated at the boundary between the bands, where the filaments are antialigned (*16*, *27*–*30*). The generation of extensile stresses occurs predominantly at the interfaces between polar domains and can drive the instability of a global nematic phase. This has been shown by a multiscale modeling approach that uses filament-based computer simulations of rods to provide parameters for a continuum polar-fluid model (*16*, *28*, *29*).

An important question arising from all the previous studies is the role of hydrodynamic interactions. It has been shown that systems with extensile force dipoles (corresponding to pusher-type hydrodynamics) destabilize the nematic order, whereas contractile force dipoles (corresponding to puller-type hydrodynamics) stabilize it (*1*, *19*, *23*, *24*). These hydrodynamic stresses are at the origin of the active current in active nematic theory. However, in the absence of hydrodynamic interactions, active forces can also render the nematic phase unstable, like in experiments (*31*) and simulations (*32*) of granular ellipsoids and in a phenomenological continuum model for overdamped active nematics (*33*). These different approaches raise several pertinent questions: What is the relevance of polarity in active gels? Are the individual filaments the appropriate fundamental units of modeling, or should it rather be the filament bundles? Are the predictions of active gel theory and of polar active filament descriptions consistent? Which advantages can a more microscopic approach provide? What is the role of filament flexibility?

Here, we study the emergent structures, persistent motion, and instability of the aligned nematic phase in mixtures of polar semiflexible filaments and molecular motors. Two-dimensional Langevin dynamics simulations are designed to mimic experimental systems of filament-motor suspensions confined at oil-water interfaces (*9*, *12*, *13*). Our modeling approach bridges length scales from nanometer-sized molecular motors to micrometer-long semiflexible filaments. The results presented here show that the average motor-induced force on antiparallel filaments is a robust measure for the activity in the system and that hydrodynamic interactions are not necessary to describe the system dynamics. In an initially globally nematic filament suspension, the generated active stresses induce filament sorting in polar bands, a buckling-type instability, the formation of nematic-type topological defects, and—at steady state—the emergence of complex flow patterns. The dynamics of individual filaments is characterized by active Brownian particle–like motion. The active-force dependence of several inherent length scales, like the domain size, the length of quasi-ballistic motion, and the effective persistence length of single filaments, all follow a universal inverse square root dependence. Furthermore, we show how active nematic behavior emerges with increasing filament concentration. Where applicable, the predictions from our “polar active nematics” description agree well with those of active nematic theory.

The simulations presented here are based on a model of semiflexible filaments of contour length *L* and persistence length ℓ_{p} in two dimensions. The filaments consist of *n*_{s} beads of diameter σ connected by stiff harmonic springs with rest length *a*_{0} = *L*/(*n*_{s} − 1), where *a*_{0} = σ. The filament area fraction *n*_{f}, or the box size *L*_{box}. Molecular motors are modeled by harmonic springs with rest length *r*_{0} and spring constant *k*_{m}. They attach to neighboring filaments with rate Γ_{att}. Motors walk in the direction of the filament polarity with step length *a*_{0}. The step rate is proportional to the probability *t* is the motor time step. Motors detach when they reach the end of a filament, when they encounter a motor already bound, or when their length exceeds a threshold *r*_{off}. In the following, we use rescaled quantities: filament aspect ratio *k*_{B}*T* is the thermal energy), filament friction *t* is the simulation time step), box size *D*_{0}. For details, see Materials and Methods.

## RESULTS

### Activity drives the formation of polar domains

The microscopic origin of filament motion is the dynamic cross-linking of *n*_{f} polar filaments by *n*_{m} molecular motors. The motors walk in the direction of polarity **p** on the filaments and induce sliding forces in direction −**p**. If two filaments are polar aligned and two consecutive motor steps occur on the two different filaments, the motors induce no net filament motion (see Fig. 1A). However, if the filaments are antialigned, the motors get stretched and net filament motion results (see Fig. 1B). The *n*_{m} motors are classified as

When motors are added to an initially globally nematic suspension of filaments with random left-right filament orientation, the motor-induced sliding forces first lead to a rapid sorting of filaments into narrow polar bands with antiparallel alignment at the domain boundaries (see Fig. 1C). In time, these polar bands coarsen. When the activity is large enough, a buckling instability leads to disordered configurations of polar domains with topological defects (see Fig. 1D). If the activity is too small or the filaments are too stiff, the steady state consists of stable parallel bands (see fig. S1). Movie S1 illustrates the polarity-sorting process in more detail. Movie S2 shows the entire dynamical evolution from the initial nematic state to the stationary state. Because in this system both nematic alignment and polar sorting play an important role, we call it an “active polar nematic.”

Steady-state filament dynamics emerges from a complex interplay between various filament and motor properties. For example, the number of antiparallel motors is an important factor driving the filament dynamics and the suspension structure. We quantify the activity by the motor force per filament*r*_{m}/*r*_{0} − 1 of the antiparallel motors. In steady state, the motor force is balanced by an effective medium friction with friction coefficient *v*_{∥} is defined as *v*_{∥} = lim_{τ → 0}〈 − **p**(*t*) · (**r**(*t* + τ) − **r**(*t*))/τ〉* _{t}*, where

**r**(

*t*) and

**p**(

*t*) are the position and the unit tangent vector of the filament center, respectively, and τ is the lag time (see fig. S2). Note that the active force

*n*

_{m}, the persistence length ℓ

_{p}, and the motor stepping rate

*v*

_{∥}

*L*/

*D*

_{0}increase linearly with the activity in the system

_{p}(see fig. S3C) but does depend on the width of the bands (see fig. S3D). For wide bands, the velocity decays to zero over a distance between one and two filament lengths as is also confirmed by an explicit calculation of the velocity correlation length from the velocity correlation function as outlined in figs. S6 to S8.

### Single-filament motion and active diffusion

We first characterize single-filament motion in dense systems in the active steady states by calculating the filament orientational autocorrelation functions and filament center-of-mass mean squared displacements (MSDs). The latter can be measured experimentally by fluorescently labeled tracer filaments or by using tracer particles that can be tracked. For nonzero *D*_{act} that increases with increasing _{R}, which also depends on *n*_{m}. The active rotation time τ_{R} is shown in Fig. 2B both as a function of _{R} decreases with increasing *34*). This implies

### Buckling polar bands

The strength of the active force does determine not only the average filament velocity but also structure, size, and stability of the domains. Polar bands are stable for weak active forces and large persistence lengths, whereas they become unstable and buckle at a particular wavelength λ for sufficiently large active forces. The time evolution of the motor-filament mixture in Fig. 3 (A to C) shows first polarity-sorted bands and then progressive bending and breaking of bands (see also movie S3). The wavelength λ of the instability displays a square root dependence on the filament persistence length ℓ_{p} (see Fig. 3D). The assumption that the active force in all these systems evolves independent of the persistence length (see fig. S5) suggests an Euler buckling–type instability with buckling force *F*_{b} ∝ *E*/λ^{2}, where *E* ∝ ℓ_{p} is the (effective) elastic modulus (*35*). This relation applies under the assumption that the effective elasticity *E* is only weakly affected by the activity. Together with

### Dynamics of topological defects

The system of polar filaments and motor proteins shares features with both active nematics and polar active fluids, although it is clearly distinguishable from both of them. The high densities favor nematic ordering for which the polarity is not important, whereas the two-filament pair interaction mediated by the motors has a pronounced polar character. The latter leads to local polar order in the domains. However, the characteristic topological defects that appear in polar active fluids, +1 or −1 defects (*14*), are never observed here. Instead, the defect structures in the dynamic disordered phase are *C*_{1} symmetry, and the other has *C*_{3} symmetry (see Fig. 4, B and C). Note that for the *C*_{1} defect in Fig. 4B, there are additional (unmarked) boundaries between antiparallel domains in the lower left and right corners where the motion is originated.

The f motion of the *l*_{def} between defects on the force, as shown in Fig. 4E. A similar scaling is found for other related length scales in the system. The exponential decay of the parallel and perpendicular spatial orientational correlation functions *l*_{⊥} and *l*_{∥}, which are estimations for width and length of the domains. Last, *v*_{∥}τ_{R} is proportional to the distance that a filament moves along the domain boundary, before changing direction. The activity dependence of these three lengths is shown in Fig. 4F and follows the same scaling as *l*_{def}. Note that both τ_{R} and *v*_{∥} are averaged values over a distribution of faster (antiparallel) and slower (parallel) filaments. Although these three lengths provide different quantitative estimations of the domain size, their scaling is consistent with an inverse square root dependence on the activity.

### From dilute to dense filament systems

In passive lyotropic liquid-crystalline systems, particle density and shape determine both structural and dynamical properties. The isotropic-nematic (IN) transition in a system of hard rods with aspect ratio *L*/σ = 20 at ϕ ≈ 0.2 (*36*) is a lower bound for the IN transition for our semiflexible filaments. To investigate the effect of filament density, we show in Fig. 5 (A to D) the simulation snapshots for various filament densities, below and above the IN transition in the passive systems, and different *12*). Although the structures observed strongly depend on ϕ, the (parallel) velocities still show a unique linear scaling with the force density *v*_{∥}τ_{R} show a universal scaling (see Fig. 5G). The chosen scaling variables (*v*_{∥}*L*/*D*_{0}) are independent of the filament length, as indicated by three simulations for filaments twice as long as all other simulations.

Last, we study the effect of activity and concentration on the effective filament stiffness or effective persistence length

A simple argument for this behavior can be extracted from the properties of single tangentially driven self-propelled filaments (*37*). Under the assumption that filament motion is predominantly along their contour (“railway motion”), the rotational diffusion *D*_{R} of the end-to-end tangent vector (the change of the contour) can be calculated explicitly and only depends on the active velocity and the persistence length, *v*_{‖}τ_{R} with the active force (see fig. S11).

## DISCUSSION

Our microscopic model of active gels—based on semiflexible polar filaments cross-linked by molecular motors—shows that for large enough motor concentration, an initially nematic arrangement of filaments is unstable and evolves via polarity sorting, band formation, and buckling into a stationary but highly dynamic polar-domain structure with persistent defect formation and annihilation. We find universal scaling of the domain sizes with the active force determined by the number of antiparallel motors and their extension. The predictions of our model are in good qualitative agreement with recent experiments. While polarity has received little attention for a long time, a recent in vitro experiment provides evidence for the existence of polar bands (*38*). Furthermore, the instability of the polar bands in our model is characterized by a unique length scale, in agreement with recent in vitro observations (*39*). In the disordered phase with polar domains, filaments display a short-time active ballistic motion followed by active diffusion at later times with a diffusion constant that increases with activity, in agreement with results of a model of growing and breaking filaments (*9*). For polar bands, similar active ballistic motion has been found in simulations of stiff rods (*28*).

Although polarity is important in driving the filament dynamics, the formation of *15*). The individual motion of *13*) because it modifies the ratio of bend and splay elastic moduli of the nematic phase and thereby determines the filament structure in the vicinity of the defects (*26*). We have shown that it also has a pronounced effect on the instability, domain size, and defect density. In addition, the effective filament persistence length is found to decrease with increasing activity, in agreement with experiments (*13*) and simulations (*26*).

Our model can readily be extended to three-dimensional systems, to study mixtures of stiff and (semi)flexible filaments, as well as of models with an increased level of complexity like the inclusion of passive cross-linkers or hydrodynamic interactions. Because of the simplicity of the underlying mechanisms, we hope that our results can also help in the design of novel engineered active biomimetic materials. The Supplementary Material accompanies this paper at www.scienceadvances.org/.

## MATERIALS AND METHODS

Suspensions of *n*_{f} semiflexible filaments and *n*_{m} motors are studied using Langevin dynamics (*40*) in two dimensions with periodic boundary conditions. Semiflexible filaments of length contour *L* are modeled as discrete chains of *n*_{s} beads of mass *m* with position vectors *q* = {1. …*n*_{f}} and *i* = {1. …*n*_{s}}. The *n*_{s} beads define *n*_{s} − 1 segments of length *a*_{0} = *L*/(*n*_{s} − 1) with unit tangent vectors *a*_{0} and spring constant **t**(0) · **t**(*s*)⟩ = exp [ − *s*/ℓ_{p}] = exp [ − (*d* − 1)*k*_{B}*Ts*/2κ] with dimensionality *d* and arc length *s* along the contour, such that the bare persistence length is ℓ_{p}/*L* = 2κ/*Lk*_{B}*T* (*41*). Nonbonded inter- and intramolecular excluded volume and attractive interactions are described by a Lennard-Jones potential *V*_{LJ} = ε[(σ/*r*)^{12} − 2(σ/*r*)^{6}], where *r* denotes the inter- or intramolecular bead-bead separation, ε denotes the interaction strength, σ denotes the interaction diameter (i.e., the beads are not overlapping) and *r*_{c} = 3σ denotes the potential cutoff radius. In all our simulations, *k*_{B}*T*/ε = 10, i.e., attractions between filaments are negligible. These parameters lead to an effective Barker-Henderson diameter (*42*), which is 14% less than that for a system with a purely repulsive Weeks-Chandler-Andersen (WCA) potential at ε/*k*_{B}*T* = 1.0 and 4% less than the WCA potential with ε/*k*_{B}*T* = 0.1. The filament beads also act as binding sites for molecular motors. These are modeled as harmonic springs *V*_{m} = *k*_{m}(*r* − *r*_{0})^{2}/2 with spring constant *k*_{m} and equilibrium length *r*_{0} = σ that walk on the filaments (*16*, *27*, *29*, *30*). Note that here, the discretization length *a*_{0} directly sets the motor step size. When the separation between two empty sites on nearby filaments is smaller than *r*_{att}, a free (unbound) motor can bind at a rate *p*_{att} = 1 in all simulations. When bound, the *n*_{s}) in discrete steps of size *a*_{0} at a rate *t* is the motor time step, and Δ*E* is the energy difference before and after the motor has stepped. Motors detach when one leg reaches the end of the filament, when the motor is stretched beyond a length *r*_{off} = 3.25 *r*_{0}, or when a motor encounters an occupied site on a filament. After detaching, the free motors form a bath with homogeneous concentration *L*_{box} is the box size. Different choices can be made for different events that lead to motor-filament dissociation. What matters for the filament dynamics is the total number of parallel and antiparallel motors on the filaments.

The Langevin equation of motion for each bead *i* on filament *q**t* using the integrator proposed in (*40*). To satisfy the fluctuation-dissipation theorem, the random forces are Gaussian distributed with mean zero and variance 2*k*_{B}*T*γδ(*t* − *t*^{′}). The integration time step δ*t* and friction γ were chosen such that δ*t*/*m* = 0.005 and γδ*t*/*m* ≥ 0.005. For these parameters, the (passive) center-of-mass motion of a single filament is diffusive for center-of-mass displacements larger than a fraction of the length of a filament, i.e., for center-of-mass diffusion, the dynamics is essentially overdamped. Note that the simulation time increases linearly with the friction constant. Simulations typically contain 900 filaments of 21 beads and up to 3200 motors.

The coupling of filaments by molecular motors leads to sliding and binding forces that depend on the relative orientation of the filaments (see Fig. 1, A and B), where antiparallel motors (coupling two antiparallel filaments) exert larger forces than parallel motors (coupling two parallel filaments). Moreover, the activity induced by dimeric motors (one motor arm is grafted) tetrameric motors (both motor arms move simultaneously) is larger (*30*). The motor model here is a hybrid of the two, and both motor arms can move with equal probability but only one at a time.

Rotational diffusion is measured through the filament orientation correlation function, where the orientation is defined as the eigenvector corresponding to the largest eigenvalue of the moment of inertia tensor. Defects are traced using the method outlined in (*25*).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/30/eaaw9975/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank R. G. Winkler for the helpful conversations.

**Funding:**The authors acknowledge that they received no funding in support of this research.

**Author contributions:**G.A.V., T.A., and G.G. conceived the research. G.A.V. wrote all the simulation codes, performed the simulations, and analyzed the data. All authors discussed the results and wrote the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).