Simulations suggest robust microtubule attachment of kinesin and dynein in antagonistic pairs

Intracellular transport is propelled by kinesin and cytoplasmic dynein motors that carry membrane-bound vesicles and organelles bidirectionally along microtubule tracks. Much is known about these motors at the molecular scale, but many questions remain regarding how kinesin and dynein cooperate and compete during bidirectional cargo transport at the cellular level. The goal of the present study was to use a stochastic stepping model constructed by using published load-dependent properties of kinesin-1 and dynein-dynactin-BicD2 (DDB) to identify specific motor properties that determine the speed, directionality, and transport dynamics of a cargo carried by one kinesin and one dynein motor. Model performance was evaluated by comparing simulations to recently published experiments of kinesin-DDB pairs connected by complementary oligonucleotide linkers. Plotting the instantaneous velocity distributions from kinesin-DDB experiments revealed a single peak centered around zero velocity. In contrast, velocity distributions from simulations displayed a central peak around 100 nm/s, along with two side peaks corresponding to the unloaded kinesin and DDB velocities. We hypothesized that frequent motor detachment events and relatively slow motor reattachment rates resulted in periods in which only one motor is attached. To investigate this hypothesis, we varied specific model parameters and compared the resulting instantaneous velocity distributions, and we confirmed this systematic investigation using a machine learning approach that minimized the residual sum of squares between the experimental and simulation velocity distributions. The experimental data were best recapitulated by a model in which the kinesin and dynein stall forces are matched, motor detachment rates are independent of load, and the kinesin-1 reattachment rate is 50 s-1. These results provide new insights into motor dynamics during bidirectional transport and put forth hypotheses that can be tested by future experiments. Statement of significance: Bidirectional transport of vesicles along microtubules is vital for cellular function, particularly in the highly elongated axons and dendrites of neurons, and transport defects are linked to neurodegenerative diseases. For developing future therapeutic strategies, a better understanding is needed for how motors, cargo adapters, and accessory proteins coordinate their activities to transport cargo to their proper cellular locations. We approached this problem by simulating how antagonistic kinesin and dynein motors compete in pairs. We constrain our simulations by recent experimental results and conclude that the motors spend nearly all their time attached to the microtubule and competing against one another. This behavior is not predicted by existing single-molecule experiments and thus provides new insights into bidirectional transport.

one motor is attached. To investigate this hypothesis, we varied specific model parameters and compared 23 the resulting instantaneous velocity distributions, and we confirmed this systematic investigation using a 24 machine learning approach that minimized the residual sum of squares between the experimental and 25 simulation velocity distributions. The experimental data were best recapitulated by a model in which the 26 kinesin and dynein stall forces are matched, motor detachment rates are independent of load, and the 27 kinesin-1 reattachment rate is 50 s -1 . These results provide new insights into motor dynamics during 28 bidirectional transport and put forth hypotheses that can be tested by future experiments. 29 30 Statement of significance: 31 Bidirectional transport of vesicles along microtubules is vital for cellular function, particularly in the highly 32 elongated axons and dendrites of neurons, and transport defects are linked to neurodegenerative 33 diseases. For developing future therapeutic strategies, a better understanding is needed for how motors, 34 cargo adapters, and accessory proteins coordinate their activities to transport cargo to their proper 35 cellular locations. We approached this problem by simulating how antagonistic kinesin and dynein motors 36 compete in pairs. We constrain our simulations by recent experimental results and conclude that the 37 motors spend nearly all their time attached to the microtubule and competing against one another. This 38 behavior is not predicted by existing single-molecule experiments and thus provides new insights into 39 bidirectional transport. 40 41 Introduction: where R is a uniformly distributed random number in range 0 to 1. In a system with N possible events, the 99 rate of any event occurring equals the sum of the rate constants for all possible events. Thus, at each time 100 point, a two-stage process (the 'direct method' of Gillespie (41)) was used to determine the time and 101 identity of the next transition. First, a random number ! is generated and used to compute the time to 102 the next event, i as follows: 103 In the second stage, a new round of random number generation is used to determine which of all possible 104 events occurs; here, the probability of any event occurring is proportional to its rate constant (27,41). We 105 split this stage into two steps, with the goal of modularizing the code and enabling expansion to multiple 106 motors bound to the same cargo. We define j = 1,…M motors in the simulation (for our  simulations, M = 2). For each motor, j, there are i = 1,…N possible transitions, where N = 4 in most cases 108 for our model (forward step, backward step, attachment and detachment). In the first step, a uniformly 109 distributed random number, % is generated and used to choose which of the M motors present will make 110 a transition, and in the second step a random number & is used to choose which of the N possible 111 transitions for that motor will occur. Formally, using M motors, the probability of motor j, which has N 112 possible transitions acting is: 113 From N possible events on motor j, the probability of event i with rate constant kjj occurring is: 114 To identify the next event, the unit interval is divided into segments 115 corresponding to each possible transition: *0, , 1-. A random number is then 116 used to determine which transition will occur. This method is repeated until all motors detach from the 117 microtubule or a maximum run length or run duration is reached. 118 In the simulations, force is defined as positive in the plus-end direction and negative in the minus-end 119 direction; thus, for kinesin-1 negative forces are hindering loads and positive forces are assisting loads. 120 The force, F, applied to each motor is defined as the motor-cargo distance, Δ , multiplied by stiffness of 121 the motor, -./.0 : 122 In our previous model (27), we used )"1 = 0.3 pN/nm and 223 = 0.065 pN/nm based on published 123 work (42)(43)(44). For investigating the role of stiffness in the present work, we simplified the model by 124 incorporating an identical stiffness, -./.0 = 0.1068 pN/nm for each motor, which sets the cargo at the 125 midpoint between two motors. 126 At every timepoint, the position of the virtual cargo is set by computing a force balance between the two 127 motors when both motors are attached, and by setting it to the position of the attached motor if only one 128 motor is attached. Thus, following a motor step or detachment, the time for the system to reach 129 mechanical equilibrium is assumed to be negligible. This instantaneous force balance can be justified by 130 considering the relaxation time of the ~30 nm quantum dot attached to the Kin-DDB complex in the 131 experiments. The drag coefficient gQdot of a microsphere with radius rQdot = 15 nm in an aqueous solution 132 (viscosity coefficient h = 10 -9 pN×s/nm 2 ) is given by the Stokes' law as gQdot = 6phrQdot = 2.8 ´ 10 -7 pN×s/nm 133 (45). Following a motor step, we can calculate the characteristic relaxation time constant for the particle 134 to exponentially relax back to an equilibrium position as = 45./ / -./.0 , where -./.0 = 0.1068 135 pN/nm. This calculation yields a 3-microsecond relaxation time constant, which is negligible compared to 136 the motor stepping and attachment/detachment rates in the msec range or slower. 137

