Towards the cellular-scale simulation of motor-driven cytoskeletal assemblies

The cytoskeleton -- a collection of polymeric filaments, molecular motors, and crosslinkers -- is a foundational example of active matter, and in the cell assembles into organelles that guide basic biological functions. Simulation of cytoskeletal assemblies is an important tool for modeling cellular processes and understanding their surprising material properties. Here we present aLENS, a novel computational framework to surmount the limits of conventional simulation methods. We model molecular motors with crosslinking kinetics that adhere to a thermodynamic energy landscape, and integrate the system dynamics while efficiently and stably enforcing hard-body repulsion between filaments -- molecular potentials are entirely avoided in imposing steric constraints. Utilizing parallel computing, we simulate different mixtures of tens to hundreds of thousands of cytoskeletal filaments and crosslinking motors, recapitulating self-emergent phenomena such as bundle formation and buckling, and elucidating how motor type, thermal fluctuations, internal stresses, and confinement determine the evolution of active matter aggregates.


Introduction
Living systems are built hierarchically, where smaller structures assemble themselves into larger functional ones. Such organization is fundamental to life, where it is seen across scales from molecules to organelles to cells to tissues to organisms. An example is the cellular cytoskeleton, made up of polymer filaments (and other accessory proteins) crosslinked by motor proteins that exert forces by walking processively along filaments (Howard et al. (2001)). Cytoskeletal assemblies such as the cortex, mitotic spindle, and cilia and flagella, underlie cell polarity, division, and movement (Bornens (2008); Barnhart et al. (2015); McIntosh (2016); Pollard and O'Shaughnessy (2019)). Cytoskeletal components have been reconstituted outside of cells to study self-organization (Nedelec et al. (1997); Foster et al. (2015)) and to create new active materials (DeCamp et al. (2015)). Understanding how cytoskeletal structures assemble from their molecular components remains challenging, in part because of the variety of motors and crosslinkers with different behavior. Improved understanding of the cytoskeleton would allow us to predict how molecular perturbations change cell behavior and to design new complex and adaptive materials (Li and Gundersen (2008); Fletcher and Mullins (2010); Needleman and Dogic (2017)).
Here we describe aLENS, a framework of computational methods and software designed to more efficiently and accurately simulate large cytoskeletal systems (Fig. 1). Since motor proteins must bind, crosslink, and unbind from filaments to evolve such systems, aLENS simulates motors as traversing a (well-defined) free energy landscape Lamson et al. (2021). This prevents artificial energy flux during crosslinking and maintains detailed balance in the passive limit. As motors crosslink filaments, the spacing between filaments is on the order of the length of motor proteins (10-100 nm) (Fig. 1A), comparable to the filament diameter. Therefore, steric interactions between filaments occur frequently and must be treated carefully to avoid unphysical filament overlap, stress and deformation (Fig. 1B). Most other cytoskeletal simulation methods implement a repulsive pairwise potential between filaments, but this requires a small timestep for hard potentials because of the instability of timestepping methods (Heyes and Melrose (1993)). Therefore, potential-based models limit simulations to short timescales. To circumvent this limitation, here we utilize our recently developed constraint method to enforce hard-core repulsion between particles (Anitescu et al. (1996); Yan et al. (2019)). We further develop constraint-based modeling by introducing a related method to treat stiff spring forces due to crosslinking motors. Both steric interactions and crosslinking forces are incorporated in a unified implicit solver. This approach ensures numerical stability of the method and allows for timesteps two or more orders of magnitude larger than currently available. Additionally, aLENS is parallelized with OpenMP and MPI to reach length and timescales comparable to to those of experiments ( Fig. 5 and 7).
As an illustration of aLENS, Fig. 1C (and movie video1.mp4) shows a simulation of 3200 microtubules within a spherical volume driven by 9600 motors that, when bound, walk to the microtubule minus-end (modeling the activity of dynein). Though the microtubules are initially unorganized (C1), the combination of motor crosslinking and walking causes the microtubule minusends to contract into the center of a large aster (C2). The motor-driven steric interactions between filaments, however, eventually fragment this into smaller asters and bottle-brush-like structures (C3,C4). This simulation displays the complex interplay between steric and crosslinking forces in determining the dynamics and steady state configurations of cytoskeletal materials.

Methodology
In this work we model filaments as rigid spherocylinders. (While not presented here, flexible filaments can be modeled within our framework as segmented, jointed filaments; See Appendix 8.) Crosslinking motors are modeled as Hookean spring tethers connecting two binding domains referred to as heads, with steric interactions between motors neglected.
As outlined below, our algorithm performs 3 tasks sequentially at every timestep: motor diffusion and stepping, motor binding and unbinding, and filament movement. The major computational challenges arise in task 2, computing binding and unbinding while maintaining realistic macroscopic statistics, and in task 3, updating filament position while overcoming stiffness constraints and maintaining steric exclusion. The timestep is determined by the shortest characteristic timescale in the simulated system (filament collision, motor binding/unbinding kinetics, and filament motion). All other degrees of freedom (e.g., internal conformational changes of motor binding heads) are assumed to occur on shorter timescales. Figure 1. A: aLENS simulates dynamics of rigid filaments crosslinked and driven by motors, thermal fluctuations, and steric interactions. Motors bind to, unbind from, and walk along filaments. B: To achieve high efficiency, aLENS computes motor forces implicitly, and steric interactions through a novel geometric constraint method that avoids filament overlaps. C1-C3: Example simulation of microtubules organized into asters by minus-end-directed motors. The 300 s Brownian simulation contains 3200 microtubules, each 1 µm long, inside a sphere of radius 3 µm. The initial position of each microtubule is random and the half of each filament on the minus-end is colored pink. Three end-pausing dynein motors are fixed at the minus-end of each microtubule and walk toward the minus-end of any microtubule they crosslink. After initial contraction into a single large aster, strong steric interactions in the aster center break up the system into several smaller asters and a bottle-brush structure. C4: Motors are highly concentrated at the centers of asters.

Crosslinking motor diffusion and stepping
Each unbound motor executes Brownian motion independently. Each bound motor updates information on the filament to which it is attached, following filament movement in the previous timestep. During the motor movement step, singly bound motors move Δ and doubly bound motors move Δ along the filaments. Here is the motor stepping velocity that depends on force on the motor head (Gao et al. (2015a)): ( proj ) = max 0, min(1, 1 + proj ∕ stall ) , where proj is the projection of tether force along filament in the stepping direction. As typically found experimentally, this stepping model means that if proj is assisting stepping, the velocity saturates at ; while for proj hindering stepping, stepping is halted when proj = − stall .

Crosslinker binding and unbinding
In filament networks, the spatial variation of unbound and bound motors is integral to network self-organization. For example, crosslinking proteins concentrate in volumes with high filament densities, producing ripening effects as passive crosslinkers are depleted from the bulk (Weirich et al. (2017)) (e.g. see Fig. 1C). Furthermore, if motors or crosslinkers bind, unbind, or diffuse at rates not set by free energy barriers, the system's energy and/or entropy can be artificially elevated or lowered, changing the system dynamics and steady-state configuration. Entropic forces bundle and increase overlaps among crosslinked filaments (Lansky et al. (2015); Gaska et al. (2020)), and free-energy-dependent binding kinetics contribute to organization of cortical microtubules (Allard et al. (2010)) and induce actin bundling (Yang et al. (2006)). Ad-hoc models, like those that attach crosslinking motors to filaments at a fixed length or randomly sample a uniform distribution to set the binding length, are unlikely to recover the force or final configuration of bundled filaments. For example, if passive crosslinkers only bind in a non-stretched configuration, they will not generate entropic forces that drive bundle overlap, as seen experimentally (Lansky et al. (2015)). Further, if crosslinkers are modeled as binding with a uniform length distribution and zero tether rest length, the contractile stress of networks will be overestimated, condensing filament networks with greater rapidity.
The assemblies of filaments/motors are assumed to explore an underlying free energy landscape, where all 'fast' degrees of freedom can be subsumed into an effective free energy that depends only on filament and crosslinking motor degrees of freedom. We require that our model correctly recapitulates the distribution and chemical kinetics of crosslinking proteins in the passive limit, i.e., when = 0 for the bound velocity of motor heads. We achieve this with a kinetic Monte Carlo procedure in which motor protein binding and unbinding events are modeled as stochastic processes. Transition rates recover the correct limiting (equilibrium) distribution by imposing detailed balance (Appendix 3). That is, we model binding and unbinding as passive processes, but it is in principle possible that certain such processes consume chemical energy.
To enforce the macroscopic thermodynamic statistics, including correct equilibrium bound-  )), we explicitly model each crosslinker as a Hookean spring connecting two binding heads labeled as or . Each crosslinker has 4 possible states: both heads unbound ( ), either or singly bound ( or ), or both heads (doubly) bound ( ). For each timestep Δ , we first calculate the rates ( ) at which each head ( and ) transitions from their current state to a new binding state (i.e. for the transitions ⇌ ( , ) ⇌ ). The transition probabilities are modeled as inhomogeneous Poisson processes with the cumulative probability function The transitions ⇌ ( , ) do not stretch or compress the tether and so do not depend on tether deformation energy. However, the transitions ( , ) ⇌ do account for tether deformation energy (Table 1).

