# Unraveling Single-Particle Trajectories Confined in Tubular Networks

Yunhao Sun∗ Department of Physics, University of California, San Diego, San
Diego, California 92093  Zexi Yu∗ Department of Physics, University of
California, San Diego, San Diego, California 92093  Christopher J. Obara
Janelia Research Campus, Howard Hughes Medical Institute, Ashburn, Virginia
20147  Keshav Mittal  Department of Physics, University of California, San
Diego, San Diego, California 92093  Jennifer Lippincott-Schwarz  Janelia
Research Campus, Howard Hughes Medical Institute, Ashburn, Virginia 20147
Elena F. Koslover<EMAIL_ADDRESS>Department of Physics, University of
California, San Diego, San Diego, California 92093

(August 27, 2024)

###### Abstract

The analysis of single particle trajectories plays an important role in
elucidating dynamics within complex environments such as those found in living
cells. However, the characterization of intracellular particle motion is often
confounded by confinement of the particles within non-trivial subcellular
geometries. Here, we focus specifically on the case of particles undergoing
Brownian motion within a tubular network, as found in some cellular
organelles. An unraveling algorithm is developed to uncouple particle motion
from the confining network structure, allowing for an accurate extraction of
the diffusion coefficient, as well as differentiating between Brownian and
fractional Langevin motion. We validate the algorithm with simulated
trajectories and then highlight its application to an example system:
analyzing the motion of membrane proteins confined in the tubules of the
peripheral endoplasmic reticulum in mammalian cells. We show that these
proteins undergo diffusive motion with a well-characterized diffusivity. Our
algorithm provides a generally applicable approach for disentangling geometric
morphology and particle dynamics in networked architectures.

##  I Introduction

Particle tracking experiments have a long history in the field of soft matter
physics, where they are used to analyze the material properties of complex
fluids and the dynamic behavior of active media (reviewed in Wirtz (2009);
Squires and Mason (2010)). In recent years, tracking of particles inside
living cells has been extensively employed to elucidate the dynamics of
cellular components ranging from single molecules Elf  _et al._ (2007); Parry
_et al._ (2014); Holcman  _et al._ (2018) to vesicular organelles Guo  _et
al._ (2014); Gal  _et al._ (2013); Koslover  _et al._ (2016); Lin  _et al._
(2016). Quantification of in vivo particle trajectories can be leveraged to
identify state transitions and organelle interactions Lin  _et al._ (2016);
Chen  _et al._ (2015); Bálint  _et al._ (2013), to explore the rheology of
intracellular fluidsWeihs  _et al._ (2006); Gal  _et al._ (2013), and to
establish the underlying physical forces that drive particle motion Weber  _et
al._ (2012a); Guo  _et al._ (2014); Fodor  _et al._ (2015).

A classic analysis approach computes the mean square displacement (MSD) of
particle trajectories. In a purely thermal homogeneous system, an MSD that
scales linearly with time indicates diffusive motion through a viscous fluid,
with the prefactor establishing the particle diffusivity. An MSD that scales
subdiffusively ($\sim t^{\alpha}$, $\alpha<1$) can instead indicate a
viscoelastic medium with characteristic power-law scaling $\alpha$
Kollmannsberger and Fabry (2011). This analysis has been employed in a number
of cellular systems, establishing diffusive behavior in the case of some
protein-sized particles Wang  _et al._ (2006); Li  _et al._ (2015); Etoc  _et
al._ (2018), and subdiffusive dynamics for vesicles and similar-sized
exogenous probes in the cytoplasm Kollmannsberger and Fabry (2011); Yamada
_et al._ (2000); Hoffman  _et al._ (2006); Weihs  _et al._ (2006). An
alternate recent approach focuses instead on the velocity auto-correlation
function, with negative correlations that decay as a power-law in time taken
to be a sign of viscoelastic rheology Weber  _et al._ (2012b); Sabri  _et al._
(2020); Lucas  _et al._ (2014). Still other studies analyze the distribution
of individual step sizes over short time-intervals Wang  _et al._ (2012); Gal
_et al._ (2013); Lampo  _et al._ (2017).

All of these approaches are confounded by a number of complications in the
intracellular environment Mogre  _et al._ (2020). First and foremost, the
prevalence of active forces with many different correlation timescales implies
that individual particle trajectories cannot be directly related to the
rheology of the medium Guo  _et al._ (2014). Furthermore, heterogeneity in the
intracellular medium indicates the existence of spatially varying
diffusivities Duits  _et al._ (2009); Li  _et al._ (2015) and is thought to be
responsible for the broad non-Gaussian distribution of step sizes observed for
many particles Wang  _et al._ (2012); Lampo  _et al._ (2017); Chechkin  _et
al._ (2017). Nevertheless, extracting the effectively diffusive or
subdiffusive behavior of intracellular particles, albeit over a finite range
of timescales, can give insight into functionally important consequences such
as particle search times, kinetics, or interaction frequencies.

In certain cases it is possible to disentangle the underlying particle
dynamics and the effects of forces, flows, or confinement from the measured
particle trajectories. For example, when particles are driven by slowly-
varying fluid flows overlaid on top of diffusive behavior, the raw MSD appears
super-diffusive. However, the diffusive behavior can be extracted by
subtracting out a smoothed trajectory and analyzing the resulting MSD curves
with an appropriate rescaling Koslover  _et al._ (2016). In other cases,
separate measurement of the confounding factors is necessary. For instance,
active microrheology measurements of cytoplasmic material properties enable
‘passive’ particle trajectories to then be used for extraction of the spectrum
of forces driving particle motion Guo  _et al._ (2014). Mapping out the
underlying surface on which a particle is confined and measuring distances
along that surface has been shown to enable accurate characterization of
diffusing particle dynamics Adler  _et al._ (2019).

Confinement effects are a common source of complication in analyzing the
trajectories of intracellular particles. Confinement in a finite region can
result in apparently subdiffusive MSD curves, on time scales up to an order of
magnitude shorter than the typical time to traverse the confining region
Saxton (1995). Similarly, velocity autocorrelation functions for particles in
confinement exhibit negative peaks that increase when the velocity is measured
across longer time intervals Weber  _et al._ (2012b). In this manuscript, we
focus on how to decouple confinement from the underlying particle dynamics in
the case of particles moving along tubular networks.

Eukaryotic cells contain a number of structures that result in confinement of
particles at different length scales. The extent of the cell itself provides
an upper limit for confinement. Furthermore, many particles are embedded in
the membrane or interior of organelles that offer confinement on a subcellular
scale. One example of interest is the reticulated mitochondrial network formed
in many cell types. These networks exhibit varying connectivity that can be
tuned by genetic perturbation of the proteins responsible for mitochondrial
fusion and fission Viana  _et al._ (2020). Another organelle with a networked
morphology is the endoplasmic reticulum (ER), whose functional roles include
lipid distribution, calcium buffering, and protein processing, sorting, and
quality control Schwarz and Blower (2016); Perkins and Allan (2021). The ER
forms an interconnected system of tubules and stacked sheets, with a
continuous lumen and membrane, that spans throughout the cell Westrate  _et
al._ (2015); Friedman and Voeltz (2011). Away from the cell nucleus, the ER
can be approximated as a tubular network with largely three-way junctions,
through which secretory proteins must move in order to encounter the point-
like ER exit sites Chen  _et al._ (2013); Griffing (2010); Barlowe and
Helenius (2016).

Mathematical models of transport within the ER and mitochondrial networks
indicate that the network morphology has the potential to alter search rates
and kinetics for proteins to find each other and specific targets within the
network Brown  _et al._ (2020); Viana  _et al._ (2020); Scott  _et al._
(2021). However, examining the interplay between network architecture and
kinetics requires an assumed physical model for particle motion (eg: diffusive
Brown  _et al._ (2020); Viana  _et al._ (2020), subdiffusive Stadler  _et al._
(2018), or locally persistent Holcman  _et al._ (2018); Dora and Holcman
(2020)) and an accurate parameterization of that model. Thus, precise
empirical quantification of particle dynamics within reticulated networks is
of growing interest for elucidating the structure-function relationship of
these critical cellular organelles.

Most prior studies on the movement of ER membrane and luminal proteins have
focused on bulk measurements of protein spread, including FRAP analysis Dayel
_et al._ (1999); Siggia  _et al._ (2000); Nehls  _et al._ (2000) and the
spatiotemporal quantification of spreading from a locally photoactivated
region Holcman  _et al._ (2018); Konno  _et al._ (2021). However, it has
recently become possible to track the movement of individual proteins through
the ER tubules, enabling new observations of dynamics within the ER. For
example, single particle trajectories of ER luminal proteins have been found
to exhibit unexpectedly fast motion along individual tubes, followed by
trapping at junctions – an effect that has been attributed to putative luminal
flows over short timescales Holcman  _et al._ (2018). ER membrane protein
trajectories, however, do not appear subject to these rapid motions.
Nevertheless, the analysis of such trajectories is complicated by their
confinement within the peripheral ER network structure.

In this work we present a novel method for analysis of diffusive particle
trajectories that are confined within a spatial network, consisting of tubules
connected by narrow junctions. Namely, we describe how such trajectories can
be ‘unraveled’ to establish a properly sampled underlying trajectory that
describes the motion of an identical particle on an infinite line. This
unraveling process allows the accurate extraction of a diffusion coefficient
for the particle motion by removing the confounding effects of the network
structure itself. The procedure is generally applicable to analysis of any
locally one-dimensional trajectories where the confining network architecture
can be separately imaged. Thus, for example, it can be applied to particles
moving along a single finite-length tube, or along a spine-studded morphology
such as that found in neuronal dendrites Santamaria  _et al._ (2006); Sartori
_et al._ (2020).

We validate the method using simulations of particles diffusing on a network,
and demonstrate that the technique can differentiate between particles
undergoing Brownian versus fractional Langevin motion (a common model for
subdiffusion due to viscoelastic rheology Mogre  _et al._ (2020); Deng and
Barkai (2009)). The technique is then applied to the analysis of membrane
protein trajectories in the ER, showing that these trajectories are indeed
consistent with diffusive motion and providing an estimate of their
diffusivity.

Overall, the work presented here enables quantitative analysis of particle
trajectories confined on network structures, allowing for mathematically
accurate decoupling of confinement effects from the underlying particle
dynamics.

Figure 1: Simulation of diffusive particle trajectories on a network. (a)
Schematic of diffusion simulation and unraveling algorithm. Top: example
single-step trajectories starting at $x_{0}$ on an infinite line; Bottom:
sample corresponding trajectories along the network. Trajectories that pass
the node can land on either the same or a different edge. (b) Simulations of
particle trajectories using the discrete time propagation algorithm (dots)
reproduce exact analytic distributions (lines) on a triskelion network.
Densities are shown for a particle starting at the center node, propagating
for a total time $t=0.8\ell^{2}/D$. (c) Diffusive particle mean squared
displacement depends on the structure of the confining network. Dimensionless
MSD is shown for particles on a triskelion geometry (blue), a honeycomb
network (magenta), a honecomb with 30% of its edges removed (red), and a
section of a peripheral ER network extracted from a COS7 cell (green). Length
units are scaled by the average network edge length ($\ell$) and time units by
$\ell^{2}/D$.

##  II Results

###  II.1 Simulating diffusive particles on a network

We begin by simulating the behavior of diffusive particles on a tubular
network. Individual narrow tubules are treated as effectively one-dimensional
segments, joined together at point-like junctions. The network structure is
thus idealized as a set of nodes connected by edges with well-defined length:
$\ell_{m}$ for the $m^{\text{th}}$ edge. The edges are not necessarily
straight, so that $\ell_{m}$ can be longer than the spatial distance between
two adjacent nodes.

Several algorithms are possible to simulate the trajectory of a diffusive
particle on such a network. In recent work, we demonstrated an exact kinetic
Monte Carlo algorithm that leverages analytical propagator functions obtained
by solving the diffusion equation in the neighborhood of each node Scott  _et
al._ (2021). By sampling transition times from these functions, the particle
can be propagated between adjacent nodes as a continuous-time random walk,
with no discretization artifacts. Here, we take an alternate approach, which
starts from a diffusive particle trajectory $z(t)$ on an infinite line and
maps it onto a corresponding trajectory $x(t)$ on the network. In a subsequent
section, we demonstrate how this mapping can be reversed to extract the
underlying diffusive dynamics of a particle whose motion is observed on a
network of known structure.

