Stochastic simulations of cargo transport by processive molecular motors
Christian B. Korn, Stefan Klumpp, Reinhard Lipowsky, Ulrich S. Schwarz
aa r X i v : . [ q - b i o . S C ] D ec Stochastic simulations of cargo transportby processive molecular motors
Christian B. Korn, Stefan Klumpp,
2, 3
Reinhard Lipowsky, and Ulrich S. Schwarz
1, 4 University of Heidelberg, Bioquant 0013,Im Neuenheimer Feld 267, D-69120 Heidelberg, Germany Center for Theoretical Biological Physics,University of California at San Diego,9500 Gilman Drive, La Jolla, CA 92093-0374, USA Max Planck Institute for Colloids and Interfaces,Science Park Golm, D-14424 Potsdam, Germany University of Heidelberg, Institute for Theoretical Physics,Philosophenweg 19, D-69120 Heidelberg, Germany
Abstract
We use stochastic computer simulations to study the transport of a spherical cargo particle alonga microtubule-like track on a planar substrate by several kinesin-like processive motors. Our newlydeveloped adhesive motor dynamics algorithm combines the numerical integration of a Langevinequation for the motion of a sphere with kinetic rules for the molecular motors. The Langevinpart includes diffusive motion, the action of the pulling motors, and hydrodynamic interactionsbetween sphere and wall. The kinetic rules for the motors include binding to and unbindingfrom the filament as well as active motor steps. We find that the simulated mean transportlength increases exponentially with the number of bound motors, in good agreement with earlierresults. The number of motors in binding range to the motor track fluctuates in time with aPoissonian distribution, both for springs and cables being used as models for the linker mechanics.Cooperativity in the sense of equal load sharing only occurs for high values for viscosity andattachment time.
PACS numbers: 82.39.-k,87.10.Mn,87.15.hj . INTRODUCTION Molecular motors play a key role for the generation of movement and force in cellularsystems [1]. In general there are two fundamentally different classes of molecular motors.Non-processive motors like myosin II motors in skeletal muscle bind to their tracks onlyfor relatively short times. In order to generate movement and force, they therefore have tooperate in sufficiently large numbers. Processive motors like kinesin remain attached to theirtracks for a relatively long time and therefore are able to transport cargo over reasonabledistances. Indeed processive cytoskeletal motors predominantly act as transport engines forcargo particles, including vesicles, small organelles, nuclei or viruses. For example, kinesin-1motors make an average of 100 steps of size 8 nm along a microtubule before detaching fromthe microtubule [2, 3], and therefore reach typical run lengths of micrometers.However, for intracellular transport even processive motors tend to function in ensemblesof several motors, with typical motor numbers in the range of 1-10 [4]. The cooperationof several motors is required, for example, when processes like extrusion of lipid tethersrequire a certain level of force that exceeds the force generated by a single motor [5, 6]. Thecooperative action of several processive motors is also required to achieve sufficiently longrun length for cargo transport [7], as transport distances within cells are typically of theorder of the cell size, larger than the micron single motor run length. In this context the mostprominent example is axonal transport, as axons can extend over many centimeters [8, 9].Another level of complexity of transport within cells is obtained by the simultaneous presenceof different motor species on the same cargo, which can lead to bidirectional movements andswitching between different types of tracks [10, 11], and by exchange of components of themotor complex with the cytoplasm.Cargo transport by molecular motors can be reconstituted in vitro using so-called beadassays in which motor molecules are firmly attached to spherical beads that flow in aqueoussolution in a chamber. On the bottom wall of the chamber microtubules are mounted alongwhich the beads can be transported [1, 12, 13, 14]. This assay has been used extensively tostudy transport by a single motor over the last decade [12], but recently several groups haveadapted it for the quantitative characterization of transport by several motors [14, 15]. Ifseveral motors on the cargo can bind to the microtubule, then the transport process continuesuntil all motors simultaneously unbind from the microtubule. Based on a theoretical model2or cooperative transport by several processive motors, it was recently predicted that themean transport distance increases essentially exponentially with the number of availablemotors [7]. Indeed these predictions are in good agreement with experimental data [14, 16].However, both the theoretical approach and the experiments do not allow us to investigatethe details of this transport process. A major limitation of the bead assays for transport byseveral motors is that the number of motors per bead varies from bead to bead and that onlythe average number of motors per bead is known [14, 15]. In addition, even if the numberof motors on the bead was known, the number of motors in binding range would still be afluctuating quantity. Recently two kinesin motors have been elastically coupled by a DNAscaffold and the resulting transport has been analyzed in quantitative detail [16]. However,it is experimentally very challenging to extend this approach to higher numbers of motors.One key property of transport by molecular motors is the load force dependence of thetransport velocity. For transport by single kinesins, the velocity decreases approximatelylinearly with increasing load and stalls at a load of about 6 pN [3]. Thus, when the cargohas to be transported against a large force, the speed of a single motor is slowed down.However, if several motors simultaneously pull the cargo, they could share the total load.This cooperativity lets them pull the cargo faster. Assuming equal load sharing, one can showthat in the limit of large viscous load force the cargo velocity is expected to be proportionalto the number of pulling motors [7]. Indeed, this is one explanation that was proposedby Hill et al. to give plausibility to their results from in vivo experiments, which showedthat motor-pulled vesicles move at speeds of integer multiples of a certain velocity [9]. Ingeneral, however, one expects that the total load is not equally shared by the set of pullingmotors. The force experienced by each motor will depend on its relative position along thetrack and can be expected to fluctuate due to the stochasticity of the motor steps [17]. Inaddition, for a spherical cargo particle curvature effects are expected to play a role of theway force is transmitted to the different motors. Because of its small size, the cargo particleis perpetually subject to thermal fluctuations. This diffusive particle motion is also expectedto affect the load distribution and depends on the exact height of the sphere above the walldue to the hydrodynamic interactions.In order to investigate these effects, here we introduce an algorithm that allows us tosimulate the transport of a spherical particle by kinesin-like motors along a straight filamentthat is mounted to a plain wall. Binding and unbinding of the motor to the filament3an be described in the same theoretical framework as the reaction dynamics of receptor-ligand bonds [18, 19, 20]. Similarly as receptors bind very specifically to certain ligands,conventional kinesin binds only to certain sites on the microtubule. Thus from the theoreticalpoint of view a spherical particle covered with motors binding to tracks on the substrate isequivalent to a receptor-covered cell binding to a ligand-covered substrate. This situation isreminescent of rolling adhesion, the phenomenon that in the vasculature different cell types(mainly white blood cells, but also cancer cells, stem cells or malaria-infected red bloodcells) bind to the vessel walls under transport conditions [21]. Different approaches havebeen developed to understand the combination of transport and receptor-ligand kinetics inrolling adhesion. Among these Hammer and coworkers developed an algorithm that combineshydrodynamic interactions with reaction kinetics for receptor-ligand bonds [22]. Recently,we introduced a new version of this algorithm that also includes diffusive motion of thespherical particle [23]. Here, we further extend our algorithm to include the active steppingof motors ( adhesive motor dynamics ). Simulation experiments with this algorithm provideaccess to experimentally hidden observables like the number of actually pulling motors,the relative position of the motors to each other, and load distributions. The influence ofthermal fluctuations on the motion of the cargo particle is also influenced by the propertiesof the molecules that link the cargo to the microtubule. In our simulations various polymermodels can be implemented and their influence can be tested directly. In general our methodmakes it possible to probe the effects of various microscopic models for motor mechanics onmacroscopic observables that are directly accessible to experiments.The organization of the article is as follows. In the first part, Sec. II, we explain our modelin detail. This is based on a Langevin equation that allows us to calculate the position andorientation of a spherical cargo particle as a function of time. In addition, we include rulesthat model the reaction kinetics of the molecular motors being attached to the cargo andcomment on the different kinds of friction involved. We then explain how theoretical resultsfor the dependence of the mean run length of a cargo particle on the number of availablemotors previously obtained in the framework of a master equations can be compared to thesituation where only the total number of motors attached to the cargo sphere is known.We also briefly comment on the implementation of our simulations. In the results part,Sec. III, we first measure the mean run length and the mean number of pulling motors at lowviscous drag and find good agreement with earlier results. We then present measurements4f quantities that are not accessible in earlier approaches, including the dynamics of thenumber of motors on the cargo that are in binding range to the microtubule. Finally,we consider cargo transport in the high viscosity regime and investigate how the load isdistributed among the pulling motors. We find that cooperativity by load sharing stronglydepends on appropriate life times of bound motors. In the closing part, Sec. IV, we discussto what extend our simulations connect theoretical modeling with experimental findings.Furthermore, we give an outlook on further possible applications of the adhesive motordynamics algorithm introduced here.
II. MODEL AND COMPUTATIONAL METHODSA. Bead dynamics
In experiments using bead assays for studying the collective transport behavior of kinesinmotors one notes the presence of three very different length scales, namely the chamberdimension, the bead size and the molecular dimensions of kinesin and microtubule, respec-tively. The chamber dimension is macroscopic. The typical radius R of beads in assaychambers is in the micrometer range [13]. The kinesin molecules with which they are cov-ered have a resting length l of about 80 nm [24]. Kinesins walk along microtubules which arelong hollow cylindrical filaments made from 13 parallel protofilaments and with a diameterof about h MT = 24 nm [25]. Thus the chamber dimensions are large compared to the beadradius, which in turn is large compared to the motors and their tracks. This separation oflength scales allows us to model the microtubule as a line of binding sites covering the walland means that the dominant hydrodynamic interaction is the one between the sphericalcargo and the wall. For sufficiently small motor density, this separation of length scales alsoimplies that we have to consider only one lane of binding sites. In the following we thereforeconsider a rigid sphere of radius R moving above a planar wall with an embedded line ofbinding sites as a simple model system for the cargo transport by molecular motors along afilament.For small objects like microspheres typical values of the Reynolds number are muchsmaller than one and inertia can be neglected (overdamped regime). Therefore, the hydro-dynamic interaction between the sphere and the wall is described by the Stokes equation.5hroughout this paper we consider vanishing external flow around the bead. Directionalmotion of the sphere arises from the pulling forces exerted by the motors. In addition thebead is subject to thermal fluctuations that are ubiquitous for microscopic objects. Tra-jectories of the bead are therefore described by an appropriate Langevin equation. For thesake of a concise notation we introduce the six-dimensional state vector X which includesboth the three translational and the three rotational degrees of freedom of the sphere. Thetranslational degrees of freedom of X refer to the center of mass of the sphere with respectto some reference frame (cf. Fig. 1). The rotational part of X denotes the angles by whichthe coordinate system fixed to the sphere is rotated relative to the reference frame [26].Similarly, F denotes a combined six-dimensional force/torque vector.With this notation at hand the appropriate Langevin equation reads [27, 28]˙ X = M F + k B T ∇ M + g It . (1)Here, M is the position-dependent 6 × M depends on the height of the sphere above the wall in such a waythat the mobility is zero when the sphere touches the wall [29]. Thus, the hydrodynamicinteraction between the sphere and the wall is completely included in the configurationdependence of the mobility matrix. The last term in Eq. (1) is a Gaussian white noise termwith h g It i = 0 , h g It g It ′ i = 2 k B T M δ ( t − t ′ ) , (2)with Boltzmann’s constant k B . The second equation represents the fluctuation-dissipationtheorem illustrating that the noise is multiplicative due to the position-dependance of M .Thus we also have to define in which sense the noise in Eq. (1) shall be interpreted. Asusual for physical processes modeled in the limit of vanishing correlation time we choosethe Stratonovich interpretation [30]. However, Eq. (1) is written in the Itˆo version markedby the super-index I for the noise term. The Itˆo version provides a suitable base for thenumerical integration of the Langevin equation using a simple Euler scheme. The gradientterm in Eq. (1) is the combined result of using the Itˆo version of the noise and a term thatcompensates a spurious drift term arising from the no-slip boundary conditions [26, 27].6 (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) PSfrag replacements v F s v m F t h M T F m F m,x v m v m l | F m,x | U xyz χ h R R ˆ r FIG. 1: A single sphere of radius R and surface separation h from a planar wall. The translationalcoordinates R of the sphere are given relative to a reference frame that is fixed to the wall. Thesphere is pulled by one molecular motor that is attached to the surface at position ˆ r measured withrespect to the center of the sphere. The bead is subject to the motor force F m with x -component F m,x . In addition, an external force F t acts parallel to the filament, typically arising from anoptical trap. The force unbalance between F m,x and F t leads to the bead velocity U . The motorwith resting length l is firmly attached to the bead and can bind to and unbind from a MT andmoves with velocity v m . χ denotes the angle between motor and MT. For our numerical simulations we discretize Eq. (1) in time and use an Euler algorithmwhich is of first order in the time step ∆ t ∆ X t = F ∆ t + k B T ∇ M ∆ t + g (∆ t ) + O (∆ t ) (3)where g (∆ t ) has the same statistical properties as g It from above. In order to compute theposition-dependent mobility matrix M we use a scheme presented by Jones et al. [29, 31]that provides accurate results for all values of the height z . A detailed description of thecomplete algorithm including the translational and rotational update of the sphere can befound in Ref. [26]. B. Motor dynamics
In our model the spherical cargo particle is uniformly covered with N tot motor proteins.These molecules are firmly attached to the sphere at their foot domains. Consequently, N tot is a constant in time. The opposite ends of these molecules (their head domains) canreversibly attach to the microtubule (MT), which is modeled as a line of equally spaced7inding sites for the motors covering the wall. A motor that is bound to the MT exerts aforce and a torque to the cargo particle. If r h and r f are the positions of the head and footdomains of the motor, respectively, then the force by the motor F m is given by: F m = ˆ r m F m , F m = F ( r m ) , ˆ r m = r h − r f k r h − r f k , r m = k r h − r f k , (4)where the absolute value F m of the motor force is given by the force extension relation F ( x ) for the motor protein. Throughout this article we consider two variants of the forceextension curve for the polymeric tail of the motor. The harmonic spring model reads F ( x ) = κ ( x − l ) . (5)This means, a force is needed for both compression and extension of the motor protein.Actually, it was found that kinesin exhibits a non-linear force extension relation [32], withthe spring constant varying between κ = 0 . · − N/m for small extensions and κ =0 . · − N/m for larger extensions [32, 33]. For extensions close to the contour lengththe molecule becomes infinitely stiff [34]. For small extensions, however, the harmonicapproximation works well. Alternatively, we consider the cable model F ( x ) = κ ( x − l )Θ( x − l ) , Θ( x ) = , x > , else . (6)In the cable model force is only built up when the motor is extended above its resting length l . In the compression mode, i. e. when the actual motor length is less then the restinglength, no force exists. The cable model can be seen as the simplest model for a flexiblepolymer.Besides the force each motor attached to the MT also exerts a torque on the cargo particle.With ˆ r being the position of the motor foot relative to the center of the sphere (cf. Fig. 1),this torque reads T m = ˆ r × F m . (7)The combined force/torque vector F = ( F m , T m ) T enters the Langevin algorithm Eq. (3).8n addition to the firm attachment of the motors to the sphere, each motor can in principlereversibly bind and unbind to the MT. We model these processes as simple Poissonian rateprocesses in similar a fashion as it is done for modeling of formation and rupture of receptorligand bonds in cell adhesion (e. g., [22]). The motor head is allowed to rotate freely about itspoint of fixation on the cargo. The head is therefore located on a spherical shell with radiusgiven by the motor resting length l . However, in contrast to the anchorage point on the cargoparticle the head position of the motor is not explicitely resolved by the algorithm. In orderto model the binding process we introduce a capture length r . A motor head can then bindto the MT with binding rate π ad whenever the spherical shell of radius l and thickness r around the motor’s anchorage point has some overlap with a non-occupied binding site on theMT. The MT’s binding sites are identified by the tubulin building block of length δ = 8 nm,so we choose r = δ/
2. Note that binding rate and binding range are complementaryquantities and that a more detailed modeling of the binding process would have to yieldappropriate values for both quantities. If the overlap criterion is fulfilled within a timeinterval ∆ t , the probability for binding p on within this time step is p on = 1 − exp( − π ad ∆ t ).With a standard Monte-Carlo technique it is then decided whether binding occurs or not:a random number rand is drawn from the uniform distribution on the interval (0 ,
1) and inthe case p on > rand binding occurs. If the overlap criterion is fulfilled for several bindingsites, then using the Monte-Carlo technique it is first decided whether binding occurs andthen one of the possible binding sites is randomly chosen.Each motor bound to the MT can unbind with escape rate ǫ . In single motor experimentsit was found that the escape rate increases with increasing motor force [35]. This forcedependence can be described by the Bell equation [7, 38] ǫ = ǫ exp( F m /F d ) , (8)with the unstressed escape rate ǫ and the detachment force F d . Because the details of forcedmotor unbinding from a filament are not known, here we make the simple assumption thatthe unbinding pathway is oriented in the direction of the tether. We set ǫ = ǫ wheneverthe motor is compressed.The major conceptual difference between a motor connecting a sphere with a MT anda receptor-ligand bond is that a motor can actively step forward from one binding site to9arameter typical value meaning reference R µ m bead radius ǫ − unstressed escape rate [7] π ad s − binding rate [6, 7] F d κ − . . . − . . . − N/m motor spring constant [33, 34, 36] δ v µ m/s maximum motor velocity [7] λ s := v /δ
125 s − forward step rate r δ/ F s . . . . . . l h MT
24 nm microtubule diameter [25]
TABLE I: Parameters used for adhesive motor dynamics. For ambient temperature we use T =293 K, for viscosity η = 1 mPa s (if not otherwise stated). If a range is given, then figure in boldface denotes the value used in the numerical simulations. the next with step length δ (the length of the tubulin unit). The mean velocity v of anunloaded kinesin motor is about v = 1 µ m/s depending on ATP concentration [3]. If themotor protein is mechanically loaded with force opposing the walking direction, the motorvelocity v m is decreased. For a single kinesin molecule that is attached to a bead on which atrap force F t pulls, the velocity was found experimentally to decrease approximately linearly[3, 39]: v m = v (cid:18) − F t F s (cid:19) , < F t < F s , (9)with the stall force F s and the trap force F t acting antiparallel to the walking direction.Different experiments have reported stall forces between 5 and 8 pN. Changes in this rangeare not essential for our results and therefore we use the intermediate value F b = 6 pN .If the force is higher than the stall force, kinesin motors walk backwards with a very lowvelocity [40], which we will neglect in the following. Finally, for assisting forces, i.e. if themotor is pulled forward, the effect of force is relatively small [40, 41]. In order to derivean expression similar to Eq. (9) for our model, we have to identify the proper term thatreplaces F t in Eq. (9). First, we rewrite Eq. (9) as v m = µ m ( F s − F t ), with some internal motor mobility coefficient µ m := v /F s . This version of Eq. (9) allows us to interpret the10 Sfrag replacements v F s v m F t F m F m,x v m v m | F m,x | U FIG. 2: Force velocity relation for a single motor: velocity v m as a function of load | F m,x | accordingto Eq. (10) with maximum velocity v and stall force F s . motor head as an over-damped (Stokesian) particle that constantly pulls with the force F s against some external load F t resulting in the effective velocity v m . According to Eq. (4) amotor pulls with force F m on the bead, so we can identify the “load” to be − F m,x , where F m,x is the x -component of F m and the minus sign accounts for Newton’s third law ( actio =reactio ). Thus, we obtain the following piecewise linear force velocity relation for the singlemotor (see also Fig. 2): v m = v F m · e x ≤ − | F m,x | F s if 0 < F m · e x < F s , F m · e x ≥ F s (10)where e x is the walking direction of the motor, see Fig. 1. Thus, if the motor pulls antiparallelto its walking direction on the bead, it walks with its maximum speed v . If it is loaded withforce exceeding the stall force F s , it stops. For intermediate loadings the velocity decreaseslinearly with load force. Eq. (10) defines the mean velocity of a motor in the presence ofloading force. We note that our algorithm also allows us to implement more complicatedforce-velocity relations and is not restricted to the piecewise linear force-velocity relation.It is used here because we do not focus on a specific kinesin motor and because it is easyto implement in the computer simulations. In practice the motor walks with discrete stepsof length δ . In the algorithm we account for this by defining a step rate λ s := v m /δ . Thedecision for a step during a time interval ∆ t is then made with the same Monte-Carlotechnique used to model the binding and unbinding process of the motor head. A step isrejected if the next binding site is already occupied by another motor (mutual exclusion)1142]. C. Bead versus motor friction
In principle, the velocity of a single motor v m pulling the sphere and the component inwalking direction of the sphere’s velocity U (see Fig. 1) are not the same. For an externalforce F t (against walking direction) in the pN range acting on the sphere and F m,x beingthe motor force in walking direction, the bead velocity U is given by U = µ ttxx ( F m,x − F t ),.Here, µ ttxx denotes the corresponding component of the mobility matrix M of the sphere(cf. Eq. (1)) evaluated at the height of the sphere’s center with super indices tt referringto the translational sector of the matrix and sub indices xx referring to the responses in x -direction. On the other hand from Eq. (10) it follows that the motor head moves withvelocity v m = µ m ( F s − F m,x ). Only in the stationary state of motion the two speeds areequal, U = v m , and we obtain the force with which the motor pulls (in walking direction)on the bead: F m,x = µ ttxx µ m + µ ttxx F t + µ m µ m + µ ttxx F s . (11)Thus, if the internal friction of the motor is large compared to the viscous friction of thesphere, i. e., 1 /µ m ≫ /µ ttxx , the second term in Eq. (11) can be neglected and one has F m,x ≈ F t . That means, only the trap force pulls on the motor. If µ m ≈ µ ttxx both terms inEq. (11) are of the same order of magnitude. Then, both external load F t and the frictionforce on the bead will influence the motor velocity. Experimentally, these prediction canbe checked in bead assays by varying the viscosities of the medium (e. g., by adding sugarlike dextran or Ficoll [43]). Numerically, we can vary η in the adhesive motor dynamicsalgorithm.When several motors are simultaneously pulling, they can cooperate by sharing the load.Assuming the case that n motors are attached to the MT which equally share the total load,then the force in x -direction exerted on the bead is nF m,x ( n ) with F m,x being again theindividual motor force. In the stationary state with U = v m we have (with external load12 t = 0): F m,x ( n ) = µ m nµ ttxx + µ m F s . (12)The mean bead velocity U ( n ) with n equally pulling motors is U ( n ) = nµ ttxx F m,x ( n ). Thus,in general, the velocity of the beads will increase with increasing number n of pulling motors.However under typical experimental conditions in vitro, where bead movements are probedin aqueous solutions, the internal motor friction 1 /µ m dominates over the viscous frictionof the bead, 1 /µ ttxx , and the velocity becomes independent of the number of motors. Onlyif the bead friction becomes comparable to the internal motor friction, the velocity exhibitsan appreciable dependence on the motor number. This dependence can be illustrated byconsidering the ratio of U ( n ) and U (1): U ( n ) U (1) = nµ ttxx + nµ m nµ ttxx + µ m = n/µ m + n/µ ttxx n/µ m + 1 /µ ttxx . (13)In the opposite limit, n/µ m ≪ /µ ttxx , i. e., when the viscous bead friction is very large anddominates over the internal motor friction, Eq. (13), leads to U ( n ) ≈ nU (1), and the beadvelocity increases linearly with n [7]. D. Vertical forces
We note that, although we are mainly interested in the x -component of the motor force F m , i.e. the component parallel to the microtubule, which enters the force–velocity relation,the motor force also has a component perpendicular to the microtubule, see Fig. 1. Thisforce component tends to pull the bead towards the microtubule and thus to the surface,whenever the bead is connected to the microtubules by a motor. This force is balanced bythe microtubule repelling the bead. Additionally, if the bead touches the filament or thewall, diffusion can only move the bead away from the wall. In case of the full spring model,compressed motors can also contribute a repulsion of the bead from the wall. If viscousfriction is strong, the normal component of the force arises mainly from the microtubule. Inthat case, the distance h between the bead and the microtubule is very small, h h i ≈
0. Forsmall viscous friction, thermal fluctuations play a major role and lead to non-zero distancesbetween the bead and the microtubule, as discussed further in Sec. III B.13ll parameters used for the adhesive motor dynamics simulations together with typicalvalues are summarized in Tab. I. For the numerical simulations we non-dimensionalize allquantities using R for the length scale, 1 /ǫ for the time scale and the detachment force F d as force scale. E. Master equation approach
When no external load is applied a motor walks on average a time 1 /ǫ before it detachesand the cargo particle might diffuse away from the MT. When several motors on the cargocan bind to the MT the mean run length dramatically increases as was previously shownwith a master equation approach [7]. For the sake of later comparison to our simulationswe briefly summarize some of these results in the following.Let P i be the probability that i motors are simultaneously bound to the MT with i = 0 , . . . , N m and N m being the maximum number of motors that can bind to the MTsimultaneously. Assuming the system of N m motors is in a stationary state and the totalprobability of having i = 0 , . . . , N m motors bound is conserved then the probability fluxfrom one state to a neighboring state is zero. This means that the probability P i of having i bound motors can be calculated by equating forward and reverse fluxes( N m − i ) π ad P i = ( i + 1) ǫP i +1 , i = 0 , . . . , N m − , (14)where it is assumed that the escape rate ǫ is a constant with respect to time. The solutionto Eq. (14) is given by P i = N m i ǫ N m − i π i ( ǫ + π ) N m , i = 0 , . . . , N m . (15)The probability that i motors are simultaneously pulling under the condition that at leastone motor is pulling is P i / (1 − P ) for i = 1 , . . . , N m . Then, the mean number of boundmotors N b (given that at least one motor is bound) is [7, Eq. [13]] N b = N m X i =1 iP i − P = ( π ad /ǫ ) [1 + π ad /ǫ ] N m − [1 + π ad /ǫ ] N m − N m . (16)14 Sfrag replacements l ha b FIG. 3: Illustration of the area fraction a b of the sphere (cut and placed besides the sphere) onwhich motor proteins can reach the MT (thin cylinder). a b depends in a geometrical fashion onthe minimum distance h between the sphere and the MT and on the resting length l . The effective unbinding rate ǫ eff , i. e., the rate with which the system reaches the unboundstate, is determined from ǫ eff (1 − P ) = π ad P . This quantity can also be identified with theinverse of the mean first passage time for reaching the unbound state, when starting withone motor bound [7]. If the medium viscosity is small, i. e., similar to that of water, and noexternal force is pulling on the bead, we assume that the velocity of the bead U does notdepend on the number of pulling motors. The mean run length, that is the mean distancethe cargo is transported by the motors in the case that initially one motor was bound, isthen the product of mean velocity U and mean lifetime (1 /ǫ eff ) [7, Eq. [14]]: h ∆ x b i = Uǫ eff = UN m π ad (cid:20)(cid:16) π ad ǫ (cid:17) N m − (cid:21) . (17)For kinesin-like motors at small external load with π ad ≫ ǫ this expression can be ap-proximated by h ∆ x b i ≈ ( U/N m ǫ )( π ad /ǫ ) N m − , i. e., the mean run length grows essentiallyexponentially with N m . In the stationary state the bead velocity U and the motor velocity v m are equal. For no external load and small viscous friction on the bead one can approximate ǫ ≈ ǫ in Eq. (16) and Eq. (17). F. Mean run length for a spherical cargo particle
In contrast to the master equation model in which one fixes the maximal number ofbound motors N m , in the computer simulations only the total number N tot of motors on thespherical cargo particle is fixed. A similar situation arises in experiments where only the total15mount of molecular motors on the sphere is measured [14] (but not in experiments withdefined multi-motor complexes [16]). If A b is the area on the sphere’s surface that includesall points being less than l apart from the MT (cf. Fig. 3), we expect on average n b = N tot a b motors to be close enough to the MT for binding, with the reduced area a b := A b / (4 πR ).While the master equation model assumes a fixed N m , in the simulations the maximumnumber of motors that can simultaneously bind to the MT is a fluctuating quantity aboutthe mean value n b . In the simulations the motors are uniformly distributed on the cargo,thus the probability distribution function P ( k ) for placing k motors inside the above definedarea fraction a b is a binomial distribution. As we have l ≪ R , a b is small and P ( k ) is wellapproximated by the Poissonian probability distribution function P ( k, n b ) = n kb k ! exp( − n b ) , n b = N tot a b . (18)Thus, given a fixed number N tot of motors on the sphere the number of motors that areinitially in binding range to the MT might be different from run to run. In addition,because of thermal fluctuations and torques that the motors may exert on the cargo particlethe relative orientation of the sphere changes during a simulation run. With the change oforientation also the number of motors in binding range to the MT is not a constant quantityduring one run.In order to make simulation results for the mean run length and the mean number ofbound motors comparable with the theoretical predictions Eq. (17) and Eq. (16), respec-tively, we have to average over different N m . Neglecting fluctuations in the number ofmotors that are in binding range to the MT during one simulation run we perform the aver-age with respect to the Poisson distribution, Eq. (18). Averaging the mean walking distance h ∆ x b i ( N m ) from Eq. (17) over all N m with weighting factors given by Eq. (18) we obtainthe following expression ( n b = a b N tot ): hh ∆ x b ii P oisson = N m,max X N m =1 n N m − b ( N m − UN m π ad "(cid:18) π ad ǫ (cid:19) N m − N m,max X N m =1 n N m − b ( N m − . (19)We note that during the initialization of each simulation run we place one motor on thelower apex of the sphere and then distribute the other N tot − P ( N m − , n b ) denotes the probability of having in total N m motors inside the areafraction a b . Furthermore, we introduced a cutoff N m,max of maximal possible motors ( N tot is obviously an upper limit for N m,max ). In the limit N m,max → ∞ we have hh ∆ x b ii P oisson = U ( e π ad n b /ǫ − / ( π ad n b ).In a similar way we can calculate the Poisson-averaged mean number of bound motors h N b i P oisson . For this it is important to include the correct weighting factor [44]. Froman ensemble of many simulation runs those with larger run length contribute more thanthose with smaller run length. For N simulation runs the mean number of bound motorsis obtained as h N b i sim = P i t i n ( i ) / P i t i , where t i is a period of time during which n ( i )motors are bound and the sum is over all such periods of time. Assuming the bead velocity U to be a constant, the time periods t i can also be replaced by the run lengths ∆ x b,i during t i . Picking out all simulation runs with a fixed N m , their contribution to the sum is themean number of bound motors N b times the total run length of beads with given N m . Thelatter is the mean run length h ∆ x b i N m times the number of simulation runs with the given N m (for sufficiently large N ). Clearly, the fraction of runs with given N m is the probability P ( N m − , n b ) introduced in Eq. (18). Consequently, we obtain again with a truncation ofthe sum at some N m,max ≤ N tot h N b i P oisson = N m,max X N m =1 ¯ P ( N m ) N b ( N m ) = P N m,max N m =1 P ( N m − , n b ) h ∆ x b i N m N b ( N m ) P N m,max N m =1 P ( N m − , n b ) h ∆ x b i N m . (20)Here we introduced the probability ¯ P ( n ) for having n motors in binding range to the MTwhen picking out some cargo particle from a large ensemble of spheres. Explicitely, thisprobability is given by¯ P ( n ) = Uπ ad ((1 + π ad /ǫ ) n − n n − b n ! e − n b P N m,max i =1 Uπ ad ((1 + π ad /ǫ ) i − n i − b i ! e − n b = ((1 + π ad /ǫ ) n − n n − b /n ! P N m,max i =1 ((1 + π ad /ǫ ) i − n i − b /i ! . (21) G. Computer simulations
We use the following procedure for the computer simulations. In each simulation run thesphere is covered with N tot motors. Initially, one motor, located at the lowest point of thesphere, is attached to the microtubule such that the distance of closest approach h between17he sphere and the microtubule is given by the resting length of the motor, i. e., h = l . Theother ( N tot −
1) motors are uniformly distributed on the sphere’s surface (cf. Sec. B). Whenthe motor starts walking, it pulls the sphere closer to the MT because there is a z -componentin the force exerted on the sphere by the motor stalk (which is strained after the first step).Then, other motors can bind to the MT. The system needs some time to reach a stationarystate of motion, so initially the motor velocity v m and the bead velocity U are not the same(for reasons of comparison a fixed initial position is necessary; other initial positions havebeen tested but initialization effects were always visible). In principle a simulation run lastsuntil no motor is bound. For computational reasons, each run is stopped after 2 · s(which is rarely reached for the parameters under consideration). For each run quantitieslike the mean number of bound motors N b , the walking distance ∆ x b and the mean distanceof closest approach h h i between sphere and MT are recorded. III. RESULTSA. Single motor simulations
In Sec. II B we defined a force-velocity relation for the single motor. In this section weperform simulation runs with a single motor, i. e., N tot = 1, and measure the effective forcevelocity relation by tracking the position of the sphere. Inserting Eq. (11) into Eq. (10)provides a prediction for the velocity of a bead subject to one pulling motor and an externaltrap force F t . In Fig. 4a this prediction is shown for three different values of viscositytogether with the actual measured velocities during the simulations. More precisely, wemeasured the mean velocity of the bead and the motor obtained from a large number ofsimulation runs (to avoid effects resulting from the initial conditions we first allowed therelative position/orientation of bead and motor to “equilibrate” before starting the actualmeasurement). The mean velocity is then given as the total (summed up over all simulationruns) run length divided by the total walking time. The good agreement between thenumerical results and the theoretical predictions provides a favorable test to the algorithm.At η = 1 mPa s (the viscosity of water), friction of the bead has almost no influence onthe walking speed. At hundred times larger viscosities, however, bead friction reduces themotor speed to almost half of its maximum value already at zero external load. Although18 PSfrag replacements U , v m i n µ m / s trap force F t in pN η = 1 mPas η = 10 mPas η = 100 mPas PSfrag replacements U , v m i n µ m / s trap force F t in pN (a) (b) FIG. 4: (a) Measured force velocity relation of a single motor (with l = 80 nm) pulling a sphereof radius R = 1 µ m for three different viscosities η = 1 , ,
100 mPa s. Shown is the relationaccording to Eq. (9), the actual measured force velocity relation of the motor head and the beadcenter, respectively, and the theoretical prediction according to Eq. (10) and Eq. (11). (b) Themeasured force velocity relation for η = 1 mPa s is shown where in Eq. (10) not | F m,x | but k F m k isused. The dotted line emphasizes the linear decrease of the velocity. The negative velocity of thebead at large F t results from thermal fluctuations. Fluctuations against walking direction increasethe escape probability. In case of escape they cannot be compensated by fluctuations in walkingdirection. (Numerical parameters: ∆ t = 10 − , number of runs N = 2 · − · .) the velocities of the motor and the bead are expected to be equal, Fig. 4a shows that themotor is slightly faster than the bead. This is a result of the discrete steps of the motorand can be considered as a numerical artefact: at the moment the motor steps forward themotor stalk is slightly more stretched (loaded) than before the step, therefore, the escapeprobability is increased. The result of unbinding at the next time step would then be that thebead moved a distance δ less than the motor. For loads close to the stall force the observedvelocity is slightly larger than the prediction, which is a result of thermal fluctuations of thebead in combination with the stepwise linearity of the force velocity relation: a fluctuationin walking direction slightly reduces the load on the motor, thus increasing the step rate,whereas fluctuations against walking direction lead to zero step rate.It was observed by Block et al. that vertical forces on the bead (i. e., in z -direction) alsoreduce the velocity of the motor [41]. But the same force that leads to stall when appliedantiparallel to the walking direction has a rather weak effect on the motor velocity whenapplied in z -direction. Using the absolute value of the total force of the motor k F m k inEq. (10) instead of its x -component | F m,x | , we measure a force velocity relation as shown in19 -3 -2 -1
0 1 2 3 4 5 6step length10 -3 -2 -1
0 1 2 3 4 5 6
PSfrag replacements η = 1 mPas η = 10 mPas η = 100 mPas η = 1 mPas, χ = 0 ◦ η = 100 mPas, χ = 60 ◦ η = 1 mPas, χ = 60 ◦ η = 10 mPas, χ = 60 ◦ h ∆ x b i i n µ m trap force F t in pN FIG. 5: Mean run length h ∆ x b i of a bead pulled by a single motor as a function of an external forceon the bead F t and for three different viscosities η = 1 , ,
100 mPa s. The lines give the theoreticalpredictions according to Eq. (22) assuming an angle of 60 ◦ between the motor and the MT. Forcomparison also the theoretically predicted h ∆ x b i -curve for χ = 0 is shown (double dotted line). Fig. 4b. Again, the velocity decreases essentially linearly with applied external force, butstalls already at around F t ≈ F s / k F m k . Asthe effect of vertical loading reported in Ref. [41] seems to be much weaker than that shownin Fig. 4b, we reject this choice of force velocity relation.From the simulations carried out for Fig. 4a we can also obtain the mean run length h ∆ x b i for a single motor as a function of external load. The results are shown in Fig. 5.Using Eq. (10) and the Bell equation, Eq. (8), we obtain h ∆ x b i = v m ǫ = v ǫ − | F m,x | /F s exp( k F m k /F d ) . (22)The numerical results shown in Fig. 5 fit very well to the theoretical prediction of Eq. (22)when assuming the angle χ between the motor and the MT to be χ = 60 ◦ . The angle χ depends on the bead radius R , the resting length l [45] and the polymer characteristics ofthe motor protein, e. g., its stiffness κ . B. Run length for several motors pulling
We now measure the run length distributions and the mean run length as a function ofmotor coverage N tot . For motors modeled as springs according to Eq. (5) with two differentresting lengths l = 50 ,
65 nm the run length distributions are shown in Fig. 6a,b. For each20 f r equne cy PSfrag replacements walking distance ∆ x b in µ m N tot = 120 N tot = 200 N tot = 280 N tot = 360 N tot = 440 N tot = 520 f r equne cy PSfrag replacements walking distance ∆ x b in µ m N tot = 120 N tot = 200 N tot = 280 N tot = 360 N tot = 440 (a) (b) FIG. 6: Distribution of run lengths ∆ x b in semi-logarithmic scale for different values of motorcoverage N tot . The motor protein is modeled as a harmonic spring according to Eq. (5). (a)Resting length of the motor protein l = 50 nm. (b) Resting length of the motor protein l = 65 nm.(Numerical parameters: time step ∆ t = 10 − /ǫ , number of simulation runs N ≈ .) value of N tot the run length was measured about N = 10 times. The simulations turnout to be very costly, especially for large N tot as the mean run length increases essentiallyexponentially with the number of pulling motors (cf. Eq. (17)). From Fig. 6 we see that thelarger N tot , the more probable large run lengths are, resulting in distribution functions thatexhibit a flatter and flatter tail upon increasing N tot . Fig. 6 nicely shows that the shape ofthe distributions depends not only on the total number of motors N tot attached to the spherebut also on the resting lengths. Given the same N tot we can see that longer run lengths aremore probable for the longer resting length l = 65 nm shown in Fig. 6b. This can simplybe explained by the fact that the larger the motor proteins, the larger is the area fraction a b introduced in Fig. 3 and therefore the more motors are on average close enough to themicrotubule to bind.Fig. 7a shows the mean run length as a function of N tot as obtained by numerical simu-lations of the transport process (points with error bars). For the motor stalk three differentvalues of the resting length l = 50 , ,
80 nm are chosen and both the full-spring and thecable model are applied for the force extension relation. Similarly to what we have alreadyseen for the run length distributions in Fig. 6, the larger the resting length l the more mo-tors can simultaneously bind for given N tot , and therefore the larger is the mean run length.Furthermore, Fig. 7a also demonstrates that it makes a clear difference whether the motorstalk behaves like a full harmonic spring or a cable. If the motor protein behaves like a cable21
0 80 160 240 320 400 480 560fit for area-fraction
PSfrag replacements total number of motors N tot h ∆ x b i i n µ m l = 50 nm, spring l = 50 nm, cable l = 65 nm, spring l = 65 nm, cable l = 80 nm, spring l = 80 nm, cable PSfrag replacements total number of motors N tot N b l = 50 nm, spring l = 50 nm, cable l = 80 nm, spring l = 80 nm, cable (a) (b) FIG. 7: (a) Mean run length h ∆ x b i (data points with error bars) as a function of motors on thebead N tot obtained from adhesive motor dynamics. The lines are fits of Eq. (19) with respect tothe area fraction a b . (b) Mean number of bound motors N b (data points with error bars). Thelines are the values obtained from the Poisson-averaged mean number of bound motors h N b i P oisson in Eq. (20) using for a b the fit value from (a). (Parameters: π ad = 5 , ǫ = 1, λ s = 125, ∆ t = 10 − , N ∼ .) (semi-harmonic spring, Eq. (6)) it exhibits force only if it is stretched. The vertical compo-nent of this force always pulls the cargo towards the MT. Thus, the mean height betweenthe cargo and the MT (which determines how many motors can bind at maximum) resultsfrom the interplay between this force and thermal fluctuations of the bead. In contrast, ifthe motor also behaves like a harmonic spring when compressed, it once in a while may alsopush the cargo away from the MT. This results in less motors being close enough to theMT for binding than in the case of the cable-like behavior of the motor stalk. Consequently,given the same l and N tot , the cargo is on average transported longer distances when pulledby “cable-like motors”.In order to apply the theoretical prediction for the mean run length of a sphericalcargo particle, Eq. (19), we need to determine proper values for the three parameters a b = n b /N tot , U, N m,max . From the simulations we measure the mean velocity of the sphere U . It turns out that U is up to 15 % less than the maximum motor velocity v due togeometrical effects. Depending on the point where the motor is attached on the sphere,some motor steps may result mainly in a slight rotation of the sphere instead of transla-tional motion of the sphere’s center of mass equal to the motor step length δ . For N m,max we choose the overall measured maximum value from all simulation runs for given N tot andpolymer model of the motor. Then, we use the remaining parameter, the area fraction a b , as22 , motor-model fit value for a b measured h h i → a b ( h h i )50nm, spring Eq. (5) 0.00211 7-14 nm → → → → TABLE II: Obtained fit values for the area fraction a b for different l and the two applied polymermodels. For comparison the area fraction which is obtained from the measured mean distance h h i is also displayed. h h i is measured for fixed N tot , the left boundary of the provided intervalcorresponds to the largest N tot . a fit parameter to the numerical results. The fits are done using an implementation of theMarquardt-Levenberg algorithm from the Numerical Recipes [46]. The resulting a b valuesare summarized in Tab. II. The theoretical curves for those parameter values are shown(dashed lines) in Fig. 7a. The increase in a b for larger resting length and the cable model isin excellent accordance with the above discussed expectation. An independent estimate ofthe area fraction can be obtained by measuring the mean distance h h i between cargo andMT and calculating the area fraction a b as a b = a b ( h h i ). For comparison those values are alsogiven in Tab. II. They turn out to be about 60 % larger than the values for a b obtained fromthe fit. This indicates that the height of the sphere above the MT (and therefore also thearea fraction) is a fluctuating quantity that is not strongly peaked around some mean value.Then, because of the non-linear dependence of a b on h we have in general h a b ( h ) i 6 = a b ( h h i ).Fig. 7b shows the mean number of bound motors (the average is obtained over all N simulation runs) as a function of N tot (symbols with error bars). The lines in Fig. 7b areplots of Eq. (20) using the same parameters as for the correspondig lines in Fig. 7a. Onerecognizes that again the theoretical prediction and the measured values match quite well.This means that on the level of the mean run length and mean number of bound motorsthe Poission average that was introduced in Sec. II F works quite well, even though thenumber of motors that are in range to the MT is not a constant during one simulation run(cf. Sec. III C). The large error bars for the cable model data in comparison to the springmodel results partly from a poorer statistics (for the l = 80 nm cable simulations thenumber of runs is in the range of some hundreds only). In addition, for cable-like motors thefluctuations in the area fraction a b are much larger than for spring-like motors as repulsive23
0 2 4 6 8 10theory
PSfrag replacements mean number of bound motors N b h ∆ x b i i n µ m l = 50 nm, spring l = 50 nm, cable l = 65 nm, spring l = 65 nm, cable l = 80 nm, spring l = 80 nm, cable FIG. 8: Combination of data from Fig. 7a,b for the resting lengths l = 50 , ,
80 nm. h ∆ x b i isshown as a function of N b . For the dashed line (theory) Eq. (19) and Eq. (20) were combinedwith a truncation of the sums at N m,max = 18. (Parameters: π ad = 5 s, ǫ = 1 s, λ s = 125 s,∆ t = 10 − /ǫ , N ∼ .) spring forces stabilize the distance between the sphere and the MT. Therefore, the width ofthe distribution density of the number of bound motors is larger for the cable model thanfor the spring model. Fig. 7b also shows that for a bead radius of 1 µ m, around 80 motorshave to be attached to the bead, otherwise binding and motor stepping become unstablebecause there are less than two motors left in the binding range.In Fig. 8 the simulation results of Fig. 7a,b are combined into one plot. Here we showthe measured mean run length h ∆ x b i as a function of the mean number of bound motors N b . All curves collapse on one master curve that can be parametrized by n b = a b N tot , i. e.,the product of the fit parameter a b and the total number of motors on the sphere. Thefact that all data points turn out to lie on one master curve again demonstrates the goodapplicability of the theoretical predictions to the simulation results. The curve shown inFig. 8 has a positive curvature in the semi-logarithmic plot. This turns out to be an effectof the finite truncation of the sums in Eq. (19) and Eq. (20). C. Distribution of motors in binding range
For the evaluation of the numerical data in Sec. III B we assumed that the number ofmotors that are in binding range to the MT are constant in time and that this number isdrawn from a Poisson distribution for every individual run. We now further examine this24 f r eguen cy number of motors in range to MTPoisson PSfrag replacements N tot = 80 N tot = 120 N tot = 160 N tot = 200 N tot = 240 N tot = 440 N tot = 520 f r eguen cy number of motors in range to MTPoisson PSfrag replacements N tot = 80 N tot = 120 N tot = 160 N tot = 200 N tot = 240 N tot = 400 N tot = 440 (a) (b) FIG. 9: Histograms for the number of motors that are in binding range to the MT. Symbols referto simulation results for different values of the total number of motors on the sphere N tot . Lines arePoisson distributions with mean value that is propotional to N tot . (a) Resting length l = 50 nm,spring model, Eq. (5). (b) Resting length l = 50 nm, cable model, Eq. (6). (Parameters: π ad =5 , ǫ = 1, λ s = 125, ∆ t = 10 − /ǫ , N = 2 · .) aspect in order to see to what extend this assumption is fulfilled. First, we measure directlythe distribution of the number of motors n MT in binding range to the MT. For this wecount n MT at every numerical time step during one simulation run and repeat this for alarge ensemble of runs ( N ∼ ). Thus, the histograms (i. e., approximately the probabilitydistributions we are looking for) that are obtained in this way for the relative frequency of n MT are not based on the assumption of constant n MT for every single run.Fig. 9 shows examples of such histograms (symbols) for a series of different values of thetotal number of motor coverage N tot . For Fig. 9a we used the spring model for the motorpolymer, Eq. (5), and for Fig. 9b the cable model, Eq. (6). In both cases the resting lengthof the motor protein is l = 50 nm. In addition, Fig. 9 also displays probability distributions(solid lines) that are obtained from the Poission distribution given that at least one motoris in binding range, i. e., p ( n ) = µ n e − µ /n !(1 − e − n ), with mean value µ given by µ = µ N tot that can be parametrized by some variable µ . The parameter µ was chosen to be 0.015 and0.0165 for the spring and cable model, respectively, by matching the Poisson distributionto the simulation data. For the spring model (Fig. 9a), the fit is excellent for all values of N tot . One must note however that these distributions are not given by Eq. (18) as indicatedby the fact that the parameter µ is much larger than the area fraction a b determined inthe previous section. Instead one needs to account for the correct weighting factors from25 f r eguen cy number of motors in range to MT PSfrag replacements N tot = 80 N tot = 120 N tot = 160 N tot = 200 N tot = 240 N tot = 400 N tot = 440 FIG. 10: The relative frequencies of the number of motors that are in binding range to the MTduring a single run are shown for six different sample runs. For the motor protein the cable modelwith resting length l = 50 nm is used. (Parameters: N tot = 400, π ad = 5 s, ǫ = 1 s, λ s = 125 s,∆ t = 10 − /ǫ ). the different run lengths. If for the moment we consider the number n MT to be a constantduring one run, then the distribution Eq. (21) can be considered to be a useful estimate forthe data in Fig. 9. Taking the limit N m,max → ∞ in Eq. (21) we obtain( a b (1 + π ad /ǫ )) n MT − a n MT b e a b π ad /ǫ − e − a b n MT ! ≈ b n MT n MT ! e − b , (23)with b = a b (1 + π ad /ǫ ). Thus, the result is approximately—except for very small n MT =1 , b . With π ad /ǫ = 5 and a b = 0 . b = 0 .
013 which is very close to the value µ = 0 .
015 used in Fig. 9a.For the cable model (Fig. 9b), the Poisson fit using a single value for µ works well onlyfor the smaller values of N tot . For large N tot the data points cannot be fitted by a Poissondistribution. It rather turns out that the ratio of mean value and standard deviation becomesless than one in these cases.Finally, we shall note that despite the good agreement between the simulation data andthe the estimate from Eq. (23), which was based on the assumption that n MT is constantduring one run, n MT is in fact not constant, but a fluctuating quantity. The fluctuationsresult partly from thermal fluctuations of the height and orientation of the sphere and partlyfrom orientation changes of the sphere that are induced by motor forces. Fig. 10 shows a fewsample histograms for the frequency that n MT motors are in binding range to the MT, i.e.either bound to the MT or unbound within the binding range, during a single run. These26 f r equen cy PSfrag replacements N tot = 80 N tot = 120 N tot = 160 N tot = 200 N tot = 240 N tot = 280 N tot = 320 N tot = 360 N tot = 440 N tot = 480 N tot = 640 escape rate ǫ/ǫ FIG. 11: Measured probability distribution density for the escape rate ǫ given ǫ > ǫ . The datawas obtained for different values of N tot . For the motor proteins the full spring model, Eq. (5),with resing length l = 50 nm was used. (Parameters: π ad = 5 , ǫ = 1, λ s = 125, ∆ t = 10 − /ǫ , N = 2 · .) examples clearly show that n MT takes different values during one run. D. Escape rate distributions
Diffusive motion of the cargo sphere directly influences the length of the pulling motorproteins and therefore also the escape rate ǫ . The dependence of the escape rate on themotor length x for x > l is obtained by inserting the polymer model Eq. (6) and Eq. (5),respectively, into the Bell equation, Eq. (8). For x ≤ l the escape rate is given by ǫ = ǫ .For low viscous friction we assume x to be distributed according to a Gaussian distribution.Then, we expect the probability distribution density ˜ p ( ǫ ) for the escape rate ǫ to be givenby a log-normal distribution density˜ p ( ǫ ) = 1 bǫ exp (cid:18) F d (ln( ǫ ) − ln(¯ ǫ )) k B T κ ∗ (cid:19) . (24)Here, ¯ ǫ denotes the escape rate associated with the mean motor length in the extended state, κ ∗ is an effective spring constant that depends, e. g. on the number of pulling motors, and b is a normalization constant that is obtained from the condition Z ∞ ǫ ˜ p ( ǫ ) dǫ ! = 1 . ǫ > ǫ and different values of N tot . Itturns out that the log-normal distribution in Eq. (24) matches well to the measured data(not shown). Fitting Eq. (24) to the data for the effective spring constant κ ∗ it turns outthat κ ∗ increases with increasing N tot . This makes sense and illustrates that for severalmotors pulling the motors behave as a parallel cluster of springs.The agreement of Eq. (24) with the simulation data for the escape rate distributionssuggests that the main source of force acting on the motors and increasing the unbindingrate is due to thermal fluctuations of the micron-sized bead. This is different from what hasbeen reported in experiments with nano-scaled two-motor complexes [16]. There it has beenargued that the forces between the two motors arising from their stochastic stepping leadto an increased unbinding rate of these motors. E. Cargo transport against high viscous friction
Except for the single motor simulations of Sec. III A, all simulation data discussed sofar were obtained for a viscosity of 1 mPa s, corresponding to a water-like solution. Wementioned in Sec. II B that load sharing between several motors may lead to cooperativeeffects at high viscous friction. We now analyze this further by performing simulations atviscosities much larger than that of water (i. e., when the viscous friction on the bead iscomparable to the internal friction of the motor protein). To do this we need to measure thevelocity of the bead depending on the number of pulling motors. Because of the nature ofthe stochastic process describing the position of the cargo the instantaneous cargo velocityis however not well defined [47]. Therefore, in order to measure the cargo velocity U we haveto average over some time interval ∆¯ t . If no motors were pulling the velocity distributiondensity is given by a Gaussian, p ( U, ∆¯ t ) = r ∆¯ t πD e − U ∆¯ t/ D , (25)with diffusion constant D = k B T a µ ttxx (cf. Sec. II A). Furthermore we assume a constantheight of the sphere so that the mobility coefficient µ ttxx is a constant in time. The width ofthe distribution density, Eq. (25), is the smaller the larger ∆¯ t is. So in order to suppressfluctuation effects it seems appropriate to average over a large time interval ∆¯ t . On the other28 r e l a t i v e f r equen cy -3 -2 -1 PSfrag replacements velocity in m/s η = 1 mPa s η = 10 mPa s η = 100 mPa s m ean v e l . i n µ m / s velocity in µ m/s viscosity η in Pa s PSfrag replacements = 10 mPa svelocity in m/s m ean v e l . i n µ m / s nη = 1 m Ps η = 10 m Ps η = 100 m Ps η = 1 m Ps (theory) η = 10 m Ps (theory) η = 100 m Ps (theory) (a) (b) r e l a t i v e f r equen cy PSfrag replacements η = 100 mPa svelocity in µ m/s n = 1 n = 2 n = 3 n = 4 mean vel. in m/s n (c) FIG. 12: (a) Propability density of the cargo particle’s velocity U that is obtained by averagingover a time interval of ∆¯ t = 0 .
02 s for different values of the viscosity η . The inset shows themean velocity as a function of the η . (b) The mean velocity of the bead given that n motorsare simultaneously pulling is plotted as a function of n (symbols). For comparison the theoreticalexpectation according to Eq. (13) is plotted, too (lines). (c) For high viscosity η = 100 mPas theconditional velocity distribution density given that a certain number of motors is pulling is plotted.(The full harmonic spring model was used; parameters: l = 80 nm, N tot = 200, numerical timestep ∆ t = 10 − , other parameters as in Tab. I.) hand the number of pulling motors changes with time because of binding and unbinding.Thus, in order to measure the velocity given a certain number of motors ∆¯ t should not betoo large in order to get enough such events. Here we choose ∆¯ t = 0 .
02 s which correspondsto a typical camera resolution of 50 Hz.Fig. 12a shows the measured velocity distributions for three different values of the viscos-ity, η = 1 , ,
100 mPa s. In the inset of Fig. 12a the mean velocity is plotted as a function ofthe viscosity. It turns out that shifting the distribution Eq. (25) by the corresponding meanvelocity, the single peaked function p ( U, ∆¯ t ) fits qualitatively well to the distribution shown29n Fig. 12a, especially the dependence of the width of the distributions on η is correctlypredicted by Eq. (25). Also a decrease of the mean velocity of the cargo particle is observedwith increasing viscosity resulting from the increased frictional load. It must be emphasizedthat only a single peak is observed in Fig. 12a, even though the bead velocity is expectedto depend on the number of pulling motors, which should lead to multiple distinct peaks[7]. One reason that we do not observe multiple peaks is the small value of the time interval∆¯ t , which leads to broad peaks with a peak width governed by diffusion of the bead (cf.Eq. (25)) and thus makes it impossible to separate different peaks. Using larger values of∆¯ t leads to smaller peak widths, but also to poorer statistics as less measurement points areobtained, so that again distinct peaks cannot be resolved. Even if we do not use a constanttime interval, but average over the variable time intervals between two subsequent changesin the number of bound motors [9], distinct peaks are very hard to separate (not shown).This does however not mean that the bead velocity is independent of the number of pullingmotors. Indeed if we plot the conditional velocity distribution calculated over all intervalsin which the bead is pulled by a certain fixed number of motors, we see a clear shift in theaverage velocity (Fig. 12c). This shift is however masked by the width of the distributionsin Fig. 12a.In Fig. 12b the average of all measured velocities given that exactly n motors are pullingis plotted as a function of n , again for the three different viscosities η = 1 , ,
100 mPa s.For η = 1 mPa s the viscous friction for the bead is about 1 /µ ttxx ≈ · − Ns/m. Theinternal friction of the motor is 1 /µ m = F s /v ≈ · − Ns/m, i. e., about two orders ofmagnitude larger than 1 /µ ttxx . According to the analysis at the end of Sec. II B we thereforeexpect that the mean velocity is independent of n if all motors equally share the load. Thenumerical data in Fig. 12b shows that the mean velocity exhibits a weak dependence on n with a maximum at about n = 4 ,
5. At the higher viscosities the numerical results showthat the mean bead velocity increases with increasing n , which indicates that the motorsshare the load. The simulation data however deviate clearly from the estimate given byEq. (13), which is indicated by the lines in Fig. 12b. This discrepancy indicates that theload is not shared equally among the motors or that only a subset of the bound motors areactually pulling the bead. Besides geometrical effects one reason why this is the case is thatthe escape rate ǫ is rather high and at high frictional load is even further increased makingthe lifetimes of motors in the pull state rather short. Then, if a new motor binds to the MT30 f r eguen cy relative load force 2 motors3 motors4 motors5 motors6 motors7 motors8 motors PSfrag replacements f r eguen cy relative load force 2 motors3 motors4 motors5 motors PSfrag replacements (a) (b) f r eguen cy relative load force PSfrag replacements f r equen cy relative load force PSfrag replacements (c) (d)
FIG. 13: Frequencies of the motor forces in walking direction relative to the total load force onthe cargo particle for different numbers of pulling motors. (a) Viscosity η = 1 mPa s, escape rate ǫ = 1 s − . (b) η = 100 mPa s, ǫ = 1 s − . (c) η = 100 mPa s, ǫ = 0 . − . (d) η = 1000 mPa s, ǫ = 0 . − . For the other parameters the same values as in Fig. 12 were used. often another motor detaches already before a stationary state is reached in which all themotors equally share the load.In order to investigate the last point in more detail we consider explicitely how the forceis typically distributed among the pulling motors. For this we count the number n of motorsattached at each numerical time step and measure the force experienced by each motor inthe direction along the microtubule. For a given number n , n such motor forces can bemeasured, F ( i ) m,x , i = 1 , . . . , n . To suppress effects from fluctuations in the overall load wethen calculate the reduced forces f n := F ( i ) m,x / P ni =1 F ( i ) m,x . Given the histograms for thesequantities measured over many time steps and simulation runs we obtain an approximationfor the probability distribution density of the relative load of the motors. Fig. 13 shows31esults for such histograms obtained at different viscosities and two different values of theunstressed escape rate ǫ . Fig. 13a shows the results for η = 1 mPa s and ǫ = 1 s − .One can see some symmetries that result from the definition of the reduced forces. Thisis especially emphasized for the case of n = 2 pulling motors. The distribution density ofthe reduced forces has a mirror symmetry about the value f = 0 .
5. Besides this artefact itturns out that for n > f = 1 / f = 1 /
3, which are rather broad and thushard to resolve. On the other hand there is still a strongly peaked maximum at f n = 0.This somewhat surprising observation is due to the rather high escape rate, which is evenincreased by the load force (cf. Eq. (8)). Therefore, the binding time of the motors isshortend. On the other hand motors that bind to the MT are initially unstressed (i. e., carryzero load) in our model. Thus they always contribute to the f n = 0 peak and may alreadyescape from the MT before the load is shared equally by all pulling motors. When theescape rate is reduced to ǫ = 0 . − as done for Fig. 13c,d clearly visible maxima around f n = 1 /n appear in the histograms that indicate that the load is equally shared by the activemotors. In Fig. 13d where we used the extremely high value of 10 mPa s for the viscosity,these peaks are very pronounced. The arrows in Fig. 13c,d indicate the relative force values1 /i, i = 2 , , . . . . It turns out that the peaks are not exactly located at these values which isagain due to the binding and unbinding process of the motors.In summary, we found that at high loads the pulling motors tend to arrange in such away that the total load is equally shared amongst them. However, for typical escape ratevalues of kinesin-like motors this process often takes more time than the lifetime of a stateof a certain number of motors bound to the MT lasts, thus preventing cooperativity in thesense of equal load sharing. 32 V. DISCUSSION AND OUTLOOK
The main purpose of this paper is to introduce a novel algorithm called adhesive motordynamics as a means to study the details of motor-mediated cargo transport. Our algo-rithm is an extension of existing adhesive dynamics algorithms developed to understand thephysics of rolling adhesion [22, 26]. Basically our method allows to simulate the motionof a sphere above a wall including hydrodynamic interactions and diffusive motion by nu-merically integrating the Langevin equation, Eq. (3). In addition, motor-specific reactionssuch as binding to the microtubule and stepping are modeled as Poisson processes and thentranslated into algorithmic rules. The parameters and properties by which the motors aremodeled are based on results of single-molecule experiments with conventional kinesin. Afirst favorable test for the algorithm was provided by measuring the force velocity–relationat different viscosities and external load forces and by comparing the results to the inputdata as done in Sec. III A.Next we measured the run length and the mean number of bound motors as a functionof the total number N tot of motors attached to the sphere. The same quantities have beenpreviously calculated based on a one-step master equation model [7]. However, this hasbeen done as a function of a fixed number of motors that are in binding range to themicrotubule. In practise and also in our simulations, this quantity varies in time. In Sec. II Fwe Poisson-averaged the theoretical predictions thus rendering it possible to compare theoryand simulation results. Using the area fraction on the cargo from which the microtubule isin binding range for the motors as a fit parameter, we found good agreement between theoryand simulations for both the mean run length and the mean number of bound motors. Notethat the latter one cannot be measured in typical bead assay experiments.We also determined the mean separation height between cargo and microtubule andfound h h i = 4 −
14 nm. Modeling the motor stalk as a cable resulted in smaller distancesthan using a full spring model for the motor stalk. A recent experimental study usingfluorescence interference contrast microscopy found that kinesin holds its cargo about 17 nmaway from the MT [37]. Our smaller distance probably results from neglecting any kindof volume extension (except binding site occupation) of the motor protein, the simplifiedforce extension relation applied to model the stalk behavior, and neglecting electrostaticrepulsions. These effects, however, could easily be included into our algorithm, e. g., by using33ard core interactions that account for the finite volume of the motor protein segments.In Sec. III C we explicitely demonstrated that the theoretical assumption of having a fixednumber of motors in binding range during one walk is not justified (cf. Fig. 10). Neverthelessthe theoretical results agree well with the simulations. This might be explainable by theobservation that averaged over many runs the distribution of motors in binding range appearsto be Poissonian (cf. Fig. 9). Thus, on the one hand fast fluctuations in the number of motorsin binding range around some mean value are not visible. On the other hand periods in whichthis number fluctuates around the same mean value can be treated as a complete run. Thus,averaging over these smaller runs (i. e., which end after the sphere was e. g. rotated visiblyand not after the last motor unbinds) has the same effect as averaging over complete runs(i. e., which end after the last motor unbinds).An interesting question is to what extend several motors can cooperate by sharing load.We have addressed this question for the case of several motors pulling a cargo particleagainst high viscous friction. One of the advantages of our algorithm is that we can measurethe velocity of the bead and at the same time also the number of simultaneously pullingmotors. Thus, in Sec. III E we tried to check whether the explanation of Ref. [9] is correct alsounder the assumptions of our model, especially for the parameter values given in Tab. I. Oursimulations show that the speed of the cargo increases with the number of pulling motors forhigh viscous friction, in agreement with experimental results [48]. Our simulations howevershow pronounced deviations from the quantitative predictions based on the assumptionthat the load is shared equally among the bound motors. Furthermore, as the averagelife time of a state with a certain number of pulling motors is rather short the differentvelocities expected for different numbers of instantaneously pulling motors were smearedout by diffusion. Similarly when directly measuring how the total load force is distributedto the different motors pulling, no equal load sharing could be observed for the escape rateof about 1 s − . We observed equal load sharing only when we used a ten-fold smaller escaperate, in order to increase the life time of the motors in the bound state.Another interesting question in this context is whether the velocity distribution exhibitsseveral maxima if the cargo is pulled against a viscous load, as observed in several experi-ments in vivo [9, 49]. For example, Hill et al. [9] found that vesicle in neurites move withconstant velocity for some period of time and then switch to another constant velocity in astep-like fashion. The distribution of velocities (measured over time intervals of the order of34 s) was found to have peaked at integer multiples of the minimal observed velocity. Thesepeaks were interpreted as corresponding to different numbers of simultaneously pulling mo-tors, which equally share the visoelastic load excerted by the cytoplasm [9] (cf. also Eq. (13)).Indeed, both an earlier model for motor cooperation [7] and our present description predictthat equal sharing of a large viscous load leads to such a velocity distribution. In our simu-lations, we could however not resolve multiple peaks, presumably because the peaks are toobroad to be resolved. The latter results from a combination of the way how we measure thevelocity and from the fast dynamics of motor unbinding as discussed in Sec. III E.As already mentioned above, the framework of our method is rather general. Thereforevarious model variations can be easily implemented and probed in simulations. Here, forexample we modeled the motor stalk by two versions of a simple harmonic spring: thecable model, which represents a molecule with an intrinsic hinge, and the spring model.More advanced force–extension relations could easily be incorporate in Eq. (4) in order toprobe the influence of more realistic polymer models on the transport process. Similarly,the force dependence in unbinding from the microtubule and stepping can be altered toexplore the impact onto macroscopic observables like the mean run length or the speedof the cargo. Furthermore, the algorithm could easily be adapted to study beads to whichclusters of motors or defined motor complexes (such as those in ref. [16]) are attached. Thus,the algorithm described in this paper can be regarded as a link between purely theoreticalmodels and data from in vitro experiments.Another interesting question for future applications of our method is how cargo transportworks against some external shear flow. Since our model is based on a hydrodynamicdescription, flow can easily be included in the dynamics of our model. For these studiesthe Langevin-equation, Eq. (1), has to be extended by additional terms accounting for theeffect of an incident shear flow [23, 28]. Then, two opposing effects exist characterized bythe step rate and the strength of the shear flow, respectively. Their interplay together withthe rates for binding and unbinding π ad and ǫ , respectively, determine whether the cargomoves in walking direction or in flow direction. Experimentally, such a setup might provideinteresting perspectives for biomimetic transport in microfluidic devices.35 cknowledgments This work was supported by the Center for Modelling and Simulation in the Biosciences(BIOMS) and the Cluster of Excellence Cellular Networks at Heidelberg. S.K. was supportedby a fellowship from Deutsche Forschungsgemeinschaft (Grant No. KL818/1-1 and 1-2)and by the NSF through the Center for Theoretical Biological Physics (Grant No. PHY-0822283).
APPENDIX A: ADHESIVE MOTOR DYNAMICS
The Langevin equation, Eq. (3), and the motor dynamics rules explained in Sec. II B areconnected by the following algorithmic rules that apply in each time step ∆ t :(i) The sphere’s position and orientation is updated according to Eq. (3) (for an explicitdescription see Ref. [26]).(ii) The positions where the motors are attached to the sphere in the flow chamber coor-dinate system are calculated.(iii) For each motor that is bound to the MT its load force is calculated. Then steppingis checked according to the stepping rate derived from Eq. (10). If the motor stepsforward the load force is again calculated as motor length/direction has changed.(iv) For each motor that is not bound to the MT binding is checked according to theprocedure explained in the main text (Sec. II B).(v) For each active motor (i. e., bound to the MT), the contribution to F D is calculated.(vi) Each motor that is bound to the MT can unbind with escape rate ǫ given by the Bellequation, Eq. (8).A motor that escaped from the MT in one time step can rebind to the MT in the next timestep according to rule (iv). The same Monte-Carlo technique that is explained in the maintext (Sec. II B) to decide whether binding occurs or not is also used for the decission on step-ping and unbinding. For the pseudo random number generator we used an implementationof the Mersenne Twister [50]. 36 PPENDIX B: MOTOR DISTRIBUTION ALGORITHM
Initially the center of the sphere is located at position (0 , , R + l + h MT ) directly above amicrotubule binding site (cf. Fig. 1). The first motor that is distributed is initially fixed atposition (0 , , l + h MT ) (relative to the chamber coordinate system) with its tail. The headis bound to the microtubule at (0 , , h MT ). Thus the initial distance between the motor andthe microtubule is given by the motor resting length l . For the distribution of the other N tot − r from the uniform distribution on (0 , π ) and r from the uniform distribution on (0 , r , arccos(1 − r )) and possible overlap to already distributed motors is checked. If noother motor is located within a ball of radius 0 . l around the just distributed motor, thenits position is kept, otherwise a new pair of random coordinates are drawn until no overlapwith other motors exists.
1. J. Howard.
Mechanics of motor proteins and the cytoskeleton . Sinauer Associates, Sunderland,Massachusetts, 2001.2. D. L. Coy, M. Wagenbach, and J. Howard. Kinesin Takes One 8-nm Step for Each ATP ThatIt Hydrolyzes.
The Journal of Biological Chemistry , 274(6):3667–3671, 1999.3. K. Visscher, M. J. Schnitzer, and S. M. Block. Single kinesin molecules studied with a molecularforce clamp.
Nature , 400:184–189, 1999.4. S. P. Gross, M. Vershinin, and G. T. Shubeita. Cargo transport: two motors are sometimesbetter than one.
Current biology , 17(12):R478–R486, June 2007.5. G. Koster, M. van Duijn, B. Hofs, and M. Dogterom. Membrane tube formation from giantvesicles by dynamic association of motor proteins.
PNAS , 100:15583–15588, 2003.6. C. Leduc, O. Campas, K. B. Zeldovich, A. Roux, P. Jolimaitre, L. Bourel-Bonnet, B. Goud,J.-F. Joanny, P. Bassereau, and J. Prost. Cooperative extraction of membrane nanotubes bymolecular motors.
PNAS , 101:17096–17101, 2004.7. S. Klumpp and R. Lipowsky. Cooperative cargo transport by several molecular motors.
PNAS , Annual review of neuroscience , 23:39–71, 2000.9. D. B. Hill, M. J. Plaza, K. Bonin, and G. Holzwarth. Fast vesicle transport in PC12 neurites:velocities and forces.
European Biophysics Journal , 33:623–632, 2004.10. S. P. Gross. Hither and yon: A review of bi-directional microtubule-based transport.
Phys.Biol. , 1:R1–R11, 2004.11. Melanie J. I. M¨uller, Stefan Klumpp, and Reinhard Lipowsky. Tug-of-war as a cooperativemechanism for bidirectional cargo transport by molecular motors.
Proceedings of the NationalAcademy of Sciences , 105(12):4609–4614, 2008.12. S. M. Block, L. S. Goldstein, and B. J. Schnapp. Bead movement by single kinesin moleculesstudied with optical tweezers.
Nature , 348(6299):348–352, November 1990.13. K. J. B¨ohm, R. Stracke, P. M¨uhlig, and E. Unger. Motor protein-driven unidirectional transportof micrometer-sized cargoes across isopolar microtubule arrays.
Nanotechnology , 12:238–244,2001.14. J. Beeg, S. Klumpp, R. Dimova, R. Serral Graci`a, E. Unger, and R. Lipowsky. Transport ofbeads by several kinesin motors.
Biophysical Journal , 94:532–541, 2008.15. Michael Vershinin, Brian C. Carter, David S. Razafsky, Stephen J. King, and Steven P. Gross.Multiple-motor based transport and its regulation by tau.
Proceedings of the National Academyof Sciences , 104(1):87–92, January 2007.16. A.R. Rogers, J.W. Driver, P.E. Constantinou, D.K. Jamison, and M.R. Diehl. Negative inter-ference dominates collective transport of kinesin motors in the absence of load.
Phys. Chem.Chem. Phys. , 11:4882–4889, 2009.17. A. Kunwar, M. Vershinin, J. Xu, and S. P. Gross. Stepping, strain gating, and an unexpectedforce-velocity curve for multiple-motor-based transport.
Curr. Biol. , 18:1173–1183, 2008.18. D. A Lauffenburger and J. J. Linderman.
Receptors . Oxford University Press, 1993.19. T. Erdmann and U. S. Schwarz. Stability of adhesion clusters under constant force.
Phys. Rev.Lett. , 92:108102, 2004.20. T. Erdmann and U. S. Schwarz. Stochastic dynamics of adhesion clusters under shared constantforce and with rebinding.
Journal of Chemical Physics , 121(18), 2004.21. T. A. Springer. Traffic signals for lymphocyte recirculation and leukocyte emigration: the ultistep paradigm. Cell , 76:301–314, 1994.22. D. A. Hammer and S. M. Apte. Simulation of cell rolling and adhesion on surfaces in shear flow:general results and analysis of selectin-mediated neutrophil adhesion.
Biophysical Journal ,63:35–57, 1992.23. C. B. Korn and U. S. Schwarz. Dynamic states of cells adhering in shear flow: From slippingto rolling.
Physical Review E , 77(4):041904, 2008.24. M. Schliwa, editor.
Molecular Motors . Wiley-VCH, 2003.25. L. Limberis, J. J. Magda, and R. J. Stewart. Polarized Alignment and Surface Immobilizationof Microtubules for Kinesin-Powered Nanodevices.
Nano Letters , 1(5):277–280, 2001.26. C. B. Korn and U. S. Schwarz. Mean first passage times for bond formation for a Brownianparticle in linear shear flow above a wall.
Journal of Chemical Physics , 126(9):095103, 2007.27. D. L. Ermak and J. A. McCammon. Brownian dynamics with hydrodynamic interactions.
Journal of Chemical Physics , 69(4):1352–1360, 1978.28. J. F. Brady and G. Bossis. Stokesian Dynamics.
Ann. Rev. Fluid Mech. , 20:111–157, 1988.29. B. Cichocki and R. B. Jones. Image representation of a spherical particle near a hard wall.
Physica A , 258:273–302, 1998.30. W. Hortshemke and R. Lefever.
Noise-Induced Transitions . Springer-Verlag, 1984.31. G. S. Perkins and R. B. Jones. Hydrodynamic interaction of a spherical-particle with a planarboundary. II. Hard-wall.
Physica A , 189:447–477, 1992.32. K. Svoboda, C. F. Schmidt, B. J. Schnapp, and S. M. Block. Direct observation of kinesinstepping by optical trapping interferometry.
Nature , 365:721–727, 1993.33. F. Gibbons, J.-F. Chauwin, M. Desp´osito, and J. V. Jos´e. A Dynamical Model of Kinesin-Microtubule Motility Assays.
Biophysical Journal , 80:2515–2526, 2001.34. T. Duke and S. Leibler. Motor Protein Mechanics: A Stochastic Model with MinimalMechanochemical Coupling.
Biophysical Journal , 71:1235–1247, 1996.35. M. J. Schnitzer, K. Visscher, and S. M. Block. Force production by single kinesin motors.
Nature Cell Biology , 2:718–723, 2000.36. A. Rohrbach, E.-L. Florin, and E. H. K. Stelzer. Photonic Force Microscopy: Simulation ofprinciples and applications.
Proceedings of SPIE , 4431:75–86, 2001.37. J. Kerssemakers, J. Howard, H. Hess, and S. Diez. The distance that kinesin-1 holds itscargo from the microtubule surface measured by fluorescence interference contrast microscopy. NAS , 103(43):15812–15817, 2006.38. G. I. Bell. Models for the Specific Adhesion of Cells to Cells.
Science , 200:618–627, 1978.39. K. Svoboda and S. M. Block. Force and velocity measured for single kinesin molecules.
Cell ,77(5):773–784, 1994.40. N. J. Carter and R. A. Cross. Mechanics of the kinesin step.
Nature , 435(7040):308–312, 2005.41. S. M. Block, C. L. Asbury, J. W. Shaevitz, and M. J. Lang. Probing the kinesin reaction cyclewith a 2D optical force clamp.
PNAS , 100(5):2351–2356, 2003.42. R. Lipowsky, S. Klumpp, and T. M. Nieuwenhuizen. Random Walks of Cytoskeletal Motorsin Open and Closed Compartments.
Physical Review Letters , 87(10):108101, 2001.43. S. Chen and T. A. Springer. Selectin receptor-ligand bonds: Formation limited by shear anddissociation governed by the Bell model.
PNAS , 98(3):950–955, 2001.44. C. Korn.
Stochastic dynamics of cell adhesion in hydrodynamic flow . PhD thesis, PotsdamUniversity, 2007.45. C. L. Asbury, A. N. Fehr, and S. M. Block. Kinesin Moves by an Asymmetric Hand-Over-HandMechanism.
Science , 302:2130–2134, 2003.46. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery.
Numerical Recipes in C .Cambridge University Press, second edition, 1994.47. C. W. Gardiner.
Handbook of Stochastic Methods . Springer-Verlag, 1985.48. A. J. Hunt, F. Gittes, and J. Howard. The force exerted by a single kinesin molecule againsta viscous load.
Biophys J , 67(2):766–781, August 1994.49. V. Levi, A. S. Serpinskaya, E. Gratton, and V. Gelfand. Organelle transport along microtubulesin xenopus melanophores: evidence for cooperation between multiple motors.
Biophysicaljournal , 90(1):318–327, January 2006.50. M. Matsumoto and T. Nishimura. Mersenne twister: A 623-dimensionally equidistributeduniform pseudorandom number generator.
ACM Trans. on Modeling and Computer Simulation ,8(1):3–30, 1998.,8(1):3–30, 1998.