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

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

Figures

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

(A) aLENS simulates dynamics of rigid filaments crosslinked and driven by motors, thermal fluctuations, and steric interactions. Motors bind to, unbind from, and walk along filaments. (B) To achieve high efficiency, aLENS computes motor forces implicitly, and steric interactions through a novel geometric constraint method that avoids filament overlaps. (C1-C3) Example simulation of microtubules organized into asters by minus-end-directed motors. The 300 s Brownian simulation contains 3200 microtubules, each 1 µm long, inside a sphere of radius 3 µm. The initial position of each microtubule is random and the half of each filament on the minus-end is colored pink. Three end-pausing dynein motors are fixed at the minus-end of each microtubule and walk toward the minus-end of any microtubule they crosslink. After initial contraction into a single large aster, strong steric interactions in the aster center break up the system into several smaller asters and a bottle-brush structure. (C4) Motors are highly concentrated at the centers of asters.

Figure 1—video 1
Contraction and break-up of simulated microtubule asters.

The simulation details are described in Figure 1.

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

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

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

The left panel shows the simulation geometry and the sampling procedure. Microtubules are randomly initialized with orientations along the +x (pink) or -x (white) directions. XCTK2 motors are colored green. Nm is the number of XCTK2 motors per microtubule. We sample the local average polarity and straining velocity by inserting planes orthogonal to the x-axis into the collected data, matching the photobleaching technique used in experimental measurement (Fürthauer et al., 2019). For every sampling plane (e.g. the blue pane in the snapshot), we choose five sample points symmetrically on this plane and draw a square sampling window with edge length 0.2µm around each sample point. For each sampling window, we compute the average polarity Px along the x-axis for all microtubules intersecting this sampling window at a given time. We then compute the velocities, averaged over 10s (a duration chosen to match the experimental timescale), of microtubules intersecting each sampling window and moving along the +x and -x directions. V+x and V-x are computed from those two groups for each sampling window. The straining velocity is computed as Vx=V+x-V-x. Therefore, for every sampling window at each sampling timestep we have a pair of data values Px,Vx. The right three panels show the joint probability distribution of (Px,Vx) computed from 900,000 sampling planes for each simulation, for Nm=10,20,30, respectively.

Strong scaling (fixed system size while increasing number of cores) efficiency of a system similar to but more than 10 times larger than that shown in Figure 6, comprising 1 million microtubules and 3 million motors.

There are in total approximately 8 million constraints per time step, which is changing from step to step because collision pairs are changing and crosslinkers are stochastically binding and unbinding. The simulation is run for 100 computing steps with 1 data-saving step and the average per-step wall-clock time is shown in the figure.

Figure 5 with 2 supplements
Results for the bundling-buckling simulation of 100,000 microtubules and 500,000 dynein motors in the periodic simulation box of 600×10×10µm.

Brownian motion of microtubules is turned off. Each dynein has one non-motile head permanently attached to a microtubule and the other motile head walks processively with maximum velocity 1µms-1. If bound, the motile head moves toward the microtubule minus-end, and detaches upon reaching it. Detailed parameters for this motor are tabulated in the Appendix 1. Every microtubule has 5 dynein motors permanently attached to randomly chosen, fixed locations along the length. The initial configuration of microtubules is randomly generated, with their orientations sampled from an isotropic distribution and centers uniformly distributed within a cylinder of length 600µm and diameter 0.3µm. The motile heads of all dynein motors are unbound initially. (A, B, C) The bundle at t=0s, 4s, and 7s. Microtubules are colored by their local nematic order parameter Slocal=32QijQij, with Qij=pipj-13δij, p being the unit orientation vector of each microtubule pointing from the minus to the plus end, and δ the Kronecker delta tensor. The average . is taken over each microtubule plus all microtubules that are directly crosslinked to it by dynein motors. (A1, B1, and C1) Zoom-in views of the small region marked by red box in A, B, and C. (C2) The same region in C1 but colored by Nd, the number of microtubules averaged over when computing Slocal. (D1 and D2) The joint probability distributions Slocal and Nd for each microtubule for the entire systems at t=0.1s, when the dyneins crosslink microtubules but microtubules barely move from initial configuration, and at t=7s, when the bundle is nematic. (E) The average trajectories (solid lines) and their standard deviation (shaded area) of left-moving and right-moving microtubules. Dashed lines show linear fits to the average trajectory after t=4s, with results VRVL250nms-1. (F) The normal stresses and the weighted average Slocal¯ over time. Due to the symmetry in the y,z directions, only their average is shown σyy,zz=12(σyy+σzz). Collision stress is positive (extensile) and crosslinker stress is negative (contractile). The weighted average Slocal¯=NdiSlocali/Ndi.

