Multiscale modeling and simulation of microtubule/motor protein assemblies
Abstract
Microtubules and motor proteins self organize into biologically important assemblies including the mitotic spindle and the centrosomal microtubule array. Outside of cells, microtubulemotor mixtures can form novel active liquidcrystalline materials driven out of equilibrium by ATPconsuming motor proteins. Microscopic motor activity causes polaritydependent interactions between motor proteins and microtubules, but how these interactions yield such largerscale dynamical behavior such as complex flows and defect dynamics is not well understood. We develop a multiscale theory for microtubulemotor systems in which Brownian dynamics simulations of polar microtubules driven by motors are used to study microscopic organization and stresses created by motormediated microtubule interactions. We identify polaritysorting and crosslink tether relaxation as two polarspecific sources of active destabilizing stress. We then develop a continuum DoiOnsager model that captures polarity sorting and the hydrodynamic flows generated by these polarspecific active stresses. In simulations of active nematic flows on immersed surfaces, the active stresses drive turbulent flow dynamics and continuous generation and annihilation of disclination defects. The dynamics follow from two instabilities, and accounting for the immersed nature of the experiment yields unambiguous characteristic length and time scales. When turning off the hydrodynamics in the DoiOnsager model, we capture formation of polar lanes as observed in the Brownian dynamics simulation.
I Introduction
Active matter, the novel class of nonequilibrium materials made up of selfdriven constituents, presents scientific challenges to our understanding of material properties and has the potential to provide new technologies such as autonomously moving and selfhealing materials. Examples of active matter include flocks of birdscavagna10 (), swarms of swimming bacteriazhang10 () or selfpropelled colloidal particlesbricard13 (), and the cellular cytoskeleton and cytoskeletal extractsnedelec97 (); surrey01 (); schaller10 (); sanchez12 (). Despite their differences in composition and length scale, these diverse systems show common features absent in equilibrium systems, including collective motion, nonequilibrium ordering transitions, and anomalous fluctuations and mechanical propertiesramaswamy10 (); marchetti13 (). Understanding and predicting the properties of active matter require new theoretical approaches and models applicable to farfromequilibrium, internally driven systems.
Mixtures of cytoskeletal filaments and motors are an important class of active matter that can be reconstituted outside the cell to form novel materials. Filaments driven into selforganized patterns such as vortices and asters are reminiscent of structures observed in cellsnedelec97 (); surrey01 (); tanakatakiguchi04 (); janson07 (); bendix08 (); pinot09 (); hentrich10 (); thoresen11 (); gordon12 (); schaller10 (). In earlier experiments, filaments were driven into static selforganized patterns such as vortices and asters, reminiscent of structures observed in vivo. In recent experiments, active networks were formed of microtubules (MTs) and synthetic multimeric kinesin1 motor complexes, with the aid of a depletantsanchez12 (); henkin14 (); keber14 (). In the presence of ATP, motor complexes can bind pairs of MTs and walk along MTs towards their plusends. When suspended in bulk, depletion interactions drove the formation of extended, highly ordered MT bundles characterized by bundle extension and fracture, and correlated with spontaneous largescale fluid flowssanchez12 (); keber14 (). When MT bundles were adsorbed onto an oilwater interface, they formed a dense, nematically ordered surface state and exhibited an active nematic phase characterized by the spontaneous generation and annihilation of disclination defect pairssanchez12 ().
Theoretical studies nakazawa96 (); kruse00 (); lee01 (); kruse03 (); liverpool03 (); sanakararaman04 (); aranson05 (); Liverpool05 (); ziebert05 (); ahmadi06 (); giomi08 (); zemel09 () have investigated aspects of these activematter systems at different scales, from the dynamics and mechanical properties of filament bundles to macroscopic behavior and stability of active suspensions. Inspired by the experiments of Sanchez et al. sanchez12 (), both Giomi et al. giomi13 (); giomi14 () and Thampi et al. TGY2013 (); thampi14a (); thampi14b (); thampi14c () have studied liquid crystal hydrodynamic models with fluid flow driven by an apolar active stress simha02 (); Voituriez05 (). In these rather general models the precise origins of the active stress driving the system are unidentified. Giomi et al. developed a theory for the speed at which defects move apart in active nematics, assuming the presence of a defect pair as an initial condition. Thampi et al. found an activityindependent velocityvelocity correlation length, as found in the bulk flow measurements of Sanchez et al., and studied defect dynamics in 2D simulations. These models reproduce qualitative features of the experiments. However, MT/motorprotein interactions are intrinsically polar, and how these polaritydependent microscopic interactions manifest themselves at meso or macroscopic scales is still unknown. Thus one theoretical challenge is how to resolve microscopic interactions between constituents in order to predict macroscopic material properties. While particlebased simulations can represent microscopic interactions in detail, computational cost typically limits crossscale studies. Continuum models are more tractable for describing dynamics at large scales, but can be difficult to connect to the microscopic dynamics quantitatively.
Here we construct a multiscale model that identifies the sources of destabilizing active stresses, and study their consequences in a largescale model GBGBS2015 (). We first perform detailed, hybrid Brownian dynamicskinetic Monte Carlo (BDkMC) simulations which incorporate excludedvolume interactions between model MTs, thermal fluctuations, explicit motors with binding/unbinding kinetics that satisfy detailed balance, and a forcevelocity relation. Active extensile stress is generated from polarity sorting of antialigned MTs, and from crosslink relaxation of polaraligned MTs. It also provides coefficients for polarspecific active stresses for a kinetic theory that incorporates polarity sorting and longrange hydrodynamic interactions, using a similar approach as that used to describe bacterial suspensions saintillan08 (); saintillan08a (); Wolgemuth08 (); subramanian09 (); Wensink12 (); ESS2013 (), where hydrodynamic instabilities lead to largescale collective motions including jets and vortices pedley92 (); koch11 (); simha02 (); saintillan08 (); saintillan08a (); sokolov09 (); kurtuldu11 (). We use this model to study actively streaming nematic states on an immersed surface, as in the Sanchez et al. experimentssanchez12 (). Numerical experiments demonstrate dynamics strikingly similar to the experiments, with largescale turbulentlike fluid flows and the persistent production and annihilation of defects. We correlate the defect dynamics with specific flow structures and with active stresses. We identify the hydrodynamic instability of nearly 1D coherent “cracks” as being source of the persistent dynamics. When turning off the induced background surface flow in the kinetic model, we capture the formation of polar lanes observed in the BDkMC simulation.
Ii The Microscopic Model
Figure 1 outlines the basic physical picture that underlies both our BDkMC simulations and the continuum kinetic model. Consider an immersed suspension of polar MTs, each with a plusend oriented director , and all of the same length and diameter (Fig. 1a). Adjacent MTs are coupled by plusend directed crosslinking motors consisting of one motor head on each MT connected by a tether that responds as a spring to stretching (Fig. 1b). The motor on each crosslink endpoint moves with a linear forcevelocity relation visscher99 (): , where is the magnitude of the crosslinking force, is the maximum translocation velocity, and is the stall force. For a nematically aligned suspension there are two basic types of MT pair interaction. For polar antialigned MTs (Fig. 1c) the motors on each end of an active crosslink move in opposite directions, stretching the tether. This creates forces on each MT that, acting against fluid drag, slide the MTs relative to each other towards their minusends. This process is termed polarity sorting nakazawa96 (). Conversely, for polaraligned MTs the motors on each end of the crosslink move in the same direction, there is little or no net sliding, and the tether pulling on the leading motor causes stretched tethers to relax (Fig. 1d).
Iii Brownian dynamicskinetic Monte Carlo model and simulations
We first perform 2D Brownian dynamicskinetic Monte Carlo (BDkMC) simulations of MTs driven by explicit motors with binding/unbinding kinetics. The main purpose is to quantify local MT pair interactions, with longranged hydrodynamics neglected due to its high computational cost. Our model is related to previous simulations of filaments with crosslinking motors nedelec07 (); head11 (); head14 (), but new in our work are algorithmic improvements for handling crosslinks and neglect of filament elasticity that allow us to more accurately treat the statistical mechanics of crosslinking motors, simulate larger systems and measure the stress tensor.
The particlebased BDkMC simulations use a simple, tractable model of active biomolecular assemblies that capture key physical features, including excluded volume interactions between filaments, attractive and sliding forces exerted by motors, and the thermodynamics and kinetics of crosslinking motor binding and unbinding. Filaments (MTs) are represented as perfectly rigid rods (discorectangles in 2D) of length and diameter that undergo Brownian dynamics. Forces and torques on the filaments occur due to motormediated forces, particleparticle repulsion, friction, and thermal forces. To simulate the Brownian motion of filaments, we adopt the computational scheme of Tao et al. tao05 (), which has been used successfully in simulations of concentrated solutions of highaspectratio particles. In this scheme, the filament centers of mass equations of motion are
(1) 
for all filaments , where the random displacement is Gaussiandistributed and anisotropic, with variance
(2) 
In the above, is Boltzmann’s constant and is the absolute temperature. Here is the inverse friction tensor
(3) 
where and are the parallel and perpendicular drag coefficients of the rod, and is the systematic (deterministic) force on particle . The equations of motion for particle reorientation are
(4) 
where is the rotational drag coefficient, is the systematic torque on particle , and the random reorientation is Gaussiandistributed, with variance
(5) 
The WeeksChandlerAndersen (WCA) potential between rods is
(6) 
where , is the minimum distance between the two finite line segments that define the filament axis, and sets the energy scale of the potential. Note that is an implicit function of the centerofmass positions and orientations of the two interacting MTs. For this value of , the typical distance of closest approach between rods is comparable to , and the thermodynamic properties closely resemble those of hard rods with aspect ratio , a model that is wellcharacterized both in 2D bates00 () and 3D bolhuis97 (); mcgrother96 ().
Because the Brownian dynamics scheme involves random particle displacements and rotations, close contacts between rods that produce large forces and torques occasionally occur, leading to instabilities in the dynamics. Such instabilities are avoided by softening the WCA potential at short distances to keep the resulting forces and torques within reasonable bounds tao05 (). At the same time, we adjust the integration timestep to ensure that pairs of interacting particles probe the softened region of the potential infrequently, so that excluded volume effects are properly accounted for.
The frictional forces are orientation dependent: translational diffusion is characterized by two diffusion constants, and , which describe diffusion perpendicular and parallel to the rod axis, respectively, and is the rotational diffusion coefficient. For spherocylinders where , the diffusion coefficients are lowen94 ()
(7) 
(8) 
and
(9) 
Here is the fluid viscosity and is the temperature in energy units. Note that is approximately a factor of two larger than .
To model motormediated interactions and activity, we implement a semigrand canonical ensemble in which a reservoir of motors is maintained in diffusive contact at a fixed chemical potential with filaments to/from which they can bind/unbind. The motors are assumed to be noninteracting both in solution and in the bound state, so the motor reservoir can be treated as an ideal solution, and there is no steric interference among bound motors. Bound motors have a free energy , where is the extension of the motor tether, which depends implicitly on the relative positions and orientations of the two filaments to which the motors is attached and on the positions of the points of attachment of the motor on the filament axes. We treat motor attachment (detachment) as a onestep process in which motors bind to (unbind from) two filaments simultaneously, and we assume a binding rate of
(10) 
and an unbinding rate of
(11) 
where is the inverse temperature in energy units. This choice of binding and unbinding rates ensures that the correct equilibrium distribution is recovered for static (nontranslocating) crosslinks, is a convenient choice from a computational standpoint, and has been used previouslylenz14 (). Given a distribution of motors bound to filaments, we compute the forces and torques exerted on MTs by differentiating with respect to the filament positions and orientations. As discussed in section II, the endpoints of bound motors translocate toward the plus ends of the MTs to which they are attached with a forcedependent velocity. Motors unbind immediately upon reaching the plus end of either of the two filaments to which they are attached.
Because the motor unbinding rate is , independent of motor tether extension, the probability that a given motor unbinds in a time interval is , and the average number of motors that unbind in is , where is the current number of bound motors. The number of motors that unbind in a time interval follows a binomial distribution,
(12) 
In one timestep we remove randomly selected motors, where is determined by sampling from the binomial distribution.
The kinetic MC procedure for motor binding is involved, because the rate of motor binding depends on motor tether extension, which in turn depends on the relative positions and orientations of the two MTs to which the motor is attached and on the positions of the points of attachment of the motor along the filament axes. To compute the relative probability and rate of motor binding to specific binding sites on a given pair of filaments we consider the statistical mechanics of the filament/motor system in the equilibrium limit of nontranslocating crosslinks. The semigrand canonical partition function of the filament/motor system is
(13) 
where is the fugacity of the motor reservoir and is the number of bound motors. Here, is the canonical partition function of a system of filaments and bound motors,
(14) 
where labels the particle positions and orientations, is the filament potential energy, including interparticle interactions and external potentials, is the singlemotor partition function, and interactions between bound motors have been neglected. The singlemotor partition function depends on the filament positions and orientations , i.e., . Substituting Eq. (14) into the grand partition function and carrying out the summation over leads to
(15) 
where
(16) 
In the limit in which the rate of motor binding and unbinding is large compared with the filament diffusion rate (adiabatic limit), plays the role of an effective motormediated filament interaction potential that depends on the chemical potential of the reservoir. Staticcrosslinkmediated interactions are generally attractive and shortranged, and bear a strong resemblance to depletiontype potentials bolhuis97 ().
The singlemotor partition function can be written as a sum of pairwise partition functions,
(17) 
where the sum ranges over all distinct pairs of filaments, and the pairwise partition function is
(18) 
Here the integration variables and parametrize the positions of motor endpoints on filaments and , respectively, is length of a motor between points specified by and , and is the linear density of binding sites on a single filament. Then we can write the effective motor potential as the sum of pairwise effective interactions,
(19) 
where
(20) 
is the effective motormediated pair potential in the adiabatic limit. Insertion of motors with the correct relative statistical weight in a kinetic MC procedure requires evaluation of the pairwise partition function (Eq. (18)) for all pairs of filaments. If the motor energy increases rapidly (e.g., quadratically) with increasing motor extension, the partition function (and the corresponding adiabatic effective potential ) falls off rapidly with increasing minimum distance between filament axes, and is nonnegligible only for pairs of filaments in close proximity. Thus, the pairwise partition function is analogous to a shortrange interaction potential, and the usual techniques for efficient handling of shortrange interactions (e.g., neighbor lists) can be applied. To efficiently evaluate the double integral in Eq. (18), note that for motors modeled as zeroequilibriumlength harmonic springs, the integrand can be expressed as a sum of bivariate normal distributions. Then reduces to a sum of cumulative bivariate normal distributions, which can be rapidly evaluated using standard numerical procedures vesely06 ().
To proceed further, we consider the statistical mechanics of the motor subsystem for fixed filament positions and orientations. The grand partition function for the motor subsystem is given by
(21) 
and the equilibrium number of bound motors for a given filament configuration is
(22) 
where is the average number of motors between filaments and ,
(23) 
Note that , so the problem of computing is equivalent to that of computing . Introducing the explicit form of quadratic potential for harmonic motors leads to
(24) 
where , and where the implicit dependence of on filament coordinates has been suppressed.
The average number of motors that bind to filaments in a time interval is
(25) 
As above, the number of motors that bind in the interval follows a Poisson distribution,
(26) 
In the kinetic MC cycle, the number of bound motors is drawn from this distribution, and motors are inserted by first selecting pairs of filaments with relative probability and then sampling from the appropriate bivariate normal distribution to choose motor endpoints that lie on the selected pair of filaments.
The overall hybrid BDkMC procedure thus consists of the following steps:

