# ToMATo: an efficient and robust clustering algorithm for high dimensional datasets. An illustration with spike sorting

Louise Martineau, Christophe Pouzat, Ségolen Geffray

► To cite this version:

Louise Martineau, Christophe Pouzat, Ségolen Geffray. ToMATo: an efficient and robust clustering algorithm for high dimensional datasets. An illustration with spike sorting. 2025. hal-05267380

HAL Id: hal-05267380

https://hal.science/hal-05267380v1

Preprint submitted on 18 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ToMATo: an efficient and robust clustering algorithm for high dimensional datasets. An illustration with spike sorting.

Louise Martineau, Christophe Pouzat, Ségolen Geffray

September 18, 2025

## Abstract

Clustering algorithms became an essential part of the neurophysiological data analysis toolbox in the last twenty five years. Many problems, from the definition of cell types/groups based on morphological, molecular and physiological data to the identification of sub-networks in fMRI data, are now routinely tackled with clustering analysis. Since the datasets to which this type of analysis is applied tend to be defined in larger and larger dimensional spaces, there is a need for efficient and robust clustering methods in high dimension. There is also a need for methods that assume as little as possible about the clusters shape and size. We report here our experience with the ToMATo (Topological Mode Analysis Tool) algorithm. It is based on a definitely deep mathematical theory (algebraic topology), but its Python based open-source implementation is easily accessible to practitioners. We applied ToMATo to a problem we know well, spike sorting. Its capability to work in the “native” space of the data (no dimension reduction is required) is remarkable, as well as its robustness with respect to outliers (superposed spikes).

## 1. Introduction

Clustering method [13, Chap. 14] applications to neurobiological data have undergone a spectacular growth in recent years. To cite just a few examples: spike sorting [20, 14]; grouping voxels with similar activities from fMRI data [1]; identifying neural activity patterns from extracellular recordings [2].
Recording and decoding the activity of multiple neurons is a major subject in con-
temporary neuroscience. Extracellular recordings with multi-electrode arrays is one of
the basic tools used to that end. The raw data produced by these recordings are almost
systematically a mixture of activities from several neurons. In order to find the number
of neurons which contributed to the recording and identify which neuron generated each
of the visible spikes, a pre-processing step called spike sorting is required.

Spike sorting is nowadays a semi-automatic process which involves many steps. In-
deed, following some initial steps (data normalization, spike detection, event construc-
tion), spike sorting boils down to a clustering problem in high dimension. It is therefore
accompanied most of the time by a dimension reduction step. This dimension reduc-
tion step is sensitive to the presence of superpositions, that is, the superposition of the
activity of two or more neurons that fired nearly simultaneously. These superpositions,
akin to outliers, lead to poor clustering results. Neuroscientists are then usually led
to perform an extra pre-processing step to remove these superpositions. This step is
not completely automated and does a much better job when supervised by an expert in
spike sorting. In addition, many clustering methods exhibit a serious performance re-
duction with increasing dimension and the current trend is clearly towards data defined
in larger and larger dimensional spaces (in spike sorting, but also for instance by con-
sidering longer time series). Since clustering methods are typically used when little is a
priori known about the data, it makes sense to use methods making very few hypothesis
about the data, that is nonparametric methods. With these two requirements in mind
we explored the performances of the ToMATo clustering method. An experienced reader
might rightfully be surprised hearing about a method that is both nonparametric and
able to cope with large dimensions. That is why we wrote that the data are defined in
a space of large dimension, making a distinction between the vector length (dimension)
defining each data point and the intrinsic dimension in which the data are living. If the
data have a large intrinsic dimension, ToMATo may fail, but if this intrinsic dimension
is small (of the order of 10), it performs very well (in our experience), even if the native
space of the data has large dimension (in our example 180).

The use of the ToMATo clustering algorithm helps to simplify and streamline this part
of the spike sorting procedure.

ToMATo (Topological Mode Analysis Tool) is a clustering method using persistent
homology, developed in 2013 by Chazal, Guibas, Oudot and Skraba [7]. It seems little
known in neuroscience but we show in this article how it can be extremely helpful and
effective using as an example an application to spike sorting. Indeed it enables to reduce
the number of steps typically involved in this procedure while providing valuable results.
Also, four main reasons to use the ToMATo algorithm for spike sorting are the following:

1. ToMATo works without dimension reduction as a prior step.

2. ToMATo is robust to superpositions.
3. ToMATo provides a very easy way of choosing the right number of clusters, solving thereby a significant problem in clustering.

4. ToMATo runs fast and requires tuning only a few parameters.

Very importantly, the ToMATo algorithm is implemented and well-documented in Gudhi<fnref n="1" />, a generic open source C++ library for topological data analysis, with a Python interface. It is thus very easy to use even for non specialists of persistent homology.

In section 2, we present the ToMATo algorithm and explain how to use it in practice. This section is very detailed since we are convinced that the user of a method should have a clear understanding of how it is designed and how it works (ToMATo is not a black-box!). In section 3, we show applications of this algorithm to spike sorting on simulated and real data, demonstrating the spectacular performances of this approach. Appendix A provides a high level description of a complete spike sorting procedure, in order to help the reader to see where ToMATo brings key improvements.

# 2. The ToMATo algorithm

## 2.1. Introduction

The ToMATo (Topological Mode Analysis Tool) algorithm is a mode-seeking algorithm. The general idea of mode-seeking algorithms for clustering consists in seeking peaks in the observation density $f$, and in assigning observations falling under the same peak to the same cluster. Indeed, if points are sampled under $f$, there should be a cluster of points corresponding to each peak of $f$. This type of algorithms is supposed to be able to find clusters of any shape, as opposed to algorithms such as kmeans which work well only on convex clusters. A major problem classically arises in mode-seeking: it can be very sensitive to small perturbations of $f$. It turns out that in practice we have only access to an approximation $\hat{f}$ of the true density $f$ and that the peaks of $\hat{f}$ do not in general coincide with the ones of $f$.

Several strategies can be considered to address this issue. The innovative approach of ToMATo resides in the use of persistent homology theory, and thanks to persistent homology a notion of peak prominence is introduced, such that prominent peaks of $\hat{f}$ correspond to prominent peaks of $f$. Clusters found by mode-seeking are merged together so that the final clusters correspond only to prominent peaks of the true density $f$, and not to some spurious, noise induced, peaks of $\hat{f}$.

More precisely, ToMATo is based on a graph mode-seeking algorithm, introduced in [15]. It only uses a density estimate at the data points, and performs mode-seeking on an auxiliary structure, a neighbourhood graph, instead of performing mode-seeking

<footnote n="1">¹https://gudhi.inria.fr/python/latest/clustering.html.</footnote>
directly on a density estimate. This graph-based mode-seeking method allows for making effective computations. It still suffers from the same drawback as any mode-seeking method: it is likely to find too many clusters by taking into account noise induced peaks. To recover some stability, ToMATo therefore combines this computationally effective graph-based mode-seeking, with a merging step based on peak prominence.

## 2.2. Ideas behind ToMATo in the continuous setting

To give insight and understand the ideas behind ToMATo, we first give an overview of what happens in the continuous setting, that is, in the theoretical case where we work directly on (smooth) functions. We start with some reminders on mode-seeking, and persistent homology. In the following we denote a generic function by $g$ and we reserve the notation $f$ for functions that are densities.

### 2.2.1. Gradient ascent for mode-seeking

