Nano-Particles Carried by Multiple Dynein Motors Self-Regulate Their Number of Actively Participating Motors

Intra-cellular active transport by native cargos is ubiquitous. We investigate the motion of spherical nano-particles (NPs) grafted with flexible polymers that end with a nuclear localization signal peptide. This peptide allows the recruitment of several mammalian dynein motors from cytoplasmic extracts. To determine how motor–motor interactions influenced motility on the single microtubule level, we conducted bead-motility assays incorporating surface adsorbed microtubules and combined them with model simulations that were based on the properties of a single dynein. The experimental and simulation results revealed long time trajectories: when the number of NP-ligated motors Nm increased, run-times and run-lengths were enhanced and mean velocities were somewhat decreased. Moreover, the dependence of the velocity on run-time followed a universal curve, regardless of the system composition. Model simulations also demonstrated left- and right-handed helical motion and revealed self-regulation of the number of microtubule-bound, actively transporting dynein motors. This number was stochastic along trajectories and was distributed mainly between one, two, and three motors, regardless of Nm. We propose that this self-regulation allows our synthetic NPs to achieve persistent motion that is associated with major helicity. Such a helical motion might affect obstacle bypassing, which can influence active transport efficiency when facing the crowded environment of the cell.


Introduction
Active, motor-protein-mediated transport is crucial for the intracellular conveyance of a large variety of cargos in eukaryotes. Notably, microtubule-associated (MT-associated) motor proteins-dynein and kinesin-play a cardinal role, fueling a variety of vital biological processes [1][2][3]. While dynein is responsible for transport towards the cell center, members of the kinesin family are mostly responsible for transport towards the cell periphery [4][5][6][7]. The dynamic interplay between these two classes of motion orchestrates the subcellular arrangement of organelles, e.g., mitochondria [8] and Golgi complexes [1]. In addition, different types of viruses have evolved to harness the dynein machinery for the efficient targeting of the nucleus, where the infection process of the host cell occurs. HIV and herpes-simplex virus, for instance, express nuclear localization signal (NLS) peptides [7] that recruit dynein from the cytoplasm [9][10][11][12][13]. Likewise, adenoviruses use other ligands (e.g., hexon) to engage the active transport mechanisms [14,15]. Recently, preliminary theoretical work has been performed to mimic these viruses in rationally designed cargos that can be used for drug delivery applications [16].
While knowledge of single dynein motility has increased extensively in the past decade [17][18][19][20][21][22][23][24], much less is known about the collective motility of multiple motors that carry a single cargo. Several observations have led to the prevailing belief that native cargos are carried by more than one motor [25][26][27][28][29]. Studies of supercoiled DNA plasmids, enriched by NFκB receptors, showed robust nuclear temporal localization [30], supporting multiple dynein recruitment mediated by NLS binding. Complexes of fluorescently labeled single-strand DNA, such as the VirE2 protein of Agrobacterium and the VirE2 itself, contain putative NLS regions that can recruit dynein [31].
In the absence of motor-motor elastic coupling, one finds (theoretically) an exponential increase of the run-length with an increasing number of cargo-associated motors [16,49,50]. However, motor-motor elastic coupling [48], combined with excluded volume interaction between motors, may strongly affect motion and lead to different characteristics from those of the single-motor behavior. Moreover, these theoretical studies do not account for essential features of dynein stepping, for example, they do not include the detailed locations of the dynein binding sites over the MT surface. In this spirit, as mentioned above, the single motor stepping was recently revisited to describe this two dimensional (2D) stepping of a single yeast dynein [40]. The 2D stepping model correctly accounts for the measured longitudinal step-size distribution and predicts a broad angular distribution of steps with a small right-handed bias. Elucidating the role of motor-motor coupling in realistic (biological) multi-motor complexes is crucial for understanding the behavior of actively transported native cargos.
Motivated by a previous theoretical work [16], we present here a combined experimental-theoretical research on rationally designed particles that can serve as a model particle for native multi-motor complexes, on the one hand, and allows a quantitative examination of their transport characteristics and their association with motor-motor coupling, on the other hand. Experimentally, our strategy is based on grafting a spherical nano-particle (NP) with a prescribed grafting density of Biotin-polyethylene-glycol-thiol (Biotin-PEG-thiol) molecules and end-linking a controlled fraction of the Biotin-PEG-thiols with a single NLS ( Figure 1A). We expected, therefore, that under exposure to a cell extract (CE), the NLS peptide would first recruit an α-importin protein and then bind a β-importin protein, which, subsequently, would recruit dynein. Thus, the Biotin-PEG-thiol, which connects the dynein and the NP, serves as a spacer polymer with variable flexibility (depending on its contour length, i.e., molecular weight). By stretching, these spacers allow for motors originating from the sides of the NP to readily reach the MT surface (Figure 2A,B), thereby, increasing the number of motors participating in the transport. It follows that the structural properties of the NP, such as its Biotin-PEG-thiol grafting density, linked NLS fraction, Biotin-PEG-thiol contour length, NP size, and CE concentration, will, together, govern the number of dynein motors that bind to the MT simultaneously, thereby, potentially controlling the motility of the NP.  Table 1, Supplementary Material Section 1, and Figure S2). Error bars indicate the standard deviations for three experiments. (C) Western blot (WB) analysis results, demonstrating the recruitment of α2-importin (a); mammalian dynein motors (b); and dynactin to the NPs after incubation in a Hela cells extract. Group 1 refers to NPs coated with Biotin-PEG-thiol, group 2 refers to PEG-NLS coated NPs, and group 3 refers to Hela cell extracts without NPs. (D) WB analysis results, showing that kinesin-1 does not bind the PEG-NLS coated NP. Group 1 refers to PEG-NLS coated NPs, and group 2 refers to Hela cell extracts without NPs. Table 1. Summary of the systems studied and their basic properties: NLS and CE concentrations, estimated mean number 〈 〉 of PEG-NLS per NP, estimated mean anchoring distance * , and fractions of NPs performing distinct modes of motion. An identical cell extract was used in systems I and II, and another extract was used in systems III and IV. The bare NP mean radius is 20 nm.  . The grey dots correspond to extrapolated values that were calculated from the fit of <N> vs. <NLS> (see Table 1, Supplementary Material Section S1, and Figure S2). Our NPs are quite different from the cargo designs that have been previously reported [37][38][39]. First, previously reported cargos do not involve the native protein assembly for the recruitment of dynein from the cytoplasm. Second, their size is one order of magni-tude larger than the NPs used in this study. (Note that the radius of a native cargo size is typically in the range of 20, 100 nm, similar to the size of the NP used here). Importantly, in this radius range, the drag force on the NP is entirely negligible. (Using the NP translational Stokes drag coefficient, γ t = 6πηR, and a velocity of v = 4 µm/s, leads to a drag force f = γ t v ∼ 10 −2 pN, for η = 10 mPa s as an upper bound for the CE viscosity, which is much smaller than the typical motor force of 4 pN.). This implies that the single motor will move at its free-load velocity and without a load-sharing effect [52]; such an effect will occur only for micron-sized NPs, for which the drag force becomes significant [38,53]. Third, the cargos used in previous studies [37,38] cannot control the mechanical coupling between motors. As such, they are less suitable for a fundamental experimental-theoretical study of a model multi-motor cargo, which could lead to an understanding of the motion of native cargos. We also note that, to the best of our knowledge, previous transport studies of dynein-recruiting cargos that use NLS ligands and are exposed to CEs have been performed only in 3-dimensional MT networks [7], which did not allow the characterization of the motion on individual MT tracks.
In this work, we performed-using our NP constructs-bead-motility assays with total internal reflection fluorescence (TIRF) microscopy. We used single-particle tracking algorithms to extract the trajectories of the NPs, and we measured several properties of these trajectories: (i) modes of motion, such as directional motion and hops between crossing MTs; (ii) run-lengths; (iii) run-times; and (iv) off-axis steps. To elucidate the factors that control NP motility, we modeled the active transport of the NP by including several competing processes, such as the binding to and the unbinding from the MT surface and the stepping kinetics of individual dynein motor proteins. These processes take place on the curved 2D microtubule surface and are influenced by both the elastic coupling between the motors (via the spacer polymers) and the excluded-volume interaction between motors. We used Monte-Carlo (MC) simulations to describe these dynamics. We highlight the unique features of the multi-motor carried NP by comparing between the experiment and theorybased simulations, which mutually support each other. Moreover, this comparison between the experimental and theoretical results provides insights into the mechanism underlying NP motion and, in particular, into the role of multiple motor action during transport.

