Linking path and filament persistence lengths of microtubules gliding over kinesin

Microtubules and kinesin motor proteins are involved in intracellular transports in living cells. Such intracellular material transport systems can be reconstructed for utilisation in synthetic environments, and they are called molecular shuttles driven by kinesin motors. The performance of the molecular shuttles depends on the nature of their trajectories, which can be characterized by the path persistence length of microtubules. It has been theoretically predicted that the path persistence length should be equal to the filament persistence length of the microtubules, where the filament persistence length is a measure of microtubule flexural stiffness. However, previous experiments have shown that there is a significant discrepancy between the path and filament persistence lengths. Here, we showed how this discrepancy arises by using computer simulation. By simulating molecular shuttle movements under external forces, the discrepancy between the path and filament persistence lengths was reproduced as observed in experiments. Our close investigations of molecular shuttle movements revealed that the part of the microtubules bent due to the external force was extended more than it was assumed in the theory. By considering the extended length, we could elucidate the discrepancy. The insights obtained here are expected to lead to better control of molecular shuttle movements.

www.nature.com/scientificreports/ In this study, we report how the discrepancy between L p,path and L p,filament of MTs arises using computer simulation. A simulation study has distinct advantages over an experimental one. First, since in the simulation, L p,filament is a preset parameter, comparison between L p,path and L p,filament and its interpretation are straightforward and avoid complexities associated with MT preparations in experiments. Second, the precise locations of bound kinesins are easily obtained in the simulation while in experiments such information is only available under limited conditions 29 . Using these advantages, we test a hypothesis that treating kinesins as flexible anchors may alter the path persistence length with a simulation of MTs gliding over kinesins under external forces. While the path persistence length can be calculated from trajectories without external forces, the use of external forces enables explicit calculation of the length of the fluctuating part of MTs as that of the bent part because the front part of the MT is either freely suspended or loosely anchored by kinesin motors. In the following, we first check if our simulation can reproduce the path persistence length reported from experiments. Then, we observe the length of the bent part in the MTs during translocation and find that the length is more extended than that assumed in the theory. This observation resolves the discrepancy between L p,path and L p,filament of MTs.

Results
Microtubule trajectories under external force field. Following previous experiments with controlled external force 18,30 , we simulated MT movements under external forces (Fig. 1b). MTs initially moving upward were gradually biased towards the downstream of the applied force density (Fig. 2a, Supplementary Video 1). Above a certain force density depending on the surface motor density (for example, 3.0 pN/µm for 10 µm −2 ), MTs showed occasional dissociations from the kinesin-coated surface, revealed by discontinuities along trajectories. Below the force densities inducing occasional dissociation, MTs followed the trajectories of their leading tips. In this study, we focused only on paths without such occasional dissociations.
For quantitative analysis, 30 MT trajectories under various external forces were simulated. Even under the same applied force density, MT trajectories showed some variance, due to thermal fluctuation (Fig. 2b, Supplementary Video 2). The curvature of MT trajectories became pronounced when increasing the applied force density. To eliminate the variance among MTs and proceed with quantitative analysis, the 30 MT trajectories were averaged. The averaged trajectory at each external force was smooth and had a convex shape, representing progressive alignments of MT trajectories along the applied force fields. The curvature of the averaged trajectory increased with the applied force density (Fig. 2b). The averaged trajectories were compared with the following equation based on the elastic theory of bent rods obtained by van den Heuvel et al. 18 : where k B is Boltzmann constant, T is absolute temperature, f is external force density, d is the average length of the bent part of MTs. L p is either path or filament persistence length as we specify later. The averaged simulated path was fitted with Eq. (1) (Fig. 2c, Fig. S1). There were systematic deviations especially at low motor density and high force density, presumably because of detachments of kinesin motors from MTs by the external force. The radius of the curvature, R 0 , was inversely proportional to the applied force density (Fig. 2d). The proportionality depended on the surface motor density.