Displace each motor endpoint by along the MT to which it is attached with translocation velocity given by the forcevelocity relation.

Determine the number of motors that unbind in the time interval by drawing from a binomial distribution (Eq. (12)), and remove this number of motors at random.

Compute average number of bound motors for all pairs of MTs (Eq. (23)) and determine the number of motors that bind in the time interval by drawing from a Poisson distribution (Eq. (26)). Randomly select pairs of MTs with relative probability , and insert a motor between each selected pair of MTs by sampling from a bivariate normal distribution.
The properties of the model depend on seven dimensionless parameters (tables 1 and 2): the MT aspect ratio , the MT packing fraction , the range of motor mediated interaction , the motor concentration , the motor run length , the motor stall force , and the Peclet number (the ratio of translocation and diffusion rates) . With current methods, it becomes more computationally expensive to simulate systems with MTs of high aspect ratio (e.g., ). The computation time scales approximately as . If doubles, the linear dimension of the the box in the longitudinal direction must be doubled to study the same number of rods. We use square boxes to avoid any loss of information upon nematic director reorientation. Therefore the number of rods scales as . Longer rods also have slower dynamics, because the translational and rotational mobilities go as to leading order. Therefore the timescale to reach steady state scales approximately linearly in . We present here results of simulations with for which we performed simulations of relatively large systems for long times over a wide range of parameters. A more limited investigation of longer rods reveal qualitatively similar behavior.
Parameter values  

