# Hidden Markov modeling of single particle diffusion with stochastic
tethering

Amit Federbush  Department of Condensed Matter Physics, Tel Aviv University
The Center for Physics and Chemistry of Living Systems, Tel Aviv University
Amit Moscovich  Department of Statistics and Operations Research, Tel Aviv
University  Yohai Bar-Sinai  Department of Condensed Matter Physics, Tel Aviv
University  The Center for Physics and Chemistry of Living Systems, Tel Aviv
University

###### Abstract

The statistics of the diffusive motion of particles often serve as an
experimental proxy for their interaction with the environment. However,
inferring the physical properties from the observed trajectories is
challenging. Inspired by a recent experiment, here we analyze the problem of
particles undergoing two-dimensional Brownian motion with transient tethering
to the surface. We model the problem as a Hidden Markov Model where the
physical position is observed and the tethering state is hidden. We develop an
alternating maximization algorithm to infer the hidden state of the particle
and estimate the physical parameters of the system. The crux of our method is
a saddle-point-like approximation, which involves finding the most likely
sequence of hidden states and estimating the physical parameters from it.
Extensive numerical tests demonstrate that our algorithm reliably finds the
model parameters, and is insensitive to the initial guess. We discuss the
different regimes of physical parameters and the algorithm’s performance in
these regimes. We also provide a free software implementation of our
algorithm.

##  I Introduction

Since the early days of statistical mechanics, the statistics of the
stochastic motion of mesoscopic particles are an important experimental probe
for their microscopic properties. Most prominently, the Gaussian statistics of
Brownian motion provided experimental proof of the atomic nature of matter [1,
2, 3] and was used to measure Avogadro’s number [4]. To this day, new models
are actively developed to explain deviations from purely Brownian statistics
in biological and colloidal systems [5, 6]. Of specific interest are systems
undergoing anomalous diffusion which exhibit a Fickian behavior, i.e. a mean-
squared displacement that is linear with time where the displacement
statistics are non-Gaussian [5, 6, 7]. These deviations from Gaussianity can
serve as an accessible probe for various experiments [8, 9, 10, 11, 12, 13].

In the analysis of such systems, two key sources for non-Gaussianity are
considered. First, particles may undergo several different types of diffusion
modes, stochastically switching between them [14, 15, 16, 17, 18, 19]. Second,
particles may be transiently confined to a small region [20, 21, 22, 23, 24,
25, 9, 26, 27, 28, 29, 30, 31, 32]. In this work, inspired by the experimental
system of Chakraborty et. al. [9], we focus on the latter, and specifically on
two-dimensional (2D) diffusion with transient tethering to the underlying
surface. Fig. 1 depicts an example trajectory from such an experiment. In such
experiments, colloids or nano-particles are coated with molecules of interest
(typically, peptides) and undergo 2D diffusion on a surface coated with a
different molecule. The interaction between the molecules leads to stochastic
transient tethering of the particles to the surface. The experiment aims to
extract information regarding the interaction between the peptides from the
frequency of these tethering and untethering events.

However, identifying these events may be challenging, since the tethering to
the surface is not directly observed, but rather needs to be inferred from the
observed trajectories. Several computational methods have been suggested to
tackle this in the past [22, 24, 31]. In this work, we present a simple and
computationally efficient algorithm to infer tethering and untethering events
from observed trajectories, to estimate the tethering and untethering rates,
the diffusion coefficient, and the effective confinement area of the
interaction potential. The algorithm is developed under the assumption that
the dynamics follow normal diffusion with Poissonian tethering/untethering
events, but can be easily extended to other types of diffusion and
tethering/untethering rates.

The structure of this manuscript is as follows: In Section II we describe the
stochastic hidden Markov model which we use to model the problem, discuss the
different time scales that the model introduces, and define the regimes of
model parameters where our methods are valid. In Section III we describe our
algorithm and its underlying approximations, and in Section IV we present our
numerical results on synthetic data. Finally, in Section V we discuss the
results and possible generalizations. The full source code of the implemented
algorithm is available on GitHub [33].

Figure 1: A snapshot of an experimental microscopy video of peptide-coated
microparticles, in a system similar to that of Chakraborty et. al. [9]. The
colored curve represents the top-right particle’s recorded trajectory in the
320 seconds preceding the snapshot. The trajectory is colored based on time,
as indicated by the timeline at the bottom. During the time segment
$t\in[68,294]$ the particle is confined to a small region of space (shown in
green) due to tethering. Image courtesy of Amandeep Sekhon, Roy Beck, and Yael
Roichman.

##  II Stochastic Model

Our Markov model describes a particle alternating between a freely diffusing
state and a tethered state in 2D. The transition between the states is modeled
as a standard two-state continuous-time Markov chain with characteristic times
$\tau_{0}$ and $\tau_{1}$, which are the inverses of the average rates of
tethering and untethering events, respectively. We denote the state of the
particle by $S(t)$ where $S=0$ corresponds to the free state and $S=1$ to the
tethered state. In the free state, the particle undergoes standard Brownian
motion, while in the tethered state, it is also confined to a harmonic
potential centered at the tether point $X^{*}(t)$. Explicitly, the position of
the particle $X(t)$ follows an Ornstein-Uhlenbeck process with a state-
dependent potential:

$\displaystyle\dot{X}(t)$
$\displaystyle=-\frac{k}{\gamma}(X(t)-X^{*}(t))S(t)+\sqrt{2D}\xi(t),$ (1)

where $k$ is the spring constant, $\gamma$ is the friction coefficient,
$S(t)\in\\{0,1\\}$ indicates the free/tethered state, $D$ is the diffusion
coefficient, and $\xi$ is a standard two-dimensional white noise
$\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$.
Observe that only the ratio $k/\gamma$ plays a role in the model, rather than
$k$ and $\gamma$ individually. We assume the particle is tethered to the point
at which the transition $S=0\to 1$ occurs. That is, if the particle becomes
tethered at time $t_{1}$, and remains tethered until untethering at $t_{2}$,
then $X^{*}(t)$ satisfies

$\displaystyle X^{*}(t)=X(t_{1})\qquad\forall t\in[t_{1},\ t_{2}).$ (2)

To summarize, the particle dynamics are modeled by a continuous-time Markov
process,

$\displaystyle F(t):=\big{(}X(t),\ S(t),\ X^{*}(t)\big{)}.$ (3)

The model is specified by 4 parameters: $\tau_{0},\tau_{1},D$ and the ratio
$k/\gamma$. It is more convenient to work with the variable $A=(D\gamma)/k$
which is the characteristic area that the particle explores in the tethered
state. For notational purposes, we group the model parameters as
$\Theta:=(\tau_{0},\tau_{1},D,A)$. Lastly, we note that
$\frac{A}{D}=\frac{\gamma}{k}$ has units of time. Its meaning is discussed
below.

The transition probabilities of $S$ are Poissonian. In the untethered state
$S(t)=0$, the particle undergoes classic Brownian motion. Hence its position
at $t+\Delta t$ is normally distributed around $X(t)$,