Reproduction of experimental path persistence length.
To check if our simulation reproduced the previous experimental results 18,30 , we calculated L p,path from Eq. (2). For this, as R 0 was obtained for various 1/f in the previous section, we needed to know the length of the bending part, d . By following the theoretical assumption 22 , we took d as the MT tip length between the MT leading end and the foremost kinesin (Fig. 3a). A representative time-evolution of the tip length of a MT is shown in Fig. 3b. The saw-toothed profile was plausible based on the following behaviour. As the foremost kinesin moved upward the MT plus end, the tip length was almost linearly elongated. Upon binding of a new kinesin at the tip, the tip length suddenly shortened. Reflect- www.nature.com/scientificreports/ ing random encounters between an MT and kinesins, the distribution of the tip length was an exponential one ( Fig. 3c), as considered in a previous study 18 . We calculated the average tip length of MT with various motor densities, σ , and applied force density, f (Fig. 3d). The average tip length d decreased with the increase of the motor density and did not change significantly with the applied force densities.
By substituting the obtained average tip length as d in Eq. (2), L p,path was calculated (Fig. 3e). We found that L p,path of MT was around 0.2-0.5 mm, which is much shorter than L p,filament of 5 mm that went into the simulation as a fixed parameter (Fig. 3e). L p,path obtained here is comparable with the experimentally reported values of 0.1-0.5 mm [16][17][18][19] . This simulation result showed that L p,path is not necessarily equal to L p,filament . It should be noted that the discrepancy between L p,path and L p,filament was reproduced without assuming the length-dependent rigidity of MT 24,25 . Length of the fluctuating part in microtubules during translocation. To test the theoretical assumption that the length of the bending part of MTs is equal to the tip length, we calculated the angular fluctuation of MT segments. Figure 4a shows the variance of angular fluctuations of MT segments as functions of the contour distance from MT leading ends (for details, see Fig. S2). Angular fluctuations were significantly larger www.nature.com/scientificreports/ at both ends of MTs than those at the central parts of MTs. Since the tip length was 0.25-0.6 µm (Fig. 4b), the length of the part showing significant fluctuation was longer than the tip length, indicating that the theoretical assumption may not hold. To see if the extended length of fluctuating parts is plausible, we deduced the length of the bending part of MTs by using the preset value of L p,filament and using the obtained R 0 . By substituting into Eq. (2) the obtained values of R 0 and other preset parameters including L p,filament and solving for d , we calculated the length of the bent part to be 1.1-1.7 μm. The calculated length of the bent part was around three times larger than the previously calculated tip length defined (Fig. 4b). The calculated length of the bent part showed good agreement with the range of the fluctuating part obtained by the fluctuation analysis (Fig. 4a).
These simulation results were caused by kinesins being flexible so that the fluctuating part of MT was extended more than the length between MT leading tip and the foremost kinesin as schematically shown in Fig. 4c. Since the average separation between neighboring kinesins on MTs was 0.35-0.67 μm ( Fig. 4d; for details, see Supplementary Information), on average, 2-3 kinesins bound to MTs within the bent part of MTs.

Discussion
In this study, we showed that L p,path of gliding MTs over kinesin-coated surfaces could significantly differ from the MT L p,filament , as opposed to the previous theoretical prediction. This difference arose because the fluctuating part of the MT during its gliding over a kinesin-coated surface extended much longer than the tip length, as previously suggested by Agayan et al. 19 . This happened because the bound kinesins did not act as rigid anchors while the theory assumed that kinesins behave as rigid anchors. Our simulation results showed that the significant difference between L p,path and L p,filament could arise without assuming the length-dependent MT bending stiffness 24,25 . Thus, the difference between the two persistence lengths can be explained without assuming a length-dependent MT bending stiffness. Measurements of MT L p,filament from gliding MT trajectories 17,18 based on the theoretical prediction that L p,path is equal to L p,filament , may need to be reconsidered.
We showed that L p,path resulted from not only the mechanical properties of MT but also the surface density of kinesins. As L p,path of gliding MTs over the kinesin-coated surface is a vital parameter for determining the performance of MS applications, the insights obtained here will be useful in designing elements for devices 31

