The cytoskeleton generates force and reorganizes to perform important cellular processes [1], including cell motility [2,3], cytokinesis [4], and chromosome segregation in mitosis [5]. The cytoskeleton is made of polymer filaments, molecular motors, and associated proteins. The two best-studied cytoskeletal filaments are actin and microtubules [1]. It remains incompletely understood how diverse cytoskeletal structures dynamically assemble and generate force of pN to nN [1,2].
Force generation and reorganization in the cytoskeleton depend on the activity of crosslinking motor proteins that align and slide pairs of filaments (Fig. 1. Reorganization of actin networks by myosin motors is important for muscle contraction [6][7][8], cell crawling and shape change [9][10][11], and cytokinesis [4,12]. Microtubule sliding by crosslinking kinesin and dynein motors contributes to mitotic spindle assembly [5,[13][14][15][16], chromosome segregation [17][18][19][20], cytoplasmic stirring in Drosophila oocytes [21], and beating of cilia and flagella [22][23][24].
a e-mail: alamson@flatironinstitute.org (corresponding author)
Filament-motor interactions produce diverse cellular structures and dynamics, but linking molecular properties of motors to larger-scale assembly behavior remains challenging. Crosslinking motors vary in binding affinity, speed, processivity, and force-velocity relation. These same ingredients can be reconstituted and show dynamic self-organization into asters or contractile bundles [25][26][27], active liquid crystals [28][29][30][31], or other structures [32][33][34]. Even in reconstituted systems, our ability to predict and control dynamics and self-organization is limited.
Improved theory and simulation of cytoskeletal assemblies with crosslinking motors would allow better prediction of both cellular and reconstituted systems. Currently, few mesoscale modeling methods for filamentmotor systems are available between explicit particle simulations and continuum hydrodynamic theory. Explicit motor simulations have several existing software tools, including Cytosim [35], MEDYAN [36], and AFINES [37], and others [38]. Explicit motor simulations are straightforward to extend to include, for example, a new force-velocity relation or motor cooperativity. However, the cost of explicit particle simulations scales linearly or quadratically with the number of particles (depending on the type of interactions), making simulation of large systems challenging. Continuum models of coarse-grained fields can be computationally tractable and predict macroscopic behavior [39][40][41][42][43][44][45][46]. Current continuum models invoke symmetry considerations to determine the structure of the model without reference to an underlying microscopic mechanisms [39,[47][48][49], or simplify a microscopic model by making assumptions about the physics of motor [43,[50][51][52][53][54][55][56][57]. Furthermore, previous continuum theories have coarse-grained the filament distribution, with simplifying assumptions about the motor distribution. This presents an opportunity to better understand how the distribution of motors evolves and affects filament motion. Further development of mesoscale modeling techniques focusing on crosslinking motors could help bridge the gap between detailed explicit particle models and continuum theories.
To develop mesoscale modeling tools, we focus on the fundamental unit of a crosslinked filament network: two filaments with crosslinking motors that translate and rotate the filaments. We study three different model representations in a coarse-graining hierarchy and compare computational cost and accuracy. For explicit motors, we extend previous work that uses Brownian dynamics and kinetic Monte Carlo simulation to handle filament motion and binding kinetics [43,56,[58][59][60][61][62]. At the first level of coarse-graining, we average over discrete bound motors to compute the continuum meanfield motor density (MFMD) between filaments, and evolve this density according to a first-order Fokker-Planck equation [58]. This requires computing the solution to a single partial differential equation (PDE) for each filament pair, rather than separately tracking each individual motor. The MFMD determines the force and torque on each filament needed to evolve its position and orientation. At the second level of coarse-graining, we expand the MFMD in moments to derive a system of ordinary differential equations (ODEs) for the time evolution of the moments. While the moment expansion does not close, an approximate treatment of filament motion can be modeled by low-order moments. To compare these three model implementations, we consider test cases of parallel, antiparallel, and perpendicular filaments. Under the same initial conditions, the three model implementations give similar results on average. Remarkably, the reduced moment expansion achieves a computational cost that is 10 3 -10 6 lower than the other models, suggesting a route to computationally tractable large-scale simulations.
We consider a pair of rigid, inextensible filaments that move and reorient under the force and torque applied by crosslinking motors. Filaments move in three dimensions, experience viscous drag, and are constrained to prevent overlaps. Motors bind to and unbind from the filaments consistent with detailed balance in binding. Crosslinking motors walk with a force-dependent velocity toward filament plus ends and unbind when they reach the ends. We investigate models at three levels: an explicit motor model where motors are represented with a discrete density, a continuum mean-field motor density (MFMD) model, and a moment expansion model.
We model filament motion using Brownian dynamics, balancing the force applied by motors against viscous drag and constraint forces. Because the force that induces Brownian motion is typically smaller than that due to motors, we neglect Brownian noise [65].
Filaments translate according to the force-balance equation
where r i is the center of filament i with mobility matrix M i acted on by forces F n,i . The mobility matrix for a perfectly rigid rod in a viscous medium is
where I is the identity matrix and γ ,i and γ ⊥,i are the parallel and perpendicular drag coefficients with respect to the filament orientation ûi . Cytoskeletal filaments with length L i and diameter D fil typically have a large aspect ratio L i /D fil 1, so we approximate the drag coefficients using slender body theory [66].
The torque-balance equation is ui
where T n,i are the torques acting on filament i and γ θ,i is the rotational drag coefficient about the center of filament i.
The force and torque exerted by crosslinking motors depend on where motors are attached, the motor tether extension, and the relative position and orientation of filaments. Given the crosslinking motor distribution along the filaments ψ i,j (s i , s j ), where s i is the bound motor head position on filament i, the total crosslinking force and torque exerted by filament i on filament j are
where f i,j (s i , s j ) is the force exerted on filament j by the crosslinking motor attached at s i and s j (Fig. 1E). For brevity, we use subscripts on variables such as f i,j to indicate that these are functions of the relative position and orientation of filaments i and j. Our three model implementations all models use Eqs. ( 4) and ( 5) to compute the force and torque that evolve filament position and orientation but differ in how the computation of ψ i,j . We constrain the motion of filaments to prevent overlap, which avoids numerical instabilities introduced by a hard potential between filaments. To implement the constraint, we construct a vector ûmin that is perpendicular to both infinite carrier lines defined by ûi and ûj and parallel to the vector of closest approach between these lines. The vector ûmin is used to define two normal planes that confine the filaments, leading to the modified force and torque
Note that for filaments lying in the same confining plane and |û i • ûj | < 1, ûmin = 0 and our constraints break down. However, if only the first condition is satisfied, i.e., (anti)parallel filaments, T i,j = 0 and F i,j is parallel to ûi and ûj . After computing the force and torque, we numerically integrate Eqs. ( 1) and (3) to update filament position and orientation.
In our model motors bind and unbind, crosslink between two filaments, exert force and torque when crosslinking, and walk with a force-dependent velocity. Typically motor proteins diffuse in solution until they are near a filament, then stochastically bind to that filament. Once one head binds, the other head can bind to a second filament, forming a crosslink, or the motor can unbind. Crosslinking motors can unbind to a state with one head bound, or can unbind completely from both filaments. We consider an infinite reservoir of unbound motor proteins. The diffusion of motors in solution is fast relative to the motion of filaments, so we assume the motor reservoir has uniform, constant concentration. We neglect steric interactions between motors. This approximation holds for filaments sparsely pop-ulated with motors and motors that do not cluster on filaments or in solution.
Motors crosslinking filaments have a potential energy U i,j (s i , s j ) (Fig. 1). The energy depends on the motor head separation vector h i,j (s i , s j ) = r j + s j ûj -(r i + s i ûi ) that gives the motor tether extension
where r i,j = r jr i and r 2 i,j = r i,j • r i,j (Fig. 1e). The bound motor heads walk with a speed v i,j that depends on the force component on the motor head parallel to the walking direction, ûi • f j,i [67]. This projected force is used to determine the motor speed via the force-velocity relation, as discussed below. This model is based on processive microtubule motors such as kinesin and dynein, but a similar model has been used for myosin minifilaments [36,37].
In the explicit motor model individual bound motors are modeled, allowing fluctuations in bound motor number and binding kinetics that recover the correct equilibrium distribution of crosslinking proteins in the limit of no motor walking (Fig. 2a) [43,56,[59][60][61][62].
A motor diffuses in solution until one of its heads bind to a filament; we model this by an infinite reservoir of unbound motors with a uniform and constant concentration c o . Filaments have a linear binding site density , and the binding site has an association constant K a (units of μM -1 ). First motor head binding has rate
where L tot = i L i is the total length of filaments and k o,S is the bare (force-independent) unbinding rate for singly bound heads. All binding locations have equal binding probability. Singly bound motors unbind at rate k off,S = k o,S .
A motor with one head bound crosslink to another filament, which may stretch or compress its tether. This makes crosslinking kinetics force dependent; our models satisfy detailed balance in binding, so we recover the thermal equilibrium Boltzmann distribution in the limit of passive crosslinkers. Motor motion shifts the crosslinking distribution away from equilibrium. Motor unbinding rate can depend on the force applied to bound heads [68][69][70][71][72][73]. Previous work shows how this force dependence can be included while maintaining detailed balance in binding [62,74,75]. For simplicity, here we include the force dependence in the binding rate only and discuss possible implications below. With one head bound to filament i at position s i , the free motor head binds to filament j at position s j with a probability proportional to a Boltzmann factor of binding energy
with β = (k B T ) -1 (Fig. 3a). Here, S → C denotes the motor's transition from a single head bound (S) to crosslinking (C). The total binding rate is computed by integration over all binding positions on filament j
where k o,C is the bare (force-independent) unbinding rate for a crosslinking motor, K E is the crosslinking association constant. The unbound motor head explores a volume V bind centered about the bound head, computed as the integral of the unbound head's position weighted by the Boltzmann factor
Beyond the cutoff radius R cut,C , the integrand becomes small, enabling the use of a lookup table (Appendix B).
The probability distribution of binding position depends on the Boltzmann factor. We recover the proper binding distribution through inverse transformation sampling of Eq. ( 11) (Appendix B.2). As discussed above, the unbinding rate of a single head of a motor crosslinking two filaments is assumed to be force-independent,
Force-dependent unbinding affects the density of motor proteins most when stretched [70]; larger motor stretch occurs when external force is applied against the force generated by motors. Therefore, sliding filaments slowed only by drag, like those in active nematics, will be less affected by force-dependent unbinding than stationary filaments or jammed filaments like microtubules in mitotic spindles. We can include force-dependent unbinding in the explicit motor and MFMD model but not in the moment expansion model (Sect. 5). We chose the time step small enough that individual motors undergo only one transition per time step (Appendix A). The motor force-velocity relation is
where f stall is the motor stall force (Fig. 3b).
The bound motor distribution is
123 where δ(s i ) is the Dirac delta function and N i,j is the total number of motors crosslinking filaments i and j.
Here, s n and s n are the attached positions of the heads of the nth crosslinking motor. Motors with one head bound to filament i have a distribution
where N i is the number of one head bound motors on filament i. Only motors crosslinking exert forces between filament pairs, but χ i and χ j are needed to calculate the evolution of ψ i,j .
Under typical experimental conditions, there can be tens to thousands of crosslinking motors between a filament pair. Motor force and torque fluctuations occur because of stochastic motor binding and unbinding. As the number of motors increases, the standard deviation relative to the mean decreases as 1/ √ N . For our explicit motor model, antiparallel filaments with an average of 14 motors bound show a standard deviation in bound motor number of 27% of the mean. This shows that the fluctuations are quite significant for order 10 motors. The 1/ √ N scaling predicts that for an average of 1000 motors, the standard deviation would be only 3.2% of the mean. The force and torque scale similarly. Therefore, for large motor number, we may use the average motor distribution to derive a mean-field motor density (MFMD) to accurately describe force and torque on filaments by motors. We can then evolve the MFMD instead of explicit motors (Fig. 2b). We previously showed that the average steady-state density of crosslinking motors between stationary parallel filaments agreed well with a solution to a multidimensional Fokker-Planck equation (FPE) [58]. Here, we expand this approach to model crosslinking motor density between filaments in three dimensions, allow filament motion, and study time-dependent behavior of coupled systems of motors and filaments.
For a one-step binding model, the MFMD evolves according to
with motor velocity v i,j , motor crosslinking rate k on , and unbinding rate k off . To satisfy detailed balance in binding, we use the rates k on = 2k o ce -βUi,j (si,sj ) and k off = 2k o , with the effective concentration c (units nm -2 ) [58]. The factors of two occur because there are two ways a motor can crosslink. To numerically solve the hyperbolic Eq. ( 17), we use a first-order accurate upwind difference method (Appendix C).
The mean-field motor density model differs from the explicit model in that motors with one head bound are not modeled explicitly. To properly compare the different binding models, we establish a mapping of parameters between these two models (Appendix D), which gives
Some model parameters are difficult to measure directly. For example, the association constant K E may differ from K a if proteins change their molecular conformation when bound. We discuss an approach to estimate such parameters in Appendix E.
If filaments move slowly compared to the timescale of motor rearrangement, then a quasi-steady state approximation can be used. In the quasi-steady limit, the force and torque on filaments are computed from the steady-state MFMD [61]. The quasi-steady approximation is computationally efficient compared to numerical integration of the time-dependent PDE. A steady-state solution also provides a convenient route to compare our model implementations. At steady state, Eq. ( 17) becomes
Here, we choose functional forms of U i,j and v i,j consistent with previous models [36,37,58,60,61]. Motors have a potential energy
2 determined by the tether spring constant k cl and tether length h cl (Fig. 3a), which implies a motor crosslinking filaments i and j exerts a force on
on filament j. The force-velocity relation of a motor head attached to filament i while the other head is bound to j follows Eq. ( 14). Here, we assume motors that reach filament ends walk off, i.e., no end pausing. A semi-analytic steady-state solution can be derived for antiparallel filaments when motor tethers have zero length (h cl = 0) because the FPE is symmetric under the transformation i → j. For zero-tether-length motors to mimic their non-zero-length counterparts, we modify the zero-length motor's spring constant so both types of motors stall at the same extension h i,j = h stall . This implies k cl h stall = k cl (h stallh cl ) = f stall with the solution
where h stall = f stall /k cl . Note this choice changes the binding dynamics, because the potential energy is now larger for larger motor extension (Fig. 3A).
To find the steady-state solution, note that r i,j • ûi , r j,i • ûj = 0 and ûi • ûj = -1 for antiparallel filaments with centers aligned. Therefore, h i,j =
and v j,i depend exclusively on the sum of s i and s j , we make the change of variables ξ = s i +s j in Eq. ( 19) to find
There are three regions of solution determined by the force-velocity relation Eq. ( 14): ξ ≤ 0, 0 ≤ ξ ≤ h stall , and h stall < ξ. For ξ ≤ 0, v i,j = v j,i = v o and Eq. ( 21) becomes
where l o = v o /k o is the motor run length. This is solved with an integrating factor, giving
Applying the boundary condition ψ i,j (-L 2 , -L 2 ) = ψ i,j (-L) = 0, we remove the last term in Eq. ( 23) and rewrite the Gaussian integral as
For 0 ≤ ξ ≤ h stall , the velocity v i,j = v j,i = 1 -ξ h stall . Equation ( 21) becomes
Solving with an integrating factor, we find
We match the solution for ψ i,j (0) to Eq. ( 23) to enforce continuity. The exponential term in Eq. ( 26) can be approximated by a series expansion or integrated numerically. Here, we use numerical integration. For ξ > h stall , the velocity and velocity derivatives are zero, so
Since the motor velocity is zero at ξ = h stall , motors do not walk from ξ < h stall to ξ > h stall . A nonzero MFMD exists for ξ > h stall only if motors bind at these lengths. This appears as an integrable discontinuity at ξ = h stall .
A series expansion or reduced representation of a continuous distribution can lower the computational cost of solving a system's time evolution [57,76,77]. Here, we use low-order moments of the MFMD to calculate motor number, mean and standard deviations of motor head distribution, and filament motion. The moments of ψ i,j are
where k, l are nonnegative integers. The moments are symmetric under exchange of both filaments and powers so that μ k,l i,j = μ l,k j,i . The zeroth moment μ 0,0 i,j = N i,j is the total number of motors bound to the two filaments, and the first moments μ 1,0 i,j , μ 0,1 i,j are proportional to the mean motor head position along each filament
Ni,j . The first two second moments determine the standard deviation of motor head density
The symmetric second moment term μ 1,1 i,j determines the covariance of motor head position
The positional means, standard deviations, and covariance are used to reconstruct an approximate MFMD for visualization using a bivariate normal distribution (Fig. 2c, Videos 1-6).
Using the approximation of zero-length tethers as in Sect. (4.1) above, f i,j is a linear function of s i and s j . In this case, filament motion can be computed from low-order moments using Eqs. ( 4) and ( 5):
and
Substituting Eqs. ( 31) and (32) into Eq. ( 1) and (3) shows that only moments up to second order are needed to compute filament motion from crosslinking motors. Thus, motor and filament evolution can be written as a system of ODEs that depend on the dynamical evolution of the moments. This dynamical evolution is computed by taking the time derivative of Eq. ( 28) and substituting in the FPE (17)
However, this coupled system of equations for the moment time evolution does not close. Because the piecewise motor force-velocity relation is not linear, moments depend on higher-order moments recursively. Also, filament ends introduce boundary terms that prevent closure. Despite this, in certain limits a truncated moment expansion shows good agreement with the explicit and MFMD models.
We first introduce a linear approximation to the force-velocity relation (Fig. 3b)
This approximation is valid for h stall 1/k cl β, in which case motors do not bind beyond their stall stretch. We also require that v o 2k o 1/k cl β, ensuring that motors pulled towards the plus ends with ûi •f i,j > 0 move quickly into a regime -f stall < ûi •f i,j < 0, where the linear and piecewise force-velocity functions agree.
We substitute the linearized force-velocity function from Eq. ( 34) into the MFMD Eq. ( 17) to obtain
where κ = v o k cl /f stall is the rate at which motors reach their stall force. Integrating Eq. ( 35) directly returns the zeroth-moment equation
where we have defined q k,l i,j = Li Lj s k i s l j e -βUi,j ds i ds j and
with q k,l i,j representing source terms. Here, B l j (s i ) is a moment of the MFMD integrated over s j that is a function of s i , but in practice B l j only appears in the equations evaluated at filament endpoints, and so captures behavior of the motor density at filament ends. Therefore, we refer to the B l j (s i ) as boundary terms. To show this, we define the notation [A(s
The general moment evolution obtained by integrating Eq. (33) with Eq. ( 35) is
The boundary terms in square brackets contain moments and B l j an order higher than ∂μ k,l i,j /∂t. In Appendix G, we write the analogous time evolution for the B l j , and show that it does not close. Therefore, the moment evolution equations do not close.
To close the system of equations, we set the boundary terms to zero. Physically, this means we neglect motor unbinding from filament plus ends. If motors pause at plus ends, this approximation will lead to significant error. However, if motor unbinding is relatively rapid (including at filament plus ends), this is a good approximation. To explore the impact of not including these boundary terms, below we quantify the discrepancy between this model and the explicit motor and MFMD models. Neglecting boundary terms truncates the system of equations at second order, because only terms up to second order are needed to calculate force and torque on filaments.
We evolve equations (1,3,38) using solver_ivp in the scipy.integrate library [78]. This code uses the LSODA integrator, an Adams/BDF integration method that automatically detects stiffness, from the Fortran ODEPACK library [79]. The source terms q k,l i,j are analytically integrated in one dimension and then numerically integrated using the quad method also from scipy.integrate (Appendix F).
To test the degree of agreement between explicit motor and mean-field models, we first selected parameters based on microtubules and kinesin-5 motor proteins because they are relatively well-studied cytoskeletal proteins [81,82,84,85] (Table 1). We studied three characteristic sets of initial filament pair position and orientation: antiparallel, parallel, and perpendicular (Fig. 3, Video 1-3), and compared both stationary and moving filaments. We choose an initial condition with no motors bound to filaments, in order to observe the effects of time evolution of the motor density. For stationary filaments, we found good agreement for all three models. For moving filaments, we found qualitative agreement but fluctuations in motor dynamics and different end boundary conditions contributed to quantitative differences in filament motion. We measured the computational cost for stationary antiparallel filaments and found that the moment expansion model can give a dramatic improvement in performance.
When filaments are held stationary, motor density reaches or fluctuates around a steady-state solution (Fig. 4). To compare with the mean-field models, we averaged 48 realizations of each explicit motor simulation; the results agreed within error with the mean-field models (Fig. 4b-d). This agreement between models demonstrates that the mean-field models capture the average behavior of our explicit model.
Beyond the steady state, we characterize the evolution of motor number, force, and torque (Fig. 4e-m). In all configurations, the crosslinking motor number in the explicit motor model lags that of the MFMD and moment expansion models (Fig. 4e-g). The crosslinking rate in the two-step binding algorithm depends on the density of motors with one head bound, resulting in a slower approach to steady state.
For antiparallel filaments, force generation increases with crosslinking motor number (Fig. 4h, Video 1) because motors walk in opposite directions, causing the motor tether to stretch and generate force. If free to move, these antiparallel filaments would slide. No average sliding would occur for parallel filaments, and the small number of crosslinking proteins for perpendicular filaments results in small relative force (Fig. 4i,j). The average explicit motor motor torque in the ẑdirection shows significant fluctuations about the mean (Fig. 4k-m). Because motor torque increases for motors for one head bound and crosslinking motors. h-j Motor force in the x-direction (solid lines) and ŷ-direction (dotted lines). Individual explicit motor runs are represented as blue for both directions. k-m Motor torque in the ẑdirection from filament i on j. Full explicit motor model range not shown to better see average. n-p Steady-state motor probability density as a function of motor extension for semi-analytic (black), explicit motor, and MFMD models. Motor minimum extension is set by the separation of filaments at closest point of approach, 25 nm farther from the filament centers, the torque fluctuations increase with filament length.
We compared the steady-state distribution of motor extension for both explicit motor and MFMD models (Fig. 4n-p). (Note that the moment expansion loses this information in coarse-graining.) The distribution of motors crosslinking antiparallel filaments has two peaks (Fig. 4n). The larger peak represents the most probable binding distance Δy, and the second peak corresponds to motors near their stall extension h = Δy 2 + h 2 stall . The shape of the distribution results from motor kinetics, walking, and stalling. Motors on parallel filaments show a peak at Δy (Fig. 4o, Video 2), but no second peak because the motor heads walk in the same direction with similar speed. For motors crosslinking perpendicular filaments, the extension distribution is singly peaked and broader than for parallel filaments (Fig. 4p, Video 3). This occurs because the parallel force component on perpendicular filaments increases more gradually as the motors extend, causing a more gradual decrease in motor speed. This broad distribution indicates a larger average force per motor for perpendicular filaments compared to aligned filaments.
Here, we consider the same three filament starting configurations and allow filament motion (Fig. 5, Videos 4-6). The final filament position and orientation are comparable for the explicit motor and MFMD models, while the moment expansion model overestimates the range of filament translation and rotation (Fig. 5b,c; note that filament rotation only occurs for the perpendicular initial configuration).
To compare motor activity between models over the whole simulation, we calculated the total work done by motors. We numerically integrate both filaments using the trapezoid rule [87],
where θ j is the angle the vector ûj rotates through over the simulation. The infinitesimal vector dθ i = θi dθ i where θi = ûi × ui
Total work computed for the mean-field models is within error of the explicit motor model (Fig. 5d). We note that the explicit motor model produces greater total work because fluctuations in motor binding cause fluctuations in sliding direction which generate larger work. Motors generate rotational work only for initially perpendicular filaments, due to the constraints. The magnitude of the rotational work is relatively small because filaments rotate slowly (due to high rotational drag and low motor torque), and this slower velocity produces less work in the overdamped limit.
The crosslinking motor number depends on the filament overlap length, which changes as filaments move (Fig. 5e-j). The crosslinking motor number in the explicit motor model lags the mean-field models initially due to differences in binding, but becomes comparable after the initial transient. As antiparallel filaments slide apart, their overlap decreases so fewer motors crosslink, while crosslinking motors continue to unbind at a constant rate. However, the overlap length has little effect on the number of motors with one head bound (Fig. 5e,h). The dynamics of motor number for parallel stationary and moving filaments are nearly identical because there is negligible sliding. (Fig. 5f,i). Moving perpendicular filaments maintain a similar overlap length to stationary perpendicular filaments, leading to an approximately constant motor number, until the plus-ends move close together (Fig. 5g,j). Then, motors continue to bind but immediately walk off, producing little force or torque.
The motor force between antiparallel filaments rapidly reaches a force plateau which persists until the antiparallel overlap length is small enough that motor binding is negligible (Fig. 5k). The nearly constant force implies that motor extension decreases as the number of crosslinking motors increases to give a constant sliding speed (Fig. 5n, Video 4). This steady-state force is an order of magnitude smaller than the stall force (Table 1). The moment expansion model shows a slower decrease in force as the overlap approaches zero compared to the MFMD model (Fig. 5h). This is a consequence of our neglect of boundary terms, which physically means neglecting motor dissociation at filament ends. This unphysical slow force decrease drives filaments beyond the zero overlap configuration to larger than expected separation (Fig. 5b).
Parallel filaments remain with their centers aligned on average because sampling the full distribution of motor crosslinking extension generates restoring force for any fluctuations away from full overlap (Eqs. 11, 13). Neither the MFMD nor the moment expansion models produce a net force, but in the explicit motor model fluctuations in motor number and binding lead to force and position fluctuations (Fig. 5f,i,l,Video 5). For perpendicular filaments, the small number of crosslinking motors results in large force fluctuations in the explicit motor model (Fig. 5j). The mean-field models show a rapid increase to half the maximum force of the antiparallel configuration followed by a decrease as the filaments align parallel (Fig. 5k,m). The lag caused by the two-step binding model is more apparent here because the explicit lower motor number means filaments move more slowly into the parallel configuration where binding is favored (Video 6).
To compare the accuracy and computational cost of our models, we focus on stationary antiparallel fil- aments because we can compare to the semi-analytic solution. Antiparallel filaments are also the main configuration in which motors generate extensile force, important for mitotic spindle assembly and dynamics in active nematics. We vary the time step Δt and MFMD grid spacing Δs and compare the error with the semianalytic solution. The central-processing unit (CPU) time measures the computational cost as a function of simulation parameters.
The solution error is the average magnitude of the deviation of the steady-state numerical solution from ψ i,j of Eqs. ( 23), (26), and (27),
where ψ i,j is either the average explicit motor distribution (over 48 simulations) or the MFMD distribution.
The size of the time step Δt does not change the error of explicit motor or MFMD simulations (Fig. 6a), because the steady-state solution is time independent. The number of calculations increases linearly with the number of time steps N t /Δt, making the CPU time approximately inversely proportional to Δt. The MFMD error scales near-linearly with grid spacing Δs as expected for a first-order upwind difference method (Fig. 6b). The CPU time scales approximately as Δs -2 , proportional to the number of grid points N grid ∝ Δs -2 .
Explicit motor simulations have a cost that is linear in the motor number, but the cost is constant for the MFMD and moment expansion models (Fig. 6c). Fewer explicit motor simulations (24 realizations) were needed to achieve sufficient statistics. We also note that at higher concentration, the mean-field models return results closer to those of the explicit model because stochastic fluctuations average out. The explicit motor model has a cost linear in filament length (due to the larger number of bound motors on longer filaments), while for the MFMD model it is quadratic (Fig. 6d). The cost of the moment expansion model is length independent.
To improve modeling methods for cytoskeletal filaments crosslinked by motors (Fig. 1), we studied crosslinked filament pairs and compared an explicit motor model to two levels of coarse-grained mean-field motor models (Fig. 2). The explicit motor model uses Brownian dynamics and kinetic Monte Carlo to describe individual motor binding and unbinding, motion, and force generation. In the first level of coarse graining, we average over individual motors and solve a PDE for the mean-field motor density (MFMD). To further coarse grain, we compute a moment expansion of the MFMD and solve a system of ODEs for the motor moments and filament motion.
We compared the model implementations for filaments that are initially antiparallel, parallel, or perpendicular (Fig. 3). When filaments are held stationary, the motor distribution reaches a steady state with similar average motor distribution, force, and torque for the three implementations (Fig. 4). The explicit motor simulations showed significant fluctuations that by construction are not present in the mean-field models. Interestingly, we found that a significant portion of crosslinking motors on antiparallel filaments do not reach their stall force for our parameter set.
When filaments move, the final filament separation is similar for the explicit motor and MFMD models, although the moment expansion model overestimates the range of displacement and reorientation as a result of neglecting boundary terms (Fig. 5). The dynamics of bound motor number, force, and torque were similar for the MFMD and moment expansion models. Motor fluctuations in the explicit motor model lead to greater overall work done by motors.
To compare computational cost across the model implementations, we studied stationary filaments and motors at steady state (Fig. 6). Both mean-field models have a simulation time independent of motor concentration, potentially making them faster than explicit models for systems with many motors. The moment expansion model's CPU time is also independent of filament length, which could make it particularly efficient for systems with long filaments. Overall, the moment expansion model was 10 3 -10 6 faster than the other models. This method could therefore be useful for simulating bulk active filament networks.
Future work could address the simplifying assumptions and approximations made in the moment expansion model. An improved treatment of boundary terms may improve the computation of filament motion.
Incorporating additional motor physics into the moment expansion model, such as non-zero-length motors, forcedependent detachment, and steric interactions between motors could improve its ability to simulate microscopic motor behavior at the mesoscale, bridging current explicit motor and continuous active network theories.Implementing the moment expansion model in systems of many filaments is of interest for testing whether the improvements in computational cost we identify are present in larger systems.
The online version contains supplementary material available at https://doi.org/ 10.1140/epje/s10189-021-00042-9.
Our kinetic Monte Carlo algorithm assumes that multiple binding/unbinding events do not occur in the same time step Δt. As Δt becomes large relative to the kinetic rates, this approximation fails. A time step is appropriate if the maximum probability of two events occurring in Δt satisfies max{P (C(Δt) ∪ B(t )|A(0))} < δ (42) for a tolerance δ, where A, B, and C denote motor bound states (including unbound, single head bound, and crosslinking (
While no analytic solution exists, t max can be numerically computed.
There are four unique processes that must be considered with a two-step binding process with unbound (U), single head bound (S), and crosslinking (C) states:
Equation (11) gives the transition probability of a singly bound motor crosslinking as an integral of a Boltzmann factor. If h cl = 0, kon,C is functionally similar to an error function. However, to model non-zero-length tethers, we numerically integrate Eq. (11). Rather than directly numerically integrating at each time step, we precompute a lookup table.
The cumulative distribution function (CDF) of Eq. ( 11), is a function hi,j. All other variables in the integral are constant for a given motor species. We reduce the CDF dimensionality by considering the lab position of each bound motor head and an infinite carrier line defined by the position and orientation of the unbound filament. Binding is then determined by the minimum distance r ⊥ between the bound motor head position and the filament ends [s-, s+] on the carrier line.
The carrier line CDF is
e -βU (r ⊥ ,s ) ds , (44) allowing us to write the crosslinking rate as
We notice that e -βU i,j is symmetric in s, so CDF(r ⊥ , s) -CDF(r ⊥ , 0) is anti-symmetric. Therefore, instead of integrating from negative infinity, we use CDF (r ⊥ , s) = sgn(s) s 0 e -βU (r ⊥ ,s ) ds (46) and ( 45) to find the crosslinking rate.
We find the values of Eq. ( 46) by Gauss-Konrad integration. The accuracy desired sets the maximum values for s and r ⊥ . The integrand is always positive for real values of s and r ⊥ , so the CDF asymptotes for large values of either variable. The maximum of the integral is the point when the Boltzmann factor drops to the accuracy limit δ. Therefore, the lookup table domain is
Given a specified grid spacing Δs, Δr, the memory required for the lookup table scales as (smax/Δs) × (r ⊥,max /Δr).
Since the lookup table is not a continuous function, we interpolate values between discrete grid points. The 2D linear interpolation for input values of r ⊥ and s is
where CDFm,n = CDF(mΔr, nΔs) are the lookup table values at m and n if r ⊥ lies within mΔr and (m + 1)Δr and s lies within nΔs and (n + 1)Δs.
When a motor head binds to a filament, the binding position probability distribution function (PDF) is defined by the Boltzmann factor. We sample the PDF by using the lookup table. To transform a uniform random variable X to random variable Y with an arbitrary PDFY ,
Since the lookup table holds the CDF values and given a random number from a uniform distribution, we apply a combination of search and interpolation to quickly find the corresponding random number from the PDF. The algorithm is as follows 1. Sample a uniform random number X ∈ [0, CDFmax]. Note that the maximum value does not need to be 1. 2. Given r ⊥ , locate index m such that mΔr ≤ r ⊥ ≤ (m + 1)Δr 3. Use m to find the set of indices {n-, n+} such that CDFm,n -≤ X ≤ CDFm,n -+1 and CDFm+1,n + ≤ X ≤ CDFm+1,n + +1. 4. Use the CDF values to interpolate the binding locations s-, s+ corresponding to the perpendicular distances r-= mΔr and r+ = (m + 1)Δr. For example,
Note that sis not necessarily less than s+. 5. Find s by interpolating the across the lookup table grid with respect to r
While this algorithm succeeds in most circumstance, the low slope of the CDF at large values of s can cause errors. For example, if the lookup table has the form of Fig. 7 anda protein is located at a perpendicular distance of r ⊥ = 35 nm, given a random number of X = 10 3 , no value for swill be found since CDF(30, smax) < 10 3 . To correct for this, we solve for s using a binary search algorithm.
The binary search algorithm is as follows 1. Determine if CDFm,n max or CDFm+1,n max is less than X. If CDFm,n max < X, set s-= smax. If CDFm+1,n max < X, set s+ = smax. 2. Find other s± using the inverted lookup table and Eq. ( 50) or (51). 3. Find the average of sand s+. 4. Use the lookup table interpolation algorithm to find the CDF(r ⊥ , savg). 5. If CDF(r ⊥ , savg) > X set the larger of the two s± values to savg. Otherwise, set the smaller of the two to savg. 6. Repeat steps 3-5 until | CDF(r ⊥ , savg)-X| < δ for some desired tolerance δ.
This process converges at a rate O(log 2 (δsmax)).
We approximate the solution ψi,j(si, sj, t) by discretizing the solution in time and space ψi,j(si, sj, t) → ψ m,n,k i,j
= ψi,j(mΔs, nΔs, kΔt) (53) for ψ m,n,k i,j ∈ R (M i +1)×(M j +1)×k , where Mi is the number of discretized points along filament i. Additional boundary points for m, n = 0 are added.
We use forward Euler time-stepping so our discrete differential operator for time is
).
(
To solve the hyperbolic FPE (17), we use a first-order accurate upwind method [88]. The differential operator for si becomes ∂ψi
).
(
Note this only holds for the indices 0 < m and 0 < n. The matrix representation for Eq. ( 55) is
where cm and dm are chosen to satisfy the boundary conditions. We choose the notation m,n for this matrix. To differentiate along sj, we use the identity ψ m,n,k i,j = ψ n,m,k j,i , apply m,n on the matrix, and then convert back,
which in index notation is ψ m,a i,j ( T ) n,a . For brevity, we use the notation ψ m,a i,j ( T ) n,a = ψ m,a i,j a,n . The discretized Fokker-Planck Eq. ( 17) is then
where U m,n,k i,j and v m,n,k i,j are the discretized potential and velocity at time kΔt. Note that U m,n,k i,j
In cases where the flux of the motors ∂(v i,j ψ i,j ) ∂s i is known at the boundaries, we construct to satisfy the requirements. When filaments are in solution, there is zero flux from the minus ends, so all cm = 0. In our simulations, motors walk of filament ends with out pausing, so dM i -1 = -1 and dM i = 1 with all other dm = 0. Although not modeled in this paper, some biological motors end pause at filament plus ends. To model this, dM i -1 = -1 and every other dm = 0.
To relate binding parameters of the one-step and multi-step binding models, we use that at steady state, the motor distribution ψi,j should be equivalent for both models. Since we only compare binding kinetics, we simplify the Fokker-Planck equation to keep only the binding terms: in Eq. ( 17), we set vi,j = vj,i = 0, ∂ψi,j ∂t = 2koce -βU i,j -2koψi,j.
(
123
Eur. Phys. J. E (2021) 44 :45
The steady-state solution is ψi,j = ce -βU i,j , (60) which is a Boltzmann factor multiplied by an effective concentration.
The multi-step binding model can be written ∂ψi,j(si, sj) ∂t = K E ko,C (χi + χj)e -βU i,j -2ko,C ψi,j, (61) ∂χi(si) ∂t = coKa ko,Sko,Sχi
where χi is the mean-field density of motors with one head bound to filament i (cf. Eq. 16). We define K E = KE/V bind and solve for the steady state, giving (
The equations for χi and χj have the forms (72) After plugging Eq. ( 71) into (72) (
This can be rearranged into the form
which implies that Y (t) and X(s) each satisfy a Fredholm equation of the second kind. Both A and B are continuous given K(s, t) = e -βU i,j (s,t) , so the Fredholm equations of the second kind have unique solutions. By inspection, the solution to Eqs. ( 67) and ( 68) is X(s) = Y (t) = D. When we substitute this solution in Eqs. ( 65) and ( 66), we find χi = χj = Kaco and ψi,j = 2 KaK E coe -βU i,j .
(
The experimental parameters for motor binding are not always independently measured. If all but one binding parameters are known, then the unknown parameter can be found from Eq. ( 18) and the ratio of the number of motors with one head bound and number of motors crosslinking.
As an example, suppose we wish to find KE. The number of motors with one head bound is NS = coKa L, where L is the filament length. In vitro experiments [86] can measure the crosslinking motors number N d . Integrating Eq. ( 75), we obtain the model prediction for the number of crosslinking motors as NC = co 2 KaK E L i L j e -βU i,j dsidsj. (77) For fully parallel or antiparallel filaments of the same length with adjacent centers, the total number of motors in Eq. ( 77) is proportional to L. If L 2/βk cl , the Gaussian integral ≈ L 2π/βk cl e -βk cl r 2 ⊥ , where r ⊥ is the center-tocenter separation between filaments. The ratio of the number of crosslinking motors relative to the number motors with one head bound is
allowing us to estimate K E = ρ βk cl 2π e βk cl r 2 ⊥ .
The source terms in the moment expansion require a double integral over two filaments. To lower the numerical integration's computational cost, we find an analytic solution for either the semi-integrated term Q l j (si) or the fully integrated term q k,l i,j . The integrated source terms are (
This integral has an analytic form in terms of error functions, which can be rapidly computed. For l = 0, 1, 2, 3, we find
To generally define boundary conditions, instead of integrating over both si and si, we integrate over just one variable. This makes the boundary condition a function of a single filament attachment position. For example, the boundary terms for the first filament are (
These boundary terms are evaluated at -Li/2 and Li/2. We derive a recursion relation by integrating Eq. ( 85) over (
This shows that the boundary terms do not close. However, if the higher-order terms or their coefficients are small compared to the moments μ k,l i,j , we may take a zeroth-order approximation. We consider this approximation in Sect. 6.