Tug-of-war between two elastically coupled molecular motors: a case study on force generation and force balance

Intracellular transport is performed by molecular motors that pull cargos along cytoskeletal filaments. Many cellular cargos are observed to move bidirectionally, with fast transport in both directions. This behaviour can be understood as a stochastic tug-of-war between two teams of antagonistic motors. The first theoretical model for such a tug-of-war, the Mu¨ller–Klumpp–Lipowsky (MKL) model, was based on two simplifying assumptions: (i) both motor teams move with the same velocity in the direction of the stronger team, and (ii) this velocity matching and the associated force balance arise immediately after the rebinding of an unbound motor to the filament. In this study, we extend the MKL model by including an elastic coupling between the antagonistic motors, and by allowing the motors to perform discrete motor steps. Each motor step changes the elastic interaction forces experienced by the motors. In order to elucidate the basic concepts of force balance and force fluctuations, we focus on the simplest case of two antagonistic motors, one kinesin against one dynein. We calculate the probability distribution for the spatial separation of the motors and the dependence of this distribution on the motors’ unbinding rate. We also compute the probability distribution for the elastic interaction forces experienced by the motors, which determines the average elastic force h F i and the standard deviation of the force fluctuations around this average value. The average force h F i is found to decrease monotonically with increasing unbinding rate e 0 . The behaviour of the MKL model is recovered in the limit of small e 0 . In the opposite limit of large e 0 , h F i is found to decay to zero as 1/ e 0 . Finally, we study the limiting case with e 0 = 0 for which we determine both the force statistics and the time needed to attain the steady state. Our theoretical predictions are accessible to experimental studies of in vitro systems consisting of two antagonistic motors attached to a synthetic scaffold or crosslinked via DNA hybridization.


Introduction
Intracellular cargos such as vesicles and organelles are transported by cytoskeletal motors. 1 Conventional kinesin and cytoplasmic dynein represent two types of cytoskeletal motors that walk along microtubules in opposite directions. 2,3 Many cargos are observed to perform a bidirectional movement on the microtubules with frequent reversals. 4,5 This behaviour reflects the presence of two antagonistic motors, plus-end directed kinesin and minus-end directed dynein, on the same cargo. These motors try to pull the cargo in their preferred direction of motion, thereby performing a stochastic tug-of-war. The first theoretical model for such a stochastic tug-of-war was introduced by Müller, Klumpp, and Lipowsky (MKL) [6][7][8] and corroborated by the observations of endosome transport in amoebae 9 and fungi. 10 The MKL model was based on two simplifying assumptions: (i) all motors move with the same velocity in the direction of the stronger team and (ii) this matching of the velocities sets in as soon as the cargo is pulled by motors from both teams. Thus, for two antagonistic motors, the MKL tug-of-war state is characterized by the following properties. When both motors are simultaneously bound to the filament, they experience mutual interaction forces which are of equal magnitude and opposite direction, in accordance with Newton's third law. The absolute value of this interaction force, the so-called cargo force F ca , is determined uniquely by the characteristic force-velocity relations of the motors and the condition of velocity-matching under this force. 6,8 The cargo then moves with this generally low, but nonzero velocity v ca in the direction of the stronger motor. In the special case of equally strong motors the cargo is in a stalled state with zero velocity.
In the MKL model, the motion of the motors is described in a coarse-grained manner, averaging over the discrete steps of the individual motors. Here, we extend this model by including these discrete steps, following the theoretical approach developed in ref. 11 and 12 for two identical motors that pull in the same direction. We consider two antagonistic motors which are coupled to their common cargo via two elastic linkers. We then use the force balance between the motors and the cargo to derive an effective elastic coupling between the two motors. With each step taken by one of the motors, the effective linker is either stretched or compressed and an elastic force is induced acting on both motors. The velocity-matching condition can therefore only be reached, in general, after the motors have taken multiple steps.
When the tug-of-war involves antagonistic motor teams consisting of several identical motors, the MKL model makes the additional assumption that the overall load acting on one team is shared equally by all motors in that team. This assumption has been previously criticized to represent a mean-field approximation because it ignores fluctuations in the load sharing. [13][14][15] In order to examine the latter fluctuations, the authors in ref. [13][14][15] studied such antagonistic motor teams. In contrast to these previous studies, we focus here on the simplest case of two antagonistic motors for which load sharing does not play any role. On the one hand, the case of 1 + 1 motors is useful in order to elucidate the basic concepts of force balance and force fluctuations. In fact, as shown below, this case already implies a nontrivial force balance, with an average interaction force that depends strongly on the unbinding rate e 0 of the individual motors and decays to zero for large e 0 . On the other hand, the study of the 1 + 1 motor system allows us to perform a detailed comparison between (i) the tug-of-war of two elastically coupled motors that perform discrete steps and (ii) the tug-of-war of two antagonistic motors as described by the MKL model.
The theoretical results described below can be scrutinized by experimental in vitro studies based on recently developed protocols [16][17][18][19][20] to control the number of active motors on synthetic molecular scaffolds. Evidence for a tug-of-war mechanism between kinesin and dynein attached to such scaffolds was observed in ref. 19 and 20. Very recently, it has also been demonstrated that one can directly crosslink a single, fluo-labeled kinesin with a single, fluo-labeled dynein via DNA hybridization. 21 For such two-motor constructs, one should be able to measure the probability distribution for the spatial separation of the two motors along the filament and, in this way, directly scrutinize the predictions of our theory. Furthermore, as shown below, the distribution for the motor-motor separation also determines the probability distribution for the elastic interaction forces and, thus, the average elastic force hFi between the motors. This paper is organized as follows. First, we briefly review (i) the single motor description in Section 2.1, thereby introducing the single motor parameters used in our model, and (ii) the cargo force predicted by the MKL model in Section 2.2. Next, the force balance for a tug-of-war between two elastically coupled motors is considered and the associated state space for this process is defined in Section 2.3. We describe the system as a Markov process with (i) several transient states corresponding to different extensions of the effective elastic linker and (ii) two absorbing states defined by a single plus-or minus-end motor bound to the filament. The rates of the network are determined using single motor parameters characterising the stepping and unbinding behaviour. Section 3 reports the results on the steady state probability distributions of 2motor runs with both motors attached to the filament, the average force experienced by the motors, as well as a detailed study of the limiting case of zero unbinding rates. In the latter case, both the force statistics and the time needed to reach the steady state are determined as a function of the elastic coupling strength.

Single motor description
When a motor binds to the filament, it steps along the filament in a preferred direction. In the absence of an external force F, the motor moves with its zero-force forward velocity towards the preferred end of the filament. We use the convention that a resisting force acting as a load on the motor has a positive sign whereas an assisting force pulling the motor in its preferred stepping direction has a negative sign. This convention is used both for plus-directed and for minus-directed motors. With increasing load force, the motor velocity decreases until the force reaches the motor's stall force F s at which the motor velocity vanishes. Experimental studies provide strong evidence that the kinesin-1 motor steps backwards for load forces F 4 F s . 22,23 The explicit form of the force-velocity relation for a single motor can be obtained from fits to experimental data or provided by piecewise linear relations as in previous theoretical studies, see e.g. ref. 6, 13 and 15 and Appendix A. Here, we will use a convenient parametrization of the force-velocity relation as introduced in ref. 11 which has the form The parameters v max and v min determine the limits of v(F) for large negative and large positive values of F, respectively. The zero-force velocity v(F = 0) is given by the parameter v 0 . In the following, we use the force-velocity relation in eqn (2.1) for both motors. The parameters v 0 , v max , v min and F s can be specified in order to obtain close approximations for experimentally determined force-velocity relations, e.g. as given in ref. 22 for kinesin-1 and in ref. 24 for cytoplasmic dynein.
Another single motor property that has been measured for kinesin as a function of load force 22,23 is the forward-tobackward stepping ratio q which is defined as the number of forward steps divided by the number of backward steps as observed within a certain time interval. This ratio was found to depend exponentially on the load force and to be well fitted by 22 qðFÞ ¼ q 1ÀF=Fs 0 with q 0 ¼ 800 and F s ¼ 7 pN for kinesin: For dynein, the corresponding parameters have not been measured directly but the single motor data on yeast dynein in ref. 25  We now use the force-velocity relationship v(F) and the forward-to-backward stepping ratio q(F) together with the step size l to define the forward stepping rate of a single motor by aðFÞ vðFÞ ' qðFÞ qðFÞ À 1 (2.4) and its backward stepping rate by bðFÞ vðFÞ ' which implies a/b = q and a À b = lv. 11 The force-dependence of the two stepping rates a and b is illustrated in Fig. 1 for yeast dynein. At the stall force F = F s , the two stepping rates are equal, implying that forward and backward steps are equally likely. For F 4 F s , backward steps are more likely than forward steps and both stepping rates decay monotonically with increasing force. For F o F s , on the other hand, forward steps are more likely than backward steps but the individual stepping rates are nonmonotonic as a function of F. In fact, the backward stepping rate b for dynein exhibits a pronounced maximum at the load force F = 3.21 pN, arising from the relatively small value q 0 = 4 of the forward-to-backward stepping ratio. For kinesin, on the other hand, which is characterized by a much larger value of q 0 , no such maximum of b is found. A motor bound to a filament unbinds from this filament with a constant unbinding rate in the absence of external forces. We denote this zero-force unbinding rate of a single motor by e 0 . When a force acts on the cargo, the motor-filament bond is more likely to break. Although the unbinding process is very complex on a molecular scale, it can be approximately described as an escape process of a particle in a potential well. According to Kramers' theory, 27 the force-dependence of the motor's unbinding rate is then approximately exponential and given by e(F) = e 0 exp(|F|/F d ), (2.6) where F d is the detachment force, another force scale that characterizes each motor type. Here and below, the force F represents the tangential force component acting parallel to the long axis of the filament, with the previously mentioned sign convention that F is positive when it acts against the preferred stepping direction that the motor has in the absence of force. Note that e(F) is taken to depend only on the absolute value |F|, and not on the direction of the force which implies that resisting and assisting forces increase the unbinding rate by the same amount. This simplifying assumption is not crucial, however, because, as we will see below, assisting forces are almost never generated by a tug-of-war.
In the following, we use this single motor description for two different types of motors, kinesin and dynein, which represent the best studied examples for processive plus-end and minusend directed motors. We will use the notation F AE s and F AE d for the stall and detachment forces of these two motor species. Our sign convention for F implies that both stall and detachment forces are positive. In addition, this convention also implies that a positive force acting on the plus-end directed motor points towards the minus end of the filament whereas a positive force acting on the minus-end directed motor points towards the filament's plus end. Furthermore, because of the opposite directionality, the force-velocity relationships v + (F) and v À (F) for the plus-and minus-directed motors have the form v + (F) = +v(F) and v À (F) = Àv(F) with v(F) as in eqn (2.1).

Cargo force in the MKL model
Before we consider the tug-of-war between elastically coupled motors, we will first summarize the main properties of the tugof-war in the MKL model which provides a useful reference process. As described in Appendix A, the latter process is characterized by instantaneous velocity matching between the different motors. The corresponding matching condition can be visualized by plotting the two force-velocity relations for the individual motors in the same (F,v)-diagram. 8 The intersection point of these two relations provides the matching condition for the MKL tug-of-war as illustrated in Fig. 2 for the case in   2 Matching condition for the velocities of two antagonistic motors as used in the MKL tug-of-war model: force-velocity relations as given by eqn (2.1) for a plus (green) and a minus (red) motor with zero-force velocities v + 0 4 0 and v À 0 o 0. The intersection point (F, v) = (F ca , v ca ) of the two force-velocity relations defines the velocity-matched state in which both motors move with the same velocity v ca and experience the same single motor force as provided by the cargo force F = F ca . In this example, the stall force F + s of the plus motor exceeds the stall force F À s of the minus motor which implies v ca 4 0, i.e. both the cargo and the two antagonistic motors move towards the plus end of the filament. which the stall force F + s of the plus motor exceeds the stall force F À s of the minus motor. In general, the intersection point (F, v) = (F ca , v ca ) of the two force-velocity relations defines the velocity-matched state in which the cargo and the motors move with the same velocity v ca and the two motors experience the same single motor force, namely the cargo force F = F ca . The cargo force is always positive because it is located between the stall forces F + s and F À s of the two motors, both of which are positive by definition. Thus, according to our sign convention for single motor forces, the cargo force acts as a resisting force on both motors: it represents the absolute value of the two opposing forces that the cargo exerts simultaneously onto the plus and onto the minus motor.

Elastic coupling between the motors
We now consider one kinesin and one dynein motor pulling on the same cargo particle via elastic linkers as illustrated in Fig. 3. We use the Cartesian coordinate x parallel to the filament to describe the positions of the motors and the cargo along this filament. The coordinate x is chosen to increase towards the plus end of the filament. Thus, each motor-cargo configuration is described by the positions x ki , x dy , and x ca with x dy o x ca o x ki . Furthermore, to discuss the elastic forces acting between the motors and the cargo, we will first define these forces with respect to the coordinate x which implies that we temporarily use a different sign convention for these interaction forces compared to the single motor forces. Thus, in the following paragraph, the elastic interaction forces are taken to be positive when they point towards the plus end of the filament and negative when they point towards the filament's negative end.
Elastic forces between the cargo and the motors. The linkers between the motors and the cargo are described by harmonic springs with spring constant k and rest length L J . The kinesin motor then exerts the force onto the cargo. Likewise, the dynein motor exerts the force F dy,ca = Àk(x ca À x dy À L J ) (2.8) onto the cargo. We now assume that, for given positions x ki and x dy , the elastic forces balance each other on timescales that are short compared to the timescales of the single motor transitions. This elastic force balance implies F ki,ca + F dy,ca = 0 and x ca = 1 2 (x ki + x dy ). Eliminating the cargo position x ca from the expressions in eqn (2.7) and (2.8), we obtain the forces F ki,ca = 1 2 k(x ki À x dy À 2L J ) = K(x ki À x dy À L 0 ) (2.9) and F dy,ca = 1 2 k(x ki À x dy À 2L J ) = ÀK(x ki À x dy À L 0 ), (2.10) which depend only (i) on the coordinate difference x ki À x dy of the two motor positions and (ii) on the effective spring parameters K k/2 and L 0 = 2L J . (2.11) Introducing the combined spring extension DL x ki À x dy À L 0 , (2.12) the force that the kinesin exerts on the dynein becomes F ki,dy = F ki,ca = KDL, (2.13) while the force that the dynein exerts on the kinesin has the form F dy,ki = ÀKDL, (2.14) as required by Newton's third law. Identification with single motor forces. In order to use the single motor description as described in the previous subsection, we now return to our original sign convention for the force F acting on a single motor. As a consequence, the single motor force is given by F = F ki ÀF dy,ki for kinesin, (2.15) which is positive when the force F dy,ki points towards the minus end of the filament, and by F = F dy F ki,dy for dynein, (2.16) which is positive when the force F ki,dy points towards the filament's plus end. Newton's third law as given by F ki,dy = ÀF dy,ki then assumes the simple form F ki = F dy = F. State space for tug-of-war with elastic coupling. Kinesin and dynein have the same step size l C 8 nm, see references in Table 1. We further assume that the motor pair can attain a Fig. 3 Different states ( j) with j = 0, 1, m and n of two antagonistic motors corresponding to different extensions of the elastic linkers between the motors and the cargo. The kinesin motor (blue ''heads'') and the dynein motor (green ''wheels'') prefer to move towards the (+)-and (À)-end of the microtubule, respectively. In state (0), the motor linkers are relaxed and do not generate elastic forces. When the motors perform steps leading to a state ( j) with j 4 0, the spring becomes stretched and generates the elastic force F j = jF K . As explained in the main text, the single motor forces F ki and F dy acting on kinesin and dynein are defined in such a way that Newton's third law assumes the simple form F j = F ki = F dy . In state (1), for example, both the kinesin and the dynein motor experience the single motor force F 1 = F K . In state (m), the force F m = mF K C F À s , i.e. it is comparable to the stall force F À s of the minus motor. At this stall force, the minus motor steps forward and backward with equal probability. In state (n), the force F n is close to the stall force of the plus motor, which can now step forward and backward with (almost) equal probability.
This journal is © The Royal Society of Chemistry 2017 relaxed state with DL = 0. It is then convenient to introduce the dimensionless spring extension j = DL/l with ÀJ r j r J, (2.17) where j = J represents the maximal stretching and j = ÀJ the maximal compression of the spring. Each motor behaves as a stochastic stepper with forcedependent forward and backward stepping rates a(F) and b(F), as given by eqn (2.4) and (2.5). Following the approach in ref. 11 for two identical motors, we now introduce a discrete state space with states ( j) labeled by an integer j with ÀJ r j r J. Each state ( j) corresponds to a certain extension of the elastic spring. We consider two identical motor linkers in series which implies the effective spring constant K k/2 for the elastic coupling between the two motors as in eqn (2.11). When the spring extension is increased by a single motor step with step size l, the elastic force experienced by both motors is increased by the strain force Fig. 3. The elastic force acting between the motors is then given by When one of the motors performs a step, the elastic force acting between the two motors changes monotonically from its initial value before the step to its final value after the step. As a consequence, the effective force acting during the stretching transition from ( j) to ( j + 1) is given by i.e. by the arithmetic mean of the forces acting before and after such a step. 28 Likewise, the effective force acting during the compression transition from ( j) to ( j À 1) has the form Motor kinetics and transition rates. The dynamics of the two-motor runs can be investigated by considering the stochastic process on the discrete state space corresponding to the ''stretching'' or ''compression'' of the effective spring between the motors. In state (0) the spring is relaxed. If one of the motors performs a single step, the elastic spring of the twomotor system can be stretched or compressed by l/2 and the system undergoes a transition to states (1) or (À1), respectively, see the network representation in Fig. 4. In general, the motor system undergoes a forward transition from state ( j) to state ( j + 1) if kinesin or dynein performs a forward step. The corresponding stepping rates are given by a AE ( % F j4 ) for the individual motors which implies the forward rate with ÀJ r j r J and the boundary condition o f (J) = 0, see Fig. 4. Each forward transition leads to an increased stretching or a reduced compression of the elastic coupling between the two motors. Likewise, the motor system undergoes a backward transition from state ( j) to state ( j À 1) if kinesin or dynein performs a Table 1 Values of the parameters used for kinesin-1 and two types of dynein motors: the values for ''strong'' and ''weak'' dynein correspond to yeast and mammalian cells, respectively. A star superscript indicates a parameter for which we did not find experimental data in the literature; the corresponding parameter value was set equal to the experimentally deduced value of another type of motor. For the minimal and maximal velocities of both dynein motors in the force-velocity relationship (2.1) we used the estimated values v min C 0.12v 0 and v max C 1.12v 0 , as indicated by the † symbol. For the parameters in this table the strain force is given by F K = kl/2 = 0. 8 Fig. 4 State space associated with different states of the elastic linker. The two states (1, +) and (1, À) are the absorbing states with a single plus and minus motor bound to the microtubule. States labelled by integer j with ÀJ r j r J denote the states with a stretched and compressed linker for j 4 0 and j o 0, respectively. Starting from state ( j), the plus and minus motors unbind from the filament with rates o + off ( j) and o À off ( j). Furthermore, the motors undergo a forward transition from ( j) to ( j + 1) with rate o f ( j) and a backward transition from ( j) to ( j À 1) with rate o b ( j).
with ÀJ r j r J and the boundary condition o b (ÀJ) = 0, see Fig. 4. Each backward transition leads to a reduced stretching or an increased compression of the elastic coupling between the two motors.
In state ( j), each of the motors can unbind (or detach) from the filament. The corresponding unbinding rates have the form for the two types of motors. The overall unbinding rate from state ( j) is then given by Motor parameters. In the following, we will study the tug-of-war between one plus-directed motor and one minus-directed motor. Most single motor parameters will be kept fixed and assume values as appropriate for kinesin-1 and strong dynein, see Table 1. One parameter that we will vary systematically is the zero-force unbinding rate As we will see, the elastically coupled tug-of-war is characterized, in the limit of small e 0 , by essentially the same force balance as the MKL tug-of-war. Another parameter of the motor system that will be varied systematically is the elastic coupling between the two motors as described by the strain force F K = Kl = kl/2. In order to ensure that the finite size of the state space does not affect our results, we will always use a sufficiently large state space with (2J + 1) Z 101.

Steady state properties of tug-of-war
The elastic coupling between the two antagonistic motors is only effective as long as both motors are attached to the filament and perform a 2-motor run. The latter runs are described by transitions between the states ( j) in Fig. 4 and are terminated as soon as one of the motors unbinds from the filament. After such an unbinding event, the cargo is bound to the filament by a single motor as described by the states (1, +) and (1, À) in Fig. 4. These 1-motor runs continue until the unbound motor rebinds to the filament. We will assume that the rebinding typically leads to the state (0) which is relaxed in the sense that the two motors do not experience elastic interaction forces. Thus, after rebinding, a new 2-motor run starts from the initial state (0). As described previously for the 2-motor runs of two identical motors, 11,12 the steady state probability distribution as obtained from an ensemble average over many 2-motor runs can be replaced by a time average over a concatenated 2-motor run that is obtained by redirecting all transitions to the absorbing states towards a certain initial state. For the network depicted in Fig. 4, the absorbing states are provided by (1, +) and (1, À) and the initial state is taken to be the relaxed state (0). As a result, we obtain the redefined network in Fig. 12 which has no absorbing states. The corresponding steady state probability distribution can then be calculated by solving the master equation for the redefined network, see Appendix B.
The steady state probability distribution p st j describes the frequencies with which the effective elastic spring between the two motors has the extension j. The latter extension determines the spatial separation of the two motors along the filament. The average motor-motor separation is then given by The fluctuating motor-motor separation L should be directly accessible to experimental studies when one combines the recently introduced crosslinking of one fluo-labeled kinesin and one fluo-labeled dynein via DNA hybridization 21 with advanced methods of fluorescence imaging such as FIONA. 29 In Fig. 5(a) we plot p st j for different values of the zero-force unbinding rate e 0 . We see that for low unbinding rates, e.g. e 0 = 0.01 s À1 , the probability distribution p st j shifts towards states with larger j-values. Because the elastic force F j corresponding to spring extension j is given by F j = jF K , we can transform the occupation probabilities of the states into the corresponding force distribution by a change of variables from j to F j , see Fig. 5(b). As shown in the latter figure, a decrease in the unbinding rate leads to a shift of the force distribution towards higher force values and to an average elastic force that approaches the cargo force F ca = 7 pN as obtained for velocity-matching. For higher unbinding rates, the occupation probabilities are shifted towards lower j-values, leading to a reduced motor-motor separation and indicating that the motors are likely to unbind from the filament before reaching a state with velocity matching. The relationship in (3.3) implies that the average elastic force hFi can be determined from the average motor-motor separation hLi.

Dependence of average elastic force on the unbinding rate
For two antagonistic motors coupled by an effective spring as studied here, the elastic interaction forces fluctuate and lead to steady state force distributions as shown in Fig. 5(b) for different values of the unbinding rate e 0 . In contrast, for the MKL model, the mutual interaction force is given by the constant cargo force F ca as obtained via velocity matching, see Fig. 2. Inspection of Fig. 5 we now consider the very low unbinding rate e 0 = 10 À5 s À1 and calculate the corresponding steady state probability distribution p st j . This distribution and the associated force distribution are displayed in Fig. 6 for different values of the strain force F K . Both for kinesin vs. strong dynein as depicted in Fig. 6(a1 and a2), and for kinesin vs. weak dynein in Fig. 6(b1 and b2), the cargo force F ca provides a more accurate approximation to the average elastic force hFi for smaller values of the strain force F K . This behaviour arises because smaller F K -values imply smaller changes in the elastic force induced by single motor  Table 1 for kinesin and strong dynein, which implies the strain force F K = 0.8 pN; and (b) steady state distributions for the elastic forces as obtained from p st j by a change of variables from j to F j = jF K . The dashed vertical line (red) represents the cargo force F ca of the velocity-matched model. For low unbinding rates e 0 , the average elastic force approaches this cargo force. For larger values of e 0 , the force distribution becomes broader and shifts towards lower force values. As a consequence, the average elastic force hFi becomes smaller than the cargo force F ca , see average force values in the inset. Likewise, the average spatial separation L 0 + hjil of the two motors decreases with increasing unbinding rate e 0 as follows from the distributions p st j in (a). Fig. 6 Steady state probability distributions and force distributions for the small unbinding rate e 0 = 10 À5 s À1 and for different choices of the strain force F K : (a1 and a2) kinesin against strong dynein, see the motor parameters in Table 1. As we decrease F K , the average force hFi approaches the cargo force F ca = 7 pN more accurately; and (b1 and b2) kinesin against weak dynein, see again Table 1. The cargo force now has the lower value F ca = 4.55 pN compared to the kinesin vs. strong dynein case in (a). Accordingly, the state (0) in (b1) has a higher occupation probability than the same state in (a1), and the probability for the force value F 0 = 0 is increased in (b2) compared to (a2).
steps which lead to a larger number of accessible force values and, thus, to a smoother interpolation of the discrete force distributions. In addition, the smaller cargo force for the tug-ofwar between kinesin and weak dynein compared to kinesin and strong dynein implies that the force value F 0 = 0 has a higher probability for the weak dynein case, see Fig. 6. As shown in Fig. 5, increasing the unbinding rate e 0 leads to a reduction in the average elastic force hFi induced by the motors. This unbinding rate dependence of the average elastic force is displayed in more detail in Fig. 7(a) for different choices of the strain force F K , with Fig. 7(b) magnifying the limit of small e 0 . In this limit, the deviation of the average elastic force hFi from the cargo force F ca increases for larger values of F K . The decrease of the average force hFi with increasing unbinding rate e 0 is caused by the increasing probability that one of the motors unbinds from the filament before the motors have matched their velocities in the 2-motor run and, thus, before the motors can generate forces comparable to the cargo force F ca . The double-logarithmic plot in Fig. 7(c) reveals that the average elastic force hFi decays to zero for large e 0 . It follows from eqn (3.3) that the asymptotic behavior of hFi for large e 0 is determined by the asymptotic behavior of the steady state probability distribution p st j for large e 0 . The latter behavior can be directly obtained from the master equations (B.1)-(B.4). One then finds from the local flux balance in the states ( j) that The power law behavior hFi B 1/e 0 is corroborated by the numerical data in Fig. 7(c) which are well fitted, over the accessible range of e 0 -values as given by 1 s À1 r e 0 r 10 6 s À1 , by a power-law of the form hFi B 1/e z with the effective decay exponent z. As shown in Fig. 7(d), the effective exponent z is found to depend weakly on the strain force F K for F K o 13 pN and to approach the true asymptotic value z = 1 for F K 4 13 pN.

Statistics of elastic forces for vanishing unbinding rate
We now look at the properties of the tug-of-war in the limit in which the two motors can no longer unbind from the filament, corresponding to zero-force unbinding rate e 0 = 0 which implies that o off ( j) = 0 for all j. In this case, the state space for the tug-of-war between the two motors is reduced to the states j = ÀJ,. . .,+J. For this reduced state space, the probability distribution p j (t), which starts from the initial distribution p j (0) = d j0 at time t = 0, evolves towards a steady state distribution p st j as shown in Fig. 8. The maximum of p st j is located close to the state ( j = 9) characterized by the elastic force 9F K = 7.2 pN induced by the effective spring, while the cargo force obtained from velocity matching is F ca = 7 pN.
The average elastic force hFi for the process with vanishing unbinding rate is shown in Fig. 9 for different values of the strain force F K . We observe that, regardless of the choice of F K , the average force hFi remains close to the cargo force F ca whereas its standard deviation s F increases with increasing F K . In the limit of small strain forces F K , the average force hFi approaches the cargo force F ca more accurately, in accordance with our previous results. As shown in the right inset of Fig. 9, the standard deviation of the elastic force behaves as over the whole range of F K -values considered here. For F K = 1 pN and e 0 = 0, the standard deviation is s F C 1.5 pN. Increasing the zero-force unbinding rate e 0 for fixed F K = 1 pN, the full network in Fig. 4 leads to a slight increase in the deviation s F , which then saturates at s F C 1.9 pN for an unbinding rate of e 0 = 1 s À1 (data not shown here). Using these results for the statistics of the elastic forces, we can directly conclude that the average spring extension h ji behaves as and its standard deviation These dependencies on the strain force F K are displayed in Fig. 10. The average spring extension h ji implies that the two motors have the average separation hDLi = lh ji. Thus, we conclude that, in the limit of small F K corresponding to weakly coupled motors, the average separation between the motors increases . Furthermore, the behavior of the standard deviation s j is consistent with a Gaussian probability distribution of the form for the spring extension j where c is a proportionality factor. It is instructive to compare the distribution as given by (3.14), which arises from the stochastic nature of the tug-of-war and the underlying motor activity, with the equilibrium distribution p eq j p exp[À 1 2 Kj 2 /(k B T)] (3.15) corresponding to thermal fluctuations in the harmonic spring potential 1 2 Kj 2 . Comparing the two distributions in (3.14) and (3.15), we can draw two conclusions. First, the motor activity Fig. 8 Time evolution of the probability distribution p j (t) for unbinding rate e 0 = 0 and strain force F K = 0.8 pN, the latter parameter being appropriate for strong dynein. The initial probability distribution at time t = 0 is given by p j (0) = d j0 corresponding to a relaxed spring with zero extension. As t increases, the distribution p j (t) approaches the steady state distribution p st j as indicated by the dashed red line. Second, the fluctuations around the average value h ji can be characterized by a standard deviation s j $ 1 ffiffiffiffi K p both for the thermal and for the active process.
A comparison between Fig. 7 and 9 also shows that the tug-of-war obtained for the reduced state space with e 0 = 0 is reached in a smooth manner when we consider the full state space as depicted in Fig. 4 and take the limit of small e 0 . This agreement is to be expected because the steady state probability distributions for e 0 4 0 follow from the redefined network in Fig. 12 which also becomes identical with the reduced state space for e 0 = 0.

Relaxation time for vanishing unbinding rate
For the reduced state space consisting of the states ( j) with j = ÀJ,. . .,+J as obtained for unbinding rate e 0 = 0, the elastically coupled motors eventually reach a steady state with the average elastic force hFi C F ca as illustrated in Fig. 8 for F K = F ca /8.75. We will now address the relaxation time for this process, i.e. the time it takes to actually reach this steady state. This relaxation time t re is provided by the largest non-zero eigenvalue l 2 of the transition rate matrix for the redefined state space in Fig. 12 via the relation 30 t re = À1/l 2 . (3.16) As shown in Fig. 11, the relaxation time t re for the approach towards the steady state increases strongly as we reduce the strain force F K , i.e. as we reduce the elastic coupling between the two motors and allow them to move further apart. Indeed, in the limit of small F K , the average separation of the two motors increases as hDLi = lh ji B 1/F K as follows from the behavior of h ji in Fig. 10. When we plot the t re -data in a doublelogarithmic manner, see the inset of Fig. 11, a least-squares fit leads to the relation t re B 1/F Z K with the decay exponent Z = 0.97. We thus conclude that the relaxation time t re is also inversely proportional to F K to a very good approximation. As a consequence, the relaxation time t re for two weakly coupled motors is roughly proportional to their average separation hDLi B 1/F K which diverges in the limit of small F K = lK = lk/2.
In order to obtain a well-defined relaxation time t re , we had to consider the limiting case with unbinding rate e 0 = 0. Real motors have, of course, a finite unbinding rate which implies that their 2-motor runs are terminated after a finite Fig. 10 Average spring extension h ji and standard deviation s j of these extensions as a function of strain force F K for unbinding rate e 0 = 0. Both h ji and s j decrease with increasing F K . The F K -dependence of h ji is very well described by F ca /F K (dashed red line). The standard deviation s j is proportional to 1 ffiffiffiffiffiffi F K p (full blue line in the inset). This behavior is intimately related to the behavior of the average force hFi and the associated standard deviation s F , see the main text. Fig. 11 Strain force dependence of the relaxation time t re (green) needed to attain the steady state for unbinding rate e 0 = 0. The inset displays a double-logarithmic plot of the t re -data. A least-squares fit to these data leads to t R B 1/F Z K with the exponent Z = 0.97 from which we conclude that the relaxation time t re is inversely proportional to F K to a very good approximation. For comparison, the average run time t 2 (red) for 2-motor runs with e 0 = 0.3 s À1 has also been included. For the range of F K -values considered here, the relaxation time t re is always large compared to the average run time t 2 . time and can then be characterized by the average run time . We must now distinguish different cases depending on the relative size of the average run time t 2 and the relaxation time t re . If t 2 is large compared to t re , the two motors will be characterized by an average elastic force hFi close to the cargo force F ca of the MKL model. On the other hand, if t 2 is small compared to t re , the force balance between the two motors is different because hFi is small compared to F ca . As shown in Fig. 11, the latter case applies to the unbinding rate e 0 = 0.3 s À1 for the whole range of F K -values considered in this figure. A more detailed discussion of the different timescales and the associated dynamic regimes will be given in a subsequent publication. 31

Summary and outlook
In this paper, we considered the tug-of-war between one kinesin and one dynein motor, which are coupled to a common cargo via elastic linkers. We started from the known properties of the single motors and used these properties to derive the forcedependent stepping rates for forward and backward steps of both motors. Unexpectedly, we found that the backward stepping rate of strong dynein exhibits a maximum at an intermediate load force, see Fig. 1. We then described the elastic interaction forces between the motors and the cargo by two harmonic springs which can be combined into an effective harmonic spring between the two motors. The extension j of this effective spring was used to define the state space for the tug-of-war between two elastically coupled motors as displayed in Fig. 4. A Markov process was constructed on this state space with transition rates that were derived from the single motor rates. Starting from the relaxed state (0) with spring extension j = 0, the antagonistic motors perform a 2-motor run on the reduced state space consisting of the transient states ( j) with j = ÀJ,. . .,+J until one motor unbinds from the filament and the process ends up in one of the absorbing states (1, +) or (1, À), see Fig. 4. The ensemble average over many such 2-motor runs can be obtained by computing the steady state probability distribution p st j on the redefined state space in Fig. 12. Using this distribution, we calculated the average elastic force hFi experienced by the two motors as a function of the unbinding rate e 0 , see Fig. 7. This average force approaches the cargo force F ca for small e 0 and decays to zero as 1/e 0 for large e 0 , see (3.8). Numerically we find the power law hFi B 1/e z 0 with an effective decay exponent z that depends on the strain force F K , see Fig. 7(d).
Finally, we studied the limiting case of a tug-of-war between two elastically coupled motors that cannot unbind from the filament, corresponding to zero-force unbinding rate e 0 = 0. In this case, we found that the time evolution of the probability distribution always leads to a steady state distribution for which (i) the average elastic force hFi is close to the cargo force F ca irrespective of the strain force F K and (ii) the standard deviation s F of the force fluctuations is proportional to ffiffiffiffiffiffi F K p as shown in Fig. 9. These relationships imply that the average spring extension h ji B 1/F K B 1/K and the standard deviation Fig. 10. The latter dependence is consistent with a Gaussian probability distribution for the spring extension j as described by (3.14). The behavior of the average spring extension h ji implies that the average separation hDLi = lh ji of the two antagonistic motors increases as 1/K for small K, corresponding to the weak coupling limit. Essentially the same K-dependence is found for the relaxation time t re towards the steady state probability distribution p st j , see Fig. 11. Therefore, in the weak coupling limit of small K, both the average motor separation hDLi and the relaxation time t re diverge as 1/K for e 0 = 0.
In the present article, we have focussed on the forces acting between two elastically coupled motors. As indicated in the last subsection 3.4 on the relaxation time, the resulting tug-of-war involves different timescales that define different kinetic regimes. Likewise, as far as motor unbinding is concerned, one has to distinguish spontaneous unbinding for small elastic forces from force-induced unbinding for large elastic forces. Another interesting topic is the tug-of-war between elastically coupled motor teams involving N 1 Z 1 plus-directed motors and N 2 Z 1 minus-directed motors. These more complex motor systems will be addressed in a forthcoming paper. 31 As previously mentioned, the theoretical results presented here can be scrutinized by experimental in vitro studies based on recently developed protocols [16][17][18][19][20] to control the number of active motors on synthetic DNA scaffolds. Evidence for a tug-ofwar mechanism between kinesin and dynein attached to such scaffolds was observed in ref. 19 and 20. Very recently, it has also been demonstrated that one can directly crosslink a single, fluo-labeled kinesin with a single, fluo-labeled dynein via DNA hybridization. 21 Our tug-of-war model provides detailed predictions for the transport properties of such a two-motor system. One of these properties is the probability distribution p st j for the spatial separation of the two motors as displayed in Fig. 5a, 6a1, b1, and 8. The latter distribution should be accessible to advanced methods of fluorescence imaging such as FIONA. 29 In principle, it is also possible to measure the average interaction force hFi directly via FRET-based molecular tension probes 32 that are incorporated into the linkers between the motors and the cargo. Other quantities of interest that can be used to compare our theory with experiment include the run lengths and run times of the two-motor systems.

A Review of tug-of-war with velocity matching
Here we shortly describe the tug-of-war in the MKL model, following ref. 6 and 7. Let N + and N À denote the number of plus and minus motors attached to the cargo, respectively. At any given time t 0 the cargo is pulled by n + plus and n À minus motors, where 0 r n + r N + and 0 r n À r N À . The motility state of the cargo at t 0 is then characterised by the number of bound motors to the microtubule, denoted by (n + , n À ). Assuming that (i) opposing motors act as load and (ii) identical motors share this load 6 the transition rates between adjacent motility states can be inferred from single motor binding and unbinding behaviour. The force balance on a cargo pulled by n + plus and n À minus motors at any moment is given by: F ca (n + ,n À ) n + F + = Àn À F À . (A.1) F + and F À are the forces ''felt'' by a single plus and a single minus motor, respectively. The sign of the force is chosen to be positive if it is a load on the plus motors. The force acting on a single plus motor is then given by F + = F ca (n + ,n À )/n + . This means that all plus motors feel the same load in a given cargo state (n + ,n À ). The effective unbinding rate for a plus motor is: e + (n + ,n À ) = n + e + 0 exp(|F ca (n + ,n À )|/(n + F + d )).
Note that multiplying single motor rates by the factor n + implies that the motor-motor interactions are not considered. The effective binding rate of a plus motor is similarly given by: The effective binding and unbinding rates for the minus motors can be obtained by replacing the index ''+'' by ''À'' in eqn (A.2) and (A. 3). An expression for the cargo force F ca (n + ,n À ) of the force balance condition can be determined by observing that both plus and minus motor teams match their velocities under this load. The corresponding velocity is equal to the cargo velocity: v ca (n + ,n À ) = v + (F ca (n + ,n À )/n + ) = Àv À (F ca (n + ,n À )/n À ), (A.4) where the functions v + (Á) and v À (Á) are determined by the single plus and minus motor force-velocity relations. In ref. 6 and 7 the force-velocity relation of a single motor is given by the following piecewise-linear function: From eqn (A.5) we obtain the velocity of plus motors for the state (n + ,n À ): The minus motor velocity is given by an analogous expression by replacing the index ''+'' by ''À''. The cargo force F ca (n + ,n À ) can now be determined by the velocity matching condition in eqn (A.4): Observe that when only one motor type is active, i.e. n À = 0 or n + = 0, the force acting on the cargo disappears: F ca (n + ,n À ) = 0. From eqn (A.7) and (A.4), we obtain the cargo velocity: For the empirical force-velocity relationship v(F) as given in the main text by eqn (2.1) the matching condition (A.4) can be solved numerically. In the case of two antagonistic motors the cargo velocity defined by the matching condition eqn (A.4) can be obtained graphically from the intersection point (F,v) = (F ca ,v ca ) of two force-velocity relations, as shown in Fig. 2.
For the case of ''stronger plus motors'', i.e. n + F + s 4 n À F À s , we have: where L = (1 + (n + F + s v À B )/(n À F À s v + F )) À1 . Observe that L can only have values in the interval [0,1], which implies that the cargo force F ca (n + ,n À ) ranges between the maximal values of the plus and minus stall forces n + F s+ and n À F sÀ .
In this work we only consider two opposing motors, i.e. the state (n + = 1, n À = 1). The notation used in the main text is F ca F ca (n + = 1, n À = 1) and v ca v ca (n + = 1, n À = 1), where F ca and v ca are obtained from the numerical solution of the velocity matching condition, i.e. by determining the coordinates of the intersection point of the force-velocity functions of both motors.

B Master equations and the specification of parameters
The full network in Fig. 4 can be reduced to a closed network with all transitions to the two absorbing states being redirected to the initial state (0), see Fig. 12. The steady state distribution of the closed network can now be used to obtain some quantities of interest such as the average absorption time for the full network. 11,33 The latter quantity can also be calculated recursively without constructing the closed network. 34 A single trajectory on the closed network in Fig. 12 corresponds to the concatenation of many full network trajectories, each of which starts at the initial state (0) and is eventually absorbed in the states (1, +) or (1, À). The master equations corresponding to the closed network in Fig. 12 are given by for the two boundary states with j = ÀJ and j = +J, by o off ð jÞp j ðtÞ; ðB:3Þ for the state (0), where the sum includes all states ( j) apart from (0), and The different rates that appear in these equations are defined in (2.22)-(2.25) and depend on the single motor properties described by (2.1)-(2.6). The latter relationships involve the single motor parameters v 0 , v max , v min , F s , q 0 , F d , e 0 , and l as well as the elastic coupling parameter K = k/2. We typically vary one parameter such as the unbinding rate e 0 or the elastic coupling parameter K, keeping all other parameters fixed at their values in Table 1.

C Matrix form of master equations
The master equations can be written in the matrix form if we define the (2J + 1)-dimensional column vector (p ÀJ (t),p ÀJ+1 (t),. . .,p JÀ1 (t),p J (t)) T |p(t)i (C.1) where the superscript T stands for 'transpose' and the ket notation will be used for convenience. Using the latter vector, the master equations (B.1)-(B.4) attain the compact form with the (2J + 1) Â (2J + 1) transition rate matrix W. The diagonal matrix elements W j, j are given by for j a ÀJ, j a 0, and j a J The off-diagonal matrix elements of W are given by and W j, jÀ1 = o f ( j À 1) for j Z ÀJ + 1.
For each column of the matrix W, the matrix elements sum up to zero, i.e. P j W j;k ¼ 0 for all values of k.
Most of our results are based on the steady state solutions The numerical solutions of these equations were obtained with a self-written script based on built-in subroutines of Mathematica 9.0.

D Time evolution without motor unbinding
In Section 3.3, we describe the time evolution of the probability distribution p j (t) for unbinding rate e 0 = 0, see Fig. 8. In this limiting case, the state space in Fig. 4 becomes 1-dimensional and does not contain any cycles. The steady state probability distribution p st j then satisfies the detailed balance conditions for all j with ÀJ r j r J À 1 which is equivalent to for all nonzero matrix elements. Applying the conditions (D.1) iteratively along the 1-dimensional state space, we obtain the steady state solutions where p st ÀJ is determined by the normalization condition P j p st j ¼ 1.
In general, the matrix W is not symmetric. However, if the matrix elements W i,j fulfill the detailed balance conditions (D.2), we can define the symmetric transition rate matrix W with matrix elements 41 Because the matrix W is symmetric, it has real eigenvalues l n with n = 1,2,. . .,N and N = 2J + 1. Furthermore, the right and left eigenvectors for the eigenvalue l n have the same components ũ j (l n ), i.e. X kW j;kũk l n ð Þ ¼ l nũj l n ð Þ and X jũ j l n ð ÞW j;k ¼ l nũk l n ð Þ (D. 5) and satisfy the orthonormality condition X jũ j l m ð Þũ j l n ð Þ ¼ d m;n (D.6) Because all eigenvalues l n are negative for n Z 2, the long time behavior of p j (t) is dominated by the term with n = 1 and l 1 = 0 corresponding to the steady state probabilities p st j . Thus, in order to determine the time-dependent probabilites p j (t) as given by (D.11), we subsequently computed the steady state probabilities p st j , see (D.3), the symmetric matrix W via (D.4), as well as the eigenvalues l n of W and the corresponding eigenvectors with components ũ j (l n ), using again the built-in subroutines of Mathematica 9.0. The initial distributions were set to p j (t = 0) = 1 for j = 0 and p j (t = 0) = 0 for all other values of j.

E Dependence of average elastic force on strain force
In Fig. 9 the strain force dependency of the average elastic force hFi is presented up to a strain force value of F K C 3 pN. In this regime, the average elastic force is approximately independent of the strain force and equal to the cargo force F ca of the tug-ofwar in the MKL model. In the limit of high strain forces, however, we would expect the elastically coupled model to generate hFi = 0, since none of the motors would be able to make a single step because of the high spring stiffness.
In Fig. 13 we plot the average force hFi as a function of the strain force F K within the range 3 pN r F K r 40 pN. We observe a weak maximum of hFi around a strain force value of F K C 10 pN.
For F K Z 15 pN, the average force hFi decreases monotonically towards zero. Considering the steady state distribution over the discrete state space one can interpret the F K -dependence of hFi as follows. Up to a strain force value of F K C 2F + s = 2F À s the process can perform at least one step resulting in a stretching of the linker between the motors. Recall that the transition from state (0) to state (1) occurs with the sum of forward stepping rates a + ( % F 04 ) + a À ( % F 04 ), see eqn (2.22), where the motors feel the effective force % F 04 = F K /2. For a strain force of F K = 14 pN this force is % F 04 = 7 pN = F + s = F À s , at which the motors are equally likely to perform a forward or a backward step. Without motor unbinding, and taking into account that a further transition from (1) to (2) with an effective force of % F 14 = 14 pN acting on the motors is highly unlikely, this force balance is expected to lead to steady state probabilities p st 0 C p st 1 C 0.5 and to the average elastic force hFi C p st 0 F 0 + p st 1 F 1 = 7 pN = F ca . † A further increase in the strain force F K shifts the steady state distribution towards the relaxed state (0), i.e. p st 0 4 p st 1 , resulting in a reduction of the average force hFi towards zero. The initial increase of hFi with increasing F K implies that the corresponding increase in F 1 overcompensates the decrease of p st 1 for F K o 10 pN, see the steady state distribution in Fig. 14.
In Fig. 14 we plot the time evolution of the probability distribution over the states p j (t) starting from an initial state different from (0). We fix a high strain force F K = 10 pN, corresponding to the maximum of hFi in Fig. 13, and choose the initial state to be (5), i.e. p j (0) = d j5 . As t increases we see that the time-dependent probability distribution p j (t) approaches the steady state distribution p st j (dashed red line). We see that p st j has a clear peak in state (1), but also a nonzero value in (2), where the corresponding elastic force is given by F 2 = 2F K = 20 pN. These high force value contributions from state (2) to the average force results in hFi 4 F ca , as observed in Fig. 13. Fig. 13 Average elastic force hFi of the process for zero-force unbinding rate e 0 = 0 and strain forces F K Z 3 pN. We observe that hFi reaches a maximum around a strain force value of F K = 10 pN. For larger strain forces corresponding to stiffer springs, the average force hFi induced by the motors decreases monotonically to zero.
342 | Soft Matter, 2017, 13, 328--344 This journal is © The Royal Society of Chemistry 2017 F Dependence of force distribution on the size of state space As shown in Fig. 10, the probability distribution p st j is shifted towards larger j-values and becomes broader as we decrease the strain force F K . Therefore, a fixed size of the state space, i.e. a fixed number (2J + 1) of states, acts to truncate the distribution p st j for sufficiently small values of F K . Likewise, the average elastic force hFi starts to decrease, for a fixed value of F K , as we decrease the size (2J + 1) of the state space below a certain minimal value.
These finite size effects are illustrated in Fig. 15 where the average elastic force hFi is plotted as a function of the inverse size (2J + 1) À1 of the state space. Inspection of this figure shows that hFi attains a constant value for sufficiently small (2J + 1) À1 but decays strongly for sufficiently large (2J + 1) À1 . In practice, we chose (2J + 1) Z 201 for F K r 0.1 pN and (2J + 1) Z 51 for F K r 0.5 pN.

Glossary of mathematical symbols a(F)
Force-dependent forward stepping rate b(F) Force-dependent backward stepping rate  Fig. 15 Dependence of average elastic force hFi on the inverse size (2J + 1) À1 of the state space. If the two motors are weakly coupled, the finite size of the state space acts to truncate the probability distributions. Thus, for F K = 0.1 pN, we see large finite size effects and a strong decay of hFi already for (2J + 1) À1 \ 0.005. Fig. 14 Time evolution of the probability distribution p j (t) for a strain force of F K = 10 pN. The initial probability distribution is p j (0) = d j5 . The probability distribution p j (t) approaches the steady state distribution p st j , indicated by the dashed red line, which has a maximum for spring extension j = 1.