Consider a discretized trajectory $\left\\{z_{i}\right\\}$ of a Brownian
particle with diffusivity $D$, evaluated at discrete times $t_{i}=i\Delta t$.
The increments of such a trajectory are distributed according to a normal
distribution: $\Delta z\sim\mathcal{N}(0,\sqrt{2D\Delta t})$. To map this
trajectory onto an equivalent particle diffusing over the network structure,
we use these increments $\Delta z$ for propagation on the network. In
particular, we assume that each time the diffusing particle passes a network
node, it selects uniformly at random the next edge from the ones attached to
that node. Our approach assumes that the discrete time step is short enough
that on any single step the particle can pass only the network node that is
closest to its starting position (in terms of distance along the edge).
Specifically, we assume that $|\Delta z|<\ell_{m}/2$ for the current edge $m$
on which the particle is located.

To map an unconstrained trajectory step $\Delta z$ onto the network (Fig. 1a),
consider a particle starting at position $x_{0}$ along edge $m$, and assume
that $x_{0}<\ell_{m}/2$, so that the particle is closer to the first node of
the edge. We want to find the next position $x_{1}$ for the particle, given
its sampled trajectory on the infinite line. If $\Delta z<-x_{0}$ then the
particle trajectory must have passed the node during this time step. We then
place the particle at position $x_{1}=|x_{0}+\Delta z|$ on an edge selected
uniformly at random from among all the edges attached to that node (including
the original edge). If $\Delta z>-x_{0}$ then the particle may also have
passed the node during that time step. The probability that it has done so can
be computed as

$\begin{split}p_{\text{pass}}=\exp\left[-\frac{(x_{0}+\Delta z)x_{0}}{D\Delta
t}\right],\end{split}$ (1)

(see Appendix A.1 for derivation). If it passed the node, then it can end the
time step on any of the adjacent edges with equal probability; if it did not
pass the node then it must still be on its original edge. The subsequent
position of the particle is then given by $|x_{0}+\Delta z|$ with probability
$1-p_{\text{pass}}+p_{\text{pass}}/d$ of landing on its original edge and
probability $p_{\text{pass}}/d$ of landing on each of the other adjacent
edges, where $d$ is the degree of the node. The propagation of the particle is
then repeated from its new position. If the particle originates close to the
other end of the edge ($x_{0}>\ell_{m}/2$) then an analogous calculation is
performed by considering the nearby node on the far end of edge $m$.

It should be noted that this algorithm also handles reflections off of a dead-
end node. Increments that would lead an unconstrained particle to pass such a
node simply result in the particle being placed back upon the single edge
leading to the dead-end. Such explicit reflections are typically implemented
in Brownian and fractional Brownian dynamics simulations with hard-wall
boundariesJeon and Metzler (2010).

The proposed algorithm accurately represents the behavior of a diffusive
particle over the time step $\Delta t$ so long as we can neglect any
trajectories that pass a node other than the most nearby one over the course
of a single time step. Therefore, unlike the exact algorithm described in
previous work Scott  _et al._ (2021) it requires the constraint
$\sqrt{2D\Delta t}\ll\ell_{m}$. In Fig. 1b, we demonstrate that for small
$\Delta t$ the agent-based simulations described here accurately reproduce
analytically computed particle distributions in a simple triskelion network
(see Appendix A.2).

Figure 2: Unraveling simulated trajectories to estimate diffusivity. (a)
Dimensionless MSD of particle trajectories on a triskelion (blue), a decimated
honeycomb network (red), and a segment of peripheral ER (green). All
simulations were run with $D=1$. Original MSD (from Fig. 1c) is shown with
dashed lines, MSD after unraveling (with $D=1$) shown with solid lines. (b)
MSD of trajectories on decimated honeycomb network, unraveled using different
values of the diffusivity $D$ (thick lines), compared to the expected $2Dt$
(thin lines). Dashed black line indicates values for the correct diffusivity.
All units are non-dimensionalized by the edge length and the simulated
diffusivity. (c) Residual function $G(D)$ plotted for the different
diffusivity values used in unraveling. Dashed line shows the correct $D=1$
used in the simulation, which matches well to the minimum of the residual
function.

###  II.2 Network confinement modulates particle MSD

A traditional metric for analyzing particle trajectories is to look at the
mean squared displacement (MSD), defined by

$\begin{split}\text{MSD}(t)=\left<\left|\vec{r}(t)-\vec{r}(0)\right|^{2}\right>,\end{split}$
(2)

where $\vec{r}(t)$ gives the spatial position of the particle at time $t$, and
the average is taken both over time and over an ensemble of many particles.
Particles that diffuse unimpeded on an infinite linear domain exhibit the
relationship $\text{MSD}=2Dt$. Identical behavior is observed for particles
diffusing on a network comprising an infinite lattice of tubules. The effect
of confinement in a finite domain is a well-known complication of MSD analysis
Saxton (1995), with the MSD beginning to flatten on time scales comparable to
$\sqrt{R^{2}/2D}$ for a domain of radius $R$ (Fig. 1c).

When the connectivity of the network structure is reduced below that of the
well-connected lattice, the time dependence of the mean squared displacement
flattens still further. Scaling analysis of percolation systems indicates that
particles moving on the largest connected cluster of a network above the
percolation phase transition should maintain effectively diffusive behavior,
with an effective diffusivity that decreases according to
$D_{\text{eff}}\sim(p-p_{c})^{\mu}$, where $p$ is the neighbor connection
probability, $p_{c}$ is its critical value, and the scaling exponent is
$\mu\approx 1.3$ for planar systems Stauffer and Aharony (2018). However, for
a finite network and over a finite observation time, the MSD on a lattice with
missing connectivity can exhibit behavior that looks distinctly non-diffusive.
In Fig. 1c (red curve), we show the MSD of particles on a honeycomb network
with $30\%$ of its edges removed in such a way as to maintain a single
connected component. The MSD appears to show a subdiffusive (sublinear)
scaling with time, as missing connectivity forces particles to take much
longer paths to get from one part of the network to another Brown  _et al._
(2020).

Cellular organelle networks, such as the peripheral endoplasmic reticulum and
mitochondrial network structures, also show reduced connectivity, with dead-
end nodes and heterogeneous pore sizes between connected edges. Simulating
diffusive particle trajectories on such networks indicates that the MSD is
also expected to scale sublinearly with time (Fig. 1c, green curve).

These simulations highlight the inherent difficulty of quantifying particle
diffusivity on a network by analyzing the mean squared displacement of
observed trajectories. Namely, the MSD convolves together two separate
effects: the underlying dynamics of the moving particles and the morphology of
the network within which they are confined. We thus proceed to develop a
method for separating out these two effects and extracting particle
diffusivity from their behavior on a known network structure.

###  II.3 Unraveling diffusive trajectories on networks

To unravel the particle dynamics from the network architecture, we proceed by
inverting the algorithm used to map the discrete trajectory of an unconfined
particle ($\left\\{z_{i}\right\\}$) onto a trajectory over the network (Fig.
1a). Given a set of discrete 2D (or 3D) observations of the network-confined
particle ($\left\\{\vec{r}_{i}\right\\}$), we seek to find the trajectory that
particle would have taken if allowed to move over an infinite line. To begin,
we re-express the spatial trajectory in terms of network coordinates
$\left\\{m_{i},x_{i}\right\\}$ that list the edge $m_{i}$ and position along
that edge $x_{i}\in(0,\ell_{i})$ of the particle at time point $t_{i}$. We
then consider the start and end positions of the particle at each time-step,
assuming that there is an underlying free one-dimensional trajectory such that
$z_{i},z_{i+1}$ were mapped to network positions
$(m_{i},x_{i}),(m_{i+1},x_{i+1})$ using the algorithm described in the
previous section.

Figure 3: Diffusivity estimates for different timesteps. Simulated
trajectories on a decimated honeycomb network (red) and a triskelion network
(blue) are unraveled to obtain an estimated $D$ value. Each analysis involves
100 independent trajectories of 100 time steps each. Error bars give standard
deviation from $20$ replicate analyses. All values are non-dimensionalized by
the edge length $\ell$ and the simulated diffusivity $D$.

For sufficiently small time-steps, we can again assume that the particle could
only have passed over the node closest to $x_{i}$. Without loss of generality,
we take this node to be the one located at position $0$ on edge $m_{i}$. For
simplicity of notation, we shift the unconstrained trajectory such that
$z_{i}=x_{i}$. The next observed position of the unconstrained particle must
then be $z_{i+1}=\pm x_{i+1}$. Bayes’ rule Ross (2014) allows calculation of
the probability associated with each choice, conditional on the observed
particle positions. The probability that the unconstrained trajectory ended at
$z_{i+1}=x_{i+1}$ after the time-step, given that the particle ended at
$(m_{i+1},x_{i+1})$ on the network is:

$\begin{split}p_{+}=\mathcal{P}(z_{i+1}=x_{i+1}|m_{i+1},x_{i+1})=\\\
\frac{\mathcal{P}(m_{i+1},x_{i+1}|z_{i+1}=x_{i+1})\mathcal{P}(z_{i+1}=x_{i+1})}{\mathcal{P}(m_{i+1},x_{i+1})}\end{split}$
(3)

The free particle increments are sampled from a normal distribution, so that

$\begin{split}\mathcal{P}(z_{i+1}=\pm x_{i+1})=\frac{1}{\sqrt{4\pi D\Delta
t}}e^{-\frac{(x_{i}\mp x_{i+1})^{2}}{4Dt}}.\end{split}$ (4)

If the particle moved on to a different edge ($m_{i+1}\neq m_{i}$), then the
unconstrained trajectory must have passed zero and the conditional
probabilities of the observed end-point on the network are:

$\displaystyle\mathcal{P}(m_{i+1},x_{i+1}|z_{i+1}=x_{i+1})$ $\displaystyle=$
$\displaystyle\frac{1}{d}e^{-\frac{x_{i}x_{i+1}}{D\Delta t}},$ (5)
$\displaystyle\mathcal{P}(m_{i+1},x_{i+1}|z_{i+1}=-x_{i+1})$ $\displaystyle=$
$\displaystyle\frac{1}{d}.$ (6)

The first line gives the probability that the unconstrained trajectory passed
$0$ (as in Eq. 1) before returning to the positive side. The second line
simply indicates that if the unconstrained trajectory ends on the negative
side, it must have passed zero and is thus equally likely to be mapped to any
of the adjacent edges.

For a particle that ends on the same edge as its starting point
($m_{i+1}=m_{i}$), the conditional probability given $z_{i+1}=x_{i+1}$ is
altered to

$\begin{split}\mathcal{P}(m_{i},x_{i+1}|z_{i+1}=x_{i+1})&=1+\left(\frac{1}{d}-1\right)e^{-\frac{x_{i}x_{i+1}}{D\Delta
t}},\end{split}$ (7)

which includes an additional term for unconstrained particles that never pass
the origin.

Plugging into Bayes’ rule (Eq. 3) results in the following probabilities for
the end-point of the unconstrained trajectory:

$p_{+}=\begin{cases}\frac{1}{2},&\quad m_{i+1}\neq m_{i}\\\
\frac{1}{2-d+d\exp(x_{i}x_{i+1}/Dt)},&\quad m_{i+1}=m_{i}.\end{cases}$ (8)

The conditional probability $p_{-}$ that $z_{i+1}=-x_{i+1}$ given the observed
values of $m_{i+1},x_{i+1}$ can be computed as $p_{-}=1-p_{+}$.

For the network trajectory $\left\\{m_{i},x_{i}\right\\}$, we obtain an
‘unraveled’ trajectory $\left\\{z_{i}\right\\}$, by sampling the next position
at each step according to Eq. 8, considering only the node closest to $x_{i}$
for potential passage by the unconstrained particle. The mean squared
displacement and velocity autocorrelation function for the unraveled
trajectory can then be analyzed in the usual manner. As shown in Fig. 2a, the
unraveled trajectory mean squared displacement ($\text{MSD}_{\text{ur}}$) for
simulated particles diffusing on a network regains its linear time scaling,
with the correct prefactor: $\text{MSD}_{\text{ur}}=2Dt$.