Let us recall the principle of gradient ascent mode-seeking. We highlight the fact that in practice in the discrete case, we won't need to estimate density gradients. We talk about gradient ascent here in the continuous setting, to mathematically define the intuitive notions of peaks and their ascending regions.

Intuitively a cluster should be a group of data points being under the same peak of the (unknown) density $f$. In practice, a peak can be defined as a local maximum of $\hat{f}$, which can be identified by gradient ascent. We define the ascending region of a peak, as the set of data points that converge to this peak when used as starting points of the gradient ascent procedure. The clusters will then be defined as the ascending regions of $\hat{f}$. Ascending regions is an unstable quantity, as illustrated in Figure 1.

### 2.2.2. Persistent homology of functions

Persistent homology of functions is a method rooted in both Morse theory and topological data analysis. At the core of persistent homology of a generic function $g$, is the evolution of the connected components<fnref n="2" /> of the superlevel sets of $g$, the definition of which is recalled below.

**Definition 2.1.** Superlevel set

Let $g$ be a function from $\mathbb{R}^d$ to $\mathbb{R}$ for $d \ge 1$. Let $\alpha \in \mathbb{R}$. The superlevel set of $g$ of parameter $\alpha$ is:

$$ \{x | g(x) \ge \alpha\}. $$

We keep track of the connected components of $\{x | g(x) \ge \alpha\}$, as a parameter $\alpha$ decreases from $+\infty$ to $-\infty$. We can imagine that $g$ represents an altitude and that $\alpha$

<footnote n="2">²We can also be interested in the evolution of 1-dimensional holes (the void inside a circle or a triangle), 2-dimensional holes (the void inside a sphere or a tetrahedron) and so on. For clustering, we only need to look at connected components, which can be seen as 0-dimensional holes.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: A: A density function $f$ with two peaks $p_1$ and $p_2$. Their respective ascending regions are denoted by $\mathcal{A}(p_1)$ and $\mathcal{A}(p_2)$. B: An approximation $\hat{f}$ of $f$ based on a sample is represented in red. This approximation has many peaks, only a few are highlighted together with their ascending regions.</figcaption></figure>

represents the sea level, then the connected components correspond to the surface (cross-sections) of the islands that appear as the sea level decreases. The nested sequence of sets $(\{x | g(x) \ge \alpha\})_{+\infty \ge \alpha \ge -\infty}$ is called a **filtration**. Let us take the function in Figure 2 A as an example. When $\alpha > g(p_1)$, $\{x | g(x) \ge \alpha\}$ is empty. Then, a connected component $C_1$ appears at $\alpha = g(p_1)$, the global maximum of $g$. As $\alpha$ decreases, $C_1$ grows but remains the only connected component of $\{x | g(x) \ge \alpha\}$, until $\alpha = g(p_2)$ when another connected component $C_2$ appears. At $\alpha = g(v)$, these two connected components are merged, and we say that the one that appeared last dies and becomes merged with the one that appeared first: $C_2$ becomes merged with $C_1$. Finally, as $\alpha$ gets smaller and goes to $-\infty$, only one connected component, $C_1$, remains. The parameter $\alpha$ can be seen as time (going backwards), thus we speak in terms of birth times, death times, and lifetimes of connected components. In our example $C_1$ is born at $g(p_1)$ and never dies, and $C_2$ is born at $g(p_2)$ and dies at $g(v)$. Birth times are greater than death times since time is going backwards, and we define the lifetime of a connected component as its birth time minus its death time. The result of persistent homology is presented in what is called a **persistence diagram**. Each connected component is represented by a point (birth time, death time) in the plane. More precisely these points are in the half plane, below the diagonal $\{y = x\}$, since birth time is always greater than death time. In our example, we obtain a persistence diagram with two points: $(g(p_1), -\infty)$ and $(g(p_2), g(v))$, as depicted in Figure 2 B. Let us point out that a connected component arises at each peak of the function $g$, and that at each valley (local minimum) a connected component dies. With a slight abuse of language we call $p_i$ the peak of $g$ corresponding to the local maximum $p_i$, and we define $\tau_i$ the **prominence of a peak** $p_i$, as the lifetime of the connected component that is born at $g(p_i)$. It is important to understand that in practice the result of persistent homology is just a collection of unlabeled (birth time, death time) points in
the persistence diagram. The connected component that generated a given point in the
diagram is not indicated.

<figure><img src="image_3.png" /><figcaption>Figure 2: Persistent homology of a function.</figcaption></figure>

A fundamental property of persistent homology is what is called **stability**. If a connected component has a short lifetime, this is likely because it corresponds to a small peak of the function, due to noise. Thus the important quantity in a persistence diagram is the points that are far away from the diagonal. If there are a function *g* and an approximation *ĝ* of *g* that is close to *g*, then the stability theorem of [8] states, in simplified terms, that the persistence diagram of *ĝ* and the one of *g* are also close in the sense that they may differ close to the diagonal, but they should have approximately the same number of points far away from the diagonal. Thus the number of prominent peaks is a stable quantity under perturbations of the function. The stability property is illustrated in Figure 3.

In persistent homology, all the peaks are merged during the filtration and in the end
only one remains<fnref n="3" />. ToMATo performs a modified persistent homology, where not all the
peaks are merged. The merging condition is based on peak prominence, as explained in
Subsection 2.2.3.

**Remark 2.1.** Usually in persistent homology, there is a parameter $\alpha$ that increases from $-\infty$ to $+\infty$ and we are interested in the connected components of the sublevel sets of $g$, rather than superlevel sets with a parameter that decreases. In this setting persistence diagrams are then composed of points above the diagonal and not below. The authors of [7] chose a different perspective that is more natural and appropriate in the context of mode-seeking, but the ideas and results are the same.

<footnote n="3">³Actually, if the function is not continuous then its superlevel set of parameter $\alpha$ has several connected components even as $\alpha$ goes to $-\infty$. In this case, in the end several clusters remain.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3: A: A function *g* and a “noisy” approximation *ĝ*. B: Their respective persistence diagrams. The persistence diagram of *g* consists of the 2 black points, and the persistence diagram of *ĝ* consists of the red crosses.</figcaption></figure>

2.2.3. Merging ascending regions based on peak prominence

We have recalled the bases of gradient ascent mode-seeking and persistent homology.
On one hand with gradient ascent mode-seeking, we have a definition of peaks and their
ascending regions, that is, potential clusters. On the other hand with persistent homol-
ogy, we have a tool to obtain peak prominences leading to a hierarchy of the peaks of f.
We now see how to combine these two ideas to define a mode-seeking clustering method,
that is not sensitive to noise induced peaks of the density estimate.

Let us work on the example illustrated in Figure 4. Here the density estimate $\hat{f}$ of $f$
exhibits 5 peaks $p_1, \dots, p_5$. The ascending region of a peak $p_i$ is denoted by $\mathcal{A}(p_i)$. A
classical mode-seeking algorithm would partition the space according to these ascending
regions, resulting in 5 clusters $\mathcal{A}(p_1), \dots, \mathcal{A}(p_5)$. However, $\hat{f}$ being an estimate of our
data density $f$, visually we would be tempted to think that $p_2, p_3$ and $p_5$ are due to noise,
and that the only relevant peaks are $p_1$ and $p_4$, resulting in only 2 clusters: $\mathcal{A}(p_4) \cup \mathcal{A}(p_5)$
and $\mathcal{A}(p_1) \cup \mathcal{A}(p_2) \cup \mathcal{A}(p_3)$. ToMATo thus departs from the classical mode-seeking setting
by allowing merging between ascending regions, taking advantage of persistent homology
theory. In this way we are more likely to identify only 2 clusters, corresponding to unions
of ascending regions.

