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

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.

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.

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 i with a fixed unit vector e i . 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 [25].
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 M 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 M possible attachments for a motor among filaments with probability proportional to , where d i is the shortest distance between the motor and filament i, and s is a parameter that sets the interaction length scale. This probability is normalized by the sum i b i over all filaments within a distance 3s of the motor. Then, we determine if each such interaction exerts force to move the motor (henceforth, "active") with probability b i or not ("inactive") with probability 1 − b i . We denote the number of active attachments to filament i by k i . To determine the change in position of the motor we compute the vector v = i k i e i (orange in Fig. 2B). This choice is consistent with measurements that show that motor velocities increase with head numbers [30]. We project v 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 forcevelocity 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.
The values of the simulation parameters are given in Table I. 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 [31][32][33]. We associate our unit length with 1 µm, such that the average filament length is where x(t) is the position of a motor at time t. 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 [34]. It is important to note that the time-and ensemble-averaged MSDs can exhibit different scalings when the process of interest is non-ergodic [35][36][37][38][39].
By varying the binding radius s 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 (T ) in Fig. 3B. For an ergodic system, the MSD should be constant as T varies-the properties of the motion are the same independent of the length of recording. In contrast, a decrease in the MSD with T 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 [1] and the motion of insulin granules in vivo [1,9].
The observed aging exponent depends on the size of the binding radius s, 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 M (Fig. 3C and D). However, due to the non-ergodic nature of this transport process, the  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 [40]. Analogous observations exist for trapped ions [41] and Bose-Einstein-Condensates [42].
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 ∆ = 100, 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 I, 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 P f that l 1 out of M possible attachment sites are assigned to filament 1 and that the remaining are assigned to filament 2 is binomial: The probability P a that k 1 out of l 1 possible attachments are active is also binomial: and similarly for k 2 . Since filament assignment and selection of the active attachments are independent events, we can now write for the overall probability P of having k 1 and k 2 active binding sites on filaments 1 and 2 The limits of the sum are set by l 1 ≥ k 1 , l 2 ≥ k 2 , and l 1 + l 2 = M .
Associating filament 1 with the x direction and filament 2 with the y direction, the resulting master equation for motion in a square vortex is Here, P t (x, y) is the probability to be in a certain location (x, y) in the vortex at time t, 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.
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 [43]. However, the situation differs, in that the diffusion is position-dependent (i.e., P (k 1 , k 2 ) varies in space through b 1 and b 2 ). 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 CONTRIBU-TIONS
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 [44]. 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 particletracking 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 T shows aging, with an exponent comparable to previously published data (compare [1] 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 −0.6125 ± 0.0035 and −0.5112 ± 0.0086 (mean ± SEM high fluorescence and low fluorescence groups, respectively). We tested if the means are significantly different using the unequal variances t-test. The two-tailed p-value is less than 10 −4 , 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 (N = 22, p-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 44.1 ± 2.6 nm/s and 45.0 ± 2.5 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.

V. DISCUSSION
Motor-driven processes in cells often exhibit anomalous statistics [12,[46][47][48], and these statistics can have important functional consequences [9]. 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 [49]. A study that reconstructed these epithelial microtubule networks in vitro observed a kinesin-coated bead cycling through a vortex structure [50]. The same study (and others [51]) also saw a slowdown 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 [24] 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. Understanding how motor assemblies behave on complex filament networks in cells is an outstanding challenge [52]. 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 [25], is a useful direction for future research. It will also be interesting to understand the interplay of motor transport, force transmission, and network rearrangement.  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 [45].

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

VII. 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 p of a motor at distance d i to be attached and active is Here, E 0 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 then becomes x(t + dt) = vdt. We recover aging and similar values for the MSD scaling exponent under this relaxed assumption (Fig. S2).

VIII. 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) T . Since experimental data are often limited to short trajectories, we show the resulting time-averaged MSD exponents for shorter simulation times in Fig. S3.
Additionally, we show the trajectories for different binding radii and temporal coarsegrainings ∆ (Fig. S4). When ∆ is larger than the average trapping time in a vortex, the time-averaged MSD is diffusive.  (Table 1).  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  (Table 1).

Networks
the default values. We quantify the static tracking error for brighter and dimmer objects  (Table 1).

X. 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).
v(t; ∆) = x(t + ∆) − x(t), x(t) is the position at time t, and ∆ is the lag time as in the present paper. Here we use ∆ = 1 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.