Computing the unraveled trajectories requires knowledge of the particle
diffusivity $D$ in order to correctly sample the step direction from Eq. 8.
Since one of the primary goals of the unraveling procedure is to gain an
estimate of particle diffusivity, $D$ cannot be assumed a priori. Instead, the
unraveling process is performed with a range of different diffusivities to
find a self-consistent value. Values of $D$ that are too large or too small
result in an unraveled MSD that deviates from the expected behavior (Fig. 2b).
We estimate the underlying particle diffusivity by finding the value of $D$
that minimizes the difference between the computed MSD${}_{\text{ur}}(t;D)$
and $2Dt$. Specifically, we compute normalized squared residuals on log-log
axes to assess the goodness of fit:

$\begin{split}G(D)=\frac{\sum_{i}\log^{2}\left[\text{MSD}_{\text{ur}}(t_{i}|D)/(2Dt_{i})\right]}{\sum_{i}\log^{2}\left[\text{MSD}_{\text{ur}}(t_{i}|D)/\left<\text{MSD}_{\text{ur}}\right>_{t}(D)\right]}\end{split}$
(9)

where the residuals are evaluated at logarithmically spaced time points
$t_{i}$ and $\left<\text{MSD}_{\text{ur}}\right>_{t}$ refers to the average
value over all time-points. The final diffusivity estimate is the value of $D$
that minimizes $G(D)$ (Fig. 2c). We note that the choice of goodness-of-fit
function is equivalent to $1-R^{2}$ for the standard $R^{2}$ metric used to
assess the accuracy of a linear regression, following a logarithmic transform.

Figure 3 shows how the estimated diffusivity varies with the time-step of the
trajectory. Our approach is valid only if the time-step is small enough that
individual consecutive snapshots are very unlikely to involve the particle
trajectory passing any network node other than the one nearest to the
particle. During the analysis, any steps that involve the particle jumping
onto a non-adjacent edge result in breaking the trajectory up into separate
segments. For larger step-sizes (roughly, $\Delta t\gtrsim 0.1\ell^{2}/D$),
such large steps become increasingly common. Breaking the trajectory then
results in a systematic underestimate of the diffusivity since large steps are
removed from the analysis. Thus, in order to accurately assess particle
dynamics on a network, the particle positions must be visualized rapidly
enough that the particle passes no more than one node on each step.

###  II.4 Fractional Langevin motion and velocity autocorrelations

A common application of MSD analysis for single particle trajectories is to
elucidate the relevant dynamic model that best describes their behaviorGal
_et al._ (2013). In particular, particles that exhibit an MSD scaling linearly
with time are generally assumed to be moving in an effectively diffusive
manner. By contrast, subdiffusive scaling ($\text{MSD}\sim t^{\alpha}$ for
$0<\alpha<1$) can be explained by a number of underlying physical models,
including continuous time random walks, confinement, and fractional Langevin
motion (fLM) Weber  _et al._ (2012b); Metzler and Klafter (2000). The latter
model has been shown to apply to a variety of cytoplasmic and intranuclear
particles Lucas  _et al._ (2014); Lampo  _et al._ (2017); Weber  _et al._
(2010a); Sabri  _et al._ (2020). Fractional Langevin motion is expected for
particles subject to thermal fluctuations in a viscoelastic medium, driven by
fractional Brownian forces with a power-law correlation over time that scales
as $t^{-\alpha}$Jeon and Metzler (2010); Bressloff and Newby (2013). We set
out to establish whether unraveled trajectories can be leveraged to
distinguish if particles are undergoing classical Brownian motion or fLM,
decoupling the underlying particle dynamics from the confinement effect of the
network.

To simulate the dynamics of a particle undergoing fLM on a network, we first
generate fractional Brownian forces in two dimensions (2D) using previously
established methods Deng and Barkai (2009); Kroese and Botev (2013). Particle
step increments in 2D are then computed by convolving the past history of
steps with the appropriate power-law memory kernel, followed by zeroing out of
the step component perpendicular to the current edge axis (details in Appendix
A.4). The overall approach is analogous to that used for simulating the
fractional Langevin equation in a confined geometry Jeon and Metzler (2010),
in the limit of no inertia and narrow edge confinement.

Figure 4: Unraveling of simulated fractional Langevin motion trajectories on a
decimated honeycomb network. (a) MSD of raw trajectories, with dashed line
showing expected scaling. (b) MSD of trajectories unraveled with different
input values of $D_{\alpha}$, as in Fig. 2b. Inset shows residual function
$G(D_{\alpha})$ for the unraveling attempts with different $D$. Simulations
used scaling exponent $\alpha=0.6$, with prefactor set such that individual
step sizes are comparable to the Brownian simulations in Fig. 2, with
$2D_{\alpha}\Delta t^{\alpha}=(0.1\ell)^{2}$. Length and time units are non-
dimensionalized by $\ell$ and $\tau=(\ell^{2}/D_{\alpha})^{1/\alpha}$,
respectively.

Figure 5: Estimating scaling exponent $\alpha$ from velocity autocorrelation
function of unraveled trajectories. Normalized velocity autocorrelations are
shown for unraveled trajectories of simulated particles undergoing (a)
Brownian motion and (b) fractional Langevin motion with $\alpha=0.6$. Red box
indicates first negative peak, used to extract estimated scaling exponent
$\alpha_{\text{est}}$. The individual step size is set such that $\left<\Delta
x^{2}\right>=2D_{\alpha}\Delta t^{\alpha}=(0.1\ell)^{2}$. (c) Estimated
scaling exponent for simulated fLM particles with different true values of
$\alpha$. Estimations are shown based on the velocity autocorrelation of
either raw 2D trajectories (solid) or unraveled 1D trajectories (dashed).
Simulations were carried out with two different root mean squared step sizes:
$10\%$ (blue) and $20\%$ (red) of the network edge length. All simulations
included $100$ particles tracked for $100$ time steps. Error bars are standard
deviations from $50$ independent replicates. Cyan bar indicates the range of
$\alpha_{\text{est}}$ that should be taken to imply purely diffusive motion,
obtained as average $\pm$ standard deviation for Brownian particle
simulations. Right: zoomed-in section shows that network confinement makes
Brownian particles ($\alpha=1$) appear to have a lowered estimated scaling
($\alpha_{\text{est}}<1$), whereas unraveled trajectories show the correct
Brownian exponent $\alpha_{\text{est}}=1$.

As expected, the simulated trajectories exhibit an MSD that scales as a power-
law with time, $\text{MSD}\sim t^{\alpha}$ (Fig. 4a). These trajectories can
be unraveled as described in the previous section, to remove the effect of
network confinement. However, because the unraveling procedure assumes
diffusive transport when computing node passage probabilities, the resulting
unraveled MSD curves are notably deformed (Fig. 4b). Consequently, the
goodness-of-fit metric (comparing to a linear, diffusive scaling) is
substantially higher for these trajectories, serving as an indicator of non-
diffusive transport.

An alternative approach for distinguishing fLM from diffusive trajectories on
a network involves analysis of the velocity autocorrelation function for the
unraveled trajectories. For velocities computed over time-interval $\delta$,
the autocorrelation function is defined by Weber  _et al._ (2010b)

$\begin{split}C^{(\delta)}_{v}(t)=\frac{1}{\delta^{2}}\left<[x(t+\delta)-x(t)]\cdot[x(\delta)-x(0)]\right>,\end{split}$
(10)

where the average is taken over both time and the particle ensemble. This
function exhibits a negative peak at $t=\delta$, rising back towards zero in a
polynomial fashion for $t>\delta$. The normalized negative peak
$C_{v}^{(\delta)}(\delta)/C_{v}^{(\delta)}(0)$ becomes deeper for more
subdiffusive motion (lower $\alpha$) and, for fractional Langevin motion, is
independent of the choice of window size $\delta$ in defining the
velocitySabri  _et al._ (2020); Lucas  _et al._ (2014); Weber  _et al._
(2012b). The scaling exponent $\alpha$ can be extracted directly from this
value according to the expression Weber  _et al._ (2010b):

$\begin{split}\left|\frac{C_{v}^{(\delta)}(\delta)}{C_{v}^{(\delta)}(0)}\right|=1-2^{\alpha-1}.\end{split}$
(11)

We plot the normalized velocity autocorrelation for unraveled trajectories of
Brownian particles and particles undergoing fLM in Fig. 5a,b. For Brownian
particles, the velocity correlations are flat for $t>\delta$, emphasizing that
the unraveling procedure removes the effect of confinement and dead-ends in
the network, which would be expected to yield negative peaks Weber  _et al._
(2012b). For particles undergoing fLM, the unraveling decreases the magnitude
of the negative peaks at larger step-sizes, consistent with the turning upward
of the MSD curves for these trajectories. However, the form of the velocity
autocorrelation function is still clearly non-diffusive, demonstrating that
this metric can be used to diagnose fLM versus Brownian motion on a network.

Figure 6: Unraveling of membrane protein (Halo-TA) trajectories in the
peripheral ER networks of an example COS7 cell. (a) Extracted network
structure (green) superimposed on a confocal image with fluorescently labeled
ER marker. Trajectories within $1$ sec of this snapshot are shown as colored
circles. Inset shows single trajectory, projected onto network edges
(magenta). (b) Ensemble-averaged MSD for raw trajectories (green), projected
trajectories (blue), and unraveled trajectories (red). Dashed magenta line
indicates apparent subdiffusion for particle trajectories due to confinement
in network structure. Dashed black line shows linear scaling with the
estimated diffusivity. Inset: residual function for unraveling with different
diffusivity values, with self-consistent estimate marked in red. (c)
Normalized velocity autocorrelation for unraveled trajectories, with estimated
scaling exponent $\alpha_{\text{est}}$ is obtained from the first negative
peak. Errors for $D_{\text{est}}$ and $\alpha_{\text{est}}$ are standard
deviations obtained from bootstrapping over individual unraveled trajectories.

To generate an estimate of the fLM scaling exponent, we use Eq. 11 at
$\delta=\Delta t$ (a single trajectory step). As shown in Fig. 5c, unraveled
trajectories of simulated Brownian particles on a network give an accurate
estimate of $\alpha_{\text{est}}\approx 1.0\pm 0.01$ (mean $\pm$ standard
deviation) based on the first negative peak in the velocity autocorrelation
functions. When the same analysis is applied to simulated fLM trajectories
with $\alpha=0.6$, the unraveling procedure gives a somewhat biased estimate
of the scaling exponent ($\alpha_{\text{est}}=0.54\pm 0.025$), but still
clearly differentiates between fLM and Brownian motion. Simulations of fLM
with a range of different scaling exponents (Fig. 5c) indicate that accurate
estimates of the scaling exponent can be obtained for $\alpha\gtrsim 0.7$,
when individual step sizes are set to $10\%$ of the network edge length. For
larger step sizes, this approach overestimates $\alpha$, due to the bias
towards diffusive motion introduced in the unraveling process whenever the
particle passes a node. However, particles undergoing fLM still yield an
estimated $\alpha_{\text{est}}$ that is significantly less than $1$, allowing
the trajectories to be distinguished from those of Brownian particles. By
contrast, if the velocity autocorrelation is analyzed for raw trajectories
confined on a network, then even Brownian particles ($\alpha=1$) are likely to
yield an estimated scaling value $\alpha_{\text{est}}$ that is below the
expected range (shaded region in Fig. 5c). Thus, network confinement can make
Brownian motion appear as fLM, and the unraveling procedure removes this
confounding effect.

###  II.5 Application to membrane proteins on ER tubules

As a direct application of this methodology to experimental data, we consider
single particle trajectories of a model ER transmembrane protein (a HaloTag
fused to the minimal targeting domain of Sec61b, Halo-TA), moving within the
peripheral ER tubules of COS7 cells. Widefield images of the ER network as
well as individual fluorescently labeled Halo-TA proteins were collected
simultaneously with rapid time resolution ($100$Hz). Because ER networks are
dynamic, rearranging over tens of seconds Friedman  _et al._ (2010); Lin  _et
al._ (2017); Lee and Chen (1988); Lee  _et al._ (1989), peripheral network
structures were extracted at $1$ second intervals to serve as the underlying
domain geometry for the particle trajectories. Details on imaging, image
analysis, and particle tracking methods are provided in the Methods section
(Appendix A).