$\displaystyle P\left(X(t+\Delta t)\right)=\frac{1}{4\pi D\Delta
t}\exp\left(-\tfrac{\left(X(t+\Delta t)-X(t)\right)^{2}}{4D\Delta t}\right).$
(4)

In the tethered state $S(t)=1$, the probability density function of
$X(t+\Delta t)$ is given by the solution of the Fokker-Planck equation with
strong friction in a harmonic potential around the anchor point $X^{*}(t)$
[34],

$\displaystyle P\left(X(t+\Delta t)\right)=$ (5) $\displaystyle\frac{1}{2\pi
A^{\prime}(\Delta t)}\exp\left({-\tfrac{\big{(}X(t+\Delta t)-\phi(\Delta
t)X(t)-(1-\phi(\Delta t))X^{*}(t)\big{)}^{2}}{2A^{\prime}(\Delta t)}}\right),$

where we defined two auxiliary variables:

$\displaystyle\phi(\Delta t)$ $\displaystyle=\exp{\left(-\frac{D\Delta
t}{A}\right)},$ $\displaystyle A^{\prime}(\Delta t)$
$\displaystyle=(1-\phi^{2}(\Delta t))A.$ (6)

###  II.1 Discretized dynamics

The model defined above describes the full continuous-time dynamics. In
principle, we could use it to estimate the model parameters
$\Theta=(\tau_{0},\tau_{1},D,A)$. However, the experimental setup poses two
difficulties: First, we can only measure the particle positions $X(t)$ and do
not have access to the particle states $S(t)$ and tether points $X^{*}(t)$.
This is the core challenge of the problem which we address in Section III.
Second, we only sample the process at discrete times
$t_{1},t_{2},\dots,t_{N}$, separated by a finite time resolution $\Delta
t=t_{i}-t_{i-1}$. To account for this, we define a discrete Markov process,
analogous to the continuous process, in which the tether point is constrained
to be one of the previously observed positions. This discrete Markov process
defines a sequence of random states $\\{F_{n}\\}_{n=1}^{N}$ where

$\displaystyle F_{n}:=F(t_{n})=\big{(}X(t_{n}),\ S(t_{n}),\
X^{*}(t_{n})\big{)}.$ (7)

The transition probabilities of $F_{n}$ are given by the product of the
transition probabilities of $S$ and $X$:

$\displaystyle P(F_{n+1}|F_{n})$
$\displaystyle=P(S_{n+1}|S_{n};\Theta)P(X_{n+1}|F_{n};\Theta),$ (8)

where $X^{*}_{n+1}$ is uniquely determined by the history of $\\{X_{n}\\}$ and
$\\{S_{n}\\}$, in accordance with Eq. 2:

$\displaystyle X^{*}_{n+1}=\begin{cases}X_{n+1}&S_{n+1}=1\mbox{ and
}S_{n}=0,\\\ X_{n}^{*}&\mbox{otherwise}.\end{cases}$ (9)

The state transitions are Poissonian and to leading order in $\Delta t$ read

$\displaystyle P$
$\displaystyle\left(S_{n+1}|S_{n}\right)=\begin{cases}1-\frac{\Delta
t}{\tau_{n}}&S_{n+1}=S_{n}\\\ \frac{\Delta t}{\tau_{n}}&S_{n+1}\neq
S_{n},\end{cases}$ (10)

where $\tau_{n}=\tau_{0}$ if $S_{n}=0$ and $\tau_{n}=\tau_{1}$ if $S_{n}=1$.
The distribution of the particle position in the next step follows Eq. (4) and
(5),

$\displaystyle P(X_{n+1}$ $\displaystyle|F_{n};\Theta)$ (11)
$\displaystyle=\begin{cases}\frac{1}{4\pi D\Delta t}e^{-\frac{1}{4D\Delta
t}\Big{(}X_{n+1}-X_{n}\Big{)}^{2}}&S_{n}=0,\\\ \frac{1}{2\pi
A^{\prime}}e^{-\frac{1}{2A^{\prime}}\Big{(}X_{n+1}-\phi
X_{n}-(1-\phi)X^{*}_{n}\Big{)}^{2}}&S_{n}=1.\end{cases}$

where now $\phi$ and $A^{\prime}$ are constants that depend on $A,D$ and the
time step $\Delta t$ as in Eq. 6. Furthermore, we assume that at the beginning
of the measurement the particle is either free ($S_{1}=0$) or tethered at
$X_{1}$ ($S_{1}=1,X^{*}_{1}=X_{1}$). Importantly, $X^{*}_{n}$ can only attain
values among the previously visited positions $X_{1},\dots,X_{n}$.

As usual for Markov processes, the probability of a trajectory
$\\{F_{n}\\}_{n=1}^{N}$ is the product of the transition probabilities in each
step, thus the log-probability is additive:

$\displaystyle\begin{split}\log&{P\left(\\{F_{n}\\}_{n=1}^{N};\Theta\right)}\\\
&=\log P(F_{1})+\sum_{n=1}^{N-1}{\log{P(F_{n+1}|F_{n};\Theta)}}.\end{split}$
(12)

An example trajectory generated from the discrete model is depicted in Fig. 2.

Figure 2: An example trajectory of $N=2000$ steps generated from the discrete
model with $\tau_{0}=\tau_{1}=100$, $D=1$ and $A=0.5$. The sampling time step
is $\Delta t=1$. The trajectory is color-coded according to the particle’s
state, where blue and orange depict the free and tethered states,
respectively.

###  II.2 Time scales

Figure 3: A section of the phase space of $\Delta t,\tau_{0},\tau_{1}$, in
time units of $\frac{A}{D}$. Note the logarithmic scale. The yellow region is
where our assumptions hold. The bottom horizontal boundary of the yellow
region corresponds to the inequality $\frac{A}{D}\lesssim\Delta t$ and the
diagonal boundary corresponds to the inequality $\Delta
t\ll\tau_{0},\tau_{1}$. The bottom and left regions are annotated according to
the discussion in this subsection. The markers correspond to the seven regimes
to be analyzed in Section IV.

As seen above, the problem involves multiple time scales. Three of them are
given by the model parameters and represent the underlying physics:
$\tau_{0},\tau_{1}$, and $\frac{A}{D}$ which is the equilibration time of the
Brownian particle with the harmonic potential. The other two time scales,
$\Delta t$ and the total experiment duration $T=t_{N}-t_{1}$ are properties of
the experiment. We employ several realistic working assumptions about these
time scales that greatly reduce computational complexity. First, we assume
that $\tau_{0},\tau_{1}\gg\frac{A}{D}$, which physically means that the rate
of tethering/untethering events is much smaller than the inverse of the
harmonic equilibration time of the particle. Violation of this condition means
that the particle can untether before the tethering potential has a
significant effect and thus tethering events will not be experimentally
discernible. This regime corresponds to standard Brownian motion with an
effective diffusion constant smaller than $D$. Second, we require that $\Delta
t\ll\tau_{0},\tau_{1}$, to avoid the possibility of multiple transitions of
tethering or untethering events within a single sampling interval. This
assumption is experimentally realistic since modern cameras can easily achieve
frame rates larger than $10^{3}$Hz, and in many experiments, $\tau_{0}$ and
$\tau_{1}$ are of orders of at least seconds [9, 20, 25]. In any case, if
$\Delta t\approx\tau_{0},\tau_{1}$, this would again correspond to an
essentially pure diffusive behavior in the discretized data. Third, if $\Delta
t<\frac{A}{D}$, the time discretization resolves the equilibration of a
tethered particle with its confining potential. Since our goal is only to
extract the physical parameters of the system, such resolution does not add
relevant information and only increases the computational cost. Therefore,
under-sampling the discrete dynamics to increase $\Delta t$ should not lead to
a significant loss of accuracy in the estimation of $\tau_{0},\tau_{1}$ but
would greatly reduce the search space. This is explicitly demonstrated in
Section IV. Even if the trajectories are experimentally measured with small
$\Delta t$, we can safely downsample them such that $\frac{A}{D}\lesssim\Delta
t$, or in other words $\phi\lesssim e^{-1}$, cf. Eq. 6. For concreteness, we
mention the experimental parameters of the diffusing nanoparticles system of
Chakraborty et. al. [9]. In this system, $\tau_{0},\tau_{1}\sim 1\text{s}$,
$D\sim 10\mu$m2/s, $A\sim 1\mu\text{m}^{2}$, which means a sampling interval
of $\Delta t=A/D=0.1$s is sufficient and is easily achievable experimentally.
To conclude this discussion, we assume the following separation of time
scales:

$\frac{A}{D}\lesssim\Delta t\ll\tau_{0},\tau_{1},$ (13)

The leading-order expansion of Eq. 11, which is first order in $\Delta t$ and
zeroth order in $\phi$, now reads

$\displaystyle P$
$\displaystyle\left(S_{n+1}|S_{n}\right)=\begin{cases}1-\frac{\Delta
t}{\tau_{n}}&S_{n+1}=S_{n}\\\ \frac{\Delta t}{\tau_{n}}&S_{n+1}\neq
S_{n}\end{cases}$ (14) $\displaystyle P$
$\displaystyle\left(X_{n+1}|F_{n}\right)=\begin{cases}\frac{1}{4\pi D\Delta
t}e^{-\frac{1}{4D\Delta t}\Big{(}X_{n+1}-X_{n}\Big{)}^{2}}&S_{n}=0\\\
\frac{1}{2\pi
A}e^{-\frac{1}{2A}\Big{(}X_{n+1}-X^{*}_{n}\Big{)}^{2}}&S_{n}=1,\end{cases}$

where $\tau_{n}=\tau_{0}$ if $S_{n}=0$ and $\tau_{n}=\tau_{1}$ if $S_{n}=1$.

To conclude the discussion regarding time scales, Fig. 3 illustrates the
section of phase space where our assumptions hold (in time units of
$\frac{A}{D}$). For a point in the yellow region in the figure, the farther it
is from the two boundary lines, the better is the separation of scales and the
better our assumptions hold. The marked points correspond to regimes that will
be analyzed in Section IV.

##  III Our method

Our problem is as follows: in an experiment, we can measure the observed
states $\\{X_{n}\\}$, but we do not have access to the hidden states
$\\{S_{n},X^{*}_{n}\\}$, nor to the model parameters
$\Theta=(\tau_{0},\tau_{1},D,A)$. We use the term hidden path to denote the
sequence of the hidden states $\\{S_{n}\\}_{n=1}^{N}$ from which the sequence
$\\{X^{*}_{n}\\}_{n=1}^{N}$ can be determined. The goal is to infer the model
parameters $\Theta$ from a series of measurements. To this end, we developed
an alternating maximization algorithm [35], similar to the classical
expectation-maximization (EM) algorithm [36], to estimate both the hidden
states and the parameters of our hidden Markov model.

If the hidden path is known, the problem of optimal parameter estimation is
fairly standard and is typically solved by maximizing the likelihood of the
model parameters, which according to Bayes’ rule is proportional to the
exponential of Eq. 12. Maximizing $P\left(\Theta|\\{F_{n}\\}\right)$ can be
done numerically or using analytical approximations, as described below.
However, when the hidden path is not known, the likelihood function to
consider is

$\displaystyle\mathcal{L}\left(\Theta|\\{X_{n}\\}\right)$
$\displaystyle=\frac{P(\Theta)}{P(\\{X_{n}\\})}P\left(\\{X_{n}\\}\right|\Theta)$
(Bayes)
$\displaystyle=\frac{P(\Theta)}{P(\\{X_{n}\\})}\sum_{\\{S_{n}\\}}P\left(\\{F_{n}\\}\right|\Theta),$
(15)

where the sum is over all $2^{N}$ possibilities for the hidden paths,
$P(\Theta)$ is the prior probability distribution for the model parameters,
and the evidence $P(\\{X_{n}\\})$ is a constant we may ignore [37]. Recall
that $\\{F_{n}\\}=\\{X_{n},S_{n},X^{*}_{n}\\}$ is the trajectory of both the
observed and hidden states of the particle, and the log-probability of such a
trajectory is the sum of the log-probabilities of all the steps, as in Eq. 12.
Computing the sum in Eq. 15 is intractable for typical values of $N$. However,
numerical evidence shows that most hidden paths are very unlikely and thus
have a negligible contribution to this sum. Taking a uniform prior,
$P(\Theta)={\rm const}$, we posit that

$\log{\mathcal{L}(\Theta|\\{X_{n}\\})}\sim\max_{\\{S_{n}\\}}{\log
P(\\{F_{n}\\}|\Theta)}.$ (16)

A key idea of our method is to use the RHS of Eq. (16) as a proxy for the
computationally intractable LHS. It is analogous to the saddle-point
approximation from statistical mechanics, where the integral is replaced with
the maximum of the integrand. To compute the RHS of Eq. (16), we need only to
find the most likely hidden path $\\{\hat{S}_{n}\\}_{n=1}^{N}$, given the
model parameters $\Theta$. This discrete optimization problem can be
efficiently solved using the Viterbi algorithm from dynamic programming [38,
37]. Below we briefly describe the parameter estimation method and how the
Viterbi algorithm is implemented for our model. Then we present our
alternating maximization approach for estimating the maximum likelihood model
parameters.

###  III.1 Parameter estimation

Given the most likely hidden path $\\{\hat{S}_{n}\\}$, maximizing the log-
likelihood in Eq. 16 is maximizing a sum of the log terms from Eq. 14. We use
the following maximum-likelihood estimators (MLE) to approximate the most
likely model parameters:

$\displaystyle\hat{\tau}_{0}$
$\displaystyle:=\frac{N_{00}+N_{01}}{N_{01}}\Delta
t\qquad\hat{\tau}_{1}:=\frac{N_{11}+N_{10}}{N_{10}}\Delta t$
$\displaystyle\hat{D}$ $\displaystyle:=\frac{1}{4(N_{00}+N_{01})\Delta
t}\sum_{n=1}^{N-1}{(1-\hat{S}_{n})(X_{n+1}-X_{n})^{2}}$ $\displaystyle\hat{A}$
$\displaystyle:=\frac{1}{2(N_{10}+N_{11})}\sum_{n=1}^{N-1}{\hat{S}_{n}(X_{n+1}-X^{*}_{n})^{2}},$
(17)

where $N_{ij}$ is the number of $i\to j$ transitions of $\hat{S}_{n}$ that
occur along the trajectory. The derivation of these estimators is trivial by
taking the derivative of Eq. 16 with respect to each model parameter and
equating to zero.

###  III.2 Finding the most likely sequence of states

We now describe how to find the maximum likelihood path $\\{\hat{S}_{n}\\}$,
conditioned on the model parameters $\Theta$. Recall that the tether point
$X^{*}_{n}$ can only assume one of the previously visited positions
$X_{1},X_{2},...,X_{n}$, see Eq. 9. We represent the set of all possible paths
as a directed layer graph (Fig. 4), known as a _trellis_ in the Viterbi
literature. In this graph, columns correspond to the discrete time $n$ and
rows to the tether point $X^{*}_{n}$ which can be either one of
$X_{1},\ldots,X_{n}$ or _free_ (untethered). Each vertex represents the state
of the particle at time $n$ and each path from left to right corresponds to a
specific sequence of tethered/free states.

We set the edge weights to be the logarithm of the transition probabilities
between states according to Eq. 8. With this choice of edge weights, the log-
likelihood of a path (Eq. 12) is just $\log P(F_{1})$ plus the sum of edge
weights along the corresponding path in the graph. Thus, finding the maximum-
likelihood sequence of hidden states is reduced to the problem of finding the
path of maximum weight in a directed layer graph. The latter problem is
efficiently solved using the Viterbi algorithm [38]. This algorithm scans the
trellis column by column from left to right and computes, for each vertex, the
maximum-weight path that ends at that vertex.

Figure 4: A trellis graph of the observed and hidden states, given the
observed particle positions $X_{1},...,X_{4}$. Each node represents a
different state $F_{n}=(X_{n},S_{n},X^{*}_{n})$, with the row corresponding to
the hidden state, and the column to the time $n$ and the observed state
$X_{n}$. The edges represent allowed state transitions and are weighted
according to the logarithm of Eq. 8. Each trajectory $\\{F_{n}\\}$ is given as
a path on the graph that advances from left to right.

###  III.3 The alternating maximization algorithm

Given an estimate of the hidden path $\\{S_{n}\\}$, we can apply Eq. 17 to
obtain the maximum likelihood estimate of the model parameters
$\Theta=(\tau_{0},\tau_{1},D,A)$. Conversely, given an estimate of the model
parameters $\Theta$, it is easy to find the maximum likelihood hidden path
using the Viterbi algorithm as explained in the previous subsection. Combining
these two observations naturally leads to an alternating maximization
procedure for estimating both the most likely model parameters $\hat{\Theta}$
and the most likely hidden path:

  1. 1.

Initial guess: Guess an initial value $\hat{\Theta}^{(0)}$ for the four model
parameters.

  2. 2.

Path maximization step: Conditioned on the current parameter estimate
$\hat{\Theta}^{(m)}$, apply the Viterbi algorithm to find
$\\{S_{n}^{(m)}\\}_{n=1}^{N}$, the most likely hidden path given
$\hat{\Theta}^{(m)}$,

$\displaystyle
S^{(m)}:=\operatorname*{argmax}_{\\{S_{n}\\}}\mathcal{L}\big{(}\\{S_{n}\\}|\hat{\Theta}^{(m)},\\{X_{n}\\}\big{)}.$
(18)

  3. 3.

Parameter maximization step: Use Eq. 17 to obtain the maximum likelihood
estimate of the model parameters conditioned on the current estimate of the
hidden path,

$\displaystyle\hat{\Theta}^{(m+1)}:=\hat{\Theta}\big{(}\\{S^{(m)}\\},\\{X_{n}\\}\big{)}.$
(19)

  4. 4.

Convergence: Alternate between steps 2 and 3 until
$\left|\hat{\theta}^{(m+1)}-\hat{\theta}^{(m)}\right|/\hat{\theta}^{(m)}\leq\epsilon$
for each $\hat{\theta}\in\hat{\Theta}$.

Since the sample space is discrete, convergence typically occurs exactly, i.e.
$\hat{\Theta}^{(n+1)}=\hat{\Theta}^{(n)}$. However, in our experiments we used
$\epsilon=10^{-3}$ to stop the iterations when the relative change is small.
The specific threshold $10^{-3}$ is inconsequential. Furthermore, we set the
maximum number of iterations to 20, to prevent the possibility of infinite
loops that alternate between several discrete hidden paths without satisfying
the convergence criterion (in practice, less than 10 iterations typically
suffice for convergence). We also define a criterion for divergence: if after
the parameter maximization step the estimators $\hat{\tau}_{0}$ or
$\hat{\tau}_{1}$ exceed the arbitrary threshold of $0.9T$, we say the
algorithm diverged and stop the iterations. This is done to avoid cases where
the particle is estimated to stay tethered or untethered for the entire
sampling time due to an unphysical divergence of the MLE estimators. In
practice, we test the regimes where $T\gg\tau_{0},\tau_{1}$ and most runs do
not diverge. We say the algorithm has converged if the convergence criterion
was satisfied without first triggering the divergence criterion.

The Viterbi algorithm scans each edge on the graph once, leading to an
$O(N^{2})$ complexity. However most paths are not likely, e.g. a path where at
some point $S_{n}=1$ and $(X_{n}-X^{*}_{n})^{2}\gg A$, meaning the particle
roams far from its tether point. To improve the running time, we implemented
edge pruning of the trellis graph. Specifically, we discard the paths of all
but the $q$ most likely tethered nodes at each time step (column) during the
algorithm’s execution. This reduces the complexity to $O(qN)$. In our tests,
we used $q=10$, which showed no significant impact on the results compared to
no pruning.

In Appendix A we discuss an extension of our algorithm that outputs the $K$
most likely hidden paths rather than the single best path. We generally found
that the $100$ most likely paths yield similar parameter estimates, so the
results in the rest of this paper use $K=1$.

###  III.4  Bootstrap bias correction

For reasons that will be discussed in Section IV, the MLE predictions for the
model parameters give very good estimates of $D$ and $A$, but consistently
overestimate the waiting times $\tau_{0}$ and $\tau_{1}$. This overestimation
can be largely corrected by using a parametric bootstrap procedure [39]. The
idea is simple: define the bias $\Delta$ as the expected difference between
the MLE and the true parameters:
$\Delta=\left\langle\hat{\Theta}\right\rangle-\Theta$. Assuming the bias does
not depend sensitively on $\Theta$ itself, one can estimate $\Delta$ by
simulating trajectories with parameters $\hat{\Theta}$ and rerunning the
procedure on this simulated data (for which we know the true parameters). This
is done as follows:

  1. 1.

Apply the alternating maximization algorithm described in Section III.3 to
obtain the MLE $\hat{\Theta}$.

  2. 2.

Generate $M$ random trajectories using $\hat{\Theta}$ as the “true” simulated
model parameters.

  3. 3.

Apply the alternating maximization algorithm to each simulated trajectory,
yielding the estimated model parameters $\hat{\Theta}^{\prime}_{i}$ for $1\leq
i\leq M$.

  4. 4.

The bias is estimated by the median value of the observed biases,

$\hat{\Delta}:=\operatornamewithlimits{median}_{1\leq i\leq
M}(\hat{\Theta}^{\prime}_{i}-\hat{\Theta})\ .$

  5. 5.

Finally, the bias is used to correct the original estimate. That is, the bias-
corrected estimate for the true model parameters is given by

$\displaystyle\hat{\Theta}^{B}=\hat{\Theta}-\hat{\Delta}.$ (20)

Note that the median is calculated for each parameter $\theta\in\Theta$
independently. We chose to estimate $\hat{\Delta}$ using the median since it
is robust to outliers and we found that it gave better results in mean squared
error.

##  IV  Results

We test the algorithm’s performance on synthetically generated trajectories,
whose model parameters $\Theta$ are known. The implementation of the algorithm
in Python, as well as the code to generate all figures in this manuscript, is
available at [33].

In all trajectories, we use $D=1$ and $A=1$ and a total experiment duration of
$T=10000$. Note that both $D$ and $A$ can always be set to unity by properly
choosing length and time units, so this choice is without loss of generality.
Thus, each experiment is specified by three scalars: the sampling time $\Delta
t$ and the model parameters $\tau_{0}$,$\tau_{1}$, all expressed in time units
of $\frac{A}{D}=1$. We focus on seven different regimes of $\Theta$ and
$\Delta t$, detailed in Table 1. We mention that the conditions in Eq. 13 are
barely satisfied in regimes $3$ and $5$. This is intentional as we wish to
test the algorithm slightly beyond its bounds.

Regime |  $D,A$ |  $\Delta t$ |  $\tau_{0}$ |  $\tau_{1}$ |  Acc. (%) |  $\hat{\tau}_{0}$ |  $\hat{\tau}_{1}$ |  $\hat{D}$ |  $\hat{A}$
---|---|---|---|---|---|---|---|---|---
1 |  $1$ |  $10$ |  $100$ |  $100$ |  $96\pm 2$ |  $131\pm 24$ |  $130\pm 19$ |  $1.00\pm 0.05$ |  $0.99\pm 0.05$
2 |  $1$ |  $1$ |  $100$ |  $100$ |  $94\pm 2$ |  $122\pm 21$ |  $122\pm 16$ |  $1.00\pm 0.01$ |  $0.99\pm 0.02$
3 |  $1$ |  $0.5$ |  $100$ |  $100$ |  $88\pm 4$ |  $125\pm 22$ |  $123\pm 17$ |  $1.00\pm 0.01$ |  $0.99\pm 0.02$
4 |  $1$ |  $10$ |  $50$ |  $50$ |  $93\pm 2$ |  $77\pm 11$ |  $75\pm 8$ |  $0.99\pm 0.05$ |  $0.98\pm 0.05$
5 |  $1$ |  $10$ |  $20$ |  $20$ |  $87\pm 2$ |  $47\pm 9$ |  $43\pm 5$ |  $0.97\pm 0.06$ |  $0.94\pm 0.06$
6 |  $1$ |  $10$ |  $200$ |  $50$ |  $96\pm 1$ |  $356\pm 95$ |  $79\pm 12$ |  $0.99\pm 0.04$ |  $0.97\pm 0.09$
7 |  $1$ |  $10$ |  $50$ |  $200$ |  $97\pm 2$ |  $60\pm 11$ |  $248\pm 43$ |  $1.01\pm 0.08$ |  $1.00\pm 0.04$
Table 1: True model parameters, accuracy, and estimated parameters (mean $\pm$
standard deviation) for the 7 regimes analyzed using the alternating
maximization algorithm (Section III.3) without bias-correction. The accuracy
is defined as the fraction of time steps at which the algorithm correctly
predicted both the tethered state $S$ and the tether point $X^{*}$.

###  IV.1 Stability

As described above, the algorithm requires an initial guess of the model
parameters $\hat{\Theta}^{(0)}$. We found empirically that the algorithm
converges to a single fixed point almost regardless of the initial parameters.
To demonstrate this, we show in Fig. 5 the results of the algorithm when
applied to a single observed trajectory corresponding to regime $1$, with 1000
different $\hat{\Theta}^{(0)}$’s, randomly drawn from a log-uniform
distribution spanning two decades around $\Theta$. We mention this
experiment’s wide distribution of initializations $\hat{\Theta}^{(0)}$ is
intended to illustrate the algorithm’s stability, and in a practical
situation, better priors can often be used, especially for $D$ and $A$. The
convergence for other regimes is qualitatively similar and is not shown here.
Out of the 1000 algorithm runs, 960 converged and only these runs were taken
into account in the figure. Since the outcome of the algorithm is largely
independent of the initial guess, in what follows we use the true model
parameters $\Theta$ as the initial guess $\hat{\Theta}^{(0)}$.

Figure 5: The algorithm is insensitive to the initialization $\Theta^{(0)}$:
Here we show the evolution of the estimators
$\hat{\tau}_{0},\hat{\tau}_{1},\hat{D},\hat{A}$ of the four model parameters
along the algorithm iterations over a single synthetic trajectory. The values
are normalized by the true parameter values so that a perfect prediction would
be unity in all panels. The first iteration is the initial parameter guess. We
drew 1000 initial guesses at random from a log-uniform distribution around the
true value and ran the algorithm, keeping only the 960 runs that converged.
The evolution of the estimators of each (converged) run is shown, exhibiting a
strong convergence towards a single fixed point. Note that the vertical axis
is logarithmic. The dashed red line represents the most likely parameter
value, given the true hidden path. It does not coincide, in general, with the
true model parameters.

###  IV.2 Recovering the parameters and hidden states

For each of the regimes in Table 1, we generated 1000 synthetic trajectories
and ran the algorithm with initial parameters that are equal to the true model
parameters, as discussed in the preceding subsection regarding stability. Over
$98\%$ of the runs in each regime have converged (the other runs entered a
loop and stopped after the maximum number of iterations). Runs that converged
typically did so within 3-8 iterations. We define the accuracy as the fraction
of time steps for which the algorithm predicted both the correct state $S_{n}$
and tether point $X^{*}_{n}$ (if $S_{n}=1$). The mean and standard deviation
of the accuracy over all converged runs are detailed in Table 1. It is seen
that the accuracy is fairly high in all regimes, and that it decreases as
$\Delta t$ approaches $\frac{A}{D}$ from above or $\tau_{0}$,$\tau_{1}$ from
below, consistently with the constraints of Eq. 13.

Next, we examine how well the algorithm recovers the model parameters. For
each regime, Table 1 depicts the mean and the standard deviation of the
algorithm’s estimates
$\hat{\Theta}=(\hat{\tau}_{0},\hat{\tau}_{1},\hat{D},\hat{A})$ over all
converged runs. These results exhibit two clear trends: the temporal
parameters $\tau_{0},\tau_{1}$ are consistently overestimated, while the
spatial ones $D,A$ are correctly estimated (although there is a very slight
yet consistent underestimate of $A$ in most regimes). The overestimate of the
temporal parameters ranges from $20\%$ to $120\%$ and will be discussed in
Section IV.3.

