Filamentous active matter: Band formation, bending, buckling, and defects

See allHide authors and affiliations

Science Advances  22 Jul 2020:
Vol. 6, no. 30, eaaw9975
DOI: 10.1126/sciadv.aaw9975


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.


Active systems are driven by nonthermal, energy-consuming processes that lead to very rich collective behavior (14). 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 (913).

Various theoretical modeling approaches (1422) 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, 1720). 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, 2730). 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 ns beads of diameter σ connected by stiff harmonic springs with rest length a0 = L/(ns − 1), where a0 = σ. The filament area fraction ϕ=nfnsπσ2/4Lbox2 is varied by changing the number of filaments, nf, or the box size Lbox. Molecular motors are modeled by harmonic springs with rest length r0 and spring constant km. They attach to neighboring filaments with rate Γatt. Motors walk in the direction of the filament polarity with step length a0. The step rate is proportional to the probability pm0 to move a motor arm, which sets the bare motor velocity v0=a0pm0/Δt, where Δ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 roff. In the following, we use rescaled quantities: filament aspect ratio L˜=L/σ, motor-to-filament ratio n˜m=nm/nf, persistence length ˜p=p/L, motor spring constant k˜m=kmr02/kBT (kBT is the thermal energy), filament friction γ˜=γ/kmδtt is the simulation time step), box size L˜box=Lbox/L, and time t˜=tD0/L2, with the single passive-filament translational diffusion coefficient D0. For details, see Materials and Methods.


Activity drives the formation of polar domains

The microscopic origin of filament motion is the dynamic cross-linking of nf polar filaments by nm 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 nm motors are classified as nmp “parallel motors” that connect parallel filaments in the interior of domains and nmap “antiparallel motors” that dynamically cross-link and slide filaments at the interfaces between oppositely oriented domains relative to each other.

Fig. 1 Self-organized structures and dynamics of motor-filament mixtures driven by sliding of antiparallel filaments.

(A and B) Sketches of two consecutive motor steps when filaments are oriented (A) parallel and (B) antiparallel. The initially relaxed gray motor steps toward the filament polar end, it relaxes, and it makes a second step on the other filament and relaxes again. This results in net filament motion only for antiparallel motors. (C and D) Simulation snapshots for an initially disordered nematic system. The thin black box indicates the central simulation cell. (C) Short-time band formation and (D) long-time disordered structures. Filament colors indicate their orientation, illustrated by the color axis. (E) Number of antiparallel motors as a function of the persistence length p and the bare motor velocity, measured by pm0. (F) Peclet number as a function of the active force F˜act in Eq. 1. The table shows the parameter combinations used in (E) and (F). In all cases, (ϕ,L˜) = (0.66, 20).

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 filamentF˜act=nmapk˜m(rmr01)/nf(1)with the average relative extension rm/r0 − 1 of the antiparallel motors. In steady state, the motor force is balanced by an effective medium friction with friction coefficient γeffD01; the parallel filament velocity 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 Fact is proportional to nmap, which itself depends on nm, the persistence length ℓp, and the motor stepping rate pm0 (see Fig. 1E); thus, the same activity can be achieved by different parameter combinations. The parallel velocity and, thus, the Peclet number Pe = vL/D0 increase linearly with the activity in the system vFact (see Fig. 1F). This demonstrates that Eq. 1 provides an appropriate identification of the active forces. For systems that form regular bands, the velocity profile within the bands (see fig. S3, A and B) shows that for fixed Fact, the velocity is independent of ℓ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 n˜m and intermediate times, the filaments follow essentially ballistic trajectories, where the MSDv2τ2 is a signature of this active persistent motion. At long times, the filaments lose their initial orientation due to active rotational motion, and the MSD becomes linear with time, characterized by a plateau in Fig. 2A, with an active translational diffusion coefficient Dact that increases with increasing n˜m. The filament orientational autocorrelation functions decay exponentially with an active decay time τR, which also depends on nm. The active rotation time τR is shown in Fig. 2B both as a function of n˜m and F˜act. For small n˜m, increasing the number of motors first leads to a faster rotational motion, but at n˜m1.5, the motion becomes slower again because the number of antiparallel motors per filament nmap saturates (see fig. S4); instead, the excess motors increase the cross-linking density inside the polar domains, which reduces the filament velocity. The active rotation time τR decreases with increasing F˜act (or equivalently nmap) with a power law, τRF˜act32. The active diffusion coefficient, parallel velocity, and rotational correlation times for many different parameter variations are related by Dactv2τR (see Fig. 2C), consistent with the theory of active Brownian motion (34). This implies DactF˜act.