Simulation method
The simulation method was based on our previous work 35 which we extended to include external applied force to MTs. In the following, we briefly summarize the simulation method. The trajectory of the MS was modelled as a three-dimensional movement of MTs propelled by kinesin motors. We assumed the MTs to be infinitely thin and inextensible semiflexible bead-rod polymers with a flexural rigidity of 22.0 pN µm −236 . The length of MTs was set to be 5 µm and each MT consisted of 20 rigid segments.
Microtubule movement was simulated with Brownian dynamics under the constraint of fixed segment length. In this method, a single time step consisted of the following two steps.
In the first step, the beads representing an MT were moved without considering any constraint, using the following expression: where r i is the position vector of the i-th bead consisting of a microtubule, ζ is the viscous drag coefficient, F bending,i is the restoring force of MT bending, F kinesin,i is a force exerted by bound kinesin, F ext is an external force, D is the diffusion coefficient of the bead, and ξ i is a three-dimensional random vector whose components take random values with zero mean and standard deviation of one. t was set at 0.5 × 10 -6 s to ensure numerical stability. The viscous drag coefficient used was the average of the parallel and perpendicular drag coefficients 37 : www.nature.com/scientificreports/ where η is the viscosity of water (0.001 Pa s), d is the length of the MT segment (0.25 µm), and r was the radius of MT (12.5 nm). The diffusion coefficient was calculated using D = k B T/ζ . We took this length of the MT segmentation such that taking shorter microtubule segmentation lead to negligible change of results (Fig. S4). The restoring force of MT bending was calculated from the following bending potential 37 : where EI is the flexural rigidity. Kinesin motors were randomly distributed over the allowed surface by specifying the positions of the kinesin tails (Fig. 1b). If an MT segment came close to a kinesin motor tail within a capture radius (20 nm) 22 , the kinesin motor was assumed to be bound to the MT segment, and the position of the motor head was specified on the MT segment. Once bound, the head of the bound kinesin motor moved toward the MT plus end with a force-dependent velocity expressed as where v 0 is the translational velocity without applied force, F is the component of the pulling force along the MT, and F stall is the stall force of the kinesin motors. v 0 was set at 0.8 µm/s, and F stall was set at 5 pN.
The bound kinesin acted as a linear spring between the motor head and tail with the spring constant of 100 pN/µm 38 and with an equilibrium length of zero and exerted a pulling force on the MT segment. The pulling force was divided into two forces which acted on the two beads located at either end of the MT segment where the kinesin motor was bound, under the condition that the total force and torque on the segment remained the same.
A kinesin motor bound to an MT detached when tension reached 7 pN. By following the approach taken by Gibbons et al. 39 , we neglected the spontaneous dissociation of the bound kinesin from the MT.
As the external force, a uniform force field pointing in the x-direction was applied: where, f is a force density and F ext is an external force.
In dealing with Eq. (3), we used an implicit-explicit method, where the restoring force of MT bending was implicitly calculated while other terms were explicitly calculated.
In the second step, the unconstrained movements were corrected by considering the constraints due to the segment length and the guiding tracks. To keep the segment length constant, the coordinates of the beads representing the MT {r i } as shown were subject to the following holonomic constraints: To keep the MT movement above the substrate, the position of the beads representing the MT were subjected to the following holonomic constrains: The correction was carried out with the following expression: where �r i (t + �t) is the correction term and segment,k and track,i are Lagrangian multipliers, which were determined in order for the coordinate at t + t to satisfy the constraints given by Eqs. (8) and (9), respectively. For this, we went through the calculations for the constraints one by one, cyclically, adjusting the coordinates until the constraints were satisfied with a tolerance of 10 -6 µm.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.