Preparation of Hela Cell Extracts
Hela CEs were prepared according to Fu et al. [54]. Briefly, Hela cells were grown in a DMEM medium supplemented with 10% fetal bovine serum (FBS) and 1% of penicillinstreptomycin at 37 • C and 5% CO 2 . The cells were detached by a trypsin-EDTA solution, washed with phosphate-buffered saline (PBS) (0.01 phosphate buffer, 0.0027M KCl, 0.137 M NaCl, pH 7.4), and pelleted for 5 min at room temperature (RT) and 500 g. The cells were then incubated on ice with lysis buffer (12 mM Pipes, 2 mM MgCl 2 , 1 mM EGTA, 0.4% Triton X-100, 10% protease inhibitor cocktail, pH 6.8) for 15 min. Finally, the lysates were cleared by centrifugation at 2700× g and then at 100,000× g, at 4 • C, for 10 min for each lysate. The CEs were diluted 10× with a lysis buffer without Triton and protease inhibitors, and the total protein mass concentration was determined by a Bradford assay, henceforth denoted by [CE]. Sucrose was added to the extracts (10% in mass), which were then flashfrozen in liquid nitrogen and stored at −80 • C. For the motility experiments, the extracts were diluted to the same total protein mass concentration of 3.4 mg/mL (systems I-IV), except for system V for which a [CE] of 6.8 mg/mL was used. The extracts were used within four weeks.

Nuclear Localization Signal (NLS)
The NLS sequence (PKKKRKVED) originated from an SV40 T large antigen [55]. We used N-terminally bromine (Br-Ac) modified NLS. The bromine group was followed by a GGGG sequence ("raft"). To generate a fluorescently labeled NLS peptide, we covalently attached Tetramethylrhodamine (TAMRA) to the NLS raft. Except for the UV-Vis absorption experiments, in which we used TAMRA-NLS, all experiments were performed with the non-labeled NLS.

Preparation of NPs
Bare, green-fluorescent microspheres (Bangs Labs, Indianapolis, IN, USA or Invitrogen, Waltham, MA, USA) were used to prepare the NPs. The total surface area of the bare NPs was constant throughout all the experiments, which guaranteed that the mean anchoring distance between PEG-NLS molecules was independent of the NP diameter. As a starting point, one can also use streptavidin-coated NPs (Bangs Labs, IN, USA) instead of bare microspheres. NP preparation started with the incubation in 10 µM Neutravidin (31,000, Thermo Fisher, Waltham, MA, USA) for 20 min at 25 • C. Excess Neutravidin was separated via dialysis. Next, the NPs were incubated for 20 min at 25 • C in 10 mg/mL BSA to block the remaining uncoated regions on the NPs surface. Excess BSA was separated by centrifugation (6800× g for 8 min at 25 • C). All subsequent centrifugation steps were performed under the same conditions. The NPs were then incubated with 1 mM Biotin-PEG-thiol-5kDa (PG2-BNTH-5, NANOCS, Boston, MA, USA see Figure S3) for 30 min at 25 • C. Excess Biotin-PEG-thiol was separated by two cycles of centrifugation. Then, NLS peptides were covalently attached to the NPs via their bromine group, which reacts with the thiol group on the PEG molecule. The reaction was performed at 25 • C for 20 min. Excess NLS was separated by two cycles of centrifugation. Finally, the NPs were incubated with Hela CEs at 30 • C for 20 min, allowing for the recruitment of α-and β-importins, dynactin, and mammalian dynein motors to the surface of the NPs. The NPs were separated from excess CE by centrifugation, washed with BRB80 (80 mM PIPES pH 6.8, 1 mM MgCl 2 , 1 mM EGTA) supplemented with 1 mM Mg-ATP, and centrifuged again. Finally, the NPs were resuspended in 60 µL BRB80 supplemented with 1 mM Mg-ATP and stored on ice until used. In the motility assay experiments, two different batches of CEs were used: the same CE was used in systems I and II, and another CE was used in systems III, IV, and V. To ensure full activity of the motility protein machinery, the NPs were used within 3 h.

Cryo-Electron Microscopy (Cryo-TEM)
Vitrified specimens were prepared according to a standard procedure [56]. Briefly, 2.5 µL drops of an NP solution were applied to a copper grid coated with a perforated lacy carbon 300 mesh (Ted Pella Inc., Redding, CA, USA) and blotted with a filter paper to form a thin liquid film. The blotted samples were immediately plunged into liquid ethane at its freezing point (−183 • C) using an automatic plunge freezer (EM GP, Leica Microsystems GmbH,Vienna, Austria) and transferred into liquid nitrogen for storage. The samples were analyzed using an FEI Tecnai 12 G2 TEM at 120 kV with a Gatan cryo-holder, maintained at −180 • C. Images were recorded on a CCD camera (Gatan manufacturer, Pleasanton, CA, USA) at low-dose conditions.

Dynamic Light Scattering and ζ-Potential Experiments
DLS and ζ-potential measurements were used to extract the NP hydrodynamic diameter (D h ) and charge at the various decoration steps. For these experiments, we used bare NPs with a 0.196 µm diameter. The measurements were performed on a Malvern NanoZS instrument (ZN-NanoSizer, Malvern, UK) operating with a 2 mW HeNe laser at a wavelength of 632.8 nm. The detection angle was 173 • and 17 • for the DLS and ζ-potential measurements, respectively. All measurements were conducted at 25 ± 0.05 • C; for the analysis of D h and ζ-potential, the viscosity was taken to be that of water (0.8872 cP). The intensity size (hydrodynamic diameter) distribution was extracted from the intensity auto-correlation function, which was calculated using an ALV/LSE 5003 correlator over a time window of 30 s (10 runs of 3 s) using the software CONTIN. Each measurement was repeated three times. The value of D h was averaged over three independent experiments; error bars indicate the standard deviations for these three experiments. For ζ-potential measurements, the solution was transferred to a U-tube cuvette (DTS1070, Malvern, UK), which was operated in automatic mode. The electrophoretic mobility of the NPs was measured, from which the ζ-potential value was determined by applying the Henry equation [57]. The ζ-potential values were averaged over three independent experiments with 30 runs per experiment; error bars indicate the standard deviations for these three experiments.

UV-Vis Absorption Experiments
UV-Vis absorption experiments were used to determine the mean number of bound PEG-NLS molecules per NP, N , from which the mean anchoring distance between adjacent PEG-NLS molecules, ξ * , was deduced. These experiments were conducted at a wavelength of 553 nm, using a fluorescently labeled NLS (TAMRA-NLS). Biotin-PEG-thiolgrafted NPs were incubated with increasing amounts of TAMRA-NLS and the absorption of the remaining TAMRA-NLS molecules that did not bind to the NPs (i.e., supernatant) was measured. The number of TAMRA-NLS molecules that did bind to the NPs was then deduced (see details in Supplementary Material Section S1). An extinction coefficient ε 553 = 0.0639 1 µM cm , which we determined experimentally (not shown), was used for that purpose. The values of N and ξ * were averaged over three independent experiments; error bars indicate the standard deviations for these three experiments.

Western Blotting
Western blot (WB) was used to confirm the recruitment of importins and dynein motors by the PEG-NLS coated NPs. These experiments employed NPs of 0.196 µm in diameter, [NLS] of 7 µM, and [CE] of 3.4 mg/mL. The NPs were prepared as detailed above, then pelleted at 6800× g for 10 min at 25 • C, resuspended in a 1× Laemmli Sample Buffer, and boiled for 5 min to promote the detachment of the bound proteins. The proteins were separated by electrophoresis using a 12% agarose gel and transferred to a nitrocellulose membrane. The membrane was incubated for 1 h in a blocking buffer of PBST (PBS supplemented with 0.1 v/v % Tween) and 10% (v/w) dry skim milk (Sigma-Aldrich, St. Louis, MO, USA). The membrane was washed three times with PBST for 5 min and then incubated for 1 h at 25 • C with anti-dynein (sc-13524), anti-karyoprotein α2 (sc-55538), anti-dynactin (sc-135890), and anti-kinesin-1 (sc-133184) (Santa Cruz Biotechnology, Dallas, TX, USA). All primary antibodies were diluted 1:200 (v/v) in a blocking buffer prior to use, except anti-dynactin, which was diluted to 1:50 (v/v). Next, the membrane was washed three times with PBST and incubated for 1 h at 25 • C with an anti-mouse HRP conjugated with a secondary antibody (sc-2005, Santa Cruz Biotechnology, Dallas, TX, USA) diluted 1:10,000 (v/v) in PBST supplemented with 0.5% (v/w) skim milk. To finalize the procedure, the membrane was washed three times with PBST and incubated with an ECL Western blot reagent (1,705,060, Bio-Rad, Hercules, CA, USA) for 5 min in the dark. Images were collected by chemiluminescence, using the Fusion FX imaging system (Vilber Lourmat, Collégien, France).

Motility Assay Experiments
Chamber preparation. Flow cells were prepared using a glass slide and a glass coverslip (washed with deionized water, 70% EtOH, and dried using nitrogen gas), and two stripes of warm Parafilm were placed in between them to form a channel with a width of a few millimeters. The chamber was incubated with one volume of 0.1 mg/mL Biotinylated Casein (prepared by biotinylation of k-Casein (Sigma-Aldrich, St. Louis, MO, USA) using EZ-Link (21,336, Thermo Fisher, Waltham, MA, USA) for 5 min at 23 • C and washed with three volumes of BRB80. The chamber was then incubated with one volume of 1 mg/mL Neutravidine (31,000, Thermo Fisher, Waltham, MA, USA) for 5 min at 23 • C and washed with three volumes of BRB80.

Preparation of Microtubules
Tubulin from a porcine brain, purified via three polymerization/depolymerization cycles [58], was flash-frozen in liquid nitrogen and kept at −80 • C until used. Biotinfluorescently labeled microtubules were prepared according to the protocol described in the work of Gell et al. [59]. A tubulin mix, containing 51 µM tubulin, 6 µM Biotin tubulin (T333P, Cytoskeleton, Denver, CO, USA), 3 µM Rhodamine tubulin (TL590M, Cytoskeleton, Denver, CO, USA), and 3 mM GMPCPP (NU-405S, Jena Bioscience, Jena, Germany), was mixed on ice, divided, flash-frozen in liquid nitrogen, and kept at −80 • C until used. Prior to the mix preparation, the tubulin was thawed and kept on ice for 5 min and then centrifuged for 10 min at 126,000× g, 4 • C. The supernatant was kept on ice and its concentration was determined by absorbance at 280 nm (ε 280 = 115, 000 1 M cm ). For microtubule (MT) assembly, an aliquot of the tubulin mix was thawed, diluted with BRB80 to a final concentration of 4 µM, and incubated for 3 h at 37 • C to promote MT assembly. For motility assay experiments, the MTs were diluted to 40 nM with a warm (37 • C) wash buffer (WshB) (BRB80, containing 0.02 mM paclitaxel, 15 mM glucose, 0.1 mg/mL glucose oxidase, 0.02 mg/mL catalase, and 50 mM DTT) and used on the same day. The MTs were kept at 23 • C until used.

Formation of Marked-End Microtubules
Marked-end microtubules were prepared by using bright fluorescent MT seeds, which served as the nucleation sites for the polymerization of dim fluorescent, N-Ethylmaleimide (NEM)-modified tubulin [60,61]. The NEM-modified tubulin was used to inhibit the MT seed minus-end assembly, and it was prepared by mixing 100 µM tubulin with 1 mM NEM (Sigma-Aldrich, St. Louis, MO, USA) and 0.5 mM GMPCPP in BRB80 and placing the solution on ice for 10 min. The reaction was quenched by adding 8 mM β-mercaptoethanol. The NEM-tubulin mix was incubated on ice for an additional 10 min, and then a tubulin mix, containing 3.2 µM tubulin, 4 µM Rhodamine tubulin, 0.8 µM Biotin tubulin, and 1 mM GMPCPP in BRB80, was incubated for 15 min at 37 • C to produce short "bright" MT seeds. One volume of bright MT seeds mixed with seven volumes of a "dim" tubulin mix (containing 3.6 µM tubulin, 0.4 µM Rhodamine tubulin, 0.8 µM Biotin tubulin, 3.2 µM NEM-modified tubulin, and 1 mM GMPCPP in BRB80) was incubated at 37 • C for 1 h to promote polymerization. The marked-end MTs were diluted to 40 nM with warm (37 • C) WshB and used on the same day. The MTs were kept at 23 • C until used.

Motility Assay Experiments
For these experiments, 40 nM of MTs were introduced in the flow cell and incubated for 10 min at 23 • C. The chamber was washed with three volumes of WshB to remove unbound MTs. Prior to the assay, the NPs were pelleted by centrifugation at 6800× g for 8 min at 4 • C, resuspended in 60 µL motility buffer (BRB80 containing 0.02 mM Paclitaxel, 15 mM glucose, 0.1 mg/mL glucose oxidase, 0.02 mg/mL catalase, 50 mM DTT, 0.1% methyl cellulose 4000 cP, and 10 mM Mg-ATP), and incubated on ice for 5 min. Prior to insertion in the flow cell, the NPs were brought to 23 • C. The motility assays were performed with bare NPs of 20 nm in radius.

NP Imaging Using Total Internal Reflection Fluorescence Microscope (TIRFM)
We followed the motion of the NPs by TIRF microscopy using an automated TIRF multicolor system integrated on an inverted Leica DM6000B microscope (Leica Microsystems GmbH, Wetzlar, Germany). Imaging was performed using a HCX Plan Apo 1.47 N.A. 100× TIRFM oil immersion objective and a triple-band Laser Line Violet Blue Green (VBG) emission filter (Leica Microsystems GmbH, Wetzlar, Germany). The samples were excited by TIRF illumination using 10 mW 488 nm and 20 mW 561 nm solid-state lasers. The images were captured using an Andor DU-897 EMCCD camera (Oxford Instruments, Abingdon, UK). The imaging system produced images of the NPs with apparent radii of 3 pixels (pixel size is 230 nm), which allowed the extraction of the center-of-mass coordi-nates of the NPs, NP(X, Y), and the center-of-mass coordinates of the MTs, (X, Y) MT , with a sub-pixel resolution and a 10 nm scale precision (see below).

Particle Tracking and Data Analysis
We used the interactive data language (IDL) multi-particle tracking method [62,63], implemented in MATLAB, to automatically detect distinct NP trajectories. Using this method, we extracted the 2D center-of-mass coordinates (X(t), Y(t)) per time point t (frame) of the individual NPs by fitting the NP fluorescence profiles along the X and Y directions using a Gaussian function, from which we determined the individual NP trajectory over time, NP(X, Y) t . We conducted an initial (automatic) filtering, in which we selected trajectories that included at least three steps. Next, we performed a more delicate (manual) filtering, in which we excluded NPs or NP trajectories that fell under one or more of the following categories: (i) the NP shape was not symmetric or too large, which implies NP aggregation; (ii) the NP moved in a region where the density of MTs was too dense, and thus, it was impossible to discern between individual MT tracks; and (iii) the NP moved only for a short period and then got stuck along the MT track for a long period. For each NP, the overall run-time and run-length were calculated from the accumulated travel time and distance, respectively. In some cases, where the NPs moved continuously between crossing MT tracks, we extracted the trajectory on each MT track separately. For these NPs, the overall run-time and run-length correspond to the accumulated traveled time/distance of all individual trajectories.
The NP velocity v was calculated by taking the center-of-mass position of the NP at times t, (X(t), Y(t)) NP and t + ∆t (X(t + ∆t), Y(t + ∆t)) NP , and dividing it by ∆t, that is , where X and Y refer to lab frame Cartesian coordinates. The longitudinal velocity, v x , was calculated by projecting the NP velocity v onto the MT direction (by definition, a positive direction points towards the MT minus-end), i.e., the unit vector of the MT (along its long axis),û MT The indexes i and f refer to the initial and final time points of the NP trajectory, and (X i , Y i ) MT and X f , Y f MT are the corresponding MT center-of-mass coordinates, which are extracted by fitting the cross-section intensity profiles of the MTs using a Gaussian function. Thus, the Positive v x corresponds to an NP step ∆x in the direction of the MT minus-end, whereas negative values of v x refer to an NP step towards the MT plus-end Figures 2C and S4. Similarly, the transverse motion, v y , was determined by projecting the NP velocity v onto the normal unit vector of the MT,n MT . Thus, the transverse velocity v y = v·n MT and the corresponding transverse step is given by ∆y = v y ·∆t. Note that, for each MT or MT unit vectorsû MT , there are two opposite normal unit vectors that can be assigned,n 1,MT andn 2,MT , which account for positive and negative NP transverse motions, respectively ( Figure S4). In this paper, we defined the right-handed transverse motion as positive and the left-handed transverse motion as negative (see Figure 2C for the definition of spatial orientation).
Next, we filter unrealistic velocity values as follows: First, we filter out absolute velocity values, υ = υ 2 x + υ 2 y , exceeding 4000 nm s , which represents the very rare-and not very plausible-case where the dynein steps 100 consecutive 40 nm sized steps. Second, we filter out unrealistic transverse steps. To this end, we set a maximum size of a transverse step, ∆y = 2 × R orb = 2 × 153 nm, which corresponds to the scenario where the NP moves half a circle along the MT perimeter, from left to right or vice versa. A transverse step that is larger than the said limit is omitted.
We use the mean value of transverse steps, ∆y , to extract the mean angular step, ∆φ , and its error, ∆φ e , and the mean angular velocity, ω , and its error, ω e , as follows: where SEM(∆y) is the standard error of the mean (SEM) value of ∆y, and R orb is the approximated distance between the MT and NP center ( Figures 2C and 6A, which reads: where SEM(υ x ) is the standard error of the mean value of υ x . Similarly, we can estimate the mean angular velocities for right-( ω |ω > 0) and left-( ω |ω < 0 ) handed motions and the mean angular velocity for separated minus-end directed ( ω |υ x > 0) and plus-end directed ( ω |υ x < 0 ) motions (see Supplementary Material Section S5).

Results
Below, our analysis of the experimental and model-simulation results will address the following parameters: (i) the mean anchoring distance between neighboring PEG-NLS molecules ξ * ; (ii) the mean anchoring distance between the PEG-NLS-αβ-dynein (αβ refers to the α-and β-importin complex) ξ; (iii) the number of motors, N m , that are anchored to the surface of the NP; and (iv) the number of MT-bound motors, M B , that participate in the motion at a given time. The mean number of NP-bound motors is related to ξ via N m = 4πR 2 ξ 2 , where R is the radius of the NP.

NP Preparation and Analysis
The NP synthesis process entails several consecutive steps where, at each step, a single component is added ( Figure 1A). Bare NPs, saturated with Neutravidin, were conjugated with Biotin-PEG-thiol spacers and then incubated for a short time in a solution containing an SV40 T large antigen NLS peptide at various concentrations.
The resulting NPs are characterized (Supplementary Material Section S1) by using cryo-TEM ( Figure S1), DLS, ζ-potential, and UV-Vis absorption isotherm, yielding the surface density and the mean number of grafted Biotin-PEG-thiols, N , which are endconjugated by NLS (PEG-NLS) ( Figure S2). The latter can be transformed into a mean anchoring distance between neighboring PEG-NLS molecules, ξ * . The resulting anchoring distance, ξ * , is depicted in Figure 1B against [NLS], showing that the incubation of NPs in higher [NLS] solutions leads to a shorter distance between adjacent PEG-NLSs until it reaches saturation above [NLS] = 3 µM. Note that the theoretical value of the (freepolymer) gyration radius is R g = 2.27 nm (for a Biotin-PEG-thiol of M w = 5 kDa; see Supplementary Material Section S2 [69][70][71][72][73][74][75]). Thus, since the diameter of Neutravidin is about 5 nm [76], and since the Neutravidins are closely packed on the NP surface, we may conclude that the anchored Biotin-PEG-thiol molecules are effectively in the socalled "mushroom regime" [77], ξ * > R g , wherein the polymer chains appear as isolated "mushrooms" (i.e., random coils of radius ∼ R g ) on top of the NP surface.
Next, NPs are incubated in Hela CEs, thus, allowing the recruitments of α-and β-importins, dynactin, and mammalian dynein association [54]; the recruitments of the importins, dynactin (which is required for the proper functioning of dynein [78,79]), and mammalian dynein are verified by WB (see Materials and Methods). Figure 1Ca shows the WB image after exposure to the CE, demonstrating antibody-specific binding to α2-importin. Figure 1Cb shows a WB image demonstrating dynein recruitment, and Figure 1Cc shows a WB image demonstrating dynactin recruitment. Group 1 refers to NPs coated with Biotin-PEG-thiol incubated in CE; group 2 refers to PEG-NLS coated NPs, which were also incubated in CE; and group 3 refers to a CE without NPs. Note that although β-importin binding was not verified directly, dynein recruitment requires β-importin binding, suggesting that, in the absence of NLS, the dynein machinery is not recruited to the NP surface.
An important question regarding our NP construct is whether it allows the non-specific binding of kinesin motors, which will influence the motility characteristics. To discard this possibility, we performed WB for kinesin-1 binding to the NP ( Figure 1D). The WB image does not show kinesin-1 binding to the PEG-NLS coated NPs (although it remains possible that kinesin-1 binding occurs below the WB detection threshold). This finding is expected since the bare NP is covered with a first layer of Neutravidin, a second layer of BSA, and a third layer of Biotin-PEG-thiol, which, together, are expected to efficiently passivate the NP surface.

NP Motility Assays
The NPs were washed from the excess cell extract, incubated in an ATP solution, and then injected into a flow cell in which MTs were adsorbed and immobilized on a glass surface. We investigated the different modes of motion of the NPs under different [NLS] but identical [CE], as defined in Table 1. Thus, we investigated four systems, wherein extracts from an identical batch were used in systems I and II and extracts from another identical batch were used in systems III and IV. Since individual protein concentrations may vary from one CE to another (although the total mass protein concentration, [CE], that we used was identical; see Material and Methods), for the sake of prudence, we compared the results of the motility assays only within each pair: I vs. II and III vs. IV. We also examined another system, referred to as system V, whose CE batch was identical to that of systems III and IV but whose [CE] was twice as high; however, we do not discuss this system below since its motility is strongly reduced, presumably due to molecular crowding effects (see Supplementary Material Section S3, Table S1). To follow the motion of the NPs, we used TIRFM (e.g., see Movie 1 (system II) and Movie 2 (system III)). First, we used marked minus-end MTs to confirm the expected correlation between the polarity of the MT and the motional direction of the NP, as dictated by the dynein bias to move towards the minus-end ( Figure 3A,B). Next, we used a standard particle tracking algorithm to extract the center-of-mass (CM) position of the NP with a 10 nm-scale precision (see Materials and Methods), which yielded individual time-dependent trajectories. Notably, although several NPs can potentially move on the same MT track simultaneously, we mostly observed only a single NP at a time. Notwithstanding, to prevent possible interference from nearby NPs (e.g., obstruction), we purposely chose to analyze only NPs that moved alone on a single MT. For each CE batch, we tested how [NLS] affects NP motion. Table 1 summarizes the estimated values of the mean anchoring distance * of PEG-NLS for each of the four systems, which we deduced from the mean number of grafted PEG-NLS, 〈 〉 (see Figure S2 in Supplementary Material Section 1). Note that * serves as an estimated lower bound for the mean anchoring distance of the PEG-NLS--dynein, such that should approach * (from above) as the concentration of associated protein concentration (and -importins, dynein, dynactin, etc.) in the CE increases ( was not directly measured).
Moreover, given the expected sub-millimolar concentrations of dynein [82], we assume the system is much below the saturation of the dynein binding "isotherm", such that the increase of CE proteins (and, accordingly, the dynein bulk concentration) is likely to increase dynein surface concentration, resulting in a decrease of . Thus, we expect systems I and II, which belong to the same CE batch and have the same CE , will demonstrate the same value of / * or, equivalently, the same value 〈 〉 〈 〉 = * ; a similar argument applies to systems III and IV. The value of * determines the mean number of PEG-NLS per NP via 〈 〉 = * , which implies that 〈 〉 varies between 5 and 37 for NLS that varies between 0.025 and 0.3 μM (Table 1 and Supplementary Material Section 1). To analyze their motion, we differentiated the trajectories of the NPs into plus-end and minus-end directed motion, with respect to the movement along the long-axis of the MT (henceforth defined as "longitudinal motion"), see Figure 2C. We detected three major modes of motion ( Figure 3 and Table 1): (i) NPs that moved continuously towards the MT minus-end ( Figure 3A,B); (ii) NPs that reached the end of the MT after exhibiting continuous minus-end directed motion, followed by a single backward step (i.e., towards the MT plus-end) and detachment ( Figure 3C,D). In some cases (not shown), the NP remained immobile for some time before detaching; and (iii) NPs that traversed between crossing MTs ( Figure 3E,F, Movie 3). Importantly, trajectories that involve traversing between crossing MT tracks have been previously observed by Ross et al. in single-motor assays [80,81]; however, the simultaneous binding of two dynein motors to the NP should enhance the frequency of such traversing events with the help of an intermediate configuration where one motor binds to one track and the second motor binds to the other track.

Longitudinal velocity
For each CE batch, we tested how [NLS] affects NP motion. Table 1 summarizes the estimated values of the mean anchoring distance ξ * of PEG-NLS for each of the four systems, which we deduced from the mean number of grafted PEG-NLS, N (see Figure S2 in Supplementary Material Section S1). Note that ξ * serves as an estimated lower bound for the mean anchoring distance ξ of the PEG-NLS-αβ-dynein, such that ξ should approach ξ * (from above) as the concentration of associated protein concentration (α-and β-importins, dynein, dynactin, etc.) in the CE increases (ξ was not directly mea-sured). Moreover, given the expected sub-millimolar concentrations of dynein [82], we assume the system is much below the saturation of the dynein binding "isotherm", such that the increase of CE proteins (and, accordingly, the dynein bulk concentration) is likely to increase dynein surface concentration, resulting in a decrease of ξ. Thus, we expect systems I and II, which belong to the same CE batch and have the same [CE], will demonstrate the same value of ξ/ξ * or, equivalently, the same value ; a similar argument applies to systems III and IV. The value of ξ * determines the mean number of PEG-NLS per NP via N = 4πR 2 ξ * 2 , which implies that N varies between 5 and 37 for [NLS] that varies between 0.025 and 0.3 µM (Table 1 and Supplementary Material Section S1).

Longitudinal Velocity
The mean properties of the NP trajectories are shown in Figure 4 and are summarized in Table S1 for each system. Except for system I, we found that the longitudinal velocities of the NPs significantly decreased compared with those of a single motor. Comparing each of the pairs (I to II and III to IV), we can deduce that the number of NP-attached motors increased such that N m was higher in system II than in system I and in system IV than in system III. Assuming that the mean number of motors participating in the transport M B increases with N m (as shown by our theoretical predictions in Section 3.4, below) suggests that NP velocities decrease with increasing M B , presumably due to inter-motor interactions. We emphasize that the regime of transport corresponds to vanishing drag due to the nanometer size of the NPs, implying that a load-sharing effect should not be present whatsoever; such a reduction of the velocity along with the increase in the number of participating motors is commonly observed in standard motility assays of MTs moving on kinesin/dynein-decorated glass surface [83,84]. Examining the average NP run-time, τ p , and run-length (along the MT symmetry axis), λ , for the different systems ( Figure 4B and Table S1) revealed that τ p increases with increasing [NLS] or [CE], leading to a decrease in ξ or an increase of N m , consistent with the anticipated increase in the mean number of participating motors, M B (c.f. Section 3.4). Note, however, that since λ = υ x × τ p , the dependence of the run-length on [NLS] and [CE] is non-monotonous and shows a maximum for system IV (Table S1).
Notably, about 90% of the steps in our experiments were minus-end directed (υ x > 0); Tables S1 and S2, and Figure S4 show the distributions of the temporal absolute velocity, υ, longitudinal velocity, υ x , and transverse velocity, υ y , respectively. To further characterize the NP motion, as shown in Figure 4C,D and in Table S1, we extracted the mean values of the longitudinal velocity (minus-end directed motion, υ x > 0, and plus-end directed motion, υ x < 0) and of the transverse velocity (left, υ y < 0, and right, υ y > 0). As [NLS] increases, which implies a decrease in ξ and an increase in the number of NP bound motors, N m = 4πR 2 ξ 2 , the corresponding absolute and longitudinal mean velocities decrease. Conversely, the probability of either the plus-end or the minus-end directed motions do not vary greatly (Table S1). Notably, as the [NLS] increases (i.e., moving from system I to II and from system III to IV), the plus-end directed longitudinal velocity, υ x | υ x < 0 , decreases to a much greater extent than the minus-end directed longitudinal velocity υ x | υ x > 0 ( Figure 4C). Notably, about 90% of the steps in our experiments were minus-end directed ( > 0); Tables S1,S2, and Figure S4 show the distributions of the temporal absolute velocity, , longitudinal velocity, , and transverse velocity, , respectively. To further characterize the NP motion, as shown in Figure 4C,D and in Table S1, we extracted the mean values of  Table S1). (A) Mean longitudinal, υ x , transverse, υ y , and absolute, υ = υ 2 x + υ 2 y velocities. A minus-end directed motion corresponds to υ x > 0 and a righthanded motion corresponds to υ y > 0 (see Figure 2C for spatial orientation). (B) Run-times and run-lengths: the run-time is measured from the time the NP binds to the MT until the time it unbinds from it; the run-length is the accumulated distance traveled along the MT symmetry axis; υ x > 0 and υ x < 0 define minus-end and plus-end directed motion, respectively. The observed variability in the properties of motility between system II and system III is relatively small; indeed, as can be clearly seen in Figure 4B, the SEMs of τ p and λ overlap between systems II and III. Thus, the relatively small differences in the properties of motility between the two systems can be rationalized by the expected small difference in the protein composition of the CE.
We next argue that NP motility properties are only sensitive to the number of NPassociated motors, N m . Consider the two observables υ x | υ x > 0 and τ p . In Figure 4E, we combine the data from the four different systems, I-IV, which vary either in PEG-NLS surface coverage, CE preparation, or both. In the main panel, we plot, for each NP, its υ x | υ x > 0 against τ p as a separate data point. In the inset, we average the overall NPs to obtain a single functional dependence. It is clearly visible that there is some overlap of data between the different systems. A plausible explanation is that N m is effectively identical for overlapping NPs, regardless of their associated system-a hypothesis that will be further tested against the model simulations, as presented below.

Angular/Transverse Motion
As discussed above, while most studies of dynein motors assume a purely longitudinal motion, more recent studies discovered rich transverse dynamics [37,38,52,85]. Therefore, we analyzed the NP trajectories for transverse motion as well. We first calculated the estimated transverse velocity, υ y , of the NPs ( Figure 4D) and found that it was significantly lower than the longitudinal velocity.
To gain further insight into the transverse dynamics, we mapped the transverse velocity into an angular velocity, ω = ∆φ ∆t (where ∆φ is the angular increment of a single time step, ∆t), as shown in Figure 5, Table S3 (see Materials and Methods for details). We assumed an ideal angular motion with constant orbital radius, R orb = 153 nm ( Figure 2C), which we estimated by considering the molecular dimensions of the NP and its associated ligands (see Materials and Methods and Supplementary Material Section S4). The angular motion, together with the longitudinal motion, imply that if the MT had been elevated above the surface, the NP would have performed a helical motion around the symmetry axis of the MT. Indeed, one of the NPs showed such a helical motion, which was made possible by the (unintended) elevation of the specific MT on which it moved (see Movie 2, pink arrow). Considering the mean angular velocity (associated with Figure 5A) as a representative number, we estimated the mean pitch size of the assumed helical trajectory by using H = 2π| v x |/| ω | (see Materials and Methods and Supplementary Material Section S5), which yielded the results shown in Figure 5B and Table S3. The trajectories of a single NP taken from the theoretical model simulations (presented in the next section) show large angular motion "fluctuations" of varying magnitudes (Movies 4-7 and Figure 6). To analyze these fluctuations, we split the angular motion into right-and left-handed motions (Table S3; see Materials and Methods for details), yielding the angular velocities shown in Figure S6A. We found that, indeed, the right-(〈 〉| > 0) and left-(〈 〉| < 0) handed motions occurred at similar angular velocities in all examined systems. In the Supplementary Material, we provide the correlations between the positive and negative longitudinal steps, on the one hand, and the positive and negative angular steps, on the other hand ( Figure S6B and Table S3). However, due to the large inherent fluctuations associated with the angular motion, we cannot come to any conclusions about such correlations. The possibility of such correlations, between longitudinal and angular steps, is further examined by the model simulations, as detailed below. The trajectories of a single NP taken from the theoretical model simulations (presented in the next section) show large angular motion "fluctuations" of varying magnitudes (Movies 4-7 and Figure 6). To analyze these fluctuations, we split the angular motion into right-and left-handed motions (Table S3; see Materials and Methods for details), yielding the angular velocities shown in Figure S6A. We found that, indeed, the right-( ω |ω > 0) and left-( ω |ω < 0 ) handed motions occurred at similar angular velocities in all examined systems. In the Supplementary Material, we provide the correlations between the positive and negative longitudinal steps, on the one hand, and the positive and negative angular steps, on the other hand ( Figure S6B and Table S3). However, due to the large inherent fluctuations associated with the angular motion, we cannot come to any conclusions about such correlations. The possibility of such correlations, between longitudinal and angular steps, is further examined by the model simulations, as detailed below.

Model and Simulations
To explain the experimental findings and to gain profound knowledge regarding NP motion and behavior, we constructed a model for the active transport of the NP. Since previous work on multi-motor complexes addressed a single-motor motion on a 1D microtubule track [32,33,36,42,51], it cannot describe the motion of the NP around the 2D MT surface. The single-motor stepping was recently revisited to describe the 2D stepping of a single yeast dynein on the microtubule surface; this stepping model correctly accounts for the measured longitudinal step-size distribution and, moreover, predicts a broad angular distribution of steps (of a single motor) with a small right-handed bias [40]. More recently, Elshenawy et al. [38] studied the single mammalian dynein stepping kinetics, providing the longitudinal and transverse step-size distributions. Using the raw data of Yildiz and co-workers [38,53], we adjusted our yeast dynein model [40] to describe a single mammalian dynein (see Supplementary Material Section 6.1 and Table S4).
Our NP multi-motor model is described in Supplementary Material Section 6 and Figure S7. The model entails an NP on which polymers with an identical contour length and a fixed density, corresponding to a prescribed mean spacing * , are grafted at random positions. Each free polymer-end is assigned with a dynein motor, such that = * , by definition. Thus, = √4 / , where is the NP radius and is the number of motors that are anchored to the NP surface. Dynein is assumed to bind to and unbind from the MT surface, processes that compete with its stepping kinetics (see Supplementary Material Section 6.2.3). Note that the term "unbinding" refers to the full detachment of the

Model and Simulations
To explain the experimental findings and to gain profound knowledge regarding NP motion and behavior, we constructed a model for the active transport of the NP. Since previous work on multi-motor complexes addressed a single-motor motion on a 1D microtubule track [32,33,36,42,51], it cannot describe the motion of the NP around the 2D MT surface. The single-motor stepping was recently revisited to describe the 2D stepping of a single yeast dynein on the microtubule surface; this stepping model correctly accounts for the measured longitudinal step-size distribution and, moreover, predicts a broad angular distribution of steps (of a single motor) with a small right-handed bias [40]. More recently, Elshenawy et al. [38] studied the single mammalian dynein stepping kinetics, providing the longitudinal and transverse step-size distributions. Using the raw data of Yildiz and co-workers [38,53], we adjusted our yeast dynein model [40] to describe a single mammalian dynein (see Supplementary Material Section S6.1 and Table S4).
Our NP multi-motor model is described in Supplementary Material Section S6 and Figure S7. The model entails an NP on which polymers with an identical contour length and a fixed density, corresponding to a prescribed mean spacing ξ * , are grafted at random positions. Each free polymer-end is assigned with a dynein motor, such that ξ = ξ * , by definition. Thus, ξ = √ 4π R/ √ N m , where R is the NP radius and N m is the number of motors that are anchored to the NP surface. Dynein is assumed to bind to and unbind from the MT surface, processes that compete with its stepping kinetics (see Supplementary Material Section S6.2.3). Note that the term "unbinding" refers to the full detachment of the motor from the MT, which implies that its two microtubule binding domains (MTBDs) become disconnected from it, such that the inverse of the single motor unbinding rate defines its processivity time. Likewise, "binding" refers to the attachment of a motor to the MT surface from the bulk solution. The motion of the motors on the 2D curved microtubule surface is influenced by the elastic coupling between motors (via the spacer polymers)-a coupling that influences both the motor step vector distribution and the binding-unbinding kinetics. In addition, we account for the excluded-volume interaction between motors that prevents them from stepping over each other; however, we do not account for the excluded volume interaction between different spacer polymers, which, we believe, is negligible (Supplementary Material Section S6.2.1 and Figures S8 and S9 [86][87][88]). Consistent with the nanometer scale of the NP, the drag force on the NP is negligible (see Section 1). We used Monte-Carlo (MC) simulations to describe the different competing processes, while readjusting the NP center-of-mass and rotational angle after each MC step to achieve mechanical force balance-processes that are orders of magnitude faster than the binding/unbinding and stepping processes (Supplementary Material Section S6.2.5). The model allows the investigation of the different microscopic internal states between which the NP fluctuates in time, in particular the number of MT-bound motors, M B , that participate in the motion at a given instance.

Model Simulation Results
Results were obtained for various NP configurations, which are characterized by two parameters: R and N m . Below, we discuss an NP with radius R = 20 nm and a number of motors ranging between N m = 1 to N m = 13. For each (R = 20 nm, N m ), we ran the simulations for a few hundred (identical) particles to obtain high statistical accuracy. Note that due to the very small NP drag coefficient (Supplementary Material Section S6.2.5), we used N m = 1 as a representative of the single-motor case. Figures 7 and 9, and Supplementary Material Section S7, describe the motility characteristics, on a time scale of 0.27 s, of each NP configuration, while Figures 8-10 and Supplementary Material Section S8 complement these with analyses on the timescale of a single MC step ∼ 1 ms, which is two orders of magnitude shorter. Figure 6 and Movies 4-7 show selected trajectories of single-motor, three-motor, seven-motor, and thirteen-motor NPs, which demonstrate rich, fluctuating, helical motions. Below, we investigate the influence of the number of motors on such behavior.   Table S5 see Figure S11 and Table S9 for a Δ = MC step time. (B) Mean (±SEM) run-length (black dots) and run-time (empty circles), using the data provided in Table S10. (C) Mean (±SEM) instantaneous longitudinal velocity, 〈 〉, separated for minus-end directed, > 0, and plus-end directed, < 0, directions; data are provided in Table S6 Table  S5 see Figure S11 and Table S9 for a ∆t = MC step time. (B) Mean (±SEM ) run-length (black dots) and run-time (empty circles), using the data provided in Table S10. (C) Mean (±SEM ) instantaneous longitudinal velocity, υ x , separated for minus-end directed, υ x > 0, and plus-end directed, υ x < 0, directions; data are provided in Table S6. (D) Mean (±SEM ) instantaneous minus-end directed longitudinal velocity υ x | υ x > 0 against mean run-time τ p .  Table S11). (B) Fraction of time the NPs spend in each of the -states, for different values of (data are shown in Table S12). (C) Mean (±SEM) longitudinal velocities of the different -states for different values of (data are shown in Table S13). Note that states of ≥ 5 are extremely rare for ≤ 13; therefore, in (B), the state = 5 is not shown and it appears only for ≥ 10 in (C).  Table S12).  Table S13). Note that states of M B ≥ 5 are extremely rare for N m ≤ 13; therefore, in (B), the state M B = 5 is not shown and it appears only for N m ≥ 10 in (C). Figure 7A shows the mean longitudinal velocity, υ x , against N m . We found a suppression of υ x from the single motor velocity (833 ± 4 nm/s) as N m increases, which is a clear signature of motor-motor coupling in the absence of any load-sharing effect (due to the vanishing drag). The corresponding standard deviation of the mean (STD) also varies with N m ; the maximal STD appears at N m = 4 (corresponding to ξ = 35.4 nm), whereas the minimal STD corresponds to the single-motor case. In addition, as seen in Figure 7B, the characteristic run-time, τ p , and longitudinal run-length, λ , are also affected by N m . We observed an effectively monotonous increase of τ p and λ with an increase of N m up to N m = 13. Following the experimental analysis shown in Figure 4C, we computed, separately, the mean longitudinal velocity of minus-end and plus-end directed steps ( Figure 7C) and found a monotonous decrease of the mean minus-end directed velocity with increasing N m , effectively saturating above N m 7, similar to the trend seen in the experimental results; the velocity of the plus-end directed steps was not sensitive to N m . Noting that Figure 7B,C can be regarded as a "parametric presentation" of υ x | υ x > 0 vs. τ p , with N m serving as the parameter, we present in Figure 7D the direct dependence of υ x | υ x > 0 against τ p . The similarity to Figure 4E is striking, giving extra support to our interpretation of the experimental results.  Table S7; see Table S16 and Figure S14 for a Δ = MC step time). (B,C) Mean (±SEM) angular velocity calculated separately for minus-end directed motion, 〈 〉 | > 0, (orange bars) and plus-end directed motion, 〈 〉 | < 0, (green bars), against the number of NP-bound motors, , for a time interval Δ that equals to the MC-step time (B; data are provided in Table S18 and 0.27 s (C; data are provided in Table S8).  Table S7; see Table S16 and Figure S14 for a ∆t = MC step time). (B,C) Mean (±SEM) angular velocity calculated separately for minus-end directed motion, ω | υ x > 0, (orange bars) and plus-end directed motion, ω | υ x < 0, (green bars), against the number of NP-bound motors, N m , for a time interval ∆t that equals to the MC-step time (B); data are provided in Table S18 and 0.27 s (C); data are provided in Table S8.  Table S19). (B) Mean and STD of the helical pitch size, regardless of the helical vorticity (data are provided in Table S20). (C) Comparison between the actual helical pitch size (from the observed simulated trajectories) and the helical pitch size estimated using the mean values of angular and longitudinal velocities, 〈 〉 = |〈 〉| |〈 〉| (data are provided in Table S21).

3.4.1.
Longitudinal velocity: Figure 7A shows the mean longitudinal velocity, 〈 〉, against . We found a suppression of 〈 〉 from the single motor velocity (833 ± 4 nm/s) as increases, which is a clear signature of motor-motor coupling in the absence of any load-sharing effect (due to the vanishing drag). The corresponding standard deviation of the mean (STD) also varies with ; the maximal STD appears at = 4 (corresponding to = 35.4 nm), whereas the minimal STD corresponds to the single-motor case. In addition, as seen in Figure 7B, the characteristic run-time, 〈 〉, and longitudinal run-length, 〈 〉, are also affected by . We observed an effectively monotonous increase of 〈 〉 and 〈 〉 with an increase of up to = 13. Following the experimental analysis shown in Figure 4C, we computed, separately, the mean longitudinal velocity of minus-end and plus-end directed steps ( Figure  7C) and found a monotonous decrease of the mean minus-end directed velocity with increasing , effectively saturating above ≃ 7, similar to the trend seen in the experimental results; the velocity of the plus-end directed steps was not sensitive to . Noting that Figure 7B,C can be regarded as a "parametric presentation" of 〈 〉 | > 0 vs 〈 〉, with serving as the parameter, we present in Figure 7D the direct dependence of 〈 〉 | > 0 against 〈 〉. The similarity to Figure 4E is striking, giving extra support to our interpretation of the experimental results.  Table S19). (B) Mean and STD of the helical pitch size, regardless of the helical vorticity (data are provided in Table S20). (C) Comparison between the actual helical pitch size (from the observed simulated trajectories) and the helical pitch size estimated using the mean values of angular and longitudinal velocities, H = 2π | ω | | υ x | (data are provided in Table S21).
To elucidate the connection between these results and the actual number of motors participating in the motion, in Figure 8A, we depict the mean number (per MC unit time)-over all runs-of transient MT-bound motors, M B , for each N m (see Table S11). As expected, as N m increases, the value of M B increases; yet, surprisingly, it effectively saturates at around M B = 2. We next extracted the fractions of time for different NP "states"; we define an NP state by the number of MT-bound (transporting) motors M B ( Figure 8B and Table S12). Note that although M B ranges roughly between 1 and 2, the distribution of the motor numbers is wide and shows that there is a significant contribution of the one-, two-, and three-motor states, whereas the contribution of M B -states of four motors and above is negligible. Figure 8C shows the corresponding mean velocity for each state (see also Table S13), avoiding NP temporal velocities associated with transitions between these states via binding-unbinding events. Surprisingly, the mean velocity of a particular M B -state increases with the increase in N m , rather than remaining constant. This can be rationalized by noting that an M B -state still corresponds to several microscopic configurational states. At larger N m , i.e., smaller NP anchoring distances, most of these microscopic states are associated with MT-bound motors whose NP anchoring distance ξ is smaller. This implies a reduction in polymer tension, i.e., suppression of the elastic coupling between MT-bound motors. Indeed, such an effect has been observed in previous theoretical studies [32,33,36,42] and is associated with the fact that a forward-pulling force on a lagging motor, while somewhat increasing its temporal velocity, has a much smaller effect than the reduction in temporal velocity due to a backward-pulling force on a leading motor. The present analysis confirms that the reduction in the longitudinal velocity with increasing N m , which we observed both experimentally and in the simulations, resulted from the increased number of the MT-bound motors participating in the motion, M B .
As manifested by the STD shown in Figures 7A and S11 and by the corresponding histograms in Figures S10 and S12, respectively, the width of the NP longitudinal velocity distribution varies with N m . We associate this distribution width both with the width of the velocity distribution at each M B -state (not shown here) and with the probability distribution (i.e., time fractions) of these states ( Figure 8B), which is sensitive to the value of N m . For N m ≥ 2, we unexpectedly observed a second peak around υ x = 0, which becomes more pronounced as N m increases ( Figure S10). This phenomenon is more strongly demonstrated by the equivalent histograms corresponding to a single MC time-step ( Figure S12). We associate this new peak to the growing abundance of states M B ≥ 2 with increasing N m ( Figure 8C). These states are prone to "jamming configurations", i.e., configurations wherein the "jamming events" (controlled by motor-motor excluded-volume interactions, see Supplementary Material Section S6) dominate.
From this analysis, we conclude that the increase in N m leads to a non-linear increase of M B , which, in turn, leads to longer run-times, τ p , and diminishing υ x . Since the effect on the run-time is more pronounced than on υ x , the resulting run-length, λ , is, in most cases, also enhanced. Moreover, the increase in M B leads to a slightly narrower NP longitudinal velocity distribution for N m > 4. While an increase in run-time is expected even without the inclusion of motor-motor coupling [16], the decrease in velocity is a sole consequence of the (elastic and excluded volume) motor-motor interactions. Similarly, the non-trivial change of width of the υ x distribution with increasing N m (STD in Figure 7A) reflects the competition between (i) the increase in M B fluctuations (i.e., the corresponding STD, δM B ; Table S11), which acts to increase the velocity fluctuations, (ii) the variability of the velocity fluctuations within each M B -state (i.e., δv x,M B ; Table S14), which can be attributed to motor-motor coupling, and (iii) the contributions of the unbinding events, since after each unbinding event an immediate jump of the NP position occurs to balance elastic forces; unbinding events increase with increasing N m ( Figure S13 and Table S15).