To keep only prominent peaks of $\hat{f}$, the idea of ToMATo is to fix a threshold parameter $\tau > 0$ and to perform persistent homology but with the additional condition that we merge only peaks of prominence less than $\tau$ with peaks of prominence greater than $\tau$. The choice of parameter $\tau$ is discussed in Section 2.2.4. It is equivalent to the choice of
<figure><img src="image_5.png" /><figcaption>Figure 4: A density estimate $\hat{f}$ with 5 local maxima $p_1, \dots, p_5$. The ascending regions are denoted by $\mathcal{A}(p_1), \dots, \mathcal{A}(p_5)$.</figcaption></figure>

the number of clusters.

In our example we have $\tau_1 \ge \tau_4 \ge \tau_5 \ge \tau_2 \ge \tau_3$ and if $\tau$ is chosen such that $\tau_5 \le \tau \le \tau_4$, then:

* The peak $p_1$ never dies.

* The peak $p_4$ never dies: it does not become merged with $p_1$ because it has prominence $\tau_4 \ge \tau$.

* The peaks $p_2$ and $p_3$ are merged with $p_1$.

* The peak $p_5$ is merged with $p_4$.

The final clusters are the unions of ascending regions that got merged together. Here two clusters are obtained: $\mathcal{A}(p_4) \cup \mathcal{A}(p_5)$ and $\mathcal{A}(p_1) \cup \mathcal{A}(p_2) \cup \mathcal{A}(p_3)$, as desired.

### 2.2.4. Choosing $\tau$ or equivalently the number of clusters

In practice, the number of clusters can directly be specified, instead of specifying $\tau$.

The main idea to choose the number of clusters is the following. When performing persistent homology, we obtain the peak prominences and they are represented in the form of a persistence diagram, which in our example would look like the one on Figure 4.
5. We recall that in practice the points of the diagram are not labeled. Looking at this persistence diagram, we can see two points far away from the diagonal and other points close to the diagonal, a priori due to noise. We thus decide that we want 2 clusters. This is equivalent to choosing $\tau_5 < \tau \le \tau_4$, as illustrated on Figure 5.

The ToMATo algorithm must thus be run twice. The first time it is run without specifying the number of clusters and it performs classical persistent homology. The result is a persistence diagram from which the number of points far away enough from the diagonal can be chosen as the desired number of clusters. Then, the algorithm is run a second time, specifying the desired number of clusters.

We highlight the fact that running the algorithm with $\tau = +\infty$ is just classical persistent homology: all the clusters that can be merged, are merged. On the opposite, ToMATo with $\tau = 0$ does not merge any clusters and the output clusters are the ascending regions of peaks.

<figure><img src="image_6.png" /><figcaption>Figure 5: The persistence diagram of the density estimate $\hat{f}$ of Figure 4. The parameter $\tau$ is chosen between $\tau_5$ and $\tau_4$ to obtain 2 final clusters corresponding to the peaks $p_1$ and $p_4$ of the density. Indeed, the persistence diagram points that lie below (respectively above) the line $y = x - \tau$ have lifetime greater (respectively smaller) than $\tau$.</figcaption></figure>

## 2.3. Explanation of the ToMATo algorithm

We now explain the ideas of the ToMATo algorithm in the discrete setting from a practical point of view.
<figure><img src="image_7.png" /><figcaption>Figure 6: An example of a directed tree $T$, with its root $r(T)$.</figcaption></figure>

The inputs to ToMATo are: data points $x_1, \dots, x_n$; a density estimate at these points; the distances between the data points. Let us denote by $\hat{f}_i$ the value of $\hat{f}$ at $x_i$. The idea of ToMATo is to mimick in the discrete setting the gradient ascent and ascending regions merging previously explained in the continuous setting. The computational cost is reduced by building graphs on top of the data. We therefore need a few definitions about graphs before explaining the algorithm.

### 2.3.1. Graph definitions

Let us recall that a graph is an abstract set of points, called vertices, together with a set of arcs going from one vertex to another vertex, called edges. The edge between a vertex $x_i$ and a vertex $x_j$ is denoted by $e_{ij}$. The vertex $x_i$ is then called the initial vertex of $e_{ij}$, and $x_j$ the final vertex of $e_{ij}$. A graph is said to be undirected if, for all edges $e_{ij}$ we have $e_{ij} = e_{ji}$, otherwise it is said to be directed. Finally, a graph is a binary graph if its edges $e_{ij}$ are equal to either 0 or 1. In the following, all graphs are binary.

**Definition 2.2.** *Directed tree ([15])*

Let $G$ be a directed graph. A set of edges $e_1, \dots, e_m$ is said to be a directed path from a vertex $x$ to a vertex $x'$, if $x$ is the initial vertex of $e_1$, if $x'$ is the final vertex of $e_m$, and if for $k \in [1, m-1]$ the final vertex of $e_k$ is the initial vertex of $e_{k+1}$.

A directed tree is a directed graph with a specified vertex $r$, called its root, such that:

1. Every vertex $x \neq r$ is the initial vertex of exactly one edge.

2. The vertex $r$ is the initial vertex of no edge.

3. There is no directed path from a vertex to itself (i.e. no cycles).

An example of a directed tree is depicted in Figure 6. Directed trees will be, in the discrete setting, the equivalent of the ascending regions of the continuous setting.

**Definition 2.3.** *Upper star*

Let $G$ be a graph, and assume that there is a function defined on the vertices of $G$. Let us denote by $x_1, \dots, x_m$ the vertices of $G$. For $k \in [1, m]$, we define the upper star of the vertex $x_k$ in $G$, as the set of edges connecting $x_k$ to other vertices with higher function values, along with these vertices. The set of vertices of the upper star of $x_k$ is denoted by $S_k$.
<figure><img src="image_8.png" /><figcaption>Figure 7: A well-chosen neighbourhood graph of a dataset of points in $\mathbb{R}^2$. The number $i$ corresponds to a point $x_i$, where the indexes have been reordered such that $\hat{f}_1 \ge \cdots \ge \hat{f}_n$. This is not completely realistic since a neighbourhood graph would naturally have more edges than the ones drawn.</figcaption></figure>

**Example 2.1.** *Constructing the $S_k$ sets for Fig. 7 we get: $S_1 = \emptyset$, $S_2 = \{1\}$, $S_3 = \{1\}$, $S_4 = \emptyset$, $S_5 = \{1\}$, $S_6 = \{1\}$, $S_7 = \{4, 6\}$, $S_8 = \{4\}$, $S_9 = \{4\}$, $S_{10} = \emptyset$, $S_{11} = \{10\}$, $S_{12} = \emptyset$, $S_{13} = \{10, 12\}$, $S_{14} = \{12\}$, $S_{15} = \{12\}$, $S_{16} = \{10\}$, $S_{17} = \{1\}$, $S_{18} = \{10\}$, $S_{19} = \{1, 10\}$, $S_{20} = \{1\}$.*

### 2.3.2. Principle of ToMATo

We are now ready to explain the pseudo-algorithm of ToMATo. In a preprocessing step, the ToMATo algorithm computes a neighbourhood graph $G$ of the data, which is typically a Rips graph or a $k$-nearest neighbour graph. The Rips graph with parameter $r > 0$ is the (undirected) graph where two points are connected if and only if the distance between them is less than $r$. A $k$-nearest neighbour graph is a graph where each point is connected to its $k$-nearest neighbours. The choice of $r$ or of $k$ is important. In Figure 7 a well-chosen neighbourhood graph of a dataset of points in $\mathbb{R}^2$ is represented as an illustration.

Let $G$ be a neighbourhood graph built on top of the data points. The principle of the algorithm is the following:

* Build the directed trees in $G$ having root at a local maximum of $\hat{f}$.

* Potentially merge some trees to obtain the desired number of clusters.

We emphasize the parallel with the continuous setting. Building directed trees having root at a local maximum of $\hat{f}$ is a mode-seeking step, it is the equivalent of searching peaks and their ascending regions. In the continuous case ascending regions are defined thanks to a gradient ascent, but there is no gradient needed here in the discrete setting for the tree construction. Each tree is then a potential cluster, and allowing merging between trees enables to recover some stability. The output of the algorithm is a union of directed trees, each union corresponding to a final cluster of the set of data points.
<figure><img src="image_9.png" /><figcaption>Figure 8: [Continuation of Figure 7 and Example 2.1.] The result of the ToMATo algorithm with $\tau$ chosen such that there are 2 final clusters. Each tree is circled in black, and the unions of trees circled in red correspond to the 2 final clusters.</figcaption></figure>

We explain more precisely the tree construction. A vertex $x_i$ is declared a peak if the set $S_i$ is empty. If a vertex $x_i$ is not a peak then it is attached to the tree containing argmax$_{j \in S_i}$ $\hat{f}_j$, forming a cluster. Iterating over all vertices, directed trees are thus created, the roots of which are peaks of $\hat{f}$ in the neighbourhood graph. These trees constitute the “ascending regions”. In Figure 8, each directed tree computed by the algorithm for the dataset of Figure 7 is circled. They are the potential clusters without any merging.

We now explain how to perform persistent homology on the neighbourhood graph. We directly explain the modified persistent homology with the condition that for a parameter $\tau$, we merge only peaks of prominence less than $\tau$ with peaks of prominence greater than $\tau$. Classical persistent homology is the same procedure, with $\tau = +\infty$. For each vertex $x_i$, we denote by $T_i$ the (necessarily unique) tree which includes $i$, and $r(T_i)$ the root of $T_i$. If two vertices $x_i$ and $x_k$ are in the same tree, we denote this tree by either $T_i$ or $T_k$; for instance we always have $T_i = T_{r(T_i)}$. The algorithm is the following: let $\tau > 0$, let $x_i$ be a vertex,

for $j$ in $S_i$: (we check if other trees $T_j$ can be merged with $T_i$)

if $\hat{f}_r(T_j) - \hat{f}_i < \tau$: merge $T_j$ with $T_i$

Let $T_{\max}(i)$ be, among the trees intersecting $S_i$ having a higher root than $T_i$, the tree with the highest root:

if $\hat{f}_r(T_i) - \hat{f}_i < \tau$: merge $T_i$ with $T_{\max}(i)$

(we check if $T_i$ itself can be merged with another tree).

In full words it means that if $x_i$ connects two peaks $T_i$ and $T_j$, respectively $T_i$ and $T_{\max}(i)$, then it is near a local minimum of $\hat{f}$. The prominence of the smallest peak is thus estimated respectively by $\hat{f}_r(T_j) - \hat{f}_i$ or $\hat{f}_r(T_i) - \hat{f}_i$, and if this estimated prominence is smaller than $\tau$ then the smallest peak is merged with the highest peak.

**Example 2.2.** [Continuation of example in Figures 7 and 8.]

Let us illustrate the algorithm for different values of $i$.
• Let $i = 13$:

$S_{13} = \{12, 10\}$, thus $x_{13}$ is attached to the tree $T_{10}$.

$j = 12$: if $\hat{f}_r(T_{12}) - \hat{f}_{13} = \hat{f}_{12} - \hat{f}_{13} < \tau$ then $T_{12}$ is merged with $T_{13} = T_{10}$.

$T_{\max}(13) = \emptyset$ (thus no further merging can be done).

• Let $i = 7$:

$S_7 = \{6, 4\}$, thus $x_7$ is attached to the tree $T_4$.

$j = 6$: if $\hat{f}_r(T_6) - \hat{f}_7 = \hat{f}_1 - \hat{f}_7 < \tau$ then $T_6 = T_1$ is merged with $T_7 = T_4$.

$T_{\max}(7) = T_1$.

If $\hat{f}_r(T_7) - \hat{f}_7 = \hat{f}_4 - \hat{f}_7 < \tau$ then $T_7 = T_4$ is merged with $T_1$.

If $\tau = +\infty$, the above procedure is the persistent homology of the filtration formed by the neighbourhood graphs of $\{x_i | \hat{f}_i \ge \alpha\}$, for $+\infty \ge \alpha \ge -\infty$. In Figure 8, we circle in red the 2 final clusters that we obtain when applying ToMATo with a value of $\tau$ chosen such that we obtain 2 clusters.

**Remark 2.2.** In the above procedure, trees can in fact be unions of trees.

**Remark 2.3.** If the neighbourhood graph has several connected components, then they are never merged during persistent homology. Consequently, when ToMATo is run with $\tau = +\infty$ and performs persistent homology of the graph, there are as many points at $y = -\infty$ in the diagram, as there are connected components in the graph. The number of connected components is thus also the minimal number of clusters that ToMATo can find, with any value of $\tau$.

We have thus defined in the discrete setting, relevant notions of peak, peak prominence, ascending region, and persistent homology. Therefore everything that we explained in the continuous setting adapts to the discrete setting.

We do not go into the details of the ToMATo implementation. The actual implementation is really clever and optimized, and a bit different from the pseudo-code presented in [7] and that is explained here. We only presented the general ideas, that are sufficient to understand the algorithm and its most important parameters, the latter being detailed in Subsection 2.4.

**Remark 2.4.** A few words about theoretical guarantees

The authors of [7] proved that the algorithm recovers the exact number of clusters, and that the output clusters coincide with the peaks of the density $f$. Let us present the results in a simplified way. For results 1 and 2, let us assume that the $n$ data points are i.i.d. and sampled under a density $f$ that is known.

Result 1 ([7], theorem 4.8):

If the persistence diagram of $f$ has a significant gap between points close to the diagonal and points far away from the diagonal, and if $n$ is large enough, then there exist a parameter $r$ for the Rips graph and a parameter $\tau$ such that with high probability, the number of clusters computed by the algorithm is equal to the number of peaks of $f$ of
prominence greater than $τ$.

**Result 2 ([7], theorem 4.9):**

If the persistence diagram of $f$ has a significant gap between points close to the diagonal and points far away from the diagonal, and if $n$ is large enough, then there exist a parameter $r$ for the Rips graph and a parameter $τ$ such that with high probability, for each peak of $f$ of prominence greater than $τ$, the algorithm outputs a cluster that coincides with the part of the ascending region that correspond to the top of the peak.

A third result states that with an approximation $\hat{f}$ of $f$ that is close to $f$, results 1 and 2 still hold. More precisely:

**Result 3 ([7], section 5):**

If the persistence diagram of $f$ has a significant gap between points close to the diagonal and points far away from the diagonal, if $n$ is large enough, and if $\hat{f}$ is close enough to $f$, then there exist a parameter $r$ for the Rips graph and a parameter $τ$ such that with high probability:

* The number of clusters computed by the algorithm with input $\hat{f}$, is equal to the number of peaks of $f$ of prominence greater than $τ$.

* For each peak of $f$ of prominence greater than $τ$, the algorithm with input $\hat{f}$, outputs a cluster that coincides with the part of the ascending region that correspond to the top of the peak.