Rate
Value Table 1. The transition rates between all possible states of a crosslinker ⇌ ( , ) ⇌ . ( , ) means either head or is bound but the other is unbound. All binding rates account for the linear binding density . in, ( , , ) is the length of filament with center-of-mass position and orientation inside the capture sphere with cutoff radius , relative to position of motor/crosslinker . The sum is over all possible candidate filaments . The unbound-singly bound transition ⇌ ( , ) is determined by the association constant a and the force-independent off rate o, . Similarly, the singly bound-doubly bound transition ( , ) ⇌ is determined by the association constant e and force-independent off rate o, . = 1∕( ) is the Boltzmann factor. ( ) in the in the ( , ) ⇌ transition rates refers to the tether energy of a motor ( ) = 1 2 xl − 0 2 . 0 is the free length of a motor, while is the length for computing the force when attached to filaments and at locations and : ( , , , , , ). The dimensionless factor determines the energy dependence in the unbinding rate. Both binding and unbinding rates must depend on and , such that the equilibrium constant recovers the Boltzmann factor exp[− ( )] For force-dependent binding models, the ( ) can be simply replaced by the tether force ( ). This is not used for results shown in this work, but implemented in the code.

Filament dynamics
We sought to develop a stable, large-timestep method for updating the position of filaments, subject to spring forces from crosslinking motors, steric interactions, and Brownian motion. This requires addressing two stability restrictions on the timestep Δ . The first arises in models that use a stiff repulsive pairwise potential to prevent filament overlaps. For example, the Lennard-Jones potential ∼ ( ∕ ) 12 − ( ∕ ) 6 , where is the separation between filaments, is so steeply varying that it requires small Δ for stability. As a result, soft alternatives such as a harmonic potential are often used (Nedelec and Foethke (2007)). These soft potentials allow partial filament overlaps, and may therefore lead to unphysical system dynamics and stresses (Heyes and Melrose (1993)).
We overcome these difficulties with a novel, linearized implicit Euler timestepping scheme, which extends on our previous work on enforcing non-overlap conditions (Yan et al. (2019)). This technique is inspired by constraint-based methods for granular flow (Tasora et al. (2013)). When collisions occur between filaments, the minimal distance between them attains Φ col = 0 with collision force col > 0. If not colliding, Φ col > 0 and col = 0. This mutually exclusive condition is called a complementarity constraint, written as 0 ≤ Φ col ⟂ col ≥ 0. If one crosslinking motor connects these two filaments, its length and force magnitude xl satisfy the Hookean spring model xl = − xl ( − 0 ), which is an equality constraint. We integrate the equation of motion such that these two types of constraints for all possible collisions and all crosslinking motors are satisfied. We briefly derive the method here, and all details can be found in Appendix 3. Because the method is specific to rigid particles with arbitrary shape, we shall use 'particle' and 'filament' interchangeably.
Each particle is tracked by its center location ∈ ℝ 3 in the lab frame and its orientation = [ , ] ∈ ℝ 4 as a quaternion (Delong et al. (2015)).
[ , ] are the scalar and vector parts of the quaternion, respectively. Using a quaternion to track the rotational kinematics of a rigid body is a standard computational approach due to its compact memory footprint (4 floating point numbers) and its singularity-free nature. The geometric configuration at time for all filaments can be written as a column vector with 7 entries: Similarly, we use the vectors  ,  ∈ ℝ 6 to represent the translational & angular velocities, and forces & torques of all particles, respectively. We relate  to  via a mobility matrix  ∈ ℝ 6 ×6 , dependent only upon the geometry , and relate  to via a geometric matrix : Because the biological filaments we consider mostly have lengths on the nm to µm scales and inertial effects can be ignored. In the following, the subscript refers to constraints, which includes both unilateral (with subscript ) and bilateral (with subscript ) constraints. For our problem, unilateral constraints refer to collision constraints while bilateral constraints refer to crosslinking motor constraints. The subscript refers to non-constraint. For unilateral constraints, we define the grand distance vector = Φ ,1 , Φ ,2 , ⋯ , Φ , ∈ ℝ , where each Φ , is the minimum distance between a pair of filaments. Similarly, for bilateral constraints we define the grand distance vector = ,1 , ,2 , ⋯ , , ∈ ℝ , containing the length , of the doubly bound motor . There are in total possibly colliding pairs of filaments and crosslinking motors. The force magnitude corresponding to these constraints are also written as vectors, = ,1 , ,2 , ⋯ , , ∈ ℝ and = ,1 , ,2 , ⋯ , , ∈ ℝ . The two types of constraints can be summarized as: Here and satisfy the complementarity (collision) constraints, while and satisfy the Hookean spring law. Here  ∈ ℝ × is a diagonal matrix consisting of all the stiffness constants, while 0 represents the rest length of every crosslinking motor.
Eqs. (4) and (5) define a differential-variational-inequality (DVI). This is solvable when closed by a geometric relation mapping the force magnitude and to the force vectors  and  : where  and  are sparse matrices containing the orientation norm vectors of all constraint forces (Anitescu et al. (1996);Yan et al. (2020) and Appendix 4). Next, we discretize this DVI using the linearized implicit Euler timestepping scheme with Δ = ℎ at timestep : The unknowns to be solved for at every timestep are the constraint (collision and motor tether) force magnitude , . This is a nonlinear DVI because +1 , +1 are nonlinear functions of geometry  +1 , although  +1 is linearly dependent on and . For a small timestep (ℎ → 0), this nonlinearity can be linearized by Taylor expansion, for example, +1 = + ℎ∇    . Then, this nonlinear DVI can be converted to a convex quadratic programming problem (Nocedal and Wright (2006)) (details in Appendix 4): subject to Here = [ , ] ∈ ℝ + is a column vector, and One way to understand the constraint optimization method is that the implicit temporal integration 'jumps' on a timescale that bypasses the relaxation timescales of unilateral and bilateral constraints (collisions and crosslinking motor springs). In the limit of motor tethers being infinitely stiff ( −1 → ), the quadratic term coefficient matrix is still symmetric-positive-semi-definite (SPSD) and the Eq. (8) is still convex and can be efficiently solved. Physically speaking, in this case the bilateral constraints degenerate from deformable springs to non-compliant joints.

Instantiation in a massively parallel computing environment
Our methods naturally lend themselves to high-performance parallel computing architectures. We utilize both MPI and OpenMP and use standard spatial domain decomposition to balance the number of motors and filaments across MPI processors. The motor update step samples the vicinity of every motor, where we use a parallel near-neighbor detection algorithm and update all motors in parallel. The most expensive part of the method is finding the solution to Eq. (8), because of its very large dimension, equal to the total number of close pairs of filaments plus the number of crosslinking proteins. We use a fully parallel Barzilai-Borwein Projected Gradient Descent (BBPGD) solver (Yan et al. (2019)) because the gradient ∇ = + is efficiently computed by one parallel sparse matrix-vector multiplication operation.
aLENS is written in a modular design using standard object-oriented C++ and is available on GitHub as discussed at the end of the Discussion section.

Verification and Benchmarks
To validate and benchmark aLENS, we first note that its collision handling approach has already been benchmarked for the pure-filament phase, and shown to accurately reproduce the equation of state and the isotropic-nematic liquid crystal phase transition of densely packed rigid Brownian rods (Yan et al. (2019)). This capacity to accurately compute the dense packing phase of fibers makes aLENS valuable to simulate many dense biological filament assemblies. The accurate treatment of steric interactions extends beyond other simulation methods and software, where steric interactions are often approximated by soft repulsive potentials or neglected.
We now further benchmark of aLENS by simulating mixtures of filaments and motors and directly comparing simulation results with experimental data. Although there are many parameters in our motor model, these comparisons don't involve fitting of model parameters to experimental data. Instead, we chose motor parameters as measured from experimental data (Scharrel et

Directed transport of microtubules by mixed active and inactive motors
We begin by verifying our motor model by reproducing results from experiments on directed microtubule transport (Scharrel et al. (2014)). As in the experimental system, the simulation begins with a fixed number of motors with one head attached to a fixed surface while the other head interacts with one microtubule. Some motor heads are active and can drive gliding of the microtubule, while other heads are inactive and behave as passive crosslinkers that hinder microtubule motion.
Here is the number of active motors and is the total number of motors (active and inactive). The microtubule velocity increases as ∕ increases from 0 to 1 in experiments (Scharrel et al. (2014)) and in our simulations. As shown in Fig. 2, our simulations quantitatively reproduce the experimental data. To achieve this agreement, we set the active motor velocity to 1.0 µm s −1 , so the sliding velocity at ∕ = 1 matches experiment. Apart from this one experimentally constrained velocity, there are no fitting parameters in our simulation (further motor parameters are in Appendix 2). In initial trial simulations, we found that changing the total motor number didn't noticeably affect the microtubule transport velocity. Therefore, for the results shown here we fixed = 100, similar to the experimental system. Since the transport trajectory is stable without stochastic noise, as shown in Fig. 2, there is no need to perform ensemble average to determine the transport velocity. Therefore, we ran 1 simulation for 10 s for each ratio ∕ .