Angular/Transverse Motion
Consider, next, the angular motion around the MT symmetry axis, which leads to a transverse motion on the projected base plane. As can be seen in Figure 6 and Movies 4-7, the motion is apparently composed of both left-and right-handed helices, combined with large fluctuations. However, on a very long trajectory, the net helical motion (in which the right-and left-handed helices cancel each other) might appear minor and will not reflect the true nature of the motion. Therefore, here, we require very delicate analyses that will reflect both characteristics.
The variation of the mean angular velocity ω and its STD with N m is shown in Figure 9A and in Figure S14, both demonstrating a monotonous decrease with an increase in N m . The fact that all values are positive suggests a net right-handed helical motion of the NP for all motor numbers. Combined with a variation of the mean longitudinal velocity with N m (Figure 7A), and the relationship H = 2π| v x |/| ω | describing the mean helical pitch size H , this leads to the dependence of H on N m , as shown in Figure 10C, exhibiting an increase in H for growing N m .
As discussed above, this pitch size represents the net helical motion, i.e., it includes the cancellations of left-and right-handed helices (as seen in Figure 6 and Movies 4-7). To refine the helical motion analyses, we define a helix (be it left-or right-handed) whenever the angular motion completes a full round (i.e., φ = ±2π). Note that, even within a single full round, the motion consists of large and frequent right-and left-handed fluctuations.
In Figure 10A, we show the resulting mean pitch size for left-(φ = −2π) and right-(φ = +2π) handed helices separately, and the mean pitch size regardless of the helix vorticity. Notably, the pitch sizes obtained for these three definitions are all comparable to each other and are much shorter than those deduced using the mean angular velocity ( Figure 10C)-all increase with increasing N m . To complement these results, we also present, in Figure S15 and Table S17, the mean angular velocity separately for left-and right-handed motion.