To investigate the effect of the sampling time $\Delta t$, we focus on the
three regimes $1,2,3$ which differ only by $\Delta t$ (see Table 1). The
distributions of the estimated model parameters are depicted in the top row of
Fig. 6. Note that the centers of the distributions hardly change with $\Delta
t$. This supports the claim described in Section II that the sampling time can
be increased without losing information, as long as Eq. 13 holds. We mention
that the width of the distribution for the spatial parameters $D,A$ decreases
with $\Delta t$. That is to be expected since the accuracy of the MLEs
increases with $N$, but for a fixed $T$, the number of samples $N$ is
inversely proportional to $\Delta t$.

Figure 6: Distribution of the algorithm’s estimates of physical model
parameters for regimes 1,2,3 with $\tau_{0}=\tau_{1}=100$ and $\Delta
t=10,1,0.5$, respectively. The dashed lines represent the true parameter
value. Top row: the maximum likelihood estimates of the alternating
maximization algorithm. Bottom row: the estimators after the bootstrap bias
correction procedure.

By focusing on regimes $1,4,$ and $5$, the effect of $\tau_{0},\tau_{1}$, when
they are equal, can be isolated. As $\tau_{0}$ and $\tau_{1}$ decrease, the
accuracy decreases and the relative overestimate of $\tau_{0},\tau_{1}$
increases. This is consistent with Eq. 13, since decreasing
$\tau_{0},\tau_{1}$ challenges the assumption that $\Delta
t\ll\tau_{0},\tau_{1}$. Finally, regimes $6$ and $7$ explore the case of
$\tau_{0}\neq\tau_{1}$. Note that when $\tau_{0}<\tau_{1}$, i.e. when the
particle is more likely to be tethered than free, the parameter estimates are
better, with less overestimation for the temporal parameters, compared to the
opposite case.

###  IV.3 Overestimation and bias correction

As seen in Tables 1 and 6, while the spatial parameters $D$ and $A$ are
correctly estimated to high accuracy, the temporal parameters
$\tau_{0},\tau_{1}$ are consistently overestimated, with an increasing
overestimation as $\tau_{0},\tau_{1}$ decrease compared to $\Delta t$. This
bias is due to a systematic misidentification of brief tethered or untethered
intervals. That is, if the particle tethers and then untethers (or vice versa)
over a short time window, the likelihood is dominated by the temporal terms
(regarding $S_{n}$) in Eq. 14 over the spatial terms. As a result, the
likelihood of a trajectory that stays in the same state is higher than that of
a trajectory that switches tethering states twice within this short interval,
even if the true dynamics involved such a switch.

To see this, take a short time interval of $l$ steps, in which the particle
does not move much. Compare two sequences of hidden states, one in which the
particle is tethered in the first step and untethered in the last one, and one
in which it is free throughout the time window (a similar argument holds for
the opposite scenario). Taking the difference between the log-likelihoods of
the two hidden trajectories using Eq. 14, and expanding to leading order in
$l$, we obtain

$\displaystyle\Delta\mathcal{L}=2\log\left(\frac{1-\Delta t/\tau}{\Delta
t/\tau}\right)+\mathcal{O}\left(l,\Delta t^{2}\right)$ (21)

where for simplicity we took $\tau_{0}=\tau_{1}=\tau$. For short time windows,
this tends to be positive since $\Delta t\ll\tau_{0},\tau_{1}$, with the first
term slowly diverging as $\Delta t\to 0$.

In addition, since the distribution of the free and tethered interval
durations times is exponential, it is dominated by such short intervals, which
are systematically misclassified by the MLE estimates. An example of this
phenomenon is illustrated in Fig. 7. Since the MLE estimates for
$\tau_{0},\tau_{1}$ (Eq. 17) are inversely proportional to the number of
switching events, missing brief tethering/untethering intervals leads to over-
estimation of the waiting times. We correct this bias using a parametric
bootstrap procedure, as described in Section III.4. The bottom row of Fig. 6
shows the distribution of the bias-corrected estimates $\hat{\Theta}^{B}$,
that are concentrated within a small region around the true value. Note that
although the mode of the distributions for $\tau_{0},\tau_{1}$ is slightly
biased, the means of the distributions are close to the true values, as
evident in Table 2. Importantly, the bootstrapping successfully mitigates the
bias in the rate estimation without affecting the already accurate estimates
of $D$ and $A$. In our simulations we generated $M=100$ bootstrapping
trajectories (cf. Section III.4) for each analyzed trajectory.

Regime |  $D,A$ |  $\Delta t$ |  $\tau_{0}$ |  $\tau_{1}$ |  $\hat{\tau}_{0}^{B}$ |  $\hat{\tau}_{1}^{B}$ |  $\hat{D}^{B}$ |  $\hat{A}^{B}$
---|---|---|---|---|---|---|---|---
1 |  $1$ |  $10$ |  $100$ |  $100$ |  $102\,(71-141)$ |  $100\,(73-139)$ |  $1.00\,(0.91-1.08)$ |  $1.00\,(0.91-1.08)$
2 |  $1$ |  $1$ |  $100$ |  $100$ |  $98\,(73-129)$ |  $100\,(76-138)$ |  $1.00\,(0.97-1.03)$ |  $1.00\,(0.97-1.03)$
3 |  $1$ |  $0.5$ |  $100$ |  $100$ |  $101\,(75-137)$ |  $101\,(69-133)$ |  $1.00\,(0.99-1.02)$ |  $1.00\,(0.98-1.03)$
4 |  $1$ |  $10$ |  $50$ |  $50$ |  $49\,(38-64)$ |  $49\,(36-61)$ |  $1.00\,(0.92-1.08)$ |  $0.99\,(0.87-1.07)$
5 |  $1$ |  $10$ |  $20$ |  $20$ |  $18\,(9-28)$ |  $19\,(11-26)$ |  $0.98\,(0.88-1.11)$ |  $0.98\,(0.87-1.10)$
6 |  $1$ |  $10$ |  $200$ |  $50$ |  $190\,(100-316)$ |  $51\,(32-77)$ |  $1.00\,(0.92-1.06)$ |  $0.99\,(0.85-1.16)$
7 |  $1$ |  $10$ |  $50$ |  $200$ |  $51\,(36-76)$ |  $198\,(131-295)$ |  $1.00\,(0.87-1.20)$ |  $1.00\,(0.94-1.08)$
Table 2: True model parameters vs. their bootstrap bias-corrected estimates
for the 7 regimes analyzed. The estimates display the mean and a $95\%$
confidence interval. We used $M=100$ synthetic bootstrap trajectories per
analyzed trajectory.

##  V  Discussion and outlook