Quantity  Parameter  Value  Notes 
Thermal energy  4.11 J  Room temperature  
MT length  250 nm  Chosen  
MT diameter  25 nm  Ref alberts07,  
Energy scale of the WCA potential  Refs bates00,; bolhuis97,; mcgrother96,  
Fluid viscosity  1.0 Pa s  Cytoplasmic viscosity, ref Wirtz09,  
Linear density of motor binding sites along MT  –  Appears only in dimensionless concentration  
Motor chemical potential  –  Appears only in dimensionless concentration  
Motor binding free energy  –  Appears only in dimensionless concentration  
Motor speed (zero force)  Reference 4.5 m/s, range 0.14–18 m/s  Of order 1 m/s, ref visscher99,  
Unbinding rate of motors  28.1 s  Processivity of 160 nm, ref schnitzer00,  
Stall force  1 pN  Ref visscher99,  
Motor spring constant  0.013 pN/nm  Decreased from ref Coppin95, to give appropriate range of motormediated interaction for zeroequilibriumlength springs 
Dimensionless parameter values  

Quantity  Parameter  Value  Notes 
MT packing fraction  0.54  Chosen to give nematic state at equilibrium in the absence of motors  
MT aspect ratio  10  
Motor concentration  1  Chosen to give average of 2 motors per nearby MT pair  
Range of motor interaction  Chosen to be of order 1 for a shortrange interaction  
Motor run length  Reference 0.64, range 0.2–12.8  Motorinduced active stresses are largest when is of order 1.  
Motor stall force  6  
Peclet number  Reference 0.68, range 0.02–2.7 
iii.1 Measurement
The dynamics and stresses experienced by individual MTs depends strongly on their local environment, in particular on the relative polarity of neighboring MTs. To identify subpopulations of MTs with distinct local environments, we define a local polar orientational order parameter
(27) 
where is the motor pair partition function defined above. Since falls off rapidly with increasing pair separation, only near neighbors of particle are included in the sums in Eq. (27). The polar order parameter ranges from (MT surrounded by neighbors of opposite polarity) to (MT surrounded by neighbors of the same polarity).
The osmotic stress tensor of a periodic system of interacting MTs at temperature in a dimensional volume is given by
(28) 
where the first and second terms on the righthand side represent the ideal gas and interaction contributions, respectively, is the unit tensor, and is the virial tensor,
(29) 
where the sum ranges over all interacting pairs of MTs. The angular brackets in Eqn. (28) denote an average over time. Here we have assumed that the temperature of the system is isotropic and welldefined, so that
(30) 
where is the momentum of MT and is the MT mass (here assumed the same for all MTs). Filaments have momentum based on their instantaneous movements on short timescales. This motion is in thermal equilibrium with the background fluid, connecting molecular motion to Brownian motion. While this relation is clearly true in the equilibrium case, it’s less obvious that this it holds for active MT/motor systems. However, a purely mechanical definition of osmotic pressure leads to the same expression even for nonequilibrium particle suspensions in the low Reynolds number hydrodynamic regime brady93 (), and we will assume that Eq. (28) holds in the following discussion.
The isotropic pressure is defined as
(31) 
The average extensile stress is
(32) 
where the direction corresponds to the average nematic director orientation. We further resolve the stress tensor into contributions from subpopulations of MTs, for example according to the local polar order parameter introduced above. This can be done by writing the total virial as the sum of contributions from individual MTs,
(33) 
where
(34) 
To calculate the pair extensile stresslet as a function of the local polar order , we calculate the virial per spherocylinder. At a given time point, each interaction gives an associated virial contribution for the pair. The singleMT virial contribution is taken to be half of the pair’s contribution. Contributions from forces for all interacting partners are summed up to give the virial contribution for each MT. Similarly, the local polar order parameter is calculated for each MT. Then the virial anisotropy contribution per MT in the nematic reference frame is determined based on its local polar order. After repeating for all time points, the histogram is normalized, leading to the calculation of the extensile pair stresslet per MT as a function of .
To calculate the extensile pair stresslet in bulk simulations, we consider interacting MTs only. At each time point, the total number of interactions is calculated by summing the number of pairs for which there is a nonzero force. The total parallel and antiparallel virials in the director reference frame are calculated. Any interactions between pairs with contribute to the polaraligned virial, and the remainder contribute to the antialigned virial. This measurements is time averaged and the extensile pair stresslet calculated by dividing the average virial anisotropy by the average number of interactions.
iii.2 Extensile stress and its origins
Figure 2 illustrates the longtime behavior of MT suspensions in the BDkMC simulation model (also see video S1). Fig. 2a shows a simulation of MTs interacting only through thermal fluctuations and steric interactions (without motors). The system develops a 2D nematic state consistent with previous workbates00 (). Figure 2b shows the result of adding immobile crosslinkers with full binding/unbinding kinetics. The system shows MT bundling due to shortrange crosslinkinduced attraction. Figure 2c shows the behavior with motors. The system now shows active MT flows driven by polarity sorting, leading to the formation of polar lanes (domains of MTs with similar polarity). These polar lanes are highly dynamic and show large fluctuations. The meansquared displacement of MT position as a function of time shows diffusive behavior at longtimes in the equilibrium cases (Figs. 2a and b) and for active MTs when measured perpendicular to the average alignment direction. For motion parallel to the average alignment direction, the active MT meansquared displacement is superdiffusive and nearly ballistic at long times (Fig. 3a).
We characterized the dynamical properties of bound motors for polaraligned and antialigned MT pairs. For two MTs labeled and with orientations and and centerofmass diplacement , we define the pair’s longitudinal displacement by . For antialigned MT pairs () undergoing motordriven relative sliding, is negative when the MT pair is contracting (minusends closer than plusends), and becomes positive when the MT pair is extending (plusends closer than minusends; see Fig. 1). When crosslinks are immobile or for motors on polaraligned MTs (), the distribution of motors as a function of is symmetric (Fig. 3a). However for motors on antialigned MTs, the distribution of motors skews toward positive values of : more motors are bound during the extensile motion of the pair. This asymmetry occurs because of the translocation of the motors toward the MT plusends. This biases MT pairs toward extension, yielding an extensile stress that drives active flows (see below).
The distribution of motor extension alters significantly when crosslinks translocate (Fig. 3c). The minimum value of is approximately 1 due to excludedvolume interactions between MTs. For polaraligned pairs, the distribution is shifted toward smaller extensions than in the equilibrium case due to nonequilibrium tether relaxation, with important implications for the generation of extensile stress, as discussed below. For antialigned pairs, the distribution is shifted toward positive extension due to oppositely directed motor motion; this combination of motor extension and motion applies active forces that drive polarity sorting.
We measured the displacement distributions and average velocities of MTs along the nematic director and found that both are strongly correlated with an MT’s initial local polar environment. Defining the nematic director , we calculated MT displacement distributions in time along the projection of the local filament orientation vector onto the nematic vector: . In order to examine dynamical behavior on timescales comparable to the diffusion timescale, we grouped the MT displacements at a lag time of (chosen to clearly illustrate the different distributions) and their initial polar environment . For MTs in an initially polar environment (, the displacement distribution is approximately Gaussian with mean near zero, consistent with diffusivelike dynamics (Fig. 3d). For MTs in an initially antipolar environment (), we again find an approximately Gaussian displacement distribution, but the mean is shifted toward the MT’s minus end (Fig. 3d). This profile is consistent with drift plus diffusion dynamics. For more mixed initial environments (), we find that the dynamics are more complicated and are not likely described by a simple drift and diffusion model (Fig. 3d). MTs in initially mixed or antipolar environments exhibit significant displacements toward their minus ends due to antipolar sliding.
To further examine the polaritydependent MT movements, we measured instantaneous MT velocity component along the nematic director, at . Velocities of MTs are not constant because MTs experience relatively rapid changes in the polarity of their neightbors. MTs in initially antipolar environments tend to slow down rapidly, indicating that they move into more mixed environments, while MTs in polar or mixed environments tend to maintain their velocities for longer times. Filaments in polar environments have velocities near zero. (Fig. 3e). The instantaneous velocity depends approximately linearly on the local polar environment, as expected when filament movements are determined mainly by polarity sorting (Fig. 3f).
We measured the timeaveraged bulk stress tensor for our active particle system, and find that, over a wide range of parameters, is anisotropic with larger components in the average MT alignment direction than in the perpendicular direction. That is, since the MT alignment direction is essentially , the stress difference is positive, which corresponds to an extensile stress. Static crosslinkers or no motors (Fig. 2a, b) yield an isotropic . The stress difference can be expressed as the sum of pair interactions between nearby MTs, with each pair contributing a stresslet (with units of forcelength), prior to division by the bulk volume. We have characterized how the stresslet varies with system parameters and configurations. The average pair stresslet increases with the motor speed up to a maximum where the typical motor run length is the MT length (Fig. 3g). Increasing further leads to decreasing because the motors rapidly move to the ends of the MTs and unbind. To understand the origins of extensile stress, we studied how varies with the local polar environment (27). The stresslet is largest when is near , suggesting that polarity sorting is the dominant source of pairwise extensile stress (Fig. 3h). As increases, drops with approximate linearity, at least away from the two isolated peaks that close examination show originate through strong steric interactions of nearly parallel MTs. Nearly, but not exactly, parallel MTs experience aligning torques due to motormediated attraction; the resulting steric collisions tend to promote pair extension that increases the extensile stress for nearlyaligned pairs (relative to perfectly aligned pairs).
The extensile stress from antialigned pair interactions arises from asymmetries during polarity sorting: if an MT pair begins sliding when the two minusends touch () and slides under a force proportional to pair overlap until the two plusends meet (), then the contractile motion would perfectly balance the extensile motion and the total extensile stress would be zerokruse00 (); Liverpool05 (); lenz12 (); lenz14 (). In our simulations we observe two effects that break this symmetry. First, MTs are unlikely to begin interacting exactly when their minus ends are touching, decreasing the range of negative over which sliding occurs. Second, more motors are bound on average during extensional motion (so that ; see Fig. 3b).
We also find the surprising and counterintuitive result that remains positive even when is near 1, that is, for polaraligned pairs of MTs. This effect occurs due to an interplay between motor motion and excludedvolume interactions. We propose that the effect can be understood by considering equilibrium and nonequilibrium motor relaxation. For immobile motors, the system is at equilibrium and the stress tensor is isotropic; attractive interactions due to motors are balanced by excluded volume interactions and thermal fluctuations, and the system is at mechanical equilibrium. When motors are active, stress anisotropy becomes possible due to the nonequilibrium nature of the motor forcevelocity relation. The tether of a longitudinally stretched motor pulls back on the leading motor, slowing it, and pulls forward on the trailing motor. Hence, the motor relaxes its longitudinal extension. This effect is observable in Fig. 3c as a slight but significant shift in the distribution of motor extension toward smaller values relative to the equilibrium case. As a result, the motorinduced contractile stress along the MT alignment direction is decreased, while there is no change in the transverse stress induced by motors. This leads to a net anisotropic extensile stress in the alignment direction. In this scenario, we would predict that if the motors had a forceindependent velocity, the polaraligned extensile stress would vanish because the longitudinal motor extension would be unable to relax. We tested this prediction by studying how varies as stall force increases for simulations of perfectly aligned (unable to rotate) isolated filament pairs. We find that the extensile stress changes little with stall force for antialigned MT pairs. However, for polaraligned MT pairs the extensile stress drops as stall force increases and goes to zero for infinite stall force (which corresponds to forceindependent velocity, Fig. 3i). When effects of filament rotation are also included, the results are more subtle; we find that the interplay of filament rotation and motor activity can induce extensile stress for polaraligned pairs in bulk simulations even for infinite stall force.
While the extensile stress due to polaraligned MT pairs is typically a factor of 2–5 smaller than for antialigned pairs, when measured per pair (Fig. 3h), polarity sorting and the tendency to form polar lanes (Fig. 2c) lead to larger numbers of polaraligned MT pairs than of antialigned. In our BDkMC simulations, which lack the effect of hydrodynamics, the overall contributions of polaraligned and antialigned pairs to the extensile stress are comparable.
Iv Continuum kinetic theory
The BDkMC simulations show how polarspecific MTpair interactions give rise to extensile active stresses. To study the effect of hydrodynamic interactions and to make analytical predictions we have developed a DoiOnsager theory doi88 () similar to those used to describe the dynamics of motile rod suspensions saintillan08 (); saintillan08a (); ESS2013 (). The theory’s fluxes and active stresses arise from polaraligned and antialigned MT pair interactions produced by active motors. These stresses induce chaotic flows driven by the formation of disclination defects.
iv.1 Dynamics of polarity sorting
To coarsegrain the BDkMC simulation results and make connections with the kinetic model, we first derive a continuummechanics model to describe the MT dynamics. Here we assume the motor run length to be approximately the MT length, meaning that once bound, the motors will stay on the MTs until reaching the plus ends. As shown in Fig. 4, we consider a nematically ordered local cluster of MTs undergoing polarity sorting, with MTs pointing rightwards and MTs pointing leftwards. Let all the MTs in this cluster be coupled by active motors which create springlike forces between the MTs, and whose bound ends move at a characteristic (constant) speed toward MT plusends. For an antipolar MT pair this induces a relative sliding, each towards its negative end. The cluster is assumed small enough so that all MTs experience the same local flow field. Using Stokesian slender body theory keller76 () we can find the velocities of the left and rightward pointing MTs. For each MT, the center locates at , with the director . We assume that in the cluster there are MTs pointing leftwards (, with superscript ), and MTs pointing rightwards (, with superscript ). Each antialigned pair (say the th and the th MT) shares () motors
(35) 
where and . As shown on the right in Fig. 4, one motor locates at , and the other locates at , with initial positions and . The characteristic motor speed is constant for the antialigned pair. Hence the distance between the two motors in the tangential direction can be calculated as
(36) 
where , . When the motor is walking, it behaves like a linear spring with rigidity by exerting equal and opposite forces
(37) 
As a result, the two MTs slide past one another undergoing polarity sorting. Following slenderbody theory keller76 (), the MT speed is given by , leading to
(38) 
where , and is the fluid viscosity. We seek the timedependent solutions of the form . The coefficients and can be solved as
(39) 
leading to
(40) 
which suggests . This expression shows that the speed of each population depends on how many opposing MTs there are to pull against, with their drag as the anchor, and their relative velocity fixed at by the motor protein speed. This latter observation is in agreement with observations of antialigned sliding of MTs in the mitotic spindle yang08 ().
Next, we consider a general situation when the MTs are not perfectly aligned but with an intersection angle, i.e., where is a small angle between the and the MTs. As discussed later, at high concentration, the steric interactions align the neighbouring MTs, which makes the smallangle assumption a reasonable approximation. Similar to the perfectlyaligned case, the positions of the two motors can now be written as:
(41) 
where , and . So the relative distance becomes
(42) 
where , . The motors exert tangential force . Following the same procedure, we seek solutions of the form , yielding
(43) 
Then the relative moving speed of the two MTs becomes
(44) 
When , Eqs. (43) and (44) exactly recover the solutions in (40) for the perfectlyaligned case. To further coarse grain the above results to facilitate a continuum modeling as discussed below, we take an average in of Eq. (44) which directly yields a translational particle flux .
iv.2 Flux velocity, active stress and kinetic model
The system is described by a distribution function of MT centerofmass positions and polar orientation vectors (), evolved through a Smoluchowski equation
(45) 
which reflects conservation of particle number. Here and are MT conformational fluxes. Important macroscopic quantities for describing a polar nematic system are the local concentration , the local polarity vector , the secondmoment tensor which arises generically in capturing active stresses produced by active suspensions simha02 (), the (tracefree) order parameter tensor , with or 3 the spatial dimension, and the fourth moment .
Slenderbody theory yields the forces each rod exerts on the fluid, and hence the volumeaveraged stress batchelor70 () by polarity sorting can be calculated. If the cluster occupies a volume , the induced extra stress tensor from antialigned sorting is . Here is proportional to fluid viscosity , and with the signed distance between the centerofmasses of the and oriented subclusters. If the cluster is extending then , as would be the case if motor protein binding and unbinding kinetics biased motor densities towards the plusend of the MTs. This is seen in the BDkMC simulations (Fig. 2e), and is associated with local extensile flows similar to those of motile Pusher particles which collectively can drive macroscopic flow instabilitiessaintillan08 (); saintillan08a (); saintillan12 (). The antialigned pair stresslet strength can be derived as . When taking as , we extract the value of from the BDkMC simulations.
While active motoring between polaraligned MTs yields little MT mobility, the BDkMC simulations show that it does yield an extensile stress. However, unlike polarity sorting we lack a simple firstprinciples model of how polar interactions yield extensile stress, though the number of polar pair interactions within a cluster scales as . Given that the anti and polaraligned stresses are of the same order (Fig. 2h) we assume the form . Comparison with the BDkMC simulations suggests that .
We have generalized this simple example to a continuum model that captures polarity sorting of MTs and the dependence of the stress upon the local polarity of the MT field. The fluxes for Eq. (45) are given in dimensionless form by
(46)  
(47) 
To nondimensionalize the above equations, we assume that there are MTs in the entire computational domain of volume . At high concentration, it is useful to introduce an effective volume fraction where is the mean number density doi88 (); ESS2013 (). Further, we choose the characteristic length scale , the velocity scale , as well as the stress scale . In Eq. (46), is the background fluid flow, and the last term yields translational diffusion with constant . For nematically ordered suspensions, the term exactly reproduces the cluster velocities induced by polarity sorting given above (note that for a perfectly polar system, no polarity sorting occurs and the flux makes no contribution.) In Eq. (47), the MTs are rotated by the background flow gradient according to Jeffery’s equation Jeffery22 () while the second term arises from the MaierSaupe potential with coefficient which models torques and stresses arising from steric interactions at high concentration maier58 (); ESS2013 (). The last term yields rotational diffusion of the rod with constant . We do not account for MT rotation through interactions with the local field, as is appropriate when the MT field is nematically ordered. All constants have been made nondimensional using characteristic velocity , and a characteristic length appropriate for dense suspensions doi88 (); ESS2013 ().
Our system is closed by specifying how and are recovered from , which involves specifying the extra stress created in the fluid by activity and other sources. We assume the active stress arises separately from antialigned and polaraligned MT interactions, and construct it from and (i.e., the simplest symmetric tensors quadratic in ). In dimensionless form, the active stress tensor takes the form
(48) 
The first term (second term) captures active stress production via polarity sorting (motor relaxation) and exactly reproduces the form of () for nematically ordered suspensions. The total extra stress tensor is given by , where models extra stresses arising from flowinduced constraint forces on MTs and steric interactions ESS2013 ():
(49) 
where .
For bulk flow modeling one typically closes the system by balancing viscous and extra stresses and solving the forced Stokes equation and with velocity and pressure . This generates the background velocity and its gradient needed to evolve Eq. (45)saintillan08 (). However, this approach does not describe the streaming nematic experiments of Sanchez et al. sanchez12 (), where the active material is confined to an interface between oil and water, so surface motions are coupled to external fluid motions. To capture that coupling, we consider a flat layer of interacting MTs bound in the plane at , and immersed between two halfspaces filled with Newtonian viscous fluid (for simplicity, of the same viscosity). The activity in the MT layer generates a stress jump across the plane, and so generates a global 3D flow which is continuous at . In order to close the system, we solve the surface velocity in terms of the extra stress . To accomplish this, we first solve the (3D) velocity field of fluid flow using the Stokes equations
(50) 
where is a regular 3D spatial gradient operator. Under Fourier transform in , the above equations can be written as:
(51) 
where is a 2D wavevector, and is a 2D velocity field. When solving these equations in the upper (+) and lower () halves of the domain, we match at the MT layer through the continuous (2D) surface velocity , i.e., and . After some algebra, we obtain
(52) 
where and is a 2D unit wavevector. We further assume that the capillarity of the surface bounding the MT layer acts against concentration of MTs. We denote the liquid viscous stress as , and match the two solutions through a traction jump on the layer . Here is a 2D operator on the surface, and , arising from both the active stress due to MT activity and the stress due to a transverse pressure gradient within the MT layer which results in the background flow being incompressible in the plane (i.e., ). Then it is easy to eliminate and solve the surface flow in terms of as:
(53) 
It is useful to compare this expression to that for the 2D Stokes equation forced by a bulk stress: . The missing factor of in Eq. (53) profoundly changes the nature of system stability for the surface and 2D bulk systems. Equation (53) not only closes the system but facilitates a pseudospectral method to solve the Smoluchowski equation (45) and the fluid flow in a coupled manner.
iv.3 Flow, polarity and defects
Assuming 2D periodic boundary conditions, we have simulated our active polar nematic model over long times, using Eqs. (4549) as well as the velocitystress kernel (53). For the simulations shown here, we choose from , and fix (i.e., aspect ratio 10), , , and (estimated from the BDkMC parameters). The computation is performed on a 2D periodic domain of a square box with dimension . The governing equations are solved spectrally in a coupled manner, using the fast Fourier transform algorithm by expanding the variables in Fourier series and truncating the series after 200 – 400 modes in each spatial direction saintillan07 (); saintillan08 (); ESS2013 ().
Simulating in regions of flow instability we find persistently unsteady flows correlated with continual genesis, propagation, and annihilation of defect pairs. When we examine simulation results at late times, from initial data near uniform isotropy, we find dynamics that are complex and appear turbulent, qualitatively similar to those reported by Sanchez et al. sanchez12 () (Fig. 5). The surface velocity and vorticity show formation of jets and swirls (Fig. 5a, also see video S2). The local MT orientation is highly correlated with the flow structures, and the surface is littered with defects which propagate freely about the system (Fig. 5b, also see video S3). These defects exist in regions of small nematic order (dark blue), and are born as opposing pairs in elongated “incipient crack” regions. These are associated with surface jets, locally decreasing nematic order, and increasing curvature of director field lines. Characteristically, the defects propagate away along their central axis and have a much higher velocity than those of order. The relatively higher surface velocity in the neighborhood of a defect appears as a welllocalized jet, in the direction of defect motion, between two oppositely signed vortices.
The active force vector field is correlated with regions of rapidly changing nematic order (Fig. 5c, also see video S4). Large active force is present along an interconnected network of ridges correlated with the stringy regions of diminished nematic order and particularly with incipient cracks. Along such cracks, the active force points in the direction from which newly nucleated defects will emerge and propagate. Isolated highforce peaks correlate and move with defects, with the force pointing in the direction of their motion. Negative order defects are associated with regions of relatively low force magnitude, likely due to the local symmetry of the nematic director field.
We observe both nucleation and annihilation of defect pairs (Fig. 6). The birth and separation of a defect pair begins from an incipient crack wherein the initially smooth director field (e.g., lower arrow in Fig. 5b) morphs into singular forms in regions of low nematic order (Fig. 6a). Typically the positivelysigned defect moves away faster and roughly along its symmetry axis. Following annihilation of an oppositelycharged defect pair (Fig. 6b), the nematic order increases as the director field reknits itself into a smooth form (e.g., upper arrow in Fig. 5b). We examined how the polarity field changes near a defect and incipient crack (Fig. 6c). As the defect propagates, it leaves behind a region of increased polarity. The polarity field rapidly rotates across the incipient crack (by approximately ), and sometimes forms a shocklike structure that precedes the birth of a new defect pair. We measured the relative speed of the defect pairs (Fig. 6d). The speeds are similar to each other and on the order of the motor protein speed in our model (normalized to unity). This is consistent with experimental observations (cf. Fig. 3 of Sanchez et al. sanchez12 ()). The average fluid velocity around the defect pair is much lower than the defect speeds. Hence, as is the case for defects in more standard liquid crystalline materials, the defects here are not material structures carried along by the background surface flow.
Because our model is based on polarspecific fluxes and active stresses, the polarity field Tjhung11 (), polaritydependent active stresses and , and the local MT dynamics are coupled (Fig. 7). The polarity field develops considerable spatial variation with regions of high and low polar order (Fig. 7a, also see video S5). The two active stresses vary in strength depending on the local polarity – the polaraligned (antialigned) stress is large in regions of high (low) polar order (Fig. 7b,c). The antialigned stress yields the largest forces, by about a factor of 3 (close to the ratio ). The polarity field varies rapidly around defects, leading to gradients in the active stresses and large active force (open circles in Fig. 7). For comparison, we did another numerical test where we assume the active stress generated during motor tether relaxation is contractile (i.e., , Fig. 7df). The ratio between the antialigned and polaraligned stress is still close to the ratio . However, since the sign of the polaraligned stress changes, the two stresses exist in approximately the same regions.
To illustrate the dramatic variation of local MT fluxes with the local polarity field, we simulated the results of a photobleaching experiment in which a circular region is exposed to highintensity laser light to bleach the fluorescing molecules on the corresponding MTs AKSEW1976 () (Figs. 7gi). In a small highpolarity region (marked A in Fig. 7a), little or no polarity sorting occurs. Therefore the photobleached spot remains approximately circular (Fig. 7g top, h) and would deform due to the fluid flow over longer times. In a lowpolarity region of high nematic order (marked B in Fig. 7b), strong polarity sorting of antialigned MTs causes a photobleached spot to separate into two lobes (Fig. 7g lower, i). Each lobe mixes with unbleached surrounding MTs due to their active relative flux, showing decreased bleaching. Through the lens of our theory, this type of experiment probes the local polarity field, and hence the origins of active stress.
iv.4 Coherent structures and hydrodynamic instabilities
In our simulations, defect pairs are generated along elongated cracks that develop in regions of high polar order. To understand this instability, we consider nematicallyaligned MTs using reduced equations where particle diffusion is neglected (i.e., ) in (47). We then adopt bipolar solutions of the form , where the concentrations and and orientations are governed by: