Multiscale modeling and simulation of microtubule/motor protein assemblies

Multiscale modeling and simulation of microtubule/motor protein assemblies

Tong Gao, Robert Blackwell, Matthew A. Glaser, M. D. Betterton, Michael J. Shelley Courant Institute of Mathematical Sciences, New York University, New York, NY 10012
Department of Physics, University of Colorado, Boulder, CO 80309
July 6, 2019

Microtubules and motor proteins self organize into biologically important assemblies including the mitotic spindle and the centrosomal microtubule array. Outside of cells, microtubule-motor mixtures can form novel active liquid-crystalline materials driven out of equilibrium by ATP-consuming motor proteins. Microscopic motor activity causes polarity-dependent interactions between motor proteins and microtubules, but how these interactions yield such larger-scale dynamical behavior such as complex flows and defect dynamics is not well understood. We develop a multiscale theory for microtubule-motor systems in which Brownian dynamics simulations of polar microtubules driven by motors are used to study microscopic organization and stresses created by motor-mediated microtubule interactions. We identify polarity-sorting and crosslink tether relaxation as two polar-specific sources of active destabilizing stress. We then develop a continuum Doi-Onsager model that captures polarity sorting and the hydrodynamic flows generated by these polar-specific 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 Doi-Onsager 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 self-driven constituents, presents scientific challenges to our understanding of material properties and has the potential to provide new technologies such as autonomously moving and self-healing materials. Examples of active matter include flocks of birdscavagna10 (), swarms of swimming bacteriazhang10 () or self-propelled 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 far-from-equilibrium, 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 self-organized patterns such as vortices and asters are reminiscent of structures observed in cellsnedelec97 (); surrey01 (); tanaka-takiguchi04 (); janson07 (); bendix08 (); pinot09 (); hentrich10 (); thoresen11 (); gordon12 (); schaller10 (). In earlier experiments, filaments were driven into static self-organized 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 kinesin-1 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 plus-ends. 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 large-scale fluid flowssanchez12 (); keber14 (). When MT bundles were adsorbed onto an oil-water 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 active-matter 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 activity-independent velocity-velocity 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/motor-protein interactions are intrinsically polar, and how these polarity-dependent 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 particle-based simulations can represent microscopic interactions in detail, computational cost typically limits cross-scale 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 multi-scale model that identifies the sources of destabilizing active stresses, and study their consequences in a large-scale model GBGBS2015 (). We first perform detailed, hybrid Brownian dynamics-kinetic Monte Carlo (BD-kMC) simulations which incorporate excluded-volume interactions between model MTs, thermal fluctuations, explicit motors with binding/unbinding kinetics that satisfy detailed balance, and a force-velocity relation. Active extensile stress is generated from polarity sorting of anti-aligned MTs, and from crosslink relaxation of polar-aligned MTs. It also provides coefficients for polar-specific active stresses for a kinetic theory that incorporates polarity sorting and long-range hydrodynamic interactions, using a similar approach as that used to describe bacterial suspensions saintillan08 (); saintillan08a (); Wolgemuth08 (); subramanian09 (); Wensink12 (); ESS2013 (), where hydrodynamic instabilities lead to large-scale 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 large-scale turbulent-like 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 BD-kMC simulation.

Ii The Microscopic Model

Figure 1: (a) Schematic of a cluster of polar-aligned and anti-aligned MTs, with plus ends marked by red rings. Motors walk on neighboring MTs, and (b) exert spring-like forces with a piecewise linear force-velocity relation. (c) An anti-aligned MT pair. (d) A polar-aligned MT pair. Grey arrows characterize the magnitude of the extensile stress.

Figure 1 outlines the basic physical picture that underlies both our BD-kMC simulations and the continuum kinetic model. Consider an immersed suspension of polar MTs, each with a plus-end oriented director , and all of the same length and diameter (Fig. 1a). Adjacent MTs are coupled by plus-end 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 force-velocity 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 anti-aligned 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 minus-ends. This process is termed polarity sorting nakazawa96 (). Conversely, for polar-aligned 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 dynamics-kinetic Monte Carlo model and simulations

We first perform 2D Brownian dynamics-kinetic Monte Carlo (BD-kMC) simulations of MTs driven by explicit motors with binding/unbinding kinetics. The main purpose is to quantify local MT pair interactions, with long-ranged 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 particle-based BD-kMC 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 motor-mediated forces, particle-particle 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 high-aspect-ratio particles. In this scheme, the filament centers of mass equations of motion are