For each network structure, trajectory segments collected within $\pm 1$ sec
of the corresponding image time were mapped onto the network by projecting to
the nearest point along the edges (Fig. 6a). All projected trajectories imaged
within a given cell (corresponding to approximately $1$ minute total imaging
time) were then collected together for analysis. As seen in Fig. 2b, both the
raw and projected trajectories exhibited an apparently subdiffusive MSD
scaling with time, as expected for simulated diffusive particles on a sample
ER network (Fig. 1c). The close correspondence between raw and projected
trajectories indicates that the proteins were in fact largely confined to the
visible ER tubules, with most of their motion occuring along the tubule axis.
In Appendix B we explore further the effect of particles moving along the
circumference of the ER tubules rather than directly along the axis. We show,
using simulations, that the narrow 50 nm radius of ER tubules Schroeder  _et
al._ (2019) implies that these perpendicular particle displacements will not
significantly affect the analysis.

Projected trajectories were unraveled as described in Section II.3 to extract
the underlying free-particle dynamics. The mean squared displacement
$\text{MSD}_{\text{ur}}$ of the unraveled trajectory closely matched the
expected $2D_{\text{est}}t$ behavior (Fig. 6b) with a diffusivity of
$D_{\text{est}}=1.54\pm 0.1\mu\text{m}^{2}/\text{s}$, where the standard
deviation in the estimate is obtained by bootstrapping over individual
particle trajectories. The velocity autocorrelation function of the unraveled
trajectories is consistent with the expected behavior for diffusive particles,
showing very little correlation for $t>\delta$ (Fig. 6c). Estimating the
scaling exponent $\alpha$ from these velocity autocorrelations (as described
for simulated trajectories above) yields a value of
$\alpha_{\text{est}}=1.02\pm 0.01$, with standard deviation again obtained by
bootstrapping.

The analysis was repeated for particle trajectories and network structures
obtained from $11$ different cells (Fig. 7), with resulting diffusivity
estimates $\left<D\right>=1.45\pm 0.24\mu\text{m}^{2}/\text{s}$ (mean $\pm$
standard deviation). The extracted diffusivity is similar to values recently
estimated from bulk measurements of photoactivated protein spreading in the ER
Konno  _et al._ (2021). The scaling estimates from the velocity
autocorrelation of unravelled trajectories in each cell yielded
$\left<\alpha\right>=1.00\pm 0.02$ (standard deviation across cells),
consistent with diffusive motion of the Halo-TA ER membrane protein.

We note that the average edge length of the ER networks extracted in this
study was $1.0\pm 0.7\mu$m (mean $\pm$ standard deviation). Given our
diffusivity estimate, we would expect the unraveling procedure to be accurate
for step sizes $\Delta t\lesssim 0.07\text{sec}$, a condition which is
satisfied by the $100$Hz imaging rate.

Figure 7: Analysis of trajectories for Halo-TA ER membrane proteins from
multiple COS7 cells. (a) Estimated diffusivity $D_{\text{est}}$. (b) Estimated
scaling exponent $\alpha_{\text{est}}$. $N=13$ independent cells were
analyzed.

##  III Discussion

Traditional approaches to quantifying single particle trajectories involve
using the mean squared displacement or velocity autocorrelation function to
characterize the diffusive or subdiffusive motion of the particles. However,
for particles trapped inside the complex spatial geometry of intracellular
structures, these approaches are confounded by the effects of confinement.
Such effects can make diffusive motion appear subdiffusive and lead to
underestimation of the particle motility. In this study, we demonstrate a
methodology for decoupling diffusive particle behavior from the effect of
spatial confinement using direct observation of the structure, in the specific
case of particles confined within a tubular network.

Our approach involves unraveling a particle trajectory on a network by
sampling the hypothetical path of that particle if it were to undergo
unconstrained diffusion on an infinite line. In classic simulations of a
Brownian particle, each step of the particle represents an ensemble of
possible paths leading from the starting position to the next one. When all
likely paths stay within a single edge of the network, such a classically
sampled step remains accurate. However, if there is a non-trivial chance of
passing a network junction during that time-step, then it becomes necessary to
consider the excursions of the particle to different edges surrounding that
junction. We develop an algorithm that leverages the statistical weighting of
these excursions to ‘fold’ one-dimensional diffusive paths onto the edges
around a network junction. This serves as a mathematically accurate method for
discrete-time Brownian dynamics simulations of particles along a network of
tubules.

We then take a Bayesian approach to inverting this algorithm. Namely, starting
with an observed trajectory of discrete steps on a network, we sample the
unconfined one-dimensional trajectory that would lead to the observed particle
positions. The resulting ‘unraveled’ trajectories of simulated Brownian
particles exhibit all the features of classic Brownian motion and can be
analyzed using the usual quantification of MSD and velocity autocorrelations.
This approach makes it possible to extract the diffusivity of the particle as
well as to differentiate between diffusive and subdiffusive motion.

Simulations of particles undergoing fractional Langevin motion on network
structures indicate that applying our unraveling algorithm can accurately
differentiate such dynamics from Brownian diffusion. For sufficiently small
step-sizes, analysis of the velocity autocorrelation function for unraveled
trajectories yields an accurate estimate of the fractional scaling exponent
$\alpha$. Furthermore, although raw Brownian trajectories may appear
subdiffusive due to confinement within the network, unraveling the
trajectories allows for accurate estimates of Brownian scaling behavior
($\alpha=1$) from the velocity autocorrelation function. Hence, the unraveling
technique can be used both to diagnose particle motion as diffusive (as
opposed to fLM) and to estimate the associated diffusivity, while removing the
confounding effects of network confinement.

Our proposed approach for uncoupling diffusive transport from network
confinement is subject to certain inherent limitations. It is applicable
specifically to particles moving on a network of narrow tubules, so that
motion along individual edges can be approximated as one-dimensional. It also
requires sufficiently frequent observations of the particle position so that
each particle can be assumed to pass at most one node junction during each
time-step. This is equivalent to requiring that the particle steps be shorter
than the typical edge length. Expanding the proposed algorithm to find the
correct statistical weighting for paths that pass multiple nodes over the
course of a time-step serves as a promising avenue for future work.

As an example application of the unraveling algorithm, we analyze the
trajectories of membrane-bound proteins in the endoplasmic reticulum (ER). The
peripheral ER in many cell types forms a tubular network structure, which
hosts a variety of proteins that diffuse within its lumen and membrane. Among
its many biologically crucial functions, it is thought to serve as an
intracellular network for the sorting and delivery of ions, lipids, and
secreted proteins. Currently, the mechanism and spatial dynamics of these
phenomena are challenging to observe, since the complex structure masks the
interpretation of dynamic data collected within the organelle. Quantifying
particle motion in the ER is crucial to elucidating the impact of its
reticulated structure on the kinetics of biochemical pathways embedded within
this organelle. Such studies constitute a first step towards generating a
mechanistic understanding of ER-localized biological processes and their
regulation.

Trajectories of single transmembrane proteins on mammalian ER networks were
collected at a sampling rate of 100Hz, sufficiently rapid to ensure multiple
steps of each protein between junction passage events. Traditional analysis of
the MSD for raw particle trajectories was confounded by confinement within the
network, as indicated by an apparent subdiffusive scaling estimate. However,
application of the unraveling algorithm demonstrated that the dynamics of
these particles is consistent with classic Brownian motion, with a diffusivity
of $D\approx 1.5\mu\text{m}^{2}/\text{s}$. This estimate is compatible with
bulk measurements of protein diffusivity in the ER membrane Konno  _et al._
(2021).

The algorithm proposed here is not specific to protein motion in the ER, but
rather is a general mathematical tool for analyzing ensembles of diffusive
trajectories on networks embedded in physical space. Such spatial networks
have been analyzed in the context of roads, power grids, and venation patterns
as well as intracellular structures Barthélemy (2011). In cell biology,
mitochondrial networks Viana  _et al._ (2020) and the dendritic trees of
neurons Sartori  _et al._ (2020) provide additional examples of network
structures where diffusion can play an important role in protein transport.
Our focus here is specifically on diffusive dynamics over the network, but the
unraveling method also provides a potential tool for diagnosing when particle
motion is non-diffusive. We show that it can accurately distinguish fractional
Langevin motion from true diffusion. Future work will pursue the possibility
of distinguishing other movement patterns, such as directed flow along edges
Holcman  _et al._ (2018), or continuous-time random walks associated with
sporadic binding Malchus and Weiss (2010).

Because the unraveling approach provides a quantitative estimate of particle
diffusivity, it has potential applications in distinguishing the effect of
external perturbations on either network structure or the dynamics of the
particles themselves. For example, mutations in ER morphogen proteins can
alter the radius of ER tubules Konno  _et al._ (2021), or the density of
tubule junctions Wang  _et al._ (2016). However, the consequences of these
morphological changes on protein diffusivity in the ER remain unclear, and
could be approached by unraveling membrane and luminal protein trajectories
according to the methods proposed here. Overall, the ability to quantify
trajectories by removing the complicating effects of confinement is critical
to understanding the behavior of particles trapped in complex reticulated
geometries, in cell biology and beyond.

## Acknowledgements

This work was initiated at the QCB Cell Modeling Hackathon, supported by
Wallace Marshall and the Quantitative Cell Biology Network (sponsored by NSF
grant MCB‐1411898). The authors thank Jonathon Nixon-Abell, Zubenelgenubi
Scott, and the members of the Hackathon for helpful discussions. They also
thank Jonathon Nixon-Abell for collecting and curating some of the data used
in this study. Funding for this work was provided by NSF grant MCB-2034482
(EFK), a Cottrell Scholar Award from the Research Corporation for Science
Advancement (EFK), a Chancellor Funded Research Grant (EFK), and direct
support of the Howard Hughes Medical Institute (CJO and JL-S).

