[1]Arianna Ceccarelli

1]Mathematical Institute, University of Oxford, Woodstock Road, Oxford, OX2
6GG, UK

# Approximate solutions of a general stochastic velocity-jump process subject
to discrete-time noisy observations

0000-0002-9598-8845<EMAIL_ADDRESS>Alexander P. Browning
0000-0002-8753-1538 Ruth E. Baker  0000-0002-6304-9333 [

###### Abstract

Advances in experimental techniques allow the collection of high-space-and-
time resolution data that track individual motile entities over time. This
poses the question of how to use these data to efficiently and effectively
calibrate motion models. However, typical mathematical models often overlook
the inherent aspects of data collection, such as the discreteness and the
experimental noise of the measured locations. In this paper, we focus on
velocity-jump models suitable to describe single-agent motion in one spatial
dimension, characterised by successive Markovian transitions between a finite
network of n states, each with a specified velocity and a fixed rate of
switching to every other state. Since the problem of finding the exact
distributions of discrete-time noisy data is generally intractable, we derive
a series of approximations for the data distributions and compare them to in-
silico data generated by the models using four example network structures.
These comparisons suggest that the approximations are accurate given
sufficiently infrequent state switching, or equivalently, a sufficiently high
data collection frequency. Moreover, for infrequent switching, the PDFs
comparisons highlight the importance of accounting for the correlation between
subsequent measured locations, due to the likely permanence in the state
visited in the previous measurement. The approximate distributions computed
can be used for fast parameter inference and model selection between a range
of velocity-jump models using single-agent tracking data.

###### keywords:

Generalised velocity-jump process; continuous-time Markov chain; single-agent
tracking data; probability density function; approximate likelihood

##  1 Introduction

Mathematical modelling of motility is an area of huge interest across numerous
scientific domains, including ecology, biochemistry and cancer science.
Modelling has been used, for example, to characterize bacterial chemotaxis
[39, 38, 37, 12, 3], motion of molecular motors [13, 18, 17, 10, 22], axonal
transport [9, 46, 8, 32, 6, 5], RNA motility [26, 15, 14], cell migration [29,
19] and animal movement [31, 33, 42, 34, 24, 7, 20].

Nowadays, the spatial and temporal resolution of the data available is
increasing thanks to the continued development of experimental imaging
technologies [25]. In particular, experimental data that tracks the location
of individual motile entities over time is now readily available (see Figure 1
panels A, C, D), comprising a number of images taken at a fixed sampling rate.
These tracking data afford the potential to calibrate mathematical models and
quantify individual-level properties. However, calibration frameworks for
mathematical models that are continuous in space and time often neglect the
inherent aspects of data collection, such as fixed-time intervals and the
experimental noise of the measured locations. In this work, we describe
single-agent motion by formulating a continuous-time stochastic model, and we
provide approximations to its solution, characterising the probability density
function (PDF) of noisy location increments.

A

B

C

D

Figure 1: Examples of motion processes in biology. A. One-dimensional movement
of a kinesin molecular motor along a microtubule, taken from [16]. B. Motion
of a locust which follows a one-dimensional velocity-jump process, in which
the velocity switches between opposite directions, taken from [47]. C.
Bacterial tracks with angular reorientations, taken from [38]. D. Tracks of
fish movement in a tank, in which each dot represents a measured location in
three dimensions at 0.1 s intervals, taken from [31]. Panels A and D are
reprinted under a Creative Commons Attribution 4.0 International License.

For many systems, the behaviour of individual motile agents can be
characterised as a velocity-jump process, comprising a sequence of movements
at a constant velocity separated by instantaneous reorientations during which
the velocity is changed [43, 28]. The characteristics of velocity-jump models
are fully encoded in the re-orientation kernels, velocities and waiting time
distributions between reorientations. Random-walk in space models have been
used widely to model single-agent motion [19, 45, 4, 7, 20], and analytical
tools to estimate the model parameters and to characterise the agent motion
have been developed [36, 11, 45, 2, 41]. Alternatively, stochastic velocity-
jump models define the agent motion by defining the agent state stochastic
evolution between a network of states with specified velocities and switching
rates, and frameworks have been proposed to estimate the model parameters [14,
37], or asymptotic or mean quantities of interest [13, 21, 42].

In this paper, we focus on approximating the solutions of a general velocity-
jump model describing single-agent motion in one spatial dimension, which
entails successive transitions between $n$ states. Specifically, state
evolution follows a continuous-time Markov chain (CTMC) in which each state
specifies an agent velocity, rate of switching and probability to transition
to each of the other states upon switching (this can be seen in Figure 1B for
a type of motion that shows two alternating velocities). Additionally, we
propose a data collection model aligned with typical experiments,
characterised by a fixed time step between collected frames and a noisy
spatial resolution for the agent location. This data collection model allows
us to generate one-dimensional in-silico data tracks and extract location
increments, defined as the difference between two subsequent locations. We
then describe these data by approximating the PDF of the location increments.

First, we compute an approximation for the marginal distribution of measuring
a single location increment, and then we compute the joint distribution for a
set of subsequent location increments. The challenges of computing these PDFs
stem from the inherent noise and the experimental constraint that the agent
location is measured with a discrete time frequency, preventing the
determination of exact switching times which are necessary to directly compute
velocities and switching rates from the data. Furthermore, to assess the
validity of the approximate PDFs produced we compare them to the empirical PDF
as the switching rates vary for different model networks, as we expect them to
be accurate when switching rates are sufficiently low relative to the sampling
rate.

The paper is organised as follows. In Section 2 we formulate the $n$-state
velocity-jump model and propose a data collection model based on the intrinsic
characteristics of typical experimental tracking data, obtained with a fixed
time frequency and measurement noise, and we provide approximations for the
PDFs of single-agent tracking data. In particular, in Section 3 we provide
approximations for the PDF of a single location increment, which are based on
considering a finite number of state switches per measured interval. We
compare these approximations for four network models by plotting them against
the empirical PDF considering case in which transitions occur both
infrequently and frequently. In Section 4 we produce an approximation for the
joint PDF of a set of subsequent location increments, and we test it by
plotting it against the empirical PDF for the models previously considered. In
Section 5 we conclude that the approximate PDFs computed capture the
distribution of the data given that switching is sufficiently infrequent
compared to the time resolution of the data. More generally, the $n$-state
model presented is suitable to describe velocity-jump motion, and the
approximate PDFs provided could be used to calibrate the model using single-
agent tracking data.

##  2 Formulation of the n-state velocity-jump model

In this section, we present a one-dimensional $n$-state velocity-jump model
suitable to be calibrated with typically measured data (see, for example, the
tracking data in Figure 1A). At any point in time $t\geq 0$, the agent is in a
state $S(t)=s$, with $s\in\\{1,2,\ldots,n\\}$, and in this state it moves with
constant velocity $v_{s}\in\mathbb{R}$. We take velocities to be scalars for
simplicity, however, the model as described holds for higher-dimensional
velocity vectors, and the methods presented in this paper could be extended to
higher spatial dimensions.

The sequence of states attained is a CTMC, a continuous stochastic process in
which the time between changes of state is an exponential random variable, and
the new state attained is sampled according to the probabilities specified in
a stationary transition matrix [1, 23, 27]. In particular, the amount of time
spent in state $s$ before switching to a different state is an exponential
random variable $\tau_{s}\sim\mathrm{Exp}{(\lambda_{s})}$, for
$\lambda_{s}\neq 0$ constant. We note that
$\mathbb{E}[\tau_{s}]=1/\lambda_{s}$ is the average time spent in state $s$.

The sequence of states can be thought of as a discrete-time Markov chain
(DTMC). We denote with $p_{su}$ the probability of switching from state $s$ to
state $u$

$p_{su}:=\mathbb{P}\left(S_{k}=u\;\middle|\;S_{k-1}=s\right),$ (1)

for $k\geq 2$, where $S_{k}$ denotes the $k$-th state attained. The
probabilities $p_{su}$ are constant; thus they do not depend on $k$, and the
Markovian property holds. By construction, $p_{ss}:=0$, and
$\sum_{u=1}^{n}p_{su}:=1$. We write the transition matrix of the DTMC as

$\boldsymbol{P}:=\begin{bmatrix}0&p_{12}&p_{13}&\ldots&p_{1n}\\\
p_{21}&0&p_{23}&\ldots&p_{2n}\\\ p_{31}&p_{32}&0&\ldots&p_{3n}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
p_{n1}&p_{n2}&p_{n3}&\ldots&0\end{bmatrix}.$

For every state $u\neq s$, we define $q_{su}:=\lambda_{s}p_{su}$, and we
define $q_{ss}:=-\lambda_{s}$ such that

$\sum_{u=1}^{n}q_{su}=0,$

for all $s$. These $q_{su}$ are then the entries of the transition-rate matrix
of the CTMC

$\boldsymbol{Q}:=\begin{bmatrix}-\lambda_{1}&\lambda_{1}p_{12}&\lambda_{1}p_{13}&\ldots&\lambda_{1}p_{1n}\\\
\lambda_{2}p_{21}&-\lambda_{2}&\lambda_{2}p_{23}&\ldots&\lambda_{2}p_{2n}\\\
\lambda_{3}p_{31}&\lambda_{3}p_{32}&-\lambda_{3}&\ldots&\lambda_{3}p_{3n}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\lambda_{n}p_{n1}&\lambda_{n}p_{n2}&\lambda_{n}p_{n3}&\ldots&-\lambda_{n}\end{bmatrix}.$

Without loss of generality, starting at a time $t_{0}=0$, we define the
initial location of the agent as $x(0)=0$. Denoting $\hat{T}_{s}(t)$ for
$s=1,2,\ldots,n$ as the total amount of time spent in state $s$ in the time
interval $[0,t]$ we have

$x(t)=\sum_{s=1}^{n}v_{s}\hat{T}_{s}(t).$

We assume to have no knowledge of the agent initial state, which corresponds
to assuming that the process is in equilibrium. The initial state
$S_{0}:=S(0)$ is randomly sampled from the stationary distribution of the
CTMC, described by a vector $\boldsymbol{\pi}:=[p_{1},p_{2},\ldots,p_{n}]$
defining the probability of being in each state $s=1,2,\ldots,n$, constructed
such that $\boldsymbol{\pi}\boldsymbol{Q}=\boldsymbol{0}$. This choice of
$\boldsymbol{\pi}$ guarantees the time independence of the process, which
translates to having a constant marginal probability of being at each state
over time

$\mathbb{P}(S(t)=s)=:\mathbb{P}(s)=p_{s},$ (2)

for all $t\geq 0$. The existence and uniqueness of the vector
$\boldsymbol{\pi}$ are guaranteed by assuming that the chain is irreducible
(see Supplementary Information Section LABEL:SI:Sec:Irreducibility and [23]
Proposition 2.59 and Corollary 2.67). In this assumption, the vector
$\boldsymbol{\pi}$ can be explicitly computed by finding the kernel of
$\boldsymbol{Q}^{T}$ (see Supplementary Information Section
LABEL:SI:Sec:P(s)).

In the model outlined here, the state evolution is a CTMC, the properties of
which have been widely studied [35, 1, 23, 27]. However, existing mathematical
models often overlook aspects of data collection, such as the fixed-time
intervals and the experimental noise of the measured locations. In Section 2.1
we introduce a measurement model for data collection, designed to capture
these inherent traits of experimental setups.

Before proceeding, we first introduce four example networks that will be used
throughout this work to illustrate the validity of the approximate solutions
computed. These model networks are illustrated in Figure 2A-D, while the
parameter sets are fully specified in Supplementary Information Figure
LABEL:SI:Fig:networks. In these figures, the nodes indicate the states and the
black arrows represent a non-zero probability of switching between states,
after a time sampled from an exponential distribution. The simplest model on
the left is a two-state model (Figure 2A), in which the state switches between
a forward state (F), in which the agent has velocity $v_{F}>0$, and a backward
state (B), in which it has velocity $v_{B}<0$. The other networks presented
are obtained adding stationary states to this two-state network. The three-
state model (Figure 2B) is obtained by adding a stationary state with long
average permanence (SL), which corresponds to having a low switching rate
compared to the other states. The four-state model (Figure 2C) is obtained by
adding a stationary state with short average permanence (SS) which corresponds
to a higher switching rate compared to state SL. In the four-state model
(Figure 2D) the states F and B are not directly connected; rather, an agent
needs to enter a stationary state before changing direction. Moreover,
entering state SL implies that the following state attained must be SS, and
therefore the total time spent in a stationary phase follows a hypoexponential
distribution, as it is obtained as the sum of two exponential distributions.
Finally, the six-state network is given a cyclic structure, and it
incorporates stationary states followed by a very likely return to the
previous velocity state (F or B), denoted as pause while moving forward (PF)
and pause while moving backward (PB).

Two-state model

Three-state model

Four-state model

Six-state model

-

A

B

C

D

-

E

F

G

H

-

I

L

M

N

Figure 2: Panels A-D describe the model networks we use throughout the
manuscript to demonstrate the results obtained, more details are included in
Supplementary Information Figure LABEL:SI:Fig:networks. Panels E-N show
examples of in-silico data tracks from the model networks specified above,
with two, three, four and six states. F stands for forward state, B for
backward state, SL for stationary state with long average permanence, and SS
for stationary state with short average permanence, PF for pause while moving
forward and PB for pause while moving backward. The tracks are generated using
the sets of parameters specified in Supplementary Information Figure
LABEL:SI:Fig:networks, except the switching rates
$\boldsymbol{\lambda}=[\lambda_{1},\lambda_{2},\ldots,\lambda_{n}]$,
multiplied by 10 to obtain the tracks in panels I-N. In an experimental
setting increasing the switching rates is equivalent to increasing the time
between frames $\Delta t$.

###  2.1 Data collection model

Here, we define a data collection model that mimics the data obtained from
many experiments (see Figure 1A and D), in which measurements are noisy and
can only be obtained at discrete times, and highlight the associated
challenges. Indeed, we cannot determine the exact switching times and,
therefore, cannot directly infer velocities, switching rates and transition
probabilities from these data without using a model. The data collection model
will be used to obtain in-silico data (see Figure 2) and to validate the
approximate solutions for the model in Section 3 and Section 4.

By simulating the CTMC we generate the agent location $x(t)$, for
$t\in[0,N\Delta t]$ with $N\in\mathbb{N}$, and we obtain a set of $N+1$ data
points which are the agent locations measured with fixed time frequency
$\Delta t$ (see Figure 3). For $j\in\\{0,1,2,\ldots,N\\}$, at time
$t_{j}:=j\Delta t$, we denote the exact location of the agent as
$x_{j}:=x(t_{j})$, and the noisy location as

$y_{j}:=y(t_{j}):=x_{j}+\epsilon_{j},$

where the noise $\epsilon_{j}\sim\mathcal{N}(0,\sigma^{2})$ is normally
distributed with zero mean and variance $\sigma^{2}$, independent of $j$,
which is a standard assumption [30]. The methods presented generalise to other
noise model distributions for which the Markov property used still holds. Some
examples of in-silico tracks are shown in Figure 2. Figure 3 highlights that
since the data is noisy and discrete we cannot determine the exact switching
times.

Figure 3: Visual description of the data collection model. At time
$t_{j}=j\Delta t$ the agent is at location $x_{j}=x(j\Delta t)$. The $j$-th
exact location increment is defined as $\Delta x_{j}:=x_{j}-x_{j-1}$. The
$j$-th noisy location is defined as $y_{j}:=x_{j}+\epsilon_{j}$, where
$\epsilon_{j}$ is the measurement error, assumed to be normally distributed.
The $j$-th noisy location increment is defined as $\Delta
y_{j}:=y_{j}-y_{j-1}$.

We now present more notation that will be useful in the following sections. We
define the total simulation time as $T:=N\Delta t$. For
$j\in\\{1,2,\ldots,N\\}$, we define the exact location increment as

$\Delta x_{j}:=x_{j}-x_{j-1},$

and the noisy location increment as

$\Delta y_{j}:=y_{j}-y_{j-1}=\Delta
x_{j}+\Delta\epsilon_{j}\sim\mathcal{N}(\Delta x_{j},2\sigma^{2}),$

since
$\Delta\epsilon_{j}:=\epsilon_{j}-\epsilon_{j-1}\sim\mathcal{N}(0,2\sigma^{2})$.

We note that subsequent location increments are correlated since the previous
increments contain information on the state, and therefore on the velocity at
the beginning of the following increment (see Figure 3). In particular, if the
switching rates are small compared to the time interval between frames then
the state switches occur infrequently. Equivalently,
$\mathbb{E}[\tau_{s}]=1/\lambda_{s}\gg\Delta t$ for all $s$, and in this case
the velocity at the beginning of an increment is likely to be the same as that
at the beginning of the previous increment. We note that
$1/\lambda_{s}\gg\Delta t$ can be obtained experimentally by appropriately
choosing the time between captured data frames $\Delta t$.

The characteristics of the data collection model are considered in the
following sections to compute approximations for the PDFs of both a single
noisy location increment (Section 3) and a set of subsequent noisy location
increments obtained from a data track (Section 4).

##  3 Approximate solutions for the n-state velocity-jump model

In this section, we aim to compute the solutions for the $n$-state velocity-
jump model proposed in Section 2 to describe the distribution of measuring a
single noisy location increment $\Delta y$, denoted as $\mathbb{P}(\Delta y)$.

We can write the probability of measuring a single location increment, $\Delta
y$,

$\mathbb{P}(\Delta y)=\sum_{w=0}^{\infty}\mathbb{P}(\Delta
y\,|\,W=w)\mathbb{P}(W=w),$ (3)

where $\mathbb{P}(\Delta y\,|\,W=w)$ is the density of $\Delta y$ conditioned
on the number of switches $W$ during the given interval. The problem of
finding an exact solution for the PDF $\mathbb{P}(\Delta y)$ is intractable.
From Equation (3) we can observe that there are an infinite number of terms to
compute for which we cannot obtain a general formula. Hence, we define
approximations that consider up-to-$m$ switches per measured interval limiting
the amount of terms to compute. We rewrite Equation (3) as

$\mathbb{P}(\Delta y)=\sum_{w=0}^{m}\mathbb{P}(\Delta
y\,|\,W=w)\mathbb{P}(W=w)+\mathbb{P}(\Delta y\,|\,W>m)\mathbb{P}(W>m),$ (4)

and we define PDF approximations that assume

$\mathbb{P}(\Delta y\,|\,W>m)\approx\mathbb{P}(\Delta y\,|\,W=m).$

Using it, we obtain what we call an up-to-$m$-switch approximation
$P_{m}(\Delta y)$, given by

$P_{m}(\Delta y):=\sum_{w=0}^{m-1}\mathbb{P}(\Delta
y\,|\,W=w)\mathbb{P}(W=w)+\mathbb{P}(\Delta y\,|\,W=m)\mathbb{P}(W\geq m).$
(5)

We expect this approximation to work well either when $\mathbb{P}(\Delta
y\,|\,W>m)\approx\mathbb{P}(\Delta y\,|\,W=m)$ or when $\mathbb{P}(W>m)\approx
0$. Alternatively, we could also define a truncated PDF approximation by
considering only the terms with at most $m$ switches, using the approximation
$\mathbb{P}(W>m)\approx 0$, but it gives worse results than approximation
defined as in Equation (5). We now give some indication on the proportion of
intervals with a set number of switches, $\mathbb{P}(W=w)$, by focusing on the
two-state model and the three-state model previously introduced (see
Supplementary Information Figure LABEL:SI:Fig:networks for full details).

Two-state model

Three-state model

Figure 4: Empirical distributions of $\mathbb{P}(W=w)$, using the rates
$\boldsymbol{\lambda}=[\lambda_{1},\lambda_{2},\ldots,\lambda_{n}]$ specified
in Supplementary Information Figure LABEL:SI:Fig:networks, for the two-state
model (from Figure 2A) and three-state model (from Figure 2B).

The plots in Figure 4 illustrate the probability of having a number of
switches $w$, $\mathbb{P}(W=w)$, in an interval of fixed time length $\Delta
t$, for the two-state model and three-state model with rates
$\boldsymbol{\lambda}$ varied by the factors 1, 5 and 10. Parameters and
networks are used as specified in Supplementary Information Figure
LABEL:SI:Fig:networks. Equivalently, these plots (Figure 4) indicate how the
probability $\mathbb{P}(W=w)$ varies multiplying the measured time length
$\Delta t$ by the factors 1, 5 and 10, keeping the first set of rates
$\boldsymbol{\lambda}\times 1$.

We now approximate the PDF of a noisy location increment $\Delta y$,
$\mathbb{P}(\Delta y)$, by considering both an up-to-one-switch approximation
($P_{1}(\Delta y)$; Section 3.1) and an up-to-two-switch approximation
($P_{2}(\Delta y)$; Section 3.2).

###  3.1 Up-to-one-switch approximation for the probability distribution
function of a location increment

We now consider the up-to-one-switch approximation

$P_{1}(\Delta y):=\mathbb{P}(\Delta
y\,|\,W=0)\mathbb{P}(W=0)+\mathbb{P}(\Delta y\,|\,W=1)\mathbb{P}(W\geq 1).$
(6)

We make progress by further conditioning each term in Equation (6) on the
state at the start of the interval, denoted as $S_{1}$. Namely, for $W=0$, we
see that

$\mathbb{P}(\Delta y\,|\,W=0)\mathbb{P}(W=0)=\sum_{s=1}^{n}\mathbb{P}(\Delta
y\,|\,W=0,S_{1}=s)\mathbb{P}(W=0\,|\,S_{1}=s)\mathbb{P}(S_{1}=s),$

where $\mathbb{P}(S_{1}=s)=p_{s}$ is given by the equilibrium assumption
(Equation (2)). For $W=0$, the amount of time spent in the first state is,
trivially, given by $\Delta t$ (see Figure 5 zero-switch case). Therefore we
have

$\mathbb{P}(W=0\,|\,S_{1}=s)=\exp{(-\lambda_{s}\Delta t)},$

and

$\mathbb{P}(\Delta y\,|\,W=0,S_{1}=s)=f_{\mathcal{N}(v_{s}\Delta
t,2\sigma^{2})}(\Delta y),$

where the right hand side denotes the PDF of the normal distribution with mean
$v_{s}\Delta t$ and variance $2\sigma^{2}$. For more details on these results
see Supplementary Information Section LABEL:SI:Sec:P_1.

Figure 5: Diagrams representing the state switching times in the cases of zero
switches, one switch and two switches. The variable $\tau_{1}>0$ indicates the
time spent in state $S_{1}$, $\tau_{2}>0$ the time spent in state $S_{2}$ and
$\tau_{3}>0$ the time spent in state $S_{3}$.

To compute $\mathbb{P}(\Delta y\,|\,W=1)\mathbb{P}(W\geq 1)$, we condition on
both the first state and the second state visited within the interval, $S_{1}$
and $S_{2}$ (see one-switch case in Figure 5). This yields

$\displaystyle\mathbb{P}(\Delta y\,|\,W=1)\mathbb{P}(W\geq
1)=\sum_{\begin{subarray}{c}s_{1}\\\ s_{2}\neq s_{1}\end{subarray}}$
$\displaystyle\mathbb{P}(\Delta y\,|\,W=1,S_{1}=s_{1},S_{2}=s_{2})$
$\displaystyle\times\mathbb{P}(W\geq 1\,|\,S_{1}=s_{1},S_{2}=s_{2})$
$\displaystyle\times\mathbb{P}(S_{2}=s_{2}\,|\,S_{1}=s_{1})\mathbb{P}(S_{1}=s_{1}).$

By definition, $\mathbb{P}(S_{2}=s_{2}\,|\,S_{1}=s_{1})=p_{s_{1}s_{2}}$
(Equation (1)), and we have

$\mathbb{P}(W\geq 1\,|\,S_{1}=s_{1},S_{2}=s_{2})=\mathbb{P}(W\geq
1\,|\,S_{1}=s_{1})=1-\exp{(-\lambda_{s_{1}}\Delta t)}$

by considering that at least one switch occurs if $\tau_{1}<\Delta t$.

Computation of the PDF for $\Delta y$ if one switch occurs and given the two
states visited, denoted as

$\tilde{f}_{s_{1},s_{2}}(\Delta y):=\mathbb{P}(\Delta
y\,|\,W=1,S_{1}=s_{1},S_{2}=s_{2}),$

is significantly more involved. We proceed by considering the corresponding
PDF for the exact increment

$\tilde{g}_{s_{1},s_{2}}(\Delta x):=\mathbb{P}(\Delta
x\,|\,W=1,S_{1}=s_{1},S_{2}=s_{2}),$

and then convoluting the result with the distribution for $\Delta\epsilon$ to
obtain $\tilde{f}_{s_{1},s_{2}}(\Delta y)$.

We construct the PDF $\tilde{g}_{s_{1},s_{2}}(\Delta x)$ by first computing
the distribution of the switching time $\tau_{1}$ during the time interval
considered $[0,\Delta t]$, with distribution given by

$\displaystyle F_{s_{1},s_{2}}(t_{1}):=$ $\displaystyle\
\mathbb{P}(\tau_{1}\leq t_{1},W=1\,|\,S_{1}=s_{1},S_{2}=s_{2})$
$\displaystyle=$ $\displaystyle\ \mathbb{P}(0<\tau_{1}\leq
t_{1},\tau_{2}>\Delta t-\tau_{1}\,|\,S_{1}=s_{1},S_{2}=s_{2})$
$\displaystyle=$ $\displaystyle\ \int_{0}^{t_{1}}\int_{\Delta
t-\tau_{1}}^{\infty}f_{s_{1},s_{2}}(\tau_{1},\tau_{2})\text{d}\tau_{2}\text{d}\tau_{1},$

for $t_{1}\in(0,\Delta t]$, where
$f_{s_{1},s_{2}}(\tau_{1},\tau_{2})=f_{\text{Exp}(\lambda_{s_{1}})}(\tau_{1})\cdot
f_{\text{Exp}(\lambda_{s_{2}})}(\tau_{2})$ is the joint distribution of the
time spent in the first and second states which are independent (see two-
switch case in Figure 5). Moreover, we obtain the cumulative distribution
function (CDF) of the switching time as

$\displaystyle G_{s_{1},s_{2}}(t_{1}):=$ $\displaystyle\
\mathbb{P}(\tau_{1}\leq t_{1}\,|\,W=1,S_{1}=s_{1},S_{2}=s_{2})$
$\displaystyle=$ $\displaystyle\ \mathbb{P}(\tau_{1}\leq
t_{1}\,|\,0\leq\tau_{1}\leq\Delta t,\tau_{2}>\Delta
t-\tau_{1},S_{1}=s_{1},S_{2}=s_{2})$ $\displaystyle=$ $\displaystyle\
\frac{F_{s_{1},s_{2}}(t_{1})}{F_{s_{1},s_{2}}(\Delta t)}.$

This CDF can be differentiated to obtain the PDF of the switching time
$\tau_{1}=t_{1}$, from which we obtain the PDF of the exact location increment
$\Delta x$, $\tilde{g}_{s_{1},s_{2}}(\Delta x)$. We have that, for fixed
$S_{1}$, $S_{2}$ and switching at time $\tau_{1}$,

$\Delta x=h_{s_{1},s_{2}}(\tau_{1}):=v_{s_{1}}\tau_{1}+v_{s_{2}}(\Delta
t-\tau_{1}).$

If $v_{s_{1}}=v_{s_{2}}$, then the velocity is constant for the whole
interval; thus the exact increment is known, $\Delta x=h_{s_{1},s_{2}}(\Delta
t)$, and the distribution for the measured increment is

$\tilde{f}_{s_{1},s_{2}}(\Delta y)=f_{\mathcal{N}({v_{s_{1}}\Delta
t},2\sigma^{2})}(\Delta y),$

where the right hand side denotes the PDF of the normal distribution.
Otherwise, when $v_{s_{1}}\neq v_{s_{2}}$, the exact increment $\Delta x$ is
determined by the time of the switch $\tau_{1}$. In this case, we note that
$h_{s_{1},s_{2}}$ is a monotonic function in $\tau_{1}$, then its inverse
exists and gives

$\tau_{1}=h_{s_{1},s_{2}}^{-1}(\Delta x);$

thus the PDF of $\Delta x$ is given by

$\tilde{g}_{s_{1},s_{2}}(\Delta
x)=g_{s_{1},s_{2}}\left(h_{s_{1},s_{2}}^{-1}(\Delta
x)\right)\cdot\left|\frac{\text{d}}{\text{d}(\Delta
x)}h_{s_{1},s_{2}}^{-1}(\Delta x)\right|,$

where $g_{s_{1},s_{2}}(t_{1})$ is the PDF obtained by differentiating the CDF
$G_{s_{1},s_{2}}(t_{1})$ in $t_{1}$.

Finally, we obtain the PDF for $\Delta y$, $\tilde{f}_{s_{1},s_{2}}(\Delta
y)$, by incorporating the Gaussian noise, computing it as the distribution of
the sum of two independent random variables $\Delta x$ and $\Delta\epsilon$.
Therefore, we integrate the product of the PDF for $\Delta x$ and the PDF for
$\Delta\epsilon=\Delta y-\Delta x$ for all $\Delta x$ and obtain

$\tilde{f}_{s_{1},s_{2}}(\Delta y)=\int_{a}^{b}\tilde{g}_{s_{1},s_{2}}(\Delta
x)f_{\mathcal{N}(0,2\sigma^{2})}(\Delta y-\Delta x)\text{d}(\Delta x),$ (7)

where we define

$a=a_{s_{1},s_{2}}:=\min\\{h_{s_{1},s_{2}}(0),h_{s_{1},s_{2}}(\Delta t)\\},$

and

$b=b_{s_{1},s_{2}}:=\max\\{h_{s_{1},s_{2}}(0),h_{s_{1},s_{2}}(\Delta t)\\},$

since $\tilde{g}_{s_{1},s_{2}}(\Delta x)=0$ for $\Delta x<a$ or $\Delta x>b$.
The integral in Equation (7) can be computed explicitly; this is provided in
Supplementary Information Section LABEL:SI:Sec:P_1.

-

Four-state model

Six-state model

$\boldsymbol{\lambda}$

$10\boldsymbol{\lambda}$

Figure 6: The empirical PDF for $\Delta y$, denoted $P(\Delta y)$, (orange
histogram) is compared with its up-to-one-switch approximation $P_{1}(\Delta
y)$ (black line). The left plots are obtained from a four-state model while
the right ones from a six-state model. The top plots are obtained using the
parameters as specified in the Supplementary Information Figure
LABEL:SI:Fig:networks, while the bottom ones use the same parameters except
the rates which are multiplied by 10.

In Figure 6 we present a comparison between the empirical PDF for $\Delta y$,
denoted $P(\Delta y)$, and the up-to-one-switch approximation $P_{1}(\Delta
y)$ for the four-state model and six-state model. To assess the performance of
the approximation as the rate parameters increase for a fixed measurement
interval, we also produce results with the rate parameters multiplied by 10.
These results show that the up-to-one-switch approximation works well for
infrequent switching relative to $\Delta t$ (low rates
$\boldsymbol{\lambda}$), while the bottom plots indicate that the accuracy
decreases as the rates are increased, as expected. This is expected since the
in-silico data contains more increments with at least two switches as the
rates are increased (for $\boldsymbol{\lambda}$ versus
$10\boldsymbol{\lambda}$ for the four-state model the increments with at least
two switches are 3% versus 59%, for the six-state model they are 6% versus
72%). Figure 7 shows that also for the two-state model and the three-state
model the accuracy of $P_{1}(\Delta y)$ decreases as the rates are increased.
Finally, the Supplementary Information Figure LABEL:SI:Fig:Error_PDF_4_6 shows
the error of the up-to-one-switch PDF approximation, defined as $P_{1}(\Delta
y)-P(\Delta y)$, for the panels shown in Figure 6.

To improve the accuracy of the approximation of $\mathbb{P}(\Delta y)$ for
higher rates, one could use approximations that increase the maximum number of
switches considered. In the next section, we compute an up-to-two-switch
approximation; these steps could be adapted to produce approximations with an
increased maximum number of switches.

###  3.2 Up-to-two-switch approximation for the probability distribution
function of a location increment

We now consider the up-to-two-switch approximation defined as in Equation (5)

$\displaystyle P_{2}(\Delta y):=$ $\displaystyle\ \mathbb{P}(\Delta
y\,|\,W=0)\mathbb{P}(W=0)+\mathbb{P}(\Delta y\,|\,W=1)\mathbb{P}(W=1)$
$\displaystyle+\mathbb{P}(\Delta y\,|\,W=2)\mathbb{P}(W\geq 2).$

In Section 3.1 and in the Supplementary Information, we computed the terms
needed to obtain $\mathbb{P}(\Delta y\,|\,W=0)\mathbb{P}(W=0)$ and
$\mathbb{P}(\Delta y\,|\,W=1)\mathbb{P}(W=1)$.

In order to compute $\mathbb{P}(\Delta y\,|\,W=2)\mathbb{P}(W\geq 2)$ we
proceed as before by conditioning on all three subsequent states visited
during an interval, $S_{1}$, $S_{2}$ and $S_{3}$ (two-switch case Figure 5),
and use the Markov property and the fact that
$\mathbb{P}(S_{3}=s_{3}\,|\,S_{1}=s_{1},S_{2}=s_{2})=\mathbb{P}(S_{3}=s_{3}\,|\,S_{2}=s_{2})=p_{s_{2}s_{3}}$
(Equation (1)), $\mathbb{P}(S_{2}=s_{2}\,|\,S_{1}=s_{1})=p_{s_{1}s_{2}}$
(Equation (1)), and $\mathbb{P}(S_{1}=s_{1})=p_{s_{1}}$ (Equation (2)), to
obtain

$\displaystyle\mathbb{P}(\Delta y\,|\,W=2)\mathbb{P}(W\geq
2)=\sum_{\begin{subarray}{c}s_{1}\\\ s_{2}\neq s_{1}\\\ s_{3}\neq
s_{2}\end{subarray}}$ $\displaystyle\mathbb{P}(\Delta
y\,|\,W=2,S_{1}=s_{1},S_{2}=s_{2},S_{3}=s_{3})$
$\displaystyle\times\mathbb{P}(W\geq
2\,|\,S_{1}=s_{1},S_{2}=s_{2},S_{3}=s_{3})$
$\displaystyle\times\mathbb{P}(S_{3}=s_{3}\,|\,S_{1}=s_{1},S_{2}=s_{2})$
$\displaystyle\times\mathbb{P}(S_{2}=s_{2}\,|\,S_{1}=s_{1})\mathbb{P}(S_{1}=s_{1}).$

Hence, we only need to compute the probability of having two switches
$\mathbb{P}(W\geq 2\,|\,S_{1}=s_{1},S_{2}=s_{2},S_{3}=s_{3})$ (see
Supplementary Information Section LABEL:SI:Sec:P_2).

As for $W=1$, the PDF for $\Delta y$ if two switches occur and given the three
states visited, denoted as

$\tilde{f}_{s_{1},s_{2},s_{3}}(\Delta y):=\mathbb{P}(\Delta
y\,|\,W=2,S_{1}=s_{1},S_{2}=s_{2},S_{3}=s_{3}),$

is obtained by first considering the PDF for an exact increment

$\tilde{g}_{s_{1},s_{2},s_{3}}(\Delta x):=\mathbb{P}(\Delta
x\,|\,W=2,S_{1}=s_{1},S_{2}=s_{2},S_{3}=s_{3}),$

and then convoluting with the distribution for $\Delta\epsilon$. We proceed as
before, although note that we must deal with the case $s_{3}=s_{1}$
separately, since it leads to different time and increment distributions (see
Supplementary Information Section LABEL:SI:Sec:P_2). Here, we focus on the
case $s_{3}\neq s_{1}$, and we assume that the velocities are all distinct and
the rates are all distinct, as in the two-state model and three-state model
(Figure 2A-B).

We first determine the joint distribution of the time $\tau_{i}$ spent in each
state $S_{i}$, $i=1,2,3$, within the interval $[0,\Delta t]$ for the case of
two switches ($W=2$). We note that the number of switches is two if and only
if $0<\tau_{1}\leq\Delta t$, $0<\tau_{2}\leq\Delta t-\tau_{1}$ and
$\tau_{3}>\Delta t-\tau_{1}-\tau_{2}$. Following a similar approach to before,
we condition on the three subsequent states attained. For
$t_{1},t_{2}\in(0,\Delta t]$ with $t_{2}\leq\Delta t-t_{1}$, we compute

$\displaystyle F_{s_{1},s_{2},s_{3}}(t_{1},t_{2}):=$ $\displaystyle\
\mathbb{P}(\tau_{1}\leq t_{1},\tau_{2}\leq
t_{2},W=2\,|\,S_{1}=s_{1},S_{2}=s_{2},S_{3}=s_{3})$ $\displaystyle=$
$\displaystyle\ \mathbb{P}\left(\ \begin{aligned} &0<\tau_{1}\leq
t_{1},0<\tau_{2}\leq t_{2},\\\ &\tau_{3}>\Delta
t-\tau_{1}-\tau_{2}\end{aligned}\;\middle|\;\begin{aligned}
S_{1}=s_{1},S_{2}=s_{2},S_{3}=s_{3}\end{aligned}\right)$ $\displaystyle=$
$\displaystyle\ \int_{0}^{t_{1}}\left(\int_{0}^{t_{2}}\left(\int_{\Delta
t-\tau_{1}-\tau_{2}}^{\infty}\prod_{i=1}^{3}f_{\text{Exp}(\lambda_{s_{i}})}(\tau_{i})\text{d}\tau_{3}\right)\text{d}\tau_{2}\right)\text{d}\tau_{1},$

by integrating the joint distribution of the time spent in the first, second
and third states that are again independent. Hence, we can compute the joint
CDF of the times when the two switches occur given the three states attained
as follows

$\displaystyle G_{s_{1},s_{2},s_{3}}(t_{1},t_{2}):=$ $\displaystyle\
\mathbb{P}(\tau_{1}\leq t_{1},\tau_{2}\leq
t_{2}\,|\,W=2,S_{1}=s_{1},S_{2}=s_{2},S_{3}=s_{3})$ $\displaystyle=$
$\displaystyle\ \mathbb{P}\left(\ \begin{aligned} &\tau_{1}\leq t_{1},\\\
&\tau_{2}\leq t_{2}\end{aligned}\;\middle|\;\begin{aligned}
&0<\tau_{1}\leq\Delta t,0<\tau_{2}\leq\Delta t-\tau_{1},\tau_{3}>\Delta
t-\tau_{1}-\tau_{2},\\\
&S_{1}=s_{1},S_{2}=s_{2},S_{3}=s_{3}\end{aligned}\right)$ $\displaystyle=$
$\displaystyle\
\dfrac{F_{s_{1},s_{2},s_{3}}(t_{1},t_{2})}{\mathbb{P}\left(\begin{aligned}
&0<\tau_{1}\leq\Delta t,0<\tau_{2}\leq\Delta t-\tau_{1},\\\ &\tau_{3}>\Delta
t-\tau_{1}-\tau_{2}\end{aligned}\;\middle|\;S_{1}=s_{1},S_{2}=s_{2},S_{3}=s_{3}\right)}$
$\displaystyle=$ $\displaystyle\
\dfrac{F_{s_{1},s_{2},s_{3}}(t_{1},t_{2})}{\int_{0}^{\Delta
t}\left(\int_{0}^{\Delta t-\tau_{1}}\left(\int_{\Delta
t-\tau_{1}-\tau_{2}}^{\infty}f_{s_{1},s_{2},s_{3}}(\tau_{1},\tau_{2},\tau_{3})\text{d}\tau_{3}\right)\text{d}\tau_{2}\right)\text{d}\tau_{1}}.$

Now, we compute the PDF for an exact increment $\Delta x$, integrating over
all possible switching times

$\displaystyle\tilde{g}_{s_{1},s_{2},s_{3}}(\Delta x)=$
$\displaystyle=\int_{0}^{\Delta t}\Bigg{(}\int_{0}^{\Delta
t-t_{1}}\mathbb{P}\left(\Delta
x\,|\,\tau_{1}=t_{1},\tau_{2}=t_{2},W=2,S_{1}=s_{1},S_{2}=s_{2},S_{3}=s_{3}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\times\mathbb{P}\left(t_{1},t_{2}\,|\,W=2,S_{1}=s_{1},S_{2}=s_{2},S_{3}=s_{3}\right)\text{d}t_{2}\Bigg{)}\text{d}t_{1}$
$\displaystyle=\int_{0}^{\Delta t}\left(\int_{0}^{\Delta
t-t_{1}}\delta\Big{(}v_{s_{1}}t_{1}+v_{s_{2}}t_{2}+v_{s_{3}}(\Delta
t-t_{1}-t_{2})-\Delta
x\Big{)}g_{s_{1},s_{2},s_{3}}(t_{1},t_{2})\text{d}t_{2}\right)\text{d}t_{1},$

where $\delta(\cdot)$ denotes the Dirac delta function and
$g_{s_{1},s_{2},s_{3}}(t_{1},t_{2})$ is the PDF obtained differentiating the
CDF $G_{s_{1},s_{2},s_{3}}(t_{1},t_{2})$ in $t_{1}$ and $t_{2}$. Using the
sifting property of the Dirac delta function we obtain

$\displaystyle\tilde{g}_{s_{1},s_{2},s_{3}}(\Delta x)=\int_{0}^{\Delta
t}\frac{1}{|v_{s_{3}}-v_{s_{2}}|}\
g_{s_{1},s_{2},s_{3}}\left(t_{1},\frac{(v_{s_{1}}-v_{s_{3}})t_{1}+v_{s_{3}}\Delta
t-\Delta x}{v_{s_{3}}-v_{s_{2}}}\right)\text{d}t_{1},$

where $g_{s_{1},s_{2},s_{3}}$ is defined to be zero for $(t_{1},t_{2})$
outside of $[0,\Delta t]\times[0,\Delta t-t_{1}]$.

In order to obtain $\tilde{g}_{s_{1},s_{2},s_{3}}(\Delta x)$ we need to
compute the indefinite integral in $t_{1}$

$\displaystyle I(t_{1}):=\int_{t_{1}}\frac{1}{|v_{s_{3}}-v_{s_{2}}|}\
g_{s_{1},s_{2},s_{3}}\left(t_{1},\frac{(v_{s_{1}}-v_{s_{3}})t_{1}+v_{s_{3}}\Delta
t-\Delta x}{v_{s_{3}}-v_{s_{2}}}\right)\text{d}t_{1},$

and define the domain of integration for $t_{1}$

$A:=\left\\{t_{1}\;\middle|\;t_{1}\in[0,\Delta
t],\frac{(v_{s_{1}}-v_{s_{3}})t_{1}+v_{s_{3}}\Delta t-\Delta
x}{v_{s_{3}}-v_{s_{2}}}\in[0,\Delta t-t_{1}]\right\\}.$

We denote
$v_{\min},v_{\text{int}},v_{\max}\in\\{v_{s_{1}},v_{s_{2}},v_{s_{3}}\\}$ such
that $v_{\min}<v_{\text{int}}<v_{\max}$. We note that for $\Delta
x\in[v_{\min}\Delta t,v_{\max}\Delta t]$, $A$ is an interval,
$A:=[E_{0},E_{1}]$. In Supplementary Information Section LABEL:SI:Sec:P_2, we
obtain that $E_{0}$ and $E_{1}$ are functions of $\Delta x$, piecewise linear
in the intervals $[v_{\min}\Delta t,v_{\text{int}}\Delta t]$ and
$[v_{\text{int}}\Delta t,v_{\max}\Delta t]$. The PDF for an exact increment is
obtained as

$\tilde{g}_{s_{1},s_{2},s_{3}}(\Delta x)=I(E_{1}(\Delta x))-I(E_{0}(\Delta
x)).$

Finally, we compute the PDF of a noisy increment $\Delta y$, again by
convoluting the PDFs for $\Delta x$ and $\Delta\epsilon$, to obtain

$\tilde{f}_{s_{1},s_{2},s_{3}}(\Delta y)=\ \int_{v_{\min}\Delta
t}^{v_{\max}\Delta t}\tilde{g}_{s_{1},s_{2},s_{3}}(\Delta
x)f_{\mathcal{N}(0,2\sigma^{2})}(\Delta y-\Delta x)\text{d}(\Delta x),$

whose analytical expression is presented in Supplementary Information Section
LABEL:SI:Sec:P_2.

-

Two-state model

Three-state model

-

A

B

$\boldsymbol{\lambda}$

-

C

D

$5\boldsymbol{\lambda}$

-

E

F

$10\boldsymbol{\lambda}$

Figure 7: The empirical PDF for $\Delta y$, denoted as $P(\Delta y)$, (orange
histogram) is compared with its up-to-one-switch approximation $P_{1}(\Delta
y)$ (black line) and its up-to-two-switch approximation $P_{2}(\Delta y)$ (red
dashed line). The left plots are obtained from a two-state model while the
right ones from a three-state model. The panels A-B are obtained using the
parameters as specified in the Supplementary Information Figure
LABEL:SI:Fig:networks, while the panels C-D and E-F are obtained with the same
parameters except the rates which are multiplied by 5 and 10, respectively.

In Figure 7 we present a comparison between the empirical PDF for
$\mathbb{P}(\Delta y)$, the up-to-one-switch approximation $P_{1}(\Delta y)$,
and the up-to-two-switch approximation $P_{2}(\Delta y)$, for the two-state
model and three-state model. We do this in order to assess the improvement
made by incorporating an extra switch in the approximation. From the plots in
panels A and B we can see that the up-to-one-switch and up-to-two-switch
approximations are comparable for infrequent switching (low rates
$\boldsymbol{\lambda}$), while looking at panels C, D and F we notice a higher
accuracy given by the two-switch approximation as the rates increase. This is
expected since the probability of having more than two switches
$\mathbb{P}(W>2)$ is higher for $\boldsymbol{\lambda}\times 10$ than for
$\boldsymbol{\lambda}\times 1$ and $\boldsymbol{\lambda}\times 5$ (see
$\mathbb{P}(W=w)$ in Figure 4). In contrast, panel E in Figure 7 shows that
going from an up-to-one-switch to an up-to-two-switch approximation does not
always lead to a better approximation. In particular, we note that the
accuracy for the two-state model for rates $10\boldsymbol{\lambda}$ is not
significantly improved since the probability of having more than two switches
is still quite high ($\mathbb{P}(W>2)\approx 0.32$), suggesting that the
approximation used ($\mathbb{P}(\Delta y\,|\,W>2)\approx\mathbb{P}(\Delta
y\,|\,W=2)$) is less accurate for the two-state model in this parameter
regime. Finally, the Supplementary Information Figure
LABEL:SI:Fig:Comparison_error_PDFs_2_3 shows a comparison between the error of
the up-to-one-switch PDF approximation, defined as $P_{1}(\Delta y)-P(\Delta
y)$, and the error of up-to-two-switch approximation, defined as $P_{2}(\Delta
y)-P(\Delta y)$, in the panels shown in Figure 7.

The approximations for $\mathbb{P}(\Delta y)$ give the distribution of single
location increments. In the next section, we use the results obtained to
approximate the probability distribution of a track by considering a set of
subsequent location increments.

##  4 Approximation for the probability distribution function of a set of
subsequent location increments

Now, we aim to compute an approximation for the joint PDF of a set of $N$
subsequent location increments. We compute an up-to-one-switch approximation
using the results obtained in Section 3.1, as it is simpler and quicker
compared to the up-to-two-switch approximation. We start by computing the PDF
of $N$ exact subsequent location increments

$\boldsymbol{\Delta x}_{N}=[\Delta x_{1},\Delta x_{2},\ldots,\Delta
x_{N-1},\Delta x_{N}],$

for any $N\in\\{1,2,\ldots\\}$. Moreover, $\boldsymbol{\Delta x}_{i}$ denotes
the first $i$ elements of $\boldsymbol{\Delta x}_{N}$, while $\Delta x_{i}$
its element in position $i$. We denote the PDF of $N$ exact subsequent
location increments with $\mathbb{P}(\boldsymbol{\Delta x}_{N})$, and a
similar notation is used for $\boldsymbol{\Delta y}$.

We start by rewriting the PDF in terms of conditional distributions taking
into account that the state at the end of a measured interval corresponds to
the state at the beginning of the following measured interval (see Figure 3).
We obtain

$\mathbb{P}(\boldsymbol{\Delta x}_{N})=\mathbb{P}(\Delta
x_{N}\,|\,\boldsymbol{\Delta x}_{N-1})\ldots\mathbb{P}(\Delta
x_{3}\,|\,\boldsymbol{\Delta x}_{2})\mathbb{P}(\Delta x_{2}\,|\,\Delta
x_{1})\mathbb{P}(\Delta x_{1}),$ (8)

where $\mathbb{P}(\Delta x_{1})$ corresponds to the marginal distribution for
a single increment $\Delta x_{1}$, approximated in Section 3.1 and Section
3.2. In this section, we first compute an approximation for the PDF of a $N$
exact location increments $\mathbb{P}(\boldsymbol{\Delta x}_{N})$. We show
that for any $N\geq 2$ we can construct the joint PDF of a set of exact
subsequent location increments $\boldsymbol{\Delta x}_{N}$,
$\mathbb{P}(\boldsymbol{\Delta x}_{N})$, using a recursive method. From now
on, we denote by $S^{j}_{i}$ the $i$-th state attained and $W^{j}$ the number
of switches that occur during the $j$-th interval (as in Figure 8).

Figure 8: Diagram representing the notation used for measured intervals of
time length $\Delta t$ such that the agent location is measured at the
beginning and end of such interval. We note that the state at the end of a
measured interval coincides with the first state attained in the following
measured interval.

We compute an iterative formula for $\mathbb{P}(\Delta
x_{N}\,|\,\boldsymbol{\Delta x}_{N-1})$, which does not only depend on
$\mathbb{P}(\Delta x_{N-1}\,|\,\boldsymbol{\Delta x}_{N-2})$, but also on the
probabilities of attaining a state or having a number of switches at the
previous measured intervals. In particular, we compute an iterative formula
for $\mathbb{P}(\Delta x_{N}\,|\,\boldsymbol{\Delta x}_{N-1})$ that involves
the approximations for

$\displaystyle\mathbb{P}(\Delta x_{N}\,|\,S^{N}_{1}=s^{N}_{1}),$
$\displaystyle\ \mathbb{P}(W^{N}=0\,|\,S^{N}_{1}=s^{N}_{1}),$
$\displaystyle\mathbb{P}(W^{N}\geq 1\,|\,S^{N}_{1}=s^{N}_{1}),$
$\displaystyle\ \mathbb{P}(S^{N}_{1}=s^{N}_{1}\,|\,\boldsymbol{\Delta
x}_{N-1}),$

using $\mathbb{P}(\Delta x_{N-1}\,|\,\boldsymbol{\Delta x}_{N-2})$, and

$\displaystyle\mathbb{P}(\Delta x_{N-1}\,|\,S^{N-1}_{1}=s^{N-1}_{1}),$
$\displaystyle\ \mathbb{P}(W^{N-1}=0\,|\,S^{N-1}_{1}=s^{N-1}_{1}),$
$\displaystyle\mathbb{P}(W^{N-1}\geq 1\,|\,S^{N-1}_{1}=s^{N-1}_{1}),$
$\displaystyle\ \mathbb{P}(S^{N-1}_{1}=s^{N-1}_{1}\,|\,\boldsymbol{\Delta
x}_{N-2}).$

We note that, for $N=2$ these terms above are again the marginal PDFs, for
which we have computed approximations in the previous sections. For
simplicity, we use the up-to-one-switch approximation presented in Section 3.1
which we denote by $P_{1}$, however, the result could be extended to
incorporate more switches.

We note that conditioning on the state at the beginning of the $N$-th
interval, $S^{N}_{1}$, makes the process at the $N$-th interval independent of
$\boldsymbol{\Delta x}_{N-1}$. Thus, we obtain

$\mathbb{P}(\Delta x_{N}\,|\,S^{N}_{1}=s^{N}_{1},\boldsymbol{\Delta
x}_{N-1})=\mathbb{P}(\Delta x_{N}\,|\,S^{N}_{1}=s^{N}_{1}),$ (9)

and, similarly, we obtain

$\mathbb{P}(S^{N}_{1}=s^{N}_{1}\,|\,S^{N-1}_{1}=s^{N-1}_{1},\boldsymbol{\Delta
x}_{N-1})=\mathbb{P}(S^{N}_{1}=s^{N}_{1}\,|\,S^{N-1}_{1}=s^{N-1}_{1},\Delta
x_{N-1}).$ (10)

Using Equation (9) we write

$\displaystyle\mathbb{P}(\Delta x_{N}\,|\,\boldsymbol{\Delta x}_{N-1})=$
$\displaystyle\ \sum_{s^{N}_{1}}\mathbb{P}(\Delta
x_{N}\,|\,S^{N}_{1}=s^{N}_{1},\boldsymbol{\Delta
x}_{N-1})\mathbb{P}(S^{N}_{1}=s^{N}_{1}\,|\,\boldsymbol{\Delta x}_{N-1})$
$\displaystyle=$ $\displaystyle\ \sum_{s^{N}_{1}}\mathbb{P}(\Delta
x_{N}\,|\,S^{N}_{1}=s^{N}_{1})\mathbb{P}(S^{N}_{1}=s^{N}_{1}\,|\,\boldsymbol{\Delta
x}_{N-1}).$

Moreover, we can compute
$\mathbb{P}(S^{N}_{1}=s^{N}_{1}\,|\,\boldsymbol{\Delta x}_{N-1})$ by
considering the state at the beginning of the previous interval $S^{N-1}_{1}$,

$\displaystyle\mathbb{P}(S^{N}_{1}=s^{N}_{1}\,|\,\boldsymbol{\Delta x}_{N-1})$
$\displaystyle=\sum_{s^{N-1}_{1}}\mathbb{P}(S^{N}_{1}=s^{N}_{1}\,|\,S^{N-1}_{1}=s^{N-1}_{1},\boldsymbol{\Delta
x}_{N-1})\mathbb{P}(S^{N-1}_{1}=s^{N-1}_{1}\,|\,\boldsymbol{\Delta x}_{N-1})$
$\displaystyle=\sum_{s^{N-1}_{1}}\frac{\mathbb{P}(\Delta
x_{N-1},S^{N}_{1}=s^{N}_{1}\,|\,S^{N-1}_{1}=s^{N-1}_{1})}{\mathbb{P}(\Delta
x_{N-1}\,|\,S^{N-1}_{1}=s^{N-1}_{1})}\mathbb{P}(S^{N-1}_{1}=s^{N-1}_{1}\,|\,\boldsymbol{\Delta
x}_{N-1}),$ (11)

where we use Equation (10).

The final term in Equation (11) can be obtained using Bayes’ theorem, such
that

$\displaystyle\mathbb{P}(S^{N-1}_{1}=s^{N-1}_{1}\,|\,\boldsymbol{\Delta
x}_{N-1})=$ $\displaystyle\ \mathbb{P}(\Delta
x_{N-1}\,|\,S^{N-1}_{1}=s^{N-1}_{1},\boldsymbol{\Delta x}_{N-2})$
$\displaystyle\times\frac{\mathbb{P}(S^{N-1}_{1}=s^{N-1}_{1}\,|\,\boldsymbol{\Delta
x}_{N-2})}{\mathbb{P}(\Delta x_{N-1}\,|\,\boldsymbol{\Delta x}_{N-2})}$
$\displaystyle=$ $\displaystyle\ \mathbb{P}(\Delta
x_{N-1}\,|\,S^{N-1}_{1}=s^{N-1}_{1})\frac{\mathbb{P}(S^{N-1}_{1}=s^{N-1}_{1}\,|\,\boldsymbol{\Delta
x}_{N-2})}{\mathbb{P}(\Delta x_{N-1}\,|\,\boldsymbol{\Delta x}_{N-2})},$

where the second equality is obtained using the property in Equation (9) at
the $(N-1)$-th step. We substitute this result into Equation (11), and, by
simplifying $\mathbb{P}(\Delta x_{N-1}\,|\,S^{N-1}_{1}=s^{N-1}_{1})$, we
obtain

$\displaystyle\mathbb{P}(S^{N}_{1}=s^{N}_{1}\,|\,\boldsymbol{\Delta x}_{N-1})$
$\displaystyle=\sum_{s^{N-1}_{1}}\mathbb{P}(\Delta
x_{N-1},S^{N}_{1}=s^{N}_{1}\,|\,S^{N-1}_{1}=s^{N-1}_{1})\frac{\mathbb{P}(S^{N-1}_{1}=s^{N-1}_{1}\,|\,\boldsymbol{\Delta
x}_{N-2})}{\mathbb{P}(\Delta x_{N-1}\,|\,\boldsymbol{\Delta x}_{N-2})}.$

Here, both the numerator and denominator of the fraction are computed at the
$(N-1)$-th induction step.

The final term that we must compute is $\mathbb{P}(\Delta
x_{N-1},S^{N}_{1}=s^{N}_{1}\,|\,S^{N-1}_{1}=s^{N-1}_{1})$. To do this, we
condition on the number of switches during the $(N-1)$-th interval, $W^{N-1}$,
and obtain

$\displaystyle\mathbb{P}(\Delta
x_{N-1},S^{N}_{1}=s^{N}_{1}\,|\,S^{N-1}_{1}=s^{N-1}_{1})$
$\displaystyle=\sum_{w=0}^{\infty}\mathbb{P}(\Delta
x_{N-1}\,|\,S^{N}_{1}=s^{N}_{1},S^{N-1}_{1}=s^{N-1}_{1},W^{N-1}=w)$
$\displaystyle\qquad\quad\times\mathbb{P}(S^{N}_{1}=s^{N}_{1}\,|\,S^{N-1}_{1}=s^{N-1}_{1},W^{N-1}=w)\mathbb{P}(W^{N-1}=w\,|\,S^{N-1}_{1}=s^{N-1}_{1}).$

For simplicity, we apply the up-to-one-switch approximation in Section 3 which
yields

$\displaystyle\mathbb{P}(\Delta
x_{N-1},S^{N}_{1}=s^{N}_{1}\,|\,S^{N-1}_{1}=s^{N-1}_{1})$
$\displaystyle\approx P_{1}(\Delta
x_{N-1},S^{N}_{1}=s^{N}_{1}\,|\,S^{N-1}_{1}=s^{N-1}_{1})$
$\displaystyle=\begin{dcases}\begin{aligned} &\mathbb{P}(\Delta
x_{N-1}\,|\,S^{N}_{1}=s^{N}_{1},S^{N-1}_{1}=s^{N-1}_{1},W^{N-1}=0)\\\ &\
\times\mathbb{P}(W^{N-1}=0\,|\,S^{N-1}_{1}=s^{N-1}_{1}),\end{aligned}&\text{if
}s^{N}_{1}=s^{N-1}_{1},\\\ \begin{aligned} &\mathbb{P}(\Delta
x_{N-1}\,|\,S^{N}_{1}=s^{N}_{1},S^{N-1}_{1}=s^{N-1}_{1},W^{N-1}=1)\\\ &\
\times\mathbb{P}(S^{N-1}_{2}=s^{N}_{1}\,|\,S^{N-1}_{1}=s^{N-1}_{1})\mathbb{P}(W^{N-1}\geq
1\,|\,S^{N-1}_{1}=s^{N-1}_{1}),\end{aligned}&\text{if }s^{N}_{1}\neq
s^{N-1}_{1},\end{dcases}$

since if $W^{N-1}=0$, then $S^{N}_{1}=S^{N-1}_{1}$, and this gives a unitary
probability of keeping the same state

$\mathbb{P}(S^{N}_{1}=s^{N}_{1}\,|\,S^{N-1}_{1}=s^{N-1}_{1},W^{N-1}=0)=1.$

Otherwise, if $W^{N-1}=1$, then we have $S^{N}_{1}=S^{N-1}_{2}\neq
S^{N-1}_{1}$ and therefore

$\mathbb{P}(S^{N}_{1}=s^{N}_{1}\,|\,S^{N-1}_{1}=s^{N-1}_{1},W^{N-1}=1)=\mathbb{P}(S^{N-1}_{2}=s^{N}_{1}\,|\,S^{N-1}_{1}=s^{N-1}_{1}).$

This result can be extended to compute an approximation for the PDF of $N$
noisy subsequent location increments $\boldsymbol{\Delta y}_{N}=[\Delta
y_{1},\Delta y_{2},\ldots,\Delta y_{N-1},\Delta y_{N}].$ The PDF of $N$ noisy
subsequent location increments can be written in terms of conditional
distributions as

$\mathbb{P}(\boldsymbol{\Delta y}_{N})=\mathbb{P}(\Delta
y_{N}\,|\,\boldsymbol{\Delta y}_{N-1})\ldots\mathbb{P}(\Delta
y_{3}\,|\,\boldsymbol{\Delta y}_{2})\mathbb{P}(\Delta y_{2}\,|\,\Delta
y_{1})\mathbb{P}(\Delta y_{1}),$ (12)

where we note that $\mathbb{P}(\Delta y_{1})$ is the marginal distribution for
a single increment $\Delta y_{1}$.

The distributions of $\Delta y$ are obtained by convoluting the distributions
for $\Delta x$ with the appropriate Gaussian noise. The steps proposed in this
section are valid for the noisy increments, with only the exception of the
properties in Equation (9) and Equation (10). Equation (9) notes that the
exact location increment $\Delta x_{N}$, conditioned on the state $S_{1}^{N}$,
is independent of the previous exact location increments $\boldsymbol{\Delta
x}_{N-1}$, while Equation (9) notes that the state $S_{1}^{N}$, conditioned on
the state $S_{1}^{N-1}$ and on the increment $\Delta x_{N-1}$, is independent
of the previous exact location increments $\boldsymbol{\Delta x}_{N-2}$. These
properties hold for exact increments, but are approximations for noisy
measured increments, since the previous measured increments contain additional
information on the noise of the following measured increment. For simplicity,
when computing the PDF of a set of noisy subsequent location increments we
assume the correlation of the noise to be negligible. In practice, this
corresponds to assuming that $\Delta\epsilon_{N}:=\epsilon_{N}-\epsilon_{N-1}$
and $\Delta\epsilon_{N-1}:=\epsilon_{N-1}-\epsilon_{N-2}$ are independent. All
the steps are presented in Supplementary Information Section
LABEL:SI:Sec:P(track).

Two-state model

Two-state model

$10\boldsymbol{\lambda}$

Three-state model

Four-state model

Six-state model

Figure 9: Comparison of the empirical PDF for two subsequent measured
increments $\Delta y_{1}$ and $\Delta y_{2}$, denoted as $P(\Delta
y_{1},\Delta y_{2})$ (left column) with its up-to-one-switch approximation of
the probability of two noisy subsequent location increments ($P_{1}(\Delta
y_{1},\Delta y_{2})$) (middle column) and with the independent up-to-one-
switch approximation for each increment $P_{1}(\Delta y_{1})P_{1}(\Delta
y_{2})$ (right column). The plots are obtained in a logarithmic color scale,
using the parameters specified in Supplementary Information Figure
LABEL:SI:Fig:networks. $\Delta y_{1}$ is represented on the vertical axis,
while $\Delta y_{2}$ is on the horizontal axis.

We visualise and demonstrate the up-to-one-switch PDF approximation for a two
noisy subsequent location increments, $P_{1}(\Delta y_{1},\Delta y_{2})$, with
the empirical PDF for $\mathbb{P}(\Delta y_{1},\Delta y_{2})$, denoted as
$P(\Delta y_{1},\Delta y_{2})$, however, the result can be applied to any
number of subsequent increments. We also compare $P_{1}(\Delta y_{1},\Delta
y_{2})$ with the naive approximation $P_{1}(\Delta y_{1})P_{1}(\Delta y_{2})$
that assumes that successive increments are approximately independent, to test
the impact of neglecting these correlations. Figure 9 shows that, as expected,
$P_{1}(\Delta y_{1},\Delta y_{2})$ is a better approximation for
$\mathbb{P}(\Delta y_{1},\Delta y_{2})$ than $P_{1}(\Delta y_{1})P_{1}(\Delta
y_{2})$ for low switching rates (infrequent switches). In contrast, the panels
representing the two-state model with $\boldsymbol{\lambda}\times 10$ show no
noticeable difference between $P_{1}(\Delta y_{1},\Delta y_{2})$ and
$P_{1}(\Delta y_{1})P_{1}(\Delta y_{2})$. This suggests that for high
switching frequency the correlation between two subsequent increments is
negligible.

Comparing the two-state model $P(\Delta y_{1},\Delta y_{2})$ (left column)
with $P_{1}(\Delta y_{1},\Delta y_{2})$ (middle column) in Figure 9 we notice
peaks around $\Delta y_{1}=\Delta y_{2}=v_{1}\Delta t$ and $\Delta
y_{1}=\Delta y_{2}=v_{2}\Delta t$. Moreover, the three, four and six-state
models have a peak at zero ($0\cdot\Delta t$), since their networks include
one or more stationary states. These peaks represent a measured interval with
no change in velocity, corresponding to the noisy increment $\Delta y_{1}$,
followed by another measured interval that keeps the same velocity as before,
corresponding to the noisy increment $\Delta y_{2}$. On the other hand, the
right panels representing $P_{1}(\Delta y_{1})P_{1}(\Delta y_{2})$ include
extra peaks which correspond to the probability of having a measured interval
with no change in velocity, corresponding to the noisy increment $\Delta
y_{1}$, followed by another measured interval with no change in velocity,
corresponding to the noisy increment $\Delta y_{2}$, but each with a distinct
velocity and therefore state. Having two subsequent measured intervals with
two different states attained is not a likely scenario, as the state switch
would have to occur exactly at the time between the two measured increments
($t_{1}$) which is a zero-probability event. We also notice that these peaks
are more oval in the left panel while round in the middle (and right) panels,
which is due to the assumption used in the approximation that neglects the
correlation of the noise of subsequent increments (see Equation
(LABEL:Eq:Approx_no_noise_1)).

To summarise, Figure 9 shows that, for low switching rates (infrequent
switches), it is important to take into account the correlation of subsequent
measured increments to approximate $\mathbb{P}(\Delta y_{1},\Delta y_{2})$.
Indeed, $P_{1}(\Delta y_{1},\Delta y_{2})$ is a better approximation for
$\mathbb{P}(\Delta y_{1},\Delta y_{2})$ than $P_{1}(\Delta y_{1})P_{1}(\Delta
y_{2})$ for all model networks presented. Finally, the Supplementary
Information Figure LABEL:SI:Fig:Err_P(Delta_y_1,_Delta_y_2) compares the error
of the joint up-to-one-switch PDF approximation for two noisy subsequent
location increments, defined as $|P_{1}(\Delta y_{1},\Delta y_{2})-P(\Delta
y_{1},\Delta y_{2})|$, with the error of the product of the marginal
approximate PDFs, defined as $|P_{1}(\Delta y_{1})P_{1}(\Delta y_{2})-P(\Delta
y_{1},\Delta y_{2})|$, for the panels shown in Figure 9.

##  5 Discussion and conclusions

In this manuscript, we present a general $n$-state one-dimensional velocity-
jump model and provide approximate solutions to the PDFs of noisy location
increments, motivated by the recent improvement of techniques to collect
individual tracking data. The state evolution is a CTMC, and each state is
associated with a defined velocity, rate of switching, and probabilities of
transitioning to every other state. Single-agent tracking data is
intrinsically noisy and collected at discrete time points, which impedes the
determination of state-switching times, and therefore the direct computation
of the model parameters. Since the approximate PDFs computed could be used as
approximate likelihoods in a calibration framework, we overcome the challenge
of estimating the model parameters using single-agent tracking data.

We compute the approximate PDFs by conditioning on the states attained during
a measured interval, limiting the number to at most two and three states, as
for frequent enough image collection we expect the majority of the increments
to include few state transitions. Hence, we derive an up-to-one-switch and an
up-to-two-switch approximation for the PDF of a single increment, and compare
them to give insights on the parameter regimes at which they are valid.
Comparing these approximations to the empirical PDFs we note that both
approximations agree with the empirical distributions if the switching between
states occurs at sufficiently low frequencies, and they become less precise as
the switching rates increase. As expected, the up-to-two-switch approximation
is more accurate than the up-to-one-switch approximation, especially as the
switching rates increase. We note that the up-to-two-switch approximation was
explicitly computed for states with distinct velocities and rates, however, it
is possible to extend the up-to-two-switch approximation to models in which
different states have the same velocities or rates by distinguishing the
cases, as in the up-to-one-switch approximation. Moreover, the methods
proposed are applicable for any number of switches per interval, but the
possible sets of ordered visited states are combinatorially explosive, which
leads to a high amount of terms to compute.

We also derive an up-to-one-switch approximation for the joint PDF of a set of
subsequent location increments, to extend the results of the marginal PDF
previously obtained, considering the correlations between successive location
increments. Comparing the joint PDF approximation to the empirical PDF and to
the product of the marginal PDFs for each increment shows that the joint PDF
approximation agrees with the empirical PDF, highlighting the importance of
considering the correlation between successive location increments, not
captured by multiplying the marginal PDFs. Again, the methods are valid with
any up-to-$m$-switch marginal PDF approximation.

The approximate PDFs provided can be used to calibrate the $n$-state velocity-
jump model proposed to tracking data. Existing approaches for model
calibration are computationally expensive [40], and based on pseudo-marginal
method [44]. On the converse, our approach provides an analytically tractable
approximation of the likelihood that can be used with any likelihood-based
inference method. Moreover, we provide general results extending the works
that estimate parameters for a two-state Markov model [14, 37], or asymptotic
quantities of interest [13, 21].

The approach taken is, in general, agnostic to the formulation of the noise
model as it is based on convoluting the PDF of an exact increment with the
noise distribution. Thus, the methods can be adapted for use with more general
noise distributions, for which the Markov property used still holds. However,
it may not be trivial to convolute the exact increment distributions with
other noise distributions, but numerical approximations may be used as an
alternative approach. A future challenge is to extend the results obtained for
non-Markovian switching time distributions, or for time-dependent parameters
or self-exciting states.

This work is focused on one-dimensional motion, but it could be extended to
model two-dimensional or three-dimensional motion. In particular, the methods
proposed to approximate the PDFs of exact increments are carried over in
higher dimensions for states with a fixed velocity vector by noting that the
exact location increments become vectors, and the noise is added by
convoluting a multivariate normal distribution. In contrast, having states
with fixed velocity modulus, but varying direction would be a more complex
extension. In this case, one could include angular reorientations following a
state switch, extending the work of [14] with a set of distinct velocities,
but this may lead to complex integrals to solve.

In conclusion, in this paper we approximate the solutions to an $n$-state
velocity-jump model designed to describe single-agent motion in one dimension.
We characterised the PDFs of noisy location increments, overcoming the
challenges posed by the experimental constraints of data collection, which is
noisy and discrete. These PDF approximations have the potential to be used to
calibrate the model proposed to experimental tracking data. Moreover, this
work paves the way for providing solutions and calibration frameworks for
motion models in higher spatial dimensions or with different data
distributions.

*Acknowledgements A.C. is supported by an EPSRC/UKRI Doctoral Training Award. A.P.B. thanks the Mathematical Institute, Oxford for a Hooke Research Fellowship. R.E.B. is supported by a grant from the Simons Foundation (MP-SIP-00001828).

## Declarations

*Conflict of interest The authors declare that they have no conflict of interests.

*Code availability The code is available on GitHub in the repository
a-ceccarelli/distributions_n-state_VJ_model.

*Author contribution A.C., R.E.B. and A.P.B conceived the original idea. A.C. computed the approximations, created the code, carried out the analysis, and wrote the manuscript. A.P.B and R.E.B. gave additional inputs to the manuscript and supervised the work. All authors gave final approval for publication.

## References

  *   * Anderson [2012] Anderson WJ (2012) Continuous-time Markov chains: an applications-oriented approach. Springer Science & Business Media, New York, 10.1007/978-1-4612-3038-0
  * Berg [1993] Berg HC (1993) Random walks in biology. Princeton University Press, Princeton, 10.1515/9781400820023
  * Berg and Brown [1972] Berg HC, Brown DA (1972) Chemotaxis in Escherichia coli analysed by three-dimensional tracking. Nature 239(5374):500–504. 10.1038/239500a0
  * Bergman et al [2000] Bergman CM, Schaefer JA, Luttich S (2000) Caribou movement as a correlated random walk. Oecologia 123(3):364–374. 10.1007/s004420051023
  * Blum and Reed [1985] Blum JJ, Reed MC (1985) A model for fast axonal transport. Cell Motility 5(6):507–527. 10.1002/cm.970050607
  * Blum and Reed [1988] Blum JJ, Reed MC (1988) The transport of organelles in axons. Mathematical Biosciences 90(1-2):233–245. 10.1016/0025-5564(88)90068-5
  * Bovet and Benhamou [1988] Bovet P, Benhamou S (1988) Spatial analysis of animals’ movements using a correlated random walk model. Journal of Theoretical Biology 131(4):419–433. 10.1016/S0022-5193(88)80038-9
  * Bressloff and Newby [2013] Bressloff PC, Newby JM (2013) Stochastic models of intracellular transport. Reviews of Modern Physics 85(1):135. 10.1103/RevModPhys.85.135
  * Cho et al [2020] Cho AY, Roque VRM, Goldman C (2020) Fast and slow axonal transport: A unified approach based on cargo and molecular motor coupled dynamics. Physical Review E 102(3):032410. 10.1103/PhysRevE.102.032410
  * Clancy et al [2011] Clancy BE, Behnke-Parks WM, Andreasson JO, et al (2011) A universal pathway for kinesin stepping. Nature Structural & Molecular Biology 18(9):1020–1027. 10.1038/nsmb.2104
  * Codling and Hill [2005] Codling E, Hill N (2005) Sampling rate effects on measurements of correlated and biased random walks. Journal of Theoretical Biology 233(4):573–588. 10.1016/j.jtbi.2004.11.008
  * Erban and Othmer [2004] Erban R, Othmer HG (2004) From individual to collective behavior in bacterial chemotaxis. SIAM Journal on Applied Mathematics 65(2):361–391. 10.1137/S0036139903433232
  * Han and Fricks [2024] Han L, Fricks J (2024) A semi-Markov approach to study a group of kinesin motors. Bulletin of Mathematical Biology 86(2):15. 10.1007/s11538-023-01241-5
  * Harrison and Baker [2018] Harrison JU, Baker RE (2018) The impact of temporal sampling resolution on parameter inference for biological transport models. PLoS Computational Biology 14(6):e1006235. 10.1371/journal.pcbi.1006235
  * Harrison et al [2019] Harrison JU, Parton RM, Davis I, et al (2019) Testing models of mRNA localization reveals robustness regulated by reducing transport between cells. Biophysical Journal 117(11):2154–2165. 10.1016/j.bpj.2019.10.025
  * Heber et al [2024] Heber S, McClintock MA, Simon B, et al (2024) Tropomyosin 1-I/C coordinates kinesin-1 and dynein motors during oskar mRNA transport. Nature Structural & Molecular Biology pp 1–13. 10.1038/s41594-024-01212-x
  * Hughes et al [2011] Hughes J, Hancock WO, Fricks J (2011) A matrix computational approach to kinesin neck linker extension. Journal of Theoretical Biology 269(1):181–194. 10.1016/j.jtbi.2010.10.005
  * Hughes et al [2012] Hughes J, Hancock WO, Fricks J (2012) Kinesins with extended neck linkers: a chemomechanical model for variable-length stepping. Bulletin of Mathematical Biology 74:1066–1097. 10.1007/s11538-011-9697-6
  * Jones et al [2015] Jones PJ, Sim A, Taylor HB, et al (2015) Inference of random walk models to describe leukocyte migration. Physical Biology 12(6):066001. 10.1088/1478-3975/12/6/066001
  * Kareiva and Shigesada [1983] Kareiva P, Shigesada N (1983) Analyzing insect movement as a correlated random walk. Oecologia 56:234–238. 10.1007/BF00379695
  * Knoops and Vanderzande [2018] Knoops G, Vanderzande C (2018) Motion of kinesin in a viscoelastic medium. Physical Review E 97(5):052408. 10.1103/PhysRevE.97.052408
  * Kutys et al [2010] Kutys ML, Fricks J, Hancock WO (2010) Monte Carlo analysis of neck linker extension in kinesin molecular motors. PLoS Computational Biology 6(11):e1000980. 10.1371/journal.pcbi.1000980
  * Liggett [2010] Liggett TM (2010) Continuous time Markov processes: an introduction, vol 113. American Mathematical Society, Providence, Rhode Island, 10.1090/gsm/113
  * Medvinsky et al [2002] Medvinsky AB, Petrovskii SV, Tikhonova IA, et al (2002) Spatiotemporal complexity of plankton and fish dynamics. SIAM Review 44(3):311–370. doi.org/10.1137/S0036144502404442
  * Meijering et al [2012] Meijering E, Dzyubachyk O, Smal I (2012) Methods for cell and particle tracking. Methods in Enzymology 504:183–200. 10.1016/B978-0-12-391857-4.00009-4
  * Miles et al [2024] Miles CE, McKinley SA, Ding F, et al (2024) Inferring stochastic rates from heterogeneous snapshots of particle positions. Bulletin of Mathematical Biology 86(6):74. 10.1007/s11538-024-01301-4
  * Norris [1998] Norris JR (1998) Markov chains. Cambridge University Press, Cambridge, 10.1017/CBO9780511810633
  * Othmer et al [1988] Othmer HG, Dunbar SR, Alt W (1988) Models of dispersal in biological systems. Journal of Mathematical Biology 26(3):263–298. 10.1007/BF00277392
  * Patel et al [2018] Patel S, Fu S, Mastio J, et al (2018) Unique pattern of neutrophil migration and function during tumor progression. Nature Immunology 19(11):1236–1247. 10.1038/s41590-018-0229-5
  * Pawitan [2001] Pawitan Y (2001) In all likelihood: statistical modelling and inference using likelihood. OUP Oxford, Oxford, 10.1111/1467-9884.00369_20
  * Pike and Burman [2023] Pike TW, Burman OH (2023) Simulating individual movement in fish. Scientific Reports 13(1):14581. 10.1038/s41598-023-40420-1
  * Popovic et al [2011] Popovic L, McKinley SA, Reed MC (2011) A stochastic compartmental model for fast axonal transport. SIAM Journal on Applied Mathematics 71(4):1531–1556. 10.1137/090775385
  * Powalla et al [2022] Powalla D, Hoerner S, Cleynen O, et al (2022) A numerical approach for active fish behaviour modelling with a view toward hydropower plant assessment. Renewable Energy 188:957–966. 10.1016/j.renene.2022.02.064
  * Preisler et al [2004] Preisler HK, Ager AA, Johnson BK, et al (2004) Modeling animal movements using stochastic differential equations. Environmetrics 15(7):643–657. doi.org/10.1002/env.636
  * Ross [2014] Ross SM (2014) Introduction to probability models. Academic press, Cambridge, Massachusetts, 10.1016/C2012-0-03564-8
  * Rosser et al [2013a] Rosser G, Fletcher A, Maini P, et al (2013a) The effect of sampling rate on observed statistics in a correlated random walk. Journal of the Royal Society Interface 10(85):20130273. 10.1098/rsif.2013.0273
  * Rosser et al [2013b] Rosser G, Fletcher AG, Wilkinson DA, et al (2013b) Novel methods for analysing bacterial tracks reveal persistence in rhodobacter sphaeroides. PLoS Computational Biology 9(10):e1003276. 10.1371/journal.pcbi.1003276
  * Rosser et al [2014] Rosser G, Baker RE, Armitage JP, et al (2014) Modelling and analysis of bacterial tracks suggest an active reorientation mechanism in rhodobacter sphaeroides. Journal of The Royal Society Interface 11(97):20140320. 10.1098/rsif.2014.0320
  * Salek et al [2019] Salek MM, Carrara F, Fernandez V, et al (2019) Bacterial chemotaxis in a microfluidic T-maze reveals strong phenotypic heterogeneity in chemotactic sensitivity. Nature Communications 10(1):1877. 10.1038/s41467-019-09521-2
  * Simpson et al [2022] Simpson MJ, Baker RE, Buenzli PR, et al (2022) Reliable and efficient parameter estimation using approximate continuum limit descriptions of stochastic models. Journal of Theoretical Biology 549:111201. 10.1016/j.jtbi.2022.111201
  * Skellam [1991] Skellam J (1991) Random dispersal in theoretical populations. Bulletin of Mathematical Biology 53(1-2):135–165. 10.1007/BF02464427
  * Taylor-King et al [2015] Taylor-King JP, van Loon EE, Rosser G, et al (2015) From birds to bacteria: generalised velocity jump processes with resting states. Bulletin of Mathematical Biology 77:1213–1236. 10.1007/s11538-015-0083-7
  * Treloar et al [2011] Treloar KK, Simpson MJ, McCue SW (2011) Velocity-jump models with crowding effects. Physical Review E 84(6):061920. 10.1103/PhysRevE.84.061920
  * Warne et al [2020] Warne DJ, Baker RE, Simpson MJ (2020) A practical guide to pseudo-marginal methods for computational inference in systems biology. Journal of Theoretical Biology 496:110255. 10.1016/j.jtbi.2020.110255
  * Wu et al [2000] Wu Hi, Li BL, Springer TA, et al (2000) Modelling animal movement as a persistent random walk in two dimensions: expected magnitude of net displacement. Ecological Modelling 132(1-2):115–124. 10.1016/S0304-3800(00)00309-4
  * Xue and Jameson [2017] Xue C, Jameson G (2017) Recent mathematical models of axonal transport. Stochastic Processes, Multiscale Modeling, and Numerical Methods for Computational Cellular Biology pp 265–285. 10.1007/978-3-319-62627-7_12
  * Yates et al [2009] Yates CA, Erban R, Escudero C, et al (2009) Inherent noise can facilitate coherence in collective swarm motion. Proceedings of the National Academy of Sciences 106(14):5464–5469. 10.1073/pnas.0811195106