# A Study of Particle Motion in the Presence of Clusters

Christian Wald11footnotemark: 1 Gabriele Steidl111Institute of Mathematics,
Technische Universität Berlin, Straße des 17. Juni 136, 10623 Berlin, Germany,
{wald<EMAIL_ADDRESS>http://tu.berlin/imageanalysis.

(August 28, 2024)

###### Abstract

The motivation for this study came from the task of analysing the kinetic
behavior of single molecules in a living cell based on Single Molecule
Localization Microscopy. Given measurements of both the motion of clusters and
molecules, the main task consists in detecting if a molecule belongs to a
cluster. While the exact size of the clusters is usually unknown, upper bounds
are available. In this study, we simulate the cluster movement by a Brownian
motion and those of the particles by a Gaussian mixture model with two modes
depending on the position of the particle within or outside a cluster. We
propose various variational models to detect if a particle lies within a
cluster based on the Wasserstein and maximum mean discrepancy distances
between measures. We compare the performance of the proposed models for
simulated data.

##  1 Introduction

In this paper, we are interested in a setting where particles and particle
clusters are moving with different ,,speed” in space, which we assume to be
$\mathbb{R}^{d}$. In practice, $d=2$ or $d=3$ will be of interest. A single
particle might enter a cluster and will then be carried by the cluster, but is
also allowed to leave it at some point. A typical real-world setting is Single
Molecule Localization Microscopy (SMLM) [6] which includes high precision in
localization of both single molecules and clusters. However, the size of the
clusters is usually unknown and only upper bounds exist. The original
motivation to deal with this topic is related to the study the consistency and
morphology of protein microclusters formed in the context of so-called $T$
cell activation, which have been a topic of interest in biology over the last
15 years [9, 8]. In particular, it is not possible to see precisely whether a
single molecule is part of a cluster.

Assuming that the movement patterns of single particles and cluster follow a
Brownian motion with quite different parameters, we propose a simple
simulation of the particle-cluster behavior using a Gaussian mixture model
(GMM) with two modes which address the cluster membership of a particle.
Knowing particle trajectories the parameters of the GMM can be estimated by
the expectation-minimization (EM) algorithm [1].

Based on this model, our main task is to determine whether a particle belongs
to some cluster. We will use a variational approach to determine the
probability that a particle belongs to a cluster at a certain time. More
precisely, we will define different functionals which minimizers serve as
estimators of the above probability. These functionals consists of two
summands, which take i) the distance of the single particle to the closest
cluster and ii) the movement patterns from the GMM into account. For the
second part, we present several choices:

  * -

we relate path trajectories of particles to empirical probability measures and
use distances between probability measures such as the Wasserstein distance
and the maximum mean discrepancy (MMD) to consider their distance to the
cluster component of the GMM.

  * -

we relate the associated empirical probability measures of particle
trajectories and cluster trajectories.

We provide numerical examples on simulated data for both approaches and
different distances between measures. Here we recognized that the latter
approach obtains a better accuracy when comparing with the ground truth of the
simulated data.
Outline of the paper. In Section 2, we explain our model for simulating
particle-cluster motion and recall the EM algorithm for determining the
parameters of the associated two mode GMM if samples of particle trajectories
are available. Then, in Section 3, we introduce the functionals for
determining the probability that a particle belongs to a cluster at a given
time. In Section 4 we compare our different approaches numerically for
simulated data. Conclusions are drawn in Section 5.

##  2 Simulation of single particle versus clusters motion

We start by describing our particle-cluster movement simulation in
$\mathbb{R}^{d}$, where we will rely on $d=2$ in our numerical experiments.
Then we will see that the motion of the particle follows a GMM with two modes
depending on its position inside or outside some cluster. Given (simulated)
measurements of the particle positions at several times, the parameters of
this GMM can be estimated by the EM algorithm.

In the following, we use $\|\cdot\|$ for the Euclidean norm on
$\mathbb{R}^{d}$ and the abbreviation $[N]:=\\{1,\ldots,N\\}$.

###  2.1 Simulation

We consider $N$ clusters at time points $t\in[T]$ modeled as balls with the
same radius $r>0$ and centers $c_{n,t}$, $n\in[N]$, $t\in[T]$. The clusters
move according to a Brownian motion, i.e.,

$\Delta c_{n}(t)\coloneqq
c_{n,t+1}-c_{n,t}\sim\mathcal{N}(m_{c},\sigma^{2}_{c}).$

where we write shorthand $\mathcal{N}(m_{c},\sigma^{2}_{c})$ for
$\mathcal{N}(m_{c},\sigma^{2}_{c}\,I_{d})$. The movement of a single particle
is described by its position $p_{t}\in\mathbb{R}^{d}$, $t\in[T]$. Outside
clusters we assume a Brownian motion

$\Delta p(t)\coloneqq p_{t+1}-p_{t}\sim\mathcal{N}(m_{out},\sigma^{2}_{out}),$

where $\sigma_{out}>\sigma_{c}$. Inside a cluster, the particle is carried by
the cluster with an additional small Brownian motion
$\mathcal{N}(0,\sigma^{2}_{pc})$. This results in an overall Brownian motion
within a cluster given by

$\Delta p(t)\sim\mathcal{N}(m_{in},\sigma^{2}_{in}),\qquad m_{in}\coloneqq
m_{c},\,\sigma^{2}_{in}\coloneqq\sigma^{2}_{c}+\sigma^{2}_{pc}.$

Note that due to the small Brownian motion inside the cluster, a single
particle might also leave the cluster at some time point. In summary, the
cluster and particle motion is simulated by Algorithm 1. An illustration can
be found in figure 1

$\mathbf{Input}$ Initial particle position $p_{0}=0\in\mathbb{R}^{d}$ and
initial cluster positions $c_{0,i}$, $i\in[N]$ uniformly sampled in
$[-b,b]^{d}$ for some $b>0$. Furthermore $r$: cluster radius,
$m_{c},\sigma_{c}^{2}$: cluster motion parameters, $\sigma_{pc}^{2}$: variance
of the small particle motion inside a cluster, $m_{out},\sigma_{out}^{2}$:
particle motion parameters outside clusters

for $t=0,\ldots,T$ do

$r_{t}=\min\limits_{n\in[N]}\|p_{t}-c_{n,t}\|$

$n_{t}=\operatorname*{arg\,min}\limits_{n\in[N]}\|p_{t}-c_{n,t}\|$

if $r_{t}\leq r$ then

$\sigma=\sigma_{pc}$

$m=0$

else

$\sigma=\sigma_{out}$

$m=m_{out}$

$n_{t}=\mathsf{None}$

end if

$p_{t+1}=p_{t}+z_{t}$, where $z_{t}\sim\mathcal{N}(m,\sigma^{2})$

$c_{n,t+1}=c_{n,t}+z_{n,t}$, where
$z_{n,t}\sim\mathcal{N}(m_{c},\sigma^{2}_{c})$

if $\hat{n}\neq\mathsf{None}$ then

$p_{t+1}=p_{t+1}+(c_{n_{t},t+1}-c_{n_{t},t})$

end if

end for

Algorithm 1 Simulation of particle and cluster movement

Figure 1: Trajectory of a particle together with its closest cluster. The
cluster is only displayed if the particle is inside of it. The movement
outside and inside the cluster clearly differs. Algorithm 2 EM algorithms for
determining the GMM parameters in (1).

Let
$f(x,\mu,\sigma)=\frac{1}{2\pi\sigma}\exp\left(-\frac{1}{2}\frac{(x-\mu)^{2}}{\sigma}\right)$
be the likelyhood of a multivariate Gaussian distribution and

$l(\alpha^{r},\mu^{r},\sigma^{r})=\sum_{t=1}^{T-1}\log\left(\sum_{j\in{in,out}}\alpha^{r}_{j}f(\Delta
p(t),\mu^{r}_{j},\sigma^{r}_{j})\right)$ be the log likelihood of the GMM for
$\\{\Delta p(1),\ldots,\Delta p(T-1)\\}$.

$\mathbf{Input}$ $(\Delta p(1),\ldots,\Delta p(T-1))$, Initial values
$\alpha^{0}=(\alpha_{in}^{0},\alpha_{in}^{0})$,
$\mu^{0}=(\mu^{0}_{in},\mu^{0}_{out})$ and
$\sigma^{r}_{=}(\sigma^{0}_{in},\sigma_{out}^{0})$. Stopping tolerance $tol$

while
$l(\alpha^{r},\mu^{r},\sigma^{r})-l(\alpha^{r-1},\mu^{r-1},\sigma^{r-1})>tol$
do

$\mathbf{E}$ step:

for $j\in\\{in,\ out\\}$ do
$\beta_{i,j}^{r}=\frac{\alpha_{j}^{r}f(x_{i},\mu_{j}^{r},\sigma_{j}^{r})}{\underset{j\in\\{in,out\\}}{\sum}\alpha_{j}^{r}f(x_{i},\mu_{j}^{r},\sigma_{j}^{r})}$

end for

$\mathbf{M}$ step

