Toward the cellularscale simulation of motordriven cytoskeletal assemblies
Abstract
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 hardbody 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.
https://doi.org/10.7554/eLife.74160.sa0Introduction
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 selforganization (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 selforganization, 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 10^{4}–10^{7} or more filaments (Petry, 2016). While current simulations may reach $O({10}^{4}{10}^{5})$ filaments (Belmonte et al., 2017; Strübing et al., 2020), molecular modeling has required significant compromises in treating steric interactions and motorproteins.
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 (welldefined) 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 $\mathrm{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, potentialbased models limit simulations to short timescales. To circumvent this limitation, here we utilize our recently developed constraint method to enforce hardcore repulsion between particles (Anitescu et al., 1996; Yan et al., 2019). We further develop constraintbased 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).
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 minusend (modeling the activity of dynein). Although the microtubules are initially unorganized (C1), the combination of motor crosslinking and walking causes the microtubule minusends to contract into the center of a large aster (C2). The motordriven steric interactions between filaments, however, eventually fragment this into smaller asters and bottlebrushlike structures (C3,C4). This simulation displays the complex interplay between steric and crosslinking forces in determining the dynamics and steady state configurations of cytoskeletal materials.
Methodology
In this work, we model filaments as rigid spherocylinders. (While not presented here, flexible filaments can be modeled within our framework as segmented, jointed filaments; See Appendix 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 ${v}_{m}\mathrm{\Delta}t$ and doubly bound motors move ${v}_{F}\mathrm{\Delta}t$ along the filaments. Here, ${v}_{F}$ is the motor stepping velocity that depends on force on the motor head (Gao et al., 2015a):
where ${F}_{\mathrm{proj}}$ is the projection of tether force along filament in the stepping direction. As typically found experimentally, this stepping model means that if ${F}_{\mathrm{proj}}$ is assisting stepping, the velocity saturates at v_{m}; while for ${F}_{\mathrm{proj}}$ hindering stepping, stepping is halted when ${F}_{\mathrm{proj}}={F}_{\mathrm{stall}}$.
2. Crosslinker binding and unbinding
In filament networks, the spatial variation of unbound and bound motors is integral to network selforganization. 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 steadystate configuration. Entropic forces bundle and increase overlaps among crosslinked filaments (Lansky et al., 2015; Gaska et al., 2020), and freeenergydependent binding kinetics contribute to organization of cortical microtubules (Allard et al., 2010) and induce actin bundling (Yang et al., 2006).
Adhoc 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 nonstretched 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 ${v}_{m}=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 boundunbound 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 (${S}_{A}$ or ${S}_{B}$), or both heads (doubly) bound ($D$). For each timestep $\mathrm{\Delta}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\rightleftharpoons ({S}_{A},{S}_{B})\rightleftharpoons D$). The transition probabilities are modeled as inhomogeneous Poisson processes with the cumulative probability function
The transitions $U\rightleftharpoons ({S}_{A},{S}_{B})$ do not stretch or compress the tether and so do not depend on tether deformation energy. However, the transitions $({S}_{A},{S}_{B})\rightleftharpoons D$ do account for tether deformation energy (Table 1).
3. Filament dynamics
We sought to develop a stable, largetimestep 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 $\mathrm{\Delta}t$. The first arises in models that use a stiff repulsive pairwise potential to prevent filament overlaps. For example, the LennardJones potential $V\sim {(\sigma /r)}^{12}{(\sigma /r)}^{6}$, where $r$ is the separation between filaments, is so steeply varying that it requires small $\mathrm{\Delta}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 ${\mathrm{\ell}}_{f}$ relaxes according to ${\dot{\mathrm{\ell}}}_{f}=\lambda ({\mathrm{\ell}}_{f}{\mathrm{\ell}}_{0})$, where ${\mathrm{\ell}}_{0}$ is the preferred length and $\lambda =N{\kappa}_{\mathrm{xl}}/(4\pi \eta L/\mathrm{log}(2L/{D}_{\mathrm{fil}}))$ (Howard, 2001). Explicit timestepping schemes require $\mathrm{\Delta}t<C/\lambda$, for some constant $C$. For $N=10$ motors, tether stiffness ${\kappa}_{\mathrm{xl}}\approx 100\text{}\mathrm{pN}\text{}{\mathrm{\mu}\mathrm{m}}^{1}$, and slender body drag coefficient $4\pi \eta L/\mathrm{log}(2L/{D}_{\mathrm{fil}})\approx 0.003\text{}\mathrm{pN}\text{}\mathrm{s}\text{}{\mathrm{\mu}\mathrm{m}}^{1}$ for $1\text{}\mathrm{\mu}\mathrm{m}$long microtubules in aqueous solvent, we have $1/\lambda \approx 3\text{\xd7}{10}^{6}\text{}\mathrm{s}$.
We overcome these difficulties with a novel, linearized implicit Euler timestepping scheme, which extends on our previous work on enforcing nonoverlap conditions (Yan et al., 2019). This technique is inspired by constraintbased methods for granular flow (Tasora et al., 2013). When collisions occur between filaments, the minimal distance between them attains ${\mathrm{\Phi}}_{\mathrm{col}}=0$ with collision force ${\gamma}_{\mathrm{c}\mathrm{o}\mathrm{l}}>0$. If not colliding, ${\mathrm{\Phi}}_{\mathrm{c}\mathrm{o}\mathrm{l}}>0$ and ${\gamma}_{\mathrm{col}}=0$. This mutually exclusive condition is called a complementarity constraint, written as $0\le {\mathrm{\Phi}}_{\mathrm{col}}\u27c2{\gamma}_{\mathrm{col}}\ge 0$. If one crosslinking motor connects these two filaments, its length ${\mathrm{\ell}}_{f}$ and force magnitude ${\gamma}_{\mathrm{xl}}$ satisfy the Hookean spring model ${\gamma}_{\mathrm{xl}}={\kappa}_{\mathrm{xl}}({\mathrm{\ell}}_{f}{\mathrm{\ell}}_{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 $\mathit{x}\in {\mathbb{R}}^{3}$ in the lab frame and its orientation $\mathit{\theta}=[s,\mathit{p}]\in {\mathbb{R}}^{4}$ as a quaternion (Delong et al., 2015). $[s,\mathit{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 singularityfree nature. The geometric configuration at time $t$ for all $N$ filaments can be written as a column vector with $7N$ entries:
Similarly, we use the vectors $\mathcal{U},\mathcal{F}\in {\mathbb{R}}^{6N}$ to represent the translational & angular velocities, and forces & torques of all particles, respectively. We relate $\mathcal{U}$ to $\mathcal{F}$ via a mobility matrix $\mathcal{M}\in {\mathbb{R}}^{6N\times 6N}$, dependent only upon the geometry $\mathcal{C}$, and relate $\mathcal{U}$ to $\dot{\mathcal{C}}(t)$ via a geometric matrix $\mathcal{G}$:
Because the biological filaments we consider mostly have lengths on the $\mathrm{nm}$ to $\mathrm{\mu}\mathrm{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 nonconstraint.
For unilateral constraints, we define the grand distance vector ${\mathbf{\Phi}}_{u}={\left[{\mathrm{\Phi}}_{u,1},{\mathrm{\Phi}}_{u,2},\cdots ,{\mathrm{\Phi}}_{u,{N}_{u}}\right]}^{T}\in {\mathbb{R}}^{{N}_{u}},$ where each ${\mathrm{\Phi}}_{u,j}$ is the minimum distance between a pair of filaments. Similarly, for bilateral constraints we define the grand distance vector $\mathbf{\Phi}}_{b}={\left[{\ell}_{f,1},{\ell}_{f,2},\cdots ,{\ell}_{f,{N}_{b}}\right]}^{T}\in {\mathbb{R}}^{{N}_{b}$, containing the length ${\mathrm{\ell}}_{f,j}$ of the doubly bound motor $j$. There are in total ${N}_{u}$ possibly colliding pairs of filaments and ${N}_{b}$ crosslinking motors. The force magnitude corresponding to these constraints are also written as vectors, $\mathit{\gamma}}_{u}={\left[{\gamma}_{u,1},{\gamma}_{u,2},\cdots ,{\gamma}_{u,{N}_{u}}\right]}^{T}\in {\mathbb{R}}^{{N}_{u}$ and $\mathit{\gamma}}_{b}={\left[{\gamma}_{b,1},{\gamma}_{b,2},\cdots ,{\gamma}_{b,{N}_{b}}\right]}^{T}\in {\mathbb{R}}^{{N}_{b}$. The two types of constraints can be summarized as:
Here, $\mathbf{\Phi}}_{u$ and $\mathit{\gamma}}_{u$ satisfy the complementarity (collision) constraints, while $\mathbf{\Phi}}_{b$ and $\mathit{\gamma}}_{b$ satisfy the Hookean spring law. Here, $\mathcal{K}\in {\mathbb{R}}^{{N}_{b}\times {N}_{b}}$ is a diagonal matrix consisting of all the stiffness constants, while $\mathbf{\Phi}}_{b}^{0$ represents the rest length of every crosslinking motor.
Equations 4 and 5 define a differentialvariationalinequality (DVI). This is solvable when closed by a geometric relation mapping the force magnitude $\mathit{\gamma}}_{u$ and $\mathit{\gamma}}_{b$ to the force vectors $\mathcal{F}}_{u$ and $\mathcal{F}}_{b$:
where $\mathcal{D}}_{u$ and $\mathcal{D}}_{b$ 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 $\mathrm{\Delta}t=h$ at timestep $k$:
The unknowns to be solved for at every timestep are the constraint (collision and motor tether) force magnitude $\mathit{\gamma}}_{u}^{k},{\mathit{\gamma}}_{b}^{k$. This is a nonlinear DVI because $\mathbf{\Phi}}_{u}^{k+1$, $\mathbf{\Phi}}_{b}^{k+1$ are nonlinear functions of geometry $\mathcal{C}}^{k+1$, although $\mathcal{C}}^{k+1$ is linearly dependent on $\mathit{\gamma}}_{u}^{k$ and $\mathit{\gamma}}_{b}^{k$. For a small timestep ($h\to 0$), this nonlinearity can be linearized by Taylor expansion, for example, $\mathbf{\Phi}}_{u}^{k+1}={\mathbf{\Phi}}_{u}^{k}+h{\mathrm{\nabla}}_{\mathcal{C}}{\mathbf{\Phi}}_{u}{\mathcal{G}}^{k}{\mathcal{U}}^{k$. Then, this nonlinear DVI can be converted to a convex quadratic programming problem (Nocedal and Wright, 2006) (details in Appendix D):
Here, $\mathit{\gamma}}^{k}=[{\mathit{\gamma}}_{u}^{k},{\mathit{\gamma}}_{b}^{k}]\in {\mathbb{R}}^{{N}_{u}+{N}_{b}$ is a column vector, and
One way to understand the constraint optimization method is that the implicit temporal integration ‘jumps’ on a timescale that bypasses the relaxation timescales of unilateral and bilateral constraints (collisions and crosslinking motor springs). In the limit of motor tethers being infinitely stiff ($\mathcal{K}}^{1}\to \mathbf{0$), the quadratic term coefficient matrix $\mathit{M}$ is still symmetricpositivesemidefinite (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 noncompliant joints.
Instantiation in a massively parallel computing environment
Our methods naturally lend themselves to highperformance 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 nearneighbor 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 BarzilaiBorwein Projected Gradient Descent (BBPGD) solver (Yan et al., 2019) because the gradient $\mathrm{\nabla}f=\mathit{M}\mathit{\gamma}+\mathit{q}$ is efficiently computed by one parallel sparse matrixvector multiplication operation.
aLENS is written in a modular design using standard objectoriented 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 purefilament phase, and shown to accurately reproduce the equation of state and the isotropicnematic 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, ${N}_{A}$ is the number of active motors and $N$ is the total number of motors (active and inactive). The microtubule velocity increases as ${N}_{A}/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\text{}\mathrm{\mu}\mathrm{m}\text{}{\mathrm{s}}^{1}$, so the sliding velocity at ${N}_{A}/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 ${N}_{A}/N$.
Selfstraining 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 manymicrotubule 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 ‘selfstraining motion’, with the system interpreted as being composed to two polar microtubule gels whose interconnecting motors pulled them past one another.
We simulate this experiment using 3,000 model microtubules with $L=0.5\text{}\mathrm{\mu}\mathrm{m}$. Initially the filaments are confined in a tube of diameter $D=1\text{}\mathrm{\mu}\mathrm{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\text{}\mathrm{\mu}\mathrm{m}$. There are approximately 25 motors per microtubule according to the experimental estimates, and in our simulations we vary the motortomicrotubule number ${N}_{m}$ from 10 to 30. There is no accurate measurement for the XCTK2 motor in these experimental conditions. Therefore, we used experimental estimates of $46\text{}\mathrm{nm}\text{}{\mathrm{s}}^{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 ${V}_{x}$ is largely independent of the number of motors ${N}_{m}$ and the local average polarity ${P}_{x}$ over the range simulated.
Intuitively, the straining velocity ${V}_{x}$ is predominantly determined by the free walking velocity of the motors in limit of many crosslinkers. From our simulations, we find a straining velocity of approximately $26\text{}\mathrm{nm}\text{}{\mathrm{s}}^{1}$, close to the experimental measurement of $18.6\pm 0.9\mathrm{nm}\text{}{\mathrm{s}}^{1}$.
Largescale parallelization efficiency
Simulation of cellularscale cytoskeletal assemblies requires methods that can reach large system sizes and timescales. Therefore, we developed aLENS to efficiently utilize modern highperformance computing resources. Millions of objects and constraints can be simulated with aLENS. Figure 4 shows detailed parallel efficiency measurements for one largescale 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.
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 $\mathrm{\Delta}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 $1\text{}\mathrm{s}$ physical time per day, using a timestep size of $1.0\text{\xd7}{10}^{5}\text{}\mathrm{s}$.
Results
Here, we illustrate the ability to use aLENS to study the interplay between microscopic dynamics and macroscopic order in active cytoskeletal assemblies. The specific examples shown here are the formation and extension of a band of microtubule bundles, polarity sorting of short microtubules on a spherical shell, the development of asters with and without thermal fluctuations, and the effect of confinement on assembling microtubulemotor mixtures. For the results presented here, all simulations were conducted in solvent with viscosity $\eta =0.01\text{}\mathrm{pN}\text{}\mathrm{s}\text{}{\mathrm{\mu}\mathrm{m}}^{2}$ at room temperature, using a fixed timestep $\mathrm{\Delta}t={10}^{4}\text{}\mathrm{s}$ unless otherwise stated.
Bundle formation and buckling in a filament band
Microtubules driven by crosslinking motors can bundle; sliding of microtubules within the bundles causes them to fracture dynamically (Sanchez et al., 2012; Foster et al., 2015; Roostalu et al., 2018). We study such phenomena through a largescale simulation of 100,000 filaments modeling microtubules and 500,000 minusenddirected 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 wellaligned 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).
The macroscopic stresses and dynamics depend on the spatial organization of filaments and motordriven sliding. To characterize this, we measure the joint probability distribution of the local nematic order parameter ${S}_{\mathrm{local}}$ and the number ${N}_{d}$ of neighboring filaments crosslinked to a filament (Figure 5D). While the network contracts, the distribution of ${N}_{d}$ 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 nearperfectly nematic (${S}_{\mathrm{local}}\approx 1$), although lessordered 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 labframe $x$axis, we observe left and rightmoving filaments that speed up early in the simulation, and then maintain constant average velocities at later time ($t>4\text{}\mathrm{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 denselypacked 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 $t\ge 4\text{}\mathrm{s}$ (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 wellstudied 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/{D}_{\mathrm{fil}}=10$ are confined between two closely spaced concentric spherical shells at 40% volume fraction. The shell gap is $\mathrm{\Delta}R=0.102\text{}\mathrm{\mu}\mathrm{m}$, shorter than the filament length, with $\mathrm{\Delta}R/{D}_{\mathrm{fil}}\approx 4$ so filaments can move over each other in a restricted way. The filaments are initialized such that the nematic directors are along the meridians everywhere. 200,000 motors, modeled after kinesin5 tetramers, drive relative filament sliding (Figure 6A). Brownian motion is modeled at room temperature $300\text{}\mathrm{K}$ and timestep $\mathrm{\Delta}t$ is set to $1\text{\xd7}{10}^{5}\text{}\mathrm{s}$. Motors move toward minus ends of bound filaments at ${v}_{m}=1.0\text{}\mathrm{\mu}\mathrm{m}\text{}{\mathrm{s}}^{1}$. Once they reach the minus ends, they immediately detach.
Motors walk along the filaments, driving sliding of antiparallel filaments (Figure 6B). This leads to polaritysorted regions at the north and south ‘poles’ of the sphere, meaning that the filament orientation $\mathit{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 polaritysorted regions approximately parallel to the polarity direction. Instead, on the sphere the boundaries between polaritysorted 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 minusends (Figure 6C2, C5 and G). This accumulation is illustrated by the positive correlation between motor density $n$ and $\mathrm{\nabla}\cdot \mathit{p}$ at $t=4\text{}\mathrm{s}$ 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, $\mathrm{\nabla}\cdot \mathit{p}>0$ in regions where minusends meet minusends and vice versa in regions where plusends meet plusends. Minusend directed motors accumulate in regions with $\mathrm{\nabla}\cdot \mathit{p}>0$, while plusend motors accumulate in regions with $\mathrm{\nabla}\cdot \mathit{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 (endpausing, EP), the filament network contracts (Fig. 6E15) with volume fraction increases from 40% to 60% and eventually freezes at $t=0.27\text{}\mathrm{s}$. We observe neither substantial polarity sorting nor motor accumulation. This indicates that the ability of motors to continuously walk, without endpausing, is crucial to effective polarity sorting.
Aster formation in bulk
Aster formation is driven by motor pausing at ends of rigid filaments (endpausing). 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, minusenddirected, endpausing 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).
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=25\text{}\mathrm{nm}={D}_{\mathrm{fil}}$ and $r=78\text{}\mathrm{nm}={\mathrm{\ell}}_{0}+{D}_{\mathrm{fil}}$ 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\text{}\mathrm{\mu}\mathrm{m}$ in $1\text{}\mathrm{s}$. This allows filaments to diffuse until they are captured in regions of high motor density, such as aster centers. Furthermore, with thermal fluctuations the asters themselves diffuse, which leads to aster coalescence (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 motorinduced stress ${\mathrm{\Pi}}^{Linker}$ initially increases quickly, reaching a peak at roughly $t=4\text{}\mathrm{s}\sim 5\text{}\mathrm{s}$, similar to the behavior during bundle contraction shown above (Figure 5F), before declining. The average time required for motors to walk to filament ends, ${\tau}_{walk}=L/{v}_{m}\approx 5\text{}\mathrm{s}$, determines the initial contraction timescale. After reaching minus ends, motors pause and relax toward their equilibrium lengths. As a result, both the motor and collision stress grow in magnitude as more motors accumulate at minus ends.
Confined filamentmotor 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 highdensity crosslinked filament aggregates that coexist with a relatively low density vapor of noncrosslinked filaments. In bulk systems as shown above and in previous work, endpausing 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\text{}\mathrm{\mu}\mathrm{m}$ long filaments at a fixed packing fraction ($\varphi =0.16$), confined in two cylinders with diameters ${D}_{\mathrm{cyl}}=0.25\text{}\mathrm{\mu}\mathrm{m}$ and $0.75\text{}\mathrm{\mu}\mathrm{m}$.
For a smalldiameter cylinder where one filament length can fit across the cylinder (${D}_{\mathrm{cyl}}/L=1$, ${D}_{\mathrm{cyl}}=0.25\text{}\mathrm{\mu}\mathrm{m}$), the cylinder is too narrow for asters to form. Instead, motor sliding and endpausing drive the filaments into polaritysorted bilayers (PSBs, Figure 8A and movie Figure 8—video 1). A single polaritysorted bilayer contains a central interface of highly crosslinked filament minusends between two antiparallel polar layers of filaments (Figure 8). At steady state, the system consists of individual PSBs separated by lowdensity vapor regions containing few motors. As expected, the local nematic order parameter ${S}_{\mathrm{local}}^{x}(x)$ nearly reaches 1 within PSBs. Even the the vapor phase is close to nematic ${S}_{\mathrm{local}}^{x}\approx 0.6$ (Figure 8A4), due to the strong confinement effect.
Next we increased the diameter of the cylinder to ${D}_{\mathrm{cyl}}/L=3$ (${D}_{\mathrm{cyl}}=0.75\text{}\mathrm{\mu}\mathrm{m}$) to weaken the confinement (Figure 8B and movie Figure 8—video 2). Here, the polaritysorted bilayers are not present, because the larger cylinder diameter allows filaments to reorient and organize into bottlebrushlike 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 $\sim 2.5$ times the vapor when compared to the PSBs (Figure 8B4, red line).
Discussion
We designed aLENS to (i) model crosslinking motor kinetics conforming to an underlying free energy landscape, (ii) circumvent the timescale limitation imposed by conventional explicit timestepping methods, and (iii) efficiently utilize modern parallel computing resources to allow simulation of cellularscale systems. This efficient framework allows both modeling the individual building cytoskeletal building blocks (filaments, motors) and gathering mesoscale statistical information such as stress and order parameters from a large system. This multiscale capability will make it possible to directly compare simulations with experimental observations on mesoscopic and macroscopic scales over timescales from seconds to minutes.
The aLENS framework is not limited to a specific motor model. Because of the modular design of the motor code, the motor model can be extended to include additional physics such as forcedependent binding and unbinding rates, or even entirely replaced, say, with a passive crosslinker or other model. Dynamic instability and branching of cytoskeletal filaments can also be integrated with the constraint minimization problem, as we showed previously in modeling the divisiondriven 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 beadspring 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 coarsegrained 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 crosslinkmediated interactions play an important role.
aLENS has been opensourced on GitHub: https://github.com/flatironinstitute/aLENS (copy archived at swh:1:rev:f2dd484f82443735562ad7b480fe7ed9fc020fb0; Adam, 2022) and precompiled binary executable is available on DockerHub: https://hub.docker.com/r/wenyan4work/alens. Our GitHub documentation provides a clear roadmap for developing additional userspecific modules.
Appendix 1
Crosslinker and motor properties
Appendix 2
Crosslinker binding and unbinding
Kinetic MonteCarlo: crosslinking proteinfilament interactions
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 meanfield 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
where $\mathrm{H}$ is a head and $\mathrm{B}$ is a binding site on a filament. The association constant ${K}_{\mathrm{a}}$ of the heads to binding sites is described by the equilibrium equation
where $[X]$ defines the concentration of substance $X$. The association constant has units of inverse molarity, and relates the the on and offrate constants $k}_{\mathrm{o}\mathrm{n},\phantom{\rule{0.056em}{0ex}}S$ and $k}_{\mathrm{o}\mathrm{f}\mathrm{f},\phantom{\rule{0.056em}{0ex}}S$.
Unbound crosslinkers are modeled as diffusing points with center of mass positions ${\mathit{x}}_{o}(t)$ and diffusion constants d_{u}. The heads of a crosslinker have spatial and timedependent concentrations $[H]=c(\mathit{x},t)$. The head binding rate is the volume integral over the product of the onrate constant, binding site density, and crosslinker concentration
where $k}_{\mathrm{o}\mathrm{n},\phantom{\rule{0.056em}{0ex}}S}^{A}={K}_{\mathrm{a}}^{A}{k}_{\mathrm{o}\mathrm{f}\mathrm{f},\phantom{\rule{0.056em}{0ex}}S}^{A$ and $\u03f5$ is the linear binding site density along filaments. The lab position along the ith filament ${\mathit{x}}_{i}(s)$ is parameterized by $s$.
The binding probability in a timestep $\mathrm{\Delta}t$ is an inhomogeneous Poisson process with the cumulative probability function
We assume the tight binding limit $k}_{\mathrm{o}\mathrm{n},\phantom{\rule{0.056em}{0ex}}S}\gg {k}_{\mathrm{o}\mathrm{f}\mathrm{f},\phantom{\rule{0.056em}{0ex}}S$ and do not consider multiple binding and unbinding events of one crosslinker during a timestep $\mathrm{\Delta}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(\mathit{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 ${\mathrm{\ell}}_{o}$. Because of the tether’s stiffness, the radius of gyration of an unbound protein is assumed to be ${r}_{g}={\mathrm{\ell}}_{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, $\sqrt{6{d}_{U}\mathrm{\Delta}t}<{r}_{g}$, rotational diffusion dominates and we approximate heads to be within a sphere of radius ${r}_{c,S}={r}_{g}$ centered at $\mathit{x}}_{o$. Realistically, the head distributions can vary within this volume but because ${\mathrm{\ell}}_{o}$ is small compared to filament lengths, we approximate the head distribution as being uniform, that is, $c(\mathit{x},t)=(4\pi {r}_{c,S}^{3}/3{)}^{1}(1\mathrm{\Theta}(\mathit{x}{r}_{c,S}))$, where $\mathrm{\Theta}(x)$ is the Heaviside step function. For larger crosslinking proteins, more detailed spatial distributions may be calculated using freelyjointed or wormlike chain models.
For uniformly distributed heads, the head binding rate is
where ${L}_{in}$ is the filament $i$’s length segment within ${r}_{c,S}$. To account for cylindrical filaments with diameter ${D}_{\mathrm{fil}}$, we augment the binding radius such that ${r}_{c,S}={r}_{g}+{D}_{\mathrm{fil}}/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 ${R}_{\mathrm{o}\mathrm{n},S}(t)\approx {R}_{\mathrm{o}\mathrm{n},S}({t}_{i})$, for $t\in \left({t}_{i},{t}_{i}+\mathrm{\Delta}t\right)$.
However, it is uncommon that an unbound crosslinker or motor will diffuse less than r_{g} in a timestep, and so we must account for the protein’s translational diffusion. The diffusion equation models the mean spatial distribution of a unbound crosslinking protein’s center
which has the solution
If the characteristic diffusion length $\sqrt{{d}_{U}\mathrm{\Delta}t}\gg {r}_{g}$, 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, ${c}_{o}(\mathit{x},t)\approx c(\mathit{x},t)$. Substituting the binding rate equation (12) and the solution to the diffusion equation (16) into the integral of equation (13) gives
For straight, rigid filaments, we take the volume and time integrals while reparameterizing $\mathit{x}{\mathit{x}}_{o}{}^{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
Integrating over s_{i} gives the binding probability of one crosslinker head to a single filament. The total binding probability is then
where $N$ is the number of filaments surrounding the crosslinker head.
Calculating the binding probability from this function and ensuring that the protein unbinds so that detailedbalance is satisfied is computationally prohibitive. Instead, setting ${r}_{c,S}$ to the root mean square diffusion distant $\sqrt{6{d}_{U}\mathrm{\Delta}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
If implementing equation 14, the protein unbinds into a uniform sphere of radius ${r}_{c,S}$. This ensures crosslinkers bind to and from regions in a way that satisfies detailedbalance.
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 s_{i}, the binding rate constant ${k}_{\mathrm{on},[},D]$ for binding to a location s_{j} on filament $j$ is
where ${k}_{\mathrm{o},[},D]={k}_{\mathrm{off},[},D]({U}_{i,j}=0)$ is the unbindingrate of crosslinking proteins when no force is applied, ${K}_{\mathrm{e}}$ is a binding association constant similar to ${K}_{a}$, and ${U}_{i,j}$ is the free energy contribution from the tether
Before crosslinking, the unbound motor head explores an effective volume ${V}_{\mathrm{bind}}$ centered around the bound motor head. Not considering steric interactions with filaments, this volume is the free head’s position weighted by the Boltzmann factor integrated over all space.
We impose an integration cutoff radius ${r}_{c,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 catchbondlike behavior where proteins remain crosslinked for longer if under tension and release quicker if compressed. We replicate this behavior with the function
where λ and x_{c} are the energy factor and characteristic length specifying the behavior of the energy and forcedependent binding/unbinding, respectively. For values of ${x}_{c}<0$, you see catchbond like behavior whereas values of ${x}_{c}>0$ exhibit slip bond behavior Walcott, 2008; Edelmaier et al., 2020. This formalism can also be used to add in angle dependence.
When we include an effective energy and/or force dependence, the unbinding rate becomes
This does not change the final stored energy of either bound state but does effect the frequency at which the motors will switch between having one head bound and crosslinking.
Meanfield theory for crosslinking proteins
We expand on our previous meanfield motor density model to include motors that have dissimilar heads, diffusion and walking in singly and doubly bound states, and a timedependent 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(N1)$ crosslinking densities ${\psi}_{i,j}^{A,B}$, $2N$ singly bound densities ${\chi}_{i}^{A}$ and ${\chi}_{i}^{B}$, and an unbound density $C$ to model all crosslinking proteins between $N$ filaments. By convention,${\psi}_{i,j}^{A,B}={\psi}_{j,i}^{B,A}$
For heads $A$ and $B$, crosslinking diffusion constants ${d}_{i,j}^{A}$ and ${d}_{i,j}^{B}$, the drag speeds ${v}_{drag,i,j}^{A}$ and ${v}_{drag,i,j}^{B}$, and walking speeds ${v}_{walk,i,j}^{A}$ and ${v}_{walk,i,j}^{B}$ have been shown to depend on the force exerted on the binding heads and thus the stretch of the tether $\mathrm{\ell}({s}_{i},{s}_{j})$. No tether force acts on singly bound motor proteins but the singly bound diffusion constants ${d}_{i}^{A}$ and ${v}_{i}^{B}$ and walking speeds ${v}_{drag,i}^{A}$ and ${v}_{drag,i}^{B}$ may depend on s_{i} 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
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 $\mathit{x}$ in the lab frame and its orientation $\mathit{\theta}=[s,\mathit{p}]\in {\mathbb{R}}^{4}$ as a quaternion Delong et al., 2015. This $\mathit{p}$ is the vector component of the quaternion, not the unit orientation vector. There are other choices to specify the orientation, such as Eulerangles and rotation matrices, but we prefer quaternions for simplicity. The geometric configuration $\mathcal{C}$ for all $N$ filaments can be written as a column vector:
which is a function of time: $\mathcal{C}(t)$. The translational and angular velocity $\mathit{U},\mathbf{\Omega}$ of all filaments can also be written as a column vector:
Similarly we can write the force and torque $\mathit{F},\mathit{T}$ applied on all filaments as a column vector:
The kinematic equation of motion 35 maps $\mathcal{U}$ to $\dot{\mathcal{C}}(t)=\mathrm{\partial}\mathcal{C}/\mathrm{\partial}t$, via a geometric matrix $\mathcal{G}$.
$\mathcal{G}\in {\mathbb{R}}^{7N\times 6N}$ is a block diagonal matrix, with one $3\times 3$ and one $4\times 3$ block for each particle:
$\mathit{I}}^{3$ is the $3\times 3$ identity matrix, same for every particle. Each $\mathit{I}}^{3$ block simply corresponds to the translational motion $\dot{\mathit{x}}}_{j}={\mathit{U}}_{j$ of each particle $j$. Each $\mathbf{\Psi}}_{j}\in {\mathbb{R}}^{4\times 3$ refers to the rotational motion $\dot{\mathit{\theta}}}_{j}={\mathbf{\Psi}}_{j}{\mathbf{\Omega}}_{j$ of each particle, where for each $j$:
Here ${\u03f5}_{ikj}$ is the LeviCivita symbol for crossproduct in 3D space.
The biological filaments we consider mostly have lengths on the nm to $\mathrm{\mu}\mathrm{m}$ scales. At these scales, solvent viscosity dominates and inertia effects can be ignored, which is the socalled Stokes regime where the mobility matrix $\mathcal{M}$ maps the force $\mathcal{F}$ linearly to the velocity $\mathcal{U}$:
$\mathcal{F}$ includes collision force ${\mathcal{F}}_{C}$ between particleparticle and particlecontainer pairs, linker force between particle pairs ${\mathcal{F}}_{L}$ generated by doubly bound crosslinkers, Brownian force on each particle ${\mathcal{F}}_{B}$ generated by thermal fluctuations, and other externally applied forces ${\mathcal{F}}_{E}$ through gravity and electrostatic fields.
In principal, Equation 35 together with Equation 38 can be integrated directly because both $\mathcal{M}$ and $\mathcal{F}$ are functions of the geometry $\mathcal{C}$ and time only. However, this approach is usually impractical, because ${\mathcal{F}}_{C}$ or ${\mathcal{F}}_{L}$ is usually very stiff functions of the geometry. For example, the collision force ${\mathcal{F}}_{C}$ is usually computed by assuming a very stiff pairwise potential between filaments, such as the LennardJones 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 ${\mathcal{F}}_{C}$ in our previous work on Brownian spherocylinders Yan et al., 2019 and rigid spheres in Stokes flow Yan et al., 2020. Instead of computing ${\mathcal{F}}_{C}$ using repulsive potentials, we imposed nonoverlapping constraints on the geometry $\mathcal{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 nonconstraint, i.e., physical components that are independent of the constraints.
For unilateral constraints, we define the grand distance function ${\mathbf{\Phi}}_{u}$ between every pair of particles as a column vector:
where each ${\mathrm{\Phi}}_{u,{P}_{j}{Q}_{j}}$ is the minimal distance between particles with indices ${P}_{j}$ and ${Q}_{j}$. Similarly, we define the grand distance function ${\mathbf{\Phi}}_{b}$ for bilateral constraints:
where each ${\mathrm{\Phi}}_{b,{P}_{j}{Q}_{j}}$ is the distance between two fixed points on particles ${P}_{j}$ and ${Q}_{j}$, respectively. Physically, ${\mathrm{\Phi}}_{b,{P}_{j},{Q}_{j}}$ is simply the length of each doubly bound crosslinker. With this definition, there are in total ${N}_{u}$ unilateral and ${N}_{b}$ bilateral constraints in the system. In other words, there are in total ${N}_{u}$ possibly colliding pairs of filaments and ${N}_{b}$ doubly bound crosslinkers. Both kinds of constraints are functions of the system geometry, so we shall write them as ${\mathbf{\Phi}}_{b}(\mathcal{C})$ and ${\mathbf{\Phi}}_{u}(\mathcal{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:
For each ${\mathrm{\Phi}}_{u,{P}_{j}{Q}_{j}}$ or ${\mathrm{\Phi}}_{b,{P}_{j}{Q}_{j}}$, there is a corresponding force magnitude ${\gamma}_{u,j}$ or ${\gamma}_{b,j}$, the (normalized) direction vector ${\widehat{\mathit{e}}}_{{P}_{j}}={\widehat{\mathit{e}}}_{{Q}_{j}}$ of this force, and the location ${\mathit{y}}_{{P}_{j}}$ and ${\mathit{y}}_{{Q}_{j}}$ where this force is applied on the filament ${P}_{j}$ and ${Q}_{j}$ respectively, as shown in Appendix 3—figure 1. With norm vectors defined in this way, ${\gamma}_{u}$ or ${\gamma}_{b}$ is positive when the force is repulsive between two filaments.
For unilateral constraints ${\mathbf{\Phi}}_{u}$ and ${\mathit{\gamma}}_{u}$ satisfy this complementarity condition:
This condition means ${\mathbf{\Phi}}_{u}(\mathcal{C})$ and ${\mathit{\gamma}}_{u}$ are orthogonal to each other, and all components of ${\mathbf{\Phi}}_{u}(\mathcal{C})$ and ${\mathit{\gamma}}_{u}$ are nonnegative Yan et al., 2019.
For bilateral constraints ${\mathbf{\Phi}}_{b}$ and ${\mathit{\gamma}}_{b}$ satisfy this linear equality condition because they are modeled as Hookean springs:
$\mathcal{K}\in {\mathbb{R}}^{{N}_{b}\times {N}_{b}}$ is a diagonal matrix, with the stiffness constant κ for each spring on its diagonal $[{\kappa}_{1},{\kappa}_{2},\mathrm{\dots}]$. Obviously every constant ${\kappa}_{j}$ is positive. ${\mathbf{\Phi}}_{b}(\mathcal{C})$ and ${\mathbf{\Phi}}_{b}^{0}$ represent the current and free length of every spring.
Both unilateral and bilateral constraints change over time, as particles move and springs attach to and detach from particles.
All combined together, we reach the equation of motion with geometric constraints:
These equations are solvable when closed by a geometric relation, which maps the force magnitude ${\mathit{\gamma}}_{u}$ and ${\mathit{\gamma}}_{b}$ to the force vectors ${\mathcal{F}}_{u}$ and ${\mathcal{F}}_{b}$:
where ${\mathcal{D}}_{u}$ and ${\mathcal{D}}_{b}$ are sparse matrices containing all orientation vectors of unilateral and bilateral forces, i.e., all $\widehat{\mathit{e}}$ vectors as shown in Appendix 3—figure 1. More details about the definition of $\mathcal{D}$ can be found in the following.
Both ${\mathcal{D}}_{u}$ and ${\mathcal{D}}_{b}$ depend only on the geometry norm vectors $\widehat{{\mathit{e}}_{{P}_{j}}},{\widehat{\mathit{e}}}_{{Q}_{j}}$ and location of constraints ${\mathit{y}}_{{P}_{j}},{\mathit{y}}_{{Q}_{j}}$, together with the particle indices ${P}_{j},{Q}_{j}$, 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 ${\mathcal{D}}_{u}$ and ${\mathcal{D}}_{b}$. For example, if a particle $P$ collides with a fixed substrate, we only include $\widehat{{\mathit{e}}_{P}}$ and ${\mathit{y}}_{P}$ in ${\mathcal{D}}_{u}$, because this substrate does not appear in the mobility matrix $\mathcal{M}$.
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 $\mathrm{\Delta}t$ is limited to be tiny by the temporal stiffness of collision and doubly bound crosslinkers.
The implicit scheme is linearized to avoid expensive largescale nonlinear problems.
With timestep $\mathrm{\Delta}t=h$, Equations 45 and 46 are discretized at timestep $k$ as:
The unknowns to be solved at every timesteps are the constraint force magnitude ${\mathit{\gamma}}_{u}^{k},{\mathit{\gamma}}_{b}^{k}$. Equations 47d and e are nonlinear because ${\mathbf{\Phi}}_{u}^{k+1}$ and ${\mathbf{\Phi}}_{b}^{k+1}$ are nonlinear functions of ${\mathcal{C}}^{k+1}$. Therefore, we linearize these two terms:
Here we have also rewritten the Equation 47 e into a equivalent form, similar to Equation 47 d. The right side, ${\mathit{\gamma}}_{u}\ge 0$ and ${\mathit{\gamma}}_{b}\in \mathbb{R}$ should be understood in the componentwise sense. Equations 48 a and b are now closed and ${\mathit{\gamma}}_{u}^{k},{\mathit{\gamma}}_{b}^{k}$ 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 blockmatrix form:
where the blocks are clear from equation (48)
Here we used the fact that:
The first relation has been well known in the problem of collision constraints Anitescu et al., 1996. In this work we extend this result to bilateral constraints. A proof of this is detailed in Section Symmetry of the geometrically constrained optimization problem.
This formulation means that the coefficient matrix is SymmetricPositiveSemiDefinite (SPSD), because the mobility matrix $\mathcal{M}$ is SPD and $\frac{1}{h}{\mathcal{K}}^{1}$ is positive & diagonal:
Because of this SPSD property, solving Equations 48 is equivalent to solving a constrained quadratic programming (CQP) due to the KarushKuhnTucker condition Nocedal and Wright, 2006:
Here $\mathit{\gamma}=[{\mathit{\gamma}}_{u},{\mathit{\gamma}}_{b}]\in {\mathbb{R}}^{{N}_{u}+{N}_{b}}$ is a column vector, and
This can be conveniently understood as following. $\mathit{q}$ represent the current values of the constraint functions $\mathbf{\Phi}$ plus the (linearized) changes due to nonconstraint forces $\mathcal{F}$, such as Brownian fluctuations. $\mathit{M}$ represent the linearized relation between the unknown constraint force $\mathit{\gamma}$ and the changes of the constraint functions $\mathbf{\Phi}$.
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 BarzilaiBorwein 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 ${\mathcal{K}}^{1}\to 0$ the quadratic term matrix $\mathit{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 noncompliant 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 $\mathrm{\Delta}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 $\mathrm{\Delta}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 isotropicnematic 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 $\mathit{x}$ in the lab frame and its orientation as a unit quaternion $\mathit{\theta}=[s,\mathit{p}]\in {\mathbb{R}}^{4}$. For an arbitrary 3D vector $\mathit{Y}$ which is attached to a particle and follows the particle’s motion, its image $\mathit{y}$ in the lab frame following the particle’s rotation is:
where $\mathit{R}\in {\mathbb{R}}^{3\times 3}$ denotes the rotation matrix generated by the unit quaternion $\mathit{\theta}$.
For both unilateral and bilateral constraints, $\mathcal{D}$ has a sparse column structure:
where ${P}_{i},{Q}_{i}$ 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 ${\mathcal{D}}_{u}$ matrix for collision (unilateral) constraints is:
Because of this structure, to prove Equation 51 we only need to prove the equality ${\mathit{D}}_{PQ}={\nabla}_{\mathcal{C}}{\mathbf{\Phi}}_{PQ}\mathcal{G}$ for a pair of particles $P,Q$, as shown in Appendix 3—figure 1.
We consider two rigid particles centered at ${\mathit{x}}_{P},{\mathit{x}}_{Q}$, each has a point fixed on the body (not necessarily on the surface). ${\mathit{y}}_{P}$ and ${\mathit{y}}_{Q}$ are vectors in the lab frame from the particle centers to the points. The distance between these two points follows the rigid body motion of both particles:
where ${\mathit{R}}_{P}$ and ${\mathit{R}}_{Q}$ are the wellknown rotation matrices. ${\mathit{Y}}_{P}$ and ${\mathit{Y}}_{Q}$ are locations of those two points in their intrinsic coordinate systems. ${\mathrm{\Phi}}_{PQ}=\mathit{r}$ is simply the distance between the two points, dependent on the motion of the two rigid particles.
According to our definition, ${\mathit{D}}_{PQ}$ maps the force magnitude γ between the two particles to force and torque vectors on each particle:
$\left({\nabla}_{\mathcal{C}}\mathbf{\Phi}\right)\mathcal{G}$ can also be explicitly written as follows:
Further, we notice the symmetry of P and Q in the above equations of ${D}_{PQ}$ and ${\nabla}_{\mathcal{C}}{\mathbf{\Phi}}_{PQ}\mathcal{G}$, we only need to prove the following equality for P:
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 $\mathrm{\Phi}$ 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 ${\mathit{Y}}_{P}$ is an arbitrary vector.
The first row of Equation 62 is straightforward because
The second row can be proved as follows. We first derive some general results about quaternions and rotation matrices, dropping the subscript $P$ to simplify equations. When the particle rotates with an angular velocity $\mathit{\omega}$, the motion of $\mathit{y}$ satisfies
$\dot{\mathit{y}}$ can also be directly computed by applying the chain rule on Equation 55, because $\mathit{Y}$ is intrinsic to the particle invariant over time:
The matrix $\mathrm{\Psi}$ bridges angular velocity and quaternion by definition:
We have
This must be valid for arbitrary $\mathit{\omega}$, which is only possible when
Now for another arbitrary vector $\mathit{r}$:
Using Equation 69 we can prove the second row of Equation 62. We first calculate the derivatives of Equation 62 using dummy indices:
Multiply the matrix ${\mathrm{\Psi}}_{kl}$ on both sides:
Substitute the right side by Equation 69, we get:
This is exactly the right side of Equation 62 because by definition ${\mathit{y}}_{P}={\mathit{R}}_{P}{\mathit{Y}}_{P}$ and ${\mathit{e}}_{P}=\mathit{r}/\mathrm{\Phi}$. Therefore Equation 62 holds and the equality Equation 51 holds.
Implementation
As mentioned above, at each timestep we first update the crosslinkers and then the filaments. We implement the two steps in a fully parallelized C ++ codebase, utilizing MPI and OpenMP and scalable to hundreds of CPU cores.
In the crosslinkerupdate step, we have assumed that every crosslinker has bindingunbinding probabilities independent of other crosslinkers. Therefore, it is straightforward to parallelize this step, we only need to search the vicinity of each crosslinker to find the candidate filaments that this crosslinker may bind to. This can be conveniently accomplished by a standard near neighbor detection operation based on bounding volume hierarchy Iwasawa et al., 2016, where the search radius is determined by the maximum stretch of each crosslinker. Once the candidate filaments for each crosslinker have been found, we compute the kMC 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 ${\mathbf{\Phi}}_{b}$ in the constraint solver. If one crosslinker has changed its status from doubly bound to singly bound, the corresponding constraint is removed from ${\mathbf{\Phi}}_{b}$, and vice versa. The geometric matrix ${\mathcal{D}}_{b}$ is also updated according to the current geometry, that is, those locations ${\mathit{y}}_{P},{\mathit{y}}_{Q}$ and norm vectors ${\widehat{\mathit{e}}}_{P},{\widehat{\mathit{e}}}_{Q}$. Then, a near neighbor detection operation is performed for all filaments to determine the unilateral constraints ${\mathbf{\Phi}}_{u}$ and its geometry ${\mathcal{D}}_{u}$. 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 $\mathrm{\Delta}t$. Therefore, we include only close pairs whose minimal distance is below some threshold value $\delta}_{c$. $\delta}_{c$ is controlled by system dynamics, that is, how far each filament may move within each timestep. Empirically, we take ${\delta}_{c}$ to be the diameter of each filament.
Once the constraint problem Equation 48 has been constructed, we run a fully parallel iterative BarzilaiBorwein Projected Gradient Descent (BBPGD) solver Yan et al., 2019 to solve for constraint forces ${\mathit{\gamma}}_{b}$ and ${\mathit{\gamma}}_{u}$, together with the velocities ${\mathcal{U}}_{b}$ and ${\mathcal{U}}_{c}$ due to constraint forces ${\mathcal{F}}_{b}$ and ${\mathcal{F}}_{u}$ by solving the equivalent CQP 53. The cost of every BBPGD iteration scales as $O({N}_{u}+{N}_{b})$, 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 $\nabla f=\mathit{M}\mathit{\gamma}+\mathit{q}$. This sparse matrixvector multiplication spmv is a well optimized standard mathematical operation. The BBPGD solver is implemented using the Trilinos package for distributed linear algebra operations. Once ${\mathcal{U}}_{b}$ and ${\mathcal{U}}_{c}$ have been solved, the filament configuration is updated and then next timestep starts.
Appendix 4
Performance measurements
The bundle contractionbuckling simulation runs on 2 nodes connected by Infiniband, and each node has two AMD EPYC 7742 64Core 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 $S\rightleftharpoons D$ 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.
For the aster formation in bulk problem, each case runs on 1 node of dual Intel Xeon 14core CPUs E52680 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 kinesin5 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 kinesin5 motors become doubly bound and more collisions occur as the asters form. The increase in BBPGD steps dominates because while the dimension of $\mathit{\gamma}$ increases, the dimension of $\mathcal{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 $\mathcal{M}$ when calculating $\nabla f$ in solving Equation 53.
Appendix 5
Aster center analysis of asters formation in bulk
This section provides more details about the simulation in Section Confined filamentmotor protein assemblies of main text.
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
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 $q\approx 0.8\text{}{\mathrm{\mu}\mathrm{m}}^{1}$. This reflects the spacing between individual asters at approximately $1.2\text{}\mathrm{\mu}\mathrm{m}$. This length scale is straightforward to understand. Since we used microtubules of length $0.5\text{}\mathrm{\mu}\mathrm{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 kinesin5 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 ${({L}_{\mathrm{box}}/(2L))}^{3}\approx {10}^{3}$. This simple estimation agrees with our aster center identification results, which on average 1200 aster centers are found for each snapshot of steady state configuration.
The BMT case does not show such a significant special length scale in $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 DensityBased Spatial Clustering of Applications with Noise and is a method to identify clusters from points in space. With a given distance $\u03f5$ and a threshold ${N}_{\mathrm{min}}$ of minimal number of points, DBSCAN searches all clusters such that each cluster has no less than ${N}_{\mathrm{min}}$ points and no point in one cluster is more than distance $\u03f5$ 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 scikitlearn. Once clusters have been identified, we compute the aster centers by averaging the location of all points in each cluster.
We set $\u03f5=100\text{}\mathrm{nm}$, because according to Figure 7 in main text, the minus ends of microtubules are separated roughly $25+53\mathrm{nm}$. We also set ${N}_{\mathrm{min}}=5$.
Identify aster centers by Graph method
The entire microtubulekinesin 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 ${N}_{\mathrm{min}}=5$. We identify aster centers by computing the average location of minus ends of these connected components.
Appendix 6
Confined filamentmotor protein assemblies
This section provides more details about the simulation in Section Confined filamentmotor 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\text{}\mathrm{\mu}\mathrm{m}$ and diameter $25\text{}\mathrm{nm}$ (aspect ratio of $10\text{})$. 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 $300\text{}\mathrm{K}$ and the simulation timestep is ${10}^{4}\text{}\mathrm{s}$, with the system configuration recorded every 500 steps. Solvent viscosity is set at $0.01\text{}\mathrm{pN}\mathrm{\mu}{\mathrm{m}}^{2}\mathrm{s}$. Values for the cylinder diameter, ${D}_{cyl}\in \{0.25,0.75\}\text{}\mathrm{\mu}\mathrm{m}$ are chosen to disrupt the selfassembly of an ideal aster. Initially, microtubules were aligned along the $x$ direction (cylinder axis) such that the initial nematic order parameter, $S=\u27e8\frac{1}{2}(3{\mathrm{cos}}^{2}\theta 1)\u27e9$, was 1. Here, θ is the angle between the microtubule orientation vector and the $+x$ axis, and $\u27e8.\u27e9$ 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=\u27e8\mathrm{cos}\theta \u27e9=0$.
Structural quantification
To measure the structure of our steadystates, we compute the local packing fraction ${\varphi}_{\mathrm{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 ${D}_{cyl}$, and the height is chosen as $25\text{}\mathrm{nm}$. The local packing fraction is computed by calculating the cumulative volume of all microtubules that fall inside each bin, and then dividing by the bin volume. For simplicity, we treat the filaments as cylinders (such that there as no hemispherical caps at their ends). For the local nematic order parameter ${S}_{\mathrm{local}}^{x}(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, ${S}_{\mathrm{local}}^{x}{(x)}_{i}=\frac{1}{2}(3{\mathrm{cos}}^{2}{\theta}_{i}1)$. We weight each ${S}_{\mathrm{local}}^{x}{(x)}_{i}$ by a factor ${W}_{i}(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 doublybound 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 $\rho =\sqrt{{y}^{2}+{z}^{2}}$. Due to the nonperiodic nature in ρ, this pair distirbution function does not decay to 1.
The simulation volume is a cylinder with height $144\text{}\mathrm{\mu}\mathrm{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 polaritysorted 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 steadystate (Appendix 6—figure 1D) shows that plusends (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.
In this case, the simulation volume is a cylinder with height $20\text{}\mathrm{\mu}\mathrm{m}$. Over time, microtubules condense into a bottlebrushlike (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 halfaster. Crosslinking motor proteins are highly concentrated along the central axis of the BB. Here is a kymograph for the local microtubule packing fraction, ${\varphi}_{\mathrm{local}}(x)$ (Appendix 6—figure 2A). We show the evolution of the local nematic order parameter, ${S}_{\mathrm{local}}^{x}(x)$, in Appendix 6—figure 2B. A negative ${S}_{\mathrm{local}}^{x}$ 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\text{}$ 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 plusends (top plot), showing that this state is asterlike, there is significant density present along the X axis. This indicates that there is an accumulation of plusends 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 asterlike. High density at $\rho =0$ suggests that microtubule centers tend to be stacked in $x$.
Ideal bottlebrush state
The ideal bottlebrush 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 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
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 ${\mathrm{\ell}}_{E},{\mathrm{\ell}}_{B}$ geometrically:
With the deformed geometry, the lengths of the two springs are:
When $\alpha \to 0$, the energy $U$ of the two springs can be expanded as:
Here in the first term is simply the linear extension of both springs when α is small. The twospring system generate a equivalent extensional rigidity ${\kappa}_{B}+{\kappa}_{E}$. The second ${\alpha}^{2}$ term governs the bending energy. We can tune the five parameters ${\mathrm{\ell}}_{E}^{0},{\mathrm{\ell}}_{B}^{0},{d}_{B},{\kappa}_{E},{\kappa}_{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 biofilaments such as microtubules, we sometimes assume filaments are inextensible, i.e., ${\kappa}_{E}=\mathrm{\infty}$ and ${\mathrm{\ell}}_{E}={\mathrm{\ell}}_{E}^{0}$. In this special case, the energy of the two springs depends only on $U=\frac{1}{2}{\kappa}_{B}{({\mathrm{\ell}}_{B}{\mathrm{\ell}}_{B}^{0})}^{2}$. By imposing ${\mathrm{\ell}}_{E}={\mathrm{\ell}}_{E}^{0}$, we can solve for $b$:
Then in this case $U$ depends on ${\alpha}^{4}$ in the limit of $\alpha \to 0$. To simplify the notations of the expansion, we define:$s={\mathrm{\ell}}_{B}^{0}{\mathrm{\ell}}_{E}^{0}2{d}_{B}.$
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 $\alpha \to 0$:
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 ${\mathrm{\ell}}_{E},{\mathrm{\ell}}_{B},{d}_{B}$ only. Therefore, the first term $frac12{\kappa}_{B}{s}^{2}$ only ‘shifts’ the zeropoint 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(\alpha )\propto {\alpha}^{4}$. When $s>0$, the leading order terms are nonzero and the energy is asymptotically a quadratic function of α: $U(\alpha )\propto {\alpha}^{2}$.
Special case 2 If we further assume that ${\kappa}_{E}=\mathrm{\infty}$ and ${\mathrm{\ell}}_{E}={\mathrm{\ell}}_{E}^{0}=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:
The expansion of $U$ as $\alpha \to 0$ is also further simplifed:
Here we have the same conclusion as the previous special case, that the dependence of $U$ on α can be tuned between ${\alpha}^{4}$ and ${\alpha}^{2}$ by choosing a proper value of $s$.
Method 2: use bilateral constraints
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 ${\mathit{p}}_{1}$ and ${\mathit{p}}_{2}$ form a plane. This plane is orthogonal to a unit norm vector
For most relevant biological filaments, the bending rigidity is isotropic along different directions on a crosssection of the filament. In other words, the recovering torque is always colinear with the vector ${\mathit{p}}_{1}\times {\mathit{p}}_{2}$ and the recovering deformation is always in plane spanned by ${\mathit{p}}_{1}$ and ${\mathit{p}}_{2}$. This is important because we can simplify the deformation to inplane 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, ${\mathit{p}}_{1}\cdot {\mathit{p}}_{2}$.
Case 1:
When the angle α between ${\mathit{p}}_{1}$ and ${\mathit{p}}_{2}$ is small, we have:
Case 2:
In this form when $\alpha \to 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 ${\mathit{\omega}}_{1}$ and ${\mathit{\omega}}_{2}$ have arbitrary directions, and to the first order of $\mathrm{\Delta}t$ the orientation vectors ${\mathit{p}}_{1}$ and ${\mathit{p}}_{2}$ rotates within $\mathrm{\Delta}t$:
Then, the bending energy after this rotation is:
where we have utilized the vector triple product identity:
This means, to the first order of $\mathrm{\Delta}t$ only the component of rotation ${\mathit{\omega}}_{1}$ and ${\mathit{\omega}}_{2}$ that is inside this plane spanned by ${\mathit{p}}_{1},{\mathit{p}}_{2}$ affect the bending energy. Therefore, to the first order of $\mathrm{\Delta}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:
The rotational mobility matrix for these two rods is:
where ${M}_{1}^{R}$ and ${M}_{2}^{R}$ 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 ${U}_{B}$. More specifically:
where the scalar torque ${T}^{n+1}$ is:
where the higher order terms in $\mathrm{\Delta}t$ have been neglected. If the bending energy Case 2 is used, instead of ${T}^{n+1}\propto E{\alpha}^{n+1,3}$ we have ${T}^{n+1}\propto E{\alpha}^{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:
where $K=3E{\alpha}^{n,2}$ and the geometric matrix $\mathcal{D}$ defines the direction of torque on each rod:
The left side of Equation 94 means the motion of filament segments must satisfy the torquedeformation 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 twosegment problem analytically if the constraint optimization problem contains only Equation 94, in the absence of collisions and Hookean springs:
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\to \mathrm{\infty}$ as the bending rigidity modulus $E$ increases to infinity. In this case, $1/K\to 0$ and the above solution is still stable, and is simplified to $({\omega}_{2}{\omega}_{1})\mathrm{\Delta}t=\frac{1}{6}\alpha $. 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 opensource modeling software is hosted at GitHub.
References

A mechanochemical model explains interactions between cortical microtubules in plantsBiophysical Journal 99:1082–1090.https://doi.org/10.1016/j.bpj.2010.05.037

Formulating ThreeDimensional Contact Dynamics ProblemsMechanics of Structures and Machines 24:405–437.https://doi.org/10.1080/08905459608905271

A theory that predicts behaviors of disordered cytoskeletal networksMolecular Systems Biology 13:941.https://doi.org/10.15252/msb.20177796

Organelle positioning and cell polarityNature Reviews. Molecular Cell Biology 9:874–886.https://doi.org/10.1038/nrm2524

Prime movers: the mechanochemistry of mitotic kinesinsNature Reviews. Molecular Cell Biology 15:257–271.https://doi.org/10.1038/nrm3768

Orientational order of motile defects in active nematicsNature Materials 14:1110–1115.https://doi.org/10.1038/nmat4387

Brownian dynamics of confined rigid bodiesThe Journal of Chemical Physics 143:144107.https://doi.org/10.1063/1.4932062

Intrinsic rates and activation free energies from singlemolecule pulling experimentsPhysical Review Letters 96:108101.https://doi.org/10.1103/PhysRevLett.96.108101

Dynamic strength of molecular adhesion bondsBiophysical Journal 72:1541–1555.https://doi.org/10.1016/S00063495(97)788027

Selfstraining of actively crosslinked microtubule networksNature Physics 15:1295–1300.https://doi.org/10.1038/s4156701906421

Multiscale modeling and simulation of microtubulemotorprotein assembliesPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 92:062709.https://doi.org/10.1103/PhysRevE.92.062709

Multiscale polar theory of microtubule and motorprotein assembliesPhysical Review Letters 114:048101.https://doi.org/10.1103/PhysRevLett.114.048101

Brownian dynamics simulations of model hardsphere suspensionsJournal of NonNewtonian Fluid Mechanics 46:1–28.https://doi.org/10.1016/03770257(93)80001R

BookMechanics of Motor Proteins and the CytoskeletonSunderland: Sinauer associates.

Implementation and performance of FDPS: a framework for developing parallel particle simulation codesPublications of the Astronomical Society of Japan 68:54.https://doi.org/10.1093/pasj/psw053

Theory of Cytoskeletal Reorganization during CrossLinkerMediated Mitotic Spindle AssemblyBiophysical Journal 116:1719–1731.https://doi.org/10.1016/j.bpj.2019.03.013

Comparison of explicit and meanfield models of cytoskeletal filaments with crosslinking motorsThe European Physical Journal. E, Soft Matter 44:45.https://doi.org/10.1140/epje/s10189021000429

Beyond polymer polarity: how the cytoskeleton builds a polarized cellNature Reviews Molecular Cell Biology 9:860–873.https://doi.org/10.1038/nrm2522

Integralbased spectral method for inextensible slender fibers in Stokes flowPhysical Review Fluids 6:014102.https://doi.org/10.1103/PhysRevFluids.6.014102

MitosisCold Spring Harbor Perspectives in Biology 8:a023218.https://doi.org/10.1101/cshperspect.a023218

A fast platform for simulating semiflexible fiber suspensions applied to cell mechanicsJournal of Computational Physics 329:173–209.https://doi.org/10.1016/j.jcp.2016.10.026

Collective Langevin dynamics of flexible cytoskeletal fibersNew Journal of Physics 9:427.https://doi.org/10.1088/13672630/9/11/427

Active matter at the interface between materials science and cell biologyNature Reviews Materials 2:1–14.https://doi.org/10.1038/natrevmats.2017.48

Mechanisms of Mitotic Spindle AssemblyAnnual Review of Biochemistry 85:659–683.https://doi.org/10.1146/annurevbiochem060815014528

Molecular Mechanism of CytokinesisAnnual Review of Biochemistry 88:661–689.https://doi.org/10.1146/annurevbiochem062917012530

MEDYAN: Mechanochemical Simulations of Contraction and Polarity Alignment in Actomyosin NetworksPLOS Computational Biology 12:e1004877.https://doi.org/10.1371/journal.pcbi.1004877

Multimotor transport in a system of active and inactive kinesin1 motorsBiophysical Journal 107:365–372.https://doi.org/10.1016/j.bpj.2014.06.014

The Dynamics of Microtubule/MotorProtein Assemblies in Biology and PhysicsAnnual Review of Fluid Mechanics 48:487–506.https://doi.org/10.1146/annurevfluid010814013639

Wrinkling Instability in 3D Active NematicsNano Letters 20:6281–6288.https://doi.org/10.1021/acs.nanolett.0c01546

A compliant viscoplastic particle contact model based on differential variational inequalitiesInternational Journal of NonLinear Mechanics 53:2–12.https://doi.org/10.1016/j.ijnonlinmec.2013.01.010

ConfinementInduced SelfPumping in 3D Active FluidsPhysical Review Letters 125:268003.https://doi.org/10.1103/PhysRevLett.125.268003

The load dependence of rate constantsThe Journal of Chemical Physics 128:215101.https://doi.org/10.1063/1.2920475

Computing collision stress in assemblies of active spherocylinders: Applications of a fast and generic geometric methodThe Journal of Chemical Physics 150:064109.https://doi.org/10.1063/1.5080433

A scalable computational platform for particulate Stokes suspensionsJournal of Computational Physics 416:109524.https://doi.org/10.1016/j.jcp.2020.109524

Energetics and dynamics of constrained actin filament bundlingBiophysical Journal 90:4295–4304.https://doi.org/10.1529/biophysj.105.076968
Decision letter

Pierre SensReviewing Editor; Institut Curie, CNRS UMR168, France

Anna AkhmanovaSenior Editor; Utrecht University, Netherlands
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Decision letter after peer review:
Thank you for submitting your article "aLENS: towards the cellularscale simulation of motordriven cytoskeletal assemblies" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Anna Akhmanova as the Senior Editor. The reviewers have opted to remain anonymous.
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential revisions:
As you will see from the reports below, the referees were favourably impressed by the performance on your numerical method, which allows to stimulate a large number of filaments and motors in a reasonable time. This paper is seen as the presentation of a new tool illustrated by examples, rather than the presentation of new scientific results. The essential revisions include:
1. Software availability (all referees) To be of value, the software should be freely available and useable. This does not seem to be the case at present.
2. Benchmarking (referees 1 and 3) It is important to provide some level of validation of the algorithm, which is currently missing. This can be done by comparing some aspects of the numerical results to existing analytical models. Regarding the binding and motor activity of cross linker, some simple toy model (e.g. gliding assay for motors, bound fraction for crosslinkers) would already be very informative, but more complex properties could also be considered as in LeraRamirez and Nedelec , Cytoskeleton 2019 (Theory of antiparallel microtubule overlap stabilization by motors and diffusible crosslinkers). Regarding collective effect, validation of the results could follow approaches such as the ones employed in Gao etal. PRE 2015 (Multiscale modeling and simulation of microtubulemotorprotein assemblies). Performance benchmark compared to other available algorithm should also be discussed.
3. Modularity (referees 2 and 3). At present, only rigid filaments interacting with one type of crosslinker are presented. The extension flexible filaments and different types of interacting proteins are discussed as possibilities but are not implemented at the moment. Although this can be seen as the natural next step, this modularity is essential for the algorithm to be useful to the community.
Reviewer #1:
The study by Yan et al., developed a novel computational framework for modelling cytoskeletal cellular processes that allows for further investigations into the material properties associated with such processes. The computational methodology involves modelling cytoskeletal filaments as rigid spherocylinders while the Hookean law is employed to model crosslinkers. The computational algorithm, aLENS performs three key tasks in a sequential manner that lends itself naturally to high performance computing. To demonstrate the applicability and usefulness of aLENS, the authors present (i) selfaligning and buckling networks for a significantly large number of filaments than previously studied, and (ii) the interplay between polarity of motor walking and polarity of filaments, which seems to suggest that the ability of motors to continuously walk without endpausing is crucial to effective polarity sorting. The authors also investigated the formation of asters, which seem to form when crosslinking motors reorganise filaments so that their minus ends are clustered and held tight by paused motors.
Strengths
1. The development of an alternative novel computational framework that offers far more flexibility, applicability and is scalable across multiscales.
2. The methodology overcomes the timescale limitations imposed by conventional explicit timestepping methods that are key to modelling the dynamics of the cytoskeletal filaments and motors.
3. aLENS utilizes efficiently highperformance parallel computing resources to scale to cellular scale systems.
Weakness
1. The lack of benchmarking that is associated with algorithm comparisons for performance and robustness.
2. The lack of rigorous validation of the algorithm against suitably identified grounds truths.
3. No clear demonstration of the efficiency of the algorithm and how it compares to current conventional algorithms of this nature.
4. There are no clear comparisons between predictions of the computational algorithm and experimental data or observations.
My recommendations to authors are as follows:
(i) To consider rigorous validation of the computational methodology either by using synthetic or experimental data.
(ii) To demonstrate computational efficiency and robustness of the algorithm by comparing results to current conventional methodologies.
(iii) To demonstrate efficiency and accuracy when aLENS is compared to current conventional methodologies.
Reviewer #2:
In this article, the authors present a new method for modeling cytoskeletal networks inside cells. In particular they model the interactions between filaments (microtubules in this case but the method is not limited to this) and motors. The main contribution appears to be in the efficiency of this method. The main challenge in this area of numerical research is handling steric interactions between the large numbers of interacting cytoskeletal filaments. Essentially, hard repulsion methods require very small timescales that can make reasonable simulations infeasible while soft repulsion approximations of these can lead to numerical issues. Their method takes a mathematically different approach to addressing this issue that may substantially speed these simulations allowing users to simulate more realistically sized systems consisting of up to 10^{6}10^{7} filaments.
This is a technically strong and well written article addressing a problem of significant importance. I will note that it is difficult for me assess the value of its content to the field without seeing how this is implemented. The main value of this method is that it can do what other methods already do, but faster. My understanding is that much of the method is previously published and has been repackaged for use in this domain. In that sense, it is a tool. That is only of value if it is deployed in a manner that is well structured, documented, and at least somewhat usable and modifiable by technically capable users. Since the github link has not been included, I cannot assess this at the moment. Given that this is life sciences journal rather than a numerical methods journal, I think for publication at eLife it is critical that this be more than a methods article.
Without seeing this, I have a few questions. The computationally intensive elements of this method are in C++, but what about the more model specification oriented elements. Is this fully C++ or are you interfacing with a higher level language? Second, is this essentially an internal use code, or are you attempting to make this at least somewhat usable to other researchers. I strongly recommend the latter otherwise this will not really be an advancement over approaches such as Cytosim or MEDYAN. Faster but less usable is a losing combination.
I think this approach has strong potential, provided this issue is addressed. Below I'll note a few specific comments.
Specific comments
Line 74 – Here it is mentioned that all filaments will be considered to be strait and rigid but that this can be relaxed by jointing smaller segments together. While these filaments may have long persistence lengths relative to the cell size, their bending capacity is still vitally important to their dynamics and so if this is simple to implement, the authors should include this capability in the origination of this method.
Equation 3 – Please specify here what [s,p] are. As far as I can tell, these are not even fully defined in the SM. Also, is there really value in calling this a quaternion? Are you ever using either the geometric or algebraic properties of these? If not, I would suggest just defining it as a set of orientation variables. Also, for the general reader, it would be useful to describe why the orientation variable here is in R4 rather than R3.
SM Figures 4 / 5 – Can you clarify what you are quantifying for speed? Is this wall time / timestep or is it wall time / simulated second.
Reviewer #3:
In this article, the authors describe a new software for simulating cytoskeleton assemblies, and provide exiting examples of the software capabilities. The software offers the possibility to treat steric interactions as constraints rather than as stiff potentials, which should greatly improve performance. Thermodynamics of protein binding also seem carefully implemented.
While the software ships with great features, and fulfill a possible need in the field, there are several issues to take in consideration:
Software issues:
1 (No link to the software is provided in the article. aLens nonetheless is available on GitHub after a quick search, so this review is based on the current GitHub version)
My first concern is the difficulty to compile aLens. Installing aLens requires running an external script to alter one's environment (or significantly altering the compile file CMakeList.txt). Therefore, I was not able to compile this program, being unwilling to run a script to change my environment. Compiling of the software should be made much easier to be used by a wider, biologyoriented audience.
2 Modularity:
While the authors mention the cross linking model to have a modular design, the software itself does not appear to be modular. There is currently no possibility of having different types of proteins (or different type of filaments). This is a very strong limitation for its general use.
Scientific issues:
1 – Currently, aLens can only simulate straight filaments. The authors claim that aLens could easily as chain of short jointed segments. I would not think this easy, except possibly for the authors. However, many studies showed the paramount importance of bending rigidity in the mechanics of networks. Notably, Lenz and Gardel showed that because of this flexibility, networks tended to be globally contractile rather than extensile (filaments can buckle but not stretch). Therefore, aLens currently has access to only a manifold of the phase diagram, that may not be relevant e.g. for example 1, figure 2. Moreover, the filaments are not dynamic.
2 – While very interesting examples are provided, the authors do not provide a verification of their algorithm and implementation. While simple examples such as buckling cannot be addressed, maybe collective effects could be used as theoretical benchmarks? These could include collective effects of motors and nematic ordering, for which there exist theoretical results.
The authors therefore found success in simulating large systems with adequate thermodynamics. This seems like a notable technical advance. They were successful in providing impressive examples derived from the software. However, the lack of benchmarking (theoretical and performancewise) mean that is not currently possible to assess exactly how successful the authors were.
The lack of filament flexibility and dynamics, as well as the impossibility to use more than one type of proteins will drastically limit its use by the biological community. The lack of code modularity will limit the possibility of external developers participating.
Overall, it seems that this is a very powerful approach, but the article needs to find a message. It could be emphasizing either:
– the implementation, in which case theoretical, and maybe performance benchmarks should be provided.
– the software package itself, in which case the software should be made more usable by the community, and possibly encompass a more general family of problems (or make them at least possible to implement via code modularity).
– the scientific results: in which case theoretical benchmark should be provided for simpler systems as a verification, and more complicated results should be put in their scientific context.
Currently, while very impressive technically, the software and associated article seem to fall short in each category.
https://doi.org/10.7554/eLife.74160.sa1Author response
Essential revisions:
As you will see from the reports below, the referees were favourably impressed by the performance on your numerical method, which allows to stimulate a large number of filaments and motors in a reasonable time. This paper is seen as the presentation of a new tool illustrated by examples, rather than the presentation of new scientific results. The essential revisions include:
1. Software availability (all referees) To be of value, the software should be freely available and useable. This does not seem to be the case at present.
We agree that it is important that aLENS be freely available and useable. The source code has been, and is, available on GitHub (https://github.com/flatironinstitute/aLENS) but this was not made sufficiently clear in the original submission. Now we have added a footnote to the last section about software availability. On GitHub detailed instructions have been provided for compiling the source and, as an aid to users, we have now added an automated script in the source code to install dependency libraries. We have also packaged the precompiled executable as a Docker image, freely available on Docker Hub (https://hub.docker.com/r/wenyan4work/alens), so that any interested users can run aLENS on any operating system, without compilation.
In sum, in this revision we have provided the full information of how to obtain and compile the source code, with additional resources to use the program in more userfriendly ways, e.g., precompiled executables.
2. Benchmarking (referees 1 and 3) It is important to provide some level of validation of the algorithm, which is currently missing. This can be done by comparing some aspects of the numerical results to existing analytical models. Regarding the binding and motor activity of cross linker, some simple toy model (e.g. gliding assay for motors, bound fraction for crosslinkers) would already be very informative, but more complex properties could also be considered as in LeraRamirez and Nedelec , Cytoskeleton 2019 (Theory of antiparallel microtubule overlap stabilization by motors and diffusible crosslinkers). Regarding collective effect, validation of the results could follow approaches such as the ones employed in Gao etal. PRE 2015 (Multiscale modeling and simulation of microtubulemotorprotein assemblies). Performance benchmark compared to other available algorithm should also be discussed.
We agree with the overall comment. To address this issue, we have added a new section Verification and Benchmarking to the manuscript, providing details about verification and performance benchmarks. On verification, one of the few available groundtruth cases, centrally relevant to this oftendense system, is the isotropicnematic transition of rigid Brownian liquid crystals. We have shown in recent work [1] that we recover quantitatively the equation of state and isotropicnematic phase transition for this system using the core steric interaction solver of aLENS. Of course, this example does not include motors. To address this, we chose two experimental systems that focus on different aspects of our method.
Reviewer 3 believed that aLENS did not have the capability to handle mixed motor types. This is not the case. Thus, as our first validation we compare with the experiments in [2], focusing on a simulated gliding assay in which directed transport of microtubules is controlled by mixtures of active and inactive motors. Our new Figure 2 in the revision shows that our simulation results are in excellent quantitative agreement with experiments over the entire range of active motor number fraction, while also demonstrating the capability of aLENS in handling multiple motor types. Further, our current aLENS implementation also supports simulation of arbitrary mixtures of filaments with different radii and lengths as well as spherical particles.
We next compare with the experiments in [3], focusing on the observed large scale collective behavior there. They show that in densely packed and highly crosslinked microtubuleXCTK2 assemblies (XCTK2 is a minusend directed kinesin motor), microtubules show nearly constant motorinduced speeds with those speeds nearly independent of motor concentration and local mean polarity of the filament assemblage. This is a nontrivial phenomena that occurs in densely packed and highly crosslinked assemblies, and has also been observed in extract spindles [4]. Our new Figure 3 in the revision shows that our simulation results agree with these experiments. This confirms the theoretical interpretation given in [3], based on continuum modeling, that these highly aligned microtubule/motor systems can be seen as two interpenetrating active gels being pulled past each other by motors.
For verification, we chose to directly compare with experimental data, rather than with other simulation work, because we focus on the predictive power of our approach, and prefer to avoid debate about technical detail differences between our methods and other simulation tools. To this end, we show that our simulation methods and modeling can directly match experimental observations. We note that for all comparisons we chose motor parameters based on experimental estimates for that particular motor type, or for other members in the same motor family if the data for that particular motor is not available. In particular, we did not perform any parameter fitting process in these benchmark simulations, which nonetheless accurately reproduce the experimental measurements. These two direct comparisons with experiment cover multiple aspects of our method and implementation, including mixed motor species, longtime directed motion, densely packed assemblies, and largescale collective interactions between motors and filaments. We believe that this provides convincing verification of the simulation methods.
For the referees’ interest, we note that we are preparing a new paper where we combined aLENS simulation results with analytic modeling to reveal a selfregulating mechanism for the “selfstraining motion” of densely crosslinked nematic bundles of microtubules [3]. This lays the microscopic foundation of the continuum theory, and its generalizations, for such systems [3].
In the revision, we have also added a large scale parallelization efficiency test case to show the major advantage of aLENS. For a system with 4 million objects (1 million microtubules and 3 million motors) and 8 million constraints, we achieve almost ideal linear speedup with increasing resources up to 1536 cores, at which point one timestep takes less than 1 second wall clock time. Such largescale parallel performance has never been demonstrated by other methods in this field, to the best of our knowledge.
3. Modularity (referees 2 and 3). At present, only rigid filaments interacting with one type of crosslinker are presented. The extension flexible filaments and different types of interacting proteins are discussed as possibilities but are not implemented at the moment. Although this can be seen as the natural next step, this modularity is essential for the algorithm to be useful to the community.
Our current implementation supports an unlimited number of different species of motors, and we have provided one such example with two species of motors as the first comparative example in the new section Verification and Benchmark of the revision.
Regarding the incorporation of flexible filaments, we have added a section in the supplemental material covering two different methods for implementing bending rigidity. One method requires no modification of the current codebase as it involves only additional Hookean springs. The other incorporates the bending rigidity as a set of equality constraints, solved simultaneously with collision and doubly bound motors, within our constraintoptimization solver.
Our codebase is actually modular. It follows a standard objectoriented design, so that users can conveniently swap our motor model for their own models, or add more features including additional constraints (as above for bending rigidity), using standard C++ programming techniques for high performance. In addition to the methodology derivations in the supplemental material describing the motor model, we have also included a detailed document, along with the source code distributed via GitHub, describing the codebase implementation of the motor protein model Protein/ProteinModel.md. Our GitHub documentation provides a clear roadmap for developing additional modules.
Further development of aLENS does require some C++ programming knowledge but such expertise is widely found in academia and is certainly not limited to our group. Indeed, we have already witnessed the development of aLENS extensions in ongoing development work with other collaborating groups.
Reviewer #1 (Recommendations for the authors):
The study by Yan et al., developed a novel computational framework for modelling cytoskeletal cellular processes that allows for further investigations into the material properties associated with such processes. The computational methodology involves modelling cytoskeletal filaments as rigid spherocylinders while the Hookean law is employed to model crosslinkers. The computational algorithm, aLENS performs three key tasks in a sequential manner that lends itself naturally to high performance computing. To demonstrate the applicability and usefulness of aLENS, the authors present (i) selfaligning and buckling networks for a significantly large number of filaments than previously studied, and (ii) the interplay between polarity of motor walking and polarity of filaments, which seems to suggest that the ability of motors to continuously walk without endpausing is crucial to effective polarity sorting. The authors also investigated the formation of asters, which seem to form when crosslinking motors reorganise filaments so that their minus ends are clustered and held tight by paused motors.
Strengths
1. The development of an alternative novel computational framework that offers far more flexibility, applicability and is scalable across multiscales.
2. The methodology overcomes the timescale limitations imposed by conventional explicit timestepping methods that are key to modelling the dynamics of the cytoskeletal filaments and motors.
3. aLENS utilizes efficiently highperformance parallel computing resources to scale to cellular scale systems.
We thank the reviewer for the positive feedback.
Weakness
1. The lack of benchmarking that is associated with algorithm comparisons for performance and robustness.
2. The lack of rigorous validation of the algorithm against suitably identified grounds truths.
3. No clear demonstration of the efficiency of the algorithm and how it compares to current conventional algorithms of this nature.
4. There are no clear comparisons between predictions of the computational algorithm and experimental data or observations.
My recommendations to authors are as follows:
(i) To consider rigorous validation of the computational methodology either by using synthetic or experimental data.
(ii) To demonstrate computational efficiency and robustness of the algorithm by comparing results to current conventional methodologies.
(iii) To demonstrate efficiency and accuracy when aLENS is compared to current conventional methodologies.
We appreciate the points the reviewer is making. As the reviewer can see, we have addressed many of these issues in our response to the editors above. Here we respond in more detail to the reviewer’s comment regarding comparison with conventional algorithms. One main focus of aLENS design is to guarantee the temporal stability for systems with steric interactions and the elastic forces arising from doubly bound crosslinking motors. On the former, as a good point of comparison we will compare aLENS with the use of steep and stiff WCA potentials to prevent geometric overlaps. As an aside, we note that steric interactions are often either entirely neglected (as in AFINES [5]) or approximated by soft potentials (as in Cytosim [6]).
First of all, the original unmodified WCA potential is unstable for the smallest timestep we tried ∆t = 4.43 × 10^{−7} s, and none of the simulations can generate meaningful data, unless in the limit of zero volume fraction where the microtubules do not touch each other.
The WCA potential must be softened for practical simulations. The strategy [8] is to linearize the WCA potential for r < r_{c}, where r_{c} is a cutoff distance, but maintain the continuity of the modified force at r_{c}. This linearization strategy generates a constant repulsive force when the filaments are closer than r_{c}. Then, the timestep size must be small enough such that the chance of a pair getting closer than r_{c} is small. With this strategy, the softened WCA potential can reproduce the equation of state for the entire range of isotropicnematic transition.
In aLENS we keep the order of accuracy at first order in timestep, comparable to the conventional explicit timestepping methods using WCA potentials, but decouple the stability requirement from the accuracy requirement. The constraint solver of aLENS allows stable treatment of both collision and crosslinker forces in an unconditionally stable method.
This means if the accuracy constraint is sufficiently high then conventional methods will be stable and the two methods will have comparable (firstorder) accuracy. As shown in Author response image 1, aLENS is equally accurate as conventional methods in reproducing the equation of state at the same timestep size, but remains stable if the timestep is increased by 10 − 100×, with about 10% − 20% loss in accuracy, in accordance with its firstorder convergence. For the such systems without crosslinkers, this stability may also be achieved using conventional softened WCA potentials, if the potential is further softened with increasing timesteps. This softening strategy, however, cannot be applied to the Hookean spring forces of doubly bound crosslinkers, without significantly affecting the desired mechanical properties of the crosslinkers. In our previous work [9], we used a similar explicit timestepping method to handle the crosslinkers which limited crosslinker stiffness to ∼ 1pNµm^{−1}, which is significantly softer than realistic crosslinkers.
The linearized optimization problem is constructed such that even in the limit where the crosslinkers become infinitely stiff, i.e., noncompliant clamps or in the limit where an infinite number of crosslinkers connect a pair of filaments, the minimization problem is still convex and the solution remains stable.
This stability guarantee of aLENS is crucial in largescale Brownian dynamics simulations because a large Brownian displacement or for a geometry where many crosslinkers bind two filaments, although the probability is very low, is almost surely generated for systems with tens of thousands of objects. In conventional potential type methods, such situations require very complicated ‘failsafe’ mechanisms such as partially linearized potentials, softcore potentials, force limiters, or displacement limiters. If these ‘failsafe’ mechanisms are not properly finetuned, one unstable event between a single pair of filaments may crash the entire simulation. None of these are necessary in aLENS because aLENS is designed to be stable. When large Brownian jumps or other large stochastic crosslinker forces occur, locally the accuracy is lost for this object at this time, but this local error usually resumes to its ‘normal’ state and statistically the sampled system property is not affected. Therefore, the users of aLENS have the option to either quickly scan the parameter space with large timesteps if larger error is acceptable, or accurately scrutinize a few cases with small timesteps.
One major motivation of developing aLENS is to correctly simulate densely packed and crosslinked filament systems accurately. The editors might consider the possibility that aLENS is also providing the benchmarking computational problems that would be so very useful to this field.
Reviewer #2 (Recommendations for the authors):
In this article, the authors present a new method for modeling cytoskeletal networks inside cells. In particular they model the interactions between filaments (microtubules in this case but the method is not limited to this) and motors. The main contribution appears to be in the efficiency of this method. The main challenge in this area of numerical research is handling steric interactions between the large numbers of interacting cytoskeletal filaments. Essentially, hard repulsion methods require very small timescales that can make reasonable simulations infeasible while soft repulsion approximations of these can lead to numerical issues. Their method takes a mathematically different approach to addressing this issue that may substantially speed these simulations allowing users to simulate more realistically sized systems consisting of up to 10^{6}10^{7} filaments.
This is a technically strong and well written article addressing a problem of significant importance. I will note that it is difficult for me assess the value of its content to the field without seeing how this is implemented. The main value of this method is that it can do what other methods already do, but faster. My understanding is that much of the method is previously published and has been repackaged for use in this domain. In that sense, it is a tool. That is only of value if it is deployed in a manner that is well structured, documented, and at least somewhat usable and modifiable by technically capable users. Since the github link has not been included, I cannot assess this at the moment. Given that this is life sciences journal rather than a numerical methods journal, I think for publication at eLife it is critical that this be more than a methods article.
Without seeing this, I have a few questions. The computationally intensive elements of this method are in C++, but what about the more model specification oriented elements. Is this fully C++ or are you interfacing with a higher level language? Second, is this essentially an internal use code, or are you attempting to make this at least somewhat usable to other researchers. I strongly recommend the latter otherwise this will not really be an advancement over approaches such as Cytosim or MEDYAN. Faster but less usable is a losing combination.
I think this approach has strong potential, provided this issue is addressed. Below I'll note a few specific comments.
We thank the reviewer for this constructive feedback. We have addressed some of these issues above in our response to the editors. We do disagree with the statement that the main value in the new work is that it does what other methods already do, but faster. While it may appear at first glance that the new method is similar to others that are currently available, the choice of model and implementation are significantly different. That allows aLENS to be more accurate.
In particular, we accurately handle the steric collisions and crosslinking interactions in a unified constraint solver while maintaining very strong temporal integration stability. This allows aLENS to be more accurate for treating cytoskeletal systems. For example, in previous work studying rods in the absence of motors we demonstrated that the constraint method can accurately and fully reproduce the equationofstate and isotropicnematic phase transition in dense systems of thermally fluctuating rods. By contrast, in other currently available simulation packages steric interactions are sometimes neglected (in AFINES [5]) or approximated by soft potentials (in Cytosim [6]). These approximations are suitable for relatively dilute systems, but are less accurate for dense packing.
Moreover, our treatment of motor/crosslinker binding and unbinding correctly recapitulates the distribution and chemical kinetics of crosslinking proteins in the equilibrium limit, which is important for quantitative modeling of motor mechanochemistry. Other available packages make use of ad hoc rules for crosslinker binding and unbinding, which for some problems can give physically unrealistic results.
Finally, the method here incorporates Brownian dynamics, which is neglected in some other packages (e.g., MEDYAN).
For these reasons, we believe that in addition to the improved computational efficiency, the range of capabilities that aLENS offers are not present in any other available method.
In the revision we provide further information on the implementation and availability of aLENS. In short, it is fully implemented in C++ to ensure efficiency, following a standard objectoriented design. It is available as open source code on GitHub, and we also provide a precompiled readytorun binary executable through Docker Hub. The design is modular such that, for example, the motor model can be swapped out for other models.
Specific comments
Line 74 – Here it is mentioned that all filaments will be considered to be strait and rigid but that this can be relaxed by jointing smaller segments together. While these filaments may have long persistence lengths relative to the cell size, their bending capacity is still vitally important to their dynamics and so if this is simple to implement, the authors should include this capability in the origination of this method.
We agree with the reviewer that filament bending can be quite important in cytoskeletal systems. Accordingly, we have now included in the revised supplemental material a detailed description of how flexible filaments with bending rigidity can be implemented within the constraint solver in aLENS with two different methods. One of them requires no new additions to the codebase. We can further prove that both methods are still stable even in the limit of infinite bending rigidity, similar to the stability for infinitely stiff doubly bound springs. One of our collaborators is working on the extension and will present the method and applications in a separate publication.
Equation 3 – Please specify here what [s,p] are. As far as I can tell, these are not even fully defined in the SM. Also, is there really value in calling this a quaternion? Are you ever using either the geometric or algebraic properties of these? If not, I would suggest just defining it as a set of orientation variables. Also, for the general reader, it would be useful to describe why the orientation variable here is in R4 rather than R3.
We appreciate the reviewer’s suggestion to clarify this point. In response to this comment, we have clarified the definition in the revised main text. In particular, the use of quaternions is a standard approach in rigid body rotational kinematics. [s,p] are the scalar and vector parts of the quaternion, respectively. Representing orientations using quaternions has several crucial advantages, including its compact memory footprint (4 for a quaternion but 9 for a rotation matrix) and the important singularityfree property (the use of Euler angles may lead to the problem of gimbal locks). We therefore make use of this standard practice in the computational geometry of rigid body kinematics.
SM Figures 4 / 5 – Can you clarify what you are quantifying for speed? Is this wall time / timestep or is it wall time / simulated second.
The units are wall clock time per simulation timestep. We have clarified this in the revised figure captions.
Reviewer #3 (Recommendations for the authors):
In this article, the authors describe a new software for simulating cytoskeleton assemblies, and provide exiting examples of the software capabilities. The software offers the possibility to treat steric interactions as constraints rather than as stiff potentials, which should greatly improve performance. Thermodynamics of protein binding also seem carefully implemented.
While the software ships with great features, and fulfill a possible need in the field, there are several issues to take in consideration:
We thank the reviewer for the constructive comments.
Software issues:
1 (No link to the software is provided in the article. aLens nonetheless is available on GitHub after a quick search, so this review is based on the current GitHub version)
My first concern is the difficulty to compile aLens. Installing aLens requires running an external script to alter one's environment (or significantly altering the compile file CMakeList.txt). Therefore, I was not able to compile this program, being unwilling to run a script to change my environment. Compiling of the software should be made much easier to be used by a wider, biologyoriented audience.
We thank the reviewer for this important comment. We hope the reviewer will be satisfied with our response, which can be found in our reply to the editor’s comments.
2 Modularity:
While the authors mention the cross linking model to have a modular design, the software itself does not appear to be modular. There is currently no possibility of having different types of proteins (or different type of filaments). This is a very strong limitation for its general use.
Thank you again for this comment. And, again, we think the reviewer will be satisfied on this through our response to the editor’s comments. In short, aLENS does have a modular design, and is perfectly capable of simultaneously handling multiple motor types, as we now demonstrate in our new section on verification and validation. Further, our current aLENS implementation also supports simulation of arbitrary mixtures of filaments with different radii and lengths as well as spherical particles.
Scientific issues:
1 – Currently, aLens can only simulate straight filaments. The authors claim that aLens could easily as chain of short jointed segments. I would not think this easy, except possibly for the authors. However, many studies showed the paramount importance of bending rigidity in the mechanics of networks. Notably, Lenz and Gardel showed that because of this flexibility, networks tended to be globally contractile rather than extensile (filaments can buckle but not stretch). Therefore, aLens currently has access to only a manifold of the phase diagram, that may not be relevant e.g. for example 1, figure 2. Moreover, the filaments are not dynamic.
We agree with the reviewer that filament flexibility and polymerization dynamics are important in cytoskeletal systems. Therefore, we have included a new section in the supplemental material to show how flexible filaments with adjustable extensional and bending rigidity can be implemented in aLENS, via two different methods. The first method requires no modifications to the current codebase, while the second method incorporates bending joints as new constraints in the constrained minimization solver but requires some code extensions. Both types of method/proposed implementation resolve collisions, doubly bound motor forces, and bending rigidity in a single unified solver. Further, both of these remain stable in the limit where the extensional and bending rigidity of filaments are infinite.
While we have not yet implemented filament polymerization dynamics, this could be done in a straightforward way since the change of filament length can be completely decoupled from the geometric constraint optimization solver. The only necessary improvement to the codebase will be to address how a crosslinker interacts with a dynamic filament, i.e., how the crosslinker tracks its own binding status and location while walking along a dynamic filament. Although currently not implemented, this feature is in our future roadmap of aLENS.
2 – While very interesting examples are provided, the authors do not provide a verification of their algorithm and implementation. While simple examples such as buckling cannot be addressed, maybe collective effects could be used as theoretical benchmarks? These could include collective effects of motors and nematic ordering, for which there exist theoretical results.
The authors therefore found success in simulating large systems with adequate thermodynamics. This seems like a notable technical advance. They were successful in providing impressive examples derived from the software. However, the lack of benchmarking (theoretical and performancewise) mean that is not currently possible to assess exactly how successful the authors were.
The lack of filament flexibility and dynamics, as well as the impossibility to use more than one type of proteins will drastically limit its use by the biological community. The lack of code modularity will limit the possibility of external developers participating.
Thank you for the comments. We have addressed many of them directly and constructively in our response to the editors. In short, we have added a new Verification and Benchmarks section in the revised manuscript to verify that our simulation results are in good quantitative accord with experimental measurements for both smallscale, single microtubule systems with multiple motor types, and largescale, manyfilament systems that are driven collectively by motorfilament interactions. There you will also find a large scale parallel efficiency test, showing nearly ideal linear efficiency for 4 million objects on 1536 cores. Such scaling behavior and efficiency is not available with other methods.
Overall, it seems that this is a very powerful approach, but the article needs to find a message. It could be emphasizing either:
– the implementation, in which case theoretical, and maybe performance benchmarks should be provided.
– the software package itself, in which case the software should be made more usable by the community, and possibly encompass a more general family of problems (or make them at least possible to implement via code modularity).
– the scientific results: in which case theoretical benchmark should be provided for simpler systems as a verification, and more complicated results should be put in their scientific context.
Currently, while very impressive technically, the software and associated article seem to fall short in each category.
We thank the reviewer for their comments. We found them very helpful. We believe that in the revision we have strengthened the manuscript on all these fronts. We have provided both theoretical and performance benchmarks, improved code accessibility and usability for external users, emphasized its modularity, and used aLENS to solve problems of current scientific import.
One of our goals is to provide the cytoskeletal physics community with a new and transformative software tool which they can apply on their own problems and to which they can contribute new functionalities. That said, using aLENS assumes a certain level of programming expertise but that expertise is certainly not singular to our own group and is widely found in academia.
In summary, aLENS is designed to simulate cytoskeletal systems that are challenging to model with existing software, by enabling efficient and accurate treatment of steric interactions and crosslinking forces, chemical kinetics, and Brownian motion in a highperformance parallel implementation. This gives access to lengthscales and timescales of relevance to cell biology. The design of aLENS also allows continuous method development. New features such as bending, polymerization, and even hydrodynamics can and are being incorporated into aLENS.
References
[1] Wen Yan, Huan Zhang, and Michael J. Shelley. Computing collision stress in assemblies of active spherocylinders: Applications of a fast and generic geometric method. The Journal of Chemical Physics, 150(6):064109, 2019.
[2] Lara Scharrel, Rui Ma, René Schneider, Frank Jülicher, and Stefan Diez. Multimotor Transport in a System of Active and Inactive Kinesin1 Motors. Biophysical Journal, 107(2):365–372, 2014.
[3] Sebastian Fürthauer, Bezia Lemma, Peter J. Foster, Stephanie C. EmsMcClung, CheHang Yu, Claire E. Walczak, Zvonimir Dogic, Daniel J. Needleman, and Michael J. Shelley. Selfstraining of actively crosslinked microtubule networks. Nature Physics, 15(12):1295–1300, 2019.
[4] Daniel J. Needleman, Aaron Groen, Ryoma Ohi, Tom Maresca, Leonid Mirny, and Tim Mitchison. Fast microtubule dynamics in meiotic spindles measured by single molecule imaging: Evidence that the spindle environment does not stabilize microtubules. Molecular Biology of the Cell, 21(2):323–333, 2010. PMID: 19940016.
[5] Simon L. Freedman, Shiladitya Banerjee, Glen M. Hocky, and Aaron R. Dinner. A Versatile Framework for Simulating the Dynamic Mechanical Structure of Cytoskeletal Networks. Biophysical Journal, 113(2):448–460, 2017.
[6] Francois Nedelec and Dietrich Foethke. Collective Langevin dynamics of flexible cytoskeletal fibers. New Journal of Physics, 9(11):427–427, 2007.
[7] Peter Bolhuis and Daan Frenkel. Tracing the phase boundaries of hard spherocylinders. The Journal of Chemical Physics, 106(2):666–687, 1997.
[8] Robert Blackwell, Oliver SweezySchindler, Christopher Baldwin, Loren E. Hough, Matthew A. Glaser, and M. D. Betterton. Microscopic origins of anisotropic active stress in motordriven nematic liquid crystals. Soft Matter, 12(10):2676–2687, 2016.
[9] Peter J. Foster, Wen Yan, Sebastian Fürthauer, Michael J. Shelley, and Daniel J. Needleman. Connecting macroscopic dynamics with microscopic properties in active microtubule network contraction. New Journal of Physics, 19(12):125011, 2017.
https://doi.org/10.7554/eLife.74160.sa2Article and author information
Author details
Funding
National Science Foundation (DMR2004469)
 Michael Shelley
National Science Foundation (CMMI1762506)
 Michael Shelley
National Science Foundation (DMS1821305)
 Saad Ansari
 Matthew A Glaser
 Meredith D Betterton
National Science Foundation (ACI1532235)
 Saad Ansari
 Matthew A Glaser
 Meredith D Betterton
National Science Foundation (ACI1532236)
 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.
Acknowledgements
MJS acknowledges support from NSF grants DMR2004469 and CMMI1762506. SA, ARL, MAG, and MB acknowledge support from NSF grants DMS1821305, 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 ACI1532235 and ACI1532236), the University of Colorado Boulder, and Colorado State University. We thank Prof. Dimitrios Vavylonis for discussions on implementing flexible filaments.
Senior Editor
 Anna Akhmanova, Utrecht University, Netherlands
Reviewing Editor
 Pierre Sens, Institut Curie, CNRS UMR168, France
Publication history
 Preprint posted: September 16, 2021 (view preprint)
 Received: September 23, 2021
 Accepted: April 24, 2022
 Version of Record published: May 26, 2022 (version 1)
Copyright
© 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.
Metrics

 1,235
 Page views

 270
 Downloads

 2
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
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)
Further reading

 Computational and Systems Biology
Cycling of cosubstrates, whereby a metabolite is converted among alternate forms via different reactions, is ubiquitous in metabolism. Several cycled cosubstrates are well known as energy and electron carriers (e.g. ATP and NAD(P)H), but there are also other metabolites that act as cycled cosubstrates in different parts of central metabolism. Here, we develop a mathematical framework to analyse the effect of cosubstrate cycling on metabolic flux. In the cases of a single reaction and linear pathways, we find that cosubstrate cycling imposes an additional flux limit on a reaction, distinct to the limit imposed by the kinetics of the primary enzyme catalysing that reaction. Using analytical methods, we show that this additional limit is a function of the total pool size and turnover rate of the cycled cosubstrate. Expanding from this insight and using simulations, we show that regulation of these two parameters can allow regulation of flux dynamics in branched and coupled pathways. To support these theoretical insights, we analysed existing flux measurements and enzyme levels from the central carbon metabolism and identified several reactions that could be limited by the dynamics of cosubstrate cycling. We discuss how the limitations imposed by cosubstrate cycling provide experimentally testable hypotheses on specific metabolic phenotypes. We conclude that measuring and controlling cosubstrate dynamics is crucial for understanding and engineering metabolic fluxes in cells.

 Computational and Systems Biology
 Neuroscience
Seizure generation, propagation, and termination occur through spatiotemporal brain networks. In this paper, we demonstrate the significance of largescale brain interactions in highfrequency (80200 Hz) for identification of the epileptogenic zone (EZ) and seizure evolution. To incorporate the continuity of neural dynamics, here we have modeled brain connectivity constructed from stereoelectroencephalography (SEEG) data during seizures using multilayer networks. After introducing a new measure of brain connectivity for temporal networks, named multilayer eigenvector centrality (mlEVC), we applied a consensus hierarchical clustering on the developed model to identify the epileptogenic zone (EZ) as a cluster of nodes with distinctive brain connectivity in the ictal period. Our algorithm could successfully predict electrodes inside the resected volume as EZ for 88% of participants, who all were seizurefree for at least 12 months after surgery. Our findings illustrated significant and unique desynchronization between EZ and the rest of the brain in early to midseizure. We showed that aging and duration of epilepsy intensify this desynchronization, which can be the outcome of abnormal neuroplasticity. Additionally, we illustrated that seizures evolve with various network topologies, confirming the existence of different epileptogenic networks in each patient. Our findings suggest not only the importance of early intervention in epilepsy but the possible factor which correlates with disease severity. Moreover, by analyzing the propagation patterns of different seizures, we asserted the necessity of collecting sufficient data for identifying the epileptogenic networks.