## 2.4. ToMATo in practice with Gudhi

ToMATo is implemented in Gudhi<fnref n="4" />, and well documented. The ToMATo tutorial<fnref n="5" /> can also be consulted. The most important parameters are:

* 'density_type': the density type can be chosen among 'KDE', 'log KDE', 'DTM' or 'log DTM'. The DTM (Distance To a Measure) can be seen as an improvement over the k-NN estimator and is defined in [3]. The ToMATo parameters 'k_DTM', 'q' and 'dim' are associated with the DTM<fnref n="6" />, see [3] for details.

* 'graph_type': the graph type can be chosen among 'knn' or 'radius'. The choice 'radius' is the Rips graph; the parameter $r$ of which is to be specified by the user. The choice 'knn' is the k-nearest neighbour graph, with default $k = 10$.

* 'n clusters': the number of clusters.

<footnote n="4">⁴https://gudhi.inria.fr/python/latest/clustering.html.</footnote>

<footnote n="5">⁵https://github.com/xetaxe/ToMATo-Notebook/blob/master/Guide%20ToMATo.ipynb.</footnote>

<footnote n="6">⁶In high dimension one may need to change the 'dim' parameter. This is explained in the Gudhi documentation.</footnote>
The default parameter values are `density_type=log DTM` and `graph_type=knn` with $k = 10$. In practice the algorithm has to be run first without specifying `n clusters`, giving a persistence diagram. Then we see how many points are sufficiently far away from the diagonal, and the algorithm has to be rerun with this number as `n clusters`.

In our experience, the values that give the best results are usually the default values. By “best results”, we mean that the persistence diagram presents a clear gap between points close to the diagonal and points far away from the diagonal, allowing an un-equivocal choice of the number of clusters. ToMATo thus requires tuning only a few parameters.

# 3. Application of ToMATo in a spike sorting problem

The problem of spike sorting is the following. The raw data consists of the recording of mixed activity of multiple neurons, and the goal is to recover how many neurons contributed to the recording, and find the times at which each neuron fired a spike. The result of spike sorting is called a rasterplot. A rasterplot is a collection of discrete time series, presented in rows. Each row is a time series that corresponds to an identified neuron, and it shows the series of the times at which this neuron emitted a spike. Each discrete time series is represented as a series of bars. A rasterplot is showed in Figure 22 in Subsection 3.3 when ToMATo is performed on real data.

We apply ToMATo in the context of spike sorting on both simulated and real data, and we explain how it supersedes the usual method.

## 3.1. General outline and interest of the ToMATo method

We explain the general spike sorting outline.

Neuronal activity is recorded at different sites. In our simulated and real data there are 4 sites. We make the assumption that anytime a neuron fires a spike, the same underlying waveform with some additive auto-correlated Gaussian noise with variance 1 is recorded on each site. More precisely, there is one waveform per electrode and per neuron. The raw data, or the recorded neuronal activity, is then the mixture of the activity of several neurons.

First a preprocessing step to detect spikes and spike times, is performed. The goal is to identify times at which the recording presents a large local maximum at at least one of the 4 sites: the waveforms around such times are spike candidates. A collection of 4 waveforms corresponding to the detected spikes of a neuron on the 4 sites, is called an event. Each waveform on each site has 45 sampling points, thus each event is a 4×45 = 180 dimensional vector. At this point, this is a clustering problem. The goal is to
cluster the events, each obtained cluster being interpreted as the spikes from one neuron.

After clustering, the final step is to go back to the raw data and to the spike times,
and assign each spike time to a neuron in order to obtain a rasterplot.

In the usual spike sorting procedure, clustering cannot be performed directly. The
clustering phase, explained in full details in [21] and reexplained in Appendix A, can be
summarized with the following steps:

1. Obtaining clean events (eliminating superpositions):

A superposition occurs when two different neurons fire nearly simultaneously, lead-
ing to an event that is the superpersition of the spikes of the two neurons. In a
first step, the most obvious superpositions are eliminated, to keep only what are
called “clean” events<fnref n="7" />. This step is essential for the principal component analysis
that follows, otherwise superpositions skew the principal components.

2. Dimension reduction with **principal component analysis**:

A principal component analysis is performed on the clean events, to reduce the
dimension from 180 to a smaller integer *d*, typically *d* ≤ 10.

3. **Dynamic visualization of the data to find the number of clusters:**

Projections on the principal components of the data are visualized to help the user
choose the number of clusters. This step is not reproducible, the choice of the
number of clusters depends a lot on the user.

4. Automatic clustering:

An automatic clustering method in $\mathbb{R}^d$ is applied on the projections on the principal
components of the data, with K clusters where the integer K has been chosen in
the previous step. K-means is often used since it is a common method that gives
satisfying results, however in the case where clusters seem to be non convex in the
visualization step, other more appropriate clustering methods are used.

This procedure is laborious. Let us explain how ToMATo supersedes these four steps.

A first run of ToMATo, i.e., without specifying ‘n clusters’, provides an efficient way
to choose the right number of clusters, by looking at the persistence diagram. Thus, the
visualization of step 3 is no longer necessary. Since dimension reduction is performed
mostly to enable data visualization, and since the superposition elimination is performed
to not skew the dimension reduction, step 1 and 2 become superfluous as well. Finally,
once the number of clusters is chosen, a second run of ToMATo is an automatic clustering
method, the equivalent to step 4. ToMATo alone thus replaces the whole four-step
procedure. Let us highlight some facts:

<footnote n="7">⁷When the model will be subsequently used to classify data, superpositions have to be looked for and accounted for.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 9: The chosen standard shape.</figcaption></figure>

* For both simulated and real data, we apply ToMATo directly in high dimension (180) and it turns out that it works and runs within a reasonable time. In fact, with a dimension reduction it can be seen that our data intrinsically live at most in a 5 or 6 dimensional space, explaining why the high dimension does not raise an issue. For general data however, it may be necessary to perform a dimension reduction before applying ToMATo.

* The persistence diagrams obtained with ToMATo on our data often allow an un-equivocal choice of the number of clusters, but sometimes there can still be an ambiguity. This is discussed in Subsection 3.2.

* Our data present two main perturbation sources. The first one consists of the presence of superpositions, and the second one consists of small waveforms. In the simulations in Subsection 3.2 we find that ToMATo is extremely robust to superpositions. This fact shows that the superposition elimination of the usual spike sorting procedure becomes indeed unnecessary with ToMATo.

## 3.2. Simulations

We simulate directly the events on the 4 sites. To do so, we consider a standard spike shape that we multiply by a scale factor, called hereafter amplitude, to obtain the ideal waveform for each neuron at each site. The considered standard shape is displayed in Figure 9. Each neuron is thus characterized by 4 amplitudes. Then, some Gaussian noise is added (with constant variance equal to 1), to simulate the events. Let us recall that each event has $4 \times 45 = 180$ sampling points, it is thus a 180 dimensional vector. Finally, to create a superposition, 2 neurons are randomly selected, then the event of the second neuron is shifted from the event of the first one with a small random inter-val, such that the spikes of the two neurons are almost aligned. A superposition thus generates a different shape.
<figure><img src="image_11.png" /><figcaption>Figure 10: The 3 ideal shapes. Each shape corresponds to a neuron. For each shape, the 4 ideal waveforms on the 4 sites are concatenated horizontally.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Simulated events.</figcaption></figure>

The three main parameters for the simulations are:

1. the number of neurons,

2. the amplitudes,

3. the superposition frequency.