We developed and tested an efficient and accurate algorithm to analyze
trajectories of diffusing particles with transient tethering, to identify
tethering/untethering events, and to estimate the physical parameters of the
system. The crux of our method is using the Viterbi algorithm to find the most
likely sequence of hidden states, estimate the most likely model parameters
for that sequence, and then bias-correct the estimates using a parametric
bootstrap procedure. Our algorithm successfully recovers the model parameters
and is largely insensitive to the parameter initialization. It is applicable
when the time interval between frames $\Delta t$ is significantly shorter than
the typical tethering/untethering times and longer than, or comparable to, the
equilibration time of the particle with the tethering potential. When the
latter condition is not met, downsampling of the time series can be performed
without significant loss of accuracy in estimating the physical parameters.

Our method as presented above only applies, of course, to normal 2D diffusion
with Poissonian transient tethering. However, it can be readily generalized to
other situations. For example, it can easily be generalized to a different
spatial dimension ($d\neq 2$) by adjusting the prefactors in Eq. 11. In higher
dimensions, the aforementioned misidentification of brief tethered intervals
should be less prominent since a sequence of short steps becomes less likely
(no recurrence). Other diffusive statistics (anomalous diffusion) or non-
Poissonian tethering/untethering can also be modelled by adjusting the
transition probabilities, Eq. 14. For example, it could incorporate a drift
term or an observational error term for position samples, as in Bernstein and
Fricks’s work [22]. However, one must assume a specific functional form for
the diffusion and tethering dynamics – in the regime in which we are working,
only a few dozen tethering/untethering events per trajectory are observed.
Hence, due to the small sample size, statistical goodness-of-fit tests would
lack sufficient power to distinguish, e.g., between Poissonian and non-
Poissonian tethering.

###### Acknowledgements.

We thank Roy Beck, Yael Roichman, Indrani Chakraborty and Amandeep Sekhon for
useful discussions and for providing experimental data. YBS is supported by
ISF grant 1907/22 and by Google Gift grant. AM is supported by ISF grant
1662/22 and NSF-BSF grant 2022778.

##  Appendix A  K most likely paths

Our saddle-point approximation in Eq. 16 replaces the sum over all hidden
paths with the likelihood value of the single most likely path. An immediate
generalization is to sum over the top $K$ likeliest paths, which would improve
the estimation accuracy of $\mathcal{L}$. Given the top $K$ likeliest paths,
we can estimate the model parameters as a weighted sum of the MLEs of each
path, cf. Eq. 17 as

$\hat{\Theta}^{[K]}=\frac{\sum_{i=1}^{K}{P(\\{F_{n}\\}^{[i]}|\Theta)\hat{\Theta}(\\{F_{n}\\}^{[i]})}}{\sum_{i=1}^{K}{P(\\{F_{n}\\}^{[i]}|\Theta)}},$
(22)

where $\\{F_{n}\\}^{[i]}$ is the trajectory corresponding to the $i$th most
likely hidden path. Finding the $K$ highest likelihood paths can be easily
done by a slight modification of the Viterbi algorithm, also known as list
Viterbi [40, 41]. The modification is discarding all but the $K$ most likely
paths ending at each node, rather than all but the most likely path. Doing
this increases the computational complexity by $K^{2}$.

In our numerical tests, we found no significant improvement in the results
when increasing $K$ up to $K=100$. This is because the $K$ most likely paths
are very similar, differing from each other in just a few steps, leading to
very similar MLEs of $\Theta$. Typically, the top paths only differ by slight
perturbations of the tethering and untethering times and do not display
qualitative differences.

To demonstrate this, Fig. 7 depicts the true hidden path of a single
trajectory along with the top $K=10$ most likely hidden paths. The likelihood
of these paths can be computed either with respect to the most likely
parameters given the true hidden path (“oracle”), or the estimated
$\hat{\Theta}$ (without bootstrapping). It is seen that this choice does not
significantly change the result, and the most likely hidden paths are all
qualitatively very similar, differing by just a few steps from each other, and
are similar to the true hidden path. They miss the brief tethered interval at
time $t=93$ and the brief free interval at time $t=220$.

Figure 7: The hidden paths of a single trajectory realization with $N=300$
steps, $\tau_{0}=\tau_{1}=100$, $D=1$, $A=1$, $\Delta t=10$. The true hidden
path is at the top. In the middle are the 10 most likely hidden paths
conditioned on the true parameters, with the upper row being the most likely.
At the bottom are the 10 most likely hidden paths given the algorithm’s
estimated model parameters, again with the upper row being the most likely.
Blue intervals correspond to the free state, and orange intervals correspond
to the tethered state.

Because the results were similar for $K$ up to 100, in the results section we
presented only results pertaining to $K=1$. However, using $K>1$ top paths may
be beneficial in other regimes.