## References

  * Wirtz (2009) D. Wirtz, Ann Rev Biophys 38, 301 (2009).
  * Squires and Mason (2010) T. M. Squires and T. G. Mason, Annu Rev Fluid Mech 42, 413 (2010).
  * Elf  _et al._ (2007) J. Elf, G.-W. Li, and X. S. Xie, Science 316, 1191 (2007).
  * Parry  _et al._ (2014) B. R. Parry, I. V. Surovtsev, M. T. Cabeen, C. S. O’Hern, E. R. Dufresne, and C. Jacobs-Wagner, Cell 156, 183 (2014).
  * Holcman  _et al._ (2018) D. Holcman, P. Parutto, J. E. Chambers, M. Fantham, L. J. Young, S. J. Marciniak, C. F. Kaminski, D. Ron, and E. Avezov, Nat Cell Biol 20, 1118 (2018).
  * Guo  _et al._ (2014) M. Guo, A. J. Ehrlicher, M. H. Jensen, M. Renz, J. R. Moore, R. D. Goldman, J. Lippincott-Schwartz, F. C. Mackintosh, and D. A. Weitz, Cell 158, 822 (2014).
  * Gal  _et al._ (2013) N. Gal, D. Lechtman-Goldstein, and D. Weihs, Rheol Acta 52, 425 (2013).
  * Koslover  _et al._ (2016) E. F. Koslover, C. K. Chan, and J. A. Theriot, Biophys J 110, 700 (2016).
  * Lin  _et al._ (2016) C. Lin, M. Schuster, S. C. Guimaraes, P. Ashwin, M. Schrader, J. Metz, C. Hacker, S. J. Gurr, and G. Steinberg, Nature communications 7, 11814 (2016).
  * Chen  _et al._ (2015) K. Chen, B. Wang, and S. Granick, Nat Mater 14, 589 (2015).
  * Bálint  _et al._ (2013) Š. Bálint, I. V. Vilanova, Á. S. Álvarez, and M. Lakadamyali, P Natl Acad Sci 110, 3375 (2013).
  * Weihs  _et al._ (2006) D. Weihs, T. G. Mason, and M. A. Teitell, Biophys J 91, 4296 (2006).
  * Weber  _et al._ (2012a) S. C. Weber, A. J. Spakowitz, and J. A. Theriot, P Natl Acad Sci 109, 7338 (2012a).
  * Fodor  _et al._ (2015) É. Fodor, M. Guo, N. Gov, P. Visco, D. Weitz, and F. van Wijland, EPL (Europhysics Letters) 110, 48005 (2015).
  * Kollmannsberger and Fabry (2011) P. Kollmannsberger and B. Fabry, Annu Rev Mater Res 41, 75 (2011).
  * Wang  _et al._ (2006) Y. Wang, R. H. Austin, and E. C. Cox, Phys Rev Lett 97, 048302 (2006).
  * Li  _et al._ (2015) H. Li, S.-X. Dou, Y.-R. Liu, W. Li, P. Xie, W.-C. Wang, and P.-Y. Wang, J Am Chem Soc 137, 436 (2015).
  * Etoc  _et al._ (2018) F. Etoc, E. Balloul, C. Vicario, D. Normanno, D. Liße, A. Sittner, J. Piehler, M. Dahan, and M. Coppey, Nat Mater 17, 740 (2018).
  * Yamada  _et al._ (2000) S. Yamada, D. Wirtz, and S. C. Kuo, Biophys J 78, 1736 (2000).
  * Hoffman  _et al._ (2006) B. D. Hoffman, G. Massiera, K. M. Van Citters, and J. C. Crocker, P Natl Acad Sci 103, 10259 (2006).
  * Weber  _et al._ (2012b) S. C. Weber, M. A. Thompson, W. E. Moerner, A. J. Spakowitz, and J. A. Theriot, Biophys J 102, 2443 (2012b).
  * Sabri  _et al._ (2020) A. Sabri, X. Xu, D. Krapf, and M. Weiss, Phys Rev Lett 125, 058101 (2020).
  * Lucas  _et al._ (2014) J. S. Lucas, Y. Zhang, O. K. Dudko, and C. Murre, Cell 158, 339 (2014).
  * Wang  _et al._ (2012) B. Wang, J. Kuo, S. C. Bae, and S. Granick, Nat Mater 11, 481 (2012).
  * Lampo  _et al._ (2017) T. J. Lampo, S. Stylianidou, M. P. Backlund, P. A. Wiggins, and A. J. Spakowitz, Biophys J 112, 532 (2017).
  * Mogre  _et al._ (2020) S. S. Mogre, A. I. Brown, and E. F. Koslover, Phys Biol 17, 061003 (2020).
  * Duits  _et al._ (2009) M. H. Duits, Y. Li, S. A. Vanapalli, and F. Mugele, Phys Rev E 79, 051910 (2009).
  * Chechkin  _et al._ (2017) A. V. Chechkin, F. Seno, R. Metzler, and I. M. Sokolov, Physical Review X 7, 021002 (2017).
  * Adler  _et al._ (2019) J. Adler, I.-M. Sintorn, R. Strand, and I. Parmryd, Communications biology 2, 1 (2019).
  * Saxton (1995) M. J. Saxton, Biophys J 69, 389 (1995).
  * Viana  _et al._ (2020) M. P. Viana, A. I. Brown, I. A. Mueller, C. Goul, E. F. Koslover, and S. M. Rafelski, Cell Syst 10, 287 (2020).
  * Schwarz and Blower (2016) D. S. Schwarz and M. D. Blower, Cell Mol Life Sci 73, 79 (2016).
  * Perkins and Allan (2021) H. T. Perkins and V. Allan, Cells 10, 2341 (2021).
  * Westrate  _et al._ (2015) L. Westrate, J. Lee, W. Prinz, and G. Voeltz, Annu Rev Biochem 84, 791 (2015).
  * Friedman and Voeltz (2011) J. R. Friedman and G. K. Voeltz, Trends Cell Biol 21, 709 (2011).
  * Chen  _et al._ (2013) S. Chen, P. Novick, and S. Ferro-Novick, Curr Opin Cell Biol 25, 428 (2013).
  * Griffing (2010) L. R. Griffing, Biochem Soc T 38, 747 (2010).
  * Barlowe and Helenius (2016) C. Barlowe and A. Helenius, Annu Rev Cell Dev Bi 32, 197 (2016).
  * Brown  _et al._ (2020) A. I. Brown, L. M. Westrate, and E. F. Koslover, Sci Rep 10, 1 (2020).
  * Scott  _et al._ (2021) Z. C. Scott, A. I. Brown, S. S. Mogre, L. M. Westrate, and E. F. Koslover, Eur Phys J E 44, 1 (2021).
  * Stadler  _et al._ (2018) L. Stadler, K. Speckner, and M. Weiss, Biophys J 115, 1552 (2018).
  * Dora and Holcman (2020) M. Dora and D. Holcman, Proceedings of the Royal Society B 287, 20200493 (2020).
  * Dayel  _et al._ (1999) M. J. Dayel, E. F. Hom, and A. S. Verkman, Biophys J 76, 2843 (1999).
  * Siggia  _et al._ (2000) E. D. Siggia, J. Lippincott-Schwartz, and S. Bekiranov, Biophysical journal 79, 1761 (2000).
  * Nehls  _et al._ (2000) S. Nehls, E. L. Snapp, N. B. Cole, K. J. Zaal, A. K. Kenworthy, T. H. Roberts, J. Ellenberg, J. F. Presley, E. Siggia, and J. Lippincott-Schwartz, Nature cell biology 2, 288 (2000).
  * Konno  _et al._ (2021) T. Konno, P. Parutto, D. M. Bailey, V. Davì, C. Crapart, M. A. Awadelkareem, C. Hockings, A. Brown, K. M. Xiang, A. Agrawal, _et al._ , bioRxiv (2021).
  * Santamaria  _et al._ (2006) F. Santamaria, S. Wils, E. De Schutter, and G. J. Augustine, Neuron 52, 635 (2006).
  * Sartori  _et al._ (2020) F. Sartori, A.-S. Hafner, A. Karimi, A. Nold, Y. Fonkeu, E. M. Schuman, and T. Tchumatchenko, Cell Rep 33, 108391 (2020).
  * Deng and Barkai (2009) W. Deng and E. Barkai, Phys Rev E 79, 011112 (2009).
  * Jeon and Metzler (2010) J.-H. Jeon and R. Metzler, Phys Rev E 81, 021103 (2010).
  * Stauffer and Aharony (2018) D. Stauffer and A. Aharony,  _Introduction to percolation theory_ (CRC press, 2018).
  * Ross (2014) S. M. Ross,  _Introduction to probability models_ (Academic press, 2014).
  * Metzler and Klafter (2000) R. Metzler and J. Klafter, Physics reports 339, 1 (2000).
  * Weber  _et al._ (2010a) S. C. Weber, A. J. Spakowitz, and J. A. Theriot, Phys Rev Lett 104, 238102 (2010a).
  * Bressloff and Newby (2013) P. C. Bressloff and J. M. Newby, Rev Mod Phys 85, 135 (2013).
  * Kroese and Botev (2013) D. Kroese and Z. Botev, Analysis, modeling and simulation of complex structures 2 (2013).
  * Weber  _et al._ (2010b) S. C. Weber, J. A. Theriot, and A. J. Spakowitz, Phys Rev E 82, 011913 (2010b).
  * Friedman  _et al._ (2010) J. R. Friedman, B. M. Webster, D. N. Mastronarde, K. J. Verhey, and G. K. Voeltz, J Cell Biol 190, 363 (2010).
  * Lin  _et al._ (2017) C. Lin, R. R. White, I. Sparkes, and P. Ashwin, Biophys J 113, 214 (2017).
  * Lee and Chen (1988) C. Lee and L. B. Chen, Cell 54, 37 (1988).
  * Lee  _et al._ (1989) C. Lee, M. Ferguson, and L. B. Chen, The Journal of cell biology 109, 2045 (1989).
  * Schroeder  _et al._ (2019) L. K. Schroeder, A. E. Barentine, H. Merta, S. Schweighofer, Y. Zhang, D. Baddeley, J. Bewersdorf, and S. Bahmanyar, Journal of Cell Biology 218, 83 (2019).
  * Barthélemy (2011) M. Barthélemy, Physics Reports 499, 1 (2011).
  * Malchus and Weiss (2010) N. Malchus and M. Weiss, Biophys J 99, 1321 (2010).
  * Wang  _et al._ (2016) S. Wang, H. Tukachinsky, F. B. Romano, and T. A. Rapoport, elife 5, e18605 (2016).
  * MATLAB (2018) MATLAB,  _version 9.5 (R2018b)_ (The MathWorks Inc., Natick, Massachusetts, 2018).
  * Sengupta  _et al._ (2013) P. Sengupta, T. Jovanovic-Talisman, and J. Lippincott-Schwartz, Nature protocols 8, 345 (2013).
  * Nixon-Abell  _et al._ (2016) J. Nixon-Abell, C. J. Obara, A. V. Weigel, D. Li, W. R. Legant, C. S. Xu, H. A. Pasolli, K. Harvey, H. F. Hess, E. Betzig, _et al._ , Science 354 (2016).
  * Grimm  _et al._ (2016) J. B. Grimm, B. P. English, H. Choi, A. K. Muthusamy, B. P. Mehl, P. Dong, T. A. Brown, J. Lippincott-Schwartz, Z. Liu, T. Lionnet, _et al._ , Nature methods 13, 985 (2016).
  * Berg  _et al._ (2019) S. Berg, D. Kutra, T. Kroeger, C. N. Straehle, B. X. Kausler, C. Haubold, M. Schiegg, J. Ales, T. Beier, M. Rudy, K. Eren, J. I. Cervantes, B. Xu, F. Beuttenmueller, A. Wolny, C. Zhang, U. Koethe, F. A. Hamprecht, and A. Kreshuk, Nat Methods (2019), 10.1038/s41592-019-0582-9.
  * Tinevez  _et al._ (2017) J.-Y. Tinevez, N. Perry, J. Schindelin, G. M. Hoopes, G. D. Reynolds, E. Laplantine, S. Y. Bednarek, S. L. Shorte, and K. W. Eliceiri, Methods 115, 80 (2017).
  * Jaqaman  _et al._ (2008) K. Jaqaman, D. Loerke, M. Mettlen, H. Kuwata, S. Grinstein, S. L. Schmid, and G. Danuser, Nature methods 5, 695 (2008).

##  Appendix A Methods

All code used in this work was written in Matlab MATLAB (2018). A software
package for carrying out simulations on networks and for analyzing particle
trajectories will be provided at
https://github.com/lenafabr/unravelNetworkTraj prior to publication.

###  A.1 Algorithm for Brownian motion simulations on networks

We develop an algorithm for discrete time-step simulations of Brownian
particles on a network of one-dimensional edges connected by point-like
junctions, as discussed in Sec. II.1.

At each step, for each individual particle, we consider the particle’s
position along its current edge. Without loss of generality, we can assume
that the nearest boundary (node) of that edge is at position $0$ and the
particle itself is at position $x_{0}$ along the edge. We sample the
displacement of the particle using a normal distribution $\Delta
z\sim\mathcal{N}(0,\sqrt{2D\Delta t})$. We then treat the particle as if it
were on an infinite line and consider, of all possible paths that start at
$x_{0}$ and end at $x_{0}+\Delta z$, what is the probability that the
particle’s path crossed $0$ during that timestep. This probability
$p_{\text{pass}}$ is derived as follows.

We start with the Green’s function for unconstrained diffusion:

$\begin{split}G(x|x_{0},\Delta t)=\frac{1}{\sqrt{4\pi
Dt}}e^{-(x-x_{0})^{2}/{4Dt}},\end{split}$ (12)

which gives the spatial probability density of a particle ending its path at
position $x$ after time $\Delta t$, given that it started at position $x_{0}$
at time $0$. The distribution of first passage times for a particle starting
at $x_{0}$ to hit $0$ for the first time at time $t$ is given by

$\begin{split}J(t;x_{0})=D\left.\frac{\partial}{\partial
x}G(x|x_{0};t)\right|_{x=0}.\end{split}$ (13)

Because diffusion is Markovian (memory-less), if the particle first hits $0$
at time $t$, we know the probability density of ending up at position $x$ is
simply given by $G(x|0,\Delta t-t)$. From there, we can calculate the
conditional probability that a particle hit $0$ at some point in the path,
given that it ended at $x$, by a simple application of Bayes’ rule:

$\begin{split}p_{\text{pass}}=\frac{\int_{0}^{\Delta t}J(t;x_{0})G(x|0,\Delta
t-t)dt}{G(x|x_{0},\Delta t)}.\end{split}$ (14)

Here, the numerator is the joint probability density of hitting $0$ at some
point during the time-step and then ending the path at $x$ and the denominator
is the overall probability density of ending at $x$. Plugging in Eq. 12 and 13
gives the formula for $p_{\text{pass}}$ stated in the main text (Eq. 1).

Every path of an unconstrained particle can be mapped to a set of ‘folded’
paths for the equivalent particle on the network. Each passage of the node at
$0$ involves selecting which of the adjacent network edges will serve as the
axis along which the network-bound particle continues to move. For a memory-
less Brownian particle that hits the node, each of the edges (including the
one from whence it came) are equally likely to be selected. We know the
unconstrained particle ends its timestep at distance $|x_{0}+\Delta z|$ from
$0$. Hence, the network-confined particle is placed at this distance from the
node along the randomly selected edge, completing its time-step.

###  A.2 Distribution of diffusing particles on a triskelion

We compute explicitly the solution to the diffusion equation on a triskelion
structure (Fig. 1b) with reflecting boundaries at the tips and continuity of
the concentrations at the degree-3 node. The initial condition is assumed to
be a delta function at the junction. The approach taken is analogous to the
equivalent problem with absorbing boundaries, described in more detail in Ref.
Scott  _et al._ (2021). Specifically, after a Laplace transform in time
($t\rightarrow s$), the concentration profile $\hat{c}_{k}(x,s)$ on edge $k$
is given by

$\begin{split}\hat{c}_{k}(x,s)=\frac{1}{D\alpha}\frac{\cosh\alpha
x}{\cosh\alpha\ell_{k}}\left[\sum_{j=1}^{3}\tanh\alpha\ell_{j}\right]^{-1}\end{split}$
(15)

where $D$ is the particle diffusivity, $\alpha=\sqrt{s/D}$, and $\ell_{j}$ is
the length of the $j^{\text{th}}$ edge.

The Laplace transform is then inverted using a Bromwich integral evaluated
with the Cauchy residue theorem, yielding the time-dependent concentration
profile

$\begin{split}c(x,t)&=\sum_{p}r_{p}e^{-Du_{p}^{2}t},\end{split}$ (16)

where the poles $s_{p}=-Du_{p}^{2}$ can be found by taking the roots of the
following equation,

$\begin{split}\cos\ell_{k}u_{p}&\left[\sum_{j}\tan\ell_{j}u_{p}\right]=0,\end{split}$
(17)

and $r_{p}$ are the residues of Eq. 15 at those poles.

The resulting expressions are plotted in Fig. 1b and compared to the
distribution of $10^{7}$ simulated particles at time $t=0.8\ell^{2}/D$ (after
starting at the three-fold junction) where $\ell$ is the length of the
shortest edge.

###  A.3 Parameters for simulations and analysis

####  A.3.1 Brownian motion on networks

The simulations are tested on synthetic networks, as well as network
structures extracted from images of the ER. For comparing particle
distributions on the triskelion network (Fig. 1b), $10^{7}$ particles were
simulated with timesteps of $\Delta tD/\ell^{2}=0.01$, for $80$ steps. All
particles were started on the central node.

The synthetic honeycomb network in Fig. 1c was created by cropping a standard
honeycomb lattice with $N=15$ cells in each dimension to a circle of radius
$1$. All reported length units for simulations on synthetic networks are
normalized by the single edge length $\ell$ in the network. The decimated
honeycomb structure was obtained by removing $30\%$ of all network edges,
randomly selected in such a way that the network maintains a single connected
component. While all the example networks used in this study are two-
dimensional, the simulations and unraveling code also applies to 3D networks.

Brownian simulations for the MSD comparison were run with $100$ particles,
using dimensionless timesteps of $\Delta tD/\ell^{2}=0.01$, for up to
$N_{t}=3\times 10^{4}$ time steps. The mean squared dispacement (MSD) is
computed as $\text{MSD}=\left<|\vec{x}(t)-\vec{x}(0)|^{2}\right>$, where the
average is done over non-overlapping time-windows of a single trajectory and
over all particle trajectories.

####  A.3.2 Trajectory unraveling analysis

The algorithm for unraveling observed trajectories on a network to sample the
corresponding trajectories on an infinite line is described in Sec. II.3. We
perform the unraveling procedure with $20$ values of diffusivity in the range
$0.2-4$ (compared to the simulated value of $D=1$). Residuals are calculated
by interpolating each unravelled MSD curve onto $100$ logarithmically spaced
time-points and then applying Eq. 9. The self-consistent diffusivity is then
found by estimating the value of $D$ that yields a local minimum in the
normalized residuals curve (eg: Fig. 2c) via a cubic spline interpolation.

###  A.4 Simulation of fractional Langevin motion on networks

The long-tailed memory kernel of fractional Langevin motion implies that
individual steps of such particles are not Markovian, necessitating an
alternate approach to the simulation. We take the approach of assuming that
individual particles are undergoing fLM in two dimensions, but that rapid
reflections from the walls of a tubule prevent their perpendicular motion. Our
overall approach is analogous to that described previously Jeon and Metzler
(2010) for simulating the fractional Langevin equation under confinement, in
the limit of no inertia (as appropriate for highly viscous intracellular
fluids).

Specifically, we consider the power-law memory kernel Weber  _et al._ (2010b)

$\begin{split}K(t-t^{\prime})=\frac{(2-\alpha)(1-\alpha)}{|t-t^{\prime}|^{\alpha}}\end{split}$
(18)

and generate 2D stochastic forces $F_{x}^{(B)},F_{y}^{(B)}$ as fractional
Gaussian noise Kroese and Botev (2013) with time correlations in accordance
with the fluctuation dissipation theorem:

$\begin{split}\left<F^{(B)}_{i}(t)F^{(B)}_{j}(t^{\prime})\right>=\gamma
k_{b}TK(t-t^{\prime})\delta_{ij},\end{split}$

where $\gamma$ is the friction coefficient and $k_{b}T$ the thermal energy of
the particles. The particle motion obeys the overdamped Langevin equation:

$\begin{split}0=-\gamma\int_{0}^{t}K(t-t^{\prime})\vec{v}(t^{\prime})dt^{\prime}+\vec{F}^{(b)}(t)\end{split}$
(19)

where $\vec{v}(t^{\prime})$ is the instantaneous particle velocity. In the
absence of confinement, the MSD for such particles is described by the
following power law Weber  _et al._ (2010b):

$\begin{split}\left<x^{2}\right>&=2D_{\alpha}t^{\alpha},\\\
D_{\alpha}&=\frac{k_{b}T/\gamma}{(2-\alpha)(1-\alpha)\Gamma(\alpha+1)\Gamma(1-\alpha)}\end{split}$
(20)

To propagate forward the particles via discrete time-steps, we take the
approach described in Ref. Deng and Barkai (2009). Namely, we discretize the
integral in Eq. 19 and solve for the next trial spatial step
$\Delta\widehat{\vec{x}}_{n+1}$ based on all prior steps $\Delta\vec{x}_{n+1}$
and the fractional Brownian force $\vec{F}^{(B)}_{n+1}$.

The trial step $\Delta\widehat{\vec{x}}$ is then projected along the direction
of the current edge containing the particle (removing the perpendicular
component). If the remaining step takes the particle past a degree-1 node, it
is reflected off the node, back along the same edge. If it steps past a
degree-3 node, then it is placed randomly on one of the adjacent edges
(including the one it came from). For our simulated networks, edges are evenly
distributed around each such junction, so that multiple reflections in the
tight junctional space are expected to be equally likely to bounce the
particle to any edge. Application of this algorithm to more complex
structures, including bent degree-2 nodes and degree-3 nodes with uneven edge
distribution would require further refinement of this algorithm (left to
future work) to account for biased probabilities of entering subsequent edges.

The final saved step $\Delta{\vec{x}}$ is then computed as the two-dimensional
difference between the final particle position and its prior position at the
start of the step. This saved step is used in the computation of all future
particle displacements, allowing the particle to maintain a memory of its
prior speed and direction of motion.

We note that this simulation approach gives the expected mean squared
displacement of $2D_{\alpha}t^{\alpha}$ when particles are placed along a
single long edge or along a fully connected honeycomb lattice.

###  A.5 Experimental methods

####  A.5.1 Cell culture, plating, and transfection

COS7 cells were purchased from ATCC and maintained in phenol red-free
Dulbecco’s modified Eagle medium supplemented with 10% (v/v) FBS, 2mM
L-glutamine, 100U/ml penicillin and $100\mu\text{g/ml}$ streptomycin at 37°C
and 5% CO2. All experiments were performed within 40 passages of the initial
thaw, and passaging was performed using phenol red-free trypsin (Corning).

High tolerance, 25mm Number 1.5 coverslips were purchased from Warner
scientific and precleaned with a modified version of a previously described
protocol Sengupta  _et al._ (2013). Briefly, the coverslips were sonicated for
12 hours in 0.1% Hellmenex™ (Sigma), followed by five washes in 300ml of
distilled water, followed by an additional 12 hour sonication in distilled
water and an additional round of five washes. Coverslips were then ethanol
sterilized in 200 proof ethanol and allowed to air dry in a clean tissue
culture hood. After cleaning, coverslips were stored in an airtight container
until ready for use.

Coverslips were pre-coated with $500\mu\text{g/ml}$ phenol red-free matrigel
(Corning). Cells were seeded to achieve 60% confluency at the time of imaging.
Transfections were performed after letting the cells adhere to the coated
glass for at least 12 hours, using Fugene6 (Promega) according to the
manufacturer’s protocol. Each coverslip was transfected with 750 ng of PrSS-
mEmerald-KDEL to label the ER structure in the 488 channel, and with 250ng of
HaloTag-Sec61b-TA for tracking on the second camera. The HaloTag construct is
a minimal targeting domain from Sec61b fused with a flexible linker and the
HaloTag, as this construct has been shown to be biochemically inert within the
ERNixon-Abell  _et al._ (2016).

####  A.5.2 Labeling conditions

Immediately before imaging, each sample was labeled with 10nM PA-JF646 Grimm
_et al._ (2016) in OptiMEM (ThermoFisher) for 1 minute. Following staining,
cells were immediately washed at least 5 times with 10 ml of PBS under
continuous aspiration, taking care not to let the cells contact the air. This
was followed by an additional wash in 10 ml of pre-warmed, complete medium.
Samples were then left 10 minutes in an incubator to let the cells recover
before moving immediately to the microscope.

####  A.5.3 Microscope and imaging conditions

Dual-color imaging was performed using a customized inverted Nikon Ti-E
microscope outfitted with a live imaging stage to maintain temperature, CO2
level, and relative humidity during imaging (Tokai Hit). The sample was
illuminated with two fiber-coupled 488nm and 642nm lasers (Agilent
Technologies) introduced into the system with a conventional rear-mount TIRF
illuminator. Imaging was performed with the angle of incidence of excitation
light manually adjusted beneath the critical angle as needed to produce the
most even illumination possible in the ER. If necessary, a small amount of
405nm light was introduced to increase the rate of photoconversion in sparsely
labeled cells, but in practice this was rarely needed. The 488 illumination
power was kept beneath 100 $\mu$W total in the back aperture, in order to
avoid undesirable activation of the photoconvertible dye.

Fluorescence emission light was collected with 100x $\alpha$-Plan-Apochromat
1.49 NA oil objective (Nikon Instruments). Emission light was split into dual
paths using a 565LP or 585LP dichroic mirror in a MultiCam optical splitter
(Cairn Research), and the two channels were cleaned up by passing the light
through 525/50 and 647LP filters (Chroma) placed before the camera. Final
signal was focused back onto synchronized dual iXon3 electron mutliplying
charged coupled device cameras (EM-CCD, DU-897; Andor Technology). Imaging was
performed with 5ms exposure times for 60-90 seconds at a time, and the timing
of each frame was monitored using an oscilloscope directly coupled into the
system (mean frame rate $\approx 95$Hz).

###  A.6 Image analysis for experimental data

####  A.6.1 Image preparation and preprocessing

The ER is relatively stable at the level of diffraction-limited imaging over
the time scale of a second Friedman  _et al._ (2010). Thus, we performed 10
frame (110 msec) median filter on the channel for the ER structure to increase
the signal to noise and minimize the necessary 488 radiation. Filtering was
performed for every frame, but for segmentation and downstream analysis, the
structure was analyzed approximately every second (100 frames).

####  A.6.2 Segmentation

Segmentation and analysis of the tubular ER structure was performed using a
two-step process. First, the intensity in the filtered image was made uniform
through time with a simple ratio bleach correction. The location of the ER
within the image was identified using an interactive pixel classification
workflow in Ilastik Berg  _et al._ (2019). Once the predictions for the ER
label were judged to be of sufficiently good quality, the image was made
binary using a simple threshold on the label probability exported from
ilastik.

####  A.6.3 Extracting network structure

For each binary segmented image of the ER structure, we skeletonize using the
bwmorph subroutine in MATLAB, which also identifies junction pixels.
Neighboring junction pixels are grouped into a single node located at their
center of mass. We then use the bwtraceboundary subroutine to trace out the
skeletonized edges starting from each node, until a neighboring node is
reached. This is repeated until the full skeleton has been categorized into
nodes and edges connecting specific node pairs. The path of the edge between
each pair of connected node is smoothed using cubic splines. Network
construction and manipulations is carried out using Matlab code provided at:
https://github.com/lenafabr/networktools.