for $j\in\\{in,\ out\\}$ do $\displaystyle\alpha_{j}^{r+1}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\beta_{i,j}^{r};$
$\displaystyle\mu_{j}^{r+1},\sigma_{j}^{r+1}=\underset{\mu_{j},\sigma_{j}}{\text{argmax}}\sum_{i=1}^{n}\beta_{i,j}^{r}\log\left(f(x_{i},\mu_{j},\sigma_{j})\right)$

end for

end while

###  2.2 GMM for particle motion

Using the above simulation, the single particle motion in each time step
follows a GMM with two modes, i.e.,

$\Delta
p(t)\sim\alpha_{out}\mathcal{N}(m_{out},\sigma^{2}_{out})+\alpha_{in}\mathcal{N}(m_{in},\sigma^{2}_{in})$
(1)

where $\alpha_{in},\alpha_{out}\geq 0$ with $\alpha_{in}+\alpha_{out}=1$
denote the probabilities that a particle is inside or outside a cluster at
time $t$. Given measured samples of $p_{t}$, $t\in[T]$ from a certain number
of experiments, we can estimate the parameters
$\alpha_{\iota},\mu_{\iota},\sigma_{\iota}$, $\iota\in\\{in,out\\}$ by the EM
algorithm [1] outlined in Algorithm 2.

##  3 Cluster membership models

Given a series of (simulated) cluster and particle positions $c_{n,t}$ and
$p_{t}$, $n\in[N]$, $t\in[T]$, the main task is to determine whether a
particle lies within a cluster at a certain time $t$. The problem is even more
difficult, since we usually do not know the clusters size $r$, but have
usually only an upper bound $R$.

In this section, we intend to estimate a probability vector
$e\coloneqq(e_{t})_{t\in[T]}\in[0,1]^{T}$ for a particle to be inside some
cluster at time $t$ by solving a minimization problem of the form

$\operatorname*{arg\,min}_{e\in[0,1]^{T}}\
\beta_{r}L_{r}(e)+\beta_{s}L_{s}(e),\qquad\beta_{r},\beta_{s}\geq 0.$

Here $L_{r}$ takes care of the distance of the particle to the closest cluster
and $L_{s}$ depends on the local similarity of the movement of the single
particle to those of the closest cluster. In the next subsections, we will
propose suitable choices for $L_{r}$ and $L_{s}$. We construct $L_{r}$ and
$L_{s}$ such that they are convex in $e$ which enables us to use simple
gradient descent optimization methods. More specifically we use the adaptive
moment estimation optimizer ADAM [5]. Since in our numerical experiments we
obtained convergence there was no need for more sophisticated optimization
algorithms.

###  3.1 Radius function $L_{r}$

In general, we do not know the radius $r$ of a cluster, but can estimate an
upper bound such that $r\leq R$. Let
$\mathcal{C}_{t}\coloneqq\\{c_{1,t},\ldots,c_{N,t}\\}$ denote the set of
position of all clusters at time $t$. We define an rough estimate for $e_{t}$
by

$\bar{e}_{t}\coloneqq\begin{cases}1\text{ if
}\text{dist}(p_{t},\mathcal{C}_{t})\leq R,\\\ 0\text{ otherwise}.\end{cases}$

Then we set

$L_{r}(e)\coloneqq\sum_{t=1}^{T}(e_{t}-\bar{e}_{t})^{2}F\left(\text{dist}(p_{t},\mathcal{C}_{t})\right),$

where $F$ is the following function which encodes the confidence in our
estimation $\bar{e}$:

$\displaystyle F(s)\coloneqq\frac{f(s-R)}{f(s-R)+f\left(1-(s-R)\right)}$ (2)

and

$f(s)\coloneqq\begin{cases}e^{-\frac{1}{s}}\text{ for }s>0,\\\ 0\text{
otherwise}.\end{cases}$

The function $F$ has the following desirable properties:

  * i)

$F(s)=0$ for $s\leq R$ where we don’t know if the particle is inside a cluster
and $F(s)=1$ for $s\geq R+1$ where we can be certain the particle not inside a
cluster,

  * ii)

$F$ is continuous and monotone increasing.

The function is plotted in Figure 2 for $R=0.7$. Since
$(e_{t}-\bar{e}_{t})^{2}$ is convex in $e_{t}$ and
$F\left(\text{dist}(p_{t},\mathcal{C}_{t})\right)$ is a constant also
$L_{r}(e)$ is convex in $e$.

###  3.2 Similarity functions $L_{s}$

Next we propose two different kinds of local similarity functions. To this
end, we consider trajectories of particles as empirical probability measures.
More precisely, we define for fixed $k\in\mathbb{N}$, the sliding window
measure of the particles at time $t\in[T]$ by

$\mu_{p,t}\coloneqq\frac{1}{k_{t}+K_{t}+1}\sum_{i=k_{t}}^{K_{t}}\delta_{p_{t+i+1}-p_{t+i}},$
(3)

