Effects of Obstacles on the Dynamics of Kinesins, Including Velocity and Run Length, Predicted by a Model of Two Dimensional Motion

Kinesins are molecular motors which walk along microtubules by moving their heads to different binding sites. The motion of kinesin is realized by a conformational change in the structure of the kinesin molecule and by a diffusion of one of its two heads. In this study, a novel model is developed to account for the 2D diffusion of kinesin heads to several neighboring binding sites (near the surface of microtubules). To determine the direction of the next step of a kinesin molecule, this model considers the extension in the neck linkers of kinesin and the dynamic behavior of the coiled-coil structure of the kinesin neck. Also, the mechanical interference between kinesins and obstacles anchored on the microtubules is characterized. The model predicts that both the kinesin velocity and run length (i.e., the walking distance before detaching from the microtubule) are reduced by static obstacles. The run length is decreased more significantly by static obstacles than the velocity. Moreover, our model is able to predict the motion of kinesin when other (several) motors also move along the same microtubule. Furthermore, it suggests that the effect of mechanical interaction/interference between motors is much weaker than the effect of static obstacles. Our newly developed model can be used to address unanswered questions regarding degraded transport caused by the presence of excessive tau proteins on microtubules.


Introduction
Cells use various motor proteins for active transport along the cytoskeleton. Among these proteins, kinesin-1 is responsible for the anterograde transport along microtubules (MTs). The two identical heads of kinesin-1 are connected to their neck linkers (NLs), which are composed of fourteen amino acids (AAs). NLs are wrapped around each other at the neck forming a coiled-coil structure. The neck is attached to a cargo linker, whose other end is connected to a tail domain where a cargo binds. MTs consist of numerous tubulins where the kinesin heads can successively bind in a stepping process. Kinesins repeat one mechanochemical cycle per step. At every cycle, one head moves to a different binding site approximately 16 nm away from the previous one. Thus, the cargo attached to the kinesin moves approximately 8 nm per step. The time required to complete one mechanochemical cycle is determined by the interaction of kinesin heads and adenosine triphosphate (ATP) [1][2][3][4][5][6]. The average period of the step is about 10 ms at high ATP concentrations. Hence, the velocity of kinesins is about 800 nm/s, and varies with the load acting on the cargo [4,[7][8][9]. Kinesins have a run length of about 1 μm, i.e., they can move about 1 μm before they are released from the MTs [10][11][12][13][14].
These characteristics of kinesins can be changed by other proteins bound on the MTs. Recent experiments suggest that other proteins attached to the MT surface can act as obstacles for kinesins [15][16][17][18][19][20][21]. Kinesin motion can be affected also by other nearby kinesins. Leduc et al. [15] observed that the velocity of kinesins decreases dramatically when the number of motors on a MT exceeds a critical value. Furthermore, nonmotile kinesin mutants bound on MTs decrease the velocity and the run length of other walking kinesins [17]. Moreover, the motion of kinesin is inhibited also by MT associated proteins such as tau. The effects of tau proteins on kinesin and dynein have been compared by Dixit et al. [18], where it was revealed that the dependency of kinesins on tau proteins bound to the MT is much stronger than that of dynein. Experiments with cargoes inside cells also suggest that tau proteins decrease the transport distance of molecular motors [20].
Several models have been proposed to describe the motion of kinesin heads [22][23][24][25][26][27], focusing mainly on the one dimensional motion of the kinesin head. However, several experiments with static obstacles suggest that kinesins can bypass obstacles by executing 2D movements along the MT surface [17,19,28]. To capture this 2D stepping motion, we first develop a new model to capture the diffusive motion of kinesin heads in the absence of obstacles. Various equations (i.e., Fokker-Plank equation), models (i.e., a worm-like chain model (WLC) [29]), and experimental results (i.e., unfolding of the coiled-coil structure [30]) are coupled in our model. Next, a deterministic model and stochastic model are developed by using the diffusion model to predict the motion of kinesin in the presence of obstacles.

Model
The walk cycle of kinesins is realized by their structural changes and by diffusion of its heads. It begins from a state where ATP is not bound on the fixed head of the kinesin, as described in Fig 1b1. The structural change from Fig 1b1 to 1b2 is caused by the binding of ATP to the fixed head. This binding creates the bonds between the fixed head and its NL. After the structural change, the free head diffuses until it reaches one of the eight neighboring binding sites (1)-(8), as shown in Fig 1b3. The diffusive motion of the free head is important because it determines the direction of the kinesin walk.
The motion caused by the structural changes is considered in our study by using the results of a previous study [31]. Specifically, the head which is not bound to the MT moves about 6 nm along the MT axis in a very short time immediately following a structural change. The subsequent motion (i.e., diffusion of the free head) is captured by using the developed methods explained below.

Forces acting on kinesin structure
Forces acting on kinesin structure need to be calculated in order to predict the diffusive motion of the free head. Because the free head is connected to the fixed head through two NLs which are nonlinear springs, as shown in Fig 2a, the springs have to be stretched for the free head to reach the binding sites.
The forces acting on the free head (i.e., F x and F y in Fig 2b) are calculated by considering three different forces: forces, F DCNL , acting on the docked NL, forces, F UDNL , acting on the undocked NL, and forces, F c , transferred from the cargo to the neck via the cargo linker. For a given position of the free head, the neck is located at the position where the vector summation of these forces is zero. F DCNL and F UDNL are calculated by using the WLC [29,32] as where k B is Boltzmann constant, T (= 300 K) is the absolute temperature, and ℓ p (= 0.6 nm) is the persistence length. ℓ c is the contour length of the NLs which can be calculated as N AA ℓ a , where ℓ a (= 0.4 nm) is the contour length between two adjacent AAs. N AA denotes the number of AAs in the NLs. N AA is equal to N AA, 0 + 3.5 n uw , where N AA, 0 is the number of AAs in the NLs when every bond in the neck is intact. n uw is the number of unwound turns of the neck.
Because the unfolding of one turn increases the number of AAs by 3 or 4 for each NL, the coefficient 3.5 is used when calculating N AA with n uw . N AA, 0 is 14 (AA 325-338) for the undocked NL and 4 (AA 335-338) for the docked NL because AA 334 of the docked NL is tightly bound to the fixed head with two backbone hydrogen bonds [33]. The neck consists of two coils shows the components of the kinesin structure and the procedure for its walking motion. The kinesin has one pair of identical heads and NLs. x represents the direction of MT-axis, and y is the tangential direction of the MT. xyplane represents the outer surface of the MT. The cylinder with bold outlines denotes the fixed head on the MT, and the cylinder with thin outlines is the free head. The neighboring binding sites of the kinesin are shown as gray ellipses. The plus and minus signs denote the polarity of the MT.
which are connected to each NL, as shown in Fig 1. The coils are coiled with each other for about 10 turns, and that coiled-coil structure is maintained by hydrophobic interactions [34]. In experiments performed by Bornschlogl et al. [30], it is observed that each turn of the neck can be unwound when stretching forces of about 10 pN are applied. This behavior changes the relation between the force acted in the NL and its length when the force reaches the unwinding force values, as shown in Fig 2c. The force is also affected by ℓ DC which is the distance between AA 324 and AA 334 of the NL in a condition when they are docked to the fixed head. This length determines the position of AA 334 of the docked NL, as shown in Fig 2a. Thus, the magnitude of the force in the x and y directions (i.e., F x and F y ) acting on the free head is symmetric about x = ℓ DC , as shown in Fig 2d. The refolding of the neck is also considered by assuming that the force-displacement relation for stretching is the same when releasing that stretching. When calculating F x and F y , the state is first considered where every chemical bond between two coils of the neck is connected (i.e., n uw = 0). The moment caused by F c on the neck is negligible because F c acts at the center of the neck. The two coils align along a line which is indicated with the dotted line in Fig 2b, so that the net moment caused by F UDNL and F DCNL is zero. However, both coils are stretched by the force F UW along the dotted line in Fig 2b. The magnitude of F UW can be calculated using F c , F DNL , and F UDNL . The bonds of the neck can be broken or intact depending on F UW . If F UW is less than the force required to unwind the first bond of the coiled-coil structure, then F x and F y are determined by F UDNL . If the calculated F UW is larger than the force to unwind the first bond, the simulation step is performed again, unwinding the first turn of the neck (i.e., n uw = 1). This procedure is repeated until F UW becomes less than the force required to unwind the subsequent turn.

Diffusive motion of kinesin head
To determine the direction of the kinesin step, the diffusive motion of the free head has to be considered. However, calculating that diffusive motion at every walk cycle requires intensive computations. To reduce computation effort dramatically, the following method is used. First, the probability of the direction of step is calculated by solving the Fokker-Plank equation on the position of the free head. Then, a random number is generated at every step of kinesins. Next, the direction of each step is determined by comparing the random number to the calculated probability. Because a solution of the Fokker-Plank equation can be used at every step of kinesins, this method enables us to determine the stochastic two dimensional walking motion of kinesin with higher computational efficiency.
The probability density for the position of the free head is calculated over time by solving the following Fokker-Plank equation (in a domain depicted in Fig 1a with where x and y are coordinates of the free head with respect to the fixed head along the MT-axis and along the tangential direction of the MT, respectively. The software COMSOL is used to solve this partial differential equation with the finite element method. p is the probability density function of the position of the free head. D is the diffusion coefficient, and γ is the drag coefficient of the free head in water. The head is assumed cylindrical with a length of 7 nm and a diameter of 4.5 nm. These dimensions of the head are obtained from a previous study on the crystal structure of kinesin heads [35]. γ is calculated to be 5.625 × 10 −8 g/s by using the method of Swanson et al. [36] for the drag coefficient for a cylinder moving at low Reynolds numbers. Then, the diffusion coefficient can be calculated as 73.6 μm 2 /s for a temperature of 300 K. Note that the dependency of γ and D on the direction of motion is negligible for such cylinder dimensions [36]. The effect of the motion along the radial direction of the MT is assumed to be negligible. The probability of stepping to the left sites (i.e., site (4) or (6) in Fig 3) can be different to the probability of walking to the right sites (i.e., site (5) or (8) in Fig 3) because there is a mismatch in the lattice of MTs due to a variable number of protofilaments. However, Yildiz et al. [37] observed that the probability of kinesin to move left is the same as the probability to move right. This observation suggests that the effects of the mismatch in the lattice of MTs are not significant. To confirm that this effect is small, the distance from the binding sites to the center of the fixed head is calculated. For a MT diameter of 26.5 nm and a number of 13 protofilaments, the angle of the mismatch is approximately 5.49°(= atan(8/(26.5 × π))). For a helical MT structure, the distances to the two diagonal sites (left and right) are 10.74 nm and 9.77 nm, respectively, as shown in Fig 3a. In the absence of a mismatch in the lattice array, these distances are both 10.25 nm. Because the change in the distance due to a lattice mismatch is not considerable (i.e., only about 5%), isotropic diffusivity can be used in this study.
Two types of boundaries are considered, as shown in Fig 4a. Absorbing boundary conditions (i.e., p = 0) are applied at the binding sites. Reflective boundary conditions (i.e., no flux along the normal direction to the boundary) are used for the other boundaries. Fig 4b shows one example of the spatial probability density for the free head. The probability that the free head binds to a certain binding site is obtained by integrating over time the flux _ P of the probability flowing out through the boundaries of the site (which is shown in Fig 4c).

Binding with tilted posture
To consider the effects of the affinity of the kinesin head to the binding sites, the positions of the specific AAs of the free head can be examined. First, AA 142-145, AA 273-281, AA 238 and AA 255 are responsible for binding on the MT [33]. Thus, it is assumed that the head binds to a binding site when the geometric center of these AAs reaches the binding site. Second, AA 324, which connects the NL to the head, is located about 0.9 nm away from that geometric center along the x direction and about 1.8 nm away along the y direction, as shown in Fig 5a. Due to this structure, when the free head is located at a diagonal site with a tilted posture, the length of the NLs to reach that site is shorter than the length when the free head reach the site with no rotation, as shown in Fig 5b. Thus, it is assumed that the binding to the diagonal site with tilted posture is favorable to the kinesin free head. The parameter θ d represents the allowable tilting angle. We further hypothesize that if the free head is tilted with an angle larger than θ d , then the affinity between the binding site and the AAs of the free head which interacts with the binding sites is not strong enough to cause binding. When the geometric center of AAs responsible for binding is located at the diagonal binding site, and the tilted angle is less than θ d , the interaction between the head and the site is strong. Then, the head binds to the site. Subsequently, the head is aligned along the MT by that interaction. To mathematically model this behavior, the positions of the absorbing boundaries in the domain are calculated as follows. First, we locate the geometric center of the AAs responsible for binding of the free head at a certain binding site (e.g., site (8)), as shown in Fig 5b. Then, we rotate the free head with respect to the geometric center by an angle θ d . Then, the position of the binding is determined as the relative position of AA 324 of the free head with respect to the position of AA 334 of the fixed head. For example, when θ d is 47°, the distances between the diagonal site and the fixed head decreases by 14.7% compared to the same distance when θ d is 0°. Fig 5c shows the position of binding sites obtained using this method. Note that this behavior is only applicable to the diagonal sites because there is not enough space for the free head to tilt when the head is near the side, forward, or backward binding sites due to the geometric interference between the free head and the fixed head.

Comparison to experimental data
The diffusion model has three parameters; ℓ DC , θ d , and d side , where d side is the distance between two adjacent binding sites along the tangential direction of the MT. The values of these three parameters are determined by calculating the probability (P b, i ) of the free head to bind at each neighboring binding site and by comparing it with experimental data. Yildiz et al. [37] tracked the motion of kinesin heads by labeling them with quantum dots. Their results show that 13% of the kinesin steps involve motion along the tangential direction of the MT, and 70% of that lateral motion occurs together with the motion toward the plus end of the MT. Also, the measurements of Nishiyama et al. [38] show that the probability of backward steps increases exponentially over the resisting load acting on the cargo. According to their observations, the probability ratio of forward and backward steps is about 1.3 when the resisting load is about 7 pN. A set of parameters (ℓ DC = 2.9 nm, d side = 6.4 nm, and θ d = 47°) satisfies the probability of step observed in both experiments. The probabilities of forward, backward and sideway step in the absence of resisting loads are presented in Fig 6a. ℓ DC and d side are similar to the measured values of the previous studies. The value of ℓ DC in the crystal structure of kinesin [39] is about 3 nm. Also, d side of 6.4 nm can be obtained when the diameter of the MT is 26.5 nm. This diameter value is in the range of the actual diameters of MTs observed in the cell [40,41].

Effects of load on the probability of step
The external load acting on the cargo affects the probabilities of kinesin step. When an external load is applied to the cargo, it is transferred to the cargo linker. That force, F c , changes the force acting on the free head (i.e., F x and F y ). As a result, the probabilities of sideway and backward

Effects of unwinding the neck and binding with a tilted posture
Both the unwinding of the neck and the binding with tilted posture are necessary to obtain the probability of the direction of a step measured in experiments [37]. To check the effects of unwinding the neck on the diffusion, P b, i is obtained using the model which allows the binding with tilted posture but does not allow neck unwinding. U off + T on in Table 1 represents this case. The equation used to calculate F NL is the same as Eq 1 except that the value of N AA is fixed as N AA, 0 . Similarly, U on + T off represents the case of constraints which do not allow unwinding of the neck but allow binding with tilted posture. U on + T on represents the case where both unwinding of the neck and binding with a tilted posture are allowed. The binding probability presented in Table 1 suggests that both behaviors are necessary to obtain values similar to those P b, i measured experimentally [37].

Number of binding sites occupied by kinesin
The number of binding sites occupied by kinesin molecules depends on the chemical state of its head. The kinesin head has a strong affinity to the MT when it has ATP or no nucleotide [42]. The head has a low affinity to the MT if it has adenosine diphosphate (ADP). When one head has ATP and the other head has no nucleotide, as described in Fig 1b3, both heads are strongly bound to the MT. Thus, two sites are occupied by the kinesin in this state, as shown in Fig 7a. When one head has no nucleotides or ATP and the other has ADP, as shown in Fig 1b1 or  1b2, one head diffuses around the (other) head which is strongly bound to the MT [43,44]. For these states, the strongly bound head (i.e., the head with no nucleotide or with ATP) occupies a The probability P b of backward steps is not indicated because its value is very small. (b) shows the changes in the probabilities of forward (P fw = P b, 6 + P b, 7 + P b, 8 ), sideway (P sd = P b, 4 + P b, 5 ), and backward (P bw = P b, 1 + P b, 2 + P b, 3 ) steps over the external load.

Results
The motion of kinesin is affected by the interaction between obstacles and kinesins and by the number of the obstacles on the MT. Our model is used to characterize the interaction by locating obstacles near a kinesin molecule. By using experiments, where the velocity and run length of kinesins are observed in the presence of immobile kinesins, the value of parameters regarding the interaction between kinesin molecules are obtained. Then, the velocity and run length of kinesins are calculated over the number of obstacles. Furthermore, the characterization of interaction between kinesins and obstacles enables us to predict the motion of kinesin when several types of molecular motors move along the same MT. In addition, the effects of the length of NLs and the probability of large steps (i.e. two binding sites away from the fixed head) are studied.

Motion of kinesin in front of a single obstacle
To characterize the interaction between kinesins and their neighboring obstacles, the size of the region which the kinesin head cannot reach due to the obstacle (R obs ), the number of binding sites occupied by a single obstacle (m obs ), and the degree of interference (A int ) between the kinesin and the obstacle are considered. When a binding site is occupied by obstacles, the free head is not able to reach it. Thus, a reflective boundary is located at that binding site. The shape of the reflective boundary is assumed to be circular, as shown in Fig 8a1. Next, the degree of interference (A int ) between obstacles and the free head is calculated by using the probability density of the free head to reach the boundary as where t 0 and t f are the initial and final time, respectively, when solving Eq 2. C obs is the reflective boundary of obstacles located near the kinesin. Then, the unbinding probability of kinesin per step can be calculated by using A int as where P 0 ub is the unbinding probability per step when every neighboring binding site is not occupied by obstacles. P 0 ub are obtained from our previous model [45]. ΔP ub, obs is the change in the unbinding probability due to obstacles. Note that the kinesin dissociates from the MT mostly when one of its heads is diffusing [10,45]. Thus, it is assumed that the increase of the unbinding probability due to obstacles during that diffusion is also very large compared to ΔP ub, obs for other kinesin states. α obs is the parameter of the model which calculates the probability to unbind when the free head interacts with the obstacles.
If the site ahead of the fixed head is occupied by a the small obstacle, the kinesin is likely to bypass it by moving forward to a diagonal site. However, for a large obstacle, the probability to move to one of the the side sites is larger than that of moving to a forward diagonal site, as shown in Fig 8a2. Thus, the kinesin needs a larger number of steps (n step ) to bypass an obstacle when R obs is large, as shown in Fig 8b2. n step also increases with m obs exponentially. Note that A int is maximum when the forward site is occupied by an obstacle, as shown in Fig 8c2. This means that the interference of obstacles is most significant when they are located in front of kinesins.

Static obstacles
The parameters describing the effects of obstacles (i.e., R obs , m obs , and α obs ) depend on the type of obstacle. In this study, immobile kinesins which act as static obstacles for walking kinesins are characterized by using previously reported experimental data [17]. In those experiments, the number of obstacles (i.e., immobile kinesins) attached to the MT was 8% of the maximum number of obstacles (obtained when the MT is fully coated with the obstacles). The velocity and run length were observed to decrease to 81% and 57% compared to their values in the absence of obstacles. The values of R obs and m obs are obtained by using the changes in the kinesin velocity. The observed decrease of the run length is used to calculate the value of α obs .
Velocity in the presence of static obstacles. The motion of the heads of the immobile kinesins is the same with that of the walking kinesin when one of its head is bound to the MT  7b2, 7b3 and 7b4, respectively. The behavior of the unbound head, which can share binding sites with other kinesins, is also considered as described in Fig 7c1. To estimate the number of immobile kinesins present in the experiment, we randomly distributed them until the MT was fully coated with immobile kinesins. When the MT is saturated with immobile kinesins, the molar ratio ρ of the immobile kinesins and tubulin dimers is calculated as 0.4352 for m obs = 3, 0.3572 for m obs = 5, and 0.1843 for m obs = 9. When we compared the velocity predicted by the model to the velocity observed in the experiments, the molar ratio ρ of 0.08 is used for m obs = 1, 0.0346 (= 8% × 0.4352) is used for m obs = 3, ρ of 0.0286 (= 8% × 0.3572) is used for m obs = 5, and ρ of 0.0147 (= 8% × 0.1843) is used for m obs = 9.
To obtain the value of aforementioned parameters, we developed a mathematical equation capable of calculating the (average) velocity and run length for various numbers of obstacles by using the behavior of kinesin in front of single obstacles, as explained previously. Specifically, velocities can be calculated as where V(ρ) is the velocity for a molar ratio ρ. V 0 is the velocity in the absence of obstacles. The derivation of this equation is provided in S1 File. The coefficient a k is determined by n step , which is calculated using the diffusion model. Table 2 shows the velocities obtained for R obs = 4 and 5 nm and m obs = 1, 3, 5, and 9. The values of these parameters were determined as R obs = 5 nm and m obs = 3 for the immobile kinesin obstacles. Note that m obs = 3 is also used for the walking kinesin when one of its head is fixed and the other is not bound. Run length in the presence of static obstacles. The run length of the deterministic model can be calculated as where RL(ρ) and RL 0 are the run length in the presence and absence of obstacles, respectively. The values of b k depend on α obs . The derivation of this equation is provided in S1 File.
To obtain precise values of the parameters, the time and spatial resolutions used in experiments were applied to our stochastic model. The mathematical models in Eqs 5 and 6 are based on the assumption that the time and spatial resolution is sufficiently small. However, a finite resolution is used in the experiments. Thus, the velocity and run length are obtained again by  [45]. Details on the stochastic model are described in S1 File. The stochastic model with R obs = 5 nm and m obs = 3 provides velocities similar to those obtained using Eq 5. The value of α obs which leads to a good match with experimentally observed run lengths is of 0.044 nm/ms. It is worthy to note the effect of obstacles on the run length and velocity of kinesin. The run length decreases as the unbinding probability is increased due to obstacles. The effect of unbinding probability on the velocity depends on the concentration of obstacles anchored on the MT. When the concentration of obstacles on the MT is small, the velocity does not change. However, if the number of molecules is considerable, both the size of obstacles and the degree of interference (A int ) should be considered. If the size of obstacles (which is considered by using R obs ) is small, the probability regarding the direction of step does not change significantly. Thus, the change in velocity is very small. However, if R obs is large and A int is small, the kinesins have more chances to walk sideways. In this case, the velocity can be reduced over the concentration of obstacles. If both R obs and A int are large, both the probability to walk sideway and the unbinding probability increase. As a result, the kinesin is likely to unbind from the MT near obstacles before walking to side binding sites. Thus, the decrease in the velocity is not considerable.

Moving obstacles
Since the state of immobile kinesins is the same with the state of walking kinesins when they wait for ATP to bind, the interactions between moving kinesins and immobile kinesins, which is characterized in this study, can be used to consider the interaction between moving motor proteins. The minus-end directed motors, such as dyneins, are considered in our model simply by reversing the direction of the steps of kinesin. (K + M F ) represents the situation when the motion of kinesin is disturbed by immobile kinesins. Also, kinesins can walk while surrounded by other kinesins (K + M + ) or in the presence of minus-end directed motors (K + M − ). These cases are depicted in Fig 9a. The changes of the kinesin motion by surrounding motors are related to the probability of the kinesin to interact with other motors. That probability is high for (K + M F ), and is small for (K + M + ) due to the motility of the surrounding motors. Thus, Effects of Obstacles on the Dynamics of Kinesins Predicted by a Model of Two Dimensional Motion (K + M + ) has the smallest reduction in velocity and run length, whereas K + M F causes the most significant decreases in velocity and run length, as shown in Fig 9b and 9c.

Length of NLs
The proposed model can be used to study the effects of the numbers of AAs at the NLs. This changes the probabilities of sideway and diagonal stepping and the interference with obstacles. These quantities are predicted by using our model. For example, Hoeprich et al. [46] predicted that the probability of stepping to each binding sites is slightly altered depending on the number of AAs in the NLs. Our model provides consistent results with their findings when the number of AAs in the NL is changed from 14 to 17, i.e., changed to the number of AAs in the NLs of kinesin-2. Our model predicts that when NLs of kinesin-2 are used, the probabilities of sideway and diagonal stepping increase by 1% compared to the same probabilities for kinesin-1.
Additionally, according to our model, when a single obstacle of R obs = 5 nm is located in front of a kinesin (i.e., site (7)), the value of A int is 0.46 ms/nm for kinesin-2, which is smaller than the respective one (i.e., 0.58 ms/nm) for kinesin-1. Hence, the interference between the free head of kinesin-2 and obstacles is not as strong as the interference between kinesin-1 and obstacles. This means that kinesin-2 is less likely to unbind from the MT in the presence of obstacles than kinesin-1. Therefore, the model predictions are consistent with the measurements of Hoeprich et al. [46] who observed that the run length of kinesin-1 decreased due to the presence of tau proteins, whereas and the run length of kinesin-2 did not decrease significantly.

Probability of large step
Another use of our model includes studying calculating the probability of larger steps (i.e., probability of stepping to distant binding sites). The domain, which is shown in Fig 4a, is extended to consider P b, i of twenty four neighboring binding sites. Also, the forces acting on the head (i.e., F x (x, y) and F y (x, y)) are calculated in this extended domain. Next, the diffusion of the free head is obtained, as shown in Fig 10. The result suggests that binding to sites far from the fixed head is possible but the probability of such events is very small. The sum of the probabilities (= P b, 9 + P b, 10 + P b, 11 + . . . + P b, 24 ) is about 1.3 × 10 −8 .

Discussion
A diffusion model was developed in this study to predict two dimensional motion of kinesin which is necessary to bypass obstacles. Next, the interference between the walking kinesin and immobile kinesins was characterized by comparing the kinesin motion obtained from the diffusion model with experimentally measured velocities and run lengths. By using that interaction, the effects of other motor proteins on kinesins are predicted when they move on the same MT. If several kinesins share one MT, their velocities are almost not affected by other kinesins on the MT. The deceleration of kinesins by minus-end directed motors is also not considerable. These results suggest that cells can handle several intracellular transport processes effectively with a small number of MTs. However, both the velocity and run length decrease considerably in the presence of immobile kinesins. As a result, the malfunction of a few motors could degrade the whole transport system.
The diffusion model developed to capture the interaction between kinesins (or between kinesins and other motor proteins) can be also applied to other types of obstacles, such as tau proteins which are attached to the MT. The run length of kinesins has been observed to decrease for an excessive density of tau proteins in several experiments [18,20,47,48]. Also, experiments have shown that the motion of single kinesins (with no cargo) is not slowed down by tau proteins. However, the transport velocity of cargoes can be modulated by tau proteins. Large cargoes such as mitochondria are transported by several kinesins to overcome large viscous forces in the cytoplasm. It is observed that the unbinding rate of kinesins increases and the rebinding rate of kinesins to the MTs decreases in the presence of tau proteins [18,[47][48][49]. Thus, the number of kinesins which effectively transport a large cargo is decreased by tau proteins, leading to slower transport of large cargo. Also, this decelerated cargo transport and the attached kinesins can behave as obstacles for other kinesin molecules. As a result, the entire transport system cannot operate properly.
At this point, it is worth to note that the effects of obstacles on the motion of kinesin is important because several neurodegenerative diseases such as Alzheimer's, Parkinson's, and Huntington's disease have been reported to be related to an abnormal intracellular transport [50][51][52]. The negative effects of obstacles on the intracellular transport have been observed also in vivo. Specifically, the transport of various types of cargoes was decelerated or failed due to an excessive number of tau proteins [53][54][55]. Thus, the characterization of the motion of kinesin in the presence of tau proteins is necessary to reveal quantitative aspects of the relation between neurodegenerative diseases and tau proteins. Our diffusion model can be effectively used to obtain this kind of information and thus shed light on the relation between excessive tau proteins and degradation of axonal transport.
Supporting Information S1 File. Deterministic and stochastic model for the motion of kinesin in the presence of obstacles. Fig A describes several situations where a kinesin molecule encounters obstacles. Table A presents the parameters to calculate velocity using the deterministic model. The values of parameters of the mechanistic model are presented in Table B. (ZIP)