####  A.6.4 Particle trajectories

Single molecule localization and tracking was performed using the TrackMate
plugin in Fiji Tinevez  _et al._ (2017); Jaqaman  _et al._ (2008). Linking
parameters were experimentally selected for each data set to minimize visible
linkage artifacts as determined by eye. Datasets were then projected onto the
simultaneously collected structure of the ER and manually curated to remove
trajectory linkages that were close in 2D but far from one another in the
underlying organelle.

For each network structure, obtained at $\approx 1.08$ sec ($100$ frame)
intervals, we identify particle trajectory segments that fall within $\pm 100$
frames of the timepoint corresponding to that structure. Those trajectories
are then projected onto the network structure by finding the nearest point
along the edge paths to each particle position. Because the time intervals
surrounding each network structure overlap with each other, each trajectory
segment enters the analysis twice by projection onto the structure before and
after its time-point.

Whenever the original particle position is more than $2$px ($\sim 0.3\mu$m)
away from the nearest network edge, the projected point is removed and the
trajectory is broken into separate segments. The trajectories are also broken
whenever a particle appears to step onto a non-adjacent edge (bypassing more
than one node) within a single time-step.

The projection step allows us to re-express a particle trajectory in terms of
the edge index and position along the edge contour at each time step.
Projected trajectories of length at least $10$ time-steps were kept for
analysis and unravelled using the same algorithm as for simulated data. An
average of 380 projected trajectories per cell, with average length 48
timesteps ($0.51$ sec) were analyzed.

##  Appendix B Effect of motion around tubule circumference

Figure 8: Effect of particle motion around tubule circumference on diffusivity
estimates. (a) Sample trajectories of particles on tubules of radius $100$nm.
(b) Estimated diffusivity $D_{\text{est}}$ normalized by simulated diffusivity
$D$, for particles moving along the network in (a), with different tubule
radii. Particles are placed at an arbitrary position around the circumference
of the tubule at each step.

Our simulations and data analysis rely on the assumption that the edges of the
network are very narrow, and can be treated as effectively one-dimensional.
However, in a realistic ER tubular network, membrane proteins are confined to
the circumference of the tubule and may thus move perpendicular to the edge
axis, as well as along the edge. Typical ER tubules have a radius of
approximately $50$nm Schroeder  _et al._ (2019).

Here, we use simulations to examine how displacements along the tubules’
circumference (orthogonal to the edge paths) might impact the diffusivity
estimation. In particular, we assume tubules are narrow enough that particles
can explore the entire circumference over an individual time-step. This is a
valid assumption for ER membrane proteins, with approximate diffusivity
$1.5\mu\text{m}^{2}/\text{s}$, where the timescale to cover the circumference
is approximately $0.008$sec. Thus, we assume at each step that the
circumferential position of the particle (around the cross-section
perpendicular to the edge path) is selected uniformly, with no correlation to
the previous steps.

When the particle position is projected in two dimensions, the displacement
orthogonal to the edge axis is in the range $-r\leq x\leq r$, where $r$ is the
tubule radius. This displacement obeys the cumulative distribution function:

$f(x)=\frac{1}{2}+\frac{1}{\pi}\arcsin{\frac{x}{r}}$ (21)

We simulated particle trajectories on the network shown in Fig. 8a for tubule
radii $r$ up to $300$nm. Each of the trajectories was then analyzed in a
manner analogous to experimental data, by projecting the particle onto the
nearest point along the network edges. After applying the unraveling
procedure, we see that the estimated diffusivity is not substantially affected
by the circumferential displacements for tubules of radius below $100$nm.
Tubules in the peripheral ER fall well within this criterion, and we thus
expect our assumption of infinitely narrow tubules to not yield substantial
bias in the analysis.


Recall that the atypical dominant weights for $\mathfrak{q}_{2}$ are of the
form $s(\varepsilon_{1}-\varepsilon_{2})$ for $s\in\frac{1}{2}\mathbb{N}$.

Proposition.  Take $\mathfrak{g}=\mathfrak{q}_{2}$ with
$\Delta^{+}=\\{\alpha\\}$. For $\lambda\in P^{+}(\mathfrak{q}_{2})$ one has