for all filaments , where the random displacement is Gaussian-distributed and anisotropic, with variance


In the above, is Boltzmann’s constant and is the absolute temperature. Here is the inverse friction tensor


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


where is the rotational drag coefficient, is the systematic torque on particle , and the random reorientation is Gaussian-distributed, with variance


The Weeks-Chandler-Andersen (WCA) potential between rods is


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 center-of-mass 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 well-characterized 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 ()




Here is the fluid viscosity and is the temperature in energy units. Note that is approximately a factor of two larger than .

To model motor-mediated interactions and activity, we implement a semi-grand 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 one-step process in which motors bind to (unbind from) two filaments simultaneously, and we assume a binding rate of


and an unbinding rate of


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 (non-translocating) 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 force-dependent 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,


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 non-translocating crosslinks. The semi-grand canonical partition function of the filament/motor system is


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,


where labels the particle positions and orientations, is the filament potential energy, including interparticle interactions and external potentials, is the single-motor partition function, and interactions between bound motors have been neglected. The single-motor 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




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 motor-mediated filament interaction potential that depends on the chemical potential of the reservoir. Static-crosslink-mediated interactions are generally attractive and short-ranged, and bear a strong resemblance to depletion-type potentials bolhuis97 ().

The single-motor partition function can be written as a sum of pairwise partition functions,


where the sum ranges over all distinct pairs of filaments, and the pairwise partition function is


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,




is the effective motor-mediated 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 non-negligible only for pairs of filaments in close proximity. Thus, the pairwise partition function is analogous to a short-range interaction potential, and the usual techniques for efficient handling of short-range interactions (e.g., neighbor lists) can be applied. To efficiently evaluate the double integral in Eq. (18), note that for motors modeled as zero-equilibrium-length 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


and the equilibrium number of bound motors for a given filament configuration is


where is the average number of motors between filaments and ,


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


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


As above, the number of motors that bind in the interval follows a Poisson distribution,


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 BD-kMC procedure thus consists of the following steps:

  1. Compute forces and torques on MTs, and evolve MT positions and orientations forward in time according to the Brownian dynamics equations of motion (Eqs. (1) and (4)).

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

  3. 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.

  4. 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 motor-mediated interaction for zero-equilibrium-length springs
Table 1: Parameter values of the BD-kMC simulation.
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 short-range interaction
Motor run length Reference 0.64, range 0.2–12.8 Motor-induced active stresses are largest when is of order 1.
Motor stall force 6
Peclet number Reference 0.68, range 0.02–2.7
Table 2: Dimensionless groups of the BD-kMC simulation.

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 sub-populations of MTs with distinct local environments, we define a local polar orientational order parameter


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


where the first and second terms on the right-hand side represent the ideal gas and interaction contributions, respectively, is the unit tensor, and is the virial tensor,


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 well-defined, so that


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 time-scales. 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


The average extensile stress is


where the direction corresponds to the average nematic director orientation. We further resolve the stress tensor into contributions from sub-populations 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,




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 single-MT 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 polar-aligned virial, and the remainder contribute to the anti-aligned 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: Snapshots of the BD-kMC particle simulations. Insets are zoomed views with motors explicitly shown in white. (a) System with no motors, illustrating the 2D nematic state. (b) An equilibrium system with static crosslinkers exhibits MT bundling due to short-range crosslink-induced attraction. (c) An active system with motors exhibits active flows and formation of polar lanes.

Figure 2 illustrates the long-time behavior of MT suspensions in the BD-kMC 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 short-range crosslink-induced 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 mean-squared displacement of MT position as a function of time shows diffusive behavior at long-times 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 mean-squared displacement is superdiffusive and nearly ballistic at long times (Fig. 3a).

We characterized the dynamical properties of bound motors for polar-aligned and anti-aligned MT pairs. For two MTs labeled and with orientations and and center-of-mass diplacement , we define the pair’s longitudinal displacement by . For anti-aligned MT pairs () undergoing motor-driven relative sliding, is negative when the MT pair is contracting (minus-ends closer than plus-ends), and becomes positive when the MT pair is extending (plus-ends closer than minus-ends; see Fig. 1). When crosslinks are immobile or for motors on polar-aligned MTs (), the distribution of motors as a function of is symmetric (Fig. 3a). However for motors on anti-aligned 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 plus-ends. This biases MT pairs toward extension, yielding an extensile stress that drives active flows (see below).