Fig. 2 Single-filament motion and active diffusion for various propulsion forces.

(A) Filament MSD divided by the lag time for different motor concentrations n˜m. The lag time τ is normalized with the passive single-filament rotation time τR0, and the MSD with the filament length L2. The parameters are (ϕ,˜p,L˜) = (0.66, 3.4, 20). (B) Normalized active rotation times τR as function of n˜m and F˜act. Solid lines are a guide to the eye, parameters as in Fig. 1F with (ϕ,L˜) = (0.66, 20). (C) Normalized active diffusion constants obtained from MSD curves as a function of the filament measured v2τR, parameters as in Fig. 1F with (ϕ,L˜) = (0.66, 20).

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 FbE2, 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 Fbnmap, this provides the scaling λ2p/nmap. The investigation of the stability of the system with simulations at various combinations of n˜m and ˜p allows us to determine the phase diagram shown in Fig. 3E. Here, the stability limit of the polarity-sorted bands increases with increasing persistence length, which nicely agrees with the predicted linear dependence of the buckling force on p.

Fig. 3 Onset of domain formation: Euler-like buckling instability of polar bands.

(A to C) Snapshots of the time evolution of an initially nematic system with (ϕ,˜p,L˜)=(0.66,10,40), for t/τR0=0.1,0.2,0.3, respectively. (D) Wavelength of the instability of polar-sorted bands as a function for varying persistence length and (ϕ,L˜,n˜m)=(0.66,40,1). (E) Stability phase diagram for various active forces and persistence lengths with parameters (ϕ,L˜)=(0.66,20). The solid line separates different phases and indicates the stability limit of the banded phase.

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 +12 and 12 topological defects, as shown in Fig. 4, which are characteristic for active nematics. Therefore, the polar character of the filament interaction seems only of secondary importance. A defect pair is formed when a polar band buckles by extensional forces, such that the convex side forms a +12 defect, while the concave side forms a 12 defect (see Fig. 4A). These are not the standard defects of active nematics, because the three domains that meet at a 12 defect display polar order, which implies that there can be active forces where domains with antiparallel polar order meet. In particular, our polar filament model gives rise to two types of 12 defects with different orientations of the polar domains around the defect core. One has C1 symmetry, and the other has C3 symmetry (see Fig. 4, B and C). Note that for the C1 defect in Fig. 4B, there are additional (unmarked) boundaries between antiparallel domains in the lower left and right corners where the motion is originated.

Fig. 4 Densities and MSDs of topological +1/2 and −1/2 defects.

(A) Snapshot of a pair of 12 (green) and +12 (magenta) defects. (B) Polar and nematic order near a 12 defect with C1 symmetry. (C) Polar and nematic order near a 12 defect with C3 symmetry. In (A) to (C), filament colors indicate the polar angle of filament orientation. In (B) and (C), the yellow lines indicate boundaries between polar domains. (D) MSD of the filament center-of-mass (red), 12 defect (green), and +12 defect (magenta). (E) Average distance between defects as a function of the active force. (F) Estimation of the domain size with three different approaches. Solid lines show the 1/F˜act dependence and symbols simulations with datasets from the table in Fig. 1 with (ϕ,L˜)=(0.66,20).

The f motion of the +12 defects and subdiffusive motion of the 12 defects, as follows from the defects MSDs in Fig. 4D and illustrated in movie S4. The annihilation of two defects is shown in movie S5. However, we do not find any differences in the dynamic behavior of the two types of 12 defects. The defect density depends linearly on the activity, which is reflected by the inverse square root dependence of the distance ldef 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 Ω̂(r) and Ω̂(r) (see fig. S10) provides the correlation lengths 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 ldef. 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 n˜m, chosen such that nmap/nf is comparable so that the systems have similar active force densities. We observe a gradual change from small local bundles held together by motors at low densities to a dense disordered nematic dominated by excluded-volume interactions similarly as in experiments (12). Although the structures observed strongly depend on ϕ, the (parallel) velocities still show a unique linear scaling with the force density vPeF˜act (see Fig. 5E). However, the active rotation times shown in Fig. 5F strongly depend on the density, especially for small densities and low activities, where the rotation time in the active systems is similar to the passive case. Nevertheless, for large activities and large densities ϕ > 0.2, again, a unique scaling behavior τRF˜32 emerges (see also Fig. 2B). In this regime, also the active diffusion coefficient Dactv2τR and domain size vτR show a universal scaling (see Fig. 5G). The chosen scaling variables (F˜act, τR/τR0, and vL/D0) are independent of the filament length, as indicated by three simulations for filaments twice as long as all other simulations.

