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

  1. Wen Yan  Is a corresponding author
  2. Saad Ansari
  3. Adam Lamson
  4. Matthew A Glaser
  5. Robert Blackwell
  6. Meredith D Betterton
  7. Michael Shelley  Is a corresponding author
  1. Center for Computational Biology, Flatiron Institute, United States
  2. Department of Physics, University of Colorado Boulder, United States
  3. Department of Molecular, Cellular, and Developmental Biology, University of Colorado Boulder, United States
  4. Courant Institute, New York University, United States


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 Living Ensemble Simulator), a novel computational framework designed 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 tens to hundreds of thousands of cytoskeletal filaments and crosslinking motors, recapitulating emergent phenomena such as bundle formation and buckling. This simulation framework can help elucidate how motor type, thermal fluctuations, internal stresses, and confinement determine the evolution of cytoskeletal active matter.

Editor's evaluation

This article presents a new method for simulating cytoskeletal dynamics inside cells. This is an important problem in the life sciences, and the numerical methods and derived results described in the paper seem very promising to facilitate computational modelling of cell dynamics. Although the user- friendliness of the software can still be improved, the method will be of interest to a broad community of biologists and biophysicists.


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, 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 (Nédélec 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).

Computational modeling of the cytoskeleton has elucidated principles of self-organization, suggested hypotheses for experimental test, and helped interpret results of experiments (Gao et al., 2015b; Rincon et al., 2017; Bun et al., 2018; Saintillan et al., 2018; Varghese et al., 2020). Several software packages for cytoskeletal modeling are currently available, including Cytosim (Nedelec and Foethke, 2007), MEDYAN (Popov et al., 2016), AFINES (Freedman et al., 2017), and CyLaKS (Fiorenza et al., 2021). A challenge for molecular simulation is the large size of cytoskeletal systems, typically 104–107 or more filaments (Petry, 2016). While current simulations may reach O(104-105) filaments (Belmonte et al., 2017; Strübing et al., 2020), molecular modeling has required significant compromises in treating steric interactions and motor-proteins.

Here, we describe aLENS, a framework of computational methods and software designed to more efficiently and accurately simulate large cytoskeletal systems (Figure 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) (Figure 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 (Figure 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 those of experiments (Figure 5 and 7).

Figure 1 with 1 supplement see all
The computational model and demonstration of aLENS.

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

As an illustration of aLENS, Figure 1C (and movie Figure 1—video 1) 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). Although the microtubules are initially unorganized (C1), the combination of motor crosslinking and walking causes the microtubule minus-ends 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.


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 H.) 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 three 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.

1. 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 vmΔt and doubly bound motors move vFΔt along the filaments. Here, vF is the motor stepping velocity that depends on force on the motor head (Gao et al., 2015a):

(1) vF(Fproj)=vmmax(0,min(1,1+Fproj/Fstall)),

where Fproj is the projection of tether force along filament in the stepping direction. As typically found experimentally, this stepping model means that if Fproj is assisting stepping, the velocity saturates at vm; while for Fproj hindering stepping, stepping is halted when Fproj=-Fstall.

2. 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 Figure 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, that is, when vm=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 C). 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-unbound concentrations and distributions (Appendix C) (Gao et al., 2015a; Lamson et al., 2019; Allard et al., 2010), we explicitly model each crosslinker as a Hookean spring connecting two binding heads labeled as A or B. Each crosslinker has four possible states: both heads unbound (U), either A or B singly bound (SA or SB), or both heads (doubly) bound (D). For each timestep Δt, we first calculate the rates R(t) at which each head (A and B) transitions from their current state to a new binding state (i.e. for the transitions U(SA,SB)D). The transition probabilities are modeled as inhomogeneous Poisson processes with the cumulative probability function

(2) P(Δt)=1exp(0ΔtR(t)dt)=1exp(R(0)Δt+O(Δt2)).

The transitions U(SA,SB) do not stretch or compress the tether and so do not depend on tether deformation energy. However, the transitions (SA,SB)D do account for tether deformation energy (Table 1).

Table 1
The transition rates between all possible states of a crosslinker U(SA, SB)D.

(SA, SB) means either head A or B is bound but the other is unbound. All binding rates account for the linear binding density ϵ is the length of filament i with center-of-mass position xi and orientation pi inside the capture sphere with cutoff radius rc,S relative to position of motor/crosslinker x. The sum is over all possible candidate filaments i. The unbound-singly bound transition U(SA,SB) is determined by the association constant Ka and the force-independent off rate ko,S. Similarly, the singly bound-doubly bound transition (SA,SB)D is determined by the association constant Ke and force-independent off rate ko,D is the Boltzmann factor. E() in the in the (SA,SB)D transition rates refers to the tether energy of a motor E()=12κxl(f0)2. 0 is the free length of a motor, while f is the length for computing the force when attached to filaments i and j at locations si and sj: f(si,sj,xi,pi,xj,pj). The dimensionless factor λ determines the energy dependence in the unbinding rate. Both binding and unbinding rates must depend on λ and ko,d such that the equilibrium constant recovers the Boltzmann factor exp[-βE(f)] For force-dependent binding models, the E() can be simply replaced by the tether force F(). This is not used for results shown in this work, but implemented in the code.


3. 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 Δt. The first arises in models that use a stiff repulsive pairwise potential to prevent filament overlaps. For example, the Lennard-Jones potential V(σ/r)12-(σ/r)6, where r is the separation between filaments, is so steeply varying that it requires small Δt 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).

The second stability restriction arises from the fast relaxation times of crosslinking motors. When crosslinkers connect two parallel filaments, the spring tether length f relaxes according to ˙f=-λ(f-0), where 0 is the preferred length and λ=Nκxl/(4πηL/log(2L/Dfil)) (Howard, 2001). Explicit timestepping schemes require Δt<C/λ, for some constant C. For N=10 motors, tether stiffness κxl100pNµm-1, and slender body drag coefficient 4πηL/log(2L/Dfil)0.003pNsµm-1 for 1µm-long microtubules in aqueous solvent, we have 1/λ3×10-6s.

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γcol0. If one crosslinking motor connects these two filaments, its length f and force magnitude γxl satisfy the Hookean spring model γxl=-κxl(f-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 C. 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 xR3 in the lab frame and its orientation θ=[s,p]R4 as a quaternion (Delong et al., 2015). [s,p] 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 t for all N filaments can be written as a column vector with 7N entries:

(3) C(t)=[x1,θ1,,xN,θN]TR7N.

Similarly, we use the vectors U,FR6N to represent the translational & angular velocities, and forces & torques of all particles, respectively. We relate U to F via a mobility matrix MR6N×6N, dependent only upon the geometry C, and relate U to C˙(t) via a geometric matrix G:

(4) C˙(t)=GU,U=MF,

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 c refers to constraints, which includes both unilateral (with subscript u) and bilateral (with subscript b) constraints. For our problem, unilateral constraints refer to collision constraints while bilateral constraints refer to crosslinking motor constraints. The subscript nc refers to non-constraint.

For unilateral constraints, we define the grand distance vector Φu=[Φu,1,Φu,2,,Φu,Nu]TRNu, where each Φu,j is the minimum distance between a pair of filaments. Similarly, for bilateral constraints we define the grand distance vector Φb=[f,1,f,2,,f,Nb]TRNb, containing the length f,j of the doubly bound motor j. There are in total Nu possibly colliding pairs of filaments and Nb crosslinking motors. The force magnitude corresponding to these constraints are also written as vectors, γu=[γu,1,γu,2,,γu,Nu]TRNu and γb=[γb,1,γb,2,,γb,Nb]TRNb. The two types of constraints can be summarized as:

(5) 0Φu(C)γu0,K[Φb(C)Φb0]=γb.

Here, Φu and γu satisfy the complementarity (collision) constraints, while Φb and γb satisfy the Hookean spring law. Here, KRNb×Nb is a diagonal matrix consisting of all the stiffness constants, while Φb0 represents the rest length of every crosslinking motor.

Equations 4 and 5 define a differential-variational-inequality (DVI). This is solvable when closed by a geometric relation mapping the force magnitude γu and γb to the force vectors Fu and Fb:

(6) Fu=Duγu,Fb=Dbγb,

where Du and Db are sparse matrices containing the orientation norm vectors of all constraint forces (Anitescu et al., 1996; Yan et al., 2020 and Appendix D). Next, we discretize this DVI using the linearized implicit Euler timestepping scheme with Δt=h at timestep k:

(7a) 1h(Ck+1Ck)=GkUk,Uk=Mk(Fuk+Fbk+Fnck),
(7b) Fuk=Dukγuk,Fbk=Dbkγbk,
(7c) 0Φuk+1γuk0,
(7d) Kk[Φbk+1Φb0]=γbk.

The unknowns to be solved for at every timestep are the constraint (collision and motor tether) force magnitude γuk,γbk. This is a nonlinear DVI because Φuk+1, Φbk+1 are nonlinear functions of geometry Ck+1, although Ck+1 is linearly dependent on γuk and γbk. For a small timestep (h0), this nonlinearity can be linearized by Taylor expansion, for example, Φuk+1=Φuk+hCΦuGkUk. Then, this nonlinear DVI can be converted to a convex quadratic programming problem (Nocedal and Wright, 2006) (details in Appendix D):

(8a) minγf(γk)=12γk,TMkγk+qk,Tγ,
(8b) subject to [INu×Nu0]γk0.

Here, γk=[γuk,γbk]RNu+Nb is a column vector, and

(9) Mk=[Duk,TDbk,T]Mk[DukDbk]+[0001hKk,1],q=[1hΦuk+Duk,TMkFnck1h(ΦbkΦb0)+Dbk,TMkFnck].

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 (K10), the quadratic term coefficient matrix M is still symmetric-positive-semi-definite (SPSD) and the Equation 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 Equation 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 f=Mγ+q 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 al., 2014; Fürthauer et al., 2019) or estimate them based on similar motor proteins (Cross and McAinsh, 2014).

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, NA is the number of active motors and N is the total number of motors (active and inactive). The microtubule velocity increases as NA/N increases from 0 to 1 in experiments (Scharrel et al., 2014) and in our simulations. As shown in Figure 2, our simulations quantitatively reproduce the experimental data. To achieve this agreement, we set the active motor velocity to 1.0µms-1, so the sliding velocity at NA/N=1 matches experiment. Apart from this one experimentally constrained velocity, there are no fitting parameters in our simulation (further motor parameters are in Appendix B). In initial trial simulations, we found that changing the total motor number N didn’t noticeably affect the microtubule transport velocity. Therefore, for the results shown here we fixed N=100, similar to the experimental system. Since the transport trajectory is stable without stochastic noise, as shown in Figure 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 NA/N.

Directed transport velocity and displacement of microtubules driven by mixed active and inactive Kinesin-1 motors.

The total number of active and inactive motors is fixed at N=100 for all simulations. NA is the number of active motors. Left panel: comparison of microtubule velocities as a function of NA/N from aLENS simulations (blue diamonds) with from the reference experiment (orange circles) (Scharrel et al., 2014). Right panel: displacement vs. time of the transported microtubule obtained from simulation for several values of NA/N. The free walking velocity of active motors was set to 1.0µms-1 to match the experimental sliding velocity at NA/N=1. There are no other fitting parameters. All motor parameters are estimates based on experiments on Kinesin-1 (Scharrel et al., 2014) or similar motor proteins (Cross and McAinsh, 2014).

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 3,000 model microtubules with L=0.5µm. Initially the filaments are confined in a tube of diameter D=1µm, randomly initialized with their orientations along the +x (pink) and -x (white) directions, and packed at about 30% volume fraction. The simulated system is periodic along the x 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 Nm from 10 to 30. There is no accurate measurement for the XCTK2 motor in these experimental conditions. Therefore, we used experimental estimates of 46nms-1 for the walking speed of NCD motors (Furuta and Toyoshima, 2008). To approximate the experimental measurement of velocity that used line photobleaching (Fürthauer et al., 2019), we sample the local polarity and straining velocity using virtual sampling planes, as shown in the left panel of Figure 3. As in Fürthauer et al., 2019, Figure 3 shows that the straining velocity Vx is largely independent of the number of motors Nm and the local average polarity Px over the range simulated.

Sampled microtubule straining motion velocity vs local polarity in actively crosslinked microtubule network.

The left panel shows the simulation geometry and the sampling procedure. Microtubules are randomly initialized with orientations along the +x (pink) or -x (white) directions. XCTK2 motors are colored green. Nm is the number of XCTK2 motors per microtubule. We sample the local average polarity and straining velocity by inserting planes orthogonal to the x-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 Px along the x-axis for all microtubules intersecting this sampling window at a given time. We then compute the velocities, averaged over 10s (a duration chosen to match the experimental timescale), of microtubules intersecting each sampling window and moving along the +x and -x directions. V+x and V-x are computed from those two groups for each sampling window. The straining velocity is computed as Vx=V+x-V-x. Therefore, for every sampling window at each sampling timestep we have a pair of data values Px,Vx. The right three panels show the joint probability distribution of (Px,Vx) computed from 900,000 sampling planes for each simulation, for Nm=10,20,30, respectively.

Intuitively, the straining velocity Vx 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 26nms-1, close to the experimental measurement of 18.6±0.9nms-1.

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. Figure 4 shows detailed parallel efficiency measurements for one large-scale test case, similar to that in Figure 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.5 GHz. 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.

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

We achieve nearly ideal linear speed up as the number of cores increases (Figure 4). At 1536 cores, the efficiency remains at 93% and each timestep takes less than 1 s, 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 Δt that is one or two orders of magnitude larger than conventional pairwise potential methods. For the system simulated in Figure 4, aLENS can reach 1s physical time per day, using a timestep size of 1.0×10-5s.


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 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.01pNsµm-2 at room temperature, using a fixed timestep Δt=10-4s 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 (Figure 5). Motor crosslinking drives contraction of initially disordered, bundled filaments (Figure 5A and B). Aligning steric and crosslinking forces drive the system into a series of well-aligned bundles spanning several filament lengths (Figure 5C, see Figure 5—video 1 and Figure 5—video 2). The motors slide filaments parallel to each other, generating macroscopic extensile motion. Later, the extended network buckles and fractures (Figure 5C).

Figure 5 with 2 supplements see all
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µms-1. If bound, the motile head moves toward the microtubule minus-end, and detaches upon reaching it. Detailed parameters for this motor are tabulated in the Appendix 1. 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, C) The bundle at t=0s, 4s, and 7s. Microtubules are colored by their local nematic order parameter Slocal=32QijQij, with Qij=pipj-13δij, p 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 Nd, the number of microtubules averaged over when computing Slocal. (D1 and D2) The joint probability distributions Slocal and Nd for each microtubule for the entire systems at t=0.1s, when the dyneins crosslink microtubules but microtubules barely move from initial configuration, and at t=7s, when the bundle is nematic. (E) The average trajectories (solid lines) and their standard deviation (shaded area) of left-moving and right-moving microtubules. Dashed lines show linear fits to the average trajectory after t=4s, with results VRVL250nms-1. (F) The normal stresses and the weighted average Slocal¯ over time. Due to the symmetry in the y,z directions, only their average is shown σyy,zz=12(σyy+σzz). Collision stress is positive (extensile) and crosslinker stress is negative (contractile). The weighted average Slocal¯=NdiSlocali/Ndi.

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 Slocal and the number Nd of neighboring filaments crosslinked to a filament (Figure 5D). While the network contracts, the distribution of Nd does not 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 near-perfectly nematic (Slocal1), although less-ordered regions occur between aligned bundles of different orientations (Figure 5C1 and D2).

Inside the bundles, filament sliding by motors leads to transport along the local nematic director. Projecting filament trajectories onto the lab-frame x-axis, we observe left- and right-moving filaments that speed up early in the simulation, and then maintain constant average velocities at later time (t>4 s in Figure 5E), as filaments align due to steric and motor forces (Figure 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 (Figure 5F). However, the fixed simulation box size hinders the networks’ elongation, causing the bundles aligned with x-axis to buckle due to the net extensile stress (Figure 5F, see Figure 5—video 1 and Figure 5—video 2). In contrast, bundles not aligned with the x-axis are not constrained and so evolve into straight spikes. This misalignment of bundles is seen as a small net stress in the y,z-directions for t4s (Figure 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 (Figure 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 L/Dfil=10 are confined between two closely spaced concentric spherical shells at 40% volume fraction. The shell gap is ΔR=0.102µm, shorter than the filament length, with ΔR/Dfil4 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 (Figure 6A). Brownian motion is modeled at room temperature 300K and timestep Δt is set to 1×10-5s. Motors move toward minus ends of bound filaments at vm=1.0µms-1. Once they reach the minus ends, they immediately detach.

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 rin=5µm and rout=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 ±eθ at each point, where eθ is the polar basis norm vector of spherical coordinate system. The pure filament system is relaxed for 1s to resolve the overlaps in the initial configuration. Afterwards at t=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 t=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 t=0s to 1s. Trajectories are colored by time. It is clear that filaments move along the meridians. (C1-C5) Configuration and statistics at t=4s. (C1) The filaments and motors. Motors clearly concentrate in some areas. (C2) The motor number density, that is, number of motors per 1µm3. (C3) The nematic director field (shown as black bars) and the nematic order parameter S. (C4) The filament volume fraction. (C5) The divergence of polarity field p non-dimensionalized by filament length, that is, change of mean polarity per filament length. (D) The development of the correlation between motor number density n/nave and the polarity divergence field, at different times of the simulation. Clearly high n/nave are correlated with positive polarity divergence. (E1-E5) Configuration and statistics at t=0.27s 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 n/nave and the local filament volume fraction. For the polarity sorting case at t=4s 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.

Motors walk along the filaments, driving sliding of antiparallel filaments (Figure 6B). This leads to polarity-sorted regions at the north and south ‘poles’ of the sphere, meaning that the filament orientation p on average points toward the poles. Filaments with reversed initial polarity are transported to the equatorial region (Figure 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 (Figure 6C1).

Motors also accumulate in some regions according to the filament polarity (Figure 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 (Figure 6C2, C5 and G). This accumulation is illustrated by the positive correlation between motor density n and p at t=4s in Figure 6D. Furthermore, motor accumulation regions appear to show slightly lower filament volume fraction (Figure 6C2 and C4), as shown in Figure 6F. These correlations can be understood through the behavior of crosslinking motors near filament ends (Figure 6G). Once polarity sorted regions of filaments form, as the blue arrows represent, p>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 p>0, while plus-end motors accumulate in regions with p<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 (end-pausing, EP), the filament network contracts (Fig. 6E1-5) with volume fraction increases from 40% to 60% and eventually freezes at t=0.27s. 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.

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 (Figure 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 (Figure 7A). In one version of the model, we included thermal fluctuations that drive filament motion (Figure 7D and movie Figure 7—video 1), while in the other thermal fluctuations of filaments were neglected (Figure 7E and movie Figure 7—video 2). 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 (Figure 7D vs E).

Figure 7 with 2 supplements see all
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 B. (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 t=35s. 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 and E1-2 (colored by the spring force). Negative values mean the crosslink forces are contractile (attractive). (B) The radial distribution function (RDF) g(r) for the minus ends of all filaments at t=0.1s (dashed lines) and t=35s (solid lines). The first peak of g(r) at r=25nm corresponds to close contacts between filaments. The second peak of g(r) at r=78nm=25nm+53nm corresponds to the minus ends of filaments crosslinked by motors whose rest length is 53nm. 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: Π=13Trσ. The collision pressure ΠCol is positive (extensile), and the motor pressure ΠLinker 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 t=4s.

These differences are clear in the radial distribution function of filament minus ends, which are clustered by motors paused at filament ends (Figure 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 r=25nm=Dfil and r=78nm=0+Dfil which correspond to scale on which filaments bind to or are crosslinked by motors, respectively (Figure 7B, D2 and 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 1s. 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 (Figure 7D1). These observations and estimated lengthscale are quantitatively confirmed by analyzing the static structure factor of aster centers (details in Appendix F), which shows that the athermal simulation has approximately three times more asters than the thermal case (Figure 7D vs E).

The differences in the dynamics of aster formation are also reflected in stress measurements (Figure 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 ΠLinker initially increases quickly, reaching a peak at roughly t=4s5s, similar to the behavior during bundle contraction shown above (Figure 5F), before declining. The average time required for motors to walk to filament ends, τwalk=L/vm5s, 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. Although 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 because crosslinking motors pull filament ends together. A confining cylindrical boundary strongly modifies the conformation of these aggregated structures (Figure 8). These simulations used 0.25µm long filaments at a fixed packing fraction (ϕ=0.16), confined in two cylinders with diameters Dcyl=0.25µm and 0.75µm.

Figure 8 with 2 supplements see all
Results for the confined filament-motor protein assembly simulations with 9,216 filaments modeling microtubules and 27,648 motors modeling crosslinking kinesin-like proteins at a cylinder diameter of Dcyl=0.25µm and 0.75µm.

Initially, 0.25µm long filaments are uniformly distributed and aligned along the x-axis, with equal numbers oriented in the +x and -x directions. Crosslinking motor proteins are initially unbound and distributed uniformly as well. A and B: Snapshots of the simulation with Dcyl=0.25µm and 0.75µm at t=58s and t=120s. A1 and B1: All 9,216 simulated filaments. In A1, the the cylinder is too long to be displayed contiguously, therefore a stacked representation is shown. The filaments are colored by the value of cosθ where θ is the angle between the filament direction vector p (oriented from the minus-end to the plus-end) and the positive x-axis (pointing to the right). A2 and B2: Zoomed-in view of the filaments in the boxed regions in A1 and B1. A3 and B3: Doubly-bound motors in the boxed regions, colored by their binding force. Negative values represent contractile force while positive values indicate extensile force. A4 and B4: The local packing fraction (red line) and the local nematic order parameter (blue line), Slocalx(x)=iN(x)Wi(x)Slocalx(x)i where a filament i contributes Slocalx(x)i=12(3cos2θi-1) to the local order at x. Filament contributions are weighted by Wi(x) and summed over all filaments at x. Line plots represent an average over 1s for the snapshots in A2 and B2. Detailed parameters and calculations for the crosslinking motor proteins are presented in Appendix G.

For a small-diameter cylinder where one filament length can fit across the cylinder (Dcyl/L=1, Dcyl=0.25µm), the cylinder is too narrow for asters to form. Instead, motor sliding and end-pausing drive the filaments into polarity-sorted bilayers (PSBs, Figure 8A and movie Figure 8—video 1). A single polarity-sorted bilayer contains a central interface of highly crosslinked filament minus-ends between two antiparallel polar layers of filaments (Figure 8). 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 Slocalx(x) nearly reaches 1 within PSBs. Even the the vapor phase is close to nematic Slocalx0.6 (Figure 8A4), due to the strong confinement effect.

Next we increased the diameter of the cylinder to Dcyl/L=3 (Dcyl=0.75µm) to weaken the confinement (Figure 8B and movie Figure 8—video 2). Here, the polarity-sorted bilayers are not present, because the larger cylinder diameter allows filaments to reorient and organize into bottle-brush-like 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 (Figure 8B2). Motors become highly concentrated along the line defects at the center of the cylinder (Figure 8B3). The radial hedgehog structure of BBs is evidenced by a negative local nematic order parameter (Figure 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 (Figure 8B4, red line).


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 force-dependent 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 al., 2019). Long and flexible polymers can be simulated by chaining short and rigid segments together with flexible connections (Appendix H), even with nonlocal interactions mediated by hydrodynamics, electrostatics, or other fields (Shelley, 2016; Nazockdast et al., 2017; Maxian et al., 2021). For example, in ongoing work we have used aLENS to simulate chromatin in the nucleus as a bead-spring chain moving through the nucleoplasmic fluid, and confined by the nuclear envelope.

Recent years have seen considerable innovation in computational approaches to cytoskeletal modeling, implemented in powerful simulation packages including Cytosim (Nedelec and Foethke, 2007), MEDYAN (Popov et al., 2016), and AFINES (Freedman et al., 2017). These packages utilize a variety of coarse-grained representations of cytoskeletal elements and numerical simulation schemes, with the diversity of approaches in part reflecting the diversity of cytoskeletal systems and phenomena of interest. aLENS brings a powerful set of new capabilities to the table, significantly expanding the range of accessible time and length scales in simulations of systems in which excluded volume and crosslink-mediated interactions play an important role.

aLENS has been open-sourced on GitHub: (copy archived at swh:1:rev:f2dd484f82443735562ad7b480fe7ed9fc020fb0; Adam, 2022) and precompiled binary executable is available on DockerHub: Our GitHub documentation provides a clear roadmap for developing additional user-specific modules.

Appendix 1

Crosslinker and motor properties

Appendix 1—table 1
The parameters of the two springs controlling extension and bending, respectively.
End-pausingTrue or FalseND
One head fixedTrue or FalseND
λenergy factorND
λP,APparallel to anti-parallel factorND
0free lengthµm
rccapture radiusµm
κHookean spring constantpNµm-1
Fstallstall forcepN
dUunbound diffusivityµm2s-1
ϵbinding site densityµm-1
vmmax walking velocityµms-1
Kaassociation constant (US)(μmol/L)-1
ko,Soff-rate constant (US)s-1
Keeffective association constant (SD)ND
ko,Dforce-independent off-rate constant (SD)s-1
dSsingly bound head diffusivityµm2s-1
dDdoubly bound head diffusivityµm2s-1
vSsingly bound walking velocityµms-1
xcforce-dependent unbinding lengthµm
Appendix 1—table 2
Properties of crosslinkers used in the main text.

ND means dimensionless. Parameters given as an array [a,b] means the two values are used for each each of a crosslinker, respectively. Kinesin-5 parameters are adapted from Blackwell et al., 2017. Dynein parameters are adapted from Foster et al., 2017. Kinesin-1 parameters are adapted from Scharrel et al., 2014.

ParameterKinesin-5DyneinKinesin-1Inactivated Kinesin-1
One head fixedFalseTrueTrueTrue

Appendix 2

Crosslinker binding and unbinding

Kinetic Monte-Carlo: crosslinking protein-filament interactions

Appendix 2—figure 1
Labels and definition of kinetic rates for crosslinking proteins binding to filaments (green) implemented in the kinetic Monte Carlo algorithm.

Crosslinking proteins (blue) exist in three different states: neither head attached to a filament (unbound), bound with one head attached to a filament (singly bound), and crosslinking two filaments (doubly bound). Motors and crosslinkers may have different rates for separate binding heads (A,B).

Our molecular model simulates distinct filaments and crosslinking proteins (crosslinking motor proteins, passive crosslinkers, etc.). This model includes fluctuations in bound protein number and binding kinetics that recovers the equilibrium distribution of static crosslinking proteins Gao et al., 2015a; Blackwell et al., 2017; Rincon et al., 2017; Lamson et al., 2019; Edelmaier et al., 2020. Modeled crosslinking proteins in solution bind to one filament and then crosslink two filaments (Appendix 2—figure 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 C.2). All binding and unbinding rate calculations are summarized in Appendix 2—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

(10) H+B=HB,

where H is a head and B is a binding site on a filament. The association constant Ka of the heads to binding sites is described by the equilibrium equation

(11) Ka=[HB][H][B]=kon,Skoff,S,

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 kon,S and koff,S.

Unbound crosslinkers are modeled as diffusing points with center of mass positions xo(t) and diffusion constants du. The heads of a crosslinker have spatial- and time-dependent concentrations [H]=c(x,t). The head binding rate is the volume integral over the product of the on-rate constant, binding site density, and crosslinker concentration

(12) Ron,SA(t)=kon,SAiLidsdx3ϵδ3(xxi(s))c(x,t),

where kon,SA=KaAkoff,SA and ϵ is the linear binding site density along filaments. The lab position along the ith filament xi(s) is parameterized by s.

The binding probability in a timestep Δt is an in-homogeneous Poisson process with the cumulative probability function

(13) Pon,S(Δt)=1exp(0ΔtdtRon,S(t)).

We assume the tight binding limit kon,Skoff,S and do not consider multiple binding and unbinding events of one crosslinker during a timestep Δt. The average number of multiple events may be calculated from binding parameters and the timestep allowing one to set a probability threshold (Lamson et al., 2021). Heads of the same crosslinking protein are forbidden to be bound to the same filament at the same time.

To describe c(x,t) during a timestep, we first consider a crosslinking protein with two heads connected by a flexible but relatively stiff polymer tether with length o. Because of the tether’s stiffness, the radius of gyration of an unbound protein is assumed to be rg=o/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, that is, 6dUΔt<rg, rotational diffusion dominates and we approximate heads to be within a sphere of radius rc,S=rg centered at xo. Realistically, the head distributions can vary within this volume but because o is small compared to filament lengths, we approximate the head distribution as being uniform, that is, c(x,t)=(4πrc,S3/3)1(1Θ(|x|rc,S)), where Θ(x) 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

(14) Ron,S(t)=3ϵkon,S4π(rc,S)3iLin,i(t),

where Lin is the filament i’s length segment within rc,S. To account for cylindrical filaments with diameter Dfil, we augment the binding radius such that rc,S=rg+Dfil/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 Ron,S(t)Ron,S(ti), for t(ti,ti+Δt).

However, it is uncommon that an unbound crosslinker or motor will diffuse less than rg 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

(15) co(x,t)t=dU2co(x,t),

which has the solution

(16) co(x,t)=1(4πdU)3/2exp[|xxo|24dUt].

If the characteristic diffusion length dUΔtrg, then equation (14) underestimates binding (Appendix 2—figure 2). The large diffusion distance also allows us to approximate the head distribution as matching the protein’s spatial distribution, co(x,t)c(x,t). Substituting the binding rate equation (12) and the solution to the diffusion equation (16) into the integral of equation (13) gives

(17) 0ΔtdtRon,S(t)=i0Δtdtkon,Sϵ(4πdUt)3/2Lidsidx3δ3(xxi(si))exp[|xxo|24dUt].

For straight, rigid filaments, we take the volume and time integrals while reparameterizing |xxo|2 by the crosslinker’s perpendicular h and parallel s distances from a filament segment’s center. This gives the linear binding probability density for a filament

(18) pon,S(h,si,Δt)=0ΔtdtRon,Ssi=Kaϵko,S4πdU(1hi2+si2erfc[hi2+si24dUΔt]).

Integrating over si gives the binding probability of one crosslinker head to a single filament. The total binding probability is then

(19) Pon,S(Δt)=1exp(iNLidsipon,S(hi,si,Δt)),

where N is the number of filaments surrounding the crosslinker head.

Appendix 2—figure 2
Comparison of initially unbound passive crosslinkers binding to a 1µm filament with binding radii set to a crosslinker’s radius of gyration versus a binding radius dUΔt.

(A, D) Number of singly bound crosslinkers over time as the unbound diffusion constant dU (A) and time step Δt (D) vary while binding radius remains unchanged rc,S=(o+Dfil)/2. Red lines mark the steady-state number of singly bound crosslinkers for a homogeneous reservoir calculated from equations (26)-(30). (B, E) Same as A and D but binding radius scales as the root mean square of diffused distance in a time step rc,S=6dUΔt. (C, F) Comparison of the steady-state number of singly bound crosslinkers as a function of dU (C) and Δt (F) for both definitions of rc,S. Simulation parameters: periodic box length =2µm, filament length L=1µm, linear binding site density ϵ=27µm-1, crosslinker number N=4000, crosslinker length o=50nm, association constant Ka=90.9(μmol/L)-1, unbinding rate ko,S=5s1. Unless otherwise stated unbound diffusion constant dU=1µm2/s and timestep Δt=0.0001s

Calculating the binding probability from this function and ensuring that the protein unbinds so that detailed-balance is satisfied is computationally prohibitive. Instead, setting rc,S to the root mean square diffusion distant 6dUΔt 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 (Appendix 2—figure 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

(20) koff,S=ko,S.

If implementing equation 14, the protein unbinds into a uniform sphere of radius rc,S. 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 i at position si, the binding rate constant kon,[,D] for binding to a location sj on filament j is

(21) kon,D(si,sj)=Keko,DeβUi,j(si,sj),

where ko,[,D]=koff,[,D](Ui,j=0) is the unbinding-rate of crosslinking proteins when no force is applied, Ke is a binding association constant similar to Ka, and Ui,j is the free energy contribution from the tether

(22) Ui,j=κxl2((si,sj)-o-Dfil)2

Before crosslinking, the unbound motor head explores an effective volume Vbind 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.

(23) Vbind=eβUi,jdr3=4π0rc,DeβUi,jr2dr.

We impose an integration cutoff radius rc,D 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) Lamson et al., 2021.

Theory and experimental evidence suggests that binding rates depend not only on energy but also force Evans and Ritchie, 1997; Dudko et al., 2006; Walcott, 2008; Guo et al., 2019. This allows for catch-bond-like behavior where proteins remain crosslinked for longer if under tension and release quicker if compressed. We replicate this behavior with the function

(24) fF(si,sj)=κxl(λ2((si,sj)-o)2+xc((si,sj)-o))

where λ and xc are the energy factor and characteristic length specifying the behavior of the energy- and force-dependent binding/unbinding, respectively. For values of xc<0, you see catch-bond like behavior whereas values of xc>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

(25) koff,D(si,sj)=ko,DeβfF(si,sj).

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.

Appendix 2—table 1
The transition rates between all possible states of a crosslinker U(SA, SB)D.

(SA, SB) means either head A or B is bound but the other is unbound. All binding rates account for the linear binding density ϵ is the length of filament i with center-of-mass position xi and orientation pi inside the capture sphere with cutoff radius rc,s relative to position of motor/crosslinker x. The sum is over all possible candidate filaments i. The unbound-singly bound transition U(SA,SB) is determined by the association constant Ka and the force-independent off rate ko,s. Similarly, the singly bound-doubly bound transition (SA,SB)D is determined by the association constant Ke and force-independent off rate ko,d with an additional factor Vbind, the effective volume explored by the unattached head while the motor/crosslinker is singly bound. Energy dependence in the (SA,SB)D transition rates is imposed by the Boltzmann factor that is a function of β=1/(kBT), the tether length of the motor/crosslinker attached to filaments i and j at locations si and sj(si,sj,xi,pi,xj,pj), the characteristic length of the tether not under load o, and the filament diameter Dfil. The dimensionless factor λ determines the energy dependence in the unbinding rate while the xc 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 N(N-1) crosslinking densities ψi,jA,B, 2N singly bound densities χiA and χiB, and an unbound density C to model all crosslinking proteins between N filaments. By convention,ψi,jA,B=ψj,iB,A

(26) ψi,jA,Bt+si[di,jAψi,jA,Bsi+(vdrag,i,jA+vwalk,i,jA)ψi,jA,B]+sj[di,jBψi,jA,Bsj+(vdrag,i,jB+vwalk,i,jB)ψi,jA,B]=ϵ(kon,i,jAχiB+kon,i,jBχjA)(koff,i,jA+koff,i,jB)ψi,jA,B,
(27) ψi,jB,At+si[di,jBψi,jB,Asi+(vdrag,i,jB+vwalk,i,jB)ψi,jB,A]+sj[di,jAψi,jB,Asj+(vdrag,i,jA+vwalk,i,jA)ψi,jB,A]=ϵ(kon,i,jBχiA+kon,i,jAχjB)(koff,i,jA+koff,i,jB)ψi,jB,A,
(28) χiAt+si[diAχiAsi+vwalk,iAχiA]=ϵkon,SACkoff,sAχiA+jLjdSj(koff,i,jBψi,jA,Bϵkon,i,jBχiA),
(29) χiBt+si[diBχiBsi+vwalk,iBχiB]=ϵkon,SBCkoff,sBχiB+jLjdSj(koff,i,jAψi,jB,Aϵkon,i,jAχiB),
(30) Ct=iLidsi[koff,SAχiA+koff,SBVϵ(kon,SA+kon,SB)C].

For heads A and B, crosslinking diffusion constants di,jA and di,jB, the drag speeds vdrag,i,jA and vdrag,i,jB, and walking speeds vwalk,i,jA and vwalk,i,jB have been shown to depend on the force exerted on the binding heads and thus the stretch of the tether (si,sj). No tether force acts on singly bound motor proteins but the singly bound diffusion constants diA and viB and walking speeds vdrag,iA and vdrag,iB may depend on si through some other physical mechanism such as crowding or state of the filament’s lattice. The total number of crosslinking proteins of the system is

(31) N=ijLidsiLjdsjψi,jA,B+iLidsi(χiA+χiB)+CV

and is constant in time.

Appendix 3

Filament dynamics

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 x in the lab frame and its orientation θ=[s,p]R4 as a quaternion Delong et al., 2015. This p 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 C for all N filaments can be written as a column vector:

(32) C=[x1,θ1,,xN,θN]TR7N,

which is a function of time: C(t). The translational and angular velocity U,Ω of all filaments can also be written as a column vector:

(33) U=[U1,Ω1,,UN,ΩN]TR6N.

Similarly we can write the force and torque F,T applied on all filaments as a column vector:

(34) F=[F1,T1,,FN,TN]TR6N.

The kinematic equation of motion 35 maps U to C˙(t)=C/t, via a geometric matrix G.

(35) C˙(t)=GU.

GR7N×6N is a block diagonal matrix, with one 3×3 and one 4×3 block for each particle:

(36) G=[I3Ψ1I3Ψ2].

I3 is the 3×3 identity matrix, same for every particle. Each I3 block simply corresponds to the translational motion x˙j=Uj of each particle j. Each ΨjR4×3 refers to the rotational motion θ˙j=ΨjΩj of each particle, where for each j:

(37) Ψ(θ)=12[pTsIP],Pij=ϵikjpk.

Here ϵikj 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 M maps the force F linearly to the velocity U:

(38) U=MF,F=FC+FL+FB+FE.

F includes collision force FC between particle-particle and particle-container pairs, linker force between particle pairs FL generated by doubly bound crosslinkers, Brownian force on each particle FB generated by thermal fluctuations, and other externally applied forces FE through gravity and electrostatic fields.

In principal, Equation 35 together with Equation 38 can be integrated directly because both M and F are functions of the geometry C and time only. However, this approach is usually impractical, because FC or FL is usually very stiff functions of the geometry. For example, the collision force FC 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 FC in our previous work on Brownian spherocylinders Yan et al., 2019 and rigid spheres in Stokes flow Yan et al., 2020. Instead of computing FC using repulsive potentials, we imposed non-overlapping constraints on the geometry C while integrating Equation 35.

Equation of motion with geometric constraints

In the following, the subscript c refers to constraints, which includes both unilateral (with subscript u) and bilateral (with subscript b) 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 nc refers to non-constraint, i.e., physical components that are independent of the constraints.

For unilateral constraints, we define the grand distance function Φu between every pair of particles as a column vector:

(39) Φu=[Φu,P1Q1,Φu,P2Q2,,Φu,PNuQNu]TNu,

where each Φu,PjQj is the minimal distance between particles with indices Pj and Qj. Similarly, we define the grand distance function Φb for bilateral constraints:

(40) Φb=[Φb,P1Q1,Φb,P2Q2,,Φb,PNbQNb]TNb,

where each Φb,PjQj is the distance between two fixed points on particles Pj and Qj, respectively. Physically, Φb,Pj,Qj is simply the length of each doubly bound crosslinker. With this definition, there are in total Nu unilateral and Nb bilateral constraints in the system. In other words, there are in total Nu possibly colliding pairs of filaments and Nb doubly bound crosslinkers. Both kinds of constraints are functions of the system geometry, so we shall write them as Φb(C) and Φu(C) 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:

(41) γu=[γu,1,γu,2,,γu,Nu]TNu,
(42) γb=[γb,1,γb,2,,γb,Nb]TNb.

For each Φu,PjQj or Φb,PjQj, there is a corresponding force magnitude γu,j or γb,j, the (normalized) direction vector e^Pj=-e^Qj of this force, and the location yPj and yQj where this force is applied on the filament Pj and Qj respectively, as shown in Appendix 3—figure 1. With norm vectors defined in this way, γu or γb is positive when the force is repulsive between two filaments.

For unilateral constraints Φu and γu satisfy this complementarity condition:

(43) 0Φu(C)γu0

This condition means Φu(C) and γu are orthogonal to each other, and all components of Φu(C) and γu are non-negative Yan et al., 2019.

For bilateral constraints Φb and γb satisfy this linear equality condition because they are modeled as Hookean springs:

(44) K[Φb(C)-Φb0]=-γb.

KNb×Nb is a diagonal matrix, with the stiffness constant κ for each spring on its diagonal [κ1,κ2,]. Obviously every constant κj is positive. Φb(C) and Φb0 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:

(45a) C˙(t)=G(C)U,
(45b) U=MF=M(Fu+Fb+Fnc),
(45c) 0Φu(C)γu0,
(45d) K[Φb(C)-Φb0]=-γb.

These equations are solvable when closed by a geometric relation, which maps the force magnitude γu and γb to the force vectors Fu and Fb:

(46) Fu=Duγu,Fb=Dbγb,

where Du and Db are sparse matrices containing all orientation vectors of unilateral and bilateral forces, i.e., all e^ vectors as shown in Appendix 3—figure 1. More details about the definition of D can be found in the following.

Both Du and Db depend only on the geometry norm vectors ePj^,e^Qj and location of constraints yPj,yQj, together with the particle indices Pj,Qj, 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 Du and Db. For example, if a particle P collides with a fixed substrate, we only include eP^ and yP in Du, because this substrate does not appear in the mobility matrix M.

Appendix 3—figure 1
The geometry for a pair of rigid particles.

The distance between two marked points Φ=|r|, where r=xP+yP-xQ-yQ.

Temporal discretization and convex quadratic programming

Equations 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., 2020; Yan et al., 2019, 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 Δt 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 Δt=h, Equations 45 and 46 are discretized at timestep k as:

(47a) 1h(Ck+1-Ck)=GkUk,
(47b) Uk=Mk(Fuk+Fbk+Fnck),
(47c) Fuk=Dukγuk,Fbk=Dbkγbk,
(47d) 0Φuk+1γuk0,
(47e) Kk[Φbk+1-Φb0,k]=-γbk.

The unknowns to be solved at every timesteps are the constraint force magnitude γuk,γbk. Equations 47d and e are nonlinear because Φuk+1 and Φbk+1 are nonlinear functions of Ck+1. Therefore, we linearize these two terms:

(48a) 0Φuk+hCΦukGkMk[Fnck+Dukγuk+Dbkγbk]γuk0,
(48b) 0=ΦbkΦb0+hCΦbkGkMk[Fnck+Dukγuk+Dbkγbk]+K1γbkγbkR.

Here we have also rewritten the Equation 47 e into a equivalent form, similar to Equation 47 d. The right side, γu0 and γb should be understood in the component-wise sense. Equations 48 a and b are now closed and γuk,γbk can be solved. We shall drop the superscript k in the following derivations because we shall repeat this solution process at every timestep.

Then equation (48) can be written in the block-matrix form:

(49) 00=[AD]+[BCEF][γuγb][γuγb]0

where the blocks are clear from equation (48)

(50a) A=1hΦu+DuTMFnc
(50b) B=CΦuGkMDu=DuTMDu
(50c) C=CΦuGMDb=DuTMDb
(50d) D=1h(Φb-Φb0)+DbTMFnc
(50e) E=CΦbGMDu=DbTMDu
(50f) F=CΦGMDb=DbTMDb+1hK-1

Here we used the fact that:

(51a) CΦuG=DuT,
(51b) CΦbG=DbT.

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 M is SPD and 1hK-1 is positive & diagonal:

(52) [BCEF]=[DuTDbT]M[DuDb]+[0001hK-1]

Because of this SPSD property, solving Equations 48 is equivalent to solving a constrained quadratic programming (CQP) due to the Karush-Kuhn-Tucker condition Nocedal and Wright, 2006:

(53a) minγf(γ)=12γTMγ+γTq,
(53b) subject to [INu×Nu0]γ0.

Here γ=[γu,γb]Nu+Nb is a column vector, and

(54) M=[BCEF],q=[AD].

This can be conveniently understood as following. q represent the current values of the constraint functions Φ plus the (linearized) changes due to non-constraint forces F, such as Brownian fluctuations. M 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., 2020; Yan et al., 2019 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 K-10 the quadratic term matrix M is still SPSD and the Equation 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 Equations 48 our geometric constraint method has some inevitable numerical errors in imposing both types of constraints for any finite timestep size Δt=h. In other words, there may be some slight residual overlaps between filaments even if Equation 48 are exactly solved. Such residual overlaps converge to zero as the timestep size Δt=h 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 Equation 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 x in the lab frame and its orientation as a unit quaternion θ=[s,p]4. For an arbitrary 3D vector Y which is attached to a particle and follows the particle’s motion, its image y in the lab frame following the particle’s rotation is:

(55) y=RY.

where R3×3 denotes the rotation matrix generated by the unit quaternion θ.

For both unilateral and bilateral constraints, D has a sparse column structure:

(56) D=[DP1Q1,DQ2Q2,],

where Pi,Qi are particle indices for the i-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 Du matrix for collision (unilateral) constraints is:

(57) D=[D0,1,D1,3,].

Because of this structure, to prove Equation 51 we only need to prove the equality DPQ=CΦPQG for a pair of particles P,Q, as shown in Appendix 3—figure 1.

We consider two rigid particles centered at xP,xQ, each has a point fixed on the body (not necessarily on the surface). yP and yQ 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:

(58) r=xP+yP-xQ-yQ=xP+RPYP-(xQ+RQYQ),

where RP and RQ are the well-known rotation matrices. YP and YQ are locations of those two points in their intrinsic coordinate systems. ΦPQ=|r| is simply the distance between the two points, dependent on the motion of the two rigid particles.

According to our definition, DPQ maps the force magnitude γ between the two particles to force and torque vectors on each particle:

(59) DPQ=[eP,yP×eP,eQ,yQ×eQ]T,
(60) eP=r/|r|=-eQ

(CΦ)G can also be explicitly written as follows:

(61) CΦPQG=[Φ/xPΦ/θPΦ/xQΦ/θQ][I3ψPI3ψQ]

Further, we notice the symmetry of P and Q in the above equations of DPQ and CΦPQG, we only need to prove the following equality for P:

(62) [Φ/xPΦ/θP][I3ψP]=[ePyP×eP]

In Equation 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 Equation 62 when YP is an arbitrary vector.

The first row of Equation 62 is straightforward because

(63) Φ/xP=|r|/xP=r/|r|=eP

The second row can be proved as follows. We first derive some general results about quaternions and rotation matrices, dropping the subscript P to simplify equations. When the particle rotates with an angular velocity ω, the motion of y satisfies

(64) y˙=ω×y,i.e.,y˙i=yit=ϵiαβωαyβ=ϵiαβωαRβγYγ

y˙ can also be directly computed by applying the chain rule on Equation 55, because Y is intrinsic to the particle invariant over time:

(65) y˙i=Rijθkθ˙kYj

The matrix Ψ bridges angular velocity and quaternion by definition:

(66) θ˙k=Ψklωl,Ψkl4×3

We have

(67) ωlRijθkΨklYj=ϵiαβRβγYγωα

This must be valid for arbitrary ω, which is only possible when

(68) RijθkΨklYj=ϵilβRβγYγ

Now for another arbitrary vector r:

(69) riRijθkΨklYj=riϵilβRβγYγ=ϵlβiRβγYγri=[(RY)×r]l

Using Equation 69 we can prove the second row of Equation 62. We first calculate the derivatives of Equation 62 using dummy indices:

(70) Φθk=Φririθk=1Φririθk=1ΦriRijθkYj

Multiply the matrix Ψkl on both sides:

(71) ΦθkΨkl=1ΦriRijθkYjΨkl

Substitute the right side by Equation 69, we get:

(72) ΦθkΨkl=1ΦϵlβiRβγYγrl=1Φ[(RY)×r]l.

This is exactly the right side of Equation 62 because by definition yP=RPYP and eP=r/Φ. Therefore Equation 62 holds and the equality Equation 51 holds.


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 binding-unbinding 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 Φb in the constraint solver. If one crosslinker has changed its status from doubly bound to singly bound, the corresponding constraint is removed from Φb, and vice versa. The geometric matrix Db is also updated according to the current geometry, that is, those locations yP,yQ and norm vectors e^P,e^Q. Then, a near neighbor detection operation is performed for all filaments to determine the unilateral constraints Φu and its geometry Du. 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 Δt. Therefore, we include only close pairs whose minimal distance is below some threshold value δc. δc is controlled by system dynamics, that is, how far each filament may move within each timestep. Empirically, we take δc to be the diameter of each filament.

Once the constraint problem Equation 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 γb and γu, together with the velocities Ub and Uc due to constraint forces Fb and Fu by solving the equivalent CQP 53. The cost of every BBPGD iteration scales as O(Nu+Nb), that is, 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 Equation 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 f=Mγ+q. 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 Ub and Uc have been solved, the filament configuration is updated and then next timestep starts.

Appendix 4

Performance measurements

The bundle contraction-buckling simulation runs on 2 nodes connected by Infiniband, and each node has two AMD EPYC 7742 64-Core CPUs 2.25 GHz. Appendix 4—figure 1 shows the performance of the solver. Different from the aster formation case shown in Appendix 4—figure 2, computational time spent on crosslinkers is negligible. This is because as the fixed head of each dynein is permanently attached to the microtubule, we only need to update the status of the free head. Also, the free heads only experience the SD transition, which further reduces the computational cost. On the other hand, the collisions in this case is more difficult to resolve compared to the aster cases, because in nematic bundles more collisions happen and collisions may happen anywhere along the microtubule instead of only at the center of each aster. Similar to the aster case, we see that computational time for constraint solution is proportional to the number of BBPGD steps.

Appendix 4—figure 1
Performance of aLENS for the buckling simulation shown in Figure 5 of main text.

The left panel shows the wall clock time that every timestep takes. The right panel shows the number of BBPGD steps to solve the constraint optimization problem at every timestep.

For the aster formation in bulk problem, each case runs on 1 node of dual Intel Xeon 14-core CPUs E5-2680 v4 2.40 GHz. Appendix 4—figure 2 shows the performance of the solver for simulations with and without thermal fluctuations. Updating the binding states of kinesin-5 motors requires roughly the same wall clock time per timestep for the entire simulation. However, the time required to solve the constraint problem grow quickly in the initial stage. The solver cost increases mostly due to the increased number of BBPGD steps (as shown in the right panels of Appendix 4—figure 2) even though the dimension of the constraint problem Equation 53 grows as more kinesin-5 motors become doubly bound and more collisions occur as the asters form. The increase in BBPGD steps dominates because while the dimension of γ increases, the dimension of M remains constant since the number of microtubules does not change and the cost of each BBPGD step mainly depends on the cost of applying M when calculating f in solving Equation 53.

Appendix 4—figure 2
Performance of aLENS for aster formation simulations shown in Figure 7 of main text.

The left panels show the wall clock time that every timestep takes to simulate the Brownian and Non-Brownian cases. The right panels show the number of BBPGD steps to solve the constraint optimization problem at every timestep for those two cases.

Appendix 5

Aster center analysis of asters formation in bulk

This section provides more details about the simulation in Section Confined filament-motor protein assemblies of main text.

Appendix 5—figure 1
The radial distribution function g(r) and structure factor S(q) of identified aster 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 g(r) and S(q), 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 g(r). Then, we compute structure factor S(q) based on g(r) as

(73) S(q)=1+4πρ1qrsinqr[g(r)-1]dr,

because the structure of aster centers is isotropic and the orientation of q does not matter.

Appendix 5—figure 1 summarizes the results for BMT and NBMT systems. Both ‘DBSCAN’ and ‘Graph’ methods generate similar results. According to S(q), there is a clearly length scale for the NBMT case at q0.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 2L, 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 (Lbox/(2L))3103. 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 S(q), but they do show larger spacing between asters according to g(r), compared to the NBMT case. This agrees with the snapshots shown in Figure 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 Nmin of minimal number of points, DBSCAN searches all clusters such that each cluster has no less than Nmin 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 ϵ=100nm, because according to Figure 7 in main text, the minus ends of microtubules are separated roughly 25+53nm. We also set Nmin=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 Nmin=5. We identify aster centers by computing the average location of minus ends of these connected components.

Appendix 6

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 25nm (aspect ratio of 10). Crosslinking motor proteins are modeled as Hookean springs. The cylinderical axis is oriented along the +x direction, with a periodic boundary condition. The radial direction has a hard confinement boundary. System temperature is fixed at 300K and the simulation timestep is 10-4s, with the system configuration recorded every 500 steps. Solvent viscosity is set at 0.01pNµm-2s. Values for the cylinder diameter, Dcyl{0.25,0.75}µm are chosen to disrupt the self-assembly of an ideal aster. Initially, microtubules were aligned along the x direction (cylinder axis) such that the initial nematic order parameter, S=12(3cos2θ-1), was 1. Here, θ is the angle between the microtubule orientation vector and the +x axis, and . denotes an average over all microtubules. Equal numbers of microtubules are oriented in the +x and the -x direction such that the polar order parameter, P=cosθ=0.

Structural quantification

To measure the structure of our steady-states, we compute the local packing fraction ϕlocal(x), 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 +x direction. The diameter of the bins is equal to Dcyl, and the height is chosen as 25nm. 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 Slocalx(x), we find the total number of microtubules, N(x), that pass through each bin at some location x. For each microtubule i in the bin, we compute it’s individual contribution to the local nematic order parameter, Slocalx(x)i=12(3cos2θi-1). We weight each Slocalx(x)i by a factor Wi(x) that depends on the length of microtubule i 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, C(x), 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 x) 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 x and ρ=y2+z2. Due to the non-periodic nature in ρ, this pair distirbution function does not decay to 1.

Appendix 6—figure 1
Results for the confined microtubule-motor protein assembly simulations with Dcyl=0.25µm.

(A) A kymograph of the local microtubule packing fraction ϕlocal(x). Initially, crosslinking motor proteins drive contraction of the system into condensed regions that break into PSBs over time. (B) A kymograph of the local nematic order parameter Slocalx(x). (C) A kymograph of the density of the crosslinking motor proteins, C(x). Condensation of microtubules coincides with condensation of the crosslinking motor proteins. (D) Pair distribution function for microtubule plus-ends (top) and microtubule centers (bottom).

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 Appendix 6—figure 1. The local nematic order (Appendix 6—figure 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 Appendix 6—figure 1C, we show a kymograph of the local density of the crosslinking motor proteins.

The microtubule pair distribution function at steady-state (Appendix 6—figure 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.

Appendix 6—figure 2
Results for the confined microtubule-motor protein assembly simulations with Dcyl=0.75µm.

(A) A kymograph of the local microtubule packing fraction ϕlocal(x). Crosslinking motor proteins drive contraction of the system. Self organization of these regions leads to emergence of the BB-like state. (B) A kymograph of the local nematic order parameter Slocalx(x). The negative order parameter suggests that there is alignment of microtubules in the radial direction (YZ plane). (C) A kymograph of the density of the crosslinking motor proteins, C(x). (D) Pair distribution function for microtubule plus-ends (top) and microtubule centers (bottom).

In this case, the simulation volume is a cylinder with height 20µm. Over time, microtubules condense into a bottlebrush-like (BB) state with a hedgehog line defect. This consists of microtubules having a degree of alignment in the radial direction. The ends of the BB state contain a half-aster. Crosslinking motor proteins are highly concentrated along the central axis of the BB. Here is a kymograph for the local microtubule packing fraction, ϕlocal(x) (Appendix 6—figure 2A). We show the evolution of the local nematic order parameter, Slocalx(x), in Appendix 6—figure 2B. A negative Slocalx indicates a significant degree of radial alignment. Maximum radial alignment (the ideal bottlebrush state) is evidenced by a nematic order parameter value of 0.5 The condensation of microtubules is mediated by crosslinking motor proteins. In Appendix 6—figure 2C, we show a kymograph of the local density of the crosslinking motor proteins. Appendix 6—figure 2D depicts the microtubule pair distribution function. While there is a ring clearly visible for microtubule plus-ends (top plot), showing that this state is aster-like, there is significant density present along the X axis. This indicates that there is an accumulation of plus-ends throughout the line defect. Microtubule centers (bottom plot) are distributed uniformly along x while there is a decay in density along ρ. The absence of a ring indicates that this state is not aster-like. High density at ρ=0 suggests that microtubule centers tend to be stacked in x.

Ideal bottle-brush state

The ideal bottle-brush state (BB) consists of microtubules aligned in the radial direction directed away from a central line defect. A schematic and different views are shown in Appendix 6—figure 3. Microtubule orientation is indicated by the color wheel. For such a state, the local nematic order parameter along x has a value of −0.5 along the length of the BB.

Appendix 6—figure 3
The perfect bottle-brush state.

Microtubules are aligned in the XY planes such that there is a line defect along the z axis. (A) 3D view. (B) Side view. (C) Top view. Microtubule orientation is shown by the color wheel.

Appendix 7

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

Appendix 7—table 1
The parameters of the two springs controlling extension and bending, respectively.

The relation between B0 and E0 determine the equilibrium configuration of the two connected filaments. When B0E0+2dB, the straight configuration is the preferred configuration.

Rolespring stiffness constantfree length
Appendix 7—figure 1
The geometry of two short rigid straight fibers connected at a bending joint.

The separation is exaggerated to clearly show the geometry. p1 and p2 are the orientation norm vectors of the two segments. κE and κB are the stiffness of the spring for extension and the spring for bending. dB is the displacement distance from the joint rotation center. In the more detailed view of the deformed geometry of the two springs, a,b are the lengths of the two edges of the triangle. E=a2+b2+2abcosα. B=(a+dB)2+(b+dB)2+2(a+dB)(b+dB)cosα.

We can use two permanently bound springs for each joint, as shown in Appendix 7—figure 1, to implement the bending rigidity. The separations in the figure is exaggerated to show the geometry clearly. The energy of the two springs depend on their lengths E,B geometrically:

(74) U=12κE(E-E0)2+12κB(B-B0)2

With the deformed geometry, the lengths of the two springs are:

(75) E=a2+b2+2abcosα
(76) B=(a+dB)2+(b+dB)2+2(a+dB)(b+dB)cosα

When α0, the energy U of the two springs can be expanded as:

(77) U=12(κB(a+b+2dB-B0)2+κE(a+b-E0)2)+[κB(-a-dB)(b+dB)(a+b+2dB-B0)2(a+b+2dB)-abκE(a+b-E0)2(a+b)]α2+124κB((a+dB)(b+dB)(a2+dB(a+b)-ab+b2+dB2)(a+b+2dB-B0)(a+b+2dB)3+3(a+dB)2(b+dB)2(a+b+2dB)2)α4+124κE(ab(a2-ab+b2)(a+b-E0)(a+b)3+3a2b2(a+b)2)α4+O(α6).

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 κB+κE. The second α2 term governs the bending energy. We can tune the five parameters E0,B0,dB,κE,κB such that the connected segments reproduce the desired mechanical behavior of a flexible filament. Although the expansion Equation 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., κE= and E=E0. In this special case, the energy of the two springs depends only on U=12κB(B-B0)2. By imposing E=E0, we can solve for b:

(78) b=12(2a2cos(2α)-a2+2E02-2acos(α))

Then in this case U depends on α4 in the limit of α0. To simplify the notations of the expansion, we define:s=B0-E0-2dB.

The value s defines three cases of the equilibrium configuration:

  • s>0. The equilibrium configuration of the joint is a straight line, and the bending spring is compressed at equilibrium.

  • s=0. The equilibrium configuration of the joint is a straight line, and the bending spring is not compressed nor stretched at equilibrium.

  • s<0. The equilibrium configuration of the joint is bent.

For the first two cases, the equilibrium configuration is a straight line and we can expand U in the limit of α0:

(79) U=κBs22+dBκB(2a2-2aE0+E0(dB+E0))2E0(2dB+E0)sα2+dB2κB(2a2-2aE0+E0(dB+E0))28E02(2dB+E0)2α4+[3a4dB2κB2E02(2dB+E0)3+a4dB3κBE03(2dB+E0)3-a3dB2κBE0(2dB+E0)3-4a2dB2κB3(2dB+E0)3-11a2dB3κB6E0(2dB+E0)3+a4dBκB4E0(2dB+E0)3-7a2dBκBE012(2dB+E0)3+5adB3κB6(2dB+E0)3+5adB2κBE06(2dB+E0)3+adBκBE023(2dB+E0)3-dB2κBE0212(2dB+E0)3-dB3κBE012(2dB+E0)3-dB4κB24(2dB+E0)3-dBκBE0324(2dB+E0)3]sα4+O(α6).

With this form, it is clear that the bending energy is tunable with the parameter s, that is, how much the bending spring is strained in the equilibrium configuration. Note that s here is a constant determined by the lengths E,B,dB only. Therefore, the first term frac12κBs2 only ‘shifts’ the zero-point of the energy. This term does not contribute to the stretching or bending energy of the joint. When s=0, the leading order terms all vanish and U(α)α4. When s>0, the leading order terms are non-zero and the energy is asymptotically a quadratic function of α: U(α)α2.

Special case 2 If we further assume that κE= and E=E0=0, we have that a=b=0. This means the extension spring degenerates into a point joint between the two segments. In this case the energy U can be further simplified:

(80) U=12κB(-2dBcosα(dB+E0)+2dB2+2dBE0+E02+2dB+E0+s)2

The expansion of U as α0 is also further simplifed:

(81) U=κBs22+dBκBs(dB+E0)2(2dB+E0)α2+dBκB(dB+E0)(3dB(dB+E0)(2dB+E0)-s(dB2+dBE0+E02))24(2dB+E0)3α4+O(α6)

Here we have the same conclusion as the previous special case, that the dependence of U on α can be tuned between α4 and α2 by choosing a proper value of s.

Method 2: use bilateral constraints

Appendix 7—figure 2
The geometry of two short rigid straight fibers connected at a bending joint.

The separation is exaggerated to clearly show the geometry. p1 and p2 are the orientation norm vectors of the two segments. ω1 and ω2 are the rotational angular velocities. UB is the bending energy of this joint. α is the angle from p1 to p2. E is the bending rigidity modulus.

Here we briefly derive the constraint optimization formulation for handling the bending rigidity of flexible fibers with bilateral constraints. To fit in the geometric constraint formulation, we represent a long and flexible fiber as many short rigid straight fibers chained together by joints. The linear extension of each joint can be straightforwardly handled by the bilateral spring constraints as for those doubly bound motors. For the bending rigidity, we first realize that for each joint the two norm orientation vectors p1 and p2 form a plane. This plane is orthogonal to a unit norm vector

(82) T^=p1×p2|p1×p2|

For most relevant biological filaments, the bending rigidity is isotropic along different directions on a cross-section of the filament. In other words, the recovering torque is always co-linear with the vector p1×p2 and the recovering deformation is always in plane spanned by p1 and p2. This is important because we can simplify the deformation to in-plane rotations in the following derivations. Note that this plane can be different for each joint since each joint is handled as an independent constraint in our method.

There are different models of how the bending energy depends on the deformation, p1p2.

Case 1:


When the angle α between p1 and p2 is small, we have:


Case 2:


In this form when α0 the energy depends on the second order instead of the fourth order of the angle:


The following derivation and method still applies.

The two cases can be handled in the same way. In the following we derive the equations for the first case, where the second case only requires a simpler small α expansion in the derivation.

There is one more relation we can utilize to simplify the derivation. Assume that ω1 and ω2 have arbitrary directions, and to the first order of Δt the orientation vectors p1 and p2 rotates within Δt:

(83) p1p1+ω1×p1Δt
(84) p2p2+ω2×p2Δt

Then, the bending energy after this rotation is:

(85) UB=E[1-p1p2-Δt(p2(ω1×p1)+p1(ω2×p2))]2
(86) =E[1-p1p2-Δt(ω2-ω1)(p2×p1)]2,

where we have utilized the vector triple product identity:

(87) a(b×c)=b(c×a)=c(a×b).

This means, to the first order of Δt only the component of rotation ω1 and ω2 that is inside this plane spanned by p1,p2 affect the bending energy. Therefore, to the first order of Δt 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 n and n+1. We have, to the first order:

(88) αn+1=αn+(ω2n+1-ω1n+1)Δt.

The rotational mobility matrix for these two rods is:

(89) M=[M1R00M2R],

where M1R and M2R 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 UB. More specifically:

(90) ω1n+1=-M1RTn+1T^
(91) ω2n+1=M2RTn+1T^,

where the scalar torque Tn+1 is:

(92) Tn+1=-Eαn+1,3=-E[αn+(ω2n+1-ω1n+1)Δt]3
(93) =-E[αn,3+3αn,2ω2n+1Δt-3αn,2ω1n+1Δt]

where the higher order terms in Δt have been neglected. If the bending energy Case 2 is used, instead of Tn+1-Eαn+1,3 we have Tn+1-Eαn+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 Equation 90. Skipping the timestep index n, we can write the result in the same way as the bilateral Hookean spring constraints as:

(94) 0={DT[M1R00M2R]D+1K}[T]+13αn1ΔtTR,

where K=3Eαn,2 and the geometric matrix D defines the direction of torque on each rod:

(95) D=[-T^T].

The left side of Equation 94 means the motion of filament segments must satisfy the torque-deformation relation, while the right side means the torque can take any values.

Equation 94 is mathematically identical to the Hookean spring constraints and can be incorporated in the constraint minimization problem in the same way.

We can solve this two-segment problem analytically if the constraint optimization problem contains only Equation 94, in the absence of collisions and Hookean springs:

(96) (ω2-ω1)Δt=-13M1R+M2R2(M1R+M2R)+1Kα.

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, K as the bending rigidity modulus E increases to infinity. In this case, 1/K0 and the above solution is still stable, and is simplified to (ω2-ω1)Δt=-16α. This means the solution to Equation 94 has very strong temporal integration stability even when the deformation force is infinitely stiff, the same as what we discussed for the infinitely stiff Hookean spring case.

Data availability

The current manuscript is a computational study, so no data have been generated for this manuscript. The open-source modeling software is hosted at GitHub.


  1. Book
    1. Howard J
    Mechanics of Motor Proteins and the Cytoskeleton
    Sunderland: Sinauer associates.
    1. McIntosh JR
    (2016) Mitosis
    Cold Spring Harbor Perspectives in Biology 8:a023218.

Article and author information

Author details

  1. Wen Yan

    Center for Computational Biology, Flatiron Institute, New York, United States
    Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Software, Validation, Visualization, Writing - original draft, Writing - review and editing
    For correspondence
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9189-0840
  2. Saad Ansari

    Department of Physics, University of Colorado Boulder, Boulder, United States
    Data curation, Investigation, Validation, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
  3. Adam Lamson

    1. Center for Computational Biology, Flatiron Institute, New York, United States
    2. Department of Physics, University of Colorado Boulder, Boulder, United States
    Formal analysis, Investigation, Methodology, Software, Validation, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2616-2801
  4. Matthew A Glaser

    Department of Physics, University of Colorado Boulder, Boulder, United States
    Supervision, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8366-5598
  5. Robert Blackwell

    Center for Computational Biology, Flatiron Institute, New York, United States
    Data curation, Validation
    Competing interests
    No competing interests declared
  6. Meredith D Betterton

    1. Center for Computational Biology, Flatiron Institute, New York, United States
    2. Department of Physics, University of Colorado Boulder, Boulder, United States
    3. Department of Molecular, Cellular, and Developmental Biology, University of Colorado Boulder, Boulder, United States
    Conceptualization, Funding acquisition, Investigation, Project administration, Resources, Supervision, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-5430-5518
  7. Michael Shelley

    1. Center for Computational Biology, Flatiron Institute, New York, United States
    2. Courant Institute, New York University, New York, United States
    Conceptualization, Funding acquisition, Investigation, Project administration, Resources, Supervision, Writing - review and editing
    For correspondence
    Competing interests
    No competing interests declared


National Science Foundation (DMR-2004469)

  • Michael Shelley

National Science Foundation (CMMI-1762506)

  • Michael Shelley

National Science Foundation (DMS-1821305)

  • Saad Ansari
  • Matthew A Glaser
  • Meredith D Betterton

National Science Foundation (ACI-1532235)

  • Saad Ansari
  • Matthew A Glaser
  • Meredith D Betterton

National Science Foundation (ACI-1532236)

  • Saad Ansari
  • Matthew A Glaser
  • Meredith D Betterton

National Science Foundation (RGM124371A)

  • Saad Ansari
  • Matthew A Glaser
  • Meredith D Betterton

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.


MJS acknowledges support from NSF grants DMR-2004469 and CMMI-1762506. SA, ARL, MAG, and MB acknowledge support from NSF grants DMS-1821305, and NIH grant RGM124371A. This work utilized resources from the University of Colorado Boulder Research Computing Group, which is supported by the National Science Foundation (awards ACI-1532235 and ACI-1532236), the University of Colorado Boulder, and Colorado State University. We thank Prof. Dimitrios Vavylonis for discussions on implementing flexible filaments.

Version history

  1. Preprint posted: September 16, 2021 (view preprint)
  2. Received: September 23, 2021
  3. Accepted: April 24, 2022
  4. Version of Record published: May 26, 2022 (version 1)


© 2022, Yan et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.


  • 1,937
  • 383
  • 12

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Wen Yan
  2. Saad Ansari
  3. Adam Lamson
  4. Matthew A Glaser
  5. Robert Blackwell
  6. Meredith D Betterton
  7. Michael Shelley
Toward the cellular-scale simulation of motor-driven cytoskeletal assemblies
eLife 11:e74160.

Share this article

Further reading

    1. Computational and Systems Biology
    2. Medicine
    Zachary Shaffer, Roberto Romero ... Nardhy Gomez-Lopez
    Research Article


    Preterm birth is the leading cause of neonatal morbidity and mortality worldwide. Most cases of preterm birth occur spontaneously and result from preterm labor with intact (spontaneous preterm labor [sPTL]) or ruptured (preterm prelabor rupture of membranes [PPROM]) membranes. The prediction of spontaneous preterm birth (sPTB) remains underpowered due to its syndromic nature and the dearth of independent analyses of the vaginal host immune response. Thus, we conducted the largest longitudinal investigation targeting vaginal immune mediators, referred to herein as the immunoproteome, in a population at high risk for sPTB.


    Vaginal swabs were collected across gestation from pregnant women who ultimately underwent term birth, sPTL, or PPROM. Cytokines, chemokines, growth factors, and antimicrobial peptides in the samples were quantified via specific and sensitive immunoassays. Predictive models were constructed from immune mediator concentrations.


    Throughout uncomplicated gestation, the vaginal immunoproteome harbors a cytokine network with a homeostatic profile. Yet, the vaginal immunoproteome is skewed toward a pro-inflammatory state in pregnant women who ultimately experience sPTL and PPROM. Such an inflammatory profile includes increased monocyte chemoattractants, cytokines indicative of macrophage and T-cell activation, and reduced antimicrobial proteins/peptides. The vaginal immunoproteome has improved predictive value over maternal characteristics alone for identifying women at risk for early (<34 weeks) sPTB.


    The vaginal immunoproteome undergoes homeostatic changes throughout gestation and deviations from this shift are associated with sPTB. Furthermore, the vaginal immunoproteome can be leveraged as a potential biomarker for early sPTB, a subset of sPTB associated with extremely adverse neonatal outcomes.


    This research was conducted by the Perinatology Research Branch, Division of Obstetrics and Maternal-Fetal Medicine, Division of Intramural Research, Eunice Kennedy Shriver National Institute of Child Health and Human Development, National Institutes of Health, U.S. Department of Health and Human Services (NICHD/NIH/DHHS) under contract HHSN275201300006C. ALT, KRT, and NGL were supported by the Wayne State University Perinatal Initiative in Maternal, Perinatal and Child Health.

    1. Computational and Systems Biology
    2. Genetics and Genomics
    Ardalan Naseri, Degui Zhi, Shaojie Zhang
    Research Article

    Runs of homozygosity (ROH) segments, contiguous homozygous regions in a genome were traditionally linked to families and inbred populations. However, a growing literature suggests that ROHs are ubiquitous in outbred populations. Still, most existing genetic studies of ROH in populations are limited to aggregated ROH content across the genome, which does not offer the resolution for mapping causal loci. This limitation is mainly due to a lack of methods for the efficient identification of shared ROH diplotypes. Here, we present a new method, ROH-DICE, to find large ROH diplotype clusters, sufficiently long ROHs shared by a sufficient number of individuals, in large cohorts. ROH-DICE identified over 1 million ROH diplotypes that span over 100 SNPs and are shared by more than 100 UK Biobank participants. Moreover, we found significant associations of clustered ROH diplotypes across the genome with various self-reported diseases, with the strongest associations found between the extended HLA region and autoimmune disorders. We found an association between a diplotype covering the HFE gene and hemochromatosis, even though the well-known causal SNP was not directly genotyped or imputed. Using a genome-wide scan, we identified a putative association between carriers of an ROH diplotype in chromosome 4 and an increase in mortality among COVID-19 patients (P-value=1.82×10-11). In summary, our ROH-DICE method, by calling out large ROH diplotypes in a large outbred population, enables further population genetics into the demographic history of large populations. More importantly, our method enables a new genome-wide mapping approach for finding disease-causing loci with multi-marker recessive effects at a population scale.