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
Article 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.
Version 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,937
 views

 383
 downloads

 12
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
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
 Medicine
Background:
Preterm birth is the leading cause of neonatal morbidity and mortality worldwide. Most cases of preterm birth occur spontaneously and result from preterm labor with intact (spontaneous preterm labor [sPTL]) or ruptured (preterm prelabor rupture of membranes [PPROM]) membranes. The prediction of spontaneous preterm birth (sPTB) remains underpowered due to its syndromic nature and the dearth of independent analyses of the vaginal host immune response. Thus, we conducted the largest longitudinal investigation targeting vaginal immune mediators, referred to herein as the immunoproteome, in a population at high risk for sPTB.
Methods:
Vaginal swabs were collected across gestation from pregnant women who ultimately underwent term birth, sPTL, or PPROM. Cytokines, chemokines, growth factors, and antimicrobial peptides in the samples were quantified via specific and sensitive immunoassays. Predictive models were constructed from immune mediator concentrations.
Results:
Throughout uncomplicated gestation, the vaginal immunoproteome harbors a cytokine network with a homeostatic profile. Yet, the vaginal immunoproteome is skewed toward a proinflammatory state in pregnant women who ultimately experience sPTL and PPROM. Such an inflammatory profile includes increased monocyte chemoattractants, cytokines indicative of macrophage and Tcell activation, and reduced antimicrobial proteins/peptides. The vaginal immunoproteome has improved predictive value over maternal characteristics alone for identifying women at risk for early (<34 weeks) sPTB.
Conclusions:
The vaginal immunoproteome undergoes homeostatic changes throughout gestation and deviations from this shift are associated with sPTB. Furthermore, the vaginal immunoproteome can be leveraged as a potential biomarker for early sPTB, a subset of sPTB associated with extremely adverse neonatal outcomes.
Funding:
This research was conducted by the Perinatology Research Branch, Division of Obstetrics and MaternalFetal Medicine, Division of Intramural Research, Eunice Kennedy Shriver National Institute of Child Health and Human Development, National Institutes of Health, U.S. Department of Health and Human Services (NICHD/NIH/DHHS) under contract HHSN275201300006C. ALT, KRT, and NGL were supported by the Wayne State University Perinatal Initiative in Maternal, Perinatal and Child Health.

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