In the rest of Subsection 3.2, we apply ToMATo with all parameters set to default, except the 'dim' parameter that we set to 2 because of the high dimension.

### 3.2.1. A first example

To illustrate the principle of our simulations we present a first example. We consider 3 neurons, and amplitudes between 0 and 20. The 3 ideal shapes that we use, are depicted in Figure 10. We simulate events by adding Gaussian noise with constant variance equal to 1, and we add 40 percent of superpositions. Simulated events are displayed in Figure 11.

We apply ToMATo on this event collection. We first apply ToMATo without specifying `n clusters` to obtain the persistence diagram that is depicted in Figure 12 Left. Green points are points with ordinate $y = -\infty$. There are unequivocally 3 points far away from the diagonal. We then rerun ToMATo with `n clusters`= 3. To assess the quality of clustering, we compute the event medians of identified clusters and superpose them with the
<figure><img src="image_13.png" /><figcaption>Figure 12: Left: Persistence diagram. Right: Confusion matrix. The color scale goes from light blue (0) to dark blue (1).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13: Event medians of identified clusters (in red), superposed with the ideal shapes (in blue).</figcaption></figure>

ideal shapes. This is illustrated in Figure 13. We can see that the medians and the ideal shapes superpose perfectly. We also plot a confusion matrix (Figure 12 Right): the entry $(i, j)$ corresponds to the percentage of events of neuron $i$ assigned to cluster $j$. In the rest of this Subsection, only confusion matrices are displayed to summarize clustering results.

### 3.2.2. On the automatic identification of the number of neurons

Let us say a few words about the detection of the number of clusters on the persistence diagram. We implemented an automatic identification method adapted to our type of data. After performing simulations in the setting of Subsection 3.2.3, it turns out that we mostly obtain two types of diagrams. When the perturbation level (from either superpositions, or low amplitudes or both) is small, persistence diagrams look like the one in Figure 14 Left. They present a clear gap between an agglomeration of points extremely close to the diagonal, and points far away from the diagonal. When the perturbation level is important, persistence diagrams tend to look like the one in Figure 15 Left. This diagram shows one group of points on the top right of the diagram. They are not exactly points far away from the diagonal, but our experiments showed that in this case the right number of neurons is the number of points in this group.

To detect the number of neurons, we first draw a parallel line to the diagonal, with equation $y = x - \tau$ where $\tau$ is the mean of the lifetimes<fnref n="8" /> of the points. It gives a first

<footnote n="8">⁸We keep the persistent homology vocabulary. The “lifetime” of a point in a persistence diagram is an abuse of language to talk about the lifetime of the peak corresponding to this point. The lifetime of</footnote>
<figure><img src="image_15.png" /><figcaption>Figure 14: Left: A persistence diagram obtained from a simulation with 15 neurons, a superposition frequency of 0.05, amplitudes between 0 and 20. Right: The same diagram, with, in orange, the line with equation $y = x - \tau$ where $\tau$ is the mean of the lifetimes of the points. We detect 15 points below this line.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15: Left: A persistence diagram obtained from a simulation with 15 neurons, a superposition frequency of 0.9, amplitudes between 0 and 20. Right: The same diagram, with, in orange, the line with equation $y = x - \tau$ where $\tau$ is the mean of the lifetimes of the points. Points below this line are cut in two groups, blue and purple. The blue group is the rightmost group and we detect 15 points in this group.</figcaption></figure>

threshold to detect points far away from the diagonal. If the points below this line show only one group of points, we count the points in this group and take it to be the number of neurons: this is what happens in Figure 14 Right. If the points below this line show two groups of points, we identify these two groups by performing a K-means with $K = 2$, on the birth times<fnref n="9" /> of the points. Then we take the number of points in the rightmost group, as the number of clusters. This is what happens in Figure 15 Right, where we take only the blue points. It turns out that the diagrams obtained with our type of data always look like Figure 14 or Figure 15, so that this detection method works well and finds the right number of clusters when possible. All the simulations in Subsection 3.2.3 have been performed with this identification method.

<footnote>a point $(x, y)$ is thus defined as $(x - y)$.</footnote>

<footnote n="9">⁹This is still a persistent homology vocabulary. The birth time of a point $(x, y)$ is its abscissa $x$.</footnote>
### 3.2.3. Simulation results

#### Realistic simulation of locust and hippocampus data

Our simulations are motivated by real life data. We focus particularly on two cases of interest: the simulation of locust data and the simulation of hippocampus data. Locust data is characterized by the following set of parameters:

1. There are between 5 and 15 neurons.

2. The superposition frequency is smaller than 0.1, typically around 0.03.

3. Amplitudes are between 0 and 20.

This is a priori an easy setting: there are not too many neurons, not too many superpositions, and amplitudes are high. Our second study case is hippocampus data, which is more complicated. There can be up to 30 neurons, amplitudes are small and there are a lot of superpositions. More precisely for hippocampus data, we have that:

1. There are between 15 and 30 neurons.

2. The superposition frequency is between 0.1 and 0.5, typically around 0.3.

3. Amplitudes are between 0 and 10.

For both cases we simulate 50, 100 or 200 events per neuron. We apply the ToMATo method and we present a few results, that are typical of what we generally obtain.

First we show in Figure 16 the simulation results for locust data. The superposition frequency is set at 0.03, and the number of neurons vary from 5 to 15. The obtained confusion matrices have most of the weight on the diagonal: our clustering results are thus almost perfect.

In Figure 17 we show the results for the simulation of hippocampus data. The superposition frequency is set at 0.3, and the number of neurons vary from 15 to 30. For this type of data we do not recover all the neurons. This was predictable since it is always a difficult type of data to work with. In Figure 18 we plot one of the persistence diagrams, to illustrate the effect of perturbations on the identification of the number of neurons. We can see that there is no clear separation between points close to and far away from the diagonal. This is due to the small amplitudes. However, in this difficult context, identifying about 15 over 30 neurons is satisfying, especially if the neurons that are identified, are properly identified. As we can see in the confusion matrices, this is the case.

To conclude, ToMATo works perfectly for locust data, and it gives the most satisfying results as possible, for hippocampus data.
<figure><img src="image_17.png" /><figcaption>Figure 16: Confusion matrices for simulated locust data. 3 percent of superpositions, $0 \le \text{amplitudes} \le 20$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17: Confusion matrices for simulated hippocampus data. 30 percent of superpositions, $0 \le \text{amplitudes} \le 10$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 18: The persistence diagram for simulated hippocampus data of Figure 17, in the case of 29 neurons. 17 neurons are identified (points in the blue group).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 19: Simulation results for amplitudes between 0 and 20. The number of neurons vary: on the first row we take 10 neurons, on the second row we take 20 neurons and on the third row we take 30 neurons. On each row, the superposition frequency *f* increases from left to right: 0.01, 0.03, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. “Identified neurons” is abbreviated by “id. neurons”.</figcaption></figure>

## Robustness of ToMATo to perturbations

It is interesting to explore other parameter values even if they do not reflect a real life setting. In particular, it is interesting to study the robustness of ToMATo to superpositions and to low amplitudes since these are the two main perturbation sources in experimental (real) data. We present some results here. For a fixed number of neurons, we simulate events with different superposition frequencies from 0.01 to 1, once for amplitudes between 0 and 20 and once for amplitudes between 0 and 10.

In Figure 19, we show the results for amplitudes between 0 and 20. We take 10, 20 and 30 neurons, and for each fixed number of neurons we vary the superposition frequency. We can see that the confusion matrices show a lot of weight on the diagonal, even with a high superposition frequency. We can conclude that in this setting with high amplitudes, ToMATo is extremely robust to superpositions.

Simulations that give Figure 20 are performed in the same setting but for amplitudes between 0 and 10. Results are a little bit less satisfying but once again, there is no clustering method that could perfectly detect neurons with events of low amplitude. Moreover, we can observe that in this case too, that the results stay consistent when the superposition frequency increases.

From these simulations, we can conclude that:

* As any clustering method, ToMATo is not very robust to low amplitudes, but it still gives very satisfying results.
<figure><img src="image_21.png" /><figcaption>Figure 20: Simulation results for amplitudes between 0 and 10. The number of neurons vary: on the first row we take 10 neurons, on the second row we take 20 neurons and on the third row we take 30 neurons. On each row, the superposition frequency increases from left to right: 0.01, 0.03, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. "Identified neurons" is abbreviated by "id. neurons".</figcaption></figure>

* ToMATo is extremely robust to the presence of superpositions, which is very impressive.

## 3.3. Application on real data

We now apply ToMATo on a real locust dataset.

### 3.3.1. Description of the data

The data used here were recorded from the first olfactory relay, the antennal lobe, of a locust (*Schistocerca americana*). Recording setting and acquisition details are described in [22]. Measurements are performed using a tetrode, thus neuronal activity is recorded at 4 different sites.

In Figure 21 A, we show a 100 ms of data recorded at the 4 sites of a tetrode. The data were filtered between 300 Hz and 5 kHz before being digitized with a sampling frequency of 15 kHz. They were then normalized by dividing the signal amplitude at each site by a robust estimator, the median absolute deviation, of the noise deviation σ<sub>noise</sub>. In Figure 21 B, after detection of the spike candidates as sufficiently large local extrema in absolute value, cuts are made on each of the four sites, each cut having 45 sampling points. This group of four cuts determines an event, of dimension 4 × 45 = 180. The procedure is detailed in [21] and summarized in Appendix A.
<figure><img src="image_22.png" /><figcaption>Figure 21: Recordings data and event detection.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 22: Left: Persistence diagram. Right: The 2 first seconds of the final obtained rasterplot.</figcaption></figure>

### 3.3.2. Result of ToMATo

Let us apply ToMATo to this event collection. The first run of ToMATo gives the persistence diagram displayed in Figure 22 Left. There are clearly 6 points far away from the diagonal thus we run ToMATo again with `n clusters= 6`. In Figure 22 Right, we display the obtained rasterplot. Each horizontal row corresponds to one of the 6 identified neurons, and shows a series of the times at which this neuron emitted a spike.

