Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Phys Biol. Author manuscript; available in PMC 2010 April 21.
Published in final edited form as:
PMCID: PMC2858005

Robust transport by multiple motors with nonlinear force–velocity relations and stochastic load sharing


Transport by processive molecular motors plays an important role in many cell biological phenomena. In many cases, motors work together to transport cargos in the cell, so it is important to understand the mechanics of the multiple motors. Based on earlier modeling efforts, here we study effects of nonlinear force–velocity relations and stochastic load sharing on multiple motor transport. We find that when two or three motors transport the cargo, then the nonlinear and stochastic effects compensate so that the mechanical properties of the transport are robust. Similarly, the transport is insensitive to compliance of the cargo-motor links. Furthermore, the rate of movement against moderate loads is not improved by increasing the small number of motors. When the motor number is greater than 4, correlations between the motors become negligible, and the earlier analytical mean-field theory of the multiple motor transport holds. We predict that the effective diffusion of the cargo driven by the multiple motors under load increases by an order of magnitude compared to that for the single motor. Finally, our simulations predict that the stochastic effects are responsible for a significant dispersion of velocities generated by the ‘tug-of-war’ of the multiple opposing motors.

1. Introduction

Motor proteins are remarkable molecular machines able to transform chemical energy into movement and force generation [1]. They use linear polar actin and microtubule filaments as ‘tracks’ to transport vesicles and organelles in the cell [2, 3], as well as to produce stresses and deformations, most notably in mitosis [4], cytokinesis [5] and cell motility [6]. Over the last two decades, the mechanics of many single-molecular motors was understood quantitatively with the help of both experimental use of optical traps [7] and theoretical applications of statistical mechanics [8].

These studies uncovered that the mechanical behavior of the individual motors can be described by force–velocity and force–processivity relations, namely by how fast a motor moves on average against a mechanical load, and how far would the motor advance before dissociating from its track. For example, the best studied microtubule-based motor, Kinesin-1, is processive, able to move freely as fast as ~1 μm s−1 and as far as ~1 μm before detachment [9]. Against the force, however, this motor slows down [10] and dissociates sooner [11].

Experimental considerations suggest that the cell can achieve moving cargos by great distances against significant resistance by using multiple motor copies [12]. Indeed, the load would be distributed among the engaged motors, so each motor can move faster against a smaller force, and if some motors dissociate, others stay engaged moving the cargo. In vitro experiments confirm greater processivity and higher stall forces of collective Kinesin-1 [9] and cytoplasmic Dynein [13] transport. Controlling the motor numbers and forces is very difficult experimentally, so mathematical modeling is very useful in interpreting the data.

A number of earlier elaborate models of mechanochemical cycles of coupled motors [1423] predicted complex, sometimes counter-intuitive, behaviors, such as bidirectional movements and oscillations, as well as effects of strain- or stress-dependent ATP binding and hydrolysis, flexibility of the cargo and cargo-motor-links, and spatial distribution of the motors on the cargo. One recent model [24] was especially instrumental in examining the collective motor mechanics by using the mean-field approximations, namely, assuming that all motors share the load equally. Based on this assumption, the authors of [24] computed effective force–velocity and force–processivity relations for N motors.

However, there are two factors, not considered in [24], that can have a significant impact on the collective motor behavior. First, because the stalks connecting motor heads and their cargo-binding domains are flexible and timing of the motors’ steps is largely random, there are stochastic fluctuations in individual motors’ positions and resulting forces applied to the motors, as well as related correlations between the coupled motors. Second, only the simplest force–velocity relation of individual motors, such that the motor velocity decreased linearly with the hindering load, was considered in [24]. Such a linear force–velocity curve is usually assumed in theoretical models [25, 26] and was in fact observed at moderate loads for the mitotic kinesin motor Eg5 [27]. Most of the measured force–velocity relations, though, are nonlinear. Most notably, a super-linear force–velocity curve is firmly established for Kinesin-1 [10], so that the velocity is less sensitive to the force at low loads and decreases rapidly with force at loads close to the stall. Mathematically, such a force–velocity curve is concave up. Similar super-linear force–velocity relations are reported for an actin-based Myosin-V motor [28]. On the other hand, sub-linear (decreasing rapidly with the force at low loads and less sensitive to the force at loads close to the stall, which mathematically has the form of the convex up curve) force–velocity relations are known for polymerizing microtubules [29] and for Dynein [13, 30], with the caveats that the former is a very special ‘one-shot’ motor, and that other data suggest a super-linear relation for the latter [31].

We started to take these factors into account in our recent paper [32], where a very specific Kinesin-1-like force–velocity relation and non-monotonic dissociation rate as a function of load were considered for small number of motors. Here, we systematically examine the force-dependent velocity and run length of N motors characterized by the nonlinear force–velocity curves. To take into account the uneven load sharing and correlations, we use the Monte Carlo stochastic simulations. In what follows, we describe the mean-field model of [24] for the nonlinear relations and the stochastic model. Then, we report the results of analytical solutions of the mean-field model and simulations of the stochastic model. In short, we find that the nonlinear and stochastic effects combine into surprisingly robust, ‘linearized’ collective behavior of the small number of motors, and that the mean-field approximation is valid for more than four motors, in which case simple scaling of the force–velocity and force–processivity relations emerges. We also report the significant increase in the effective diffusivity of the collective transport and ‘velocity-smearing’ due to the stochastic effects in the tug-of-war between multiple opposing motors. Details of the simulations are given in the appendix.

2. Mathematical model of the collective motor transport

2.1. Model of the single motor

All ever observed force–velocity relations can be well approximated by the following mathematical expression [33]:


where ν is the load-free gliding rate of the motor, F is the load force and Fs is the stall force. If the motor moves in steps of the length d, then the motor can be described effectively with the load-dependent rate of stepping:


Following [24], in the model we use Kramers’s theory formula for the force-dependent rate of motor detachment:


where ε is the load-free dissociation rate, and Fd is the characteristic force of detachment. This simplest force dependence of the motor dissociation rate is used in tens of recent models, i.e. [34].

Equations (1) and (3) for ω = 1 describe the linear force–velocity relation and exponentially decaying run length that has been used in the published mean-field model [24]. Here we examine the collective transport for other values of parameter ω. For ω < 1, equation (1) describes a sub-linear motor characterized by the convex up force–velocity curve, while ω > 1 corresponds to a super-linear motor with the concave up force–velocity relation. The sub-linear motor velocity decreases rapidly at small forces and becomes force insensitive at loads close to the stall, while the super-linear motor velocity is insensitive to the small loads and decreases rapidly at greater loads.

2.2. Mean-field model

According to the mean-field model [24], a cargo particle is transported cooperatively by N molecular motors along a filament. The model is based on the assumptions that the motors share the load equally, so that the force F/i is applied to each of i engaged motors. The state of the system is characterized by the number i of engaged motors pulling on the cargo. In this state, the cargo has the velocity


and the number of engaged motors increases with the rate


and decreases with the rate


where π is the attachment rate per motor, and N is the total motor number on the cargo.

We use the stationary solutions of the master equation obtained in [24] expressing the probability for the system to have i engaged motors in terms of the motor parameters,


to find analytically the average velocity of the cargo,


and the average run length before detachment of the cargo,


In order to compare our results with the mean-field model [24], we use the same parameter values: ν = 1 μm s−1, ε = 1 s−1, π = 5 s−1, Fs = 6 pN and Fd = 3 pN, unless otherwise stated. A number of rate values and assumptions of the mean-field model (velocity at zero load is independent of the motor number; association rate is proportional to the motor number) were confirmed experimentally in [35].

2.3. Stochastic model

To model the multiple motor transport, we place N motors on the cargo, so that each motor head is attached to a single spot by a link (figure 1). According to recent measurements [36], the link between the motor domain and the bead is highly nonlinear: when stretched beyond the rest length, it behaves as a relatively stiff linear spring characterized by the spring constant ~0.3 pN nm−1. However, its effective compressional rigidity is very low, ~0.05 pN nm−1, i.e. the link buckles almost without resistance when compressed [36]. Thus, we model each link as a linear spring exerting restoring force when stretched beyond the rest length and not generating any force if the distance between the bead attachment point and the motor head is less than the rest length. The model is one dimensional, so all distances are measured along the motors’ track. In the simulations, we use the appropriate values of linkage stiffness 0.32 pN nm−1 measured for single kinesin-1 motors in vitro [11, 36] and rest length 0.11 μm [11]. Each dissociated motor binds to the track with the constant on-rate, and each engaged motor detaches with the rate given by equation (3) dependent on the instantaneous force applied to this motor by the elastic link. Each motor makes a forward step with the force-dependent rate (equation (2)). (Probabilities of binding, unbinding and stepping events are computed by multiplying the respective rate by the time step.) The instantaneous position of the cargo is calculated at each step from the requirement that the total force on the cargo from the load and all elastic links is equal to zero.

Figure 1
Scheme of the stochastic model: the cargo (bead) has three motors attached to it. The cargo moves forward with velocity ν and is slowed down by the backward load F. At the depicted moment, two out of three motors are bound to the track, and the ...

The Monte Carlo simulations of such a single motor correctly reproduce the deterministic force–velocity and force–processivity curves for all tried nonlinearity parameter values (ω = 0.25, 0.50, 1.00, 2.00 and 4.00). We ran the stochastic simulations for two to four motors. Initially, all motor heads were attached to the track and placed at the origin, and each simulation ended when all motors dissociated from the track. The simulations were run repeatedly for various values of the load force applied to the cargo, until reliable statistics of the average velocities and run lengths were gathered. We emphasize that in this model the motors share the load stochastically, unevenly, unlike in the mean-field approximation: a ‘leading’ motor experiences the greatest hindering load, while the motor ‘lagging behind’ has a possibility of being pulled forward. We assume that under the influence of the forward load, the motor steps forward as if it was unloaded, but the forward load has the same quantitative effect on detachment as that at the backward load. (The force value in the exponent of Kramers’s formula is the magnitude of the load applied to the motor.) Details of the simulations are described in the appendix. Also, we demonstrate in the appendix that the detailed mechanochemical cycle of individual motors is not crucial for the collective motor behavior on the scale of microns and seconds.

3. Results

3.1. Nonlinear effects in the mean-field approximation

We used the formulas of the mean-field theory [24] and nonlinear force–velocity relations to obtain the force dependence of the average collective motor velocity and run length given by equations (8) and (9), respectively. The results are shown in figure 2. The run length exponentially increases with the number of motors and exponentially decreases with the external load. As expected, the multiple motor system performance improves with increasing parameter ω: the super-linear motors have the greater run length and velocity at any given load than linear motors, and those, in turn, perform better than the sub-linear motors, simply because the time to unbinding is independent of the velocity and therefore of parameter ω, so the increased/decreased run length simply reflects the fact that motors with ω > 1/ω < 1 are faster/slower.

Figure 2
Force–processivity (A) and force–velocity (B) curves (solid lines and symbols) obtained for multiple motors using the mean-field model for the force–velocity relations characterized by (a–b1) ω = 2 and (a–b2) ...

Note that the collective force–velocity properties of the super-linear motors are improved more significantly for greater loads. Also, ‘kinks’ (points where the velocity is not a smooth function of the force) in the nonlinear force–velocity curves appear at the stall forces for one and two motors, and can be explained as follows. For example, for two motors, if the total load is less than Fs, then each motor moves continuously against either half or total load, depending on whether one or both motors are engaged. On the other hand, if the total load is greater than Fs, then either both engaged motors move continuously against half load, or, if one of the motors dissociates, the remaining motor is stalled completely.

3.2. Stochastic effects lead to worse performance and ‘linearization’ of the force–velocity relations

Simulations of the stochastic model, in which the individual motors are characterized by the same parameters as those in the mean-field model, illustrate that the stochastic effects worsen the performance of a small number of motors. The results plotted in figure 3 demonstrate that both run lengths and velocities of two or three motors are lower than those predicted by the mean-field theory for any given force and for any nonlinearity (parameter ω). Same is true for groups of four or more motors (data not shown). The stochastic model predicts especially great drop in the run length in the high load regime. We also find that for any given parameter ω, having more (three instead of two) motors on the cargo does not help much to increase the processivity in the higher load limit or velocity for any load. It is as if just one motor takes on most of the load, and additional motors help very little. Interestingly, two or three motors make velocity even lower than that of one motor for low loads, but higher at greater loads. This causes ‘linearization’ effect for the small motor number: the resulting collective force–velocity relation is closer to the linear one despite varying nonlinearity in the individual motor properties.

Figure 3
(AC): Average run lengths (a–c 1–2) and velocities (a–c 3–4) as functions of the load obtained from the simulations of the stochastic model for ω = 2 (A), ω = 1 (B), ω = 0.5 (C), for two ...

The reason for the poor performance of the multiple motor system under high load in the stochastic model has its origin in the stochastic load sharing, as illustrated in figure 4. Indeed, when one of the motors advances leaving other motors behind, this leading motor takes on a disproportionate high load and detaches with a greater rate than the other motors. This increases the load on the remaining motors and, more importantly, the load causes the cargo to rapidly retract at the moment of the leading motor dissociation, stretching the remaining motor links supporting this new excessive load. The combined effect of the increased detachment rate of links’ stretching of the remaining engaged motors and of the backward cargo excursions reduces the overall run length and average velocity of the cargo.

Figure 4
(A) Average force against which the motors step (a1), average backward excursion length (a2), average frequency of the backward excursions (a3), average separation between motor heads (a4) as functions of ω for two motors with load 1 pN applied ...

These arguments are illustrated by figure 4(B) showing separation between the two motors driving the cargo, which is greater for the super-linear motors because one of these motors can more easily advance without slowing down against moderate force. For the sub-linear motors, the separation between motor heads cannot grow too much because if one of them moves ahead of another, the leading motor is slowed down significantly. The super-linear motors step against a greater average load (figure 4(A)), so they detach faster, which in turn increases the frequency of the backward excursions (figure 4(A)). Note that the motors rarely share the load equally, but rather step either against almost zero load (lagging behind a motor), or almost maximal load (leading motor) (figure 4(C)). Note also that each of the two super-linear motors moves against the average force equal to the half-total load, while the sub-linear motors advance against a lower than half total load (figure 4(A)): in the sub-linear case, mostly the rearward motor advances, while the forward motor stays put ‘subsidizing’ the rearward motor ‘catching up’.

Interestingly, the average backward excursion length is roughly constant, independent of the nonlinearity (value of ω) (figure 4(A)). The reason is that when the motors are super-linear, the spatial separation between them is significant, and when the trailing motor detaches, the cargo does not move much, still supported by the leading motor. When the leading motor detaches, the cargo makes a long backward travel. When the motors are sub-linear, they are not separated to a great extent, so when either leading or trailing motor detaches; the remaining motor link stretches to a moderate extent. Thus, the super-linear motors detach more frequently, but the long backward cargo excursions are less frequent than that for the sub-linear motors; in the latter case, these excursions are also shorter.

3.3. Robust collective motor behavior

We have investigated the effect of the motor link stiffness on cargo transport in the stochastic model. The simulations demonstrated (data not shown) that neither run length, nor velocity is sensitive to the elasticity of the links between the motors and cargo, with one exception: the groups of stiffer motors stall at slightly higher forces. This conclusion is different from that in [32], where it was suggested that stiffer links make the collective transport more effective. This effect was due to ‘strain-gating’: for example, in the case of two motors, if the motors start from the same location, one of the motors steps against the half-load, and then takes on the greater load share helping the second motor to step against the less than half-load. This leads to a more effective transport, but only when the link stiffness is high enough for the leading motor to take on a significantly greater load. However, this mechanism only works if the leading motor does not detach too frequently. In [32], a non-monotonic force dependence of the dissociation rate was used, so that the detachment became infrequent at stall. We use the dissociation rate exponentially increasing with the load, which cancels the strain-gating effect. This result emphasizes the significance of the force dependence of the dissociation rate, which is not accurately measured. We have also investigated the effect of the detachment force, Fd, on the cargo transport in the stochastic model. Not surprisingly, motors with higher detachment force perform better, as they detach less frequently at the same loads, and more attached motors transport the cargo faster (fewer backward excursions) and to longer distances.

The mean-field model [24] can be applied to the in vivo transport when the cargo experiences the viscous load. The cargo of radius r, driven by i engaged motors through a medium with viscosity η, moves with velocity vi and experiences the viscous resistance γνi, where γ is the viscous drag given by the Stokes formula: γ = 6πηr. To account for the viscous drag on the cargo transport, parameter F in equation (4) has to be replaced by γνi. This results in a linear velocity equation for ω = 1:


the solution of which has the form


In addition, one has to replace the parameter F in equation (6) by the expression γνi:


Note that equation (11) predicts that when viscous load becomes significant, the cargo velocity is proportional to the number of motors, which was observed [37]. One can now use equations (11) and (12) to calculate the average run length and velocity of the cargo as functions of the viscosity. In the nonlinear case, after replacing the parameter F in equations (4) and (6) by the expression γνi, we obtain the generalization of equation (10):


For ω ≠ 1 equation (13) has multiple roots; however, only one root lies in the physical range (0 ≤ νiν). Using solutions of equations (12) and (13) for a given parameter set (ω, r and η), we calculated the average run length and velocity as functions of the viscosity η. Respective velocity plots are shown in figure 5(A). Both velocities and run lengths (the latter are not shown) for any given ω remain almost unaffected by increasing viscosity up to ~0.01 Pa×s. At higher viscosity, both velocities and run lengths start decreasing with increasing viscosity. The velocities are insensitive to the motor number or to the nonlinearity up to viscosity ~1 Pa×s, though more super-linear motors perform slightly better.

Figure 5
Velocity–viscosity curves obtained for multiple motors using the mean-field model (A) and stochastic model (B) for the cargo of radius r = 0.5 μm for (a–b1)ω = 2 and (a–b2)ω = 0.5. Parameters for mean-field ...

We used the stochastic model to predict average run lengths and velocities in the presence of the viscous load (simulation details are discussed in the appendix). The results, obtained from the stochastic model in the presence of viscous drag (figure 5(B)), are qualitatively similar to those in the mean-field approximation. However, interestingly, multiple motors now move slower than a single motor even at low viscosity: when the viscous load is relatively small, the motors move almost with their unloaded velocity, and so sometime the forward motor advances too far and becomes loaded by the rearward motor. This causes two motors to move slower than a single motor under a small load.

3.4. Two to four motors are significantly correlated, while five and more motors become uncorrelated

We saw that stochastic fluctuations for the small number of motors are essential, as the results of the stochastic simulations deviate significantly from the conclusions of the mean-field theory. We observed that for a greater number of motors, the differences between the simulations and mean-field formulas decrease. To quantify this effect, we first generated time series of forces experienced by motors transporting the cargo against the constant total load. One such plot for two motors with applied total load of 4 pN at ω = 0.25 is shown in figure 6(A). This plot shows that the motors are highly mechanically anti-correlated. We then examined such time series for various motor numbers and nonlinearities by computing the average correlation (normalized covariance function) between the force time series (computation details are given in the appendix). We observed that the correlation depends only on the total number of motors N for small values of ω and is almost independent of applied load. For large values of ω, the correlation increases with the applied load. Importantly, however, for all values of ω, the correlations between the motors decrease with the motor number (see figure 6(B) for ω = 0.25). Figure 6(B) shows that two to four motors are strongly anti-correlated (negative values mean anti-correlation; values close to −1 mean strong anti-correlation, close to 0 mean weak correlation). The correlations weaken as the motor number increases and become insignificant if N > 4. Therefore, the transport by more than four motors can be modeled by using the mean-field approximation. Note that the average cross-motor-correlation coefficient decreases approximately as 1/N (figure 6(B)). The simple reason is that a random rapid displacement of any one motor generates the change of load almost equally shared by other (N − 1) motors, and so intuitively, the influence of force fluctuations in one motor on another scales as ~1/N.

Figure 6
(A) Time series of the forces applied to two motors transporting the cargo against the total load of 4 pN with ω = 0.25. Other parameter values are the same as in figure 3(D). (B) Average correlation between pairs of N motors collectively transporting ...

3.5. Scaling in the limit of the great motor number

We have calculated the force–processivity and force–velocity relations for tens of motors using the mean-field approximation. The results are shown in figure 7. The following asymptotic (large N) scaling behavior is clear from these plots. There is the analytical result for the average run length exponential dependence on N at zero load [24]:

Figure 7
(A) Force–processivity curves obtained from the mean-field model for ω = 2. The critical force (force at which the average run length is 8 nm) as a function of N is shown in the inset. (BD) Force–velocity curves obtained ...

This conclusion, of course, is valid for any ω. From figure 7(A), we observe that the multiple motor run lengths are exponentially decreasing with force:


for any ω (i.e. they are independent of the nature of single motor force–velocity curves). Here, coefficient α exponentially increases with N, while the simulations show that the coefficient β is not very sensitive to either ω, or N. The critical force (defined as the load at which the motors make on average but one step before the cargo detaches) increases linearly with N and is not sensitive to ω (figure 7(A)).

The force–velocity curves for the great motor number reveal interesting scaling (figures 7(B)–(D)). It can be understood from the following simple analysis: when the motors are almost uncorrelated, they indeed share the load almost equally. Then the force–velocity relation has the form V=v[1(FFsn)w], where the average number of working motors n has to be calculated from the nonlinear algebraic equation π(Nn)=εnexp(FFdn). Approximate asymptotic solution of this equation in the not very interesting case when Fs < Fd shows that almost all motors remain attached up to the stall, and the effective force–velocity curve is only slightly lower than the simple prediction V=v[1(FFsn)w], n = Nπ/(π + ε), up to the stall. In the case when Fs > Fd, numerical solution of these algebraic equations demonstrates at load force F less than cNFd, where c ≈ 0.7–0.3 for ε/π = 0.3–0.9; the effective force–velocity relation for many motors can also be approximated with the same formula. However, if F > cN Fd, then n → 0, almost all motors detach, and the velocity of the cargo plunges almost to zero. This simple semi-analytical scaling is easily seen in the plots of figures 7(B)–(D) obtained numerically from the complex formulas of the mean-field theory.

3.6. Significant increase with the load of the diffusivity of the cargo driven by multiple motors

Much useful information is contained in the statistical fluctuations about the mean velocity of the motors’ cargo. One quantity that can be monitored as the cargo progresses is the variance of the cargos displacement about its mean. It is easy to show that this variance grows linearly with time; the proportionality coefficient is equal to twice the effective diffusion coefficient, Deff [38], which in turn can be expressed by the formula Deff = Rdν/2, where R is the quantity called the randomness parameter in the case of a single motor [39]. This parameter is equal to 1 in the case of the simple hypothetic molecular motor, ‘Poisson stepper’ [39] that makes constant spatial steps at random times distributed exponentially. In certain sense, the more complex the motors mechanochemical cycle is, the higher its randomness parameter [39]. (Strictly speaking, the parameter R in the case of multiple motors, though a very useful quantity, cannot be called the ‘randomness parameter’.)

To the best of our knowledge, the parameter R was never estimated for the multiple motors. We recorded the rate of growth of variance of two and three coupled motors at various loads and calculated the respective parameters R. The results are shown in figure 8: as expected, the parameter R of one motor, which is the Poisson stepper in our case, is equal to 1. This is true also for multiple unloaded motors, because in this case the motors carry the cargo almost independently. However, the diffusivity increases dramatically with the load; the effect is an order of magnitude. The reason is the frequent rearward excursions when one of the motors detaches; such excursions increase and become more frequent at higher loads. This effect is more pronounced for the sub-linear motors, and is also less for three than for two motors, because the latter take on greater load each and detach more frequently springing back more. This prediction can be used in principle to determine the number of motors carrying the cargo.

Figure 8
R as a function of applied load for (A)ω = 0.25 (B)ω = 1 and (C)ω = 4. Parameter values are the same as in figure 3.

3.7. Stochastic effects in the tug-of-war between multiple opposing motors

The tug-of-war between multiple motors of opposing polarity was investigated in the framework of the mean field theory in [40]. In the symmetric case (four motors with the same free velocities and stall forces), the authors of [40] found that when the motors are ‘weak’ (small stall to detachment force ratio), then the cargo is almost always stalled. On the other hand, when the motors are ‘strong’ (great stall to detachment force ratio), either one type of motors ‘wins’ or another, and the cargo alternates between two directions of movement with the free motor speed. In the intermediate cases, the cargo is either stalled or moves in any direction with nearly free motor speed. We simulated the stochastic model with the same parameters as those used to produce figure 3 in [40]. We found that the nonlinearity parameter had little effect on the cargo’s speed histogram. However, the stochastic effects are very strong: figure 9(A) illustrates that though there is indeed a significant probability of the stalled cargo, the non-zero velocities of the cargo are greatly dispersed: all non-zero velocities from 0 to the plus/minus free motor speed are almost equally probable, and the significant non-zero velocity peaks predicted by the mean field theory are almost completely smeared out by the stochastic effect (our figure 9(A) is to be compared with figure 3(A3,B3,C3) in [40]. Note that our figure 9(A) shows data for a different number of motors).

Figure 9
(A) Velocity distributions obtained from the stochastic model for the symmetric tug-of-war of N+ = N = 2 plus and minus motors with ω = 1. Motility behavior shown in (a2) was obtained using single motor parameters as in figure 3 that ...

Similarly, in the asymmetric case, when one motor type is stronger than the other, and the numbers of the opposing motors are equal, multiple peaks in the velocity distribution are predicted by the mean-field theory [40]. In contrast, we found that, again, these peaks are obliterated by the stochastic movements, and the velocity is widely distributed around these peaks (our figure 9(B) is to be compared with figure 4(A3,B3,C3) in [40]. Note that our figure 9(B) shows data for a different number of motors).

4. Discussion

In this paper, we investigated the nonlinear and stochastic effects on the collective motor transport. We show that in the mean-field approximation, the super-linear force–velocity relation improves the multiple motors’ performance, while the sub-linear relation makes it worse. However, when the stochastic load sharing is taken into account, we see that the ‘leading’ motor takes on a disproportionately great load and detaches frequently causing retractions of the cargo. This effect significantly worsens the super-linear motors’ performance at low loads. Thus, we make the following predictions to be tested in future experiments: due to the combination of the nonlinear and stochastic effects, the collective force–velocity curve for two or three motors becomes almost linear, and the rate of movement against moderate loads is not, in fact, accelerated by increasing the small number of motors. The run length, on the other hand, increases exponentially with the motor number. An additional, potentially useful, future application of our model is using a comparison between the predicted relations and data to infer the number of motors: this technique was used in [35] by applying the formula of the mean-field theory, which is not quantitatively accurate for the small motor number.

We observe that the collective motor behavior is insensitive to the compliance of the motor stalks in the physiological range. Furthermore, we determine that the velocity of the small number of coupled motors is insensitive to viscous loads if the effective viscosity of the cytoplasm is less than ~1 Pa×s. At higher viscosity, the velocity decreases linearly with the viscosity; respective relation is not sensitive to the exact form of the force–velocity curve and the motor number. We also argue that details of the exact mechanochemical cycle are not crucial for the average collective motor performance (see the appendix). In summary, our simulations predict a very robust collective motor behavior.

The model predicts that two to four coupled motors are significantly anti-correlated, so the stochastic effects for the small number of motors are significant. We find that five or more coupled motors correlate weakly, share the load almost evenly, and the mean-field approximation describes the collective motor behavior accurately. For a great number of motors, the mean-field theory makes two predictions. First, the run length increases exponentially with the motor number (this conclusion was first reached in [24]) and decreases exponentially with the load. Multiple motor force–processivity relations are insensitive to the nonlinearities of the individual motor force–velocity relation. Second, there are two regimes in the multiple motor force–velocity relations if Fd < Fs. At a significant load, proportional to the total number of motors F > cN Fd, c ~ 1, the average cargo velocity decreases rapidly. At smaller loads, the collective velocity scales in a very simple way with the motor number.

We calculated the effective diffusivity and randomness parameter characterizing the rate of growth of the displacement variance of the cargo with time and found that the diffusivity increases drastically, by the order of magnitude for two and three motors due to significant rearward excursions of the cargo upon detachment of one of the motors. For large loads, the diffusivity depends on the motor number, so in principle, this prediction can be used to infer the number of motors carrying the cargo if the variance of the displacement and the force–velocity relation are measured. Furthermore, we investigated the tug-of-war between the opposing multiple motors, and found that the nonlinearity does not affect the results, but the stochastic effects cause wide dispersion of the cargo velocity and smearing out of the multiple peaks predicted by the mean field theory. One interesting practical application of this prediction is that qualitative information about the number of opposing motors can be inferred in principle from the measurements of the cargos velocity distribution (such distributions for competing opposite bipolar kinesin and ncd motors were measured, for example, in [41]).

One of the interesting conclusions we can make from the model results is that if the ‘objective’ of the cell is to keep the velocity of the cargo unaffected by a significant load, then more than three super-linear motors have to transport the cargo. Most of the quantitative model predictions will have to wait to be tested in future experiments. Only a very preliminary attempt to measure effective force–velocity curves for two and three motors is made [42]. Just very general rules of the collective motor transport are tested experimentally by now: additive stall forces and drastically increased with the motor number travel distances are established in vitro for Kinesin-1 [9] and Dynein [13]. Note that recently the unexpected experimental finding was reported [43]: while two motors produce longer average run lengths than single kinesins, the system effectively behaves as though a single-motor attachment state dominates motility. The authors of this study proposed that negative motor interference derived from asynchronous motor stepping can explain this effect. Our model does not reproduce this observation, indicating perhaps the presence of more complex elastic properties of the motor link and mechanochemical coupling of the motor cycle with the link mechanics.

In the future, the most immediate problems that the modeling of the collective motor transport will have to address are as follows. First, by comparing the increasingly available data with computational screening of the model parameter space, we need to obtain the dissociation rate dependence on the load for individual motors, as well as to understand better the motor behavior at super-stall forces and at forces pulling the motor forward. Second, there are very puzzling differences between collective motor in vitro and in vivo behaviors: in vivo, a reduced motor number causes slight increases of the travel velocities, while the travel distances are not reduced [44], which is in stark contrast to the simple theoretical predictions and observed in vitro behavior. We believe that insight gained from modeling studies will be crucial for correct interpretation of future experimental results.


This study was funded by National Institutes of Health grant GM068952 to AM.


A.1. Simulation of the single motor

We use the Monte Carlo procedure [32, 38] to update the state (position and engaged or detached state) in increments of the time step Δt. The time step Δt is chosen to be sufficiently smaller than the fastest characteristic time (in our case, detachment of the last attached motor under a high load). We used Δt = 10−5 s that conforms with this requirement. The computational procedure is as follows.

  1. Initial condition: at t = 0, x = 0, where x is the position of the motor on the track.
  2. Updating procedure: repeat the following steps up to tmax in increments of Δt.
    1. If t > tmax go to step 3.
    2. Detachment: calculate Poff = ε(F)* Δt. First try detachment with the probability Poff. If the detachment occurs, go to step 3, else go to step (iii).
    3. Stepping: if the motor remains attached after step (ii), stepping occurs with probability Pstep = kstep(F)* Δt. After stepping, x is changed to x + d where d = 8 nm.
  3. Run length is the current value of x. Velocity is obtained by dividing the current position x by the current time t.

A.2. Simulation of the multiple motors

We put N motors on cargo, so that the motor heads are attached to the cargo via the links of the rest length l = 0.11 μm [11]. Each link exerts a restoring force when stretched beyond their rest length. The links have no compressional rigidity, i.e. they exert no force when compressed. Initially, we place the bead’s center of mass at the origin and allow all motors to attach to any discrete binding site on the track within distance l on either side of the bead. Once the motors are attached, we calculate the initial position of the bead’s center of mass so that the sum of all elastic forces applied to the bead is equal to zero; in what follows, the bead’s position is calculated at each step so that the sum of all elastic forces from the motor links has to be equal to zero.

For each time step, we visit each of the N motors and determine their tentative states (attached or detached) and positions. During the updating procedure, at each computational step, each motor’s state is updated once. If the motor is currently unattached, we allow it to attach with a probability Pon = π * Δt, determined by the ‘on-rate’ π, to any binding site on the track within distance l on either side from the bead’s center of mass.

If the motor is currently attached, a load Fi felt by the ith motor is obtained by multiplying the extension of its link Δli by the link’s stiffness k, and there are three possibilities: the motor can remain stationary, advance, or detach. Probabilities of these three events are determined from the single motor model based on the current load on the motor: (i) Poff is calculated using equation (3) irrespective of the direction of the force applied to the motor; (ii) Pstep is calculated using equation (2) for backward loads FiFs; for backward load greater than Fs, Pstep = 0; a forward load does not alter the motor cycle, so we substitute Fi = 0 for forward loads in equation (2). If the motor steps, its position xi is changed to xi + d. When we determine the tentative states and positions of all N motors, we update the states and positions of all motors simultaneously. Then, the number of engaged motors n and their locations are recorded and the bead position is updated.

Viscous load: in the presence of the viscous load, the position of the bead is determined not by the balancing of the elastic motor forces to zero, but by the viscous force that bead experiences: if the bead is subjected to the net force f, this causes it to move with velocity νdrift = f/γ. The net motion of the bead over the time interval 3t is given by the deterministic drift xdrift = νdrift × Δt. The net force f on the cargo is given by the formula f=i=1Nfi, where fi is the elastic restoring force exerted by the ith motor on the cargo, which magnitude depends on the extension of the ith link.

Solution of equation (13): the nonlinear algebraic equation (13) was solved with Mathematica using the function Solve.

Calculation of correlations: correlations between the motors’ time series were calculated using the Maltab time-series tool and function corrcoef.

A.3. Model with mechano-chemical cycle

In order to test if the details of a mechano-chemical cycle affect the average mechanical properties of the collective motor transport, we analyzed the following published kinesin model [30] in the multiple motor case. The motor cycle is based on the ATP binding to the motor (M) and subsequent hydrolysis:


Here kon(koff) is the rate constant for binding (unbinding) of ATP, and kcat is the rate constant for ATP hydrolysis. The mechanical motor step takes place synchronously with the hydrolysis and release of its products. Stationary solution of the standard Markov chain equations based on this cycle gives the Michaelis–Menten expression for the deterministic motor velocity:


Here d is the step size, F is the load and [ATP] is the ATP concentration. kcat(F) and koff are the load-dependent hydrolysis and ATP-unbinding rates, respectively. By choosing these load dependences in the form


we get a one-motor force–velocity curve that is very similar to the one given by equation (1). In these expressions, kB T is the thermal energy, and dl is the characteristic molecular scale. At saturating ATP concentrations, motor velocity predicted by this model


becomes identical to that given by equation (1), if dkcat(0) = ν is the unloaded velocity. To model the load-dependent unbinding of the motor, we assume that it detaches with rate given by equation (3) from any state.

We ran Monte Carlo simulations of this mechano-chemical model simply by adding the transitions between the chemical steps of each of the motors attached to the bead to the mechanical steps described previously (mechanical steps took place synchronously with the hydrolysis steps). The probability of each of these transitions was computed by multiplying the respective instantaneous load-dependent rate by the time step duration. The parameter values used for the simulations were the same as described above; in addition, kon = 2 × 106 M−1 × s−1, koff(0) = 55 s−1, dl = 1.6 nm, T = 300 K, kcat(0) = 125 s−1 and [ATP] = 5 mM.

Results obtained from the multiple motor simulations using this mechano-chemical model were not significantly different from those obtained from the purely mechanical model. There was no difference in the average run-lengths and velocities under low loads. For significant loads, maximal differences in the average run-lengths were 0.004 μm for two motors, and 0.006 μm for three motors, respectively. Similarly, maximal differences in the average velocities under significant loads were 0.02 μm s−1 for two motors, and 0.04 μm s−1 for three motors, respectively.


1. Vale RD. The molecular motor toolbox for intracellular transport. Cell. 2003;112:467–70. [PubMed]
2. Rogers SL, Gelfand VI. Membrane trafficking, organelle transport, and the cytoskeleton. Curr Opin Cell Biol. 2000;12:57–62. [PubMed]
3. Welte MA, Gross SP, Postner M, Block SM, Wieschaus EF. Developmental regulation of vesicle transport in Drosophila embryos: forces and kinetics. Cell. 1998;92:547. [PubMed]
4. Walczak CE, Heald R. Mechanisms of mitotic spindle assembly and function. Int Rev Cytol. 2008;265:111–8. [PubMed]
5. Glotzer M. The molecular requirements for cytokinesis. Science. 2005;307:1735–9. [PubMed]
6. Schliwa M, Woehlke G. Molecular motors. Nature. 2003;422:759–65. [PubMed]
7. Moffitt JR, Chemla YR, Smith SB, Bustamante C. Recent advances in optical tweezers. Annu Rev Biochem. 2008;77:205–8. [PubMed]
8. Kolomeisky AB, Fisher ME. Molecular motors: a theorist’s perspective. Annu Rev Phys Chem. 2007;58:675. [PubMed]
9. Vershinin M, Carter BC, Razafsky DS, King SJ, Gross SP. Multiple-motor based transport and its regulation by Tau. Proc Natl Acad Sci USA. 2007;104:87–92. [PubMed]
10. Schnitzer MJ, Visscher K, Block SM. Force production by single kinesin motors. Nat Cell Biol. 2000;2:718–23. [PubMed]
11. Coppin CM, Pierce DW, Hsu L, Vale RD. The load dependence of kinesin’s mechanical cycle. Proc Natl Acad Sci USA. 1997;94:8539–44. [PubMed]
12. Gross SP, Vershinin M, Shubeita GT. Cargo transport: two motors are sometimes better than one. Curr Biol. 2007;17:R478–86. [PubMed]
13. Mallik R, Petrov D, Lex SA, King SJ, Gross SP. Building complexity: an in vitro study of cytoplasmic dynein with in vivo implications. Curr Biol. 2005;15:2075. [PubMed]
14. Badoual M, Julicher F, Prost J. Bidirectional cooperative motion of molecular motors. Proc Natl Acad Sci USA. 2002;99:6696–701. [PubMed]
15. Duke T. Cooperativity of myosin molecules through strain-dependent chemistry. Philos Trans R Soc Lond B Biol Sci. 2000;355:529–38. [PMC free article] [PubMed]
16. Shu Y, Shi H. Cooperative effects on the kinetics of ATP hydrolysis in collective molecular motors. Phys Rev E. 2004;69(2 Pt 1):021912. [PubMed]
17. Leibler S, Huse DA. Porters versus rowers: a unified stochastic model of motor proteins. J Cell Biol. 1993;121:1357–8. [PMC free article] [PubMed]
18. Julicher F, Prost J. Cooperative molecular motors. Phys Rev Lett. 1995;75:2618–21. [PubMed]
19. Vermeulen KC, Stienen GJM, Schmid CF. Cooperative behavior of molecular motors. J Muscle Res Cell Motil. 2002;23:71–9. [PubMed]
20. Vilfan A, Frey E, Schwabl F. Elastically coupled molecular motors. Eur Phys J. 1998;B 3:535–46.
21. Julicher F, Ajdari A, Prost J. Modeling molecular motors. Rev Mod Phys. 1997;69:1269–81.
22. Campas O, Kafri Y, Zeldovich KB, Casademunt J, Joanny J-F. Collective dynamics of interacting molecular motors. Phys Rev Lett. 2006;97:038101. [PubMed]
23. Campas O, Leduc C, Bassereau P, Casademunt J, Joanny J-F, Prost J. Coordination of Kinesin motors pulling on fluid membranes. Biophys J. 2008;94:5009–17. [PubMed]
24. Klumpp S, Lipowsky R. Cooperative cargo transport by several molecular motors. Proc Natl Acad Sci USA. 2005;102:17284–9. [PubMed]
25. Nedelec F. Computer simulations reveal motor properties generating stable antiparallel microtubule interactions. J Cell Biol. 2002;158:1005. [PMC free article] [PubMed]
26. Wollman R, Civelekoglu-Scholey G, Scholey JM, Mogilner A. Reverse engineering of force integration during mitosis in the Drosophila embryo. Mol Syst Biol. 2008;4:195. [PMC free article] [PubMed]
27. Valentine MT, Fordyce PM, Krzysiak TC, Gilbert SP, Block SM. Individual dimers of the mitotic kinesin motor Eg5 step processively and support substantial loads in vitro. Nat Cell Biol. 2006;8:470–6. [PMC free article] [PubMed]
28. Gebhardt JCM, Clemen AE-M, Jaud J, Rief M. Myosin-V is a mechanical ratchet. Proc Natl Acad Sci USA. 2006;103:8680–5. [PubMed]
29. Janson ME, Dogterom M. Scaling of microtubule force–velocity curves obtained at different tubulin concentrations. Phys Rev Lett. 2004;92:248101. [PubMed]
30. Singh MP, Mallik R, Gross SP, Yu CC. Monte Carlo modeling of single-molecule cytoplasmic dynein. Proc Natl Acad Sci USA. 2005;102:12059–64. [PubMed]
31. Toba S, Watanabe TM, Yamaguchi-Okimoto L, Toyoshima YY, Higuchi H. Overlapping hand-over-hand mechanism of single molecular motility of cytoplasmic dynein. Proc Natl Acad Sci USA. 2006;103:5741–5. [PubMed]
32. Kunwar A, Vershinin M, Xu J, Gross SP. Stepping, strain gating, and an unexpected force–velocity curve for multiple-motor-based transport. Curr Biol. 2008;18:1173. [PMC free article] [PubMed]
33. Keren K, Pincus Z, Allen GM, Barnhart EL, Marriott G, Mogilner A, Theriot JA. Mechanism of shape determination in motile cells. Nature. 2008;453:475–80. [PMC free article] [PubMed]
34. Grill SW, Kruse K, Julicher F. Theory of mitotic spindle oscillations. Phys Rev Lett. 2005;94:108104. [PubMed]
35. Beeg J, Klumpp S, Dimova R, Gracia RS, Unger E, Lipowsky R. Transport of beads by several kinesin motors. Biophys J. 2008;94:532–41. [PMC free article] [PubMed]
36. Jeney S, Stelzer EHK, Grubmüller H, Florin E-L. Mechanical properties of single motor molecules studied by three-dimensional thermal force probing in optical tweezers. Chemphyschem. 2004;5:1150–8. [PubMed]
37. Kural C, Kim H, Syed S, Goshima G, Gelfand VI, Selvin PR. Kinesin and dynein move a peroxisome in vivo: a tug-of-war or coordinated movement? Science. 2005;308:1469–72. [PubMed]
38. Mogilner A, Elston T, Wang H-Y, Oster G. Molecular motors: theory. In: Fall CP, Marland E, Tyson J, Wagner J, editors. Joel Keizer’s Computational Cell Biology. New York: Springer; 2002. pp. 321–55.
39. Schnitzer MJ, Block SM. Statistical kinetics of processive enzymes. Cold Spring Harb Symp Quant Biol. 1995;60:793–802. [PubMed]
40. Müller MJI, Klumpp S, Lipowsky R. Tug-of-war as a cooperative mechanism for bidirectional cargo transport by molecular motors. Proc Natl Acad Sci USA. 2008;105:4609–14. [PubMed]
41. Tao L, Mogilner A, Civelekoglu-Scholey G, Wollman R, Evans J, Stahlberg H, Scholey JM. A homotetrameric kinesin-5, KLP61F, bundles microtubules and antagonizes Ncd in motility assays. Curr Biol. 2006;16:2293–302. [PubMed]
42. Shtridelman Y, Cahyuti T, Townsend B, DeWitt D, Macosko JC. Force–velocity curves of motor proteins cooperating in vivo. Cell Biochem Biophys. 2008;52:19–29. [PMC free article] [PubMed]
43. Rogers AR, Driver JW, Constantinou PE, Kenneth Jamison D, Diehl MR. Negative interference dominates collective transport of kinesin motors in the absence of load. Phys Chem Chem Phys. 2009;11:4882–9. [PubMed]
44. Shubeita GT, Tran SL, Xu J, Vershinin M, Cermelli S, Cotton SL, Welte MA, Gross SP. Consequences of motor copy number on the intracellular transport of kinesin-1-driven lipid droplets. Cell. 2008;135:1098–107. [PMC free article] [PubMed]