## References

  * Einstein [1905] A. Einstein, Über die von der molekularkinetischen theorie der wärme geforderte bewegung von in ruhenden flüssigkeiten suspendierten teilchen, Annalen der Physik 322, 549 (1905).
  * von Smoluchowski [1906] M. von Smoluchowski, Zur kinetischen theorie der brownschen molekularbewegung und der suspensionen, Annalen der Physik 326, 756 (1906).
  * Perrin [1909] J. Perrin, Mouvement brownien et réalité moléculaire, in  _Annales de Chimie et de Physique_ , Vol. 18 (1909) pp. 1–114.
  * Kappler [1931] E. Kappler, Versuche zur messung der avogadro-loschmidtschen zahl aus der brownschen bewegung einer drehwaage, Annalen der Physik 403, 233 (1931).
  * Wang  _et al._ [2009] B. Wang, S. M. Anthony, S. C. Bae, and S. Granick, Anomalous yet Brownian, Proc. Natl. Acad. Sci. USA 106, 15160 (2009).
  * Wang  _et al._ [2012] B. Wang, J. Kuo, S. C. Bae, and S. Granick, When Brownian diffusion is not Gaussian, Nat. Mater 11, 481 (2012).
  * Chechkin  _et al._ [2017] A. V. Chechkin, F. Seno, R. Metzler, and I. M. Sokolov, Brownian yet non-Gaussian diffusion: From superstatistics to subordination of diffusing diffusivities, Phys. Rev. X 7, 021002 (2017).
  * Guan  _et al._ [2014] J. Guan, B. Wang, and S. Granick, Even hard-sphere colloidal suspensions display Fickian yet non-Gaussian diffusion, ACS Nano 8, 3331 (2014).
  * Chakraborty  _et al._ [2019] I. Chakraborty, G. Rahamim, R. Avinery, Y. Roichman, and R. Beck, Nanoparticle mobility over a surface as a probe for weak transient disordered peptide–peptide interactions, Nano Lett. 19, 6524 (2019).
  * Chakraborty and Roichman [2020] I. Chakraborty and Y. Roichman, Disorder-induced Fickian, yet non-Gaussian diffusion in heterogeneous media, Phys. Rev. Res. 2, 022020(R) (2020).
  * Pastore  _et al._ [2022] R. Pastore, A. Ciarlo, G. Pesce, A. Sasso, and F. Greco, A model-system of Fickian yet non-Gaussian diffusion: light patterns in place of complex matter, Soft Matter 18, 351 (2022).
  * Rusciano  _et al._ [2022] F. Rusciano, R. Pastore, and F. Greco, Fickian non-Gaussian diffusion in glass-forming liquids, Phys. Rev. Lett. 128, 168001 (2022).
  * Ciarlo  _et al._ [2023] A. Ciarlo, R. Pastore, F. Greco, A. Sasso, and G. Pesce, Fickian yet non-Gaussian diffusion of a quasi-2D colloidal system in an optical speckle field: experiment and simulations, Sci. Rep. 13, 7408 (2023).
  * Yoshina-Ishii  _et al._ [2006] C. Yoshina-Ishii, Y.-H. M. Chan, J. M. Johnson, L. A. Kung, P. Lenz, and S. G. Boxer, Diffusive dynamics of vesicles tethered to a fluid supported bilayer by single-particle tracking, Langmuir 22, 5682 (2006).
  * Das  _et al._ [2009] R. Das, C. W. Cairo, and D. Coombs, A hidden markov model for single particle tracks quantifies dynamic interactions between lfa-1 and the actin cytoskeleton, PLoS Comput. Biol. 5, 1 (2009).
  * Ott  _et al._ [2013] M. Ott, Y. Shai, and G. Haran, Single-particle tracking reveals switching of the HIV fusion peptide between two diffusive modes in membranes, J. Phys. Chem. B 117, 13308 (2013).
  * Thapa  _et al._ [2018] S. Thapa, M. A. Lomholt, J. Krog, A. G. Cherstvy, and R. Metzler, Bayesian analysis of single-particle tracking data using the nested-sampling algorithm: maximum-likelihood model selection applied to stochastic-diffusivity data, Phys. Chem. Chem. Phys. 20, 29018 (2018).
  * Granik  _et al._ [2019] N. Granik, L. E. Weiss, E. Nehme, M. Levin, M. Chein, E. Perlson, Y. Roichman, and Y. Shechtman, Single-particle diffusion characterization by deep learning, Biophys. J. 117, 185 (2019).
  * Falcao and Coombs [2020] R. C. Falcao and D. Coombs, Diffusion analysis of single particle trajectories in a bayesian nonparametrics framework, Phys. Biol. 17, 025001 (2020).
  * Xu  _et al._ [2011] Q. Xu, L. Feng, R. Sha, N. C. Seeman, and P. M. Chaikin, Subdiffusion of a sticky particle on a surface, Phys. Rev. Lett. 106, 228102 (2011).
  * Skaug  _et al._ [2013] M. J. Skaug, J. Mabry, and D. K. Schwartz, Intermittent molecular hopping at the solid-liquid interface, Phys. Rev. Lett. 110, 256101 (2013).
  * Bernstein and Fricks [2016] J. Bernstein and J. Fricks, Analysis of single particle diffusion with transient binding using particle filtering, J. Theor. Biol. 401, 109 (2016).
  * Wang  _et al._ [2018] I.-H. Wang, C. J. Burckhardt, A. Yakimovich, and U. F. Greber, Imaging, tracking and computational analyses of virus entry and egress with the cytoskeleton, Viruses 10, 166 (2018).
  * Slator and Burroughs [2018] P. J. Slator and N. J. Burroughs, A hidden markov model for detecting confinement in single-particle tracking trajectories, Biophys. J. 115, 1741 (2018).
  * Callegari  _et al._ [2019] A. Callegari, C. Sieben, A. Benke, D. M. Suter, B. Fierz, D. Mazza, and S. Manley, Single-molecule dynamics and genome-wide transcriptomics reveal that NF-kB (p65)-DNA binding times can be decoupled from transcriptional activation, PLoS Genet. 15, e1007891 (2019).
  * Kowalek  _et al._ [2019] P. Kowalek, H. Loch-Olszewska, and J. Szwabinski, Classification of diffusion modes in single-particle tracking data: Feature-based versus deep-learning approach, Phys. Rev. E 100, 032410 (2019).
  * Jensen  _et al._ [2019] M. Jensen, Y.-Y. Wang., S. K. Lai, M. G. Forest, and S. A. McKinley, Antibody-mediated immobilization of virions in mucus, Bull. Math. Biol. 81, 4069 (2019).
  * Huseyin and Klose [2021] M. K. Huseyin and R. J. Klose, Live-cell single particle tracking of PRC1 reveals a highly dynamic system with low target site occupancy, Nat. Commun. 12, 887 (2021).
  * Doerries  _et al._ [2022a] T. J. Doerries, A. V. Chechkin, R. Schumer, and R. Metzler, Rate equations, spatial moments, and concentration profiles for mobile-immobile models with power-law and mixed waiting time distributions, Phys. Rev. E 105, 014105 (2022a).
  * Doerries  _et al._ [2022b] T. J. Doerries, A. V. Chechkin, and R. Metzler, Apparent anomalous diffusion and non-gaussian distributions in a simple mobile–immobile transport model with poissonian switching, J. R. Soc. Interface 19, 10.1098/rsif.2022.0233 (2022b).
  * Simon  _et al._ [2023] F. Simon, J.-Y. Tinevez, and S. van Teeffelen, ExTrack characterizes transition kinetics and diffusion in noisy single-particle tracks, J. Cell Biol. 222, e202208059 (2023).
  * Doerries  _et al._ [2023] T. J. Doerries, R. Metzler, and A. V. Chechkin, Emergent anomalous transport and non-gaussianity in a simple mobile–immobile model: the role of advection, New J. Phys. 25, 063009 (2023).
  * Federbush  _et al._ [2023] A. Federbush, A. Moscovich, and Y. Bar-Sinai, Github repository, https://github.com/ybs-lab/diffusion-with-tethering (2023).
  * Reichl [2016] L. E. Reichl, A modern course in statistical physics (Wiley-VCH GmbH, 2016) Chap. 7.3, 4th ed.
  * Li  _et al._ [2019] Q. Li, Z. Zhu, and G. Tang, Alternating minimizations converge to second-order optimal solutions, in  _36th Int. Conf. Mach. Learn. ICML 2019_, Vol. 97 (2019) pp. 3935–3943.
  * Dempster  _et al._ [1977] A. P. Dempster, N. M. Laird, and D. B. Rubin, Maximum likelihood from incomplete data via the EM algorithm, J. R. Stat. Soc.: Ser. B (Methodol.) 39, 1 (1977).
  * MacKay [2006] D. J. C. MacKay,  _Information Theory, Inference and Learning Algorithms_ (Cambridge University Press, 2006).
  * Viterbi [1967] A. Viterbi, Error bounds for convolutional codes and an asymptotically optimum decoding algorithm, IEEE Trans. Inf. Theory 13, 260 (1967).
  * Efron and Tibshirani [1994] B. Efron and R. Tibshirani,  _An Introduction to the Bootstrap_ (Chapman & Hall/CRC, Philadelphia, PA, 1994).
  * Seshadri and Sundberg [1994] N. Seshadri and C.-E. Sundberg, List Viterbi decoding algorithms with applications, IEEE Transactions on Communications 42, 313 (1994).
  * Roder and Hamzaoui [2006] M. Roder and R. Hamzaoui, Fast tree-trellis list Viterbi decoding, IEEE Transactions on Communications 54, 453 (2006).