Figure 3: Measurements of BD-kMC simulations. (a) Mean-squared deviation of MTs as a function of time. (b) Mean velocity of filaments along nematic director as a function of time for different initial polar environments . (c) Histogram of motor extension , broken into contributions from polar-aligned and anti-aligned pairs in the active case. (d) Histogram of motor occupancy as a function of the particle filament longitudinal displacement , broken into contributions from polar-aligned and anti-aligned pairs in the active case. (e) Histogram of filament displacements at time separation for various initial polar environments . (f) Variation of average instantaneous velocity of filaments along the nematic director in time for different initial polar environments . (g) Variation of extensile pair stresslet with motor run length , showing results from the entire bulk simulation and contributions of polar-aligned and polar-antialigned pairs. (h) Variation of extensile pair stresslet with local polar environment . (i) Variation of extensile pair stresslet with motor stall force from simulations of isolated, perfectly parallel filament pairs.

The distribution of motor extension alters significantly when crosslinks translocate (Fig. 3c). The minimum value of is approximately 1 due to excluded-volume interactions between MTs. For polar-aligned 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 anti-aligned 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 diffusive-like dynamics (Fig. 3d). For MTs in an initially anti-polar 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 anti-polar environments exhibit significant displacements toward their minus ends due to anti-polar sliding.

To further examine the polarity-dependent 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 anti-polar 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 time-averaged 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 motor-mediated attraction; the resulting steric collisions tend to promote pair extension that increases the extensile stress for nearly-aligned pairs (relative to perfectly aligned pairs).