$\operatorname{sch}_{\mathfrak{h}}L(\lambda)=\left\\{\begin{array}[]{lcl}a_{\lambda}&&\text{
if $\lambda$ is typical};\\\ \sum_{i=1}^{s}a_{i\alpha}&&\text{ if
}\lambda=s\alpha,\ s\in\mathbb{N};\\\
\sum_{i=0}^{s}a_{i+\frac{1}{2}\alpha}&&\text{ if
}\lambda=(s+\frac{1}{2})\alpha,\ s\in\mathbb{N}.\end{array}\right.$

###### Proof.

If $\lambda$ is typical, the assertion follows from Theorem 6.3.1. Consider
the case when $\mu\in P^{+}(\mathfrak{q}_{2})$ is atypical and $\mu\not=0$.
Write $K(\mu)$ for the maximal finite-dimensional quotient of
$\operatorname{Ind}_{\mathfrak{b}}^{\mathfrak{q}_{2}}C(\mu)$. It is known that
$\operatorname{sch}_{\mathfrak{h}}K(\mu)=a_{\mu}$. Further, if $\mu=s\alpha$
for $s\in\mathbb{N}$, then it is known (for example see Section 7 of [9]) that
we have short exact sequences

$0\to V_{0}\to K(\alpha)\to L(\alpha)\to 0,\ \ \ \ 0\to\Pi L((s-1)\alpha)\to
K(s\alpha)\to L(s\alpha)\to 0$

where $V_{0}$ is a nontrivial extension of $\mathbb{C}$ and $\Pi\mathbb{C}$.
On the other hand if $\mu=(s+\frac{1}{2})\alpha$ then for $s=0$ we have
$K(\alpha/2)=L(\alpha/2)$ and for $s>0$ we have short exact sequences

$0\to\Pi L((s-\frac{1}{2})\alpha)\to K((s+\frac{1}{2})\alpha)\to
L((s+\frac{1}{2})\alpha)\to 0.$

From these results we can obtain the desired formulas by induction. ∎

###  9.7. Realization of $\mathscr{K}_{-}(\mathfrak{q}_{n})$

Recall that $\mathcal{F}(\mathfrak{g})_{int}$ denotes the full subcategory of
$\mathcal{F}(\mathfrak{g})$ consisting of the modules with integral weights
and that $\mathcal{C}$ denotes the full subcategory of
$\mathcal{F}(\mathfrak{h})$ consists of modules with weights lying in
$P(\mathfrak{g})$. We write $\mathcal{C}_{int}$ for the subcategory of
$\mathcal{C}$ consisting of those modules with integral weights.

#### 9.7.1.

By Theorem 9.4.9, the map $\operatorname{sch}$ gives embeddings

$\mathscr{K}_{-}(\mathfrak{g})\hookrightarrow\mathscr{K}_{-}(\mathcal{C})^{W},\
\ \
\mathscr{K}_{-}(\mathcal{F}(\mathfrak{g})_{int})\hookrightarrow\mathscr{K}_{-}(\mathcal{C}_{int})^{W}.$

Further, we have identified the image of $\mathscr{K}_{-}(\mathfrak{g})$
(resp., $\mathscr{K}_{-}(\mathcal{F}(\mathfrak{g})_{int})$) with the
subalgebra spanned by $a_{\nu}$ with $\nu\in P^{+}(\mathfrak{g})$ (resp.,
$\nu\in P^{+}(\mathfrak{g})_{int}$).

One has $a_{0}=1$; we denote by
$\operatorname{ev}_{0}:\mathscr{K}_{-}(\mathcal{C})\to\mathbb{Z}$ the co-unit
map given by $\operatorname{ev}_{0}(a_{\nu})=\delta_{0,\nu}$ for $\nu\in
P^{+}(\mathfrak{g})$. One has

$\mathscr{K}_{-}(\mathcal{C})=\mathscr{K}_{-}(\mathcal{C}_{int})^{W}\oplus\mathscr{K}_{-}(\mathcal{C}_{nint})^{W},$

where $\mathcal{C}_{nint}$ consists of modules with weights that lie in
$P(\mathfrak{g})$ and are non-integral. For
$b\in\mathscr{K}_{-}(\mathcal{C})^{W}$ and
$b^{\prime}\in\mathscr{K}_{-}(\mathcal{C}_{nint})^{W}$ one has
$bb^{\prime}=\operatorname{ev}_{0}(b)b^{\prime}$, so
$\mathscr{K}_{-}(\mathcal{C}_{nint})$ is an ideal of
$\mathscr{K}_{-}(\mathcal{C})$. By Corollary 9.5.5, for $x\not=0$,
$\operatorname{ds}_{x}(\mathscr{K}_{-}(\mathcal{C}_{nint})^{W})=0$ and the
image of $\operatorname{ds}_{x}$ lies in $\mathscr{K}_{-}(\mathcal{C}_{int})$.
This reduces a study of $\mathscr{K}_{-}(\mathcal{C})$ to a study of
$\mathscr{K}_{-}(\mathcal{C}_{int})$.

#### 9.7.2.

The ring
$\mathscr{K}_{-}(\mathcal{F}(\mathfrak{g})_{int})\otimes_{\mathbb{Z}}\mathbb{Z}[\sqrt{-1}]$
can be realized in the following way.

Let $V$ be a free $\mathbb{Z}[\sqrt{-1}]$-module with a basis $\\{v_{i}|\
i\in\mathbb{Z}\setminus\\{0\\}\\}$. Denote by $\bigwedge V$ the external ring
of $V$. This is a $\mathbb{N}$-graded supercommutative ring; we consider
$(\bigwedge V)^{spoil}$, which has

$\left(\bigwedge
V\right)^{spoil}=\displaystyle\sum_{i=0}^{\infty}\bigwedge^{2i}V\oplus\displaystyle\sum_{i=0}^{\infty}\bigwedge^{2i+1}V\varepsilon$

where, we recall, $\varepsilon$ is a formal variable satisfying
$\varepsilon^{2}=2\varepsilon=0$. This is a $\mathbb{N}$-graded commutative
and supercommutative ring (which means that $(\bigwedge
V)^{spoil}_{i}(\bigwedge V)^{spoil}_{j}=0$ if $i,j$ are odd).

We denote by $\Xi_{int}$ the set of finite multisets $\\{a_{i}\\}_{i=1}^{s}$
with $a_{i}\in\mathbb{Z}\setminus\\{0\\}$ and $a_{i}+a_{j}\not=0$ for all
$1\leq i,j\leq s$. For $A\in\Xi_{int}$ we denote by $(\bigwedge
V)^{spoil}_{A}$ the span of the elements $v_{i_{1}}\wedge
v_{i_{2}}\wedge\ldots\wedge v_{i_{p}}$ (resp., the elements $\varepsilon
v_{i_{1}}\wedge v_{i_{2}}\wedge\ldots\wedge v_{i_{p}}$) with
$\operatorname{Core}(\\{i_{1},\ldots,i_{p}\\})=A$ if $A$ has an even (resp.,
an odd) cardinality. Clearly, this gives a grading

$\left(\bigwedge V\right)^{spoil}_{A}\left(\bigwedge
V\right)^{spoil}_{B}\subseteq\left(\bigwedge V\right)^{spoil}_{A\diamond B}$

and $(\bigwedge V)^{spoil}_{\emptyset}$ is the subring of $(\bigwedge
V)^{spoil}$ generated by $v_{i}\wedge v_{-i}$. In other words $\left(\bigwedge
V\right)^{spoil}$ is $\Xi$-graded algebra. For each $k$ we consider the ideal

$J_{k}:=\displaystyle\sum_{i=k+1}^{\infty}\left(\bigwedge
V\right)^{spoil}_{i}.$

Recall that the map $\lambda\mapsto\operatorname{Nonzero}(\lambda)$ gives a
one-to-one correspondence between $P^{+}(\mathfrak{g})_{int}$ and the subsets
$S\subset\mathbb{Z}\setminus\\{0\\}$ of cardinality at most $n$. The weights
in $I_{0}$ correspond to the subsets of even cardinality; the zero weight
correspond to $\emptyset$. We define

$\psi:\mathscr{K}_{-}(\mathcal{F}(\mathfrak{g})_{int})\otimes_{\mathbb{Z}}\mathbb{Z}[\sqrt{-1}]\overset{\sim}{\longrightarrow}(\bigwedge
V)^{spoil}/J_{n}$

by setting, for $\lambda\in I_{i}\cap P^{+}(\mathfrak{g})$,

$\psi(a_{\lambda}):=\varepsilon^{i}\frac{t(\lambda)}{|t(\lambda)|}\cdot
v_{\lambda},$

where $v_{\lambda}:=v_{i_{1}}\wedge v_{i_{2}}\wedge\ldots\wedge v_{i_{k}}$,
and $\operatorname{Nonzero}(\lambda)=\\{i_{1}>i_{2}>\ldots>i_{k}\\}$.

Combining Corollary 9.4.11, Proposition 9.3.2 and Corollary 9.5.3 we obtain

#### 9.7.3.

Theorem.

  1. (i)

The map
$\psi:\mathscr{K}_{-}(\mathcal{F}(\mathfrak{g})_{int})\otimes_{\mathbb{Z}}\mathbb{Z}[\sqrt{-1}]\overset{\sim}{\longrightarrow}\left(\bigwedge
V\right)^{spoil}/J_{n}$ is an isomorphism of $\Xi_{int}$-graded spoiled super
rings.

  2. (ii)

In particular, for $A\in\Xi_{int}$, $\psi$ restricts to an isomorphism of
vector spaces
$\psi_{A}:\mathscr{K}_{-}(\mathcal{F}(\mathfrak{g})_{A})\to\left(\bigwedge
V\right)^{spoil}_{A}/J_{n}$.

  3. (iii)

The map $\operatorname{ds}_{s}:R(\mathfrak{q}_{n})_{int}\to
R(\mathfrak{q}_{n-s})_{int}$ corresponds to the natural quotient map

$\left(\bigwedge V\right)^{spoil}/J_{n}\to\left(\bigwedge
V\right)^{spoil}/J_{n-s}$

(one has $J_{n-s}\supset J_{n}$ for $s>0$).

###### Proof.

It is clear that $\psi$ is an isomorphism of $\mathbb{Z}$-modules which
preserves the $\Xi$-grading. The result on $ds_{s}$ is also clear from the
map. Thus it remains to check that $\psi$ respects multiplication.

Let $\lambda,\mu\in P^{+}(\mathfrak{g})$ with $\lambda,\mu\in I_{0}$; the case
when either is in $I_{1}$ is much easier. If $a_{\lambda}a_{\mu}=0$ then the
result is clear, thus we assume that $a_{\lambda}a_{\mu}\neq 0$. In this case,
there exists unique elements $\lambda^{\prime}\in W\lambda$ and
$\mu^{\prime}\in W\mu$ such that $\lambda^{\prime}+\mu^{\prime}$ is dominant,
and thus we have
$a_{\lambda}a_{\mu}=(-1)^{i}a_{\lambda^{\prime}+\mu^{\prime}}$, where
$E_{\lambda^{\prime}}\otimes
E_{\mu^{\prime}}=\Pi^{i}E_{\lambda^{\prime}+\mu^{\prime}}$. We may write
$T_{\lambda^{\prime}+\mu^{\prime}}=(-1)^{j}T_{\lambda^{\prime}}T_{\mu^{\prime}}$
for some $j$, and it is clear that this $j$ also satisfies

$v_{\lambda^{\prime}+\mu^{\prime}}=(-1)^{j}v_{\lambda}v_{\mu}.$

Then as in the proof of Proposition 9.3.2, we have
$E_{\lambda^{\prime}}\otimes
E_{\mu^{\prime}}=\Pi^{j+\ell}E_{\lambda^{\prime}+\mu^{\prime}}$, where
$(-1)^{\ell}=\frac{t(\lambda)t(\mu)}{t(\lambda^{\prime}+\mu^{\prime})}\in\\{\pm
1\\}$. Therefore
$a_{\lambda}a_{\mu}=(-1)^{j}\frac{t(\lambda)t(\mu)}{t(\lambda^{\prime}+\mu^{\prime})}a_{\lambda^{\prime}+\mu^{\prime}}$.
Therefore

$\psi(a_{\lambda}a_{\mu})=(-1)^{j}\frac{t(\lambda)t(\mu)}{t(\lambda^{\prime}+\mu^{\prime}}\left(\frac{t(\lambda^{\prime}+\mu^{\prime})}{|t(\lambda^{\prime}+\mu^{\prime})|}v_{\lambda^{\prime}+\mu^{\prime}}\right)=\left(\frac{t(\lambda)}{|t(\lambda)|}v_{\lambda}\right)\left(\frac{t(\mu)}{|t(\mu)|}v_{\mu}\right)$

where we have used that
$|t(\lambda^{\prime}+\mu^{\prime})|=|t(\lambda)||t(\mu)|$. ∎

For $A\in\Xi_{int}$, we define

$\overline{A}:=\left\\{\begin{array}[]{ll}\
\frac{\\#\\{i:a_{i}>0\\}-\\#\\{i|a_{i}<0\\}}{2}\text{ mod }2&\text{ for
}|A|\text{ even};\\\ \ \frac{\\#\\{i:a_{i}>0\\}-\\#\\{i|a_{i}<0\\}-1}{2}\text{
mod }2&\text{ for }|A|\text{ odd}.\end{array}\right.$

Note that we have $\overline{A\diamond B}=\overline{A}+\overline{B}$ if $|A|$
or $|B|$ is even.

#### 9.7.4.

Lemma.  The grading

$\mathscr{K}_{-}(\mathcal{F}(\mathfrak{g})_{int})_{\overline{0}}=\bigoplus\limits_{\overline{A}=\overline{0}}\mathscr{K}_{-}(\mathcal{F}(\mathfrak{g})_{A}),\
\ \
\mathscr{K}_{-}(\mathcal{F}(\mathfrak{g})_{int})_{\overline{1}}=\bigoplus\limits_{\overline{A}=\overline{1}}\mathscr{K}_{-}(\mathcal{F}(\mathfrak{g})_{A})$

defines another super ring structure on
$\mathscr{K}_{-}(\mathcal{F}(\mathfrak{g})_{int})$; note that it does not have
the structure of a spoiled super ring under this grading.

###### Proof.

This follows from the fact that
$\mathscr{K}_{-}(\mathcal{F}(\mathfrak{g})_{int})$ is $\Xi_{int}$-graded, and
$\mathscr{K}_{-}(\mathcal{F}(\mathfrak{g})_{A})\mathscr{K}_{-}(\mathcal{F}(\mathfrak{g})_{B})=0$
if $|A|$ and $|B|$ are odd. ∎

#### 9.7.5.

Corollary.  If $|A|$ is even, then $t(\lambda)\in\mathbb{R}^{+}$, and thus the
isomorphism $\psi_{A}$ admits an integral structure. In particular we have an
isomorphism of graded rings:

$\bigoplus\limits_{\overline{A}=\overline{0}}\mathscr{K}_{-}(\mathcal{F}(\mathfrak{g})_{A})\to\bigoplus\limits_{\overline{A}=\overline{0}}\left(\bigwedge
V^{\mathbb{Z}}\right)_{A}^{spoiled}.$

Here $V^{\mathbb{Z}}$ denotes the free $\mathbb{Z}$-module with basis
$\\{v_{i}:i\in\mathbb{Z}\setminus\\{0\\}\\}$.

## References

  * [1] A. Aramova, L.L. Avramov, J. Herzog, Resolutions of monomial ideals and cohomology over exterior algebras, Trans. of the Amer. Math. Soc., Volume 352-2, (1999), 579–594.
  * [2] I. N. Bernstein, I. M. Gelfand, S. I. Gelfand, On certain category of $\mathfrak{g}$-modules, Functional Analysis and Applications, 10 (1976), no. 2, 1–-8.
  * [3] J. Brundan, Kazhdan-Lusztig polynomials and character formulae for Lie superalgebra $\mathfrak{q}_{n}$, Adv. Math. 182 (2004), 28–77.
  * [4] S. J. Cheng Supercharacters of queer Lie superalgebras, J. Math. Phys. 58, 061701 (2017).
  * [5] J.-S. Cheng, J.-H. Kwon, Finite-dimensional half-integer weight modules over queer Lie superalgebras Comm. Math. Phys., 346 (2016), 945–965.
  * [6] V. Deodhar, O. Gabber, V. Kac, Structure of some categories of representations of infinite-dimensional Lie algebras, Adv. in Math. 45 (1982), 92–116.
  * [7] M. Duflo, V. Serganova, On associated variety for Lie superalgebras, arXiv:0507198.
  * [EAS] I. Entova-Aizenbud, V. Serganova. Jacobson-Morozov lemma for algebraic supergroups, Advances in Mathematics 398, 108240, (2022).
  * [8] A. Frisk, Typical blocks of the category $\mathcal{O}$ for the queer Lie superalgebra, J. Algebra and Applications, 6, No. 5 (2006).
  * [9] M. Gorelik, Shapovalov determinants for $Q$-type Lie superalgebras, IMRP (2006), Art. Id. 96895, 1–71.
  * [10] M. Gorelik, On the ghost centre of Lie superalgebras, Annales de l’Institut Fourier, Tome 50 (2000) no. 6, pp. 1745–1764.
  * [11] M. Gorelik, Depths and cores in the light of DS-functors, arXiv:2010.05721.
  * [12] M. Gorelik, Bipartite extension graphs and the Dulfo–Serganova funtor, arXiv: 2010.12817.
  * [13] M. Gorelik, On modified extension graphs of a fixed atypicality, preprint.
  * [14] M. Gorelik, On Duflo-Serganova functor for the queer Lie superalgebras, preprint.
  * [15] M. Gorelik, T. Heidersdorf, Gruson-Serganova character formulas and the Gruson-Serganova cohomology functor, arXiv:2104.12634.
  * [16] M. Gorelik, A. Sherman, On the Duflo-Serganova functor for the queer Lie superalgebra, arXiv:2204.05048.
  * [17] N.I. Ivanova, A.L. Onishchik, Parabolic subalgebras and gradings of reductive Lie superalgebras, Journal of Mathematical Sciences, Vol. 152, No. 1, 2008.
  * [18] V. G. Kac, D. Kazhdan, Structure of representations with highest weight of infinite-dimensional Lie algebras, Adv. in Math. 34 (1979), no. 1, 97-108.
  * [19] N. Grantcharov, V. Serganova Extension quiver for Lie superalgebra $\mathfrak{q}(3)$, SIGMA Symmetry Integrability Geom. Methods Appl. 16 (2020), Paper No. 141.
  * [20] C. Hoyt, S. Reif, Grothendieck rings for Lie superalgebras and the Duflo-Serganova functor, Algebra Number Theory 12 (2018), no. 9, 2167–-2184.
  * [21] M.S. Im, S. Reif, V. Serganova Grothendieck rings of periplectic Lie superalgebras, arXiv:1906.01948.
  * [22] V. G. Kac, Lie superalgebras, Adv. in Math., 26, no. 1 (1977), 8–96.
  * [23] V. G. Kac, D. H. Peterson, Infinite flag varieties and conjugacy theorems, Proc. Nat. Acad. Sci. U.S.A. 80 (1983), 1778–-1782.
  * [24] V. Mazorchuk, Parabolic category O for classical Lie superalgebras. Advances in Lie superalgebras, 149–166, Springer INdAM Ser., 7, Springer, Cham, 2014.
  * [25] S. Reif, Grothendieck rings of queer Lie superalgebras, arxiv:2107.02219.
  * [26] V. Serganova On the superdimension of an irreducible representation of a basic classical Lie superalgebra, in Supersymmetry in mathematics and physics, 253-–273, Lecture Notes in Math., 2027, Springer, Heidelberg, (2011).
  * [27] V. Serganova, Finite-dimensional representation of algebraic supergroups, ICM
  * [28] V. Serganova, Quasireductive supergroups, Contemporary Mathematics (2011), Vol. 544, 141-159
  * [29] A. Sergeev, The centre of the enveloping algebra for Lie superalgebra $Q(n,\mathbb{C})$, Lett. Math. Phys. 7 (1983), no. 3, 177–179.
  * [30] A. Sergeev, A. Veselov, Grothendieck rings of basic classical Lie superalgebras, Annals of Mathematics, Vol. 173 No. 2 (2011), 663-703.