Longitudinal and Angular/Transverse Motion Are Correlated
To gain further insight into the complex motion of the NP, we determined whether the longitudinal and angular motions are correlated. In Figure 9B, we dissect the angular velocity for forward and backward steps associated with MC time-step (i.e., positive and negative longitudinal velocities, respectively). Both angular velocity types, associated with either forward or backward steps, show a decrease with increasing motor number N m . In addition, the results show that the mean angular velocity is greater (more positive) for backward steps, with the single motor NP showing the highest value. This finding implies that when a backward step is performed, the NP is likely to move to the right. For comparison, we also plot, in Figure 9C, the same distributions on the experimental time interval (0.27 s), showing that these correlations are suppressed, highlighting the importance of the time interval for making the basis for comparison between different results, be it experimental or theoretical.

Binomial Distribution of the Number of NP-Bound Motors, N m
For a site-independent binding process, we may assume that the binomial distribution of N m (between the different NPs) holds. Using this distribution, in Supplementary Material Section S9, we computed some of the reported motility characteristics for a varying N m . For the mean longitudinal and transverse velocities calculated on the MC time scale, this analysis shows a small increase in their values (relative to the values for a deterministic motor number, N m ), especially for small motor numbers (Figures S16 and S17); this phenomenon is due to the contribution from NPs with very few motors (N m = 1, 2, 3), whose velocities are always higher ( Figure 7A) and, thereby, contribute more to the mean. However, the trends that we deduced for deterministic number of motors (Figures 7 and 9) are not altered. Importantly, when moving to the experimental time scale (0.27 s), the difference is almost non-discernable ( Figures S16 and S17). Moreover, the run-times and run-lengths (for which the selected time interval is irrelevant) are even less affected by the N m fluctuations, probably due to the contributions from NPs with a high N m , for which these values are much higher.