DDB state-switching model 173
In a subset of simulations, we incorporated DDB state switching based on work from Feng et al. showing 174 that DDB switches between processive, diffusive and stuck states during movement (34). To simulate this 175 behavior, three motility states for DDB were integrated into the model, as follows. In the processive state, 176 DDB can move forward, backward, or detach; in the stuck state, DDB cannot move or detach from the 177 microtubule; and in the diffusive state, DDB offers no resistance to kinesin movement and can also detach 178 from the microtubule. From every state, the motor can transition to either of the remaining states, with 179 transition rates as follows: the processive-to-stuck switching rate was 1.0 ,! with reverse rate of 1.8 ,! ; 180 the stuck-to-diffusive switching rate was 0.07 ,! with a reverse rate of 0.33 ,! ; and the diffusive-to-181 processive switching rate was 3.9 ,! with a reverse rate of 0.23 ,! (34). In the state-switching model 182 simulations, these transitions were added to the stepping and attachment/detachment transition 183 possibilities for DDB. For the majority of the simulations, this state switching was turned off. 184

Data processing 185
After running the simulations, the raw data were processed to match the experimental conditions (34), 186 as follows. To match the 20 fps frame rate (34), cargo position was averaged over 50 ms windows. To 187 account for uncertainties in fitting the experimental point spread function, a normally distributed error 188 with an 8 nm standard deviation was added to cargo position at each time point. The instantaneous 189 velocity at each timepoint was calculated by taking three-point slope of positions 50 ms before and after 190 each point. 191

Parameter optimization 192
As an alternate approach for parameter identification, we used a Bayesian Optimization method for 193 minimizing the residual sum of squares (RSS) between the experimental and model instantaneous velocity 194 distributions. The instantaneous velocity probability density functions from -2000 nm/s to 2000 nm/s with 195 10 nm/s bin widths were calculated for both experimental data and simulation results, where " is the 196 PDF of instantaneous velocity i. For each parameter set, a velocity probability density function was 197 generated from 1000 trajectories having a maximal run time of 50 s. The Matlab function 'bayesopt' was 198 used with 100 iterations to identify parameter sets that minimized the RSS function: 199