where $k_{t}\coloneqq\max\\{-k,1-t\\}$ and
$K_{t}\coloneqq\min\\{k-1,T-t-1\\}$. In our numerical experiments $k$ will be
chosen as $k=6$. Then $\mu_{p,t}$ can be considered as kind of local particle
movement. Further, we suppose that the parameters $m_{\iota},\sigma_{\iota}$,
$\iota\in\\{in,out\\}$ of the GMM (1) were estimated by the EM algorithm.

Let
$\text{d}:\mathcal{P}(\mathbb{R}^{d})\times\mathcal{P}(\mathbb{R}^{d})\to\mathbb{R}_{\geq
0}$ denote a distance function between probability measures. In this paper, we
will compare settings using the Wasserstein-2 distance and the MMD distance
which are explained in Appendix A.

We model the similarity of $\mu_{t}$ to the movement within clusters in two
different ways. The first one takes the GMM parameters $m_{in},\sigma_{in}$
into account, while the second one uses a sliding window measure of the
nearest cluster centers.

#### 3.2.1 Similarity between $\mu_{p,t}$ and
$\mathcal{N}(\mu_{in},\sigma^{2}_{in})$

For a threshold $h$ which we will determine below, we introduce a rough
estimate of $e_{t}$ by

$\displaystyle\hat{e}_{t}\coloneqq\begin{cases}1\ \text{if}\
\text{d}(\mu_{p,t},\mathcal{N}(m_{in},\sigma^{2}_{in}))\leq h,\\\ 0\
\text{otherwise}.\end{cases}$ (4)

Then we define

$L_{s}(e)\coloneqq\sum_{t=k+1}^{T-k}(e_{t}-\hat{e}_{t})^{2}G\left(\text{d}\left(\mu_{p,t},\mathcal{N}(m_{in},\sigma_{in})\right)\right),$
(5)

where

$G(s)\coloneqq\min\\{1,g(s)\\}$ (6)

and

$g(s)\coloneqq\frac{1}{h}^{2}s^{2}-\frac{2}{h}s+1.$

Here $g$ is a weight function that encodes the confidence of the estimates
$\hat{e}_{t}$. The confidence should be small if
$\text{d}(\mu_{p,t},\mathcal{N}(m_{in},\sigma_{in}))$ is close to the
threshold. This is true for $g$ since we have, see also figure 2,

  * i)

$0\leq g(s)\leq 1$ for all $s\in\mathbb{R}$

  * i)

$g(s)=0$ for $s=h$

  * ii)

$g(0)=1$ i.e. if $\text{d}(\mu_{p,t},\mathcal{N}(m_{in},\sigma_{in}))=0$ which
means the single particle moves exactly like the clusters

  * iii)

$g(2h)=1$ i.e. we are confident that the particle is not inside a cluster when
its movement differs significantly from that of clusters.

$-0.5$$0.5$$1$$1.5$$2$$-0.5$$0.5$$1$$1.5$$x$$y$

$-0.5$$0.5$$1$$1.5$$2$$-0.5$$0.5$$1$$1.5$$x$$y$

Figure 2: Left: Function $F$ from (2) for $R=0.7$. Right: Function $G$ from
(6) for a threshold of $h=0.5$

Next, we propose a choice for the threshold $h$. For independent random
vectors $X^{\iota}_{-k},\ldots,X^{\iota}_{k-1}$, $\iota\in\\{in,out\\}$
distributed as $\mathcal{N}(m_{\iota},\sigma^{2}_{\iota})$, set
$\mathbf{X}^{\iota}\coloneqq(X^{\iota}_{-k},\ldots,X^{\iota}_{k-1})$, and
define random measures by

$\bm{\mu}_{\mathbf{X}^{\iota}}\coloneqq\frac{1}{2k}\sum_{i=-k}^{k-1}\delta_{X_{i}^{\iota}},\qquad\iota\in\\{in,out\\}.$

Let the random variables $X^{\iota}_{l}$ be defined on a measure space
$\Omega,\Sigma,\alpha$. Then
$\Omega:\omega\to\text{d}\left(\bm{\mu}_{\mathbf{X}^{\iota}},\mathcal{N}(m_{in},\sigma^{2}_{in})\right)$
is a random variable. With

$\displaystyle E_{\iota}$
$\displaystyle\coloneqq\mathbb{E}_{\mathbf{X}^{\iota}}\left[\text{d}\left(\bm{\mu}_{\mathbf{X}^{\iota}},\mathcal{N}(m_{in},\sigma^{2}_{in})\right)\right]=\int_{\Omega}\text{d}\left(\bm{\mu}_{\mathbf{X}^{\iota}},\mathcal{N}(m_{in},\sigma^{2}_{in})\right)\mathrm{d}\alpha,$
(7)