**Remark 3.1.** With the usual method (presented in https://c_pouzat.gitlab.io/spike-sorting-the-diy-way/ and summarized in Appendix A), with the same dataset it is not clear whether one should choose 5 or 6 clusters. With the ToMATo method there is no ambiguity.

The code and the data used in this article are available at https://gitlab.math.unistra.fr/martineau/tomato_for_spike_sorting.

## 4. Conclusion

The ToMATo clustering algorithm (Sec. 2) is part of a recent branch of applied mathematics, topological data analysis (TDA) [6, 9]; a branch that grew out of a rather sophisticated and abstract domain of mathematics, algebraic topology. Despite of this abstract origin, the key ideas on which ToMATo is built can (and should, we think) be grasped by any serious practitioner of neurophysiological data analysis. The introduction
and use of the concept of **peak prominence** (Sec. 2.2.2) gives rise to a mode-seeking algo-
rithm that is **robust with respect to noise**. The use of a graph built from the estimated
density values at the data points, **and only at the data points**, leads to an algorithm that
is computationally very efficient (Sec. 2.3). Serious mathematical studies have moreover
provided theoretical guarantees (Sec. 2.3.2) on the capabilities of ToMATo to recover
the right clusters. These considerations convinced us that an exploration of ToMATo in
a neurophysiological data analysis context was a worthy endeavor.

As a test case, we chose a subject we know reasonably well: spike-sorting (Sec. 3).
More precisely, the part of spike sorting where clustering algorithms play a key role is
the determination of the number of “good” neurons—a good neuron is a neuron whose
spikes can be reliably identified—, together with their waveform/template (Sec. 3.1 and
Appendix A). For our test, both simulated (Sec. 3.2.3) and real data (Sec. 3.3) were
used. The very important preliminary question, the determination of the number of clus-
ters/neurons, was discussed twice; first in the general presentation of ToMATo (Sec. 2.2.3
and 2.2.4), next in the specific context of spike sorting (Sec. 3.2.2). We considered two
broad simulation scenarios : “locust”, the easy case, with a large signal to noise ratio
(SNR), few neurons, leading to few superpositions; and “hippocampus”, with a low SNR
and many neurons, leading to frequent superpositions. We showed that ToMATo does
not require a preliminary dimension reduction and is able to reliably identify the “good”
neurons **despite of a potentially large number of superposed events**. These findings
were confirmed using our “usual” real dataset.

We hope that this report has convinced our readers that the ToMATo clustering
algorithm provides an attractive alternative to more traditional algorithms (kmeans,
Gaussian mixture models, etc.). ToMATo is moreover just a tiny part of a mature,
comprehensive and well documented topological data analysis C++ library: Gudhi. The
library is interfaced with Python (and R) and can be readily tried by anyone willing to
spend a few time reading the relevant documentation. Don't hesitate!

Acknowledgments

The authors would like to thank Marc Glisse from INRIA for carefully answering our
questions about the ToMATo implementation. This work was supported by the Agence
Nationale de la Recherche (ANR): project ANR-22-CE45-0027 SIMBADNESTICOST.

**A. A spike sorting outline**

Spike sorting aims at extracting from “raw data” (nowadays continuous recordings from
several channels/electrodes, like the 4 channels of Fig. 23A) sequences of spike times
emitted by “identified” neurons [10, 20, 16]. The raw data are typically a mixture of
waveforms/motifs from different neurons with an added independent recording noise.
The data generation model that is most of the time tacitly assumed, but sometimes spelled out [24], is that the neurons generate dependent marked point processes [4], where the mark of each neuron on each channel/electrode is a waveform/motif; when two or more neurons generate events whose time separation is smaller than the duration of their waveforms, these waveforms are summed on their overlapping regions. A white or “colored” noise is then added to the realizations of these marked point processes. We wrote that the raw data are typically **continuous** recordings, but this descriptive statement should not mislead the mathematically oriented reader: the data are in fact **sampled**—that is, measurements are performed at a fixed frequency (15 kHz for the data of Fig. 23), implying that the data are intrinsically discrete, they just “look” continuous—and they are filtered before being sampled, a high-pass filter (300 Hz for the data of Fig. 23) removes low frequency oscillations and a low-pass filter (5 kHz for the data of Fig. 23) ensuring that the requirement of Nyquist-Shannon sampling theorem [25] are met (the sampling frequency must be at least twice as large as the low-pass filter cut-off frequency). Most of the (really) many spike sorting algorithms that have been proposed can be decomposed in 3 main stages illustrated by the 3 rows on Fig. 23:

1. Spike detection (Fig. 23A, in red), and event sample construction (Fig. 23B and C).

2. Estimation of the number of neurons and of their waveforms (Fig. 23D, E and F).

3. Event assignment and superposition resolution (Fig. 23G₁, G₂, G₃).

The present manuscript focuses on the second point. For completeness we state below the main operations we usually apply to the data in order to provide a clear “reference” procedure with which our new approach can be compared (referring to Fig. 23)<fnref n="10" />:

A Events detection (putative action potentials) based on extrema exceeding a threshold.

B Cuts / windows, one on each site, of “well-chosen” length (here 45 sampling points) around the detected extremes, this collection of four cuts (as we have here four sites) constitutes **an event** (our event space is here $\mathbb{R}^{180}$ as we have 4×45 amplitudes per event).

C The first 200 detected events aligned on their valley (events made of superpositions are displayed in red).

D **Dimension reduction**, here the projection of the sample on a plane of the subspace defined by the first three principal components.

E **Clustering** with the k-means method and 10 centers.

<footnote n="10">¹⁰The data and Python codes leading to this figure are publicly available: https://c_pouzat.gitlab.io/spike-sorting-the-diy-way/.</footnote>
<figure><img src="image_24.png" /><figcaption>Figure 23: Spike sorting as a succession of “simple” tasks (see text and https://c_pouzat.gitlab.io/spike-sorting-the-diy-way/ for details). These images are from https://c_pouzat.gitlab.io/spike-sorting-the-diy-way/.</figcaption></figure>
F The motifs (centers of the clusters defined in the previous step) corresponding to the 5 “largest” neurons (the 5 different colors) on each of the 4 sites.

G1 Return to the raw data (black trace) and attribution of a motif to each local extremum generating a prediction (red trace).

G2 In black the difference between the black and red traces of G1, a detection of local extrema is performed again and the closest motif is assigned to each extremum, giving rise to a new prediction (red trace).

G3 The difference between the black and red traces of G2, we continue this “peeling” procedure until there is nothing left identifiable to any of the motifs of the collection.

This somewhat long list should be taken as **an example** of what is done when doing spike sorting, since there are usually several options/approaches at each step. We use for instance PCA for dimension reduction [11], but ICA can be (and is) used. The “peeling procedure” of the third row [19, 21], a sophisticated form of **template matching**, is just one way of solving the superposition problem [23]. The fact that the data are sampled gives rise to a **sampling jitter** [17] that must be dealt with in order to get good superposition resolution. This sampling jitter together with the presence of superpositions makes spike sorting a **non trivial clustering problem**. As explained in the Introduction, the superpositions must be eliminated in an ad-hoc way in order to have useful principal components when dimension reduction is used. We are not dealing here with the added difficulties met in some datasets where the waveform of a given neuron depends on the history of its discharge [5]; using closely spaced recording electrodes and working with the amplitude ratios on the different site is usually the best way to proceed [18, 12].
References

[1] Hadeel K. Aljobouri, Hussain A. Jaber, Orhan M. Koçak, Oktay Algin, and Ilyas Çankaya. Clustering fmri data with a robust unsupervised learning algorithm for neuroscience data mining. *Journal of Neuroscience Methods*, 299:45-54, 2018.

[2] Michael J. Berry and Gašper Tkačik. Clustering of neural activity: A design principle for population codes. *Frontiers in Computational Neuroscience*, 14, 2020.

[3] Gérard Biau, Frédéric Chazal, David Cohen-Steiner, Luc Devroye, and Carlos Rodríguez. A weighted k-nearest neighbor density estimate for geometric inference. *ejs*, 5(none):204-237, January 2011.

[4] Pierre Brémaud. *Point Process Calculus in Time and Space: An Introduction with Applications*. Probability Theory and Stochastic Modelling. Springer International Publishing, 2020.

[5] William H Calvin. Some simple spike separation techniques for simultaneously recorded neurons. *Electroencephalography and Clinical Neurophysiology*, 34(1):94-96, 1973.

[6] Gunnar Carlsson. Topology and data. *Bulletin of the American Mathematical Society*, January 2009.

[7] Frédéric Chazal, Leonidas J. Guibas, Steve Y. Oudot, and Primoz Skraba. Persistence-Based Clustering in Riemannian Manifolds. *J. ACM*, 60(6):1-38, November 2013.

[8] David Cohen-Steiner, Herbert Edelsbrunner, and John Harer. Stability of Persistence Diagrams. *Discrete Comput. Geom.*, 37(1):103-120, January 2007.

[9] Edelsbrunner, Letscher, and Zomorodian. Topological Persistence and Simplification. *Discrete Comput. Geom.*, 28(4):511-533, November 2002.

[10] Gaute T Einevoll, Felix Franke, Espen Hagen, Christophe Pouzat, and Kenneth D Harris. Towards reliable spike-train recordings from thousands of neurons with multiple electrodes. *Current Opinion in Neurobiology*, 22(1):11-17, 2012. Neurotechnology.

[11] E.M. GLASER and W.B. MARKS. On-line separation of interleaved neuronal pulse sequences. In KURT ENSLEIN, editor, *Data Acquisition and Processing in Biology and Medicine*, pages 137-156. Pergamon, 1968.

[12] Charles M. Gray, Pedro E. Maldonado, Mathew Wilson, and Bruce McNaughton. Tetrodes markedly improve the reliability and yield of multiple single-unit isolation from multi-unit recordings in cat striate cortex. *Journal of Neuroscience Methods*, 63(1):43-54, 1995.
[13] Trevor Hastie, Robert Tibshirani, and Jerome Friedman. *The Elements of Statistical Learning*. Springer, 2009.

[14] Bálint Király and Balázs Hangya. Navigating the statistical minefield of model selection and clustering in neuroscience. *eNeuro*, 9(4), 2022.

[15] Koontz, Narendra, and Fukunaga. A Graph-Theoretic Approach to Nonparametric Cluster Analysis. *IEEE Trans. Comput.*, C-25(9):936-944, September 1976.

[16] Arthur Leblois and Christophe Pouzat. *Neurobiology of Motor Control*, chapter Multi-Unit Recording. Wiley, 2017.

[17] Kevin C. McGill and Leslie J. Dorfman. High-resolution alignment of sampled waveforms. *IEEE Trans Biomed Eng*, 31(6):462-468, 1984.

[18] Bruce L. McNaughton, John O'Keefe, and Carol A. Barnes. The stereotrode: A new technique for simultaneous isolation of several single units in the central nervous system from multiple unit records. *Journal of Neuroscience Methods*, 8(4):391-397, 1983.

[19] Jonathan W. Pillow, Jonathon Shlens, E. J. Chichilnisky, and Eero P. Simoncelli. A model-based spike sorting algorithm for removing correlation artifacts in multi-neuron recordings. *PLOS ONE*, 8(5):1-14, 05 2013.

[20] Christophe Pouzat. *Brain-Computer Interfaces 1: Methods and Perspectives*, chapter Analysis of Extracellular Recordings. iSTE/Wiley, 2016.

[21] Christophe Pouzat and Georgios Is Detorakis. Spysort: Neuronal spike sorting with python. Proceedings of the 7th European Conference on Python in Science (EuroSciPy 2014), pages 27-34, 2014.

[22] Christophe Pouzat, Ofer Mazor, and Gilles Laurent. Using noise signature to optimize spike-sorting and to assess neuronal classification quality. *J. Neurosci. Methods*, 122(1):43-57, December 2002.

[23] V.J Prochazka, B Conrad, and F Sindermann. A neuroelectric signal recognition system. *Electroencephalography and Clinical Neurophysiology*, 32(1):95-97, 1972.

[24] William M. Roberts. Optimal recognition of neuronal waveforms. *Biological Cybernetics*, 35:73-80, 1979.

[25] M. Unser. Sampling-50 years after Shannon. *Proceedings of the IEEE*, 88(4):569-587, 2000.