Discussion
A comparison between experiments and simulations requires the consideration of two issues: (i) knowledge of the relative NP motor coverage θ = ξ * ξ 2 , which corresponds to the mean number of NP-bound motors via N m = N × θ, where N is the mean number of PEG-NLS polymers per NP. We can reasonably assume that regardless of the [NLS], an identical CE batch and [CE] imply identical θ, which allows a qualitative comparison between different systems, as discussed below; (ii) the expected (wide) binomial distribution of N m between different NPs. We have concluded from simulations that, on the experimental time scale (0.27 s), the effect of the binomial distribution is relatively minor. Hence, we no longer emphasize the difference between (simulation) results for deterministic N m and those corresponding to N m using the binomial distribution.
The comparison between experiment and theory shows several similarities regarding the dependence of the longitudinal motion on the number of NP-bound motors. Note that, in the experiment, we estimated the number of PEG-NLS, N , for the different [NLS] used for the motility assays (Tables 1 and S1). An increase in N , when we use the same CE batch and [CE] (i.e., fixed θ), implies a proportional increase in N m . Both the model simulations and the experiment show a decrease in the mean longitudinal velocity with increasing N m ( Figures 4A and 7A and Tables S1 and S5). In the model simulations, the mean longitudinal velocity v x gradually levels off when N m is increased, rather than continuing to diminish strongly ( Figure 7A). This phenomenon occurs mainly due to the "self-regulation" of the mean number of transporting motors, M B , which effectively saturates at the value of two ( Figure 8A). If M B would have continued to increase with increasing N m , v x should have continued to decrease strongly due to motor-motor coupling. A very similar effect can be inferred from the experimental results ( Figure 4A and Table S1). A comparison of v x between system I (815 nm s , Table S1) and the model simulation results (833 nm/s in Table  S5) suggests that the mean number of NP-bound motors in system I should be N m 1. Thus, we can assume (consistent with the estimated values of N 5, 10, 10, 35) that N m 1, 2, 2, 7 for systems I, II, III, IV, respectively. First, we consider systems I and II. The drop in the simulation value of v x from N m 1 to N m 2 is about 7%, whereas the experiment yields a 6% drop. Similarly, we observe a drop of about 26% in the simulation value of v x , from N m 2 to N m 7, whereas the experiment yields a 32% drop, which can be regarded as a satisfactory agreement. This semi-quantitative agreement provides indirect evidence that M B is "self-regulated" in the experiment, as it is in the simulations.
As discussed in the Results section, in Figure 4E, we presented the experimentally observed dependence of υ x | υ x > 0 against τ p . We suggested that this behavior could emerge from a parametric dependence on the experimentally hidden variable N m . Since N m is a controlled parameter in the simulations, in Figure 7D, we combined the data presented in Figure 7B,C (for υ x | υ x > 0 and τ p against N m ) to obtain the prediction of the simulations for the direct dependence. The similarity between Figures 4E and 7D is evident, providing further support to our interpretation. In addition, the broad scatter of the single NP data, shown in Figure 4E, can be explained by the expected broad (binomial) distribution of N m discussed above; this is due to the small number of N m associated with the nanoscale size of the NPs.
One of the most peculiar features of the angular motion, seen only on the MC timescale of simulations ( Figure 9B and Table S18), was its correlation with the direction of the longitudinal motion (i.e., minus-end or plus-end motion along the MT long-axis). We found that plus-end directed steps dictate (mostly) large right-handed steps. However, for much longer time intervals in the simulation (e.g., the experimental timescale 0.27 s), these correlations disappear ( Figure 9C) and, consistently, they do not appear in the experimental results. Moreover, this characteristic also appears as a single-motor property, and it diminishes as more motors become available (increasing N m ). In fact, our single mammalian dynein stepping model, which builds on the raw data of Yildiz and co-workers [38], shows that these correlations appear strongly in the single-motor stepping. (Notably, revisiting the results of our single yeast dynein stepping model, we found similar correlations [40]). Thus, further investigation of this peculiar motility feature is required to determine whether it is evident in multi-motor complexes.
Our model simulations did not include two main experimental setup constraints and boundary conditions. First, they permitted the free motion of the NP around the MT surface (similar to the bridge-like setup of Yildiz and co-workers [37]), unlike our experimental setup in which the motion might be impaired due to the presence of the impenetrable glass surface. Second, they do not include MT ends (i.e., MTs are considered infinite) and MT junctions (i.e., only isolated MT tracks are considered), which will be studied elsewhere.
It is important to compare our experimental-theoretical studies to previous results in the literature. Previous gliding assays using kinesin motors [83] found that, as the density of motors increases, the gliding velocity of the MTs decreases, which is in qualitative agreement with our findings. This phenomenon is expected when the load-sharing effect, trending for the opposite, is negligible. However, quantitative comparisons with gilding assays are, in general, not appropriate since the moving MT in a gliding assay is effectively a 1D "cargo" (as studied theoretically in [32]), whose geometry facilitates the fast binding of dynein to the MT if the unbinding of an internal motor occurs; binding and unbinding processes occur mainly at the MT ends-binding at the plus-end and unbinding at the minus-end-as the MT moves. To the best of our knowledge, there is no other systematic experimental study of the effects of motor crowding on spherical cargo velocity; a velocity reduction due to motor crowding was indeed found in previous, albeit more simplified, theoretical studies [16,49,50].