The extensile stress from anti-aligned pair interactions arises from asymmetries during polarity sorting: if an MT pair begins sliding when the two minus-ends touch () and slides under a force proportional to pair overlap until the two plus-ends 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 polar-aligned pairs of MTs. This effect occurs due to an interplay between motor motion and excluded-volume 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 force-velocity 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 motor-induced 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 force-independent velocity, the polar-aligned 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 anti-aligned MT pairs. However, for polar-aligned MT pairs the extensile stress drops as stall force increases and goes to zero for infinite stall force (which corresponds to force-independent 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 polar-aligned pairs in bulk simulations even for infinite stall force.

While the extensile stress due to polar-aligned MT pairs is typically a factor of 2–5 smaller than for anti-aligned pairs, when measured per pair (Fig. 3h), polarity sorting and the tendency to form polar lanes (Fig. 2c) lead to larger numbers of polar-aligned MT pairs than of anti-aligned. In our BD-kMC simulations, which lack the effect of hydrodynamics, the overall contributions of polar-aligned and anti-aligned pairs to the extensile stress are comparable.

Iv Continuum kinetic theory

The BD-kMC simulations show how polar-specific MT-pair interactions give rise to extensile active stresses. To study the effect of hydrodynamic interactions and to make analytical predictions we have developed a Doi-Onsager 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 polar-aligned and anti-aligned 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

Figure 4: Schematic for a cluster of MTs undergoing polarity sorting. The plus-ends are marked by red rings. On the right: an anti-aligned pair of the th and the th MTs.

To coarsegrain the BD-kMC simulation results and make connections with the kinetic model, we first derive a continuum-mechanics 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 spring-like forces between the MTs, and whose bound ends move at a characteristic (constant) speed toward MT plus-ends. For an anti-polar 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 anti-aligned pair (say the th and the th MT) shares () motors


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 anti-aligned pair. Hence the distance between the two motors in the tangential direction can be calculated as


where , . When the motor is walking, it behaves like a linear spring with rigidity by exerting equal and opposite forces


As a result, the two MTs slide past one another undergoing polarity sorting. Following slender-body theory keller76 (), the MT speed is given by , leading to


where , and is the fluid viscosity. We seek the time-dependent solutions of the form . The coefficients and can be solved as


leading to


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 anti-aligned 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 small-angle assumption a reasonable approximation. Similar to the perfectly-aligned case, the positions of the two motors can now be written as:


where , and . So the relative distance becomes


where , . The motors exert tangential force . Following the same procedure, we seek solutions of the form , yielding


Then the relative moving speed of the two MTs becomes


When , Eqs. (43) and (44) exactly recover the solutions in (40) for the perfectly-aligned 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 center-of-mass positions and polar orientation vectors (), evolved through a Smoluchowski equation


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 second-moment tensor which arises generically in capturing active stresses produced by active suspensions simha02 (), the (trace-free) order parameter tensor , with or 3 the spatial dimension, and the fourth moment .

Slender-body theory yields the forces each rod exerts on the fluid, and hence the volume-averaged stress batchelor70 () by polarity sorting can be calculated. If the cluster occupies a volume , the induced extra stress tensor from anti-aligned sorting is . Here is proportional to fluid viscosity , and with the signed distance between the center-of-masses 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 plus-end of the MTs. This is seen in the BD-kMC 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 anti-aligned pair stresslet strength can be derived as . When taking as , we extract the value of from the BD-kMC simulations.

While active motoring between polar-aligned MTs yields little MT mobility, the BD-kMC simulations show that it does yield an extensile stress. However, unlike polarity sorting we lack a simple first-principles 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 polar-aligned stresses are of the same order (Fig. 2h) we assume the form . Comparison with the BD-kMC 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


To non-dimensionalize 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 Maier-Saupe 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 anti-aligned and polar-aligned 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


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 flow-induced constraint forces on MTs and steric interactions ESS2013 ():


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 half-spaces 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


where is a regular 3D spatial gradient operator. Under Fourier transform in , the above equations can be written as:


where is a 2D wave-vector, 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


where and is a 2D unit wave-vector. 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:


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 pseudo-spectral method to solve the Smoluchowski equation (45) and the fluid flow in a coupled manner.

iv.3 Flow, polarity and defects

Figure 5: Snapshots of streaming MT nematics on a liquid-liquid interface. The active stress magnitudes are chosen as and . (a) The background fluid velocity vector field superposed upon the colormap of the associated vorticity. (b) The nematic director field superposed on the colormap of the scalar order parameter (twice the positive eigenvalue of the tensor ). Disclination defects of order appear in localized regions of low order. Two defects are marked by an open circle () and a square (). The arrow at right marks a pair of annihilating defects, while the arrow at left identifies an “incipient crack” from which a defect pair is about to emerge. Here is a calculated characteristic length between the cracks. Principal eigenvalues of the active stress due to polarity sorting (, c) and motor relaxation (, d).

Assuming 2D periodic boundary conditions, we have simulated our active polar nematic model over long times, using Eqs. (45-49) as well as the velocity-stress kernel (53). For the simulations shown here, we choose from , and fix (i.e., aspect ratio 10), , , and (estimated from the BD-kMC 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 well-localized 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 high-force 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.

Figure 6: Time sequential snapshots of the nematic director field for nucleation (a) and annihilation (b) of defect pairs, where and are fixed. The (dimensionless) time spacing between frames is 5. (c) Polarity field associated with a motile defect and an incipient crack on the bottom. (d) Relative speed of the two oppositely charged defects, as well as the mean flow speed near this defect pair, as function of dimensionless time . In (a), (b) and (c), the color shows the scalar order parameter, plotted with the same scale as Fig. 5b.
Figure 7: Dynamics of the polarity field, the polarity-dependent active stresses, and the predicted dynamics of a photobleaching experiment. (a-c) Results from the simulation shown in Fig. 5. (a) The polarity vector field superimposed upon its magnitude (the local polar order). Circular areas labeled A and B mark regions of high and low polarity, respectively. (b,c) Polarity-dependent active stress magnitudes, showing principal eigenvalues of the active stresses due to polarity sorting (, b) and motor relaxation (, c). In (a-c), the stress magnitudes are chosen as and . For comparison, (d-f) shows the polarity field and the polarity-dependent active stress fields when choosing and . In (a-f), positions of -order defects are marked by open circles. (g) Schematic of predicted dynamics for a bleached spot of high nematic order in a region of high polar order (area A), and in a region of low polar order (area B). Arrows represent MTs with arrowheads denoting plus ends. In panels (h) and (i) these predictions are borne out by simulations of photobleached spots in areas A and B, respectively.

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 positively-signed defect moves away faster and roughly along its symmetry axis. Following annihilation of an oppositely-charged 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 shock-like 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 polar-specific fluxes and active stresses, the polarity field Tjhung11 (), polarity-dependent 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 polar-aligned (anti-aligned) stress is large in regions of high (low) polar order (Fig. 7b,c). The anti-aligned 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. 7d-f). The ratio between the anti-aligned and polar-aligned stress is still close to the ratio . However, since the sign of the polar-aligned 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 high-intensity laser light to bleach the fluorescing molecules on the corresponding MTs AKSEW1976 () (Figs. 7g-i). In a small high-polarity 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 low-polarity region of high nematic order (marked B in Fig. 7b), strong polarity sorting of anti-aligned 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 nematically-aligned 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: