# Persistent homology of spike trains with the Curto-Itskov filtration: stability and applications

Louise Martineau, Christophe Pouzat, Ségolen Geffray

► To cite this version:

Louise Martineau, Christophe Pouzat, Ségolen Geffray. Persistent homology of spike trains with the Curto-Itskov filtration: stability and applications. 2026. hal-05519731

HAL Id: hal-05519731

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

Preprint submitted on 19 Feb 2026

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 CC BY 4.0 - Attribution - International License
# Persistent homology of spike trains with the Curto-Itskov filtration: stability and applications

Louise Martineau, Christophe Pouzat, Ségolen Geffray

February 19, 2026

## Abstract

Following the work of [18], we define a filtration for spike trains based on the frequencies of cofiring neurons. This filtration, called the Curto-Itskov filtration, allows to define the persistence diagram of a spike train. We then introduce a distance on the space of spike trains and prove the stability of persistence diagrams, with respect to this distance. Finally, we illustrate the behaviour of the Curto-Itskov filtration with simulations, and apply it on a real dataset with a clustering goal.

## 1 Introduction

Actual neural networks are nowadays experimentally mostly studied with *extracellular recordings*. These recordings lead, after a rather heavy processing called spike-sorting, to sequences of action potentials, or *spike trains*, of individual neurons [31, 25]. The harsh reality of these recordings is that a *very tiny fraction* of the neurons making a local network gets recorded. This has obvious consequences for the spike trains interpretation: some unobserved neurons can form synapses on several of the observed ones (that is typically the case when the unobserved neuron is an inhibitory one [26], but this can happen with a unobserved excitatory neuron as well); such common input can lead to coordinated pauses (if the unobserved is inhibitory) or spikes (if the unobserved neuron is excitatory) in groups of observed neurons. A local network is moreover never isolated from the rest of the brain, at least in awake, behaving animals; this means that neurons outside of the recorded area can give inputs to the observed neurons. This is particularly the case when the local network under investigation is an early sensory relay receiving inputs from receptors. Then again, common unobserved inputs are expected, generating also potentially coordinated spikes among the observed neurons. Such coordinated spiking (or pausing) phenomena of more than two neurons are badly handled by
classical spike train analysis methods that focus on pairwise interactions, either using "cross-correlograms" [30, 5], or stochastic intensity models [4, 13]. This problem has been acknowledged for a long time by users of cross-correlograms and more recently by users of stochastic intensity models, and *partial* solutions have been proposed [28, 23]. On the other hand, the concept of *cell assemblies*, that are (or would be) materialized by groups of neurons cofiring within a short time window, is now pervasive in neuroscience [7]. There is therefore a need for data analysis methods able to identify such groups/assemblies from experimental data [8]. To that end, we propose here to adapt, to spike train data, tools from *topological data analysis* (TDA), that were introduced for the analysis of fMRI data [18].

Persistent homology is a tool from the field of topological data analysis, that has become more and more popular in neuroscience throughout the last decade (reviewed by [17]). The general persistent homology pipeline is described in Figure 1. The first step is to turn the data into a filtration of simplicial complexes, and the second step is to take the homology of this filtration. The final result is called a persistence diagram. This diagram gives interesting information about the data, that complements information arising from other data analysis techniques.

The usual framework of persistent homology in topological data analysis is the study of point cloud data. In this framework there is a stability property stating that if two point clouds are close then their persistence diagrams are close too, with appropriate distances for both point clouds and persistence diagrams. This stability theorem is fundamental for a statistical approach of topological data analysis, since it provides a noise robustness guarantee.

The study of point cloud data with persistent homology has been covered in the literature([21]). Thus even with other types of raw data, a topological analysis can be carried out by building point clouds on top of data, if it is relevant in the context. Sometimes however, it can be more appropriate to directly adapt the persistent homology pipeline to the encountered data. In neuroscience, scientists often deal with collections of binary time series. Such time series appear in neurophysiology in the study of neuronal activity: in this context they are called spike trains. Binary time series can also appear when thresholding fMRI data, as in the article [18] that introduced the Curto-Itskov filtration. In the persistent homology pipeline, step 2 is never affected: taking the homology of a filtration of simplicial complexes is a procedure that does not need any tuning<fnref n="1" />. Step 1 however, i.e. the filtration choice, can — and must — be adjusted to the data type and to the information we want to extract. The authors of [18] introduced a filtration for binary time series, based on appearance frequencies of variables that are all equal to 1 at the same time. In this article we exploit this idea in the context of spike trains, and we prove a stability property for this filtration, i.e. if two spikes trains are close for a

<footnote n="1">¹This is not entirely true: we must choose the coefficient field. In this article we consider homology with $\mathbb{F}_2$ coefficients.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Persistent homology pipeline.</figcaption></figure>

distance which we define, then their persistence diagrams are close too.

The paper is organised as follows. In Section 2 we recall the outlines of persistent homology. In Section 3, we define the Curto-Itskov filtration, based on frequencies in a spike train. In Section 4, we define a distance on the space of spike trains and prove a stability theorem. We present simulations in Section 5, and finally in Section 6 we apply the Curto-Itskov filtration on a real dataset.

# 2 Persistent homology reminders

In this section we recall basic facts about persistent homology. The usual framework of persistent homology of point cloud data is also introduced, because it is a benchmark and it is the natural context to apprehend persistent homology for the first time.

## 2.1 Ideas of persistent homology

### 2.1.1 Simplicial complex and simplicial homology

Topology describes some properties about shapes. More precisely, it describes properties that are preserved under continuous deformations. Homology consists of discrete topological descriptors of a shape. It counts the number of connected components, the number of 1-dimensional holes (the interior of a triangle), of 2-dimensional holes (the interior of a sphere), and so on. A “shape”, in the most general setting, is a topological space. A simpler notion of shape is a simplicial complex. Simply put, a simplicial complex is a collection of points, segments, filled triangles, filled tetrahedrons and so on. An example is given in Figure 2. Homology of simplicial complexes is easy to compute, which makes these objects at the core of topological data analysis. Let us give the definitions.

#### Simplicial complex

**Definition 2.1. Simplex**

Let $k, d$ be nonnegative integers. A $k$-simplex $\sigma$ is the convex hull of $(k + 1)$ affinely independent points in $\mathbb{R}^d$.

Let $\sigma$ be a $k$-simplex. The points that span $\sigma$ are called the vertices of $\sigma$, the number $k$ is called the dimension of $\sigma$. If $\sigma$ is spanned by points $v_1, \dots, v_{k+1}$, we denote
$\sigma = [v_1, \dots, v_{k+1}]$. Any simplex spanned by a subset of the vertices of $\sigma$, is called a face of $\sigma$.

**Definition 2.2.** *Simplicial complex*

A simplicial complex $K$ in $\mathbb{R}^d$ is a collection of simplices in $\mathbb{R}^d$ such that,

* Every face of a simplex of $K$ is in $K$.

* The intersection of any two simplices of $K$ is a face of each of them.

<figure><img src="image_3.png" /><figcaption>Figure 2: A simplicial complex.</figcaption></figure>

Given a simplicial complex $K$ in $\mathbb{R}^d$, the geometric realization $|K|$ of $K$ is defined as the subset of $\mathbb{R}^d$ that is the union of the simplices of $K$. The topology of $K$ is the topology induced on $|K|$ by the standard topology in $\mathbb{R}^d$.

A simplicial complex is characterized by the list of its simplices. It can thus be treated in a strictly combinatorial way without thinking of a geometric realisation in a space $\mathbb{R}^d$ for some dimension $d$. This leads to the definition of an abstract simplicial complex.

**Definition 2.3.** *Abstract simplicial complex*

An abstract simplicial complex is a collection $\mathcal{K}$ of finite non-empty sets, such that if $A$ is an element of $\mathcal{K}$, so is every nonempty subset of $A$.

If $\mathcal{K}$ is an abstract simplicial complex, its vertex set is defined as the union of one-point elements of $\mathcal{K}$. Finally, let us introduce the notion of isomorphic simplicial complexes.

**Definition 2.4.** *Isomorphic simplicial complexes*

Two abstract simplicial complexes $\mathcal{K}$ and $\mathcal{K}'$ are isomorphic if there is a bijective map $f$ between the vertex set of $\mathcal{K}$ and the vertex set of $\mathcal{K}'$, such that

$$ \{v_1, \dots, v_{k+1}\} \in \mathcal{K} \text{ if and only if } \{f(v_1), \dots, f(v_{k+1})\} \in \mathcal{K}'. $$

**Homology groups**

Homology groups is the rigourous notion to count $k$-dimensional holes of a shape, here a simplicial complex. We consider homology groups with coefficients in the field $\mathbb{F}_2 = \{0, 1\}$. The interested reader can consult the first chapter of [29] for a more detailed presentation on simplicial homology.

In the following, let $\mathcal{K}$ be a simplicial complex and let $k$ be a nonnegative integer.
**Definition 2.5.** k-chain

A k-chain on K is a formal sum with coefficients in $\mathbb{F}_2$ of k-simplices of K, i.e. a k-chain is a sum $\sum_{i=1}^p \lambda_i \sigma_i$ with $p \in \mathbb{N}^*$, $\lambda_i \in \mathbb{F}_2$, and $\sigma_i$ k-simplices of K.

The sum of two k-chains is defined by summing up the coefficients, and the multiplication of a k-chain by a scalar $\mu \in \mathbb{F}_2$ by multiplying the coefficients by $\mu$. This leads to the following definition.

**Definition 2.6.** Space of k-chains

The space of k-chains on K, denoted by $C_k(K)$, is defined as the set of k-chains of K with the two operations defined above.

The space $C_k(K)$ is a vector space over $\mathbb{F}_2$, and it is spanned by the k-simplices of K.

**Definition 2.7.** Boundary operator

The boundary operator $\partial_k : C_k(K) \to C_{k-1}(K)$ is defined as follows. If $\sigma = [v_1, \dots, v_{k+1}]$ is a k-simplex, we set

$$ \partial_k(\sigma) = \sum_{i=1}^{k+1} [v_1, \dots, \check{v}_i, \dots, v_{k+1}], $$

where $[v_1, \dots, \check{v}_i, \dots, v_{k+1}]$ is the simplex spanned by $\{v_1, \dots, v_{k+1}\} \setminus v_i$.

**Lemma 2.1.** [29]

$$ \partial_{k-1} \circ \partial_k = 0. $$

**Definition 2.8.** k-cycles and k-boundaries

The kernel of $\partial_k$ is called the space of k-cycles and denoted $Z_k(K)$. The image of $\partial_{k+1}$ is called the space of k-boundaries and is denoted by $B_k(K)$.

**Definition 2.9.** Homology groups and Betti numbers

The quotient vector space $H_k(K) = Z_k(K)/B_k(K)$ is called the k-th homology group of K. The k-th Betti number of K is defined as the dimension of $H_k(K)$, and is denoted by $\beta_k(K)$.

The k-th Betti number represents the number of linearly independent k-dimensional holes of K (a 0-dimensional hole is a connected component). Indeed, a cycle that is not a boundary, is a good notion of hole.