Conclusions
In this paper, we report a combined experimental and theoretical study of NPs that are carried by multiple mammalian dynein motor-proteins. We studied different motility characteristics of these multi-motor NPs and compared them to the known and wellstudied case of the single dynein [6,[17][18][19][20]38,39,85,[89][90][91]. We focused on a single NP size (R = 20 nm), although the complexity of motor-motor interactions suggests that motility properties would be strongly affected by varying the NP size; this investigation will take place in future studies.
Although the number of transporting (MT-bound) motors, M B , was found to increase with the increasing number of NP-bound motors, N m , it remains at a mean value of ≤ 2 with dominating states M B = 1, 2, 3. Thereby, the NP can achieve the following motility features: (i) longer run-times and run-lengths than those of a single motor (simulations and experiment); (ii) substantial angular velocity (simulations and experiment), resulting in pronounced helical motion (simulations), yet with a pitch size longer than the single-motor one; (iii) significant longitudinal velocity (simulations and experiment), which is somewhat reduced compared with the single-motor velocity; (iv) plus-end directed motional intervals that are strongly correlated with large right-handed motion (simulations), especially during the temporal state of a single transporting motor; and (v) a maximum width of the longitudinal velocity distribution at an intermediate number of NP-bound motors (simulations and experiment).
Some native cargos are rigid, such as the (cores of) HIV [13], herpes-simplex virus [92], and adenoviruses [14,15], and super-coiled plasmids [93]; others are relatively flexible, as in the case of cytoplasmic lipid granules (i.e., native liposomes). Our current model simulations, as well as previous theoretical results [32,33,42] and the present motility assays of a rigid NP, demonstrate that, for rigid cargos, a maximal transport efficiency can be achieved when flexible linkers mediate between the rigid NP body and the motor proteins, allowing the NP to optimally use the viable motors. We, thus, conjecture that rigid native cargos use a similar mechanism for their active intracellular transport. For instance, it is possible that the abundance of disordered loops in the NFkB structure gives it extra flexibility for the enhancement of the super-coiled plasmid motility [93].
If used intracellularly, the ability of the NP to make long trajectories enhances its active transport towards the nucleus, which is relevant both to the rigid core viruses mentioned above and, possibly, to drug delivery applications using nano cargos [16,94,95]. In addition, by maintaining its helical motion-and, in particular, by making large right-handed steps when stepping backward-we posit that the NP is capable of bypassing obstacles on the MT tracks or in their vicinity. Indeed, this suggestion (albeit without mentioning the backward-right-handed correlated motion) has already been put forward [37] and was further observed in (surface-free) bead motility assays [96]. We plan to further investigate these hypotheses in future work.

Data Availability Statement:
The data presented in this study are available on request from the corresponding authors.