Transport of synaptic vesicles is modulated by vesicular reversals and stationary cargo clusters

Neuronal cargo, such as membrane-bound organelles and macromolecular complexes, are transported by molecular motor proteins on microtubule tracks along neuronal processes. Individual motor proteins predominantly move in a unidirectional manner, depending on the intrinsic polarity of microtubules. However, membrane-bound organelles, such as mitochondria, endosomes, lysosomes and synaptic vesicles, exhibit bidirectional transport when viewed over long timescales, indicating that both anterograde and retrograde motors participate in transporting neuronal cargo (Maday et al., 2014). Although endosomes (Loubéry et al., 2008) and mitochondria (Narayanareddy et al., 2014) exhibit diverse motility characteristics based on their maturation state and size, respectively, much less is known about diversity in synaptic vesicle motility in vivo.

Several studies have found that stationary vesicle clusters (SCs) of neuronal cargo are a prominent feature of axonal transport (Che et al., 2016; Ganguly et al., 2015; Kang et al., 2008; Lasiecka et al., 2010; Sood et al., 2018). These SCs can modulate cargo transport, causing them to pause at these locations (Sood et al., 2018) or slow down in their vicinity (Che et al., 2016). SCs of synaptic vesicles appear to be associated with actin-rich regions (Sood et al., 2018) and microtubule ends in vivo (Che et al., 2016; Yogev et al., 2016), and can mobilise in response to neuronal activity (Sood et al., 2018). As SCs can both sequester and release vesicles in a context-dependent manner, their effect on cargo transport is not well understood. Furthermore, it is unclear how neurons modulate the relative fractions of stationary and moving cargo within axons, and whether such modulation is necessary to maintain healthy neuronal function. As it is experimentally infeasible to specifically perturb only SCs or moving vesicles without affecting several other features of axonal transport, we attempted to address these questions by constructing simulation models.

Previous computational studies have examined individual obstacles to transport, such as cargo–cargo interactions (Lai et al., 2018) and microtubule ends (Gramlich et al., 2021). One such theoretical study used a simple lattice model with paused states to study the formation of cargo pileups and axonal swellings around microtubule tracks in diseased neurons (Lai et al., 2018). They showed that higher frequencies of reversals and track switching events reduce cargo pileups, but did not comment on the effect of these perturbations on net cargo flow (Lai et al., 2018). Another study used a single rate constant model, in conjunction with experiments conducted in Caenorhabditis elegans motor neurons, to demonstrate that pausing of cargo at microtubule ends is regulated by the nature of the motor and the direction of transport (Gramlich et al., 2021). They propose that reversals likely contribute very little to how vesicles navigate microtubule ends (Gramlich et al., 2021). Thus, the mechanisms regulating SC formation and how moving vesicles navigate such obstacles is still unclear. We attempted to address this using our model, which accounts for a diversity of obstacles, similar to what is observed in healthy neuronal processes in vivo.

We hypothesised that SCs limit the transport of vesicles and that vesicle motility characteristics promote cargo flow by modulating the proportion of SCs formed along neuronal processes. This study combined experimental observations of in vivo pre-synaptic vesicle transport in the touch receptor neurons (TRNs) of C. elegans with kinetic Monte Carlo simulations of a model constructed to describe our experimental observations. Our model accounts for cargo crowding owing to the presence of both static and dynamic obstacles to transport. Vesicles can pause at microtubule ends, at sites of stalled cargo such as mitochondria, and upon encountering other moving or stationary vesicles using the same tracks for transport. The model further allows for microtubule track switching and bidirectional transport (reversals) of cargo. The model parameters were benchmarked extensively to our experimental observations. Using this model, we addressed (1) whether vesicular motility characteristics promote net anterograde flow of cargo in crowded environments and (2) how SCs modulate the transport of vesicles.

Our simulations and experimental data together demonstrate that there is a reciprocal relationship between SC formation and net anterograde cargo transport. The ability of moving vesicles to switch between microtubule tracks and reverse their direction of motion helps maintain net anterograde cargo flow in crowded environments, with reversals playing a much more significant role than track switching. The model further predicts that increasing the occurrence of reversals at the sites of obstacles significantly improves net anterograde transport, suggesting that reversals can aid in navigation of crowded geometries. Our simulations support the view that SCs function as dynamic reservoirs of cargo vesicles, and that vesicle motility characteristics such as reversals regulate cargo transport by modulating the proportion of SCs along the neuronal process. These predictions from simulations were validated by experiments conducted on a C. elegans tauopathy model, with implications for the progression of neurodegenerative diseases.

As we wanted to investigate the role of vesicular motility characteristics in navigating crowded environments in vivo, it was important to examine their motion properties in detail. Thus, we characterised the transport properties of the precursors of synaptic vesicles (pre-SVs) in the processes of C. elegans TRNs in vivo. The pre-SVs were fluorescently labelled by tagging GFP to the synaptic vesicle protein RAB-3 (Bounoutas et al., 2009; Kumar et al., 2010) and their motion along the neuronal process was tracked using live imaging. Moving pre-SVs predominantly exhibited unidirectional motion in the anterograde and retrograde directions (Fig. 1A). Anterogradely-moving pre-SVs comprised two distinct populations, one with a very high pause frequency (>0.15 s^–1), constituting 16.67% of all anterograde vesicles, and the other with a lower pause frequency (<0.15 s^–1) (Fig. 1B). We categorised anterogradely moving pre-SVs exhibiting <0.15 pauses s^–1 as ‘smooth anterograde’ (SmA), and pre-SVs showing >0.15 pauses s^–1 as ‘step anterograde’ (StA) (Fig. 1B). Retrogradely moving pre-SVs (R) predominantly exhibited very low pause frequencies (<0.15 s^–1) along their trajectories, with <2.5% of all retrograde vesicles exhibiting pause frequencies >0.15 s^–1 (Fig. 1C). As the proportion of retrogradely moving vesicles exhibiting pause frequencies of >0.15 s^–1 was very low, we considered the retrogradely moving vesicles to exhibit predominantly smooth motion.

Analysis of the motion properties of these distinct populations showed that StA vesicles had a significantly higher pause duration fraction (Fig. 1F) and significantly lower segment run lengths compared to those of the SmA and R populations of moving pre-SVs (Fig. 1D). The motion of SmA pre-SVs was processive (Fig. 1A), with an average segment run length of 1.13±0.75 μm (±s.d.) (Fig. 1D), a maximum run length of 7.87 μm and a pause duration fraction of 0.3±1.8% (mean±s.d.) (Fig. 1F). StA motion, on the contrary, was characterised by a higher pause duration fraction (7±3.6%) (Fig. 1F) and shorter segment run lengths (0.7±0.48 μm) compared to those of SmA (Fig. 1D). StA vesicles also showed reduced net displacement over time compared to that of SmA and R vesicles (Fig. 1D). R vesicle motion was characterised by an average segment run length of 1.63±1.53 μm (Fig. 1D), a maximum run length of 14.32 μm and a low pause duration fraction (0.6±2.3%) (Fig. 1F). In summary, StA pre-SVs exhibited shorter segment run lengths, more frequent and longer pauses, and lower mean velocities compared to those of SmA and R pre-SVs.

Although pre-SVs predominantly exhibited unidirectional transport in either the anterograde or retrograde directions, we observed that a small proportion of anterograde and retrogradely moving pre-SVs changed their direction of motion, termed a ‘reversal’ (Fig. 1A). On average, 13.68% (1.368±0.47717 per ten RAB-3 pre-SVs) of all moving RAB-3 pre-SVs reversed along the posterior lateral mechanosensory (PLM) process (Fig. 1G). Additionally, a population of pre-SVs were observed to remain stationary for a broad distribution of timescales, with several being stalled for over 2–3 min. Based on the duration for which pre-SVs remained stationary, we categorised SCs as short-lived SCs (between 15 and 45 s) and long-lived SCs (>45 s) (described in detail in the Materials and Methods). We observed that 70% of all SCs were long-lived (Fig. 1H), with a higher density along the neuronal process compared to that of short-lived SCs (Fig. 1I).

These stationary pre-SV clusters are a general feature of neuronal cargo transport and can act as crowded locations along the process (Sood et al., 2018). We propose that one or more of the vesicular motility characteristics described above might help maintain cargo transport in such crowded locations. As it is experimentally infeasible to perturb specific populations of vesicles without disrupting vesicle transport in general, we addressed this hypothesis using simulation models. We used the vesicle transport properties calculated above to inform our model described below.

We developed a coarse-grained model to examine the relationship between pre-SV motility characteristics and SC formation. We did not encode individual molecular processes such as motor attachment and/or detachment from microtubules or to pre-SVs as independent rates in our model. The rates incorporated in the model were based on experimentally derived frequencies of observable events. All these parameters involved multiple molecular processes, such as motor attachment and/or detachment to the vesicle surface or to microtubules, and tug-of-war- or coordination-based mechanisms for bidirectional transport. We did not include each of these processes as individual parameters in our simulations as experimental data on the properties of UNC-104 (a kinesin-3 motor) and the dynein motor complex that are functional in the PLM neuron are not available. Additionally, we did not explicitly model the three-dimensional diffusion of pre-SVs detached from the track as diffusion does not appear to play a significant role in cargo flow along the narrow confines of the PLM neuronal process. This was evident from our fluorescence recovery after photobleaching (FRAP) experiments, wherein all the fluorescence signal recovery observed in the photobleached region over a timescale of 3 min was associated with active vesicular transport (Fig. S2A,B).

We approximated the PLM neuronal process as a set of parallel tracks representing the peripheral microtubules, along which most cargo transport likely occurs (Cueva et al., 2007). Individual tracks comprise a one-dimensional lattice of 1000 sites, spanning a width of 8 nm each (Fig. 2A). Each particle in the system represents a pre-SV motor complex in one of the three motion states observed in vivo (Fig. 1A). Their relative proportions were benchmarked to experimentally observed ratios, as described in detail in the Materials and Methods. These pre-SVs, depending on their motion state, can move in the anterograde or retrograde directions (Fig. 2A) at rates (see Table S2) corresponding to the respective mean velocities of cargo observed in vivo (Fig. 1D). The model accounts for both static and dynamic obstacles to transport. Moving pre-SVs can stall upon encountering other moving or stationary vesicles on the same track, which act as dynamic obstacles to transport (Fig. 2A,D). Pre-SVs also stall when they encounter static obstacles such as microtubule ends or large stationary mitochondria (Fig. 2A,D). The frequency of dynamic obstacles to transport depends on the vesicle density, which was benchmarked to experiments (Fig. S2). The distribution of static obstacles encoded in the model was in keeping with reported experimental measurements, as described in detail in the Materials and Methods. A stalled pre-SV can change its position by track switching, or its state by reversing its direction of motion (Fig. 2B,C). Track switching (Fig. S3) and reversal rates were benchmarked to experiments using multiple metrics, as described in the Materials and Methods. The model was simulated using a kinetic Monte Carlo algorithm starting with a random placement of particles. The system was evolved until it reached the steady state, following which the locations and motion states of all pre-SVs and measurements of net anterograde current and proportions of long-lived SCs were recorded at intervals corresponding to rates of 5 fps used in experiments.

We did not explicitly encode cargo pausing as an independent parameter in the model, unlike previous studies (Lai et al., 2018; Gramlich et al., 2021). This allows cargo pausing to emerge dynamically based on the assumptions in our model, and we found that this was sufficient to predict the distribution of SC lifetimes observed in experiments. Our model was extensively benchmarked to observable features of pre-SV transport in the C. elegans PLM neuron. This allowed us to examine how multiple kinds of obstacles influence vesicle transport in healthy neurons, which has not been addressed before. This further provided us with a system with which to test and validate the predictions from our model, with links to neuronal function and the physiology of the animal.

In order to test the predictive power of the model, we examined its ability to predict changes to an unseen experimental dataset in response to a novel perturbation. We simulated axotomy and compared the change in vesicle density profiles as predicted from simulations with experimental observations of axotomy conducted in animals carrying a different fluorescently labelled synaptic vesicle marker. We ablated the PLM neurons of transgenic animals in which pre-SVs were fluorescently labelled by tagging eGFP to the transmembrane synaptic vesicle protein SNB-1 (Fig. 3A). As a consequence of continuous transport, pre-SVs began to pile up at the proximal cut site (process end closer to the neuronal cell body) following axotomy (Fig. 3A). We considered axotomy in our simulations as being the equivalent of introducing a permanent block to vesicle motion across all microtubules, separated from the distal process by a physical distance of 0.8 μm (100 sites), and removing all vesicles and microtubules that lay within the ablated region.

We specifically monitored the accumulation of pre-SVs at the proximal cut site 5 min and 10 min post axotomy (Fig. 3A). A visual comparison of vesicle accumulation profiles that were predicted by simulations with experimental observations indicated that our simulation model was accurately able to capture the vesicle density profile observed in experiments 5 min and 10 min post axotomy (Fig. 3B,C; Fig. S4A,B). A quantitative comparison showed that simulations reproduced the experimentally determined exponential decay in vesicle density with increasing distance from the cut site, and produced decay constant values that were similar to experimental measurements (Fig. 3B,C; Fig. S4A–D). The simulations were also able to accurately predict the change in the vesicle density profile at the proximal cut site observed between 5 min and 10 min post axotomy in experiments (Fig. 3D). Another key observation from our experiments and simulations was the increase in the width of vesicle accumulation at the proximal cut site between 5 and 10 min post axotomy (Fig. 3A; Fig. S4F). We found that the simulations reproduced the fold change in the extent of vesicle accumulation observed in experiments between 5 and 10 min post axotomy (Fig. 3E; Fig. S4E). In summary, we found that the predictions from the model were quantitatively similar to our experimental observations. This suggests that the model is an accurate representation of synaptic vesicle transport in the PLM TRN process and can be used to make predictions on the effect of perturbations in specific motility characteristics on pre-SV transport.

We examined the effect of eliminating track switching and reversals on SC lifetimes and net anterograde current (Fig. S5A–D). As reported in prior theoretical studies (Lai et al., 2018), we found that a reduction in the track-switching rate caused an increase in SC formation and decrease in net anterograde current (Fig. S5C,D). Interestingly, our simulations showed that removing reversals had a much more severe impact on pre-SV transport compared to the impact of removing track switching. For the case in which reversals were completely disallowed, the model predicted that all moving pre-SVs eventually come to a halt to form a permanent block of long-lived SCs (Fig. 4A; Movie 1) and the net anterograde current falls to zero (Fig. S5D). The observed increase in SC formation upon reduction in reversal rates is consistent with prior simulation work on axonal cargo accumulations (Lai et al., 2018). Starting from such a jammed configuration, allowing for vesicular reversals dissolved the permanent block, restoring net anterograde transport (Movie 2). Thus, our simulations predict that vesicular reversals play an important role in maintaining cargo flow and limiting SC formation.

We next asked whether net anterograde current and SC lifetimes are sensitive to the number of reversals occurring in the system. The schematic shown in Fig. 4A describes the state interconversion rates that are altered to reduce the total number of reversals in the simulations (Fig. S5E). The representative kymographs (Fig. 4A) reflect the changes observed in axonal transport upon reduction in the rate of reversals. We found that when the rates of reversals (StA to R and R to StA) were reduced to one-third, one-sixth and one-ninth the regular rate (Fig. S5E), there was a steady and consistent increase in the proportion of long-lived SCs (Fig. 4B,C) and a decrease in net anterograde transport (Fig. 4B,C).

To validate the predictions from our simulations, we examined a tauopathy model of C. elegans, bkIs10, which expresses the IN4R isoform of human tau (encoded by MAPT) carrying the V337M frontotemporal dementia with Parkinsonism chromosome 17 type (FTDP-17) mutation. C. elegans with bkIs10 form tau aggregates along the neuronal process, have uncoordinated movements and show age-dependent loss of axons (Kraemer et al., 2003).

We found that at younger stages [larval stage (L) 4 and 1-day stage], the stationary cargo density (Fig. 5A), proportion of long-lived SCs (Fig. 5B) and net anterograde current (Fig. 5C) of pre-SVs in bkIs10 animals were comparable to those of wild-type animals (Movies 3 and 4). Among 3-day-old and 5-day-old adults, bkIs10 animals had a significantly higher density of SCs compared to that of wild-type animals (Fig. 5A; Movies 5 and 6). Thus, the tauopathy model at older stages presents a system that is more densely crowded with SCs than in the wild type, which can be used to investigate the correlation between net current, the proportion of long-lived SCs and the reversal rate in vivo. We observed that upon increased crowding with age, bkIs10 animals showed a significant increase in the proportion of long-lived SCs (Fig. 5B) and a slight, but consistent, decrease in net anterograde current (Fig. 5C). The reciprocal relationship between long-lived SCs and net current predicted by simulations (Fig. 4C) was validated by our experimental observations (Fig. 5B,C,E).

Upon increased crowding with age, bkIs10 animals showed a consistent reduction in reversal rate as compared to that of wild-type animals (Fig. 5D). Like net current, the reversal rate exhibited a reciprocal trend with the proportion of long-lived SCs (Fig. 5E), suggesting that increased crowding is associated with a reduction in net anterograde current, reduced reversals and an increased proportion of long-lived SCs. This is consistent with the predictions from our model. Interestingly, our experiments showed that the decrease in reversal rate and increase in the proportion of long-lived SCs (at the 3-day adult stage) preceded the decrease in net anterograde current (at the 5-day adult stage) (Fig. 5B–D; Movies 5 and 6). This further corroborates the prediction from simulations that a reduction in reversal rates leads to decreased net anterograde transport.

As reversals play an important role in maintaining net anterograde transport and limiting SC formation, it is interesting to investigate whether reversals can serve as a useful mechanism to navigate obstacles to transport. A previous study conducted on DA9 motor neurons of C. elegans reported that the frequency of reversals occurring at microtubule ends is too low for it to contribute significantly to navigating microtubule ends (Gramlich et al., 2021). We asked whether increasing the numbers of reversals occurring at sites of obstacles could improve net anterograde transport.

Our experimental analysis showed that reversals occurred both at and away from SCs (crowded locations) along the PLM neuronal process (Fig. 6A). Nearly 40% of all the reversal events occurred at sites of SCs, with the rest occurring away from these sites (Fig. S5H). The average ratio of reversals occurring at SCs to away from SCs was 1:1.4 (Fig. S5I). Simulations run at a ratio of 1:1 for reversals occurring at and away from obstacles showed similar values for the fraction of area occupied by SCs over time (Fig. S5K) and net current (Fig. S5K) compared to the values obtained from experiments.

We simulated the following limits: (1) reversals occurring predominantly at the sites of obstacles, (2) reversals occurring uniformly along the neuronal process and (3) reversals occurring predominantly away from the sites of obstacles. The ratio of reversals occurring at and away from obstacles was varied as 100:1, 10:1, 1:10 and 1:100. In this system, net anterograde currents were the largest when reversals were allowed to occur primarily at the sites of obstacles, as shown for 100:1 (Fig. 6B). The net current observed in this configuration was nearly double the value observed at 1:1 (Fig. 6B). As this ratio was decreased, the net anterograde current reduced, whereas the proportion of long-lived SCs increased (Fig. 6C).

This prediction indicates that reversals can play a significant role in navigating crowded locations within the neuronal process. Our experimental observations, however, indicated that reversals occurred at a low frequency (Fig. 1G) and were distributed uniformly along the process (Fig. S5I). Thus, the PLM neuron is not optimised for maximum anterograde transport of synaptic vesicles. This suggests that a system that allows for the occurrence of reversals is sufficient to maintain net anterograde cargo transport at levels required for proper neuronal function. Furthermore, as we cannot currently perturb or trigger reversals specifically at crowded locations experimentally, this provides us with a testable prediction for future experiments as the molecular mechanisms governing bidirectional transport is a field of active research.

The narrow and highly confined processes of C. elegans TRNs present specific obstacles to cargo transport. Stationary vesicle clusters, which are a prominent feature of axonal transport in vivo, are crowded locations (Sood et al., 2018). At these locations, moving vesicles stall for timescales ranging from a few seconds to several minutes (Sood et al., 2018). As crowded locations pose impediments to transport by causing vesicle accumulation, potentially choking off steady vesicular flux, it appears reasonable that neurons might have evolved strategies to cope with the consequences of crowding in a dense environment. Increased SCs of vesicles and low vesicular flux are prominent features of a number of neurodegenerative diseases. Understanding the in vivo mechanisms underlying the formation and dissolution of such SCs might provide insights into disease states (Gunawardena and Goldstein, 2001; Kreiter et al., 2018).

Previous modelling studies of intracellular transport in crowded environments involved studying the effect of molecular crowding arising from interactions between kinesin-1 motors on their motion properties (Rank and Frey, 2018), strategies used by motor proteins to avoid motor density-induced traffic jams (Leduc et al., 2012) and traffic jams in pathological conditions, for instance, at regions of microtubule polarity mismatch (Kuznetsov, 2010; Kuznetsov and Avramenko, 2009a) or reduced microtubule densities (Kuznetsov and Avramenko, 2009b). Studies looking at the influence of obstacles to cargo transport have examined individual obstacles such as microtubule ends in healthy neurons (Gramlich et al., 2021) and off-track cargo pileups in axonal swellings of diseased and/or congested neurons (Lai et al., 2018).

Prior to studying specific disease states, it is important to examine cargo transport in healthy neuronal processes. Our model examines multiple obstacles to transport in healthy neurons, which allows us to investigate how vesicle clusters might form and persist for a broad range of timescales while still allowing for steady cargo transport. Although detachment from microtubules and diffusion is one possible solution to bypass crowded locations, diffusion-based processes are expected to play a subordinate role in narrow, cytoskeleton-dense neuronal processes (Fig. S2A,B). It has been previously reported that SCs, although long-lived, are dynamic, allowing for mobilisation of vesicles from such clusters (Sood et al., 2018). Thus, in narrow and dense neuronal processes, moving vesicles likely employ active transport-dependent mechanisms to maintain transport in crowded locations.

We constructed a simulation model for axonal transport of pre-SVs, which allowed us to investigate the contribution of various parameters to net anterograde cargo flow. The simulation model was benchmarked to the experimentally derived ratio of motion states and average velocities of moving pre-SVs in C. elegans TRNs. Furthermore, we systematically varied multiple simulation parameters, narrowing them to a set of conditions that produced transport behaviour that was consistent with in vivo observations, using different metrics and assays. The simulations included two active transport-dependent mechanisms to navigate crowded locations, track switching and reversals.

Reversals provide both a temporal and spatial relaxation of vesicular crowding, allowing a vesicle that has undergone reversals to return to the same region of the axon or microtubule track later. Reversals also aid in redistributing cargo to avoid pileups at crowded locations and reduce the frequency of vesicle–vesicle interactions (Wu et al., 1998; Sabharwal and Koushika, 2019; Lai et al., 2018). Reducing the rate of reversals in our simulations caused a significant increase in the proportion of long-lived SCs and a reduction in net anterograde flow (Fig. 4B,C). We also observed that SC formation and net anterograde flow were sensitive to the locations of reversals (Fig. 6B,C). In both these perturbations, there was a consistent reciprocal trend observed between long-lived SC formation and net anterograde cargo flow (Figs 4C and 6C). This suggests that SCs act as dynamic reservoirs for moving vesicles, wherein the capture of moving vesicles leads to an increase in cluster formation and the mobilisation of vesicles from these clusters leads to increased cargo flow.

The reciprocal relationship between long-lived SCs and net anterograde current predicted by simulations was validated by experiments (Fig. 5E). In a C. elegans model for tauopathy, we observed that as the animal ages, the PLM neuronal process became increasingly more crowded than for the wild type (Fig. 5A). As crowding increased, there was an increase in the proportion of long-lived SCs (Fig. 5B) and a decrease in net anterograde current (Fig. 5C). This observation corroborates the predictions from our simulations, indicating that there exists a trade-off between SC formation and cargo flow. SCs in healthy neurons likely function to limit cargo transport. Thus, it is important for neurons to optimise these two parameters, as an imbalance between the two is a characteristic feature of several neurodegenerative disease states (Gunawardena and Goldstein, 2001; Baloh et al., 2007; Bilsland et al., 2010; Collard et al., 1995; Ebneth et al., 1998; Mandelkow et al., 2003; Pigino et al., 2003; Sorbara et al., 2014; Stamer et al., 2002).

Our experimental observations further corroborate the prediction from simulations that reversal rates might modulate net cargo flow. As the neuron became increasingly crowded with age, there was a steady decline in reversal rate observed in the tauopathy models compared to that observed in the wild type (Fig. 5D), and there was a reciprocal trend observed between the reversal rate and proportion of long-lived SCs (Fig. 5E). This followed the same trend observed for net anterograde current (Fig. 5C,E). However, we found that the reversal rate started to decline significantly at the 3-day adult stage (Fig. 5D), whereas the net current declined later at the 5-day adult stage (Fig. 5C). Additionally, wild-type animals independently showed a decline in stationary cargo density with age and a corresponding increase in reversal rate (Fig. 5A,D). These observations, coupled with the predictions from simulations, suggest that the lack of modulation in reversal rates in the tauopathy model causes the decrease in net anterograde flow and increased long-lived SC formation, thereby contributing to disease progression (Mandelkow et al., 2003; Kreiter et al., 2018; Roy et al., 2005; Gunawardena and Goldstein, 2001; Sabharwal and Koushika, 2019).

An important prediction from our model is that increasing the number of reversals at sites of obstacles causes a significant increase in net anterograde current (Fig. 6B). This suggests that reversals can be a powerful means to navigating crowded locations. The ability to trigger reversals specifically at such locations could help improve cargo transport in congested or diseased neurons. In summary, long-lived SCs can function as dynamic reservoirs that modulate net anterograde cargo flow, and reversals regulate the relative proportions of SCs and moving vesicles. An important implication of our findings is that the age-dependent decline in axonal transport observed in neurodegenerative disease models could be driven by changes in vesicle motility characteristics. Studying these changes might provide a new window into a class of neurodegenerative diseases and their associated therapeutics.

All C. elegans strains used in this study were maintained at 20°C on nutrient growth media (NGM) media seeded with Escherichia coli OP50, as described in the standard protocol (Brenner, 1974). The animals were primarily imaged at the L4 stage. In some experiments, animals were imaged at the 1-day, 3-day and 5-day adult stages; this is indicated in the relevant sections. All animals used for imaging came from non-contaminated, growing plates. The strains used in this study are listed in Table S1.

Live single worms were mounted on glass slides with 5% agar pads and anaesthetised using 5 mM tetramisole (Sigma-Aldrich, St. Louis, MO, USA) prepared in M9 buffer [3 g KH₂PO₄, 6 g Na₂HPO₄, 5 g NaCl, 1 ml 1 M MgSO₄, made up to 1 l with H₂O; sterilized by autoclaving (Stiernagle, 2006)]. All time-lapse imaging was conducted on the Yokogawa CSU-X1-A3 spinning disc, using the Hamamatsu ImagEM C9100-13/14 EMCCD camera integrated with an Olympus IX83 microscope by PerkinElmer. Samples were imaged using a 100×1.63 NA objective. Some imaging was also done on the Zeiss LSM 880 using an 63×1.4 NA objective. Imaging was performed at 5 fps and 3–5 min movies were acquired for all genotypes.

All image panels used for representation and analysis of time-lapse movies were generated using Fiji-ImageJ v1.52p. Experimental kymographs were generated using the MultipleKymograph plugin. Plugins were downloaded from the NIH website (a modified version of the plugin called KymoResliceWide is available at https://imagej.net/plugins/kymoreslicewide) or http://www.emblheidelberg.de/eamnet/html/bodykymograph.html.

In kymographs, cargo moving in the retrograde direction (towards the cell body) and anterograde direction (away from the cell body) appear as sloped lines, whereas stationary cargo appear as vertical lines. A cargo is counted as moving if it has been displaced by at least 3 pixels in successive time frames. For all direct comparisons with simulations, experimental kymographs were analysed in 8 μm sections as all simulations spanned 8 μm regions.

Moving vesicles were manually annotated on experimental kymographs (Fig. S1A), with individual trajectories traced as accurately as possible, using the ‘Segmented Line’ feature in Fiji. Motion parameters such as pause duration (Fig. S1B), mean velocity (Fig. S1C) and segment run length (Fig. S1C) were calculated using an in-house macro (‘Macro for motion properties.ijm’) written in Fiji. The codes for the calculation of motion properties can be found on GitHub (https://github.com/amruta2612/Fiji-macros). Segment run length was defined as the distance moved by the cargo between pauses (Fig. S1C), using a conversion factor of 0.13 μm/pixel. Pause time was defined as the length of time for which a cargo stayed stationary between two consecutive runs (Fig. S1B). The ‘Measure’ function was used to obtain the length of these lines in pixels, which was then converted to units of time using a conversion factor of 1 pixel=0.2 s (corresponding to 5 fps). Pause duration fraction was defined as the ratio of time spent by a vesicle in the paused state divided by the total duration of motion of the vesicle (Fig. S1B). The mean velocity of a moving vesicle was calculated by dividing the net displacement of the vesicle by the total duration of its motion (Fig. S1C).

All vertical lines on the kymograph represent stationary or paused cargo. The entire length for which each particle remained stationary was annotated using the ‘Segmented Line’ feature in Fiji (ImageJ). The length of each line was calculated in pixels using the ‘Measure’ tool of ImageJ, which was then used to calculate the time for which the vesicle remained stationary (in seconds) based on the frame rate at which the movie was captured. A vesicle was counted as stationary if it did not move for more than 15 s (for synaptic vesicles). This cut-off was based on the average pause times of moving vesicles and any vesicle was said to be stationary if it remained paused for at least five times longer than its average pause time. All stationary vesicles were further binned as short lived or long lived based on the duration for which they remained stationary. For this, we tracked two types of stationary vesicles in our kymographs: (1) stationary vesicles, for which the end and origin were within the time window of the entire kymograph analysed, and (2) stationary vesicles for which the end or origin were not within the time window of the kymograph analysed and thus the exact duration for which they remained stationary was not known. All vesicles under type 1 that remained stationary for 15–45 s were binned into the short-lived category and those beyond 45 s (45 s–1 min or 1–5 min) were binned into the long-lived category. All vesicles under type 2 were binned into the long-lived category as long as they met the minimum cut-off of 15 s. All stationary cargo were annotated manually and were counted and binned into the relevant categories using a macro (‘Macro for SC lifetime.ijm’) written in ImageJ (https://github.com/amruta2612/Fiji-macros).

Stationary cargo were identified as described in the previous section. In kymographs, we marked each SC using a rectangular region of interest (ROI), with the length representing the lifetime of the SC (in seconds) and the width representing the distance occupied by the SC along the axonal process (in micrometres). The area of this rectangular ROI thereby represents the total distance along the axon that was occupied by the respective stationary cargo over time. In order to quantify and compare stationary cargo site occupancy over time, the sum of the areas of all rectangular ROIs was plotted as a percentage of the total area of the kymograph.

Sloped lines in the kymograph indicate moving vesicles. The minimum cut-off used for counting a sloped line as a moving cargo was a displacement of 3 pixels in successive time frames. All vesicles moving away from the cell body were counted as anterograde flux, whereas the ones moving towards the cell body were counted as retrograde flux. Total flux was calculated as the sum of anterograde and retrograde flux, whereas net current was calculated as the difference between anterograde and retrograde flux normalised to time, expressed as number of vesicles/(8 μm×s). All moving vesicles were annotated manually and were counted and binned into the relevant categories using a macro (‘Macro for flux.ijm’) written in ImageJ (https://github.com/amruta2612/Fiji-macros).

Reversals refer to changes in the direction of motion of moving vesicles from either anterograde to retrograde direction or vice versa. An anterogradely moving vesicle that changed its direction to retrograde along its trajectory was counted as an anterograde reversal. Similarly, any vesicle moving in the retrograde direction switching to the anterograde direction was counted as a retrograde reversal. Reversal data for all experiments is represented as a fraction of moving vesicles in a particular movie, i.e. the number of reversals per ten moving vesicles. For anterograde reversal rate calculation, the number of anterograde reversals was divided by anterograde flux; for retrograde reversal rate calculation, the number of retrograde reversals was divided by retrograde flux; and for total reversal rate, the total number of reversals was divided by total flux.

The locations of reversals were divided into two categories: (1) at sites of stationary vesicles – all reversal events that occurred on encountering a stationary vesicle or within the proximity of 2 pixels (0.26 μm) on either side of stationary vesicles were counted in this category; and (2) away from stationary vesicles – all reversal events that occurred at a distance >2 pixels (0.26 μm) from stationary vesicles were counted in this category. These data were then represented as the percentage of reversals taking place at or away from an SC. As stationary vesicles were bright, these cut-offs were chosen keeping in mind the diffraction limits of the microscope and accuracy of measuring the distance from stationary vesicles.

For each categorical variable, such as the proportion of long-lived SCs, net current and rate of reversals, all data points from both wild type (jsIs821) and the tauopathy model (bkIs10;jsIs821) were pooled together and sorted in ascending order of values. The maximum and minimum values were noted and each data point was transformed using min-max normalisation: ⁠, such that all the transformed data points lay between 0 and 1. This transformation was applied independently to data for each categorical variable and the final transformed data points were color coded according to their respective genotypes.

In order to quantitatively determine the difference between two distributions, we calculated the first Wasserstein distance or the Earth Mover's distance (EMD). We implemented this by converting both distributions to their empirical cumulative distributions or by transforming the distributions using min-max normalisation to scale between 0 and 1. We then used the ‘emdist’ package in R (https://CRAN.R-project.org/package=emdist) to compute the EMD value associated with the pair of distributions being compared. A higher EMD value represents a greater difference between the two distributions.

Photobleaching experiments were performed on an Olympus FV3000 confocal microscope using a 40×/1.3 NA oil objective, in which a 10 μm region of the PLM neuronal process near the cell body was imaged using a 488 nm laser for 1 min before bleaching, and for 3 min immediately after bleaching, to monitor fluorescence recovery over time. The rectangular region spanned the entire width of the neuronal process; hence, recovery could only occur lengthwise. The analysis of bleach recovery was conducted as follows. In each movie, the 10 μm region that was photobleached was marked as the ‘bleached region’, a 10 μm region excluding the neuronal process (chosen arbitrarily) was selected as the background, and a 10 μm region of the neuronal process excluding the bleached region (chosen arbitrarily) was selected as the ‘unbleached region’. The mean fluorescence intensity across each 10 μm region was monitored over time. The background-corrected intensity for each timepoint was calculated as ⁠. The background-corrected intensity was then smoothed using the Whittaker–Henderson smoothing function with λ=1600. The smoothed intensity profiles were internally normalised to the mean pre-bleach intensity. All metrics such as the half-time (t_1/2) and recovery ratios were calculated from the smoothed normalised intensity profiles.

Axotomy experiments were performed using a 355 nm (ultraviolet) pulsed nanosecond laser (Nano S60-30 Nd:YAG Litron laser). Live animals were anaesthetised using 5 mM tetramisole (Sigma Aldrich) in M9 buffer on 5% agarose pads. Slides were imaged on an Olympus IX73 inverted epifluorescence microscope using a UPLSAPO 100×1.4 NA oil immersion objective, with an X-cite Mercury vapour short arc lamp at 100% iris setting, and an Andor EMCCD camera operated at a gain of 50 and an exposure of 200 ms. The analysis for vesicle accumulation following axotomy was conducted as follows. For each animal or movie, static images were denoised using the ‘Gaussian blur’ feature with a Σ(radius)=8 in Fiji/ImageJ. In the denoised images, the fluorescence intensity over a 1 μm region at the proximal and distal cut sites was monitored at 5 mins and 10 min post axotomy, using the rectangular ROI feature to measure pixel intensities. An additional rectangular ROI section was used to measure the average background intensity value for each timepoint post axotomy. The adjusted intensity values over a 1 μm region for each timepoint were calculated using the formula ⁠. After obtaining the adjusted intensity values at 5 min and 10 min post axotomy for each animal, the data points were scaled using min-max normalisation such that the highest intensity corresponded to a value of 1 and the lowest intensity corresponded to 0. The normalised intensity values were averaged across animals and the average normalised intensity values were fitted to an exponential decay curve as I(x)=I₀+I₁[exp(−x/λ)], where I(x) represents the intensity at a distance of x micrometres from the cut site, I₀ represents the baseline intensity value, I₁ represents the scaling factor and λ is the decay constant.

Analysis was conducted in R (http://www.R-project.org/) and figures were produced using the package ggplot2 (Wickham, 2009). Breaks in graph axes, wherever required, were introduced using the package ggbreak (Xu et al., 2021). For comparisons between experimental and simulation data, estimation statistics were computed using the ‘dabestr’ package (https://CRAN.R-project.org/package=dabestr) and visualised using Gardner–Altman plots (Gardner and Altman, 1986).

The PLM neuronal process is diametrically 400–500 nm wide, and transverse electron microscopy sections of the neuronal process show that it is densely packed with parallel microtubule bundles that are evenly spaced apart from each other at an average distance of ∼9.7 nm (Cueva et al., 2007). The average diameter of synaptic vesicles in C. elegans neurons is ∼30–40 nm (Nonet, 1999; Choudhary et al., 2017). This means that there is not much space for cargo to traverse between microtubules. Thus, most cargo transport must occur along the periphery of the microtubule bundle (space between the microtubule bundle and the neuronal membrane). This narrow geometry likely does not allow for much interaction between vesicles in three-dimensional space, and most interactions must occur in the two-dimensional space along the microtubule tracks. We modelled an 8 μm region of the PLM neuronal process as a set of ten parallel microtubule tracks, in which each microtubule track was a one-dimensional lattice of 8 nm sites that were equally spaced. Each lattice site was associated with at most a single vesicle (Fig. 2A). To accurately reflect the geometry and spatial constraints of the PLM neuronal process in the model, the peripheral microtubule tracks were placed in a concentric ring and track switching was allowed between neighbouring microtubules on either side (Fig. 2C).

In our simulations, moving particles were vesicle–motor complexes in one of the three motion states observed in vivo – SmA, StA or R (Fig. 1A). We made the reasonable assumption that SmA and R correspond to vesicles transported predominantly by anterograde and retrograde motors, respectively. The step state was taken to correspond to a comparatively mixed distribution of plus- and minus-end-directed motors, exhibiting frequent pauses owing to a tug of war between the motor proteins. Direct transitions between the smooth anterograde and smooth retrograde states are expected to be unlikely (D'Souza et al., 2022 preprint). Thus, we only allowed transitions between the SmA and StA and between R and StA separately, but not the direct SmA to R transition (Fig. 2B). The rates (Table S2) prescribed to such transitions in motion states of vesicles were termed as ‘state interconversion rates’ in our simulations.

We accounted for the interactions between vesicles and their immediate environment in the following way: vesicle motion could be blocked by the presence of another vesicle on the same track. If the centre of a vesicle was associated with a specific site, it was assumed to block three sites in front and three sites at the back from being occupied (Fig. 2A). If a vesicle was blocked from motion, it could change its position or state. Its position could be changed through track-changing events. Alternatively, it could maintain its position but change its state.

We constructed two types of stochastic models: (1) a system without obstacles and (2) a system with obstacles. In the ‘without-obstacles’ model, vesicles move on uninterrupted tracks, and are only blocked by the presence of other vesicles in front of them (Fig. 2D). The frequency of vesicular collisions is dependent on the vesicle density, which is discussed in greater detail below. A more complex version of the ‘with-obstacles’ model incorporates additional barriers to vesicle motion (Fig. 2D). Vesicle motion can be obstructed by the presence of microtubule ends. In the ‘with-obstacle model, six breaks corresponding to 12 microtubule ends (spanning 20 sites each on six different microtubules at different regions) are incorporated every 1000 sites (corresponding to an 8 μm region). The number of microtubule breaks introduced in our model (six distal ends over an 8 μm region for a ten-microtubule system) were very similar to the number of microtubule distal ends observed in serial electron microscopy sections of the PLM neuronal process, as reported in Fig. 5D in Cueva et al. (2007).

It has been observed that sites of stalled mitochondria serve as crowded locations that contribute to the formation of SCs along the PLM neuronal process (Sood et al., 2018). The density of mitochondria along the PLM neuronal process has been previously reported to be five to six mitochondria per 100 μm (Awasthi et al., 2020 preprint; Mondal et al., 2021). A prior study has also reported that the dynamics of these mitochondria are low and occur on much longer timescales than examined for the purpose of this study, with one to two mitochondria moving per hour in an 8–10 μm region (Mondal et al., 2021). In our simulations, mitochondria were modelled as stationary obstacles to transport, with one mitochondrion for every 8–10 μm in the simulations. An examination of the reported mitochondrial distribution showed that this assumption accounted for 65–70% of all observed inter-mitochondrial distances (Awasthi et al., 2020 preprint).

When vesicles reach the end of a microtubule track or encounter large physical obstructions such as mitochondria, they pause as they cannot advance on the same track. Such paused vesicles can eventually either switch to a neighbouring track or reverse their direction through state inter-conversion, with a small probability.

The hopping rates for SmA, StA and R pre-SVs were assigned based on the experimentally derived mean velocity values (Fig. 1E). StA pre-SVs exhibited ∼0.5 times lower mean velocities compared to those of SmA and R pre-SVs (Fig. 1E). R pre-SVs exhibited ∼1.4 times higher mean velocities compared to those of SmA pre-SVs (Fig. 1E). However, for the sake of simplicity, their hopping rates were kept equal to those of SmA in the simulations.

The reversal rates StA→R and R→StA were constrained to the same order of magnitude as that observed in experiments (Fig. 1G). Interconversion between SmA and StA states was only introduced in the model as a tool to modulate the SmA:StA:R proportions so that they could be benchmarked to the ratio observed in experiments (Fig. S1D). These interconversions do not have any meaning in the experimental context, so there are no measures by which to constrain these parameter values. We used trial and error to arrive at these rates by benchmarking simulations to the experimentally observed 50:40:10 ratio of SmA:StA:R (Fig. S1D).

The two types of stochastic models developed were both able to capture key features of axonal transport of cargo in vivo (Fig. S1E). We quantitatively assessed whether the model without obstacles or the model with obstacles was better able to reproduce experiments. In order to test this, we used two measures for comparison between simulations and experiments: (1) net anterograde current, which measures the net anterograde bias (per unit time) in transport of pre-SVs, and (2) fraction of area occupied by SCs over time, which measures the proportion of total sites along the neuronal process that are occupied by SCs of pre-SVs over a given time interval. We found that the model with obstacles produced net current values that were similar to experimental measurements, whereas the model without obstacles generated significantly higher net current values (Fig. S1F). The fraction of area occupied by SCs over time in the model with obstacles was similar to the experimentally derived fraction, whereas the model without obstacles showed a significantly reduced fraction compared to experiments (Fig. S1G). These observations collectively demonstrate that the model with obstacles is a better representation of synaptic vesicle distribution and transport along the PLM neuronal process in vivo.

An examination of the total number of moving pre-SVs in an 8 μm region of the PLM neuronal process showed that their density ranges from 6.25–25 vesicles/μm (Fig. S1H). Accordingly, simulations were carried out with vesicle densities of 6.25, 12.5 and >15 vesicles/μm in a system with 1000 sites spanning an 8 μm region. As vesicle density influences the frequency of cargo-based crowding, and subsequently cargo flow in the simulations, this value was chosen such that the model reproduced experimental measures using an independent assay.

FRAP has previously been used as a measure of active microtubule-based transport in both neuronal and non-neuronal cells (Ciocanel et al., 2017; Dey and Ray, 2018). Recovery occurs when unbleached fluorophores associated with vesicles from neighbouring unbleached regions transit into the bleached portion of the neuronal process. We compared the fluorescence recovery curves obtained from different vesicle densities in the simulations with experiments and used two measures to identify the vesicle density that best reproduced experimental observations: (i) t_1/2 – a measure of the time taken (in seconds) to reach half the intensity of recovery seen at 3 min post bleaching, and (2) the extent of recovery observed at 1, 2 and 3 min post bleaching.

We saw from the representative kymographs in both experiments and simulations (Fig. S2A) that the recovery in fluorescence after photobleaching comes largely from active vesicular transport, with negligible contributions from diffusive processes. We compared the fluorescence recovery curves from simulations with obstacles, run at three different vesicle densities with experimental measurements. We found that simulations run at the highest vesicle density were best able to reproduce the experimental bleach recovery curve (Fig. S2B). Quantitative comparison of (Fig. S2C) and extent of recovery at 1 min (Fig. S2D), 2 min (Fig. S2E) and 3 min (Fig. S2F) post bleaching between simulations and experiments also showed that simulations run at the highest vesicle density were the most similar to experimental results.

We further observed that the extent of recovery in simulations run at lower vesicle densities was higher than the extent of recovery at high vesicle densities (Fig. S2B,D–F). Our simulations showed that higher vesicle densities corresponded to a reduced extent of fluorescence recovery owing to lower net anterograde cargo transport (Fig. S2B,D–F). This addresses previously proposed open questions in the field on the effect of vesicle density on cargo transport kinetics (Lai et al., 2018). We saw that in addition to static physical obstacles, higher vesicle densities and vesicle–vesicle interactions in the system might contribute to SC formation. These stalled vesicles can present a barrier to the transport of unbleached vesicles from neighbouring sites. This is consistent with predictions from previous theoretical studies on the formation of cargo pileups along microtubules (Lai et al., 2018).

The probability for moving cargo to switch microtubule tracks relies on a diverse set of parameters, such as (1) the type of motor, (2) the number of motors engaged on the cargo surface, (3) the angle of intersecting microtubules, (4) the size of the cargo being transported and (5) axial separation between microtubules (Ross et al., 2008; Bálint et al., 2013; Bergman et al., 2018; Lee and Higuchi, 2019). Most of these observations are derived from in vitro studies, and the ability of moving synaptic vesicles to switch tracks along neuronal processes in vivo is poorly understood. As it is difficult to predict the extent to which the SmA, StA and R populations of vesicles differ with respect to their track-switching rates in vivo, for the purposes of this simulation, we assumed the track-switching rates to be equal for the three populations of vesicles.

The rate at which moving pre-SVs switch microtubule tracks is expected to be much lower than the rate at which pre-SVs move along the microtubule track. The vesicle hopping rates used in the simulations ranged from 40–100 hops/s. Thus, the system was simulated with track-switching rates ranging from 10⁻⁴ to 10⁻¹/s (2×10^–1, 2×10^–2, 2×10^–3 and 6×10^–4/s). Simulations run at different track-switching rates were then benchmarked to experiments using multiple measures, such as (1) FRAP, (2) net current and (3) area occupied by SCs over time (Fig. S3). Track-switching rates were benchmarked to experiments using the system with obstacles at high vesicle density, as it had already been established to better reproduce experimental measurements.

We saw that simulations with a track-switching rate of 0.02/s were best able to reproduce the experimental recovery curve (Fig. S3A), t_1/2 values (Fig. S3B), and the extent of fluorescence recovery at 1 min (Fig. S3C), 2 min (Fig. S3D) and 3 min (Fig. S3E) post bleaching. Furthermore, the track-switching rate of 0.02/s was also able to reproduce the experimental net current values (Fig. S3F) and the fraction of area occupied by SCs over time (Fig. S3G). We also observed that an increase in track-switching rates caused an increase in net anterograde current (Fig. S3F) and a decrease in the fraction of area occupied by SCs over time (Fig. S3G). In the absence of track switching (rate=0/s), the net current values were very low (Fig. S3F), whereas the fraction of area occupied by SCs over time was very high (Fig. S3G). This corroborates with prior theoretical models that show that higher rates of track switching are associated with reduced cargo pileups (Lai et al., 2018). Our simulations suggest that track switching plays an important role in maintaining net cargo flow and limiting SC formation. Owing to technical limitations, it was not possible for us to observe switching of vesicles between microtubule tracks along the PLM neuronal process in vivo. Thus, we did not experimentally investigate the role of track switching in maintaining cargo transport along the neuronal process.

As all vesicles were tracked in the simulations, it was easy to generate kymograph tracks from them. In order to make these tracks visually similar to experimental kymographs, we assigned different greyscale values and line thicknesses to vesicles depending on their motion. With this in mind, we chose the following algorithm for generating our kymographs. We plotted all vesicles, with the structure of the kymograph representing the relative proportion of stationary versus non-stationary vesicles. In our kymographs, there were two parameters for the tracks: the greyscale value and the thickness. Of all the vesicles, we identified those that were part of long-lived clusters, representing them by thicker (darker) lines in the kymograph. The greyscale was assigned in proportion to the time spent by such vesicles in SCs. To set the thickness and greyscale values, we computed the maximum time each vesicle spent at a site, assigning a line thickness proportional to the time spent in a stationary state. For the greyscale value, a lighter shade of grey was used for vesicles that showed substantial movement, whereas a darker shade of grey (RGB value: 0.1, 0.1, 0.1) was used for vesicles that were largely stationary. These choices yielded simulation kymographs that visually resembled those obtained from experiments.

We took simulation snapshots of the axon at time intervals that corresponded to the frame rate of experimental movies. We collapsed sites with the same index across microtubules, as vesicles on different tracks corresponding to the same indexed site could not be resolved separately in experiments. A site was considered to be occupied, contributing to the fluorescence signal, if the corresponding site in any of the microtubules was occupied. When a site was occupied across two adjacent snapshots, we considered it to be an SC and the duration of the respective SC was incremented by an amount equivalent to the time interval of the snapshot. The duration assigned to an SC was computed as the duration for which the SC was observed within the time window of the simulation. As in the experiments, long-lived SCs were defined as clusters that were observed to survive for longer than 45 s. All vesicles in simulations could be tracked throughout their lifetimes.

Similar to analysis of experimental data, an anterogradely moving vesicle that changed its direction to retrograde along its trajectory was counted as an anterograde reversal, whereas the retrograde-to-anterograde change was counted as retrograde reversal.

The current (equivalently the flux) was obtained by computing the time-averaged currents of different states of vesicles in the following way: choosing a link connecting two neighbouring sites, we computed the number of anterograde vesicles of both states (smooth and step) moving across it, as well as the number of retrograde vesicles, over an interval of time. Subtracting the numbers of retrograde vesicles from anterograde vesicles and dividing by the elapsed time gave us the current, which could then be computed over a large number of configurations to generate averages in steady state.

Stationary cargo site occupancy analysis for the simulation data was performed in a manner similar to that in the experiment. The rectangular ROI for each SC in the kymograph was calculated as the product of the lattice width and the time length occupied by the SC. The rectangular ROI data for all SCs in the system was then represented as the fraction of the total area of the kymograph.

FRAP experiments were simulated on a 64 μm lattice with 400, 800 and 1000 vesicles. Once the system settled into a steady state, vesicles in an 8 μm section of the lattice were tagged as ‘bleached’ – in effect, hiding them from that point on in the simulation. The recovery of unbleached vesicles in the bleached region was then followed over time. Although the bleached vesicles continued to hop, change track, change type and interact as before, they were excluded from subsequent analysis of bleach recovery and generation of kymographs.

We simulated axotomy by blocking a group of 100 sites (0.8 μm) on all contiguous microtubules. Vesicles could thus move up to but not across the axotomised region. This led to an accumulation of vesicles over time and subsequent saturation of the density at and close to the two extremal sites of the axotomised region. Exponential fits of vesicle densities in the vicinity of the axotomized region were determined and compared to the experiments.

Code for counting and binning of manually annotated experimental images are available at https://github.com/amruta2612/Fiji-macros. Codes for simulation are available at https://github.com/reshmamaiya/transportCrowdedSystem.

G.I.M. is grateful to the Department of Biological Sciences, Tata Institute of Fundamental Research, for an Adjunct Professorship. We also thank the Institute of Mathematical Sciences, Chennai, for the use of its high-performance computing facilities. We also thank Badal Singh Chauhan for generating and integrating the transgenic strain used for our axotomy experiments.