Figure 5—video 1
Contraction and buckling of a long microtubule-motor bundle.

The bottom panel is a zoom-in view to the area in a gray box in the top panel. The simulation details are described in Section Bundle formation and buckling in a filament band.

Figure 5—video 2
Motor motion and stretching during the contraction and buckling of a long microtubule-motor bundle.

This is a zoom-in view to the area in a gray box in the bottom panel in Figure 5—video 1. The simulation details are described in Section Bundle formation and buckling in a filament band.

Results for the polarity sorting simulation in a spherical shell.

Initially, 100,000 0.25 µm-long filaments modeling microtubules and 200,000 motors modeling crosslinking kinesin-like proteins are placed between two concentric spherical shells with radii rin=5µm and rout=5.102µm, to maintain the volume fraction of filaments between these two shells at 40%. Initially, all filaments are evenly distributed on the spherical shell, with their orientation randomly chosen to be either ±eθ at each point, where eθ is the polar basis norm vector of spherical coordinate system. The pure filament system is relaxed for 1s to resolve the overlaps in the initial configuration. Afterwards at t=0, 200,000 motors are added to the system homogeneously distributed between the two shells. Sample points are evenly placed to measure the statistics by averaging the volume within 0.25µm from each sample point. (A) The configuration at t=0. Filaments are colored by their polarity, while motors are colored as black dots. Only randomly selected 10% of all motors (same after) are shown in the image to illustrate the distribution. (B) Randomly selected trajectories of filaments from t=0s to 1s. Trajectories are colored by time. It is clear that filaments move along the meridians. (C1-C5) Configuration and statistics at t=4s. (C1) The filaments and motors. Motors clearly concentrate in some areas. (C2) The motor number density, that is, number of motors per 1µm3. (C3) The nematic director field (shown as black bars) and the nematic order parameter S. (C4) The filament volume fraction. (C5) The divergence of polarity field p non-dimensionalized by filament length, that is, change of mean polarity per filament length. (D) The development of the correlation between motor number density n/nave and the polarity divergence field, at different times of the simulation. Clearly high n/nave are correlated with positive polarity divergence. (E1-E5) Configuration and statistics at t=0.27s for a comparative simulation where motors have end-pausing, arranges in the same style as C1-C5. This case shows significant contraction instead of polarity sorting as filaments are pulled away from the north and south poles and the overall volume fraction significantly increases to approximately 60%. The structure becomes densely packed and does not significantly evolve further. (F) The correlation between motor number density n/nave and the local filament volume fraction. For the polarity sorting case at t=4s the motor number density correlates with low filament volume fraction. This is not seen in the end pausing (EP) case. (G) A schematic for the correlations shown in D and F.

Figure 7 with 2 supplements
Results for the aster formation simulations with Brownian motion of simulated microtubules turned on (BMT) and off (NBMT).

Initially, 40,000 0.5µm-long filaments modeling microtubules and 80,000 motors modeling crosslinking kinesin-like proteins are placed in a periodic cubic box of 10×10×10µm with uniform distribution. Filament orientations are isotropic and motors are all in the unbound state. Motors are assumed to have two minus-end-directed walking heads with symmetric properties. They are assumed to pause when they reach the minus end of filaments until detaching. Detailed parameters are tabulated in Appendix B. (A, D, and E) Simulation snapshots. Each filament is shown as a cylinder colored in half pink (minus end) and half white (plus end). (A) The initial configuration for both NMT and BNMT cases. Each motor is colored as a green dot. (D and E) The snapshot for both cases at t=35s. D1-2 and E1-2 Expanded views of a aster core for D and E. Only doubly bound motors are shown in D and E (in green color), and in D1-2 and E1-2 (colored by the spring force). Negative values mean the crosslink forces are contractile (attractive). (B) The radial distribution function (RDF) g(r) for the minus ends of all filaments at t=0.1s (dashed lines) and t=35s (solid lines). The first peak of g(r) at r=25nm corresponds to close contacts between filaments. The second peak of g(r) at r=78nm=25nm+53nm corresponds to the minus ends of filaments crosslinked by motors whose rest length is 53nm. Blue and red lines are results for the BMT and NBMT cases, respectively. (C) The collision (solid) and crosslinker (dashed) pressure for BMT (blue) and NBMT (red) cases. Pressure is defined as the trace of the stress tensor: Π=13Trσ. The collision pressure ΠCol is positive (extensile), and the motor pressure ΠLinker is negative (contractile). The inset plot shows the pressure for the NBMT case in the initial stage of the simulation. The black dashed lines mark the time t=4s.

Figure 7—video 1
Aster formation in bulk of Brownian microtubules.