Formulation of bidirectional stepping model 202
In recently published experiments, we reconstituted bidirectional cargo transport in vitro by connecting a 203 truncated kinesin motor to an activated dynein-dynactin-BicD2 (DDB) complex through complementary 204 single stranded DNA oligonucleotides (34). A quantum dot was then linked to a biotin on one end of the 205 DNA (Fig. 1A), and the position of the fluorescent cargo was tracked via total internal reflection 206 fluorescence (TIRF) microscopy at 20 frames per second. Consistent with previous work (33), the resulting 207 complexes moved slowly for long durations along immobilized microtubules, with some complexes 208 moving net plus-end and others moving net minus-end ( given in Fig. S2). The simulated cargo has neither mass or drag, and its position is updated following every 216 step to maintain force balance. The parameters for this 'basic model' are shown in Table 1  In contrast, the simulated instantaneous velocity distribution had a distinct central peak that was right 254 shifted compared to the experiments, along with distinct side peaks corresponding to the unloaded 255 velocities of DDB and kinesin-1. Using a Gaussian mixture model (GMM), the velocity distribution was fit 256 well by three normal distributions: a peak centered at 127 nm/s that accounted for 53% of the population, 257 a minor peak at -334 nm/s that accounted for 41% of the population, and another minor peak centered 258 at 587 nm/s that accounted for 6% of the population ( Table 2). The motor reattachment rate in our 259 simulations was 5 s -1 based on published experimental data for kinesin-1 (48,52,53), and because 260 equivalent experimental data are not available for dynein, we chose the same value for the DDB complex 261 (Table 1). This relatively slow reattachment rate means that when either motor detaches from the 262 microtubule, it will take on average 200 ms to reattach, giving the other motor time to move at its 263 unloaded velocity. Thus, these lateral peaks in the simulated velocity distribution can be explained by the 264 simulations having more periods when only one motor is attached than the experiments, and their 265 positions are close to the unloaded kinesin and DDB velocities. 266

270
To better understand the underlying motor dynamics that lead to the experimental Kin-DDB traces being 271 relatively smooth and the velocity distribution having a single peak centered around zero, we investigated 272 how the modeled properties of kinesin-1 and DDB contribute to the distinctive instantaneous velocity 273 distribution of our 'basic' model described above. We hypothesize that there are properties of kinesin-1 274 and/or DDB motors that differ when the motors are joined in a two-motor complex compared to in their 275

Testing the influence of dynein motor properties on bidirectional transport 290
The first hypothesis we tested was that the simulations differ from the experiments because the motile 291 properties of DDB in a bidirectionally moving Kin-DDB complex differ from the DDB properties measured 292 in single-motor optical tweezer experiments. To test this possibility, we first investigated whether either 293 including experimentally observed pausing and diffusive states of the DDB complex could better align the 294 simulations with the experiments. In virtually all published unloaded single-molecule studies, DDB exhibits 295 long processive runs, but also undergoes episodes of 1D diffusion along the microtubule and spends a 296 significant fraction of the time stuck to the microtubule in an immobilized state (32,34,37,38). Feng et al 297 found that DDB spent 65% of the time in a processive state, 31% of the time in a stuck state, and 4% of 298 the time in a diffusive state, and also quantified the switching rates between states (34). We incorporated 299 this DDB switching behavior as a series of transitions into our model. We assumed that the diffusive state 300 of DDB offered no resistance to kinesin stepping, and that in the stuck state, DDB neither moved nor 301 detached from microtubule. In the instantaneous velocity distribution, incorporating this DDB switching 302 behavior caused a decrease in the magnitude of the unloaded DDB velocity peak, but had no effect on the 303 unloaded kinesin velocity peak (Fig. 3A). The fall in the unloaded DDB peak can be explained by DDB 304 spending time in the stuck state instead of processive walking. The lack of change in the unloaded kinesin 305 peak is likely due to the fact that DDB is in the diffusive state only a small fraction of the time (Fig. 3A). 306 Thus, DDB state switching cannot account for the discrepancy between model and experimental results. 307 Due to this lack of an effect, all subsequent simulations used only a single processive state for DDB. 308 The second property we examined was the load-dependent detachment dynamics of DDB. There is 309 abundant evidence in the literature that dynein alone has catch-bond behavior, meaning that the 310 detachment rate slows down with increasing load (26,29,(54)(55)(56)(57). There is no experimental evidence to our 311 knowledge that activated dynein in a DDB complex shows catch-bond behavior; however, only one 312 published study directly addresses this question (33). Thus, we investigated whether the simulated 313 bidirectional behavior is shaped by the load-dependent detachment kinetics of DDB by using an 314 exponential dependence of the off-rate on load (Methods,Equation 8). In this formulation, a positive 315 detachment force parameter ( 5>/9;? ) corresponds to a slip-bond, a negative value corresponds to a 316 catch-bond, and a value of infinity corresponds to an ideal-bond. We found that slowing the unloaded 317 detachment rate ( 5>/9;? 6 ) or increasing the detachment force ( 5>/9;? ) decreased the peak 318 corresponding to the unloaded kinesin-1 velocity, bringing the distribution more in line with the 319 experimental ( Fig. 3B and C). When an ideal bond was simulated by setting 5>/9;? to infinity, the 320 instantaneous velocity distribution was indistinguishable to that for 5>/9;? = 15 pN (Fig. S3A). Setting 321 5>/9;? to -3 pN, corresponding to a catch-bond resulted in no change from the ideal-bond model (Fig.  322 S3A). 323 The third property we examined was the DDB stall force, Fstall. In the literature, stall forces for isolated 324 dynein were measured to be 1 pN by a number of labs (29), but activated dynein in a DDB complex was 325 later shown to have a stall force of 3.6 pN (33). In our basic model, the kinesin stall force is twice that of 326 DDB (Table 1) and the central peak in the instantaneous velocity distribution is shifted towards the plus-327 end by nearly 150 nm/s relative to the experimental peak ( Table 2; Fig. 2C); thus, it is reasonable to expect 328 that bringing the stall forces more into alignment should shift the central peak. We found this to be the 329 case ( Fig. 3D) -increasing the DDB stall force to 6 pN shifted the central peak leftward and increasing it 330 to 8 pN to match the kinesin stall force brought the central peak nearly into alignment with the 331 experimental peak. As a final investigation, we examined the sensitivity of the simulations to the DDB 332 reattachment rate by increasing 9//9;? from 5 s -1 to 50 s -1 . This modification had virtually no effect on 333 the simulated velocity distribution (Fig. S3B). This lack of an effect likely stems from the slow DDB 334 unloaded detachment rate of 0.1 s -1 -because DDB detachments are relatively infrequent, the unbound 335 episodes are rare and shortening their duration has little effect. 336 To summarize the DDB parameter investigations, removing the load-dependence of DDB detachment by 337 modeling DDB as an ideal-bond nearly eliminated the kinesin peak in the instantaneous velocity 338 distribution, and increasing the DDB stall force shifted the central velocity peak leftward, closer to the 339 experimental peak. In contrast, incorporating a switching model for DDB or increasing the DDB 340 reattachment rate had only minimal effects on the velocity distribution. Next, to explore other 341 modifications of the model that could diminish the DDB velocity peak, we examined the sensitivity of the 342 model simulations to changes in the kinesin detachment and reattachment rates. Testing the influence of kinesin-1 motor properties on bidirectional transport 354 The second hypothesis we tested was that the discrepancy between simulations and experiments is due 355 to differences in the motor properties of kinesin when operating in an antagonistic motor pair compared 356 to single-bead optical tweezer experiments with isolated kinesin motors. These differences could be due 357 to, for instance, the different geometries of being attached to a trapped bead versus being tightly 358 connected to dynein through a short DNA (31,58,59), and they could in principle affect both motor on-359 and off-rates. The first investigation was to determine whether decreasing the detachment rate of kinesin 360 by reducing 5>/9;? 6 or increasing 5>/9;? reduced the magnitude of the DDB velocity peak in the 361 instantaneous velocity distribution. Decreasing 5>/9;? 6 for kinesin-1 substantially decreased the 362 magnitude of the unloaded DDB velocity peak and increased the magnitude of the central peak, although 363 it also increased the unloaded kinesin velocity peak (Fig. 4A). Increasing 5>/9;? to reduce the sensitivity 364 of detachment to load had a similar effect but to a lesser degree (Fig. 4B), with the effect plateauing 365 around an 5>/9;? value of 35 pN. Thus, decreasing the kinesin detachment rate diminished the minus-366 end velocity peak, bringing the simulations closer to the experimental results, but neither of these changes 367 had a significant effect on the unloaded kinesin peak. 368 The other way to decrease the fraction of the time Kin-DDB complexes are moving solely by DDB is to 369 increase the kinesin-1 reattachment rate. In the 'basic' model, we set a constant reattachment rate, 370 0>9//9;? , of 5 s -1 for both motors, based on the literature (28,53). To determine whether the kinesin-1 371 reattachment rates were contributing to the simulation-experiment mismatch, we tested several larger 372 values of 0>9//9;? for each motor. In contrast to DDB where altering the reattachment rate had little 373 effect (Fig. 3D), increasing 0>9//9;? for kinesin-1 both strongly diminished the fraction of time the 374 complex moves at the unloaded DDB speed and increased the fraction of time that both motors are 375 engaged, resulting in an enhanced cargo velocity peak near zero (Fig. 4C). 376 The final kinesin parameter we investigated was the kinesin stall force, </9== . In our 'basic' model, the 377 kinesin-1 stall force was set to 8 pN and the DDB stall force was set to 3.8 pN based on published optical 378 tweezer experiments (33,46,47). Because of this discrepancy, it is perhaps not surprising that when both 379 motors are engaged, the kinesin directionality dominates. Thus, we investigated kinesin stall forces 380 between 4 and 8 pN, which are in the range of published studies (29,46,(60)(61)(62)(63)(64), and found that weaker 381 kinesin stall forces cause the central velocity peak to shift closer to zero, better matching the experiments 382 (Fig. 4D). 383 In summary, a major discrepancy between the experimental and simulated velocity distributions is the 384 peak in the simulations near the unloaded DDB velocity. We found that this minus-end velocity peak could 385 be diminished by decreasing the kinesin-1 unloaded detachment rate or the sensitivity of detachment to 386 load, or by increasing the kinesin-1 reattachment rate. Additionally, decreasing the kinesin-1 stall force 387 resulted in a leftward shift of the dominant peak in the velocity distribution, bringing the simulations in 388 better alignment with experiments. 389

Tuning parameters in an improved bidirectional transport model 398
The simulation results to this point showed that altering motor attachment and detachment rates can 399 decrease the weight of the instantaneous velocity distribution corresponding to the two unloaded motor 400 velocities, and that bringing the kinesin and DDB stall forces into closer alignment can shift the location of 401 the central peak. Informed by these single parameter adjustments, we next formulated an updated model 402 that incorporated multiple parameter adjustments, with the goal of more closely matching the simulated 403 and experimental instantaneous velocity distributions. We chose not to modify the 5>/9;? 6 parameter for 404 either kinesin-1 or DDB because single-molecule TIRF measurements from many labs have consistently 405 measured unloaded detachment rates (equal to velocity divided by run length) near 1 s -1 for kinesin-1 and 406 0.1 s -1 for DDB (29,40,46,47). Our adjustments to kinesin-1 were guided by recent work by Pyrpassopoulos 407 et al. (58), who showed that kinesin-1 KIF5B under purely horizontal loads, achieved through a three-bead 408 optical tweezer geometry, has a ~6 pN stall force, a 1.1 s engaged duration (matching the unloaded run 409 duration) (58), and a component of its reattachment rate at 100 s -1 (31,58). Thus, we set the kinesin-1 stall 410 force to 6 pN and applied an ideal-bond detachment model by setting 5>/9;? ( ) = 5>/9;? 6 . We also 411 increased the kinesin-1 reattachment rate to 50 s -1 ; a 100 rate s -1 was not used because there was little 412 improvement beyond 50 s -1 and because in the experiments there were components slower than 100 s -1 413 (58). There is no comparable three-bead study for DDB, but we hypothesize that removing the vertical 414 force component may have similar effects on dynein. Thus, we also incorporated a load-independent 415 detachment rate (ideal bond) for DDB. 416 The parameter adjustments were incorporated into an upgraded model that we named the 'better-fit' 417 model (Table S2). Simulations showed that in the 'better-fit' model the plus-end side peak was eliminated, 418 the minus-end side peak was diminished, and the central peak was substantially enhanced but was still 419 right-shifted compared to the experimental data (Fig. 5A). Starting from this improved model, we tested 420 further parameter adjustments, with the goal of shifting the central velocity peak to the left and 421 eliminating the minus-end side peak. We first tested the DDB stall force and found that increasing it to 6 422 pN to match that of kinesin-1 substantially shifted the main peak leftward. Next, we varied the DDB 423 backward stepping rate and found that decreasing it to 3 s -1 or 1 s -1 also substantially shifted the main 424 peak leftward (Fig. 5C and S3C). 425 The final property we tested was the motor stiffness, which affects how rapidly motors ramp up their 426 force as they step in opposite directions. We found that increasing motor stiffness eliminated the fastest 427 minus-end velocities in the distribution, and that this weight was shifted to the central peak, bringing it 428 closer to the height of the experimental peak ( Fig. 5D and S3D). Upon closer examination, we surmised 429 that this shift results from eliminating a 'recoil effect' occurring at lower stiffness values, in which the 430 kinesin detaches and the cargo rapidly moves toward the minus-end to become centered on the engaged 431 DDB motor. A 6 pN force is expected to stretch the spring that connects each motor to its cargo by 60 nm. 432 A recoil of this magnitude over the 0.1 s window used for the velocity distributions corresponds to a -600 433 nm/s velocity, and this is the component that is eliminated when motor stiffness is increased (Fig. 5D). In 434 summary, in the 'better-fit' model either increasing the DDB stall force or decreasing the DDB 435 backstepping rate shifted the central velocity peak leftward to better align with the experimental 436 distribution, and increasing the motor stiffness eliminated the remaining minus-end velocity side peak.

Using machine learning for parameter sensitivity exploration and optimization 448
Having narrowed the parameter choices by iterative tuning, we next applied an automated parameter 449 optimization method to identify parameter sets that best match the simulations to the experimental 450 velocity distributions. To do this we used the Matlab Bayesian Optimization function 'bayesopt', which 451 identifies an optimal parameter set based on an objective function. We defined the objective function as 452 the residual sum of squares (RSS) between the experimental and simulation velocity distributions, and we 453 defined upper and lower search bounds for each parameter (Table 3). In the first round, we allowed seven 454 parameters to be optimized: the detachment forces and reattachment rates for both motors, along with 455 the motor stiffness, DDB stall force and DDB backstepping rate. Because of the interdependency of the 456 parameters, we found that repeated optimization runs achieved similar RSS values using different 457 parameter sets. Hence, in Table 3 we report the range of optimal parameter values for 10 independent 458 optimization runs, with large ranges connoting parameters that either co-vary with other parameters or 459 do not strongly affect the final velocity distribution, and narrow ranges connoting parameters that most 460 strongly determine model performance. As a final step, we fixed the detachment forces and reattachment 461 rates and performed a second round of optimization on the remaining DDB stall force, DDB backstepping 462 rate, and motor stiffness parameters. The results of this final optimization were similar to the first, but 463 the optimal ranges were narrowed (Table 3). 464 The parameter optimization results can be summarized as follows. First, the optimal detachment forces 465 for both motors were relatively high, consistent with the ideal-bond model we implemented in our 466 'better-fit' model. Therefore, in our final 'best-fit' model, we set detachment to be independent of load 467 for both motors. Second, consistent with the 'better-fit' model, the kinesin reattachment rate converged 468 in a range around 50 s -1 for kinesin and around 5 s -1 or slower for DDB. Thus, we chose those values for 469 the reattachment rates. Third, the DDB stall force converged near the 6 pN stall force of kinesin-1, 470 supporting the observation in Fig. 5B that matching the stall forces helps to align the central velocity peak 471 with the experiments. Fourth, the DDB backstepping rate converged in a range of 1-5 s -1 . This rate is slower 472 than the experimental 15 s -1 (33), and reiterates the finding in Fig. 5C that slowing DDB backstepping shifts 473 the central velocity peak leftward. Based on the optimal fit with the reduced parameter set, we set the 474 DDB backstepping rate to 1 s -1 (see Discussion). Finally, the motor stiffness converged at 0.2 pN/nm, 475 consistent with our finding that a minimum stiffnesses in this range is needed to avoid the 'recoil' 476 phenomenon that contributes unwanted minus-end velocities. The final parameter choices that were 477 used for the 'best-fit' model are shown in Table 3, and the entire 'best-fit' parameter set is provided in 478    (Fig. 2A). Using a mean squared displacement analysis, the apparent diffusion 489 coefficient dropped from 11,800 % ,! in the 'basic' model to 1272 nm 2 s -1 in 'best-fit' model, which 490 was close to the 996 % ,! for the experimental data (Fig. S1B). More importantly, in the instantaneous 491 velocity distribution (Fig. 6B), the unloaded velocity peaks, which were a significant discrepancy in the 492 'basic' model, were nearly eliminated, and the central peak was shifted to near zero to match the 493 experimental data. One metric of this agreement is that in the 'best-fit' model, 95% of instantaneous 494 velocity values fall between -349 to +299 nm/s, similar to the -394 to +371 nm/s range for the 495 experimental data. Fitting the 'best-fit' velocity distribution to a Gaussian mixture model identified a 496 central peak centered at -8 nm/s and comprising 93% of weight, along with two side peaks centered at 497 452 nm/s and -318 nm/s, comprising 0.5% and 6.5% of the weight, respectively (Table 2). Finally, the 498 residual sum of squares (RSS) between experimental and simulation data dropped from 3.93*10 -5 in the 499 'basic' model to 2.76*10 -6 in 'best-fit' model (Table 2). Other features were also improved; for instance, 500 the simulated traces for the 'best fit' model also had longer durations, more closely matching experiments 501 (Fig. S4). Also, for the 'best-fit' model the distribution of velocities computed over each trace matched the 502 experimental values (Fig. 6C). 503 515 Transport of vesicles in axons and dendrites involves varying numbers and types of motors, as well as 516 regulation at multiple levels. Hence, to understand how kinesin and dynein motors attached to the same 517 cargo compete and coordinate during transport, the field has turned to reconstituting the in vitro motility 518 of pairs of kinesin and activated dynein (1,7,38,53,65). This approach allows an examination of how 519 differing mechanochemical properties of specific motor isotypes translate into effective movement 520 against an antagonistic partner. Consistent behavior of kinesin-dynein pairs that are observed across 521 different labs (20,32-34) include the following: 1) Kin-DDB complexes remain bound to microtubules for 522 ~tens of seconds (33,34), considerably longer than kinesin-1 or DDB alone; 2) trajectories include episodes, 523 lasting from seconds to tens of seconds, where the velocity of the complex is ~10-fold slower than the 524 unloaded velocities of the respective motors; and 3) cargo trajectories are quite smooth and show few if 525 any directional switches. In principle, it should be possible to recapitulate these motility behaviors using 526 Kin-DDB simulations that incorporate established parameters for kinesin and dynein behavior in isolation. 527 However, we find this not to be the case. Instead, simulations of Kin-DDB transport in our 'basic' model 528 show frequent directional switches and significant proportions of time when the complexes move at the 529 unloaded velocity of the motors (Fig. 2B). The discrepancies between the simulations and experiments 530 were clearly seen by comparing the instantaneous velocity distributions (Fig. 2C). The experimental data 531 were centered in a wide peak around zero velocity, whereas the simulated data had additional peaks 532 corresponding to the unloaded kinesin-1 and DDB velocities, along with a central peak centered around a 533 slow plus-end speed. Thus, we focused our in silico 'experiments' on identifying parameter adjustments 534 that reconciled the simulations with the experiments. 535 The first modification that helped to align the simulation results with the experiments was to make the 536 kinesin detachment rate independent of load in the range of the stall forces used. Although single-bead 537 optical tweezer experiments have found that kinesin detachment rates increase exponentially with force 538 (58,59,66), recent theoretical and experimental work has suggested that forces perpendicular to the 539 microtubule inherent to the geometry of these experiments may be contributing to the detachment rate. 540 In support of this, experiments using a three-bead assay that almost fully eliminate these vertical forces 541 found that the kinesin-1 detachment rate was approximately independent of load (58). Incorporating this 542 'ideal-bond' for kinesin into our model strongly diminished the minus-end side peak in the velocity 543 distribution; thus, our simulations support the hypothesis that against forces oriented purely parallel to 544 the microtubule, the kinesin detachment is independent of load. We also examined a catch-bond model 545 for kinesin, and found a further increase in the amplitude of the central velocity peak with increasing catch 546 bond properties (Fig. S5). Thus, our results are consistent with possibility that kinesin has catch-bond 547 behavior under some circumstances, but an idea-bond for kinesin was sufficient in our simulations to 548 match the experimental data. 549 We also found that incorporating a load-independent detachment rate for DDB diminished the fast plus-550 end velocities seen in the simulations, better matching the experimental data. This result can be 551 understood simply as DDB remaining bound to the microtubule against the pulling forces of kinesin, thus 552 preventing durations in which DDB is detached and kinesin is moving at its unloaded speed. A DDB ideal-553 bond model is supported by experiments with activator-free dynein and with isolated vesicles, which 554 found evidence that dynein has catch-bond behavior under some conditions (29,57), though we note that 555 a DDB optical trapping study by Elshenawy et al. found an exponential dependence of DDB detachment 556 with load (33). However, it remains a possibility that, as has been shown for kinesin-1 (30,33,59), vertical 557 forces inherent to the single-bead optical trapping assay may contribute to detachment of DDB. Thus, our 558 simulations put forth the testable hypothesis that when pulling exclusively against hindering loads 559 oriented parallel to the microtubule, DDB detachment is independent of load over the forces generated 560 by a single kinesin motor. As with kinesin, our results were also consistent with DDB catch bond behavior 561 ( Fig. S3A), but an ideal bond was sufficient to match the simulations to the experiments. 562 Incorporating a fast kinesin-1 reattachment rate was another important modification that helped to bring 563 the simulations into alignment with the experimental data. In particular, the fast kinesin-1 reattachment 564 rate diminished the minus-end velocity side peak, which can be understood as minimizing the time that 565 DDB moves at its unloaded velocity. Our basic model used a kreattach value of 5 s -1 , which was initially 566 determined in a study that measured motor-driven deformations of giant unilamellar vesicles (53), was 567 later supported by experiments that used DNA to connect two kinesins (48), and which has been 568 employed in a number of modeling studies (e.g. (28)). However, three recent optical tweezer studies 569 found that against hindering loads, kinesin-1 can slip backward in multiple 8 nm intervals and rapidly 570 reengage with the microtubule at rates consistent with our 50 s -1 reattachment rate (58,67,68). These 571 experiments suggest that against a hindering loads oriented parallel to the microtubule, kinesin-1 enters 572 a weak-binding state in which it slips backward, and then rapidly reestablishes a strong-binding state. We 573 also note that this 50 s -1 reattachment rate is still below the 125 s -1 reattachment rate predicted from 574 multiplying the bimolecular on-rate of kinesin in solution by the effective tubulin concentration when two 575 motors are connected through a DNA linker (48). Thus, our simulations put forth the testable hypothesis 576 that when pulling against hindering loads oriented parallel to the microtubule, the kinesin-1 reattachment 577 rate is 50 s -1 or faster. A possible mechanism to explain this phenomenon is that if a kinesin slips or 578 detaches against a load oriented parallel to the microtubule, the motor is presented with multiple sites of 579 reattachment as it slides backward. Further, it is possible that if the motor reattaches against a hindering 580 load, the force enhances the rate of ADP release, which maximizes the probability that the motor enters 581 a strongly-bound state and restarts its walking cycle. 582 The model adjustment that most influenced the position of the central peak in the instantaneous velocity 583 distribution was altering the relative stall forces of DDB and kinesin-1, highlighting that the antagonistic 584 motors rapidly establish a "draw" where each is pulling at near its maximum force. The paired 585 modification of reducing the kinesin-1 stall force from 8 pN down to 6 pN and increasing the DDB stall 586 force from 3.6 pN up to 6 pN shifted the central peak leftward by roughly ~140 nm/s to match the 587 experimental peak (Fig. 6B). A large body of experiments support a kinesin-1 stall force in the range of 6 588 pN (29,46,(60)(61)(62)(63)(64), justifying this modification. For DDB, one justification for increasing the stall force is 589 that the lower value we started with was taken from a study using single-bead tweezer geometry (33), 590 and it is possible that in the geometry of these DDB-kinesin pairs where forces are oriented solely parallel 591 to the microtubule, the DDB stall forces are closer to the ~6 pN of kinesin-1 (Fig. 3D). Related to the stall 592 force, we also found that slowing the DDB backstepping rate shifted the central velocity peak leftward. 593 This shift makes sense in cases where the kinesin stall force is greater than the DDB stall force, meaning 594 that DDB is under superstall forces and stepping backward at 15 s -1 (Fig. 5C and S3C). However, we note 595 that when stall forces were matched, the DDB backstepping rate did not have a strong effect on the 596 velocity distribution (Fig. S6). Hence, although single-bead optical tweezer experiments measured a DDB 597 backstepping rate of 15 s -1 (33), our simulations put forward the testable hypothesis that the DDB 598 backstepping rate may be smaller in the absence of any vertical forces. 599 Another model adjustment that had a surprising effect was increasing the stiffness of the linkages 600 connecting the two motors to the cargo. When model and experiments were compared by analyzing 601 instantaneous velocities, a 'recoil' effect becomes apparent following detachment of one of the motors 602 (usually kinesin in our simulations). Thus, the simulations provide a lower limit of 0.2 pN/nm for the motor 603 stiffness and make the testable prediction that if a more compliant DNA linkage were used to connect the 604 motors, the instantaneous velocity distribution would reveal apparent fast velocities due to large 605 displacement recoil events. Additionally, a more compliant linkage between the motors would be relevant 606 to the transport of large (0.1 -1 micron) vesicles in cells, where motor forces are sufficient to deform the 607 membrane. 608 Despite the large body of single-molecule data describing their motor properties, it remains challenging 609 to predict the bidirectional dynamics that will result from antagonistic kinesin and dynein motor pairs. In 610 the present work we find that by using experimental data to constrain model simulations, mechanistic 611 insights can be generated into how kinesin and DDB operate in antagonistic pairs. Our simulations support 612 a model (Fig. 6D) in which: 1) Over the range of forces generated by the motors, the detachment rate of 613 both kinesin-1 and DDB is insensitive to load; 2) In the two-motor geometry examined, the kinesin-1 and 614 DDB stall forces are similar; 3) kinesin-1 reattaches to the microtubule at 50 s -1 in the geometry of these 615 motor pairs; and 4) when modeled as linear springs, motor stiffness is at least 0.2 pN/nm. These model-616 generated hypotheses generate testable predictions for future single-molecule experiments. Furthermore, 617 by identifying motor parameters that are the strongest determinants of bidirectional transport, this work 618 provides a framework to interpret how diverse kinesins, different activating dynein adapters, and motor-619 cargo stiffness will affect the resulting bidirectional transport of intracellular cargo.