**Example 2.1.** It can be checked that the simplicial complex K represented in Figure 2 has the following Betti numbers: $\beta_0(K) = 3$, $\beta_1(K) = 3$, $\beta_2(K) = 1$, and $\beta_k(K) = 0$ for every $k > 2$. It coincides with the fact that K has 3 connected components, 3 holes of dimension 1 (the interior of the 3 empty triangles), and 1 hole of dimension 2 (the interior of the empty tetrahedron).

Homology groups of more general topological spaces can also be defined. They have the same interpretation than homology groups of simplicial complexes, that is the Betti numbers represent the number of k-dimensional holes. In the rest of the article, “k-dimensional holes” are sometimes referred to as “homological features of dimension k”.
**Definition 2.10.** *Simplicial map*

Given $K$ and $K'$ two simplicial complexes, we say that $f$ is a simplicial map from $K$ to $K'$, if $f: |K| \rightarrow |K'|$ is a continuous map that maps each simplex of $K$ linearly onto a simplex of $K'$.

Finally let us state that a simplicial map between simplicial complexes induces a linear map between homology groups.

**Proposition 1.** [29]

Let $K$ and $K'$ be simplicial complexes and $f$ be a simplicial map from $K$ to $K'$. Then $f$ induces a linear map $f^*: H_k(K) \to H_k(K')$ for each $k$. Moreover, the identity simplicial map on $K$ induces the identity linear map on $H_k(K)$; and given two simplicial maps $f$ and $g$, we have $(f \circ g)^* = f^* \circ g^*$.

### 2.1.2 From data to simplicial complex

The overall goal of topological data analysis is to study topological aspects of “the shape of data”. This sentence can have or not have a natural meaning, depending on the type of data. The most basic example is when data is a point cloud, i.e. a collection of points in a metric space. Then, speaking about its “shape” has a natural meaning: for instance if points are sampled on a circle, then the collection of points visually form more or less a circle. In this example the goal of a topological analysis is to recover the topological properties of a circle, that is, one connected component and one 1-dimensional hole.

Since the homology of a simplicial complex can easily be computed, the idea of topological data analysis is to define a simplicial complex on top of the data, and then to consider its homology.

Going back to the point cloud example, the most popular way to build a simplicial complex on top of a point cloud is the following. A parameter $r > 0$ is fixed, and we consider the union of balls of radius $r$ centered at each point. A simplicial complex called the Čech complex is then defined hereafter. Examples of Čech complexes are displayed in Figure 3.

**Definition 2.11.** *Čech complex on top of a point cloud ([3])*

Let $X = \{x_1, \dots, x_n\}$ be a set of points in a metric space, and let $r$ be a nonnegative real number. The Čech complex $\text{Cech}_r(X)$ is the abstract simplicial complex whose simplices are $\{x_1, \dots, x_{k+1}\}$ such that the $(k+1)$ closed balls $\bar{B}(x_i, r)$ have a non-empty intersection.

In Figure 4, we display the union of balls centered at the points of a point cloud, for different radius values. The points are sampled on a circle: the point cloud is represented on the left (radius $r_1 = 0$). In this example we can see that for a well-chosen parameter $r$, here $r_3$ or $r_4$, the Betti numbers of the union of balls are the same as the ones of a circle. Moreover, it turns out that the Čech complex $\text{Cech}_r(X)$ of a point cloud $X = \{x_1, \dots, x_n\}$, has the same Betti numbers as the union of balls $\cup_{1 \le i \le n} \bar{B}(x_i, r)$.
<figure><img src="image_4.png" /><figcaption>Figure 3: Top: A 3-points point cloud with the union of balls centered at these points for some fixed radius (left), and the corresponding Čech complex (right). Bottom: The same point cloud with the union of balls centered at the points for a larger radius (left), and the corresponding Čech complex (right).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: A point cloud and the union of balls centered at the points, for different radius values. The radius $r_1$ is equal to 0.</figcaption></figure>

This is a consequence of the nerve theorem ([22], Corollary 4.G.3). Therefore, for a well-chosen radius, computing the homology of the Čech complex reveals the desired topological properties of the point cloud.

For types of data different than points clouds, e.g. spike trains, speaking of “topological properties of the data” may have no natural meaning. But if a simplicial complex can be defined on top of the data, its homology can be considered and we can still talk about a topological analysis of the data. In Section 3 we use the Curto-Itskov construction to build a simplicial complex on top of spike train data. This construction involves a cofiring frequency parameter $f$, in the same way that the definition of the Čech complex depends on a radius parameter $r$.

### 2.1.3 Filtration

In the point cloud example we pointed out that for a well-chosen parameter $r$ for the Čech complex, the desired topological properties of the point cloud can be found. However there is no practical way to pick up the right parameter. And like the Čech complex,
most simplicial complexes built on top of data require a paramater choice.

To circumvent this issue, the principle of persistent homology is to consider all pa-
rameter values in a multiscale manner and to keep track of the evolution of homological
features of different dimensions, as the parameter value changes. This leads to the notion
of filtration.

**Definition 2.12. Filtration**

A filtration of simplicial complexes is a sequence of nested simplicial complexes.

For instance,in Figure 4 the sequence of Čech complexes for the increasing radius val-
ues defines a filtration. Such a filtration is called a Čech filtration.

The filtration choice is the key point in persistent homology. The challenge is to define
a filtration that reflects the desired information about the data.

2.1.4 Persistence diagram

Throughout a filtration, homological features of different dimensions appear, and disappear as they become merged with one another. We say that the feature that appeared last disappears when it becomes merged with the feature that appeared first. The terms of “birth time” and “death time” are employed in this setting, and the result of persistent homology is the collection of birth and death times of homological features. This result is presented in what is called a persistence diagram: a collection of points (birth time, death time) in the plane. It is important to highlight that there is one persistence diagram per dimension: there is one persistence diagram for homological features of dimension 0, one for homological features of dimension 1, etc... In Figure 5, persistence diagrams for dimension 0 and 1 for the Čech filtration of the point cloud of Figure 4 are shown as an example.

On the one hand, long living homological features give points far away from the
diagonal in the persistence diagram, and these points are interpreted as topological in-
formation. In the persistence diagrams of Figure 5 for instance, there is one point far
away from the diagonal for dimension 0 and one point far away from the diagonal for
dimension 1. They are interpreted as one connected component and one 1-dimensional
hole, so that we recover the Betti numbers of a circle. On the other hand, short living
homological features give points close to the diagonal in the persistence diagram, and
these points are interpreted as noise. They do not bring topological information but
they can have discriminatory power in a clustering context.

The authors of [9] proved that persistence diagrams constructed from Čech filtrations
are robust to noise. Indeed if two point clouds are close, then their persistence diagrams
are close too, for an appropriate distance. This property is called stability. In Section 4,
we prove a similar property for persistence diagrams constructed from the Curto-Itskov
filtration of spike trains.
<figure><img src="image_6.png" /><figcaption>Figure 5: The persistence diagrams of the Čech filtration of the point cloud of Figure 4. The persistence diagram for homological features of dimension 0 consists of the red points, and the persistence diagram for homological features of dimension 1 consists of the blue points.</figcaption></figure>

### 2.1.5 Interpretation of persistent homology

Whatever the data type, building a filtration of simplicial complexes on top of the data allows to consider its persistent homology. However, interpretation of the information found in a persistence diagram must be done carefully.

For a point cloud, as we already mentioned, there is a desire to reveal actual topological information. Long living homological features thus have a meaningful interpretation. For other types of data, the interpretation of persistence diagrams must be taken with care. In this article we focus on spike trains: there is no natural idea that comes to mind when we want to define the shape of such data. The goal of persistent homology in this context is only to bring new information to complement other data analysis techniques. The important information in a spike train consists of the frequencies of cofiring neurons. It turns out that a filtration revealing this information can be constructed. Persistence diagrams of a spike train can thus be defined, but hasty conclusions must not be drawn from these diagrams.

## 2.2 Persistence module and persistence diagram

Let us present the algebraic definition of a persistence diagram.

### 2.2.1 Persistence module

**Definition 2.13.** *Persistence module*

A persistence module $V$ over a subset $T \subset \mathbb{R}$ is an indexed family of vector spaces $(V_t)_{t \in T}$ and a doubly indexed family of linear maps $(v_s^t : V_t \to V_s | t \le s)$ which satisfy the composition law $v_s^r \circ v_r^t = v_s^t$ whenever $t \le r \le s$; and where $v_t^r$ is the identity map on $V_t$.
A filtration of simplicial complexes $(K_t)_{t \in \mathbb{R}}$ induces a persistence module, by considering the homology groups. Indeed let $k$ be a nonnegative integer, and consider $V_t = H_k(K_t)$ for $t \in \mathbb{R}$, the homology groups of the spaces $V_t$ for a fixed dimension $k$. Let us denote $i_s^t$ the inclusion maps $K_t \subset K_s$ for $t \le s$, and take $v_s^t = (i_s^t)^*$ the induced linear maps between $H_k(K_t)$ and $H_k(K_s)$. The family of both vectors spaces $(H_k(K_t))_{t \in \mathbb{R}}$ and linear maps $((i_s^t)^*|_{t \le s})$ forms a persistence module.

Let us introduce the notion of tame persistence module.

**Definition 2.14.** *Tame persistence module*

A persistence module $\mathcal{V}$ is tame, if each vector space $V_t$ is finite-dimensional.

**2.2.2 Persistence diagram**

A persistence diagram is a way to represent the decomposition of a persistence module as a direct sum of interval modules, which is defined below.

**Definition 2.15.** *Interval module*

Let $b, d \in \mathbb{R}$, $b \le d$. We define the interval module $\mathbf{I}_{(b,d)}$ over $T \subset \mathbb{R}$, with $b, d \in T$, as the persistence module over $T$ with vector spaces $V_t = \mathbb{F}_2$ for $b \le t \le d$, $V_t = 0$ otherwise; and linear maps $v_s^t = \mathrm{Id}_{\mathbb{F}_2}$ for $b \le t \le s \le d$, $v_s^t = 0$ otherwise$^2$.

<figure><img src="image_7.png" /><figcaption>Figure 6: Interval module $\mathbf{I}_{[b,d]}$.</figcaption></figure>

An interval module $\mathbf{I}_{(b,d)}$ can be interpreted as a homological feature that appears in the filtration at index $b$ and disappears at index $d$. Hence it can also be seen as a point $(b, d) \in \mathbb{R}^2$.

The direct sum of two persistence modules $\mathcal{V}$ and $\mathcal{W}$, denoted as $\mathcal{V} \oplus \mathcal{W}$, is defined by taking the direct sum of the vector spaces $V_t \oplus W_t$, and the “direct sum” of the linear maps $v_s^r \oplus w_s^r$. Let us state the theorem of decomposition of a persistence module as a direct sum of interval modules.

**Theorem 2.1.** ([15], [35], [2])

Let $\mathcal{V}$ be a tame persistence module over $T \subset \mathbb{R}$. Then $\mathcal{V}$ can be decomposed in a unique way (up to a permutation) as a direct sum of interval modules.

<footnote>$^2$More precisely, $b$ (resp. $d$) is the infimum (resp. supremum) of the interval of indices corresponding to $\mathbb{F}_2$ vector spaces. Thus, $b$ and $d$ can be equal to $-\infty$ or $+\infty$.</footnote>
Such a persistence module $V = \bigoplus_i I_{(b_i, d_i)}$ can be seen as a union of points $(b_i, d_i) \in \mathbb{R}^2$: this is what is called a persistence diagram. A persistence diagram is defined as a multiset in $\mathbb{R}^2$, meaning that the same point can appear several times.

## 2.3 The Bottleneck distance

The Bottleneck distance is one of the distances on persistence diagrams. It was introduced in [14]. It takes into account long living features more than short living ones. We take the definition used in [12].

**Definition 2.16.** $l^\infty$ distance

Let $p_1 = (b_1, d_1)$ and $p_2 = (b_2, d_2)$ be two points of $\mathbb{R}^2$. The $l^\infty$ distance is defined as usual if $-\infty < b_1, d_1, b_2, d_2 < +\infty$, that is

$$ \|p_1 - p_2\|_{\infty} = \max(|d_2 - d_1|, |b_2 - b_1|). $$

Otherwise we take the following definitions:

$$ \bullet \ \|( -\infty, d_1 ) - ( -\infty, d_2 ) \|_{\infty} = |d_2 - d_1|, $$

$$ \bullet \ \|( b_1, +\infty ) - ( b_2, +\infty ) \|_{\infty} = |b_2 - b_1|, $$

$$ \bullet \ \| ( -\infty, +\infty ) - ( -\infty, +\infty ) \|_{\infty} = 0, $$

$$ \bullet \ \| p_1 - ( b_2, +\infty ) \|_{\infty} = \| ( b_1, +\infty ) - p_2 \|_{\infty} = +\infty, $$

$$ \bullet \ \| p_1 - ( -\infty, d_2 ) \|_{\infty} = \| ( -\infty, d_1 ) - p_2 \|_{\infty} = +\infty. $$

Let $\Delta$ be the diagonal of $\mathbb{R}^2$, that is $\Delta = \{(r, r), r \in \mathbb{R}\}$. The $l^\infty$ distance of a point $(b, d)$ with $-\infty \le b, d \le +\infty$, to the diagonal $\Delta$ is defined as:

$$ \begin{aligned} \|(b, d) - \Delta\|_{\infty} &:= \inf_{(r,r) \in \Delta} \|(b, d) - (r, r)\|_{\infty} \\ &= \frac{1}{2}|d - b|. \end{aligned} $$

Let us define the notion of $\delta$-matching between two persistence diagrams, which appears in the definition of the Bottleneck distance.

**Definition 2.17.** $\delta$-matching

A $\delta$-matching between two persistence diagrams $dgm^1$ and $dgm^2$ is a collection of pairs $M \subset dgm^1 \times dgm^2$ such that:

* for every $p_1 \in dgm^1$ there is at most one $p_2 \in dgm^2$ such that $(p_1, p_2) \in M$,

* for every $p_2 \in dgm^2$ there is at most one $p_1 \in dgm^1$ such that $(p_1, p_2) \in M$,
and,

* if $(p_1, p_2) \in M$ then $\|p_1 - p_2\|_{\infty} \le \delta$,

* if $p_1 \in \text{dgm}^1$ is unmatched then $\|p_1 - \Delta\|_{\infty} \le \delta$,

* if $p_2 \in \text{dgm}^2$ is unmatched then $\|p_2 - \Delta\|_{\infty} \le \delta$.

The Bottleneck distance is then defined as follows. An example is shown in Figure 7.

**Definition 2.18. Bottleneck distance**

The Bottleneck distance $d_B(\text{dgm}^1, \text{dgm}^2)$ between two persistence diagrams is defined as

$$d_B(\text{dgm}^1, \text{dgm}^2) = \inf \{\delta | \text{there exists a } \delta\text{-matching between } \text{dgm}^1 \text{ et } \text{dgm}^2\}.$$

<figure><img src="image_8.png" /><figcaption>Figure 7: The Bottleneck distance between two persistence diagrams. The diagram $\text{dgm}^1$ consists of the two black points and the diagram $\text{dgm}^2$ consists of the blue point. The $\delta$-matching realizing $d_B(\text{dgm}^1, \text{dgm}^2)$ is displayed with thick dotted lines. The Bottleneck distance $d_B(\text{dgm}^1, \text{dgm}^2)$ between the diagrams is represented in red and denoted by $d_B$.</figcaption></figure>

## 2.4 Interleaving and stability

Interleaved persistence modules intervene in Theorem 2.2. This theorem or its refinements, is the core of most stability proofs of persistent homology, for example the ones stated in Remark 2.1. We apply it in section 4 to prove the stability of the Curto-Itskov filtration.

**Definition 2.19. $\varepsilon$-interleaving**

Let $\varepsilon \ge 0$. Two persistence modules $\mathcal{U}$ and $\mathcal{V}$ are said to be $\varepsilon$-interleaved if there exist maps $\phi_t : U_t \to V_{t+\varepsilon}$ and $\psi_t : V_t \to U_{t+\varepsilon}$ defined for all $t$, such that the following diagrams commute for all eligible parameter values; that is, for all $t \le s$:
<figure><img src="image_9.png" /><figcaption>Figure 8: Interleaving diagrams.</figcaption></figure>

**Theorem 2.2.** [10], Theorem 4.4

Let $\varepsilon \ge 0$. If $U$ and $V$ are tame persistence modules that are $\varepsilon$-interleaved then the Bottleneck distance between the diagrams satisfies the bound

$$d_B(dgm(U), dgm(V)) \le \varepsilon.$$

**Remark 2.1.** Stability of persistence diagrams.

Let us give some examples of existing stability results.

The authors of [14] proved the stability theorem of persistence diagrams constructed from sublevel sets of functions. Let $X$ be a metric space and $f, g : X \to \mathbb{R}$ be two functions and let $dgm(f)$ and $dgm(g)$ be their respective persistence diagrams. Let $d_B$ be the Bottleneck distance. Then under some conditions on $X$, $f$ and $g$ we have:

$$d_B(dgm(f), dgm(g)) \le \|f - g\|_{\infty}.$$

The stability of persistence diagrams constructed from Čech filtrations with respect to the Hausdorff distance is a consequence of this theorem ([3], Corollary 11.29). The authors of [9] proved the stability with respect to the Gromov-Hausdorff distance, which is a stronger result: Let $X$ and $Y$ be two finite point clouds (more generally, two totally bounded metric spaces), let $dgm(X)$ and $dgm(Y)$ be their respective persistence diagrams obtained from the Čech filtration, let $d_{GH}$ be the Gromov-Hausdorff distance. Then

$$d_B(dgm(X), dgm(Y)) \le 2d_{GH}(X, Y).$$

In Section 4 we bound the Bottleneck distance between persistence diagrams obtained from spike trains with the Curto-Itskov filtration, by a distance between spike trains. In this way, we talk about stability of the Curto-Itskov filtration/stability of persistence diagrams for spike trains.

# 3 The Curto-Itskov filtration

We define the Curto-Itskov filtration for spike trains. First, let us set what we mean by spike train.
**Definition 3.1.** *Spike train*

A spike train is a collection of binary time series. Each time series corresponds to a recorded neuron, and at each time, the series at a time $t_j$ includes a 1 if the neuron emits a spike at $t_j$, and a 0 otherwise.

A spike train example with 3 neurons is displayed in Figure 9. The term “spike train” is usually employed for a unique binary time series corresponding to one neuron, but for convenience in this article this term is employed even for a collection of such time series.

As explained in Section 2, to build a filtration we must think about what information we want to focus on. For spike train data, an interesting information is the cofiring frequencies of any set of neurons.

## 3.1 Frequencies in a spike train

Let $\mathcal{S}$ be a spike train. Let us denote by $N_1, \dots, N_p$ the neurons of $\mathcal{S}$. Let us fix a time interval length, denoted as $\delta t$.

For each integer $1 \le m \le p$, we define the occurence of a set a $m$ neurons $(N_{i_1}, \dots, N_{i_m})$, as the number of co-spikes of $(N_{i_1}, \dots, N_{i_m})$ in a window of length $\delta t$. More presicely, first we cut the spike train in consecutive windows of length $\delta t$; let us say that we obtain $\ell$ windows. The occurence of $(N_{i_1}, \dots, N_{i_m})$ is then defined as the number of these windows within which each neuron of $(N_{i_1}, \dots, N_{i_m})$ has at least one spike. In practice, we choose $\delta t$ such that in the considered biological context, it is smaller than the smallest typical inter-spike interval and greater than the largest typical cospike interval. Finally, the frequency of $(N_{i_1}, \dots, N_{i_m})$ is defined as its occurence, divided by $\ell$. This frequency is denoted by $\mathbf{fr}(N_{i_1}, \dots, N_{i_m})$. To describe the cofiring behaviour of all sets of neurons of the spike train $\mathcal{S}$, let us denote the frequency vector of $\mathcal{S}$ as $F_S(N_1, \dots, N_p)$, and define it as the vector that includes the frequencies of all sets of neurons, in the following order:

$$F_S(N_1, \dots, N_p) = (\mathbf{fr}(N_i)_{1 \le i \le p}, \mathbf{fr}(N_i, N_j)_{1 \le i < j \le p}, \mathbf{fr}(N_i, N_j, N_\ell)_{1 \le i < j < \ell \le p}, \dots).$$

**Remark 3.1.** Let us highlight some facts about this definition, that are important to understand the distance on spike trains that is defined in Subsection 4.1.

1. Sets of neurons of frequency 0 are taken into account in the frequency vector.

2. It is possible to define $F_S(N_1, \dots, N_p)$ for every $p' > p$, by adding zeros in the vector $F_S(N_1, \dots, N_p)$. This will allow to compare frequency vectors of spike trains that have different number of neurons.

**Example 3.1.** As a toy example, let us consider the spike train $\mathcal{S}^1$ of Figure 9. Let us call the neurons $N_1, N_2, N_3$ by numbering them vertically from top to bottom. Imagine that for a given $\delta t$ we obtain 12 consecutive windows of length $\delta t$. Then $F_S^1$ the frequency vector
of $\mathcal{S}^1$, includes, in this order, the frequencies of $N_1, N_2, N_3, (N_1, N_2), (N_1, N_3), (N_2, N_3)$ and $(N_1, N_2, N_3)$. Thus it is given by:

$$F_{\mathcal{S}^1}(N_1, N_2, N_3) = \left( \frac{7}{12}, \frac{8}{12}, \frac{8}{12}, \frac{3}{12}, \frac{4}{12}, \frac{6}{12}, \frac{2}{12} \right) \\ \approx (0.58, 0.67, 0.67, 0.25, 0.33, 0.5, 0.17).$$

<figure><img src="image_10.png" /><figcaption>Figure 9: A spike train example.</figcaption></figure>

## 3.2 The Curto–Itskov filtration

The Curto-Itskov filtration was introduced in the article [18], itself inspired from the work of Curto and Itskov in [16]. Our definition is the same than the one in [18], with the modifications that we consider frequencies and not occurences, and that we consider the frequencies of neurons that cofire within a small time interval rather than at the exact same time.

**Definition 3.2.**

For each real number $t \in [-1, 0]$, let $K_t$ be the abstract simplicial complex defined by:

$$\{N_{i_1}, \dots, N_{i_m}\} \in K_t \text{ if and only if } \mathbf{fr}(N_{i_1}, \dots, N_{i_m}) \geq -t.$$

The sequence of $K_t$'s for increasing values $-1 \le t \le 0$, forms a filtration of simplicial complexes (this is easily verified). This filtration is called the Curto-Itskov filtration. We formally chose to index the Curto-Itskov filtration by negative numbers, in order to have a filtration indexed by an increasing parameter. A more natural way to understand the definition is to think of the values $f = -t$, $f \in [0, 1]$ as frequency levels: for $f \in [0, 1]$, we have

$$\{N_{i_1}, \dots, N_{i_m}\} \in K_{-f} \text{ if and only if } \mathbf{fr}(N_{i_1}, \dots, N_{i_m}) \geq f.$$
**Remark 3.2.** The authors of [18] chose the name “Curto-Itskov filtration” because Curto and Itskov in [16], defined the following simplicial complex K:

$$ \{N_{i_1}, \dots, N_{i_m}\} \in K \text{ if and only if } \mathrm{occ}(N_{i_1}, \dots, N_{i_m}) \ge 1, $$

where $\mathrm{occ}(N_{i_1}, \dots, N_{i_m})$ is the number of cospikes of $(N_{i_1}, \dots, N_{i_m})$ in the whole spike train. Curto and Itskov studied the homology of this simplicial complex K, for place cells spike trains.

**Example 3.2.** The Curto-Itskov filtration of the spike train $S^1$ of Example 3.1, is represented in Figure 10.

<figure><img src="image_11.png" /><figcaption>Figure 10: The Curto-Itskov filtration of the spike train $S^1$ of Example 3.1.</figcaption></figure>

For every frequency level f such that $t = -f < -8/12$, the corresponding simplicial complex $K_t$ is empty, and for every frequency level f such that $t = -f \ge -2/12$, $K_t = K_{-2/12}$ is the filled triangle with all its edges and vertices.

The Curto-Itskov filtration of simplicial complexes induces a persistence module, by considering the homology groups. Moreover we have the following fact.

**Fact 3.1.** A persistence module built from the Curto-Itskov filtration is tame. Indeed we always consider finitely many neurons, so every simplicial complex in the filtration is finite. Every homology group in the persistence module is thus finite-dimensional.

Theorem 2.1 thus ensures that the persistence diagram of a spike train is well-defined when we consider the Curto-Itskov filtration.

**Example 3.3.** Continuation of Examples 3.1 and 3.2

In Figure 11 we display the persistence diagram of the Curto-Itskov filtration of the spike train $S^1$ of Example 3.1.

# 4 Stability of the Curto-Itskov filtration

Stability is a fundamental theoretical guarantee in persistent homology. In our setting, we want to bound the Bottleneck distance between the persistence diagrams, by some distance between the spike trains.

In this section we define a meaningful distance on the space of spike trains, and deduce a stability property.
<figure><img src="image_12.png" /><figcaption>Figure 11: The persistence diagram of the Curto-Itskov filtration of the spike train $S^1$ of Example 3.1.</figcaption></figure>

## 4.1 Distance on spike trains

First, let us point out that until now, in a spike train the neurons are always labeled arbitrarily. In fact, the neurons are interchangeable: the quantity of interest in a spike train is the frequencies of sets of cofiring neurons, whatever the neuron labels. The distance that is defined takes this fact into account by allowing all neuron permutations.

Let $S^1$ and $S^2$ be two spike trains. Let $p$ be the maximum between the numbers of neurons in $S^1$ and the number of neurons in $S^2$. The neurons for both spike trains are denoted by $N_1, \dots, N_p$, even if there are not necessarily the same neurons in $S^1$ and $S^2$.

To define the distance between $S^1$ and $S^2$, the main idea is to take the $\ell^\infty$ distance between their frequency vectors. To solve the problem of the neuron arbitrary ordering, the infimum over all neuron permutations is taken.

**Definition 4.1.** Distance between spike trains

The distance between $S^1$ and $S^2$, denoted as $d(S^1, S^2)$, is defined as

$$d(S^1, S^2) := \min_{\pi \in \mathcal{S}_p} \|F_{S^1}(N_1, \dots, N_p) - F_{S^2}(N_{\pi(1)}, \dots, N_{\pi(p)})\|_{\infty},$$

where $\mathcal{S}_p$ denotes the permutation group of $\{1, \dots, p\}$.

Because of Fact 2 of Remark 3.1, $F_{S^1}(N_1, \dots, N_p)$ and $F_{S^2}(N_1, \dots, N_p)$ are well-defined. The distance $d$ satisfies almost all distance axioms. Indeed, deriving from the properties of the infinity norm for vectors, it can be verified that $d$ is non negative, symmetric and satisfies the triangle inequality. The point separation axiom is not satisfied, indeed two different spike trains can have 0 distance if, up to neuron permutation, they have the same frequency vector. This is not a problem since for us, two such spike trains are considered identical.

**Remark 4.1.** This definition serves theoretical purposes. The actual computation of the distance is thus not discussed.
**Example 4.1.** Let $S^1$ be the spike train of Example 3.1, and let $S^2$ be the spike train displayed in Figure 12. Imagine again that for a given $\delta t$ we obtain 12 consecutive windows of length $\delta t$. Denoting the 4 neurons by $N_1, N_2, N_3, N_4$, the frequency vector of $S^2$ includes, in this order, the frequencies of:

$N_1, N_2, N_3, N_4,$
$(N_1, N_2), (N_1, N_3), (N_1, N_4), (N_2, N_3), (N_2, N_4), (N_3, N_4),$ 
$(N_1, N_2, N_3), (N_1, N_2, N_4), (N_1, N_3, N_4), (N_2, N_3, N_4), \text{ and}$ 
$(N_1, N_2, N_3, N_4).$

It is given by:

$$F_{S^2}(N_1, N_2, N_3, N_4) = \left( \frac{8}{12}, \frac{8}{12}, \frac{1}{12}, \frac{6}{12}, \frac{6}{12}, \frac{1}{12}, \frac{4}{12}, \frac{1}{12}, \frac{3}{12}, 0, \frac{1}{12}, \frac{3}{12}, 0, 0, 0 \right).$$

Let us compute the distance between $S^1$ and $S^2$. The frequency vector $F_{S^1}(N_1, N_2, N_3, N_4)$ can be defined by adding an additional neuron $N_4$ with 0 frequency. This is then easily verified that the permutation $\pi$ that minimizes

$$\|F_{S^1}(N_1, N_2, N_3, N_4) - F_{S^2}(N_{\pi(1)}, N_{\pi(2)}, N_{\pi(3)}, N_{\pi(4)})\|_{\infty},$$

is the permutation $\pi$ that maps 4 to 3, 3 to 1, 1 to 4, and 2 to 2, giving $d(S^1, S^2) = 1/12$.

<figure><img src="image_13.png" /><figcaption>Figure 12: Another spike train example.</figcaption></figure>

## 4.2 Stability theorem

Let us state the stability theorem of the Curto-Itskov filtration.

**Theorem 4.1.** *Stability of the Curto-Itskov filtration*

Let $S^1$ and $S^2$ be two spike trains. Let $dgm(S^1)$ and $dgm(S^2)$ be the corresponding persistence diagrams built from the Curto-Itskov filtration. Then

$$d_B(dgm(S^1), dgm(S^2)) \le d(S^1, S^2).$$
To prove this theorem, we use Theorem 2.2. Our proof is mostly an adaptation of the stability proof of persistent homology of sublevel sets of functions, of [14] and [3] Proposition 11.26.

Proof. of Theorem 4.1

Let $\mathcal{S}^1$ and $\mathcal{S}^2$ be two spike trains, let $(K_t^1)_{-1\le t\le 0}$ and $(K_t^2)_{-1\le t\le 0}$ be their respective Curto-Itskov filtrations, and let $\mathbb{V}^1$ and $\mathbb{V}^2$ be the induced persistence modules, for a fixed dimension $k$. First, as stated in Fact 3.1, these persistence modules are tame. Let $\varepsilon = d(\mathcal{S}^1, \mathcal{S}^2)$. Let us now prove that $\mathbb{V}^1$ and $\mathbb{V}^2$ and $\varepsilon$-interleaved.

Let $p$ be the maximum between the numbers on neurons in $\mathcal{S}^1$ and the number of neurons in $\mathcal{S}^2$, and let us call $N_1, \dots, N_p$ the neurons for both spike trains.

Let $\pi \in \mathfrak{S}_p$ be a permutation that achieves the minimum

$$ \min_{\pi \in \mathfrak{S}_p} \| F_{\mathcal{S}^1}(N_1, \dots, N_p) - F_{\mathcal{S}^2}(N_{\pi(1)}, \dots, N_{\pi(p)}) \|_{\infty}. $$

Then, we have $\varepsilon = \|F_{\mathcal{S}^1}(N_1, \dots, N_p) - F_{\mathcal{S}^2}(N_{\pi(1)}, \dots, N_{\pi(p)})\|_{\infty}$ for this permutation $\pi$.

Let $f \in [0, 1]$ be a frequency level, and $t = -f$. Let $\{N_{i_1}, \dots, N_{i_m}\}$ be a $(m-1)$-simplex in $K_t^1$. Let us show that $\{N_{\pi(i_1)}, \dots, N_{\pi(i_m)}\}$ is included in $K_{t+\varepsilon}^2$.

Let us note that the frequency of $(N_{i_1}, \dots, N_{i_m})$ in $\mathcal{S}^1$, is the $j$-th coordinate of the vector $F_{\mathcal{S}^1}(N_1, \dots, N_p)$ for some $j$, and the frequency of $(N_{\pi(i_1)}, \dots, N_{\pi(i_m)})$ in $\mathcal{S}^2$ is the $j$-th coordinate of $F_{\mathcal{S}^2}(N_{\pi(1)}, \dots, N_{\pi(p)})$ for the same $j$. Since $\{N_{i_1}, \dots, N_{i_m}\} \in K_t^1$, by definition the frequency of $(N_{i_1}, \dots, N_{i_m})$ in $\mathcal{S}^1$ is greater than $f = -t$. It implies that the frequency of $(N_{\pi(i_1)}, \dots, N_{\pi(i_m)})$ in $\mathcal{S}^2$ is greater than $f - \varepsilon = -(t + \varepsilon)$. Indeed, for two vectors $v$ and $w$ and an index $j$ we have

$$ \begin{aligned} v_j &= v_j - w_j + w_j \\ &\ge -\|w - v\|_{\infty} + w_j, \end{aligned} $$

Thus if $w_j \ge f$ then $v_j > f - \|w - v\|_{\infty}$.

Taking $v = F_{\mathcal{S}^1}(N_1, \dots, N_p)$, $w = F_{\mathcal{S}^2}(N_{\pi(1)}, \dots, N_{\pi(K)})$, $\varepsilon = \|v - w\|_{\infty}$, and the index $j$ introduced above, this gives the fact that the frequency of $(N_{\pi(i_1)}, \dots, N_{\pi(i_m)})$ in $\mathcal{S}^2$ is greater than $f - \varepsilon = -(t + \varepsilon)$. In other words, $\{N_{\pi(i_1)}, \dots, N_{\pi(i_m)}\}$ is included in $K_{t+\varepsilon}^2$.

We showed that if a simplex $\{N_{i_1}, \dots, N_{i_m}\}$ is in $K_t^1$ then $\{N_{\pi(i_1)}, \dots, N_{\pi(i_m)}\}$ is in $K_{t+\varepsilon}^2$. Since $\{N_{\pi(i_1)}, \dots, N_{\pi(i_m)}\}$ and $\{N_{i_1}, \dots, N_{i_m}\}$ are isomorphic simplicial complexes, the map between $K_t^1$ and $K_{t+\varepsilon}^2$ which maps each simplex $\{N_{i_1}, \dots, N_{i_m}\} \in K_t^1$ to the simplex $\{N_{\pi(i_1)}, \dots, N_{\pi(i_m)}\} \in K_{t+\varepsilon}^2$, acts as an inclusion. Up to a neuron permutation, it is as if $K_t^1$ was included in $K_{t+\varepsilon}^2$. This “inclusion” map induces a linear map between the homology groups $H_k(K_t^1)$ and $H_k(K_{t+\varepsilon}^2)$. Let us call $(\phi_f)_{f \in \mathbb{N}}$ the collection of these linear maps. Symetrically, take $(\psi_f)_{f \in \mathbb{N}}$ the maps induced by “inclusion” of $K_t^2$ in $K_{t+\varepsilon}^1$.
It is immediate that $\phi$ and $\psi$ provide an $\varepsilon$-interleaving between $V^1$ and $V^2$. Indeed, they are induced by “inclusion”, and so are by definition the maps of the persistence modules $V^1$ and $V^2$. Thus, every commutation property needed to have an interleaving, is satisfied.

The persistence modules $V^1$ and $V^2$ are thus tame and $\varepsilon$-interleaved for the value $\varepsilon = d(S^1, S^2)$. Applying Theorem 2.2 concludes the proof. $\square$

# 5 Simulations

In this section we apply the Curto-Itskov filtration on simulated spike trains. The goal of these simulations is to apprehend the effect of the presence of sets of co-firing neurons on the persistence diagram.

## 5.1 Simulation principle

We simulate spike trains that are similar to the real data that is used in Section 6, i.e. extracellular recordings of the first olfactory relay of the locust *Schistocerca americana*. In this real dataset, there are 7 neurons, the original time window is 1/15 milliseconds and the recording time length is between 10 and 15 minutes. To apply the Curto-Itskov filtration, taking $\delta t = 10$ milliseconds is an appropriate choice for this context (this is close to the observed smallest inter-spike intervals of any given neuron).

We thus simulate spike trains with 7 neurons, for a “recording length” of 15 minutes. Since the goal is to apply the Curto-Itskov filtration with $\delta t = 10$ milliseconds, there is no need to simulate an actual whole spike train with 13 500 000 time steps of 1/15 millisecond. Instead, we directly simulate spikes and cospikes in 90 000 windows, which corresponds to the number of windows of 10 milliseconds in a 15 minutes recording.

The simulation principle is the following. For each neuron labeled from 0 to 6, we respectively take the firing rates as $\frac{1000}{90000}$, $\frac{2000}{90000}$, $\frac{2500}{90000}$, $\frac{3000}{90000}$, $\frac{4000}{90000}$, $\frac{4500}{90000}$ and $\frac{5000}{90000}$. These firing rates are in the range of those observed in the real datasets. We first let each neuron fire with respect to its firing rate, then we potentially reinforce cofirings of some sets of neurons. To reinforce the coactivity of a given set of $m$ neurons ($N_{i_1}, \dots, N_{i_m}$), we proceed as follows. For each window, we pick a reference neuron among ($N_{i_1}, \dots, N_{i_m}$) that fires in this window, say $N_{i_k}$. Then we have each of the other neurons of the set spike in this window, with probability 0.4. Finally, to be a bit more realistic we add some inhibition. To do so, we choose a fixed set of neurons receiving an input from a common *unobserved* inhibitory neuron$^3$, here (2, 3, 4, 5), and in each window we remove the possible spike of each of these neurons, with probability $\frac{10000}{90000}$. We chose (2, 3, 4, 5) as inhibited neurons because this gives results that best coincide with the real data.

<footnote>$^3$Inhibitory neurons cannot be recorded in this preparation with the recording technique used in these experiments.</footnote>
**Remark 5.1.** It is important to note that we do not simulate neuronal activity at the anatomic graph level, but we only simulate spike trains. We remind that in real datasets, neuronal activity is recorded with electrodes and the obtained signal is a mixture of the activity of several neurons. Then, with a spike sorting procedure$^4$, the neurons that contributed to the recording are identified as well as their spike times. The final result is a spike train. Thus, the way we simulate co-firing activity of sets of neurons directly on the spike train is artificial. It is different from simulating a neuronal population with co-firing sets of neurons from a biophysically based or intensity based model [20] and working on the obtained spike trains. Our simulations only produce spike trains that could have obtained if sets of neurons are really coactive at the graph level (directly coactive, or coactive via another common neuron that we do not observe). In the following subsections we often talk about “sets of co-firing neurons” for the sake of simplicity, but there must be no confusion: we only mean sets of neurons for which we have reinforced the cospikes in the spike train simulations.

After simulating the spike trains, we apply the Curto-Itskov filtration with $\delta t = 10$ milliseconds and we apply persistent homology with coefficients in $\mathbb{F}_2$.

**Remark 5.2.** About the field choice

For both simulated and real data, we tried and used different finite field choices when applying persistent homology, but they provide the same results. Thus we only present the results with $\mathbb{F}_2$.

# 5.2 First results and interpretation

## 5.2.1 Introduction

We present here the result of our simulations. We study the effect of coactive sets of neurons on the persistence diagrams obtained with the Curto-Itskov filtration. The code for simulations and application to real data is available at https://gitlab.math.unis-tra.fr/martineau/persistent-homology-of-spike-trains-with-the-curto-itskov-filtration. For each simulation we only show one persistence diagram, but this is representative of the several simulations that we carried out for each case, and that are presented in the notebook.

First let us emphasize that there is always one persistence diagram per dimension, but if there is no confusion we speak of “the persistence diagram” to talk about the collection of persistence diagrams of each dimension.

Let us remind that it is possible to have multiple points with the same coordinates, but visually we can only see one point. We should thus display the point coordinates of a diagram, in addition to displaying it, to have access to complete information.

<footnote n="4">⁴See: https://c_pouzat.gitlab.io/spike-sorting-the-diy-way/.</footnote>
We also draw attention on the way we interpret the persistence diagrams obtained with the Curto-Itskov filtration. The usual interpretation of a persistence diagram, with the Čech filtration for instance, considers that points far away from the diagonal correspond to topological structure and points close to the diagonal are due to noise, and that we can thus recover the topology of a point cloud with the points far away from the diagonal. In our context, as mentioned in Subsection 2.1.5, we cannot interpret the diagrams in this way. Therefore, to evaluate the actual effect of reinforcing the activity of a set of neurons, the main point is to compare a given diagram to the diagrams obtained in the case where all neurons are independent.

Finally, let us note that to see a homological feature of dimension $k \in \mathbb{N}$ in the Curto-Itskov filtration, there must be sets of $(k+1)$ neurons that cofire within a window of length $\delta t$. For instance, a 1-dimensional hole appears if a collection of edges in the simplicial complex forms an empty triangle, or an empty square, an empty pentagon etc. Moreover, let us note that with seven neurons, we can have homological features of dimension 5 at the most.

**Remark 5.3.**

We draw attention on the interpretation of the absence of k-dimensional holes. It can be due to two reasons :

* either there is no set of $(k+1)$ co-firing neurons, meaning that all the frequency of sets of $(k+1)$ neurons are null;

* or there are sets of $(k+1)$ cofiring neurons, but either there are not enough of them to form a k-dimensional hole, or there are enough but they are arranged in the simplicial complex in a way that they do not form a k-dimensional hole.

### 5.2.2 Case of independent neurons

Spike trains are simulated with only independent neurons. Each neuron fires with respect to its own firing rate, without interacting with other neurons, but inhibition is still present as before. The obtained diagrams look like the one displayed in Figure 13 A. In Figure 14 we also plot the point coordinates, as well as a zoom on dimensions 0 and 1.

## Results

### Points in dimension 0

If we label the neurons from 0 to 6 and denote by $f_0, \dots, f_6$ their firing rates, then at each negative frequency value ($-f_i$) a connected component (the 0-simplex corresponding to neuron $i$) appears in the filtration. The connected components in the filtration thus correspond to the different neurons, hence the presence of 7 points in the persistence diagram of dimension 0. Throughout the filtration connected components merge and in the end there is only one connected component that remains until infinity: there is one
<figure><img src="image_14.png" /><figcaption>Figure 13: **A:** Persistence diagram of a spike train with independent neurons. **B:** Persistence diagram of a spike train with one pair of coactive neurons: (1, 2). **C:** Persistence diagram of a spike train with two disjoint pairs of coactive neurons: (0, 1) and (2, 6). **D:** Persistence diagram of a spike train with a set of three coactive neurons: (1, 2, 3). **E:** Persistence diagram of a spike train with 6 coactive pairs: (1, 0), (0, 4), (4, 3), (2, 3), (3, 6), (6, 5). **F:** Persistence diagram of a spike train with 8 coactive pairs: (1, 0), (0, 4), (4, 3), (2, 3), (3, 6), (6, 5), (1, 2), (4, 5). **G:** Persistence diagram of a spike train with the set of all of the 7 neurons, being coactive.</figcaption></figure>
point at $y = +\infty$ in the diagram.

This remark holds for all the persistence diagrams in our simulations, since we always
have 7 neurons.

**Structure**

In the diagrams of Figure 13 A and Figure 14 we can see 15 points in dimension 1, meaning that 1-dimensional holes appear throughout the filtration. There are also 15 points in dimension 2, with 14 points at $y = +\infty$. As mentioned in Subsection 2.1.5, there is no direct geometric interpretation of the persistence diagram in our context, but the presence of many points in dimension 2 including points at infinity means that there is some topological structure.

It turns out that all the diagrams in our simulations also present a lot a topological
structure, with homological features at least up to dimension 2.

**Importance of the frequency vector**

To have complete information, one must look at the frequency vector as well as the
persistence diagram. Indeed, here, for instance, there is no 3-dimensional hole, but by
only looking at the diagram we do not know whether it is because there are no sets of 4
cofiring neurons, or because there are sets of 4 cofiring neurons but they just do not form
a 3-dimensional hole. If we look at the frequency vector we obtain this information:
here it turns out that here we have one set of 4 neurons that has a nonzero frequency
(we do not show the frequency vector here).

**5.2.3 Reinforcing one pair of neurons**

Here, the coactivity of one pair of neurons is reinforced as explained in Subsection 5.1.
We reinforce the pair (1, 2) (we recall that the neurons are labeled from 0 to 6). We
show the resulting diagram<fnref n="5" /> in Figure 13 B.

**Results**

Compared to the case of independent neurons, the diagram in Figure 13 B clearly
shows one point of dimension 0 lower than the others. When zooming, we can see that
there are five points with a death time around -0.0025 as in the independent case, and
there is one point that dies around -0.018. It corresponds to a connected component
that died earlier. This is explained by the fact that reinforcing one pair of neurons, makes
an edge between the two neurons appear early in the filtration. Thus, two connected
components, the 0-simplices corresponding to the two neurons of the pair, are merged
when the edge appears. It makes a connected component die: indeed we recall that in

<footnote n="5">⁵In the notebook you can zoom on the diagram and display all point coordinates but we do not show these results here.</footnote>
<figure><img src="image_15.png" /><figcaption>Figure 14: Complement plots to the independent case displayed in Figure 13 A. Top: Point coordinates of the persistence diagram. Bottom: A zoom on persistence diagrams of dimension 0 and 1 of Figure 13 A.</figcaption></figure>

persistent homology when two homological features are merged, we say that the one that appeared last in the filtration, dies. Thus reinforcing the coactivity of a pair of neurons, shifts a point of dimension 0 to the bottom of the points of the independent case.

### 5.2.4 Reinforcing two pairs of neurons

We reinforce the coactivity of two disjoint pairs of neurons: (0, 1) and (2, 6). We show the resulting persistence diagram in Figure 13 C.

#### Results

The most significative effect is the fact that there are two points in dimension 0 that have a death time lower than the typical death value in the independent case. The two reinforced pairs indeed make two edges appear early in the filtration thus making two connected components die early.

### 5.2.5 Reinforcing a set of three neurons

We reinforce the coactivity of a set of three neurons : (1, 2, 3). A typical diagram is displayed in Figure 13 D.

#### Results

* We can see two points in dimension 0 with early death times. Indeed, reinforcing (1, 2, 3) also reinforces the pairs (1, 2), (2, 3) and (1, 3), thus making two connected
<figure><img src="image_16.png" /><figcaption>Figure 15: The 2 reinforced structures of Subsubsection 5.2.6. This will be further referred to as “Figure 15 left” and “Figure 15 right”.</figcaption></figure>

components die early (the three pairs form a triangle so adding the third edge does not make a third connected component die).

* A point in dimension 1 is lower and shifted to the left, compared to the independent case. It means that a 1-dimensional hole is born and dead earlier than in the independent case. The early birth time is explained by the fact that we indirectly reinforced 3 pairs that form an empty triangle, making a 1-dimensional hole appear early. This empty triangle is filled early since we reinforced the coactivity of the set of its three vertices, hence the early death time.

* There are homological features in dimension 3. More precisely there are three points, including one at $y = \infty$. It can be explained by the fact that reinforcing three neurons among seven, has the consequence that quite a lot of structure appears.

### 5.2.6 Reinforcing pairs that form or do not form 1-dimensional holes

In this part we reinforce some sets of neurons to create a certain structure, and evaluate the impact on the persistence diagram. We run two experiments. First we reinforce the coactivity of six pairs: (1,0), (0,4), (4,3), (2,3), (3,6), (6,5). Then, we reinforce these same six pairs and add two other ones: (1,2) and (4,5). In the first case, we thus reinforce the structure drawn in Figure 15, left: we have six edges, that do not form any higher-dimensional structure. In the second case we reinforce the structure of Figure 15, right: it is the same as the left one but two 1-dimensional holes are now present. Let us highlight again that we do not reinforce these structures at the anatomic graph level. Here we only reinforce cospikes of some pairs of neurons such that in the filtration, the structures of Figure 15 emerge. The resulting diagrams are displayed in Figure 13 E and F.

## Results
* In both cases there are 6 points in dimension 0 that have earlier death times than in the independent case. This is expected since in the first simulation we reinforced 6 pairs. In the second simulation, adding two other pairs does not affect more the connected components. Indeed since there are only 7 neurons, reinforcing 6 or more pairs makes all connected components die early to form only one connected component that lives forever.

* The points in dimension 0 also have earlier birth times than in the independent case, they are shifted to the left. This is explained by the fact that reinforcing a pair, also increases the firing rates of the neurons of the pair, and here some neurons belong to several points of the reinforced structure. Thus it increases a lot the firing rates of these neurons, resulting in an early birth time.

* In Figure 13 F, there are four points in dimension 1 that are shifted to the left. We mainly reinforce two holes of dimension 1, making them appear early in the filtration hence the shift to the left. We must have indirectly reinforced a bit other structures too, explaining the four points. However, if we think in terms of points far away from the diagonal that represent significative structures, here there are two points that are really far away from the diagonal. In a more classical persistent homology interpretation, we can say that there are two significant holes of dimension 1, and indeed it coincides with the structure that we reinforce.

* In Figure 13 E and Figure 13 F, there are homological features of dimension up to 4. In both cases we reinforce more sets of neurons than in the previous simulations thus it creates higher dimensional features. However, these features have low lifetime : they are really close to the diagonal, meaning that they are not to be interpreted as significant structure. Moreover, there is no point at $y = +\infty$ in dimension 2, as opposed to the case of independent neurons (and also to all the previous simulations). It seems that reinforcing a really specific structure like we do here, makes the structure of the independent case disappear and we only see information about the structure that we reinforce. In the previous simulations this phenomenon does not happen, likely because reinforcing only one or two pairs or a set of three neurons must not generate enough structure to kill the original structure of the independent case.

### 5.2.7 Reinforcing the set of the seven neurons

Finally we simulate the case where all of the seven neurons are coactive together. This is unrealistic but we only want to appreciate the impact of such a coactivity on the persistence diagram. We display the resulting diagram in Figure 13 G. First let us note that we indirectly reinforce the coactivity of sets of two, three, four, five and six neurons. We thus expect to see homological features of all possible dimensions appear in the filtration. However, we do not expect any significant hole of any dimension to emerge, since all possible simplices are reinforced.
## Results

We can see homological features up to dimension 5. These features have small lifetimes and there are no points at $y = +\infty$ (except of course in dimension 0). There is thus no significant homological feature; and once again we lost the structure of the independent case. We can also say that at least the points of dimension 0, 1 and 2 are lower than in the independent case, showing the reinforced sets of 2 neurons, 3 neurons, etc...

### 5.2.8 Conclusion

The previous simulations showed that reinforcing the coactivity of some sets of neurons has a direct and very visual impact on the persistence diagrams.

## 5.3 Clustering the simulated diagrams

In Subsection 5.2, we run different simulations and we looked at the resulting persistence diagrams. In this subsection we want to assess if we can distinguish diagrams from different simulations with an automatic clustering method.

To cluster persistence diagrams, what is called a Betti curve is computed. The definition is the following:

**Definition 5.1.** Betti curve

The Betti curve $\beta$ of a persistence diagram dgm (for a given dimension), is defined as:

$$ \beta(t) = \sum_{x=(b,d) \in dgm} \mathbf{I}(b \le t \le d), \quad t \in \mathbb{R}. $$

In other words, the Betti curve is a function from $\mathbb{R}$ to $\mathbb{R}$ that counts at each time $t$, the number of homological features that are alive at time $t$.

**Remark 5.4.** In the following, we remove points at infinity to compute Betti curves. More precisely we compute $\beta(t) = \sum_{x=(b,d) \in dgm, d<\infty} \mathbf{I}(b \le t \le d)$, $t \in \mathbb{R}$.

Using Betti curves is really common in persistence homology. Other transformations of persistence diagrams are also possible, but we only compute Betti curves here since they are sufficient to discriminate the diagrams.

**Clustering the seven diagrams of Subsection 5.2**

For each of the seven cases of Subsection 5.2, five simulations are run. We thus obtain 35 persistence diagrams. We compute the Betti curves of dimension 0 (meaning, the Betti curves of the persistence diagrams in dimension 0). The corresponding results are displayed in Figure 16. Visually the seven groups are clearly distinguished. Then we apply the K-means algorithm to this collection of curves. We obtain a perfect clustering: each Betti curve being correctly labeled.
<figure><img src="image_17.png" /><figcaption>Figure 16: Betti curves in dimension 0 of the seven simulation cases of Subsection 5.2.</figcaption></figure>

**Remark 5.5.** Using Betti curves of dimension 1 also gives a perfect clustering.

# 6 Application to actual locust spike trains

In this section, we apply the Curto-Itskov filtration to a locust dataset, which is publicly available [32] and whose spike sorting is fully detailed at https://christophe-pouza.t.github.io/zenodo-locust-datasets-analysis/.

## 6.1 Data description

The dataset used in the sequel, locust20010214, consists of spike trains obtained from extracellular recordings of the first olfactory relay, the antennal lobe, of the locust, *Schistocerca americana*. Two neuronal types are found in this brain region [24]: circa 800 excitatory *projection neurons*, and circa 300 inhibitory *local neurons*. The projection neurons generate “classical”, sodium dependent, action potentials, while the local neurons generate calcium dependent active responses that *do not propagate* [27]. With extracellular recordings the projection neuron spikes are easily observed, *but the local neuron activity is not*. Both cell types receive inputs from the olfactory receptor neurons located on the insect antenna. Odor stimulation can lead to: a direct excitation of the projection neurons; an indirect inhibition of the projection neurons via interposed local neurons; a combination of both in succession (see Fig. 17). Given the network structure [24] and the recording constraints—neither the receptors, nor the local neurons are observed—, we can expect some co-firing of more than two observed projections neurons (because three or more of these neurons receive a common excitatory input from the receptors), as well as some co-pauses (because several neuron receive a common inhibitory input).
In order to follow the analysis presented next, the reader needs to know some details of the experiment. The dataset contains both recordings of spontaneous activity, and of odor responses. Several odors were used : cis-3-hexen-1-ol (the main component of freshly cut grass odor), citral (the main component of lemon odor), mint (a mixture of several compounds), octanol (an alcohol part of most cut green leave odors) and vanilla (also a mixture of several compounds). For each odor there were typically 25 to 30 stimulations applied, one after the other, with a 30 seconds time interval between successive stimulations<fnref n="6" />. The total recording time for a given odor was typically 12.5 to 15 minutes. The dataset includes 4 acquisition-sets<fnref n="7" /> in the spontaneous regime (that is acquisition-sets during which no odor was presented), 6 acquisition-sets with cis-3-hexen-1-ol, and 1 acquisition-set for each of the 4 other odors. We do not take into account the data corresponding to the fourth cis-3-hexen-1-ol acquisition-set, because there were only 11 odor presentations among the 25 planned<fnref n="8" />. We thus have 13 spike trains. In each spike train, 7 neurons were observed. The actual sequence of acquisition-sets was: spontaneous regime 1; spontaneous regime 2; cis-3-hexen-1-ol 1; citral ; vanilla; octanol; mint; cis-3-hexen-1-ol 2; spontaneous regime 3; spontaneous regime 4; cis-3-hexen-1-ol 3; cis-3-hexen-1-ol 4; cis-3-hexen-1-ol 5; cis-3-hexen-1-ol 6.

## 6.2 The “usual” analysis of this data set

The most common analysis performed by neurophysiologists dealing with stimulus responses is the construction of a *peristimulus time histogram* [26, p. 232]: for each neuron, the spike times of each stimulus/odor presentation are shifted such that the stimulus onset time is always the same (say, time 0); the time axis is then discretized/binned; the total number of spikes (from all the stimulus presentations of a given acquisition-set) falling in each bin are counted; the counts are divided by the binwidth and by the number of presentations in order to get an *estimated time-dependent spike frequency*. Since we are not interested here in the shape of this time-dependent spike frequency *per se*, but by its variations across neurons and odors, we use a *variance stabilized* version of the *peristimulus time histogram*, stemming from an idea of John Tukey's *hanging rootogram* [33, Sec. 9.7, p. 349]. The difference with the classical peristimulus time histogram is that, at the last step, instead of dividing each bin count, *bc*, by the binwidth and by the number of presentations, we take $\sqrt{bc} + \sqrt{bc} + 1$. This stabilizes the variance at 1 with a minimal bias when counts are close to zero [19]. The Jupyter notebook giving all the details of this procedure applied to our data can be found on a dedicated Gitlab page<fnref n="9" />.

Fig. 17 shows the stabilized peristimulus time histograms for each of the seven neurons

<footnote n="6">⁶Receptor adaptation is therefore absent and the successive responses to a given odor can be reasonably modeled as IID draws from an unknown distribution.</footnote>

<footnote n="7">⁷By “acquisition-set”, we mean a sequence of 25 to 30 successive presentations of an odor.</footnote>

<footnote n="8">⁸These experiments were done 26 years ago and the acquisition had to be stopped in order to transfer data out a filled-up hard-disk. With modern equipment this would not happen.</footnote>

<footnote n="9">⁹https://gitlab.com/c_pouzat/aspa/-/blob/master/python/locust_data_analysis/locust20010214_sta.ipynb</footnote>
and each of the spontaneous regime and odor acquisition-sets. The first column shows the
spontaneous regime. Since no odor is presented, we expect *flat* curves with tiny statistical
fluctuations of order 1. The curves for successive acquisitions should moreover overlap.
This is mostly the case, *except for neuron 6*, where a shift upward is observed for the last
two acquisitions (gray and brown curves); they come after each different odor has been
presented for at least one cycle of 25 to 30 stimulations. The second columns displays
the responses of the 5 complete acquisition-sets with cis-3-hexen-1-ol. We see that all
neurons, *except neuron 6*, exhibit essentially stable responses. Neuron 6 exhibits again
a “problematic” upward shift starting with cis-3-hexen-1-ol 2 (in blue), which comes in
real time after the mint acquisition and before spontaneous 3 and 4 (gray and brown
curves on the graph to the left). The shift increases further before stabilization with
cis-3-hexen-1-ol 3 (gray). Concomitant with this baseline shift, a strong monophasic
excitatory response develops, while the excitatory response in the first cis-3-hexen-1-ol
acquisition and in all the other odors is very weak (right column). This is most likely
due to a contamination (at the level of spike sorting) by another neuron, as suggested by
a deteriorated refractory period for this neuron starting with the second cis-3-hexen-1-ol
acquisition (see the detailed account of our sorting procedure referred to at the beginning
of this section). The acquisition-sets starting from cis-3-hexen-1-ol 2 should therefore
either be ignored for this neuron, or great care should be taken when interpreting the late
responses. The third column displays the responses to the *first* cis-3-hexen-1-ol (black),
to citral (blue), to vanilla (gray), to octanol (brown) and mint (orange). We see that
neuron 1 (top row) exhibits similar responses (excitation followed by inhibition) to all
the odors, except to octanol (brown), where the inhibition is absent. Neuron 2 exhibits
also mostly a biphasic response (inhibition followed by excitation), but the excitation
part in the response to octanol (brown) is tiny. Neuron 3 shows at most a very weak
inhibition. Neuron 4 is the most constant across odors with a monophasic inhibitory
response. Neuron 5 exhibits a strong monophasic inhibition for all the odors, except
for citral (blue) where a biphasic (inhibition-excitation) response is observed. Neuron
6, before it becomes “problematic”, exhibits a consistent weak excitatory response, like
neuron 7.

In summary, there is a clear difference between the spontaneous regime and odor
responses. Differences between odor responses are absent or weak (neuron 1 for octanol,
neuron 5 for citral). Neuron 6 exhibits a “contamination” starting with the second cis-
3-hexen-1-ol acquisition-set and should be dealt with caution, or not considered at all,
when analyzing the end of the experiment.

## 6.3 Our proposed Curto-Itskov filtration

When we apply the Curto-Itskov filtration with $\delta t = 10$ milliseconds, followed by persistent homology, we get Fig. 18.

We now proceed with comparisons between the observed persistence diagrams and the
ones obtained by simulations of well defined simple settings.
<figure><img src="image_18.png" /><figcaption>Figure 17: Variance stabilized peristimulus time histograms of the seven neurons (one per row, neuron 1 at the top and neuron 7 at the bottom). First column: the 4 spontaneous regime acquisitions, the first in black, the second in blue, the third in gray, the fourth in brown. Second column: the 5 complete cis-3-hexen-1-ol acquisitions with the colors black, blue, gray, brown and orange. Third column: the first cis-3-hexen-1-ol (black), citral (blue), vanilla (gray), octanol (brown); mint (orange). Variance is stabilized at 1. The stimulus onset is at 0. The ordinate range is always [0, 30].</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 18: Persistence diagrams of spontaneous activity and odour stimulations, of the locust olfactory pathway. C3H stands for cis-3-hexen-1-ol.</figcaption></figure>
### 6.3.1 Comparison to the independent case

First, we compare these diagrams to the independent case with seven neurons, to assess if in real data there exists coactivity between some sets of neurons. To do that, for each odour (or spontaneous activity), 1000 spike trains are simulated with the same firing rates as the real spike train corresponding to this odour, and with independent neurons. We apply the Curto-Itskov filtration and we evaluate if the points on the diagram from real data depart from the points in the corresponding simulated independent case. This is done in dimension 0 and in dimension 1.

Figure 19 shows the results for dimension 0. For each odour, we plot in black the 6 points of the real diagram (we omit the point at $y = \infty$) and in red the $6 \times 1000$ points of the 1000 simulated diagrams (again, the points at infinity are not represented). We apply a K-means clustering with 6 clusters and draw the boundary of the convex hull of each cluster. We can see that the black points are often at the left or at the bottom of the points of the simulated diagrams. This shows that some pairs of neurons are coactive in the real data, as we observed in the simulations of Section 5.

Figure 20 shows the results for dimension 1. The points of the real diagram are plotted in black and the points of the simulated diagrams in blue. This time we draw the convex hull of all the blue points. We can see that some black points are at the left of the blue points, meaning that there are holes of dimension 1 appearing early in the filtration. Since holes of dimension 1 are formed by edges, it shows again that some pairs are coactive in the real data.

The diagrams of the real data thus depart from the diagrams in the independent case. It shows that neuron coactivity is present in the real spike trains.

### 6.3.2 Clustering

We want to see if the Curto-Itskov filtration can discriminate the different odours, or the odour presentation regime from the spontaneous activity. Ideally we would like to partition the spike trains in 6 clusters: one for spontaneous, one for cis-3-hexen-1-ol, and one for each of the 4 other odours.

Visually the diagrams are hard to discriminate, even if we zoom on some parts. An automatic clustering is carried out using Betti curves. In this setting, Betti curves in dimension 0 do not discriminate at all (not shown). We plot Betti curves in dimension 1, in Figure 21. Visually 3 clusters appear. We still tried and carried out an automatic clustering with 6 clusters, but it does not give a satisfying result<fnref n="10" />: it is not surprising since even visually we do not distinguish 6 clusters. We apply the K-means algorithm to the Betti curves in dimension 1, with 3 clusters. We obtain the following results.

* Cluster 1: Spontaneous 1, Spontaneous 2.

<footnote n="10">¹⁰The results are presented in the notebook.</footnote>
<figure><img src="image_20.png" /><figcaption>Figure 19: Comparison in dimension 0 between the real diagrams (in black) and 1000 simulated diagrams (in red) with independent neurons.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 20: Comparison in dimension 1 between the real diagrams (in black) and 1000 simulated diagrams (in blue) with independent neurons.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 21: Betti Curves in dimension 1 of the locust dataset.</figcaption></figure>

• Cluster 2: C3H₁, Citral, Vanilla, Octanol, Mint.

• Cluster 3: C3H₂, C3H₃, C3H₅, C3H₆, Spontaneous 3, Spontaneous 4.

This classification is consistent with the observations obtained from the “usual” psth analysis (Sec. 6.2). Cluster 1 and 2 are successful in discriminating spontaneous activity from odor responses. Cluster 3 is generated mostly by the problematic response of neuron 6. Indeed, two spontaneous sequences are misclassified and it turns out that we can detect these two problematic sequences in the first column of Figure 17 as corresponding to neuron 6. Also, the second and third columns of Figure 17 display differences between C3H and the other odours only for neuron 6. This explains the presence of this additional cluster which contains most of the C3H responses. Considering the experimental problems, this clustering is as satisfying as it can be. Moreover, our method provides a collection of curves for each dimension, from which it can be easy to distinguish visually if a clustering is appropriate and, if so, to pick up the proper number of clusters.

**Remark 6.1.** The Betti curves of dimension 0 and 2 (not shown) give less satisfying results. We also tried and used other vectorizations of persistence diagrams for our clustering purpose, like persistence images [1], landscapes [6], silhouettes [11], Atol [34]. They either give the same results as the Betti curve in dimension 1, or worse results. We also tried to cluster the diagrams with very simple summaries such that the mean of birth times, death times, lifetimes of the points (that is, of the corresponding homological features), the number of points, the number of points at infinity, but it never yields a satisfying clustering.

# 7 Conclusion

We have shown that the Curto-Itskov filtration provides a relevant way to analyze spike trains with persistent homology. Our theoretical stability result (Sec. 4) gives moreover a robustness guarantee (something quite rare for spike train analysis methods). We
checked by simulations (Sec. 5) that the Curto-Itskov filtration behaves as expected, in particular, we are confident that cofirings are detected if they are present. Its application to actual data (Sec. 6) from the first olfactory relay of an insect, leads to results fully compatible with the “usual” analysis methods. This filtration complements therefore other spike train analysis techniques. There is clearly significant work ahead: the determination of the discretizing time window used to map the spike train to frequency vectors (Sec. 5.1) will be network dependent (the projection neurons of the locust antennal lobe, exhibit particularly long minimal inter-spike intervals); the method used for assessing the significance of the observed point positions on the persistence diagrams (Sec. 6.3.1) will need refinement; the simulation codes, written in Python for now, will have to be re-written in a compiled language in order to make the approach quick enough for satisfying our experimentalist colleagues. But the codes implementing the whole procedure are already freely available, our colleagues can therefore give it a try on their data with a minimal effort.

## Acknowledgments

The authors would like to thank Bertrand Michel for carefully answering our questions. This work was supported by the Agence Nationale de la Recherche (ANR) project ANR-22-CE45-0027 SIMBADNESTICOST.

## References

[1] Henry Adams et al. “Persistence images: a stable vector representation of persistent homology”. In: *J. Mach. Learn. Res.* 18.1 (Jan. 2017), pp. 218–252. ISSN: 1532-4435. DOI: 10.5555/3122009.3122017.

[2] G. Azumaya. “Corrections and supplementaries to my paper concerning Krull-Remak- Schmidt’s theorem”. In: *Nagoya Math. J.* 1 (1950), pp. 117–124.

[3] Jean-Daniel Boissonnat, Frédéric Chazal, and Mariette Yvinec. *Geometric and Topological Inference*. Cambridge, England, UK: Cambridge University Press, Sept. 2018. ISBN: 978-1-10829780-6. DOI: 10.1017/9781108297806.

[4] D. R. Brillinger. “Maximum likelihood analysis of spike trains of interacting nerve cells”. In: *Biological Cybernetics* 59.3 (Aug. 1988), pp. 189–200. ISSN: 1432-0770. DOI: 10.1007/bf00318010. URL: http://dx.doi.org/10.1007/bf00318010.

[5] David R. Brillinger, Hugh L. Bryant, and José P. Segundo. “Identification of synaptic interactions”. In: *Biol. Cybern.* 22.4 (Dec. 1976), pp. 213–228. ISSN: 1432-0770. DOI: 10.1007/BF00365087.

[6] Peter Bubenik. “Statistical topological data analysis using persistence landscapes”. In: *J. Mach. Learn. Res.* 16.1 (Jan. 2015), pp. 77–102. ISSN: 1532-4435. DOI: 10.5555/2789272.2789275.
[7] György Buzsáki. "Neural Syntax: Cell Assemblies, Synapsembles, and Readers". In: *Neuron* 68.3 (Nov. 2010), pp. 362–385. ISSN: 0896-6273. DOI: 10.1016/j.neuron.2010.09.023. URL: http://dx.doi.org/10.1016/j.neuron.2010.09.023.

[8] Eduarda Gervini Zampieri Centeno et al. "A hands-on tutorial on network and topological neuroscience". In: *Brain Structure and Function* 227.3 (Feb. 2022), pp. 741–762. ISSN: 1863-2661. DOI: 10.1007/s00429-021-02435-0. URL: http://dx.doi.org/10.1007/s00429-021-02435-0.

[9] Frédéric Chazal, Vin de Silva, and Steve Oudot. "Persistence stability for geometric complexes". In: *Geom. Dedicata* 173.1 (Dec. 2014), pp. 193–214. ISSN: 1572-9168. DOI: 10.1007/s10711-013-9937-z.

[10] Frédéric Chazal et al. "Proximity of Persistence Modules and their Diagrams". In: *INRIA* (2008), p. 29. URL: https://inria.hal.science/inria-00292566v4.

[11] Frédéric Chazal et al. "Stochastic Convergence of Persistence Landscapes and Silhouettes". In: *SOCG'14: Proceedings of the thirtieth annual symposium on Computational geometry*. New York, NY, USA: Association for Computing Machinery, June 2014, pp. 474–483. ISBN: 978-1-45032594-3. DOI: 10.1145/2582112.2582128.

[12] Frédéric Chazal et al. *The Structure and Stability of Persistence Modules*. Cham, Switzerland: Springer International Publishing, 2016. ISBN: 978-3-319-42545-0. URL: https://link.springer.com/book/10.1007/978-3-319-42545-0.

[13] E. S. Chornoboy, L. P. Schramm, and A. F. Karr. "Maximum likelihood identification of neural point process systems". In: *Biological Cybernetics* 59.4–5 (Sept. 1988), pp. 265–275. ISSN: 1432-0770. DOI: 10.1007/bf00332915. URL: http://dx.doi.org/10.1007/bf00332915.

[14] David Cohen-Steiner, Herbert Edelsbrunner, and John Harer. "Stability of Persistence Diagrams". In: *Discrete Comput. Geom.* 37.1 (Jan. 2007), pp. 103–120. ISSN: 1432-0444. DOI: 10.1007/s00454-006-1276-5.

[15] William Crawley-Boevey. "Decomposition of pointwise finite-dimensional persistence modules". In: *J. Algebra Appl.* 14.05 (Aug. 2014), p. 1550066. ISSN: 0219-4988. DOI: 10.1142/S0219498815500668.

[16] Carina Curto and Vladimir Itskov. "Cell Groups Reveal Structure of Stimulus Space". In: *PLoS Comput. Biol.* 4.10 (Oct. 2008), e1000205. ISSN: 1553-7358. DOI: 10.1371/journal.pcbi.1000205.

[17] Carina Curto and Nicole Sanderson. "Topological Neuroscience: Linking Circuits to Function". In: *Annual Review of Neuroscience* 48.1 (July 2025), pp. 491–518. ISSN: 1545-4126. DOI: 10.1146/annurev-neuro-112723-034315. URL: http://dx.doi.org/10.1146/annurev-neuro-112723-034315.

[18] Steven P. Ellis and Arno Klein. "Describing high-order statistical dependence using "concurrence topology," with application to functional MRI brain data". In: *Homology, Homotopy Appl.* 16.1 (June 2014), pp. 245–264. ISSN: 1532-0081. DOI: 10.4310/HHA.2014.v16.n1.a14.
[19] Murray F. Freeman and John W. Tukey. "Transformations Related to the Angular and the Square Root". In: *The Annals of Mathematical Statistics* 21.4 (Dec. 1950), pp. 607-611. ISSN: 0003-4851. DOI: 10.1214/aoms/1177729756. URL: http://dx.doi.org/10.1214/aoms/1177729756.

[20] Antonio Galves, Eva Löcherbach, and Christophe Pouzat. *Probabilistic Spiking Neuronal Nets: Neuromathematics for the Computer Era*. Springer International Publishing, 2024. ISBN: 9783031684098. DOI: 10.1007/978-3-031-68409-8. URL: http://dx.doi.org/10.1007/978-3-031-68409-8.

[21] R. Ghrist. "Barcodes: The persistent topology of data". In: *Bulletin of the American Mathematical Society* (2007). URL: https://www semanticscholar.org/paper/Barcodes%3A-The-persistent-topology-of-data-Ghrist/cad0b ff789bb0e626245b21a47791ef76a9cd5c2.

[22] Allen Hatcher. *Algebraic topology*. eng. Cambridge: Cambridge University Press, 2001. ISBN: 978-0-521-79160-1.

[23] Stephen L Keeley et al. "Modeling statistical dependencies in multi-region spike train data". In: *Current Opinion in Neurobiology* 65 (2020). Whole-brain interactions between neural circuits, pp. 194-202. ISSN: 0959-4388. DOI: https://doi.org/10.1016/j.conb.2020.11.005. URL: https://www.sciencedirect.com/science/article/pii/S0959438820301720.

[24] Gilles Laurent. "Dynamical representation of odors by oscillating and evolving neural assemblies". In: *Trends in Neurosciences* 19.11 (Nov. 1996), pp. 489-496. ISSN: 0166-2236. DOI: 10.1016/s0166-2236(96)10054-0. URL: http://dx.doi.org/10.1016/s0166-2236(96)10054-0.

[25] Arthur Leblois and Christophe Pouzat. "Neurobiology of Motor Control". In: ed. by Scott L. Hooper and Ansgar Büschges. Wiley, 2017. Chap. Multi-Unit Recording. ISBN: 9781118873397. DOI: 10.1002/9781118873397.ch3. URL: http://dx.doi.org/10.1002/9781118873397.ch3.

[26] Liqun Luo. *Principles of Neurobiology*. Garland Science, Sept. 2020. ISBN: 9781003053972. DOI: 10.1201/9781003053972. URL: http://dx.doi.org/10.1201/9781003053972.

[27] Katrina MacLeod and Gilles Laurent. "Distinct Mechanisms for Synchronization and Temporal Patterning of Odor-Encoding Neural Assemblies". In: *Science* 274.5289 (Nov. 1996), pp. 976-979. ISSN: 1095-9203. DOI: 10.1126/science.274.5289.976. URL: http://dx.doi.org/10.1126/science.274.5289.976.

[28] George P. Moore et al. "Statistical Signs of Synaptic Interaction in Neurons". In: *Biophysical Journal* 10.9 (Sept. 1970), pp. 876-900. ISSN: 0006-3495. DOI: 10.1016/s0006-3495(70)86341-x. URL: http://dx.doi.org/10.1016/s0006-3495(70)86341-x.

[29] James R. Munkres. *Elements of algebraic topology*. eng. Reading (Mass.): Addison-Wesley, 1984. ISBN: 0-201-04586-9.
[30] Donald H. Perkel, George L. Gerstein, and George P. Moore. "Neuronal Spike Trains and Stochastic Point Processes: II. Simultaneous Spike Trains". In: *Bio-physical Journal* 7.4 (1967), pp. 419-440. ISSN: 0006-3495. DOI: https://doi.org/10.1016/S0006-3495(67)86597-4. URL: https://www.sciencedirect.com/science/article/pii/S0006349567865974.

[31] Christophe Pouzat. "Brain-Computer Interfaces 1: Methods and Perspectives". In: ed. by Maureen Clerc, Laurent Bougrain, and Fabien Lotte. iSTE/Wiley, 2016. Chap. Analysis of Extracellular Recordings.

[32] Christophe Pouzat, Ofer Mazor, and Gilles Laurent. *Extracellular recordings from the locust, Schistocerca americana, olfactory pathway*. Zenodo, July 2015. DOI: 10.5281/zenodo.21589. URL: https://doi.org/10.5281/zenodo.21589.

[33] John A. Rice. *Mathematical Statistics and Data Analysis, Third Edition*. Duxbury, 2007.

[34] Martin Royer et al. "ATOL: Measure Vectorization for Automatic Topologically-Oriented Learning". In: *International Conference on Artificial Intelligence and Statistics*. PMLR, Mar. 2021, pp. 1000-1008. URL: http://proceedings.mlr.press/v130/royer21a.html.

[35] Cary Webb. "Decomposition of graded modules". In: *Proc. Amer. Math. Soc.* 94.4 (1985), pp. 565-571. ISSN: 0002-9939. DOI: 10.1090/S0002-9939-1985-0792261-6.