Self-straining state of actively crosslinked microtubule networks
As an additional verification, we compare aLENS with results of recent experiments of Fürthauer et al. (2019) in which many-microtubule assemblies are densely packed into a nematic bundle and crosslinked by a large number of motors. In this heavily crosslinked nematic regime, microtubules are found to be transported by motors along the nematic director direction at a constant velocity in a direction determined by individual microtubule polarity. Experimentally, microtubule velocity was found to be independent of the local average polarity of the ensemble, as has been observed in extract spindles (Needleman et al. (2010)), and (over the range of experimental conditions) independent of motor density. This phenomenon of oppositely-oriented, constant velocity microtubule fluxes was referred to as 'self-straining motion', with the system interpreted as being composed to two polar microtubule gels whose inter-connecting motors pulled them past one another.
We simulate this experiment using 3000 model microtubules with = 0.5 µm. Initially the filaments are confined in a tube of diameter = 1 µm, randomly initialized with their orientations along the + (pink) and − (white) directions, and packed at about 30% volume fraction. The simulated system is periodic along the direction, with periodic tube length 3 µm. There are approximately 25 motors per microtubule according to the experimental estimates, and in our simulations we vary the motor-to-microtubule number from 10 to 30. There is no accurate measurement for the XCTK2 motor in these experimental conditions. Therefore, we used experimental estimates of 46 nm s −1 for the walking speed of NCD motors (Furuta and Toyoshima (2008) Fig. 3 shows that the straining velocity is largely independent of the number of motors and the local average polarity over the range simu-lated. Intuitively, the straining velocity is predominantly determined by the free walking velocity of the motors in limit of many cross-linkers. From our simulations, we find a straining velocity of approximately 26 nm s −1 , close to the experimental measurement of 18.6 ± 0.9 nm s −1 . is the number of XCTK2 motors per microtubule. We sample the local average polarity and straining velocity by inserting planes orthogonal to the -axis into the collected data, matching the photobleaching technique used in experimental measurement (Fürthauer et al. (2019)). For every sampling plane (e.g. the blue pane in the snapshot), we choose five sample points symmetrically on this plane and draw a square sampling window with edge length 0.2 µm around each sample point. For each sampling window, we compute the average polarity along the -axis for all microtubules intersecting this sampling window at a given time. We then compute the velocities, averaged over 10 s (a duration chosen to match the experimental timescale), of microtubules intersecting each sampling window and moving along the + and − directions. + and − are computed from those two groups for each sampling window. The straining velocity is computed as = + − − . Therefore, for every sampling window at each sampling timestep we have a pair of data values , . The right three panels show the joint probability distribution of ( , ) computed from 900,000 sampling planes for each simulation, for = 10, 20, 30, respectively.

Large scale parallelization efficiency
Simulation of cellular-scale cytoskeletal assemblies requires methods that can reach large system sizes and timescales. Therefore, we developed aLENS to efficiently utilize modern high performance computing resources. Millions of objects and constraints can be simulated with aLENS. Fig. 4 shows detailed parallel efficiency measurements for one large-scale test case, similar to that in Fig. 6, but more than 10 times larger. Here we track 1 million microtubules and 3 million motors for 100 timesteps. The performance is benchmarked on a cluster interconnected with infiniband and each node has two AMD EPYC 7742 CPUs, each having 64 cores at 2.5GHz. We launched hybrid MPI+OpenMP jobs such that each MPI rank has 16 OpenMP threads. On average at each timestep the constraint optimization solver handles approximately 8 million collision and doubly bound motor constraints. The number of constraints changes at every timestep due to a variable number of collision pairs and to stochastic binding and unbinding of motors. We achieve nearly ideal linear speed up as the number of cores increases ( Fig. 4). At 1536 cores, the efficiency remains at 93% and each timestep takes less than 1 second, making it possible to track such large systems on experimental timescales (a few seconds) within days or weeks of computing time. More importantly, the constraint optimization allows a Δ that is one or two orders of magnitude larger than conventional pairwise potential methods. For the system simulated in Fig. 4, aLENS can reach 1 s physical time per day, using a timestep size of 1.0 × 10 −5 s.

Results
Here we illustrate the ability to use aLENS to study the interplay between microscopic dynamics and macroscopic order in active cytoskeletal assemblies. The specific examples shown here are the formation and extension of a band of microtubule bundles, polarity sorting of short microtubules . Strong scaling (fixed system size while increasing number of cores) efficiency of a system similar to but more than 10 times larger than that shown in Fig. 6, comprising 1 million microtubules and 3 million motors. There are in total approximately 8 million constraints per time step, which is changing from step to step because collision pairs are changing and crosslinkers are stochastically binding and unbinding. The simulation is run for 100 computing steps with 1 data-saving step and the average per-step wall-clock time is shown in the figure. on a spherical shell, the development of asters with and without thermal fluctuations, and the effect of confinement on assembling microtubule-motor mixtures. For the results presented here, all simulations were conducted in solvent with viscosity = 0.01 pN s µm −2 at room temperature, using a fixed timestep Δ = 10 −4 s unless otherwise stated.

Bundle formation and buckling in a filament band
Microtubules driven by crosslinking motors can bundle; sliding of microtubules within the bundles causes them to fracture dynamically (Sanchez et al. (2012); Foster et al. (2015); Roostalu et al. (2018)). We study such phenomena through a large-scale simulation of 100,000 filaments modeling microtubules and 500,000 minus-end-directed motor proteins modeled after dynein; (Fig. 5). Motor crosslinking drives contraction of initially disordered, bundled filaments ( Fig. 5A and B). Aligning steric and crosslinking forces drive the system into a series of well-aligned bundles spanning several filament lengths (Fig. 5C, see movies video2.mp4 and video3.mp4). The motors slide filaments parallel to each other, generating macroscopic extensile motion. Later, the extended network buckles and fractures (Fig. 5C).
The macroscopic stresses and dynamics depend on the spatial organization of filaments and motor-driven sliding. To characterize this, we measure the joint probability distribution of the local nematic order parameter local and the number of neighboring filaments crosslinked to a filament (Fig. 5D). While the network contracts, the distribution of doesn't change significantly because the number of motors per filament and the maximum number of neighboring filaments within a densely packed structure remain roughly constant. As filaments align, they become nearperfectly nematic ( local ≈ 1), although less-ordered regions occur between aligned bundles of different orientations (Fig. 5C1, D2).
Inside the bundles, filament sliding by motors leads to transport along the local nematic director. Projecting filament trajectories onto the lab-frame -axis, we observe left-and right-moving Figure 5. Results for the bundling-buckling simulation of 100,000 microtubules and 500,000 dynein motors in the periodic simulation box of 600 × 10 × 10 µm. Brownian motion of microtubules is turned off. Each dynein has one non-motile head permanently attached to a microtubule and the other motile head walks processively with maximum velocity 1 µm s −1 . If bound, the motile head moves towards the microtubule minus-end, and detaches upon reaching it. Detailed parameters for this motor are tabulated in the Appendix 2. Every microtubule has 5 dynein motors permanently attached to randomly chosen, fixed locations along the length. The initial configuration of microtubules is randomly generated, with their orientations sampled from an isotropic distribution and centers uniformly distributed within a cylinder of length 600 µm and diameter 0.3 µm. The motile heads of all dynein motors are unbound initially. A, B, and C: The bundle at = 0 s, 4 s, and 7 s. Microtubules are colored by their local nematic order parameter , being the unit orientation vector of each microtubule pointing from the minus to the plus end, and the Kronecker delta tensor . The average ⟨.⟩ is taken over each microtubule plus all microtubules that are directly crosslinked to it by dynein motors. A1, B1, and C1: Zoom-in views of the small region marked by red box in A, B, and C. C2: The same region in C1 but colored by , the number of microtubules averaged over when computing local . D1 and D2: The joint probability distributions local and for each microtubule for the entire systems at = 0.1 s, when the dyneins crosslink microtubules but microtubules barely move from initial configuration, and at = 7 s, when the bundle is nematic. filaments that speed up early in the simulation, and then maintain constant average velocities at later time ( > 4 s in Fig. 5E), as filaments align due to steric and motor forces (Fig. 5F). Note that velocity and stresses plateau only when the nematic order saturates. The filament motions created by motors cause the densely-packed filaments to collide often, creating a net extensile stress along the bundles' axes ( Fig. 5F). However, the fixed simulation box size hinders the networks' elongation, causing the bundles aligned with -axis to buckle due to the net extensile stress (Fig. 5F, see movies video2.mp4 and video3.mp4). In contrast, bundles not aligned with the -axis are not constrained and so evolve into straight spikes. This misalignment of bundles is seen as a small net stress in the , -directions for ≥ 4 s (Fig. 5F).

Polarity sorting in a spherical shell
Crosslinking motors on antiparallel filaments drive polarity sorting, which transports filaments to regions of like polarity. This has been well-studied on a planar periodic geometry, e.g. (Gao et al. (2015b)). Here we use aLENS to examine the effect of confinement geometry on polarity sorting (Fig. 6). The geometry is designed to explore the polarity sorting phenomena where initial filament alignment occurs in a spherical geometry and significantly affects the dynamics and steady state of the system. In this simulation, 100,000 filaments with aspect ratio ∕ f il = 10 are confined between two closely spaced concentric spherical shells at 40% volume fraction. The shell gap is Δ = 0.102 µm, shorter than the filament length, with Δ ∕ f il ≈ 4 so filaments can move over each other in a restricted way. The filaments are initialized such that the nematic directors are along the meridians everywhere. 200,000 motors, modeled after kinesin-5 tetramers, drive relative filament sliding (Fig. 6A). Brownian motion is modeled at room temperature 300 K and timestep Δ is set to 1 × 10 −5 s. Motors move toward minus ends of bound filaments at = 1.0 µm s −1 . Once they reach the minus ends, they immediately detach.
Motors walk along the filaments, driving sliding of antiparallel filaments (Fig. 6B). This leads to polarity-sorted regions at the north and south "poles" of the sphere, meaning that the filament orientation on average points toward the poles. Filaments with reversed initial polarity are transported to the equatorial region (Fig. 6C1). In contrast to the planar geometry (Gao et al. (2015b)), we did not observe the formation of polar lanes with boundaries between polarity-sorted regions approximately parallel to the polarity direction. Instead, on the sphere the boundaries between polarity-sorted regions are approximately orthogonal to the polarity directions, as more clearly illustrated by plotting the polarity divergence (Fig. 6C1).
Motors also accumulate in some regions according to the filament polarity (Fig. 6C1). These motor accumulation regions are actually regions where the divergence of filament polarity field is positive, meaning areas of overlap of filament minus-ends ( Fig. 6C2, C5, G). This accumulation is illustrated by the positive correlation between motor density and ∇ ⋅ at = 4 s in Fig. 6D. Furthermore, motor accumulation regions appear to show slightly lower filament volume fraction ( Fig. 6C2 and C4), as shown in Fig. 6F. These correlations can be understood through the behavior of crosslinking motors near filament ends (Fig. 6G). Once polarity sorted regions of filaments form, as the blue arrows represent, ∇⋅ > 0 in regions where minus-ends meet minus-ends and vice versa in regions where plus-ends meet plus-ends. Minus-end directed motors accumulate in regions with ∇ ⋅ > 0, while plus-end motors accumulate in regions with ∇ ⋅ < 0. Once motors accumulate, they may attach to both minus ends and push them away such that the distance between minus ends is the length of motors. As a result, the volume fraction of filaments in that region is below average.
In contrast, if the motors stop walking but do not detach when they reach the minus ends (endpausing, EP), the filament network contracts (Fig. 6E1-5) with volume fraction increases from 40% to 60% and eventually freezes at = 0.27 s. We observe neither substantial polarity sorting nor motor accumulation. This indicates that the ability of motors to continuously walk, without end-pausing, is crucial to effective polarity sorting. Figure 6. Results for the polarity sorting simulation in a spherical shell. Initially, 100,000 0.25 µm-long filaments modeling microtubules and 200,000 motors modeling crosslinking kinesin-like proteins are placed between two concentric spherical shells with radii in = 5 µm and out = 5.102 µm, to maintain the volume fraction of filaments between these two shells at 40%. Initially, all filaments are evenly distributed on the spherical shell, with their orientation randomly chosen to be either ± at each point, where is the polar basis norm vector of spherical coordinate system. The pure filament system is relaxed for 1 s to resolve the overlaps in the initial configuration. Afterwards at = 0, 200, 000 motors are added to the system homogeneously distributed between the two shells. Sample points are evenly placed to measure the statistics by averaging the volume within 0.25 µm from each sample point. A: The configuration at = 0. Filaments are colored by their polarity, while motors are colored as black dots. Only randomly selected 10% of all motors (same after) are shown in the image to illustrate the distribution. B: Randomly selected trajectories of filaments from = 0 s to 1 s. Trajectories are colored by time. It is clear that filaments move along the meridians. C1-C5: Configuration and statistics at = 4 s. C1: The filaments and motors. Motors clearly concentrate in some areas. C2: The motor number density, i.e., number of motors per 1 µm 3 . C3: The nematic director field (shown as black bars) and the nematic order parameter . C4: The filament volume fraction. C5: The divergence of polarity field ∇ ⋅ non-dimensionalized by filament length, i.e., change of mean polarity per filament length. D: The development of the correlation between motor number density ∕ ave and the polarity divergence field, at different times of the simulation. Clearly high ∕ ave are correlated with positive polarity divergence. E1-E5: Configuration and statistics at = 0.27 s for a comparative simulation where motors have end-pausing, arranges in the same style as C1-C5. This case shows significant contraction instead of polarity sorting as filaments are pulled away from the north and south poles and the overall volume fraction significantly increases to approximately 60%. The structure becomes densely packed and does not significantly evolve further. F: The correlation between motor number density ∕ ave and the local filament volume fraction. For the polarity sorting case at = 4 s the motor number density correlates with low filament volume fraction. This is not seen in the end pausing (EP) case. G: A schematic for the correlations shown in D and F.

Aster formation in bulk
Aster formation is driven by motor pausing at ends of rigid filaments (end-pausing). Previous work has focused on how motor biophysics affects aster formation (Belmonte et al. (2017); Roostalu et al. (2018)). An additional contributor to aster formation may be thermal fluctuations, which are difficult to tune experimentally but can be easily modulated in simulations (Fig. 7). To examine this, we simulated 40,000 filaments and 80,000 processive, minus-end-directed, end-pausing motors starting from the same spatially uniform and orientationally isotropic random configuration (Fig. 7A). In one version of the model, we included thermal fluctuations that drive filament motion ( Fig. 7D and movie video4.mp4), while in the other thermal fluctuations of filaments were neglected ( Fig. 7E and movie video5.mp4). The resulting structure of the system is significantly different in the absence of filament thermal motion, showing that thermal fluctuations influence the asters' shape, structure, and ultimate spatial organization. With filament thermal motion, a number of dispersed, spherically symmetric, dense asters form. By contrast, in the absence of thermal motion the number of asters is larger and more regularly spaced, but their shape is more irregular and they contain fewer filaments (Fig. 7D vs E).
These differences are clear in the radial distribution function of filament minus ends, which are clustered by motors paused at filament ends (Fig. 7B). On large length scales, the radial distribution reflects larger and denser asters for the simulation with thermal fluctuations that drive filament movement. In simulations of both cases, two prominent peaks appear in the radial distribution funcation at small length scales = 25 nm = f il and = 78 nm = 0 + f il which correspond to scale on which filaments bind to or are crosslinked by motors, respectively (Fig. 7B,D2,E2). The relatively small peak between these two maxima correspond to filaments that are geometrically confined between two crosslinked filaments.
These differences arise from the fact that athermal filaments do not move unless driven by motors, which requires that two filaments are close enough to become crosslinked. This suggests that, at steady state, athermal aster centers are separated by twice the filament length. In contrast, with thermal motion filaments may diffuse ∼1 µm in 1 s. This allows filaments to diffuse until they are captured in regions of high motor density, such as aster centers. Furthermore, with thermal fluctuations the asters themselves diffuse, which leads to aster coalescence (Fig. 7D1). These observations and estimated lengthscale are quantitatively confirmed by analyzing the static structure factor of aster centers (details in Appendix 6), which shows that the athermal simulation has approximately 3 times more asters than the thermal case ( Fig. 7D vs E).
The differences in the dynamics of aster formation are also reflected in stress measurements (Fig. 7C), where the more crowded filament configurations of the thermal case produces a larger stress throughout the simulation. In both cases the motor-induced stress Π initially increases quickly, reaching a peak at roughly = 4 s ∼ 5 s, similar to the behavior during bundle contraction shown above (Fig. 5F), before declining. The average time required for motors to walk to filament ends, = ∕ ≈ 5 s, determines the initial contraction timescale. After reaching minus ends, motors pause and relax toward their equilibrium lengths. As a result, both the motor and collision stress grow in magnitude as more motors accumulate at minus ends.

Confined filament-motor protein assemblies
Confinement of cytoskeletal structures plays an important role in cells, where the cytoskeleton is spatially constrained by membranes, organelles, and other cellular structures. Whereas in the previous examples we studied open periodic geometry, here we show results of cylindrical confinement. The microtubule motor system is constrained inside a cylinder with periodic boundary conditions at the cylinder ends. The impermeable boundary of the cylinder surface to motors and filaments was implemented by our complementarity constraints.
Similar to the previous bulk cases, motors move filaments to create high-density crosslinked filament aggregates that coexist with a relatively low density vapor of non-crosslinked filaments. In bulk systems as shown above and in previous work, end-pausing motors drive aster formation Figure 7. Results for the aster formation simulations with Brownian motion of simulated microtubules turned on (BMT) and off (NBMT). Initially, 40,000 0.5 µm-long filaments modeling microtubules and 80,000 motors modeling crosslinking kinesin-like proteins are placed in a periodic cubic box of 10 × 10 × 10 µm with uniform distribution. Filament orientations are isotropic and motors are all in the unbound state. Motors are assumed to have two minus-end-directed walking heads with symmetric properties. They are assumed to pause when they reach the minus end of filaments until detaching. Detailed parameters are tabulated in Appendix 2. A, D, and E: Simulation snapshots. Each filament is shown as a cylinder colored in half pink (minus end) and half white (plus end). A: The initial configuration for both NMT and BNMT cases. Each motor is colored as a green dot. D and E: The snapshot for both cases at = 35 s. D1-2 and E1-2: Expanded views of a aster core for D and E. Only doubly bound motors are shown in D and E (in green color), and in D1-2 & E1-2 (colored by the spring force). Negative values mean the crosslink forces are contractile (attractive). B: The radial distribution function (RDF) ( ) for the minus ends of all filaments at = 0.1 s (dashed lines) and = 35 s (solid lines). The first peak of ( ) at = 25 nm corresponds to close contacts between filaments. The second peak of ( ) at = 78 nm = 25 nm + 53 nm corresponds to the minus ends of filaments crosslinked by motors whose rest length is 53 nm. Blue and red lines are results for the BMT and NBMT cases, respectively. C: The collision (solid) and crosslinker (dashed) pressure for BMT (blue) and NBMT (red) cases. Pressure is defined as the trace of the stress tensor: Π = 1 3 Tr . The collision pressure Π is positive (extensile), and the motor pressure Π is negative (contractile). The inset plot shows the pressure for the NBMT case in the initial stage of the simulation. The black dashed lines mark the time = 4 s. because crosslinking motors pull filament ends together. A confining cylindrical boundary strongly modifies the conformation of these aggregated structures (Fig. 8). These simulations used 0.25 µm long filaments at a fixed packing fraction ( = 0.16), confined in two cylinders with diameters cyl = 0.25 µm and 0.75 µm.
For a small-diameter cylinder where one filament length can fit across the cylinder ( cyl ∕ = 1, cyl = 0.25 µm), the cylinder is too narrow for asters to form. Instead, motor sliding and endpausing drive the filaments into polarity-sorted bilayers (PSBs, Fig. 8A and movie video6.mp4). A single polarity-sorted bilayer contains a central interface of highly-crosslinked filament minus-ends between two antiparallel polar layers of filaments (Fig. 8A2-3). At steady state, the system consists of individual PSBs separated by low-density vapor regions containing few motors. As expected, the local nematic order parameter local ( ) nearly reaches 1 within PSBs. Even the the vapor phase is close to nematic local ≈ 0.6 (Fig. 8A4), due to the strong confinement effect. Next we increased the diameter of the cylinder to cyl ∕ = 3 ( cyl = 0.75 µm) to weaken the confinement ( Fig. 8B and movie video7.mp4). Here the polarity-sorted bilayers are not present, because the larger cylinder diameter allows filaments to reorient and organize into bottle-brushlike aggregates (BBs). In the bottle brushes, filament plus ends are oriented radially outward from the cylinder axis, forming a hedgehog line defect capped by half asters (Fig. 8B2). Motors become highly concentrated along the line defects at the center of the cylinder (Fig. 8B3). The radial hedgehog structure of BBs is evidenced by a negative local nematic order parameter (Fig. 8B4, blue line). The splayed nature of the BBs produces a lower relative packing fraction of ∼ 2.5 times the vapor when compared to the PSBs (Fig. 8B4, red line).

Discussion
We designed aLENS to (i) model crosslinking motor kinetics conforming to an underlying free energy landscape, (ii) circumvent the timescale limitation imposed by conventional explicit timestepping methods, and (iii) efficiently utilize modern parallel computing resources to allow simulation of cellular-scale systems. This efficient framework allows both modeling the individual building cytoskeletal building blocks (filaments, motors) and gathering mesoscale statistical information such as stress and order parameters from a large system. This multiscale capability will make it possible to directly compare simulations with experimental observations on mesoscopic and macroscopic scales over timescales from seconds to minutes.
The aLENS framework is not limited to a specific motor model. Because of the modular design of the motor code, the motor model can be extended to include additional physics such as forcedependent binding and unbinding rates, or even entirely replaced, say, with a passive crosslinker or other model. Dynamic instability and branching of cytoskeletal filaments can also be integrated with the constraint minimization problem, as we showed previously in modeling the division-driven growth of bacterial colonies (Yan et  aLENS has been open-sourced on GitHub: https://github.com/flatironinstitute/aLENS and precompiled binary executable is available on DockerHub: https://hub.docker.com/r/wenyan4work/alens. Our GitHub documentation provides a clear roadmap for developing additional user-specific modules.

Summary of videos
Here is a list of videos for this manuscript.
Video 1 (Figure 1 Video 1): Contraction and break-up of simulated microtubule asters. The simulation details are described in Fig. 1.    (2020). Modeled crosslinking proteins in solution bind to one filament and then crosslink two filaments (Fig. 1). In dense filament networks, the spatial variation of unbound proteins play an important part in the network's reorganization. To account for inhomogeneous concentrations, we explicitly model unbound crosslinkers and develop a method that reproduces one head bound and doubly bound distributions consistent with a mean-field model (Appendix 3). All binding and unbinding rate calculations are summarized in Table 1.
Unbound crosslinking proteins rapidly diffuse in the surrounding fluid until a head binds to a filament. Heads of modeled crosslinking proteins in solution bind to filaments described by the reversible chemical reaction H + B HB, where H is a head and B is a binding site on a filament. The association constant a of the heads to binding sites is described by the equilibrium equation where [X] defines the concentration of substance X. The association constant has units of inverse molarity, and relates the the on-and off-rate constants on, and of f, . Unbound crosslinkers are modeled as diffusing points with center of mass positions ( ) and diffusion constants . The heads of a crosslinker have spatial-and time-dependent concentrations [H] = ( , ). The head binding rate is the volume integral over the product of the on-rate constant, binding site density, and crosslinker concentration on, ( ) = on, where on, = a of f , and is the linear binding site density along filaments. The lab position along the th filament ( ) is parameterized by . The binding probability in a timestep Δ is an in-homogeneous Poisson process with the cumulative probability function We assume the tight binding limit on, ≫ of f, and do not consider multiple binding and unbinding events of one crosslinker during a timestep Δ . The average number of multiple events may be calculated from binding parameters and the timestep allowing one to set a probability thresholdLamson et al. (2021). Heads of the same crosslinking protein are forbidden to be bound to the same filament at the same time.
To describe ( , ) during a timestep, we first consider a crosslinking protein with two heads connected by a flexible but relatively stiff polymer tether with length . Because of the tether's stiffness, the radius of gyration of an unbound protein is assumed to be = ∕2. The binding heads at the tether's ends move by the tether's rotation and translation. Depending on the timestep's length, either rotation or translation will dominate the evolution of the head distributions.
For most biological crosslinking proteins, the rotational diffusion is fast compared to translational diffusion. When the crosslinker's center does not diffuse far from its position at the beginning of a timestep, i.e., √ 6 Δ < , rotational diffusion dominates and we approximate heads to be within a sphere of radius , = centered at . Realistically, the head distributions can vary within this volume but because is small compared to filament lengths, we approximate the head distribution as being uniform, i.e., ( , ) = (4 3 , ∕3) −1 (1− Θ(| | − , )), where Θ( ) is the Heaviside step function. For larger crosslinking proteins, more detailed spatial distributions may be calculated using freely-jointed or worm-like chain models.
For uniformly distributed heads, the head binding rate is where is the filament 's length segment within , . To account for cylindrical filaments with diameter f il , we augment the binding radius such that , = + f il ∕2. Since this scenario exists within a regime where the crosslinker or motor does not diffuse far from its initial position in a timestep, we approximate on, ( ) ≈ on, ( ), for ∈ [ , + Δ ).
However, it is uncommon that an unbound crosslinker or motor will diffuse less than in a timestep, and so we must account for the protein's translational diffusion. The diffusion equation models the mean spatial distribution of a unbound crosslinking protein's center which has the solution If the characteristic diffusion length √ Δ ≫ , then equation (14) underestimates binding ( Fig. 2A,D). The large diffusion distance also allows us to approximate the head distribution as matching the protein's spatial distribution, ( , ) ≈ ( , ). Substituting the binding rate equation (12) and the solution to the diffusion equation (16) into the integral of equation (13) gives For straight, rigid filaments, we take the volume and time integrals while reparameterizing | − | 2 by the crosslinker's perpendicular ℎ and parallel distances from a filament segment's center. This gives the linear binding probability density for a filament Integrating over gives the binding probability of one crosslinker head to a single filament. The total binding probability is then where is the number of filaments surrounding the crosslinker head. Calculating the binding probability from this function and ensuring that the protein unbinds so that detailed-balance is satisfied is computationally prohibitive. Instead, setting , to the root mean square diffusion distant √ 6 Δ and using the rate equation (14), we obtain a useful approximation to equation (19). This is computationally efficient and mitigates the low binding rates when diffusion or times steps are large (Fig. 2). We note that the accuracy of this approximation is dependent on the length of filaments in the simulation with longer filaments reducing the error from edge effects. Future work will focus on developing methods to more accurately reproduce the above binding distribution.
Once bound, a crosslinking protein's head unbinds with a constant rate If implementing equation (14), the protein unbinds into a uniform sphere of radius , . This ensures crosslinkers bind to and from regions in a way that satisfies detailed-balance. With one head bound, a crosslinker's tether deforms to bind its other head to adjacent filaments. Deforming a tether requires energy, implying crosslinking kinetics depend on tether deformation. For passive crosslinkers and motors with rapid kinetic rates compared to stepping rate, the ratio of binding and unbinding rates to and from a position on a filament is proportional to the Boltzmann factor of the binding free energy.
With one head bound to filament at position , the binding rate constant on, for binding to a location on filament is where o, = of f, ( , = 0) is the unbinding-rate of crosslinking proteins when no force is applied, e is a binding association constant similar to , and , is the free energy contribution from the tether Before crosslinking, the unbound motor head explores an effective volume bind centered around the bound motor head. Not considering steric interactions with filaments, this volume is the free head's position weighted by the Boltzmann factor integrated over all space.
We impose an integration cutoff radius , where the integrand becomes sufficiently small, making the factor consistent with a finite lookup table Lamson et al. (2021). The binding head's positional distribution must also satisfy the Boltzmann factor. We recover the proper binding distribution through inverse transformation sampling of equation (21) where and are the energy factor and characteristic length specifying the behavior of the energy-and force-dependent binding/unbinding, respectively. For values of < 0, you see catch-bond like behavior whereas values of > 0 exhibit slip bond behavior Walcott (2008); Edelmaier et al. (2020). This formalism can also be used to add in angle dependence. When we include an effective energy and/or force dependence, the unbinding rate becomes This does not change the final stored energy of either bound state but does effect the frequency at which the motors will switch between having one head bound and crosslinking. Table 1. The transition rates between all possible states of a crosslinker ⇌ ( , ) ⇌ . ( , ) means either head or is bound but the other is unbound. All binding rates account for the linear binding density . in, ( , , , ) is the length of filament with center-of-mass position and orientation inside the capture sphere with cutoff radius , relative to position of motor/crosslinker . The sum is over all possible candidate filaments . The unbound-singly bound transition ⇌ ( , ) is determined by the association constant a and the force-independent off rate o, . Similarly, the singly bound-doubly bound transition ( , ) ⇌ is determined by the association constant e and force-independent off rate o, with an additional factor bind , the effective volume explored by the unattached head while the motor/crosslinker is singly bound. Energy dependence in the ( , ) ⇌ transition rates is imposed by the Boltzmann factor that is a function of = 1∕( ), the tether length of the motor/crosslinker attached to filaments and at locations and ( , , , , , ), the characteristic length of the tether not under load , and the filament diameter f il . The dimensionless factor determines the energy dependence in the unbinding rate while the is the characteristic length that determines the force dependence. The latter is not used in the simulations of the main text but is implemented in the code base.

Mean-field theory for crosslinking proteins
We expand on our previous mean-field motor density model to include motors that have dissimilar heads, diffusion and walking in singly and doubly bound states, and a time-dependent homogeneous concentration of unbound crosslinking proteins Lamson et al. (2021). This last addition imposes the condition that the total number of proteins when all bound and unbound states are accounted for remains constant. This requires a system of equations with ( − 1) crosslinking densities , , , 2 singly bound densities and , and an unbound density to model all crosslinking proteins between filaments. By convention, ,

Constraint quadratic programming
In the main text we discussed specifically filament. In fact, our method is applicable to rigid bodies of arbitrary shapes. Here we derive the detailed equations. The configuration of each particle is tracked by its center location in the lab frame and its orientation = [ , ] ∈ ℝ 4 as a quaternion Delong et al. (2015). This is the vector component of the quaternion, not the unit orientation vector. There are other choices to specify the orientation, such as Euler-angles and rotation matrices, but we prefer quaternions for simplicity. The geometric configuration  for all filaments can be written as a column vector: which is a function of time: ( ). The translational and angular velocity , of all filaments can also be written as a column vector: Similarly we can write the force and torque , applied on all filaments as a column vector: The kinematic equation of motion 35 maps  to( ) = ∕ , via a geometric matrix .
 ∈ ℝ 7 ×6 is a block diagonal matrix, with one 3 × 3 and one 4 × 3 block for each particle: 3 is the 3 × 3 identity matrix, same for every particle. Each 3 block simply corresponds to the translational motioṅ = of each particle . Each ∈ ℝ 4×3 refers to the rotational motioṅ = of each particle, where for each : Here is the Levi-Civita symbol for cross-product in 3D space. The biological filaments we consider mostly have lengths on the nm to µm scales. At these scales, solvent viscosity dominates and inertia effects can be ignored, which is the so-called Stokes regime where the mobility matrix  maps the force  linearly to the velocity  :  includes collision force  between particle-particle and particle-container pairs, linker force between particle pairs  generated by doubly bound crosslinkers, Brownian force on each particle  generated by thermal fluctuations, and other externally applied forces  through gravity and electrostatic fields.
In principal, Eq. (35) together with Eq. (38) can be integrated directly because both  and  are functions of the geometry  and time only. However, this approach is usually impractical, because  or  is usually very stiff functions of the geometry. For example, the collision force  is usually computed by assuming a very stiff pairwise potential between filaments, such as the Lennard-Jones or WCA potential. This stiffness poses severe limits on the stability of all explicit temporal integrators. We discussed this problem in detail for collision forces  in our previous work on Brownian spherocylinders Yan et al. (2019) and rigid spheres in Stokes flow Yan et al. (2020). Instead of computing  using repulsive potentials, we imposed non-overlapping constraints on the geometry  while integrating Eq. (35).

Equation of motion with geometric constraints
In the following, the subscript refers to constraints, which includes both unilateral (with subscript ) and bilateral (with subscript ) constraints. Unilateral constraints refer to those inequality constraints, i.e., constraints imposed from one side, while bilateral constraints refer to equality constraints. In our system, unilateral constraints come from collisions and bilateral constraints come from doubly bound crosslinkers. The subscript refers to nonconstraint, i.e., physical components that are independent of the constraints.
For unilateral constraints, we define the grand distance function between every pair of particles as a column vector: where each Φ , is the minimal distance between particles with indices and . Similarly, we define the grand distance function for bilateral constraints: where each Φ , is the distance between two fixed points on particles and , respectively. Physically, Φ , , is simply the length of each doubly bound crosslinker. With this definition, there are in total unilateral and bilateral constraints in the system. In other words, there are in total possibly colliding pairs of filaments and doubly bound crosslinkers. Both kinds of constraints are functions of the system geometry, so we shall write them as () and () in the following when necessary.
The force magnitude between all pairs of particles for unilateral and bilateral constraints can be written similarly as column vectors: For each Φ , or Φ , , there is a corresponding force magnitude , or , , the (normalized) direction vector̂ = −̂ of this force, and the location and where this force is applied on the filament and respectively, as shown in Fig. 1. With norm vectors defined in this way, or is positive when the force is repulsive between two filaments.
For unilateral constraints and satisfy this complementarity condition: This condition means () and are orthogonal to each other, and all components of For bilateral constraints and satisfy this linear equality condition because they are modeled as Hookean springs:  ∈ ℝ × is a diagonal matrix, with the stiffness constant for each spring on its diagonal 1 , 2 , … . Obviously every constant is positive.
() and 0 represent the current and free length of every spring.
Both unilateral and bilateral constraints change over time, as particles move and springs attach to and detach from particles.
All combined together, we reach the equation of motion with geometric constraints: These equations are solvable when closed by a geometric relation, which maps the force magnitude and to the force vectors  and  : where  and  are sparse matrices containing all orientation vectors of unilateral and bilateral forces, i.e., all̂ vectors as shown in Fig. 1. More details about the definition of  can be found in the following. Both  and  depend only on the geometry norm vectorŝ ,̂ and location of constraints , , together with the particle indices , , i.e., which particles appear within the vicinity of each other and which are bound to each other by springs.
Further, this constraint formulation is also applicable to the case where one constraint is not between a pair of particles but between one particle and one externally imposed confinement or boundary, for example, a flat substrate or a spherical shell. The only necessary modification in this case is to ignore one side of the collision geometry when constructing the matrix  and  . For example, if a particle collides with a fixed substrate, we only includê and in  , because this substrate does not appear in the mobility matrix .
Appendix 4 Figure 1. The geometry for a pair of rigid particles. The distance between two marked points Φ = | |, where = + − − .

Temporal discretization and convex quadratic programming
Eqs. (45) and (46) generate a differential variational inequality (DVI), which can be solved when equipped with a timestepping scheme. In this work we use the linearized implicit Euler timestepping scheme, similar to our previous work Yan et al. (2019,2020), for three reasons: • It is straightforward to integrate with both the Brownian motion and the stochastic binding and unbinding of crosslinkers into an Euler scheme. • The scheme cannot be explicit. Otherwise Δ is limited to be tiny by the temporal stiffness of collision and doubly bound crosslinkers. • The implicit scheme is linearized to avoid expensive large-scale non-linear problems.
With timestep Δ = ℎ, eqs. (45) and (46) are discretized at timestep as: The unknowns to be solved at every timesteps are the constraint force magnitude , .
Eqs. (47)d and e are nonlinear because +1 and +1 are nonlinear functions of  +1 . Therefore, we linearize these two terms: Here we have also rewritten the Eq. (47)e into a equivalent form, similar to Eq. (47)d. The right side, ≥ 0 and ∈ ℝ should be understood in the component-wise sense. Eqs. (48)a and b are now closed and , can be solved. We shall drop the superscript in the following derivations because we shall repeat this solution process at every timestep. Then eqs. (48) can be written in the block-matrix form: where the blocks are clear from eqs. (48) Here we used the fact that: The first relation has been well known in the problem of collision constraints Anitescu et al. (1996). In this work we extend this result to bilateral constraints. A proof of this is detailed in Section Symmetry of the geometrically constrained optimization problem.
This formulation means that the coefficient matrix is Symmetric-Positive-Semi-Definite (SPSD), because the mobility matrix  is SPD and 1 ℎ  −1 is positive & diagonal: Because of this SPSD property, solving Eqs. (48) is equivalent to solving a constrained quadratic programming (CQP) due to the Karush-Kuhn-Tucker condition Nocedal and Wright (2006): subject to Here = [ , ] ∈ ℝ + is a column vector, and This can be conveniently understood as following. represent the current values of the constraint functions plus the (linearized) changes due to non-constraint forces  , such as Brownian fluctuations.
represent the linearized relation between the unknown constraint force and the changes of the constraint functions .
Solving one global optimization problem at every timestep is usually expensive, because the dimension of this problem (53) can be very large in a system with many particles and constraints. However, this CQP. (53) is a class of well understood optimization problem and fast algorithms exist. We previously developed a fully parallel Barzilai-Borwein projected gradient descent (BBPGD) method Yan et al. (2019,2020) to efficiently solve this problem for unilateral constraints only. In this work we found that the same BBPGD method also works very well for the current problem.
One way to understand the constraint optimization method is that the temporal integration 'jumps' on a timescale that the relaxation timescales of unilateral and bilateral constraints (collisions and crosslinker springs) are bypassed. As a special case, in the limit of infinitely stiff springs where  −1 → the quadratic term matrix is still SPSD and the Eq. (53) is still easily solved. Physically speaking, in this case the bilateral constraints degenerate from deformable springs to non-compliant joints.
Last but not least, due to the linearization in Eqs. (48) our geometric constraint method has some inevitable numerical errors in imposing both types of constraints for any finite timestep size Δ = ℎ. In other words, there may be some slight residual overlaps between filaments even if Eq. (48) are exactly solved. Such residual overlaps converge to zero as the timestep size Δ = ℎ decreases to zero, which follows the typical first order numerical convergence. In principal, such residual overlaps due to linearization errors can be eliminated if the full nonlinear constraint problem is solved. However, the cost for a full nonlinear solution is prohibitive. Therefore, in our implementation we do not pursue the elimination of such residual overlaps. Instead, we focus on the stability of temporal integration, i.e., the temporal integration of trajectory is stable even if very large forces suddenly appear on some particles due to, for example, Brownian noise or a large number of doubly bound crosslinkers. We have also benchmarked our algorithm such that the average physical properties of the entire suspension converge to the reference values. For example, our method accurately captures the system stress, the equation of state and isotropic-nematic phase transition of rigid Brownian spherocylinders Yan et al. (2019).

Symmetry of the geometrically constrained optimization problem
We briefly prove the symmetry of Eq. (51). The derivation in this section is applicable to rigid particles with arbitrary shapes. The configuration of each particle is tracked by its center location in the lab frame and its orientation as a unit quaternion = [ , ] ∈ ℝ 4 . For an arbitrary 3D vector which is attached to a particle and follows the particle's motion, its image in the lab frame following the particle's rotation is: where ∈ ℝ 3×3 denotes the rotation matrix generated by the unit quaternion . For both unilateral and bilateral constraints,  has a sparse column structure: where , are particle indices for the -th column. For example, for a system with 4 particles 0, 1, 2, 3 and two possible collision pairs 0, 1 and 1, 3, the  matrix for collision (unilateral) constraints is: Because of this structure, to prove Eq. (51) we only need to prove the equality = ∇   for a pair of particles , , as shown in Fig. 1.
We consider two rigid particles centered at , , each has a point fixed on the body (not necessarily on the surface). and are vectors in the lab frame from the particle centers to the points. The distance between these two points follows the rigid body motion of both particles: where and are the well-known rotation matrices. and are locations of those two points in their intrinsic coordinate systems. Φ = | | is simply the distance between the two points, dependent on the motion of the two rigid particles.
According to our definition, maps the force magnitude between the two particles to force and torque vectors on each particle: ∇   can also be explicitly written as follows: Further, we notice the symmetry of P and Q in the above equations of and ∇  , we only need to prove the following equality for P: In Eq. (62) the only difference between unilateral and bilateral constraints are how the two points on particles P and Q are picked. For unilateral (collision) constraints, the two points are where the distance Φ reaches the minimal distance between the two particles.
For bilateral constraints, there is no such restriction and the two points are arbitrary. Obviously, we only need to prove this latter case, i.e., to prove Eq. (62) when is an arbitrary vector.
The first row of Eq. (62) is straightforward because The second row can be proved as follows. We first derive some general results about quaternions and rotation matrices, dropping the subscript to simplify equations. When the particle rotates with an angular velocity , the motion of satisfieṡ = × , i.e.,̇ = = = can also be directly computed by applying the chain rule on Eq. (55), because is intrinsic to the particle invariant over time: The matrix Ψ bridges angular velocity and quaternion by definition: We have This must be valid for arbitrary , which is only possible when Now for another arbitrary vector : Using Eq. (69) we can prove the second row of Eq. (62). We first calculate the derivatives of Eq. (62) using dummy indices: Multiply the matrix Ψ on both sides: Substitute the right side by Eq. (69), we get: This is exactly the right side of Eq. (62) because by definition = and = ∕Φ. Therefore Eq. (62) holds and the equality Eq. (51) holds.

Implementation
As mentioned above, at each timestep we first update the crosslinkers and then the filaments. We implement the two steps in a fully parallelized C++ codebase, utilizing MPI and OpenMP and scalable to hundreds of CPU cores.
In the crosslinker-update step, we have assumed that every crosslinker has bindingunbinding probabilities independent of other crosslinkers. Therefore, it is straightforward to parallelize this step, we only need to search the vicinity of each crosslinker to find the candidate filaments that this crosslinker may bind to. This can be conveniently accomplished by a standard near neighbor detection operation based on bounding volume hierarchy Iwasawa et al. (2016), where the search radius is determined by the maximum stretch of each crosslinker. Once the candidate filaments for each crosslinker have been found, we compute the k-MC probabilities using a precomputed lookup table with interpolation to speed up the numerical integration while maintaining accuracy. This step is also parallel on all CPU cores.
After the positions of crosslinkers have been updated, we update the set of bilateral constraints in the constraint solver. If one crosslinker has changed its status from doubly bound to singly bound, the corresponding constraint is removed from , and vice versa. The geometric matrix  is also updated according to the current geometry, i.e., those locations , and norm vectorŝ ,̂ . Then, a near neighbor detection operation is performed for all filaments to determine the unilateral constraints and its geometry  . If two filaments are far away from each other, there is no need to include this pair in the constraint solver because it is impossible for them to collide within this time step Δ . Therefore, we include only close pairs whose minimal distance is below some threshold value . is controlled by system dynamics, i.e., how far each filament may move within each timestep. Empirically, we take to be the diameter of each filament.
Once the constraint problem Eq. 48 has been constructed, we run a fully parallel iterative Barzilai-Borwein Projected Gradient Descent (BBPGD) solver Yan et al. (2019) to solve for constraint forces and , together with the velocities  and  due to constraint forces  and  by solving the equivalent CQP 53. The cost of every BBPGD iteration scales as ( + ), i.e., the total dimension of the linear constraint problem. The number of iterations needed depends on the complexity of the structure. For example, if all filaments are far from each other such that almost no collisions or no doubly bound crosslinkers exist, the solution converges almost immediately. If all filaments are densely packed and many doubly bound crosslinkers form between the filaments, many iterations may be necessary. Empirically, the solution of Eq. 53 converges in a few hundred BBPGD iterative steps for common biological structures such as microtubule asters or bundles. However, each iteration of BBPGD is cheap because we only need to compute ∇ = + . This sparse matrix-vector multiplication spmv is a well optimized standard mathematical operation. The BBPGD solver is implemented using the Trilinos package for distributed linear algebra operations. Once  and  have been solved, the filament configuration is updated and then next timestep starts. centers at steady state, for BMT and NBMT cases. 'DBSCAN' and 'Graph' refer to two different methods of identifying aster centers, based on spatial locations of all microtubule minus ends, and the crosslinking connectivity, respectively. 500 snapshots at simulation steady state are used to compute ( ) and ( ), for each case.
To quantify the spatial aster center distribution, we identify aster centers for each snapshot of data. For cross validation, we use two different methods to identify the aster centers: 'DBSCAN' and 'Graph'. The implementation details are discussed in the following. Once aster centers are identified, we compute the radial distribution function ( ). Then, we compute structure factor ( ) based on ( ) as because the structure of aster centers is isotropic and the orientation of does not matter. Fig. 1 summarizes the results for BMT and NBMT systems. Both 'DBSCAN' and 'Graph' methods generate similar results. According to ( ), there is a clearly length scale for the NBMT case at ≈ 0.8 µm −1 . This reflects the spacing between individual asters at approximately 1.2 µm. This length scale is straightforward to understand. Since we used microtubules of length 0.5 µm, if two aster centers are smaller than 2 , then the edge of two asters may touch or overlap, and are likely to be crosslinked by kinesin-5 motors and gradually merge into one bigger aster. With this length scale argument, we can estimate the total number of asters in the simulation box to be ( box ∕(2 )) 3 ≈ 10 3 . This simple estimation agrees with our aster center identification results, which on average 1200 aster centers are found for each snapshot of steady state configuration.
The BMT case does not show such a significant special length scale in ( ), but they do show larger spacing between asters according to ( ), compared to the NBMT case. This agrees with the snapshots shown in Fig. 7 in the main text, where asters are larger but more distributed in space. Both methods identified on average 280 aster centers for each snapshot of steady state configuration.

Identify aster centers by DBSCAN method
DBSCAN stands for Density-Based Spatial Clustering of Applications with Noise and is a method to identify clusters from points in space. With a given distance and a threshold min of minimal number of points, DBSCAN searches all clusters such that each cluster has no less than min points and no point in one cluster is more than distance separated from other points in the same cluster.
To apply DBSCAN, we first create a point cloud using the location of all microtubule minus ends in the system, and then run the algorithm using the function cluster.dbscan from the python package scikit-learn. Once clusters have been identified, we compute the aster centers by averaging the location of all points in each cluster.
We set = 100 nm, because according to Fig. 7 in main text, the minus ends of microtubules are separated roughly 25 + 53nm. We also set min = 5.

Identify aster centers by Graph method
The entire microtubule-kinesin system can be abstracted as an undirected graph, where each microtubule is a node marked by their index and each doubly bound kinesin form an edge. Then, one aster is simply abstracted as a connected component of the graph. We use the connected_components() function in the python package networkx to find all such connected components, with minimal number of microtubules min = 5. We identify aster centers by computing the average location of minus ends of these connected components.

Confined filament-motor protein assemblies
This section provides more details about the simulation in Section Confined filament-motor protein assemblies of main text. We simulate 9,216 microtubules and 27,648 crosslinking motor proteins in a cylindrical volume. Microtubules are modeled as rigid spherocylinders with length 0.25 µm and diameter 25 nm (aspect ratio of 10 ). Crosslinking motor proteins are modeled as Hookean springs. The cylindrical axis is oriented along the + direction, with a periodic boundary condition. The radial direction has a hard confinement boundary. System temperature is fixed at 300 K and the simulation timestep is 10 −4 s, with the system configuration recorded every 500 steps. Solvent viscosity is set at 0.01 pNµm −2 s. Values for the cylinder diameter, ∈ {0.25, 0.75}µm are chosen to disrupt the self-assembly of an ideal aster. Initially, microtubules were aligned along the direction (cylinder axis) such that the initial nematic order parameter, = ⟨ 1 2 (3 cos 2 − 1)⟩, was 1. Here, is the angle between the microtubule orientation vector and the + axis, and ⟨.⟩ denotes an average over all microtubules. Equal numbers of microtubules are oriented in the + and the − direction such that the polar order parameter, = ⟨cos ⟩ = 0.

Structural Quantification
To measure the structure of our steady-states, we compute the local packing fraction local ( ), local nematic order, local crosslinker density, and pair distribution functions. For the first three quantities, we start by dividing the volume into cylindrical bins with their axis in the + direction. The diameter of the bins is equal to , and the height is chosen as 25 nm. The local packing fraction is computed by calculating the cumulative volume of all microtubules that fall inside each bin, and then dividing by the bin volume. For simplicity, we treat the filaments as cylinders (such that there as no hemispherical caps at their ends). For the local nematic order parameter local ( ), we find the total number of microtubules, ( ), that pass through each bin at some location . For each microtubule in the bin, we compute it's individual contribution to the local nematic order parameter, local ( ) = 1 2 (3 cos 2 − 1). We weight each local ( ) by a factor ( ) that depends on the length of microtubule that falls inside the bin, normalized by the cumulative length of all other microtubules that traverse the bin. We calculate the local nematic order parameter as Local crosslinker density, ( ), is found by counting the number of center points of crosslinking motor proteins that fall in each bin, and then dividing by the bin volume. For this calculation, we only consider doubly-bound crosslinking motor proteins. Finally, we compute the pair distribution functions by finding the distance (using the nearest image convention in ) of all microtubules from a single reference microtubule. Repeating this for all microtubules as a reference, and averaging yields a pair distribution function. The useful dimensions here are and = √ 2 + 2 . Due to the non-periodic nature in , this pair distirbution function does not decay to 1 = 0.25 µm The simulation volume is a cylinder with height 144 µm. We measure structural properties of the system over the course of the simulation. A kymograph of the local packing fraction is shown in Fig. 1A. The local nematic order (Fig. 1B) shows that the polarity-sorted bilayers (PSBs) have a maximum order parameter equal to 1 The condensation of microtubules is mediated by the crosslinking motor proteins. In Fig. 1C, we show a kymograph of the local density of the crosslinking motor proteins.
The microtubule pair distribution function at steady-state (Fig. 1D)shows that plus-ends (top plot) are distributed in a ring. The ring radius is set by the length of a single crosslinking motor protein. There is negligible density away from the ring. In contrast to asters (that contain microtubules isotropically distributed around a core), microtubule centers (bottom plot) are distributed in vertically extended regions. Separation between these regions is determined by the sum of the microtubule length and the length of the crosslinking motor protein. The presence of three regions in this pair distribution plot is evidence for a pair of layers.

Bending Rigidity
A flexible long fiber can be implemented by connecting short rigid segments into chains.
The key is how to properly implement the force and torque induced by deformation at the rigid segment joints. There are two ways to implement this, which we shall detail in the following. The first method implements the deformation of each joint with two linear Hookean springs and requires no modification to the current codebase. The second method directly incorporates the bending rigidity as a new set of constraints in the geometric constraint minimization solver, but requires some extensions to the current codebase.

Method 1: use two Hookean springs
Role spring stiffness constant free length Bending 0 Extension 0 Appendix 8 Table 1. The parameters of the two springs controlling extension and bending, respectively. The relation between 0 and 0 determine the equilibrium configuration of the two connected filaments. When 0 ≥ 0 + 2 , the straight configuration is the preferred configuration.
Appendix 8 Figure 1. The geometry of two short rigid straight fibers connected at a bending joint. The separation is exaggerated to clearly show the geometry. 1 and 2 are the orientation norm vectors of the two segments. and are the stiffness of the spring for extension and the spring for bending.
is the displacement distance from the joint rotation center. In the more detailed view of the deformed geometry of the two springs, , are the lengths of the two edges of the triangle.
Here in the first term is simply the linear extension of both springs when is small. The two-spring system generate a equivalent extensional rigidity + . The second 2 term governs the bending energy. We can tune the five parameters 0 , 0 , , , such that the connected segments reproduce the desired mechanical behavior of a flexible filament. Although the expansion Eq. 77 is general and can be fitted to many different models by tuning the five parameters, it is too complicated to be conveniently used in an actual simulation. In the following we discuss simpler special cases which are more relevant to biological filaments.
Special case 1 When model some bio-filaments such as microtubules, we sometimes assume filaments are inextensible, i.e., = ∞ and = 0 . In this special case, the energy of the two springs depends only on = 1 2 ( − 0 ) 2 . By imposing = 0 , we can solve for : = 1 2 √ 2 √ 2 cos(2 ) − 2 + 2 0 2 − 2 cos( ) Then in this case depends on 4 in the limit of → 0. To simplify the notations of the expansion, we define: The value defines three cases of the equilibrium configuration: • > 0. The equilibrium configuration of the joint is a straight line, and the bending spring is compressed at equilibrium. • = 0. The equilibrium configuration of the joint is a straight line, and the bending spring is not compressed nor stretched at equilibrium. • < 0. The equilibrium configuration of the joint is bent.
Special case 2 If we further assume that = ∞ and = 0 = 0, we have that = = 0. This means the extension spring degenerates into a point joint between the two segments. In this case the energy can be further simplified: = 1 2 − √ 2 cos ( + 0 ) + 2 2 + 2 0 + 0 2 + 2 + 0 + 2 The expansion of as → 0 is also further simplifed: Here we have the same conclusion as the previous special case, that the dependence of on can be tuned between 4 and 2 by choosing a proper value of .

Method 2: use bilateral constraints
where we have utilized the vector triple product identity: This means, to the first order of Δ only the component of rotation 1 and 2 that is inside this plane spanned by 1 , 2 affect the bending energy. Therefore, to the first order of Δ we can simplify the bending rigidity problem inside this spanned plane, although in reality the filament segments have true 3D rotations. We denote the current and next timesteps by and + 1. We have, to the first order: The rotational mobility matrix for these two rods is: where 1 and 2 are inverse of rotational drag coefficients for those two segments. The torque generated by the joint on each segment can be calculated by the derivative of bending energy . More specifically: where the scalar torque +1 is: where the higher order terms in Δ have been neglected. If the bending energy Case 2 is used, instead of +1 ∝ − +1,3 we have +1 ∝ − +1 . We can replace the expansion accordingly and the derivation remains valid.
Combining all of the above, we are effectively integrating the dynamics of all rods while ensuring Eq. 90. Skipping the timestep index , we can write the result in the same way as the bilateral Hookean spring constraints as: where = 3 ,2 and the geometric matrix  defines the direction of torque on each rod: The left side of eq. 94 means the motion of filament segments must satisfy the torquedeformation relation, while the right side means the torque can take any values. Eq. 94 is mathematically identical to the Hookean spring constraints and can be incorporated in the constraint minimization problem in the same way.
This simply means that if a straight fiber is bent to angle , its recovering motion within each timestep is proportional to the current angle . More importantly, → ∞ as the bending