A cycling state that can lead to glassy dynamics in intracellular transport

A cycling state that can lead to glassy dynamics in intracellular transport

Monika Scholz James Franck Institute, the University of Chicago, Chicago, IL 60637 Institute for Biophysical Dynamics, the University of Chicago, Chicago, IL 60637    Stanislav Burov Department of Physics, Bar-Ilan University, Ramat-Gan, 5290002 Israel    Kimberly L. Weirich James Franck Institute, the University of Chicago, Chicago, IL 60637 Institute for Biophysical Dynamics, the University of Chicago, Chicago, IL 60637    Björn J. Scholz Enrico Fermi Institute, the University of Chicago, Chicago, IL 60637 Kavli Institute for Cosmological Physics, the University of Chicago, Chicago, IL 60637 Department of Physics, the University of Chicago, Chicago, IL 60637    S. M. Ali Tabei Physics Department, University of Northern Iowa Cedar Falls, Iowa 50614    Margaret L. Gardel James Franck Institute, the University of Chicago, Chicago, IL 60637 Institute for Biophysical Dynamics, the University of Chicago, Chicago, IL 60637 Department of Physics, the University of Chicago, Chicago, IL 60637    Aaron R. Dinner James Franck Institute, the University of Chicago, Chicago, IL 60637 Institute for Biophysical Dynamics, the University of Chicago, Chicago, IL 60637 Department of Chemistry, the University of Chicago, Chicago, IL 60637

Power-law dwell times have been observed for molecular motors in living cells, but the origins of these trapped states are not known. We introduce a minimal model of motors moving on a two-dimensional network of filaments, and simulations of its dynamics exhibit statistics comparable to those observed experimentally. Analysis of the model trajectories, as well as experimental particle tracking data, reveals a state in which motors cycle unproductively at junctions of three or more filaments. We formulate a master equation for these junction dynamics and show that the time required to escape from this vortex-like state can account for the power-law dwell times. We identify trends in the dynamics with the motor valency for further experimental validation. We demonstrate that these trends exist in individual trajectories of myosin II on an actin network. We discuss how cells could regulate intracellular transport and, in turn, biological function, by controlling their cytoskeletal network structures locally.

thanks: To whom correspondence may be addressed: dinner@uchicago.eduthanks: To whom correspondence may be addressed: dinner@uchicago.edu

I Introduction

Individual microscopic particles (beads wong2004anomalous (); wang2009anomalous () or fluorescently labeled molecules parry2014bacterial (); sako2000single (); cai2009single ()) can now be tracked in cells. These studies reveal complex dynamics zimmerman1991estimation (); konopka2006crowding (); hofling2013anomalous (). The resulting trajectories can be treated as random walks, and quantitative analysis of their statistics can provide insights into underlying mechanisms tabei2013intracellular (); Burov et al. (2013). Often, the mean square displacement (MSD) exhibits a power-law (typically sublinear) dependence on the separation in time between two observations, known as the lag time () CoxPRL (); tolic2004anomalous (); weber2010bacterial (); weigel2011ergodic (); tabei2013intracellular (); Burov et al. (2013); bressloff2013stochastic (). In certain cases tabei2013intracellular (); weigel2011ergodic (), the MSD also decays as the amount of data included in averages (the measurement time, ) increases; this trend indicates a power-law distribution of dwell times and is known as “aging” in theories of glassy dynamic cugliandolo1994evidence ().

These correlations can have important biological implications tabei2013intracellular (); manzo2015weak (). For example, a recent study shows that the anomalous dynamics observed for insulin secretory vesicles (granules) can account for the biphasic kinetics of insulin release tabei2013intracellular () without the need to invoke separate pools of granules, as previously Seino2011 (). In particular, the sustained release relies on the glassy dynamics. Glassy dynamics are often interpreted in terms of trapping in local minima of an energy landscape with an exponential or power-law distribution of depths bouchaud1995aging (); feigel1988stochastic (). However, how such a landscape could arise from typical biomolecular interactions is unclear. Crowding is insufficient, as it results in standard Brownian motion but with a reduced diffusion coefficient Dix2008 (). Because the moving vesicles are associated with molecular motors, which consume cellular energy stores (nucleotide triphosphates) for directed motion along cytoskeletal filaments, other, intrinsically nonequilibrium mechanisms of generating these statistics may exist.

Figure 1: Schematic of a cytoskeletal assembly comprised of multiple molecular motors. Drawing is to scale for actin filaments (red) and myosin motors (blue).

In the case of insulin release, the vesicles have both kinesin and dynein associated with them varadi2003kinesin (), which walk in opposite directions on microtubules vale1992directional (). More generally, many cytoskeletal assemblies in cells have multiple motors associated with them (Fig. 1). Thus it is natural to ask if glassy dynamics could arise from a tug-of-war. Mathematical models of competing motors have been formulated and show that, depending on the number of motors, their properties, and their binding affinities, different regimes of transport kinetics can be accessed muller2008tug (); newby2010local (). While the tug-of-war model can explain bidirectional transport in a diverse set of biological systems hendricks2010motor (); soppina2009tug (); ali2011myosin (); hancock2014bidirectional (), standard formulations cannot account for aging because exponential dwell times are assumed. Moreover, the tug-of-war is essentially a single-filament mechanism. Cytoskeletal networks contain geometric structures that involve multiple filaments (e.g., junctions), and these could support other dynamics.

In this paper, we investigate the dynamics of a minimal model of a motor that can make multiple attachments to a two-dimensional network of filaments. Using simulations, we show that such a model can exhibit glassy dynamics, and we discover that the long-time correlations in this model result from vortex-like trajectories that motors follow when three or more filaments cross to form a circuit. This represents a new mechanism for trapping that does not require individual motor heads stalling, and we term it the “cycling state”. We obtain average flows for idealized junction geometries from a master equation analysis and show that trapping in vortex-like cycling states can give rise to glassy, non-ergodic dynamics like those observed in experiments. We analyze experimental particle tracking data to demonstrate the presence of the cycling state in measured trajectories and show that it relates to exponents that quantify aging. Broader implications for biological function are discussed.

Figure 2: Schematic of the molecular simulation procedure. (A) The motor (blue dot) can associate with filaments (red vectors) within a defined binding radius. (B) Active attachments are assigned stochastically in a distance-dependent fashion. The velocity vector that results from each active attachment moving a step along its filament is computed. (C) The velocity vector is projected onto the filaments. (D) The motor takes a step with size proportional to the magnitude of the largest velocity projection, along the associated filament.

Ii Molecular simulations

We model filaments as randomly oriented line segments in a plane. The length of each filament is drawn from a normal distribution with a maximum length and standard deviation . We associate the polarity of filament with a fixed unit vector . The filaments are static, and thus represent experimental situations in which cytoskeletal rearrangements are slow in comparison with the period over which motor transport is measured. For simplicity, we also neglect heterogeneities in the composition of the filaments and the solution environment, which can lead to complex dynamics newby2010local ().

We are interested in cases in which many molecular motors act in concert—e.g., a vesicle with several protein motors attached or a myosin minifilament. We refer to our model of such an assembly as “a motor”. A motor is a point particle that can bind up to filaments at once and move along them as follows. We separate the binding process into two steps (Fig. 2, panels 1 and 2; see SI Text and Fig. S1 for an alternative scheme). First, we distribute the possible attachments for a motor among filaments with probability proportional to , where is the shortest distance between the motor and filament , and is a parameter that sets the interaction length scale. This probability is normalized by the sum over all filaments within a distance of the motor. Then, we determine if each such interaction exerts force to move the motor (henceforth, “active”) with probability or not (“inactive”) with probability . We denote the number of active attachments to filament by . To determine the change in position of the motor we compute the vector (orange in Fig. 2B). This choice is consistent with measurements that show that motor velocities increase with head numbers lee2015axonal (). We project onto all the filaments with at least one active attachment and add to the motor position the projection with the maximum magnitude scaled by the time step (Fig. 2C and D). This projection rule ensures that the motor moves along rather than off filaments. It also gives rise to an effective force-velocity dependence, with opposing parallel velocities canceling each other. We assume that forces that are directed orthogonally do not create a load on the motor. For the simple scenario of orthogonal filaments, this scheme thus simplifies to a step in the direction of the filament with the highest number of active binding interactions (i.e., a majority rule). The projection rule implies that binding sites do not detach under load; rather, they stay bound and contribute to the overall velocity vector. Simulations relaxing the projection rule to a simple net velocity calculation show only minor differences (Fig. S2). Simulations with stochastic selection between projections with probabilities proportional to their magnitudes also yielded similar results.

Parameter Value
Filament density per unit area 1
Filament length 5
Filament standard deviation 5
Binding radius 0.01
Total number of binding sites 50-100
Time step 0.001
Total number of steps
Trials 2000
Table 1: Simulation parameters

The values of the simulation parameters are given in Table 1. While the model is general, we choose the values to be roughly consistent with actin and myosin to ensure that we study a physically reasonable regime. To this end, we assume a myosin speed of 1 m/s, which is in the range of speeds reported from in vitro and in vivo studies of various classes of myosin finer1994single (); howard1997molecular (); Pierobon20094268 (). We associate our unit length with 1 m, such that the average filament length is 5 m. The binding radius of the motor is then about 10 nm, and a single time step of the simulation is 0.1 s. Although the properties of actual molecular motors vary substantially, our conclusions are robust to parameter choices that range over an order of magnitude (Fig. 3).

Figure 3: Scaling in molecular simulations. (A) Time-averaged MSD as a function of lag time for binding radii (gray), 0.01 (blue), 0.1 (orange) and 1(red). (B) Mean-squared displacement as a function of measurement time for a range of binding radii (colors are the same as in (A)). Each curve was rescaled to start at 1 for easier visualization. (C) and (D) Dependence of the indicated exponents on the maximum number of possible attachments, . The dashed lines show the outcome for a tug-of-war scenario, which can be obtained from our model by simulating the dynamics on two antiparallel filaments. The full lines show the exponents for the molecular simulations on a two-dimensional random network. Lag-time exponents were calculated for .

We simulate the model according to the rules above and calculate the time-averaged mean square displacement (MSD) for the resulting trajectories:


where is the position of a motor at time . We plot the time-averaged MSD as a function of lag time () in Fig. 3A. We use the time-averaged MSD to make connection with biological experiments that typically have a limited number of trajectories saxton1997single (). It is important to note that the time- and ensemble-averaged MSDs can exhibit different scalings when the process of interest is non-ergodic bouchaud1992weak (); lubelski2008nonergodicity (); he2008random (); sokolov2008viewpoint (); schulz2013aging ().

By varying the binding radius by factors of 10 from 0.001 to 1, we can tune the transport from superdiffusive to subdiffusive. When the range of interaction is very small, the motor only attaches to the closest filament. As a result, the motion is ballistic but slow (note the intercept for the gray line in Fig. 3A) because the number of active attachments is low. For intermediate interaction ranges, the motor can simultaneously bind multiple filaments, and a tug-of-war-like mechanism gives rise to subdiffusive or diffusive motion.

We plot the MSD as a function of the measurement time () in Fig. 3B. For an ergodic system, the MSD should be constant as varies—the properties of the motion are the same independent of the length of recording. In contrast, a decrease in the MSD with suggests that there are long-lived traps. The essential idea is that the traps increasingly dominate the statistics as more data are included in the averages. We observe a power-law decay (i.e., aging), consistent with prior analysis of myosin II on an in vitro actin network Burov et al. (2013) and the motion of insulin granules in vivo tabei2013intracellular (); Burov et al. (2013).

The observed aging exponent depends on the size of the binding radius , but we obtain exponents comparable to experimental values as this parameter varies over two orders of magnitude. The molecular simulations of the model also predict more trapping for increased numbers of binding sites: both the lag-time and the aging exponents decrease with larger (Fig. 3C and D). However, due to the non-ergodic nature of this transport process, the lag-time exponent depends on the length of the recording . We show the -dependence of the exponents in Fig. S3 (see also Fig. S4). In contrast to the results that we obtain for motion on a random filament network, a simple tug-of-war scenario with a motor between two antiparallel filaments yields no glassy dynamics, independent of the number of binding sites. Also, the scaling of the MSD becomes increasingly ballistic with increased numbers of binding sites (Fig. 3C and D), and the diffusion constant increases (data not shown).

Figure 4: Single-particle trajectories reveal a cycling state. (A) Representative simulation trajectories (blue and red) projected onto the filament network (gray). (B) A magnified view of a cycling trajectory (red).

Having thus captured the statistics of experiments, we sought to use the model to elucidate the microscopic motions that underlie the statistics. In this regard, we noticed that motors frequently exhibit vortex-like motions in which they steadily cycle from one filament to another at a junction (Fig. 4A). These motions persist for long times in comparison with the duration of the simulations. Cycling motions are of particular interest given that passive particles in a vortex flow in a fluid are known to exhibit power-law-distributed trapping times solomon1993observation (). Analogous observations exist for trapped ions devoe2009power () and Bose-Einstein-Condensates lundh2000hydrodynamic ().

An approximate length scale for the vortices responsible for the glassy dynamics can be inferred from the crossover in the MSD curves. When the MSD as a function of lag time switches from super-linear to linear/sublinear scaling, is sufficiently large for the MSD to include contributions from the vortices (Figs. S1 and S2). In other words, the MSD exponent decreases when the dynamics include bounded trajectories. In our conditions, the crossover occurs at , which corresponds to about 10 s based on the numbers for actin and myosin given above. In turn, for a myosin II motor and the parameters in Table 1, the length scale of vortices would be approximately 100 nm.

Iii Cycling state gives rise to power-law dwell times

To investigate whether the observed vortex-like motion can account for the anomalous statistics, we now consider an idealized geometry. Specifically, we consider filaments that meet at right angles to form a square because it simplifies the mathematics (Fig. 5). However, we emphasize that the conclusions that we draw from this analysis are general—cycling can occur whenever the unit vectors of a group of crossing filaments sum to zero. In a square loop, the motor only interacts with one or two filaments at a time, and, due to the projection rule (Fig. 2), the motor always moves along the filament with the majority of sites bound. The average resulting motion can be described by streamlines, which form a nested set of closed loops that do not cross (Fig. 5A). If the dynamics were deterministic, the streamlines would describe the motion entirely, and the motor would stay in the vortex forever. However, due to the stochastic nature of binding and unbinding, individual trajectories deviate from the streamlines, and the motor eventually leaves the vortex.

To characterize this behavior quantitatively, we derive the master equation that governs this escape. For this purpose, we need to determine how the active binding sites distribute between the two accessible filaments via the two-step procedure described in Molecular Simulations. The probability that out of possible attachment sites are assigned to filament 1 and that the remaining are assigned to filament 2 is binomial:


The probability that out of possible attachments are active is also binomial:


and similarly for . Since filament assignment and selection of the active attachments are independent events, we can now write for the overall probability of having and active binding sites on filaments 1 and 2


The limits of the sum are set by , , and .

Associating filament 1 with the direction and filament 2 with the direction, the resulting master equation for motion in a square vortex is

Here, is the probability to be in a certain location in the vortex at time , and we assume a unit time step. The first term accounts for motors that stay in place, while the second and third terms account for taking steps along filaments 1 and 2, respectively. Note that the polarity of the filaments is fixed, so that motion along each filament is unidirectional and no terms are needed to account for motion in the other direction. The vortex has absorbing boundaries where the velocity has a saddlepoint, forming a diamond shaped region.

Figure 5: Vortex model. (A) Streamlines (blue) for motion in the vicinity of an idealized square circuit of filaments (gray). (B) Survival probability in a vortex for (red) or (blue) in a vortex of size 40. The probability is initialized 3/4 of the diagonal. The binding radius changes from 1/5, 1/7.5 or 1/10 of the vortex length. (dotted, dashed or full lines). The probability decays in three phases: constant, power-law, and exponential. The length of the power-law decay decreases with increasing binding radius. (C) Numerical solution of the master equation in a square vortex. In the example shown, the initial condition is a localized probability at a single location similar to B. The vortex size is 40, with =5, a binding radius of 1/10 of the side length and a unit length step size. The colors correspond to the logarithm of the probability density of observing the motor at a specific location (red is highest) and is normalized in each subplot. The three panels show the master equation after 5, 25, and 750 iterations, respectively. The probability density spreads over time and the maximum rotates. The final image shows a pattern of localization that is stable as the integrated density continues to decrease.

We solve this master equation numerically (Fig. 5C). The survival probability in the vortex decays in three stages. At short times, the probability is close to unity, since the probability distribution needs a finite number of steps to reach the vortex boundaries. Then, we observe that the probability density also rotates within the vortex, similar to particles in the explicit simulations (see also Supplementary Movie 1). During that time, the probability decays as a power law (see Fig. 5B). The probability distribution in the vortex ultimately reaches a quasi-steady-state, when the shape of the distribution is no longer changing, and the survival probability decays as an exponential. The durations of these three stages are determined by the starting conditions and the relative size of the binding radius to the size of the vortex, as well as the number of binding sites. The power-law decay of the survival probability is longest for small binding radii (significantly smaller than the vortex) since the binding potential favors steps parallel to the filaments, and only allows very small steps orthogonal to the filaments. This relates to the results shown in Fig. 3A and B, in that there exist a range of vortex sizes in our chosen network, which lead to characteristic trapping times. Binding radii much smaller or larger than these sizes will show less aging than midrange radii, due to the existing spatial scales in the random network. Relatedly, increasing the density of filaments in the network decreases the vortex size; this leaves the scaling with unchanged while decreasing the extent of aging (Fig. S5).

The average motion in the vortex can be viewed as a combination of drift along the streamlines and diffusion orthogonal to the streamlines. From this perspective, the cycling motion is similar to that of a particle in a Rayleigh-Benard convection cell shraiman1987diffusive (). However, the situation differs, in that the diffusion is position-dependent (i.e., varies in space through and ). In other words, the master equation limits to a drift-diffusion (Langevin) form with a multiplicative, rather than an additive, noise. This form reflects the fact that, in the motor model, the number of attachments is influenced by the position relative to the filaments defining the vortex boundary.

Iv Experimental demonstration of cycling state contributions

A key prediction of the cycling state model is that the aging exponent decreases with the number of attachment sites on each motor (Fig. 3D). We can test this idea without the confounding effects of cell signaling by studying mixtures of purified actin filaments and myosin motors. The specific system that we consider comprises actin filaments bundled by a passive crosslinker fimbrin. The motors are minifilaments of skeletal muscle myosin II, which polymerizes into large assemblies with on the order of 100 motor proteins thoresen2013thick (). The actin and myosin molecules are visualized through fluorescence microscopy, as detailed in the SI Text. Single-particle trajectories were obtained by tracking, and trajectories with fewer than 30 timepoints were discarded. There are 246 such trajectories with a mean length of 166 s, with a standard deviation of 120 s. This length is long compared with typical in vitro particle-tracking studies of isolated motors on single filaments. Both the large number of heads and the density of binding sites in the filament network make it unlikely for motors to detach, which favors processivity. The mean-squared displacement as a function of measurement time shows aging, with an exponent comparable to previously published data (compare Burov et al. (2013) with Figs. S6A and B).

To test the prediction in Fig. 3D, we exploit the fact that the number of myosin molecules in each minifilament varies naturally and use fluorescence intensity as a proxy for the number of motor heads. We divide the trajectories into two groups according to the median fluorescence of the motor. There are 132 single-particle trajectories for motors with relatively high intensity and an equal number of trajectories for motors with relatively low intensity. We expect motors with higher intensity to have more heads and thus exhibit stronger aging. We observe that this is the case (Fig. 6A). Specifically, we used case resampling (a form of bootstrapping) to get a distribution of exponents from each fluorescence group. The means of these distributions were and (mean SEM high fluorescence and low fluorescence groups, respectively). We tested if the means are significantly different using the unequal variances -test. The two-tailed -value is less than 10, indicating that the difference is significant.

Additionally, we manually selected trajectories that visibly cycled (looped over the same multipixel region more than once). 29 trajectories were found to cycle, and a significant portion of these were also classified as high-fluorescence (, -value = 0.00185, Fisher’s exact test for cycling/non-cycling vs. high/low fluorescence). Because the actin filaments do not move significantly (Fig. S6), we can project trajectories onto the network, and we see that the trajectories that result in the most pronounced decay have visible cycling that coincides with filament junctions (arrows in Figs. 6B-D).

While the observations above demonstrate the existence of cycling, we sought to exclude alternative mechanisms for the statistical differences in exponents in Fig. 6A. Reasonable candidates are intrinsic differences in motor speeds and detachment. With regard to the former, we note that the average motor step size from frame to frame (i.e., the net speed for movement over the filament network) does not differ for the high and low intensity groups: the means SEM are nm/s and nm/s, respectively. If detachment were instead responsible for the anomalous dynamics, the particles should undergo simple diffusion when apparently trapped. To test for this possibility, we divided the trajectories between trapped and non-trapped periods and then analyzed their dynamics using a published measure for detecting complex dynamics at resolutions of only a few frames. We find that the motion is not simple diffusion, as detailed in the SI Text (see also Fig. S9).

Together, these results strongly support the idea that the cycling state exists for motors with multiple attachments and gives rise to power-law decay in the MSD as a function of measurement time.

Figure 6: Experimental trajectories show cycling. (A) Trajectories were divided into two groups according to the total fluorescence intensity of the motor. The group with higher fluorescence intensity (green) shows a larger aging exponent than the low-fluorescence group (yellow). (B-E) Representative trajectories of individual myosin II minifilaments overlayed onto the actin network. Arrows denote cycling events. Trajectories were obtained from myosin II minifilaments on an actin network bundled with fimbrin. Scale bar is 1.1 m; trajectories shown are imaged at 1.5 s intervals; see SI Text for details. Single-particle trajectories were obtained using the Python-based implementation of the Crocker-Grier algorithm Trackpy daniel_b_allan_2014_9971 ().

V Discussion

Motor-driven processes in cells often exhibit anomalous statistics weiss2004anomalous (); amblard1996subdiffusion (); tolic2004anomalous (); jeon2011vivo (), and these statistics can have important functional consequences tabei2013intracellular (). While existing random-walk models of aging typically start by assuming trapped states with power-law dwell times, we introduce a microscopic mechanism for how simple biologically consistent interactions can give rise to such long-lived traps. The glassy dynamics result from a vortex-like state that emerges in multiple dimensions—motors that can attach to multiple filaments simultaneously cycle unproductively at filament junctions, and the escape times from these flows are power-law-distributed over a time range set by the motor step size and the filament spacing. We demonstrate that such cycling events occur frequently in the motion of skeletal myosin II assemblies on a dense, random network of bundles of actin filaments in vitro, and we use this system to validate the predictions of the model. To the best of our knowledge, these cycling dynamics have not been appreciated previously, and we expect that their topological features will lead to rich physics beyond idealized trap models.

Microtubule structures with many intersections are observed above the basal membrane of epithelial cells, where they function in endocytic vesicle transport reilein2005self (). A study that reconstructed these epithelial microtubule networks in vitro observed a kinesin-coated bead cycling through a vortex structure ross2008kinesin (). The same study (and others osunbayo2015cargo ()) also saw a slow-down or pausing states at intersections, which is also present in our model due to the motors interacting with both filaments at the same time. While these structures are macroscopic compared to the vortices responsible for the observed glassy behavior, they show that motor associated cargo or multi-binding complexes can navigate intersections and cycles for certain geometries.

Tug-of-war models muller2008tug () contain similar molecular elements and exhibit pauses when motors are stalling or opposing forces match exactly. However, such models are essentially one-dimensional and do not give rise to power-law-distributed dwell times. Based on Figs. 3C and D, we expect the cycling state to be prevalent only when motor assemblies comprise many protein motors and directed motion dominates over thermal processes. Indeed, measured exponents for the MSD could be used to infer the size of such assemblies. However, care is needed because different filament structures favor different amounts of directed, tug-of-war, and vortex-like motions. Random networks in three dimensions have few circuits that lead to cycling, while cytoskeletal networks that are organized by specific filament binding proteins (e.g., Arp2/3), as well as quasi-two-dimensional networks that arise in cell cortices, have larger numbers of suitable junctions. Thus accurate estimates of the numbers of active motor heads requires constructing a calibration curve for the exponent for each network structure.

Understanding how motor assemblies behave on complex filament networks in cells is an outstanding challenge elting2012future (). The degree to which the cycling state contributes to dynamics in different contexts in vivo is an open question deserving further study. Cells could potentially spatiotemporally control intracellular transport by rearranging their cytoskeletal networks to favor or disfavor cycling. Understanding how this control mechanism manifests in different types of cells and tissue environments, as well as its interplay with other regulatory processes and trapping mechanisms newby2010local (), is a useful direction for future research. It will also be interesting to understand the interplay of motor transport, force transmission, and network rearrangement.

Vi Acknowledgments

We thank Samantha Stam and David Kovar’s laboratory for protein purification; Michael Murrell, Jennifer Ross, Toan Huynh, and Norbert Scherer for helpful conversations; and Glen Hocky and Shiladitya Banerjee for critical readings of the manuscript. Support was provided by the University of Chicago Materials Research Science and Engineering Center (NSF DMR-1420709) and the W. M. Keck Foundation.


  • (1) I. Wong, M. Gardel, D. Reichman, E. R. Weeks, M. Valentine, A. Bausch, and D. Weitz, Phys. Rev. Lett. 92, 178101 (2004)
  • (2) B. Wang, S. Anthony, S. Bae, and S. Granick, Proc. Natl. Acad. Sci. USA 106, 15160 (2009)
  • (3) B. R. Parry, I. V. Surovtsev, M. T. Cabeen, C. S. O’Hern, E. R. Dufresne, and C. Jacobs-Wagner, Cell 156, 183 (2014)
  • (4) Y. Sako, S. Minoghchi, and T. Yanagida, Nat. Cell Biol. 2, 168 (2000)
  • (5) D. Cai, D. P. McEwen, J. R. Martens, E. Meyhofer, and K. J. Verhey, PLoS Biol. 7, e1000216 (2009)
  • (6) S. B. Zimmerman and S. O. Trach, J. Mol. Biol. 222, 599 (1991)
  • (7) M. C. Konopka, I. A. Shkel, S. Cayley, M. T. Record, and J. C. Weisshaar, J. Bacteriol. 188, 6115 (2006)
  • (8) F. Höfling and T. Franosch, Rep. Prog. Phys. 76, 046602 (2013)
  • (9) S. A. Tabei, S. Burov, H. Y. Kim, A. Kuznetsov, T. Huynh, J. Jureller, L. H. Philipson, A. R. Dinner, and N. F. Scherer, Proc. Natl. Acad. Sci. USA 110, 4911 (2013)
  • (10) S. Burov, S. A. Tabei, T. Huynh, M. P. Murrell, L. H. Philipson, S. A. Rice, M. L. Gardel, N. F. Scherer, and A. R. Dinner, Proc. Natl. Acad. Sci. USA 110, 19689 (2013)
  • (11) I. Golding and E. C. Cox, Phys. Rev. Lett. 96, 098102 (2006)
  • (12) Tolić-Norrelykke, Iva Marija and Munteanu, Emilia-Laura and Thon, Genevieve and Oddershede, Lene and Berg-Sorensen, Kirstine, Phys. Rev. Lett. 93, 078102 (2004)
  • (13) S. Weber, A. Spakowitz, and J. Theriot, Phys. Rev. Lett. 104, 238102 (2010)
  • (14) A. Weigel, B. Simon, M. Tamkun, and D. Krapf, Proc. Natl. Acad. Sci. USA 108, 6438 (2011)
  • (15) P. C. Bressloff and J. M. Newby, Rev. Mod. Phys. 85, 135 (2013)
  • (16) L. Cugliandolo, J. Kurchan, and F. Ritort, Phys. Rev. B 49, 6331 (1994)
  • (17) C. Manzo, J. A. Torreno-Pina, P. Massignan, G. J. Lapeyre Jr, M. Lewenstein, and M. F. G. Parajo, PRX 5, 011021 (2015)
  • (18) S. Seino, T. Shibasaki, and K. Minami, J. Clin. Invest. 121, 2118 (2011)
  • (19) J.-P. Bouchaud and D. S. Dean, J. Phys. I 5, 265 (1995)
  • (20) M. Feigel’Man and V. Vinokur, J. Phys. Paris 49, 1731 (1988)
  • (21) J. A. Dix and A. S. Verkman, Annu. Rev. Biophys. 37, 247 (2008)
  • (22) A. Varadi, T. Tsuboi, L. I. Johnson-Cadwell, V. J. Allan, and G. A. Rutter, Biochem. Biophys. Res. Commun. 311, 272 (2003)
  • (23) R. D. Vale, F. Malik, and D. Brown, J. Cell Biol. 119, 1589 (1992)
  • (24) M. J. Müller, S. Klumpp, and R. Lipowsky, Proc. Natl. Acad. Sci. USA 105, 4609 (2008)
  • (25) J. Newby and P. C. Bressloff, Phys. Biol. 7, 036004 (2010)
  • (26) A. G. Hendricks, E. Perlson, J. L. Ross, H. W. Schroeder III, M. Tokito, and E. L. Holzbaur, Curr. Biol. 20, 697 (2010)
  • (27) V. Soppina, A. K. Rai, A. J. Ramaiya, P. Barak, and R. Mallik, Proc. Natl. Acad. Sci. USA 106, 19381 (2009)
  • (28) M. Y. Ali, G. G. Kennedy, D. Safer, K. M. Trybus, H. L. Sweeney, and D. M. Warshaw, Proc. Natl. Acad. Sci. USA 108, E535 (2011)
  • (29) W. O. Hancock, Nat. Rev. Mol. Cell Biol. 15, 615 (2014)
  • (30) R. H. Lee and C. S. Mitchell, J. Theor. Biol.(2015)
  • (31) J. T. Finer, R. M. Simmons, J. A. Spudich, et al., Nature 368, 113 (1994)
  • (32) J. Howard, Nature 389, 561 (1997)
  • (33) P. Pierobon, S. Achouri, S. Courty, A. R. Dunn, J. A. Spudich, M. Dahan, and G. Cappello, Biophys. J. 96, 4268 (2009)
  • (34) M. J. Saxton and K. Jacobson, Annu. Rev. Biophys. Biomol. Struct. 26, 373 (1997)
  • (35) J.-P. Bouchaud, J. Phys. I 2, 1705 (1992)
  • (36) A. Lubelski, I. M. Sokolov, and J. Klafter, Phys. Rev. Lett. 100, 250602 (2008)
  • (37) Y. He, S. Burov, R. Metzler, and E. Barkai, Phys. Rev. Lett. 101, 058101 (2008)
  • (38) I. M. Sokolov, Physics 1, 8 (2008)
  • (39) J. H. Schulz, E. Barkai, and R. Metzler, Phys. Rev. Lett. 110, 020602 (2013)
  • (40) T. Solomon, E. R. Weeks, and H. L. Swinney, Phys. Rev. Lett. 71, 3975 (1993)
  • (41) R. G. DeVoe, Phys. Rev. Lett. 102, 063001 (Feb 2009)
  • (42) E. Lundh and P. Ao, Phys. Rev. A 61, 063612 (2000)
  • (43) B. I. Shraiman, Phys. Rev. A 36, 261 (1987)
  • (44) T. Thoresen, M. Lenz, and M. L. Gardel, Biophysical Journal 104, 655 (2013)
  • (45) D. B. Allan, T. A. Caswell, and N. C. Keim, “Trackpy v0.2,” (May 2014), http://dx.doi.org/10.5281/zenodo.9971
  • (46) M. Weiss, M. Elsner, F. Kartberg, and T. Nilsson, Biophys. J. 87, 3518 (2004)
  • (47) F. Amblard, A. C. Maggs, B. Yurke, A. N. Pargellis, and S. Leibler, Phys. Rev. Lett. 77, 4470 (1996)
  • (48) J.-H. Jeon, V. Tejedor, S. Burov, E. Barkai, C. Selhuber-Unkel, K. Berg-Sørensen, L. Oddershede, and R. Metzler, Phys. Rev. Lett. 106, 048103 (2011)
  • (49) A. Reilein, S. Yamada, and W. J. Nelson, J. Cell Biol. 171, 845 (2005)
  • (50) J. L. Ross, H. Shuman, E. L. Holzbaur, and Y. E. Goldman, Biophys. J 94, 3115 (2008)
  • (51) O. Osunbayo, J. Butterfield, J. Bergman, L. Mershon, V. Rodionov, and M. Vershinin, Biophys. J. 108, 1480 (2015)
  • (52) M. W. Elting and J. A. Spudich, Developmental Cell 23, 1084 (2012)

Supplementary Information for

A cycling state that can lead to glassy dynamics in intracellular transport

Appendix A Alternative models

Here we illustrate that the results are robust for alternative formulations of the model. First, we tested a model with a one-step binding process. In this scheme, the probability of a motor at distance to be attached and active is

Here, is the maximum attractive energy that can be achieved by binding, relative to the unbound state. Fig. S1 compares the MSD results for this scheme and that in the main text.

Second, we explored alternative ways of moving the motors in response to the forces. We proposed the projection rule in the main text to ensure that motors walk parallel to tracks, as observed in many experiments. When this requirement is relaxed, the velocity at each time step can be written as , where is the number of binding sites attached to filament and is the unit vector in the direction of the filament. The position update then becomes . We recover aging and similar values for the MSD scaling exponent under this relaxed assumption (Fig. S2).

Figure S1: Mean-squared displacement for a Boltzmann binding scheme (colored curves), compared with results for the two-step attachment scheme in the main text (gray). (A) Time-averaged mean square displacement (MSD) for (green) and (blue). (B) The MSD as a function of measurement time for (orange) and (red). Other parameters are the same as those in the main text (Table 1).
Figure S2: Mean-squared displacement for the model with a net velocity (red and blue). Results for the two-step attachment scheme in the main text are shown in gray. The parameters used are the same as those in the main text (Table 1).

Appendix B Effect of nonergodicity on lag-time exponents

The scaling of the MSD with lag time is dependent on the duration of the recording (measurement time) . Since experimental data are often limited to short trajectories, we show the resulting time-averaged MSD exponents for shorter simulation times in Fig. S3.

Figure S3: MSD as a function of lag time for binding radii (gray), 0.01 (blue), 0.1 (orange), and 1(red) and duration of simulation , and (from top to bottom). The parameters used are the same as those in the main text (Table 1).

Additionally, we show the trajectories for different binding radii and temporal coarse-grainings (Fig. S4). When is larger than the average trapping time in a vortex, the time-averaged MSD is diffusive.

Figure S4: Example trajectories for binding radii , 0.01, and 0.1. The trajectories are coarse-grained by only showing points that are temporally separated by . The parameters used are the same as those in the main text (Table 1).
Figure S5: Mean-squared displacement for different filament network densities. (A) MSD as a function of lag time for network density (number of filaments per unit area) (red), (orange), (blue), and (green) and (black). (B) The MSD as a function of measurement time for the same densities as in (A). Other parameters are the same as those in the main text (Table 1).

Appendix C Experimental details

Figure S6: Mean-squared displacement for myosin II motors. (A) MSD as a function of lag time. (B) The MSD as a function of measurement time.

Networks of bundled actin were polymerized in a thin layer at a surface of a flowcell. A supported lipid bilayer passivated the surface. The supported bilayer was formed by incubating a UV-ozone cleaned borosilicate coverslip (Fisherbrand) with 1 mM vesicle suspension. Vesicles were prepared by the standard method of extrusion (200 nm and 50 nm pore membranes, Liposofast extruder, Avestin) from dried films of phospholipid (1,2-dioleoyl-sn-glycero-3-phosphocholine, Avanti Polar Lipids) suspended in buffer (140 mM NaCl, 8.5 mM NaHPO, 1.5 mM NaHPO, pH 7.5). After a complete bilayer formed, excess vesicle suspension was exchanged with actin polymerization buffer (10 mM imidazole, 50 mM KCl, 0.2 mM EGTA pH 7.5, 300 M ATP). Monomeric actin (rabbit skeletal muscle purified from acetone powder, Pel-freeze; 2.0 M unlabeled and 0.64 M labeled with the fluorophore tetramethylrhodamine-6-maleimide, Life Technologies) was added to initiate the polymerization of long, entangled actin. Depletion agent (0.3 wt %, 15 centipoise methylcellulose, Sigma) crowded actin to the surface. We used an oxygen scavenging system (50 M glucose, 0.5 vol % -mercaptoethanol, glucose oxidase, and catalase) to reduce photobleaching. After 30 min of polymerization, fimbrin (0.53 M, pombe) was added to crosslink the actin filaments into a network of bundles. Thick filaments of myosin II were polymerized in a similar manner by adding monomeric myosin II (20 nM, chicken skeletal muscle, fluorescently labeled with Alexa 647, Life Technologies) to a separate, actin-free solution. After 10 min of polymerization, myosin in ATP (2% of total sample volume) was gently mixed with the solution above the actin network, such that the final concentrations were 3.8 pM myosin and 2.3 mM ATP. Actin and myosin were imaged with a inverted microscope (Nikon Eclipse Ti-PFS) equipped with a spinning disk confocal head (CSUX, Yokogawa), 561 nm and 647 nm laser lines, 60/1.49 NA oil immersion objective (Zeiss), and a CCD camera (Coolsnap HQ2, Photometrics). Images, collected at 1.5 s intervals, began 10 min after myosin was added to the sample.

Single-particle trajectories were obtained using the Python-based implementation of the Crocker-Grier algorithm Trackpy. The settings are given in Tab. 2. The particle identification parameters were chosen such that all visible features were detected and the mask size was chosen large enough to obtain subpixel resolution (The histogram of the decimal values of the identified particles showed a flat distribution, indicating no bias within a pixel.). The linking values were set by looking at a previously manually tracked data set and using the largest displacement of that dataset as search range. All other parameters were set to the default values. We quantify the static tracking error for brighter and dimmer objects (implemented in Trackpy). The Fig. S8 shows the ensemble averaged tracking error as a function of time. There is no overall trend in the accuracy over the time of tracking in high (red) or low (blue) brightness fluorescent objects, and there is no statistical difference in the mean tracking error between the two populations (mean and standard deviation shown, unequal variances -test -value=0.0459)

Parameter Value
mask diameter 19
minimum mass 650000
noise size 3
search range 15
memory 1
Table 2: Particle tracking parameters for Trackpy.
Figure S7: Fluorescently labelled actin network. The images show snapshots of the actin network while myosin II motors are moving on it (not shown). While the frames show some photobleaching over the course of the experiments, the network does not remodel.
Figure S8: (A, B) Tracking accuracy for trajectories with high (red) or low (blue) fluorescence. The shaded area denotes the standard deviation. (C) Mean tracking accuracy for high and low fluorescence trajectories. Lines indicate standard deviations. (D) Tracking error for randomly selected trajectories over time.

Appendix D Analysis of trapped states

Here, we consider the possibility that detachment could give rise to the power-law dwell times evidenced by the decay of the MSD with the measurement time. If detachment were responsible for the anomalous dynamics, the particles should undergo simple diffusion when apparently trapped. To test for this possibility, we divided the trajectories between trapped and non-trapped periods and then analyzed their dynamics as follows.

To define the trapped portions of single-particle trajectories, we calculated the average position of each particle (motor) every 40 frames. We then compared the instantaneous positions with the averaged ones. If a point in a trajectory deviated by less than 2.0 pixel widths (182 nm) from the corresponding average position, we count that frame as trapped. The results are shown in Figure S9. Representative trapped trajectory segments are shown in red in the context of the non-trapped segments and the filament network (Figure S9A) and in a magnified view (Figure S9B).

To test for simple diffusion with a time resolution of a few frames, we employ the relative angle distribution introduced by Burov et al. Burov et al. (2013). The relative angle is defined by


, is the position at time , and is the lag time as in the present paper. Here we use frame. We build a histogram of these values for the trapped (red) trajectory segments and normalize such that it integrates to one.

The relative angle distribution for the trapped trajectory segments is shown in Figure S9C. If the motion were simple diffusion, the distribution would be flat, because there would be no directional correlation between steps of the random walk. Instead, we see that it is peaked at , which indicates that the particles are making reversals. The only way that this could happen from detachment would be if the particles were caged when off the filaments. However, this cannot be the case, as the prevalent switching between trapped and non-trapped states shows that the majority of the particles are unimpeded by obstacles. The only possibility consistent with Figure S9C is that the particles are switching direction while on the filaments, owing to the motor dynamics—this is the cycling mechanism that we propose.

Figure S9: The dynamics in apparent trapped states are not simple diffusion. (A) Representative single-particle trajectories colored by whether they are assigned to be trapped (red) or not (blue) according to the criteria in the accompanying text. The actin network is shown in gray in the background. (B) Magnified view of selected trapped trajectory segments. (C) Relative angle distribution of trapped trajectory segments for frame. Shaded area indicates standard deviation.


  • Burov et al. (2013) S. Burov, S.M.A. Tabei, T. Huynh, M.P. Murrell, L.H. Philipson, S.A. Rice, M.L. Gardel, N.F. Scherer,  and A.R. Dinner,  Proc. Natl. Acad. Sci. USA 110, 19689–19694 (2013).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description