Fig. 5 Filament density–dependent structure formation and single-filament dynamics.

(A to D) Snapshots for systems at various densities, ϕ = (0.15, 0.29, 0.44, 0.66), where nm is varied to result in similar nmap/nf0.165. (E to G) Dimensionless observables as a function of the active force for different area fractions. Symbols correspond to simulation data, dashed lines are used to guide the eye, and full lines indicate limiting power laws at large densities: (E) Peclet number, (F) active rotation time, (G) active diffusion coefficient v2τR (open squares), and domain size vτR (solid up-triangles). Colored labels in (E) refer to densities in (E) to (G). Black symbols in (G) correspond to simulations with longer filaments with (ϕ,˜p,L˜)=(40,0.66,1.7).

Last, we study the effect of activity and concentration on the effective filament stiffness or effective persistence length p*. This effective persistence length is calculated by an exponential fit of the tangent correlation function (see fig. S11D) in the motor-filament mixture at filament concentration ϕ and motor concentration n˜m. For passive systems, excluded volume interactions imply that p* increases with volume fraction ϕ. We find that for large enough activity, p* decreases like 1/F˜act for all filament concentrations (see fig. S11E).

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 DR 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, DR=v/p*. If the rotational diffusion constant is replaced by τR1, one obtains p*=vτR1/F˜, which is exactly what is found in fig. S11E. We note that for our systems with ϕ > 0.22, this scaling is identical to the one of vτR with the active force (see fig. S11).


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 +12 defects has been found, as in active gel theory, demonstrating that nematic order dominates the development of steady-state structures. The dynamics of 12 defects is found to be subdiffusive, of +12 defects superdiffusive, in agreement with active nematic theory (15). The individual motion of +12 defects in the direction of their head shows that the forces are on average extensile. Filament flexibility is a relevant property of active gels (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


Suspensions of nf semiflexible filaments and nm 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 ns beads of mass m with position vectors rqi, where q = {1. …nf} and i = {1. …ns}. The ns beads define ns − 1 segments of length a0 = L/(ns − 1) with unit tangent vectors t̂qi=(rqi+1rqi)/rqi+1rqi. The beads are connected by harmonic bonds with equilibrium length a0 and spring constant kr02/kBT=250. A bending potential Vb=κ(1t̂qi·t̂qi+1)/a0 gives rise to a tangent vector correlation function ⟨t(0) · t(s)⟩ = exp [ − s/ℓp] = exp [ − (d − 1)kBTs/2κ] with dimensionality d and arc length s along the contour, such that the bare persistence length is ℓp/L = 2κ/LkBT (41). Nonbonded inter- and intramolecular excluded volume and attractive interactions are described by a Lennard-Jones potential VLJ = ε[(σ/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 rc = 3σ denotes the potential cutoff radius. In all our simulations, kBT/ε = 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 ε/kBT = 1.0 and 4% less than the WCA potential with ε/kBT = 0.1. The filament beads also act as binding sites for molecular motors. These are modeled as harmonic springs Vm = km(rr0)2/2 with spring constant km and equilibrium length r0 = σ that walk on the filaments (16, 27, 29, 30). Note that here, the discretization length a0 directly sets the motor step size. When the separation between two empty sites on nearby filaments is smaller than ratt, a free (unbound) motor can bind at a rate (patt/δt)(nmf/nm), where nmf is the number of free motors and patt = 1 in all simulations. When bound, the nmb motors walk in the direction of filament polarity (from segment 1 → ns) in discrete steps of size a0 at a rate Γ=pm0/Δtexp(βΔE). Here, pm0 is the probability for an attempted step of a motor attached to two filaments, Δ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 roff = 3.25 r0, or when a motor encounters an occupied site on a filament. After detaching, the free motors form a bath with homogeneous concentration nmf/Lbox2, where Lbox 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 qmd2rqidt2=Fqiγdrqidt+Rqi(2)is integrated with a time step δt using the integrator proposed in (40). To satisfy the fluctuation-dissipation theorem, the random forces are Gaussian distributed with mean zero and variance 2kBTγδ(tt). 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 material for this article is available at

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.


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.

Stay Connected to Science Advances

Navigate This Article