we suggest as threshold

$h\coloneqq\frac{E_{in}+E_{out}}{2}.$

In the numerical experiments we always use the empirical expected value and
variance.

#### 3.2.2 Similarity between $\mu_{p,t}$ and closest cluster motion

We can also directly compare the local particle movement to the movement of
the closest cluster. To this end, let

$n_{t}\coloneqq\operatorname*{arg\,min}_{n\in[N]}\ \|c_{n,t}-p_{t}\|$

be the index of the closest cluster at time point $t$. For a window
$(t-k,\cdots,t+k-1),$ let $\hat{n}_{t}$ be the smallest index which occurs
most in $(n_{t-k},\ldots,n_{t+k-1})$. Then we define the empirical measure

$\mu_{c,t}\coloneqq\frac{1}{k_{t}+K_{t}+1}\sum_{i=k_{t}}^{K_{t}}\delta_{c_{\hat{n}_{t},t+i+1}-c_{\hat{n}_{t},t+i}},$

where $k_{t}$ and $K_{t}$ are given as in (3). This describes the movement of
the closest cluster at times $(t-k,\ldots,t-k+1)$. For a threshold $h$ which
we will determine below, we introduce a rough estimate of $e_{t}$ by

$\hat{e}_{t}\coloneqq\begin{cases}1\ \text{if}\
\text{d}(\mu_{p,t},\mu_{c,t})\leq h,\\\ 0\ \text{else}.\end{cases}$

Then we define

$L_{s}(e)\coloneqq\sum_{t=k+1}^{T-k}(e_{t}-\hat{e}_{t})^{2}G\left(\text{d}(\mu_{p,t},\mu_{c,t})\right),$

where the confidence function $G$ can be chosen as in (6).

Finally, we suggest a threshold $h$. We assume that we have computed an
estimate $m_{c},\sigma_{c}$ for the cluster movement, e.g. by a simple maximum
likelihood estimation. Let $X_{-k},\ldots,X_{k-1}$ be independent random
vectors distributed as $\mathcal{N}(0,\sigma^{2}_{in}-\sigma^{2}_{c})$ and
$Y_{-k},\ldots,Y_{k-1}$ be independent random vectors distributed as
$\mathcal{N}(m_{c},\sigma^{2}_{c})$. Set $Z_{i}\coloneqq Y_{i}+X_{i}$,
$i=-k,\ldots,k-1$, i.e. $Z_{i}$ models a single particle motion that belongs
to this cluster. Then we consider

$\displaystyle E_{in}$
$\displaystyle\coloneqq\mathbb{E}_{\mathbf{X},\mathbf{Y}}\left[\text{d}\left(\bm{\mu}_{\mathbf{Y}},\bm{\mu}_{\mathbf{Z}}\right)\right].$
(8)

Furthermore, let $W_{-k},\ldots,W_{k-1}$ be distributed as
$\mathcal{N}(m_{out},\sigma^{2}_{out})$ and

$\displaystyle E_{out}$
$\displaystyle\coloneqq\mathbb{E}_{\mathbf{X},\mathbf{W}}\left[\text{d}\left(\bm{\mu}_{\mathbf{Y}},\bm{\mu}_{\mathbf{W}}\right)\right],.$
(9)

Similarly as before, we define

$h\coloneqq\frac{E_{in}+E_{out}}{2}.$

######  Remark 1.

We used the estimates for the GMM in order to determine the threshold $h$. If
we are in a situation, where a GMM is not available, we can also use a
clustering method on the set
$\\{\text{d}(\mu_{p,t},\mu_{c,t}):t\in\\{1,\ldots,T\\}\\}$ to determine $h$.

##  4 Numerical results

All simulations where done in the Julia programming language. For the
Wasserstein metric we used the PythonOT for Julia package.

###  4.1 EM algorithm

Table 1 gives the result of the EM algorithm on simulated data for variuos
configurations. It shows that bigger differences in the simulation parameters
$\mu_{in},\mu_{out}$ and $\sigma_{in},\sigma_{out}$ improves the estimation.
If we increase the time step the estimation also improves. The key question is
if the estimations suffice in order to get good estimates whether a molecule
is inside a cluster in the next section.

Simulation Parameters | $\sigma_{in}$ | $\sigma_{out}$ | $\mu_{in}$ | $\mu_{out}$ | $\alpha_{in}$ | $\alpha_{out}$
---|---|---|---|---|---|---
T: $1000$, r: $0.7$, runs: $10$
Simulation parameters | 0.5 | 0.7 | (0.3,0.3) | 0 | 0.54 | 0.46
EM estimated | 0.48 | 0.71 | (0.32, 0.29) | (-0.01,-0.01) | 0.54 | 0.46
Simulation parameters | 0.5 | 1 | (0.3,0.3) | 0 | 0.58 | 0.42
EM estimated | 0.50 | 1.00 | (0.31, 0.30) | (-0.02,-0.04) | 0.58 | 0.42
Simulation parameters | 0.5 | 0.7 | 0 | 0 | 0.56 | 0.44
EM estimated | 0.45 | 0.69 | (-0.02, -0.04) | (0.24,-0.09) | 0.55 | 0.45
T: $5000$, r: $0.7$, runs: $10$
Simulation parameters | 0.5 | 0.7 | (0.3,0.3) | 0 | 0.55 | 0.45
EM estimated | 0.49 | 0.69 | (0.31, 0.30) | (0.02,0.01) | 0.52 | 0.48
Simulation parameters | 0.5 | 0.7 | 0 | 0 | 0.57 | 0.34
EM estimated | 0.51 | 0.73 | (0.00,0.00) | (-0.01,0.01) | 0.61 | 0.39
Table 1: EM estimations of $\sigma$ and $\mu$

###  4.2 Cluster membership

For every run of our experiments we simulated data according to algorithm 1
and estimated the parameters of the corresponding GMM with the EM algorithm.
Using the cluster membership models from above we obtain estimates
$e=(e_{k+1},\ldots,e_{T-k})$. Note that for $t\leq k,t>T-k$ the measure
$\mu_{p,t}$ has fewer than $2k$ points why we omitted them in our analysis.
Then let $\hat{e}=(\hat{e}_{k+1},\ldots,\hat{e}_{T-k})$ for

$\hat{e}_{t}=\begin{cases}0\text{ if }\ e_{t}<\frac{1}{2}\\\ 1\text{ else
}\end{cases}$

We view $\hat{e}$ as the solution of a classification problem and report the
accuracy i.e. the ratio of correctly classified time points. We show the mean
accuracy over the runs in table 2. Detail about the distance functions on
$\mathcal{P}_{1}(\mathbb{R}^{d})$ that we used can be found in appendix A.

#### Using the similarity between $\mu_{p,t}$ and
$\mathcal{N}(\mu_{in},\sigma^{2}_{in})$

The first column of table 2 displays the accuracy of the estimation whether
the particle is inside a cluster at a given time point. If we only use the
radius term $L_{r}$, i.e., $\beta_{r}=0$ the peformance is much worse than
using only the term $L_{s}$. Combining $L_{r}$ and $Ls$ is beneficial for all
choices of $L_{s}$. Generally all choices of $L_{s}$ perform similarly with
the Wasserstein distance based beeing sligthly better than the others. All
method perform better when the movement pattern of the single molecule differ
stronger from the ones of clusters.

$T=1000$, $r=0.7$, $R=1.2$, runs: $20$, $k=6$
---
$\sigma_{in}=0.5$, $\sigma_{out}=0.7$, $\mu_{in}=(0.3,0.3)$, $\mu_{out}=0$
|  |  $\mu_{p,t}$ and $\mathcal{N}(\mu_{in},\sigma^{2}_{in})$ |  $\mu_{p,t}$ and $\mu_{c,t}$
|  | Mean Accuracy | Mean Accuracy
$\beta_{r}$ | $\beta_{s}$ | MV | WS | MMD | MV | WS | MMD
$0$ | $1$ | 0.77 | 0.77 | 0.74 | 0.73 | 0.90 | 0.90
$1$ | $1$ | 0.84 | 0.89 | 0.87 | 0.81 | 0.90 | 0.94
$1$ | $0$ | 0.70 | 0.70 | 0.68 | 0.70 | 0.70 | 0.69
$\sigma_{in}=0.5$, $\sigma_{out}=1.0$, $\mu_{in}=(0.3,0.3)$, $\mu_{out}=0$
$0$ | $1$ | 0.83 | 0.86 | 0.79 | 0.72 | 0.90 | 0.92
$1$ | $1$ | 0.87 | 0.92 | 0.90 | 0.82 | 0.91 | 0.93
$1$ | $0$ | 0.68 | 0.67 | 0.68 | 0.68 | 0.69 | 0.70
Table 2: Accuracy of the cluster membership models. Second column: membership
model using the similarity between $\mu_{p,t}$ and
$\mathcal{N}(\mu_{in},\sigma^{2}_{in})$, see Section 3.2.1. Third column:
Accuracy of the cluster membership model using the similarity between
$\mu_{p,t}$ and the closest cluster motion, see Section 3.2.2. MV: Variance
and Expected Value distance, WS: Wasserstein-2 Distance, see appendix A

#### Using the similarity between $\mu_{p,t}$ and closest cluster motion.

As demonstrated in the third row of table 2 using the movement of the closest
cluster has a better accuracy than comparing with
$\mathcal{N}(\mu_{in},\sigma^{2}_{in})$ for all methods except for the Mean
and Variance method. This is especially evindent when we compare the
performance without the radius function $L_{r}$. For the Wasserstein and MMD
based method adding $L_{r}$ does not improve the accuracy by much possibly
because the room for improvement is low. Note that there is an inherent error
that we have to expect a time point $t$ when a particle enters a cluster but
is outside for all time points $t-l$ for $l\leq k$. This could also explain
why there is no improvement in accuracy when changing the movement of the
single particle from $\mathcal{N}(0,0.7)$ to $\mathcal{N}(0,1.0)$.

##  5 Conclusions

We studied how simulated data from a two mode GMM can be used to to determine
if a particle belongs to a certain cluster which itself moves according to a
Brownian motion. Numerical examples give evidence that the approach is suited
to solve the task. Our model can be generalized at several parts. The cluster
radius may vary between different clusters and maybe also differ over the
time. It may be interesting to consider other distributions instead of the
Gaussian ones, as e.g. heavy tailed distributions. For an EM algorithm for
Student-t distributions, see, e.g. [4]. Since our cluster membership model
that exploits the similarity of $\mu_{p,t}$ and $\mu_{c,t}$ can be also used
by determining the threshold via clustering even settings where the
distributions are not known can be tackled.

Our method shares some similarities with changepoint detection models which
use e.g. Wasserstein distance or MMD in that it is a sliding window approach
that compares local behavior by distances on the space of probability
measures, see e.g. [2, 11, 7]. However they compare a window before a certain
time point with a window after that time point whereas we compare a window
around a time point to either a fixed distribution or to the local behavior of
a cluster. A question is, if there is a practical relevance for such
examinations. Clearly, we want to apply our cluster membership models to real-
word measurements of tracked molecules and clusters as indicated in the
motivation for this work. The setup for such measurements as well as a
reliable molecule position tracking algorithm to test our model is however a
hard task on its own.

### Acknowledgement

Funding by the DFG within the SFB “Tomography Across the Scales” (STE
571/19-1, project number: 495365311) is gratefully acknowledged. Many thanks
to J. Hertrich for fruitful discussions.

## References

  * [1] A. P. Dempster, N. M. Laird, and D. B. Rubin.  Maximum likelihood from incomplete data via the em algorithm.  Journal of the Royal Statistical Society: series B (methodological), 39(1):1–22, 1977.
  * [2] K. Faber, R. Corizzo, B. Sniezynski, M. Baron, and N. Japkowicz.  Watch: Wasserstein change point detection for high-dimensional time series data.  In 2021 IEEE International Conference on Big Data (Big Data), pages 4450–4459. IEEE, 2021.
  * [3] A. Gretton, K. Borgwardt, M. Rasch, B. Schölkopf, and A. Smola.  A kernel method for the two-sample-problem.  Advances in Neural Information Processing Systems, 19, 2006.
  * [4] M. Hasanasab, J. Hertrich, F. Laus, and G. Steidl.  Alternatives of the EM algorithm for estimating the parameters of the Student-t distribution.  Numerical Algorithms, 87:77–118, 2021.
  * [5] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [6] M. Lelek, M. T. Gyparaki, G. Beliu, F. Schueder, J. Griffié, S. Manley, R. Jungmann, M. Sauer, M. Lakadamyali, and C. Zimmer.  Single-molecule localization microscopy.  Nature Reviews Methods Primers, 1(1):39, 2021.
  * [7] S. Li, Y. Xie, H. Dai, and L. Song.  Scan B-statistic for kernel change-point detection.  Sequential Analysis, 38(4):503–544, 2019.
  * [8] B. Rossboth, A. M. Arnold, H. Ta, R. Platzer, F. Kellner, J. B. Huppa, M. Brameshuber, F. Baumgart, and G. J. Schütz.  Tcrs are randomly distributed on the plasma membrane of resting antigen-experienced t cells.  Nature Immunology, 19(8):821–827, 2018.
  * [9] E. Sherman, V. Barr, and L. E. Samelson.  Super-resolution characterization of tcr-dependent signaling clusters.  Immunological Reviews, 251(1):21–35, 2013.
  * [10] C. Villani et al.  Optimal transport: old and new, volume 338.  Springer, 2009.
  * [11] S. Wei and Y. Xie.  Online kernel CUSUM for change-point detection.  arXiv preprint arXiv:2211.15070, 2022.