This is a zoom-in view to the BMT case shown in Figure 7.

Figure 7—video 2
Aster formation in bulk of Non-Brownian microtubules.

This is a zoom-in view to the NBMT case shown in Figure 7.

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

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

Figure 8—video 1
Filament-motor assembly for the Dcyl/LMT=1 case.

The bottom panel is a zoom-in view to the area in a gray box in the top panel. The simulation details are described in Section Confined filament-motor protein assemblies.

Figure 8—video 2
Filament-motor assembly for the Dcyl/LMT=3 case.

The bottom panel is a zoom-in view to the area in a gray box in the top panel. The simulation details are described in Section Confined filament-motor protein assemblies.

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

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

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

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

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

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

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

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

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

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

Appendix 5—figure 1
The radial distribution function g(r) and structure factor S(q) of identified aster centers at steady state, for BMT and NBMT cases.

‘DBSCAN’ and ‘Graph’ refer to two different methods of identifying aster centers, based on spatial locations of all microtubule minus ends, and the crosslinking connectivity, respectively. 500 snapshots at simulation steady state are used to compute g(r) and S(q), for each case.

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

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

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

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

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

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

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

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

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

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

Author response image 1
Comparison of the equation of state for 1µm microtubules at room temperature across the isotropicnematic transition range of volume fractions.

ϕcp = 0.9036 is the close packing volume fraction. The reference data is adapted from [7], generally accepted as the ground truth for rigid Brownian liquid crystal equation of state. Simulations using linearized WCA potential uses the code by [8], which has a constant repulsive force when the distance between rods is smaller than 0.812DMT, where DMT is the microtubule diameter.

Tables

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

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

ProcessRateValue
U(SA,SB)Ron,S(x)ko,S3ϵKa4πrc,S3iLin,i(x)
(SA,SB)URoff,Sko,S
(SA,SB)DRon,D(si)ko,DϵKejLjdsjexp[(1λ)βE(f(s))]
D(SA,SB)Roff,D(si,sj)ko,Dexp[λE(f)]
Appendix 1—table 1
The parameters of the two springs controlling extension and bending, respectively.
ParameterExplanationUnit
End-pausingTrue or FalseND
One head fixedTrue or FalseND
λenergy factorND
λP,APparallel to anti-parallel factorND
0free lengthµm
rccapture radiusµm
κHookean spring constantpNµm-1
Fstallstall forcepN
dUunbound diffusivityµm2s-1
ϵbinding site densityµm-1
vmmax walking velocityµms-1
Kaassociation constant (US)(μmol/L)-1
ko,Soff-rate constant (US)s-1
Keeffective association constant (SD)ND
ko,Dforce-independent off-rate constant (SD)s-1
dSsingly bound head diffusivityµm2s-1
dDdoubly bound head diffusivityµm2s-1
vSsingly bound walking velocityµms-1
xcforce-dependent unbinding lengthµm
Appendix 1—table 2
Properties of crosslinkers used in the main text.

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

ParameterKinesin-5DyneinKinesin-1Inactivated Kinesin-1
End-pausingTrueFalseFalseFalse
One head fixedFalseTrueTrueTrue
λ0.2580.50.50.5
λP,AP1111
00.0530.0400.050.05
rc0.0390.0330.0380.038
κ300.0100.0100.0100.0
Fstall5.01.07.07.0
dU1.01.01.01.0
dS00[0,10-2][0,0]
dD00[0,10-2][0,0]
ϵ1,625400400400
vm[-0.1,-0.1][0,-1.0][0,1.0][0,0]
Ka[90.9,90.9][100.0,100.0][0,10.0][0,10.0]
ko,S[0.11,0.11][0.1,0.1][0,1.0][0,0.1]
Ke[90.9,90.9][100.0,100.0][0,10.0][0,10.0]
ko,D[0.11,0.11][0.1,0.1][0,1.0][0,0.1]
Appendix 2—table 1
The transition rates between all possible states of a crosslinker U(SA, SB)D.

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

ProcessRateValue
U(SA,SB)Ron,S(x,t)3ϵKako,S4πrc,S3iLin,i(x,t)
(SA,SB)URoff,Sko,S
(SA,SB)DRon,D(si,t)ϵKeko,DVbindjLjdsjexp[βκxl(1λ2(oDfil)2xc(oDfil))]
D(SA,SB)Roff,D(si,sj,t)ko,Dexp[βκxl(λ2(oDfil)2+xc(oDfil))]
Appendix 7—table 1
The parameters of the two springs controlling extension and bending, respectively.

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

Rolespring stiffness constantfree length
BendingκBB0
ExtensionκEE0

Download links

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

Downloads (link to download the article as PDF)

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

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

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