##  Appendix A Distances between probability measures

#### Wasserstein Distance.

A popular distance that emerged from optimal transport is the Wasserstein-$p$
distance for $p\geq 1$ which is defined as follows see e.g. [10] chapter 6.
Let $\mu,\nu\in\mathcal{P}_{p}(\mathbb{R}^{d})$ measures of finite $p$ moment
i.e. $\int_{\mathbb{R}^{d}}\|x\|^{p}\mathrm{d}\mu(x)<\infty$. Let
$\Gamma(\mu,\nu)\subset\mathcal{P}(\mathbb{R}^{d}\times\mathbb{R}^{d})$ be the
set of $\alpha\in\mathcal{P}(\mathbb{R}^{d}\times\mathbb{R}^{d})$ s.t.
$\pi^{1}_{\sharp}\alpha=\mu$ and $\pi^{2}_{\sharp}\alpha=\nu$ where $\pi^{i}$
are the projections onto the $i$ component. Then the Wasserstein-$p$ distance
is defined as

$\mathcal{W}_{p}(\mu,\nu)=\underset{\alpha\in\Gamma(\mu,\nu)}{\mathsf{min}}\
\left(\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\|x-y\|^{p}\mathsf{d}\alpha(x,y)\right)^{\frac{1}{p}}.$

$\mathcal{W}_{p}$ is a metric on the space of probability measures
$\mathcal{P}_{p}(X)$ with finite $p$-th moment. Note that finite empirical
measures which are the only ones relevant for our application have finite $p$
moment for all $p\geq 1$.

#### Maximum Mean Discrepancy (MMD).

Another important distance on $\mathcal{P}(\mathbb{R}^{d})$ is the MMD [3].
Let $K:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$ be a kernel. Then for
two probability measures $\mu,\nu$ we define

$\displaystyle\mathsf{MMD}^{2}(\mu,\nu)=$
$\displaystyle\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}K(x,y)\mathrm{d}\mu(x)\mathrm{d}\mu(y)-2\int_{R^{d}\times\mathbb{R}^{d}}K(x,y)\mathrm{d}\mu(x)\mathrm{d}\nu(y)$
$\displaystyle+\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}K(x,y)\mathrm{d}\nu(x)\mathrm{d}\nu(y).$

We will use as kernel the Riesz kernel $K(x,y)=-\|x-y\|$ for which the MMD
defines a metric on $\mathcal{P}_{1}(\mathbb{R}^{d})$.

#### Mean-Variance (MV)

Finally, we have also used a simple similarity function based on the variance
and expected value of probability measures. For $\mu,\nu$, we compare

$\displaystyle\mathbb{E}_{\mu},\mathbb{E}_{\nu}$
$\displaystyle\mathsf{cov}(\mu),\mathsf{cov}(\nu)$

where $\mathbb{E}_{\mu}=\int_{\mathbb{R}^{d}}x\mathrm{d}\mu\in\mathbb{R}^{d}$
and
$\mathsf{cov}(\mu)=\int_{\mathbb{R}^{d}}(x-\mathbb{E}_{\mu})(x-\mathbb{E}_{\mu})^{T}\mathrm{d}\mu$.
If we assume the covariance matrix to be $\sigma\mathsf{Id}_{d}$ for
$\sigma\in\mathbb{R}$ we also can look at

$\sqrt{\text{det}(\text{cov}(\mu))}.$

We define

$\displaystyle d_{m}(\mu,\nu)$
$\displaystyle=\|\mathbb{E}_{\mu}[x]-\mathbb{E}_{\nu}[x]\|_{2}$ (10)
$\displaystyle d_{v}(\mu,\nu)$
$\displaystyle=|\sqrt{\mathsf{det}(\mathsf{cov}(\mu))}-\sqrt{\mathsf{det}(\mathsf{cov}(\nu))}|.$
(11)

Here we have two functions $L_{s}$ comparing single molecule and cluster
movement denoted by $L_{s,mean}$ and $L_{s,var}$ coming from $d_{m}$ and
$d_{v}$, see (10). Suppose we estimated
$\hat{\mu}_{in},\hat{\sigma}_{in},\hat{\mu}_{out},\hat{\sigma}_{out}$ via EM.
We then define the threshold $h$ in (4) for $L_{s,mean}$ by
$h=\frac{\|\hat{\mu}_{in}-\hat{\mu}_{out}\|_{2}}{2}$. For $L_{s,var}$ we
define $h$ by $h=\frac{\hat{\sigma}^{2}_{in}+\hat{\sigma}^{2}_{out}}{2}$. We
then define the $L_{s}$ in the mean-variance case by

$L_{s}=L_{s,mean}+L_{s,var}.$