Geometry of spiking patterns in early visual cortex: a topological data analytic approach

In the brain, spiking patterns live in a high-dimensional space of neurons and time. Thus, determining the intrinsic structure of this space presents a theoretical and experimental challenge. To address this challenge, we introduce a new framework for applying topological data analysis (TDA) to spike train data and use it to determine the geometry of spiking patterns in the visual cortex. Key to our approach is a parametrized family of distances based on the timing of spikes that quantifies the dissimilarity between neuronal responses. We applied TDA to visually driven single-unit and multiple single-unit spiking activity in macaque V1 and V2. TDA across timescales reveals a common geometry for spiking patterns in V1 and V2 which, among simple models, is most similar to that of a low-dimensional space endowed with Euclidean or hyperbolic geometry with modest curvature. Remarkably, the inferred geometry depends on timescale and is clearest for the timescales that are important for encoding contrast, orientation and spatial correlations.


Introduction
The broad goal of understanding the nature and function of neural activity entails the intermediate step of characterizing the space in which this activity lives. This is a challenging problem, even for small populations of neurons because of the dimensional explosion and data limitations: it is simply not feasible to obtain sufficient experimental data to exhaustively sample the space of all possible spiking patterns [1]. For this reason, indirect methods are required. A promising strategy, applicable in a wide variety of contexts, is the persistent homology approach to topological data analysis (TDA) [2][3][4][5][6]. With TDA, an abstract notion of distance is used to compare observations, such as samples of neural activity. Based on these pairwise distances, a sequence of network graphs can be constructed, with successive graphs built by connecting the nodes at increasingly greater distances. These graphs are then associated with higher dimensional topological objects called simplicial complexes, which are characterized via their Betti numbers, describing (intuitively speaking) the number of disconnected components, the number of holes, the number of 'bubbles', etc. The way that the Betti numbers change as a function of the distance criterion between the samples then provides a topological characterization of the space.
Here, we use this approach to analyse the activity of clusters of cortical neurons, with a specific focus on patterns of spiking activity in time. This requires confronting another issue, but one that the persistent homology approach is also well suited to handle: even for a single neuron, activity is not a scalar quantity, but rather, a sample of a point process. Thus, it is possible to use the TDA machinery in a novel way, by applying it to a measure of distances between neural responses that is sensitive to temporal patterns within and across neurons: the family of spike distances introduced by Victor and Purpura [7]. While these distances do not correspond to Euclidean distances [8], they do satisfy the triangle inequality, which suffices for the TDA machinery.
While our approach shares the overall goal of a geometric characterization of neural activity with many other studies, there are some important differences in the questions we ask. Below we highlight the two main distinctions.
One central question is, what is the dimensionality of the response space? In previous studies [9,10], dimensionality is defined as the dimension of the space that contains response trajectories, where a response trajectory is the evolution of average firing rates over time. Thus, these previous approaches would not distinguish between spiking activity that was well described by a Poisson process and spiking activity whose patterning required additional dimensions to describe. Our construction of the response space is quite different: it is sensitive to the patterns of spikes on individual trials. We find that there are substantial differences between the geometry of the recorded neural activity and their Poisson-like equivalents. Moreover, we find that the influence of spike patterning on geometry is maximal at the timescales that are important for neural coding of visual attributes such as contrast, orientation and spatial patterns [11].
Many recent studies have shown that the geometry of the responses 'untangles' the key parameters of the stimulus, an important and novel finding for complex, naturalistic domains [12][13][14]. This is not a focus of the current study, as an untangling is to be anticipated for the stimulus set studied here: a family of artificial visual textures parametrized by local image statistics [15]. V1 and V2 neurons are selectively tuned to these statistics (see later), so decoding of the population will necessarily identify distinct directions (i.e. distinct population vectors) that capture stimulus parameters.
Rather, our focus is on whether the intrinsic geometry of neural activity patterns depends on the stimulus parameters being encoded. To this end, our visual stimuli contained a range of types of statistical structures, including textures that varied in first-, second-, third-and fourth-order local spatial correlations (see §4) known to be present in natural images [16] and to be salient for human observers [17]. First-and second-order statistics control luminance and orientation content, respectively, and thus drive neurons in both V1 and V2, while third-and fourth-order statistics control aspects of local form that are primarily detected by V2 neurons [18]. Here, we separately analyse responses to many examples of textures that explore each of these statistics. Thus, we can determine whether the geometry of neural activity patterns depends on the kind of visual structure that is encoded, or perhaps on the cortical area (V1 versus V2), or alternatively has a more universal behaviour-the latter suggesting that the spike patterning dynamics are endogenously driven.
Applying TDA to spike distances carries with it a theoretical departure from most previous applications of the persistent homology approach to neural data. Provided that mean responses depend continuously on stimulus parameters, Euclidean distances between mean firing rates or distances derived from correlations [3,19,20], which are effectively dot products in a vector space, will always lead to a manifold (although not necessarily one of low dimension). The reason is that these distances will map neural activity to a smooth and connected subset of a vector space. Spike metrics have a discrete component and consequently are fundamentally not Riemannian, so response geometry will typically not constitute a manifold-but, as we show, manifold structure is not needed to interpret the results of TDA.
Finally, we find that a novel way of passing from distance measures to graphs-linking nodes according to greatest dissimilarity rather than greatest similarity (using a 'decreasing filtration' rather than the standard 'increasing filtration')sharpens the distinction between the experimental data and candidate geometrical models.

Results
The starting point of our analysis is 28 datasets of single-unit spiking activity recorded from areas V1 and V2 of the visual cortex of anaesthetized macaques during stimulation with visual patterns drawn from a high-dimensional space of visual textures [15]. Each dataset is a recording at a separate recording site or in a different macaque, and at each site, the stimulus space is sampled along its 10 axes, corresponding to 10 independent and visually salient kinds of spatial correlation (figure 1). We sample each axis at four points: two levels of positive correlation and two levels of negative correlation (see §4 for details). Sixty-four perceptually similar examples of each of these 40 texture types (10 axes, four sample points per axis) were presented in a random order, with each stimulus lasting 320 ms. This sequence was repeated four times at each recording site. We refer to the neuronal responses to the 64 perceptually similar stimuli obtained in this way in one repeat as a 'collection'; thus, there were 160 collections (40 texture types, four repeats) in each dataset. For our analysis, we selected the 80 collections with the smallest number of empty responses from each dataset. Depending on the dataset, the responses consisted of the spike trains of one or more isolated neurons. If more than four neurons were simultaneously recorded, we considered only the four neurons with the highest firing rates. We then applied TDA to characterize the topology of the firing patterns in each of the collections and summarize the resulting characterizations to identify their consistent features.
The TDA method we used follows very closely the framework proposed by [3], but we started from a different way of endowing a set of neuronal responses with a notion of distance. Specifically, we used the multi-neuron Victor-Purpura distance [11,21] to define the distances between all pairs of neuronal responses. This choice to employ a distance between spike trains rather than the correlations (based for example on firing rate) between neuronal responses to quantify (dis)similarity is motivated by several considerations. The main consideration is that we wanted to focus on spiking patterns in neuronal data that are relatively sparse (typically only a few spikes in the 320 ms response period). Information about spiking patterns would be lost in averages over multiple stimulus presentations, or after a smoothing procedure to estimate time-varying firing rates from single trials. Secondarily, the use of the Victor-Purpura distance allows us to build on previous studies of the visual cortex that employ other data analytical techniques. We note that the generality of the TDA approach extends seamlessly to metrics such as the Victor-Purpura distance, as it does not assume a vector space embedding.
The Victor-Purpura distance between spike trains (see §4) is defined as the minimum cost of transforming one spike train into the other via addition or deletion of spikes, shift of spike times or change in the neuron of origin of the spikes. Each editing move is associated with a cost, which depends on a parameter q, controlling the relevant timescale for shifts of spike times or on a parameter k, controlling the sensitivity to the neuron of origin of each spike. The family of distances defined in this way ranges from measures of dissimilarity that ignore spike timing altogether (q = 0), measuring distance based only on spike count, to measures that consider spike timing meaningful at an arbitrarily high level of precision (q → ∞). Similarly, k = 0 corresponds to ignoring the neuron of origin, while k = 1 corresponds to considering a change in the neuron of origin equivalent to the insertion of a spike. We carried out TDA for a grid of values of these two parameters: q = 1, 2, 5, 10, 20, 50, 100, 200 (s −1 ) and k = 0, 1. The resulting grid covered the range found relevant to neural coding in previous studies [11,21].
For each choice of parameters of the Victor-Purpura distance, we computed topological summaries via Betti curves [3] (figure 2) and characterized relevant features of the topological space by computing mean Betti curves averaged over the 80 collections drawn from of each dataset. The procedure for construction of the Betti curves is detailed in §4 and summarized here. We first create a sequence of graphs from each collection; in each graph, the nodes are the individual responses in the collection, and the edges are determined by the dissimilarities (figure 2a). As in the standard implementation of TDA, each graph is built by connecting nodes whose dissimilarities are less than a given ceiling. Thus, the graph sequence begins with a ceiling of 0 (all nodes disconnected) and terminates at a ceiling equal to the maximum distance (all nodes connected), a graph sequence known as the increasing filtration, and used in the clique topology method of [3]. Here, we also create a graph sequence by adding in the edges in decreasing order (the decreasing filtration), like in [22]. In each filtration, as the edges of the graph fill in, its topology can be characterized by the number of enclosed 'holes' of different dimensions: the number of one-dimensional tunnels, two-dimensional voids and three-dimensional 'cavities', known as Betti numbers and denoted by β 1 , β 2 and β 3 , respectively. These are computed at each step of the filtration, and thus, the Betti numbers are functions of the edge density ρ, the fraction of potential edges that have been filled in: β 1 (ρ), β 2 (ρ) and β 3 (ρ). These quantities, which we computed over the range from ρ = 0 to 0.6 and then averaged over the 80 collections in a dataset, are shown in figure 2b.
Note that while the Betti curves are sensitive to the geometry and topology of the metric space of data points [3,5,6], they are invariant under monotonic transformations of the distance [3]. This means that they are sensitive only to the relative ordering of all distances.
Below, we will focus on the mean Betti curves from the 80 collections with the highest firing rates at each recording site, but first we investigated whether the Betti curves had a substantial dependence on the two critical experimental variables that distinguished the collections: whether the recordings were in V1 versus V2, and whether the spatial structure of the stimulus was low order (first and second order) versus high order (third or fourth order) (see figure 1, §4, and [17] for details and definitions). The loworder correlation structure is extracted by linear receptive fields and is already present in V1; third-and fourth-order correlation structure is primarily extracted in V2 [18].
To compare the topological characterizations of the collections of responses, we summarized the information in the Betti curves by two main features: the integral of the curves [3], referred to as integrated Betti values, and the position of  figure S3) do not depend on whether the visually salient structure of the textures is driven by first-and second-order spatial correlations, which lead to differences in luminance and spatial frequency that is already manifest in the cortical input, or by third-and fourth-order spatial correlations, which require intracortical calculations to extract [18].
Note that this analysis does not seek to determine whether the responses to one kind of texture differ from responses to another (whether they lie in different parts of a combined response space, e.g. by virtue of different firing rate profiles), but rather to analyse the intrinsic geometry of each collection via methods that focus on the firing pattern. Our results up to this point suggest that this geometry is largely independent of the recording area, or the stimuli that drive the responses. Given the similarity of the results of TDA across recording areas and texture types, we pool these results in the following further analysis and now focus on the mean Betti curves over the 80 collections within each of the datasets.
To assess which aspects of the spike train data might be responsible for the shape of these mean Betti curves (e.g. figure 2b), we synthesized surrogate spike train data of four different types by perturbing specific aspects of the original experimental data. The description of the surrogate data is detailed in §4 and summarized here.
1. Uniform resampling of spike times (U). The spike times of each response are resampled uniformly in the response interval (0-320 ms). 2. Exchange resampling of spike times between collections (EB).
The spike times of each response are resampled uniformly, without replacement, in the set of spike times of the 80 collections of the dataset. and k (sensitivity to neuron of origin) of the Victor-Purpura distance. The two top (respectively, bottom) rows correspond to the increasing (respectively decreasing) filtration method. The odd (respectively even) rows correspond to k = 0 (respectively k = 1). In the increasing filtration plots, the Betti curves for β 3 are close to zero.
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220677 3. Exchange resampling of spike times within collections (EW). As given in 2, but with the resampling restricted to the set of spike times of the collection to which the response belongs. 4. Poisson-generated spike trains (P). Each response is replaced with a sample of a Poisson process with the same firing rate.
In all cases, the spike times are resampled independently for each neuron that participates in a response. For each experimental dataset, we generated 20 surrogate spike datasets for each of the aforementioned four types, and computed Betti curves for these surrogates just as the curves were computed for the experimental data. An example of this analysis is shown in figure 5. As is typical across the datasets, the Betti curves associated with the experimental data were often quite distinct from those of the surrogate sets. These differences were most prominent in the mid-range values of q (5 to 50 s −1 ) and also more prominent for the decreasing filtration than for the increasing royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220677 filtration. In addition, the Poisson surrogates usually show a greater divergence from the real data than the other surrogates. One possible reason for this is that, unlike the other surrogates, the Poisson surrogates do not preserve the number of spikes in each response.
To determine the consistency of these features across the 28 datasets, we used the integrated Betti values and the centres of mass of the Betti curves. Thus, for every fixed Betti number (β 1 , β 2 or β 3 ), we considered the difference between the integrated Betti value of each surrogate (mean value over the 20 computations) and the integrated Betti value of the experimental data, and averaged them over the 80 collections from each dataset (see figure 6 for β 1 and electronic supplementary material, figures S4 and S5 for β 2 and increasing, k = 0  figure 5, the behaviour of the experimental data departs from the behaviour of the surrogates, especially for the mid-range values of q (5 to 50 s −1 ). This holds both when the neuron of origin is ignored (k = 0) and when it is considered relevant (k = 1), and it is seen for both increasing and decreasing filtrations. All surrogates behave in a way that is distinct from the experimental data. As shown in figure 6, three of the surrogates (U, EB, EW) behave similarly to each other, while the Poisson surrogate (P) deviates from the data more extensively. We confirmed these observations using the centres of mass of the Betti curves in place of the integrated Betti values (see statistical tests below). Figure 7 and electronic supplementary material, figures S6 and S7 assess the statistical significance of these observations for β 1 -β 3 via the two-sample Kolmogorov-Smirnov (KS) test applied to the distributions of integrated Betti values. Even for the EW surrogate, which perturbs the recorded spike trains the least, the KS test statistic takes on large values, especially for mid-range values of q, rejecting (at the 0.05 significance level) the null hypothesis that the EW surrogate and the experimental data samples come from the same distribution. More generally, the analysis shows that the surrogates for U, EB and EW diverge in a similar, timescale-dependent manner from the experimental data, while the Poisson-generated data are always significantly different from both the experimental data and the other surrogates. These observations are confirmed by the KS test applied in the same way to the distributions of centres of mass of the Betti curves (electronic supplementary material, figures S8-S10).
In sum, figures 3-7 and electronic supplementary material, figures S2-S10 demonstrate that the neural activity elicited by images with a range of statistical properties, and recorded in either V1 or V2, have a common topological structure that is distinct from that of surrogates with matching spike counts and spike timing distributions. These differences are present over a broad range of temporal scales and are most prominent at a temporal scale of 5-50 s −1 , i.e. 20-200 ms.
The analysis so far shows that the observed firing patterns depart from that of surrogates that randomize spike times but maintain mean firing rates in various ways. This departure is manifest in the Betti curves, but the Betti curves are only an indirect reflection of the structure of the response space. To gain insight into this topological structure, we compared the measured Betti curves with those obtained from randomly assigned distances or by sampling points in simple geometric spaces [3,4,6]. Specifically, we computed the Betti curves associated with: (i) a random symmetric 64 × 64 matrix with zeros on the diagonal [3], (ii) a sample of 64 random points in a unit (hyper)cube within a d-dimensional Euclidean space [3], and (iii) a sample of 64 random points in the hyperbolic ball model of the d-dimensional hyperbolic space [6]. For Euclidean and hyperbolic spaces, we considered dimensions d from 1 to 15. royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220677 For the hyperbolic models, we also varied the effective curvature. As given in [6], we implemented this by keeping intrinsic (Gaussian) curvature fixed at −1 but changing typical distances between the points (i.e. varying from 1 to 5 the maximum radius R max of the hyperbolic ball model, see §4). We therefore refer to models with low R max as hyperbolic spaces with modest curvature. For each random or geometric model, we computed the Betti curves of 300 sets of samples. We then compared these distributions with the Betti curves selected from our data via integrated Betti values and centres of mass.
We restricted the analysis to collections with at most one empty response, to have 64 data points (and to exactly match the geometric models), and we randomly chose one such collection for every dataset that had collections meeting this criterion. This yielded Betti curves from 18 different datasets.
To determine compatibility with the geometric models, we first tallied how many of the experimental integrated Betti values were within 3 standard deviations of the mean of the    figure S11). While none of these models were a good fit to the data, the Euclidean and low-curvature hyperbolic models   figure S11 shows that compatibility is reduced when considering centres of mass instead of integrated Betti values, in particular for k = 1. Compatibility is maximal for low dimensions (d = 3 to 5) in the range we considered, consistently across the geometric models. In all cases, Betti curves associated with the experimental data were found to be incompatible with the random symmetric matrix model. Interestingly, across the geometric models, compatibility of integrated Betti values between the experimental data and the models ( figure 8) is maximized in a consistent region of the parameter space of the Victor-Purpura distance. These regions correspond approximately to q ¼ 5 to 20 s À1 for k = 0, and q ¼ 1 to 10 s À1 for k = 1, while for large values of q, the compatibility is very low in all cases.

Geometry of spiking patterns
Understanding the properties of biological neural networks requires examining the structure of spike trains-the sequences and patterns of action potentials generated by groups of neurons in the cortex during the processing of information. We show here that characterizing the response space in which these patterns live through TDA provides insight into the intrinsic structure of the activity patterns generated in V1 and V2 in the macaque monkey during robust visual stimulation. The influence of timescale on this structure suggests that local networks in V1 and V2, and perhaps the interactions between these two cortical areas, play a major role in shaping the geometry of the spike response space, and in fact dominate the influence of the stimuli. While no simple model such as a circle, torus or sphere [4] recapitulates the space, across a range of timescales the spaces extracted by TDA are most similar to low-dimensional models of Euclidean geometry or hyperbolic geometry with a modest amount of negative curvature (figure 8 and electronic supplementary material, figure  S11). Moreover, these spaces appear to be largely independent of the visual stimuli used in this study (figure 4 and electronic supplementary material, figure S3) or the cortical region of the recordings (figure 3 and electronic supplementary material, figure S2). Note that while the spike pattern structures we identify are most consistent with Euclidean geometry or hyperbolic geometry with modest curvature, this geometry arises in a way that is fundamentally different from the manifolds and subspaces described earlier in §1. The latter are derived from a vector state space of instantaneous spike rates [23], or calcium activity across a population [24], and the time evolution of these signals can be interpreted as trajectories across these subspaces. Here, we begin with a family of metric state spaces, each defined by a timescale for comparing individual spike trains, and through TDA identify a family of spaces that each capture a unique spike pattern structure. These spike pattern structures exist in time as entities. Furthermore, vector space approaches assume that the actual timing of the spikes underlying the firing rates in the state space are largely irrelevant and can be effectively modelled as a Poisson point process, either homogeneous or inhomogeneous [4], with no temporal structure [25]. In our approach, we are specifically interested in the spike timing in individual and joint responses across groups of neurons. Implicit in our results is the demonstration that TDA provides a way to characterize neuronal response spaces avoiding the limitations of vector space methods. TDA begins with the pairwise distances between instances of population activity. While several previous applications of TDA to neural spiking data have assumed a vector space structure or used a vector space embedding to compute these distances [3,4,19,20], there is no fundamental reason to do so. Rather, the primary requirement is that the distance between two samples of neural activity corresponds to their dissimilarity. In [26], TDA methods combined with nonparametric (dis)similarity measures between spike trains were applied to synthetic spike train data for regimes classification in artificial neural networks. Here, we apply TDA in conjunction with a parametrized family of distances that have been shown to capture the meaningful dissimilarities between spike trains in several contexts [11,[27][28][29], and that are not Euclidean [8].

Spiking responses driven by visual textures
Recognizing that neural activity is the result of an interaction of the inputs to the population and its intrinsic network properties, we sought to identify characteristics of the spiking response space that might apply to a wide range of visual stimuli. To this end, we used stimuli consisting of mathematically defined visual textures [17,18] of known relevance to natural vision [30]. Critically, this stimulus set included some textures whose visually salient structure (first-and second-order spatial correlations) could be extracted by simple centre-surround operations in the retina, and other textures whose structure (third-and fourth-order spatial correlations) can be extracted only by cortical computations, primarily in V2 [18]. As we show (figures 3-4 and electronic supplementary material, figures S2 and S3), the characteristics we identify are independent of this distinction; that is, these characteristics are shared by V1, where only first-and second-order spatial correlations are robust drivers of activity, and V2, where many neurons are sensitive to higher-order correlations [18]. During the data collection phase, we matched the orientation and spatial scale of the texture stimuli to the receptive field properties of the well-isolated single units in the tetrode recordings. This allowed us to collect robust spiking activity even under propofol anaesthesia/ sufentanil analgesia [18]. These responses included strongly driven modulations of firing rate driven by many of the texture types, as well as trial-to-trial variability, thus facilitating wide exploration of the response space. We further promoted the exploration of the response space through the use of 64 distinct exemplars of each texture type. Finally, to ensure sensitivity to the details of spiking patterns, we based the analysis on individual spike trains rather than their averages and used a distance that was sensitive to the timing of individual spikes rather than just overall spike count [4] or firing rate envelope [31].

The importance of spiking pattern analysis in neuroscience
There is a considerable body of work spanning more than three decades, both experimental and theoretical, dedicated to elucidating the nature and function of spiking patterns in single neurons and across populations of neurons [32][33][34].
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220677 Several themes have emerged from this work. One is that the spiking patterns of individual neurons and groups of neurons can encode sensory information in the timing of their spikes [11,35]. Second, the degree of synchrony of spiking patterns and the generation of discrete sequences of spikes across neurons may play a role in transferring signals through synaptic pathways [36,37]. Third, spiking patterns can arise through phase-locking of spikes to slower oscillatory signals in neural networks, a mechanism that is thought to mediate communication across cortical areas [38]. Fourth, spiking patterns in individual neurons may arise from joint activity in neural assemblies [39]. Finally, plasticity and learning in neural networks are dependent on the precise timing of pre-and post-synaptic spikes across the synapses in the network [40,41]. Together, these themes emphasize the critical role spiking patterns play in brain function. However, the studies within this body of work use a wide range of approaches; of necessity, each provides only a partial picture of the spiking patterns of interest. Our study here is the first to apply TDA to the analysis of spiking patterns in the visual cortex. This study identifies low-dimensional structure in the spiking response space and characteristic timescales, reflecting the presence of spiking patterns in the visual cortex generated during robust visual stimulation. The patterns do not encode the visual stimuli; rather, like the low-dimensional spaces that organize during spontaneous spiking activity [42], they may reflect intrinsic states in visual cortical processing.

Timescales and neurons of origin in spike patterns
Our approach recognizes that the timing of individual spikes is potentially important, but a priori, it is agnostic as to what the relevant timescales are. This viewpoint is implemented by applying TDA to a family of distances controlled by a parameter q (in units of s −1 ) that sets the resolution with which spike timing influences the measure of dissimilarity. Applying TDA to the resulting family of distances for q = 1 to 200 s −1 allows us to focus on a range of timescales from 1 s down to 5 ms. The family of distances has a second parameter, k, which controls the importance of the neuron of origin in determining dissimilarity. For k = 0, the neuron of origin is irrelevant; for k = 1, changing the neuron of origin has unit cost. These extremes proved useful in analysing multi-neuronal coding of spatial phase [21], demonstrating maximal information near k = 1; here, the dependence of response space geometry on k was relatively small. Importantly, the distinction between experimental data and the random surrogates (U, EB, EW, P) depends systematically on the temporal resolution q of the distance between spike trains (figures 6 and 7 and electronic supplementary material, figures S4-S10). The distinction is greatest for intermediate values of q = 5-50 s −1 , corresponding to resolutions of 200-20 ms). For low values of q (less than 5), the distinction is generally lost, as the Victor-Purpura distance progressively disregards spike timing; for higher values (greater than 100), it is also generally weaker, indicating that at these timescales (with interspike intervals less than 10 ms), the systematic effect on geometry is less pronounced. Although the behaviour of Poisson surrogates statistically diverges from the experimental data for the decreasing filtration method and high values of q (figures 5 and 7 and electronic supplementary material, figures S6-S10), we observe (figures 6 and 7 and electronic supplementary material, S4-S10) that the distinction is clearly visible across the whole range of parameters and for both filtration methods. One can interpret the dependence between the timescale parameter q of the Victor-Purpura distance and the difference between experimental data and the random surrogates as the temporal scale of distinctive geometries that structure the spiking response space. It is notable that the timescales which most clearly reveal the geometry of the response space correspond to the timescales that are most informative for carrying visual information about contrast, spatial frequency, orientation and texture-both as identified by analysis methods based on the Victor-Purpura distance [11] and by unrelated approaches [43]. Because the timescales at which the ongoing activity has the most distinctive geometric structure are similar to the timescales that are most informative about visual features, it is reasonable to hypothesize that this match facilitates transmitting visual information.
One way to interpret the results with the random surrogate data is to consider the degrees of freedom of the spiking activity. The Betti curves of the surrogates typically have higher values than the experimental data for intermediate values of q (figure 6 and electronic supplementary material, figures S4 and S5). This may indicate that the response space for the real spiking activity is more constrained than the synthetic spiking activity. In this range (q = 5-50 s −1 ), the random manipulation of the spike data in different ways consistently opens up holes in the response space, leading to higher Betti values (figures 5 and 6 and electronic supplementary material, figures S4 and S5). This result suggests that the neural response space produced during visual stimulation has fewer discontinuities (is more structured) than the response spaces produced by the surrogates.

Comparison with simple geometric models
The comparison with random and geometric models highlights complementary aspects of the geometry of the space of sampled spike trains endowed with Victor-Purpura distances. Upon defining a notion of compatibility that accounts for the Betti curves of β 1 -β 3 and both increasing and decreasing filtration, we observed that even if none of the considered models consistently fits the data, Euclidean geometry and hyperbolic geometry with moderate curvature (R max = 1) are more compatible with the geometry of the spike train responses (figure 8 and electronic supplementary material, figure S11). Compatibility is concentrated in low dimensions of the geometric models (d = [3][4][5]. This observation is consistent with other geometrical descriptions of neural population activity [44]. In these descriptions, the instantaneous firing rates of a large number of neurons are found to occupy a low-dimensional manifold within a high-dimensional spiking response space. While in the study described here, we analyse recordings from at most four neurons, the distances we consider between spike trains are intrinsic and do not depend on the choice of an embedding, and as a consequence, potentially capture aspects of the neural activity manifold. However, recordings of larger numbers of neurons are needed to demonstrate that this is the case. Figure 8 also shows that the compatibility between the data and the geometric models assessed via royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220677 integrated Betti values systematically depends on the timescale q of the Victor-Purpura distance between spike trains and is maximized for mid-range values q = 5-20 s −1 when the neuron of origin of each spike is disregarded (k = 0), and for low-range values q = 1-5 s −1 when the origin of each spike is accounted for (k = 1). When assessed via centres of mass, however, compatibility is reduced for k = 0 and very low for k = 1 (electronic supplementary material, figure S11), even though a systematic dependence on q is still observable.

Considerations regarding the experimental procedure
The data analysed here were collected from monkeys under propofol anaesthesia/sufentanil analgesia and neuromuscular blockade. As a consequence, fixational eye movements could not have contributed to the fluctuations in spiking activity we observed in V1 and V2, but conversely, the natural dynamics of the visual input due to such eye movements is only partially approximated by the transient mode of stimulus presentation that we used. Another caveat is that the anaesthesia and opiate analgesia may have produced noise correlations in the neural activity that contributed to the topological structure we extracted through TDA. While it is impossible to rule out any impact of anaesthesia and analgesia on our results, there are two reasons that it is probably minor. First, as mentioned earlier, the timescales associated with the most consistent and distinctive state-space structure observed here corresponded to the timescales that are most important for carrying visual information in neurons in the visual cortex of the awake macaque [11]. Second, available evidence suggests that the effect of drug-induced changes in network state would be more likely to dilute any underlying structure, than to create it. Specifically, Ecker et al. [45] compared the variability of the spiking activity in V1 in awake monkeys with that seen in monkeys under sufentanil, inferring the presence of a sufentanil-related state variable. This state variable fluctuated with a timecourse that varied from 50 to 1650 ms. If similar sufentanil-related fluctuations were present during our recording sessions, the impact would have been distributed randomly across the many 320 ms spiking response samples.

Methodological innovations
In summary, we introduce a new framework for applying TDA to spike train data and use it to analyse patterns of spiking response in macaque visual cortex. There are two main methodological innovations. First, in contrast to most previous applications of TDA to neural data, we do not assume that the spike trains have a vector space embedding that induces distances and correlation measures between them; rather, the topological analysis is directly applied to dissimilarity measures (e.g. distances) that result from considering spike trains to be sequences of events-in this case, the Victor-Purpura spike train distances, which are non-Euclidean. A second innovation of the approach is the filtration-the sequence of simplicial complexes derived from the dissimilarities that are used to compute the Betti curves. Typically, an increasing filtration is used for TDA [2,5]: graphs are progressively filled in for pairs of points at greater and greater dissimilarities; here, we show that the decreasing filtration can reveal a clearer picture of the data's geometry and topology. Betti curves, especially when averaged over several computations, capture the statistical distribution of the persistent topological features (tunnels, voids, etc.) across a filtration, which serves as a topological descriptor of the underlying metric spaces. In contrast to the usual filtration of graphs by increasing weights and their associated clique complexes, the decreasing filtration is obtained by considering, in reverse order, the independence complex of each graph (i.e. the clique complex of the complement of the graph). Betti curves of low dimension (β 1 -β 3 ) of the decreasing filtration therefore carry information on the arrangements of high-order cliques of the original graphs, which is not related to the homology of the increasing filtration. Because the analysis is carried out for a sequence of distances parametrized by timescale, we are able to identify the range in which the geometric structure of the spiking response space is most distinctive: the range 5-50 s −1 , i.e. 20 to 200 ms. As noted earlier, this timescale corresponds to the temporal precision that is most informative for decoding visual information from spike trains. This matching of the timescale of response space geometry and the timescale of neural coding may be a general feature of brain networks, and we speculate that networks in other domains (e.g. motor planning, learning, decision-making, etc.) will behave similarly.

Experiments
All procedures followed the guidelines provided by the US National Institutes of Health and Weill Cornell Medical College Animal Care and Use Committee. Full details concerning the physiological preparation and multi-tetrode single-unit recordings can be found in [46] and [18]. Detailed descriptions of the visual stimuli, their generation and their display during the experimental sessions are given in [18], which also details how single-unit activity was characterized, and how the time-series of neural firing events, including the procedures utilized for spike sorting, were extracted from the multi-tetrode recordings. These methods are summarized in electronic supplementary material, materials and methods.

Visual stimuli and stimulation protocols
The visual stimuli used in this study are 16 × 16 black and white chequerboard-like patterns drawn from a mathematically defined stimulus space. The stimulus space has 10 coordinate axes specifying the type of the multi-point correlations in each 2 × 2 neighbourhood of 'checks' in the patterns, and the coordinates along these axes define the strength of the correlation (the fidelity with which the multi-point correlation is rendered across each pattern), and the correlation's polarity [17,18]. The 10 coordinate axes can be partitioned into four classes according to the order of the multi-point correlations [17]: first-, second-, third-and fourth-order correlations (see figure 1). In this study, the patterns were drawn from four coordinate values (2 magnitudes × 2 polarities) along each of the 10 axes in the stimulus space. The size of the 16 × 16 patterns, their position on the visual display screen, the orientation of the patterns and the check-sizes in the patterns were chosen to optimize stimulation of isolated clusters of neurons based on the neurons' receptive field properties as determined by the responses to sinusoidal gratings. Further details on this point are given in [18]. The topological data analysis presented here was performed on a subset of the neural responses obtained during the experimental runs. The experimental runs used (at least) 64 samples from the four coordinate values along each of the 10 axes in the stimulus space. In addition, 64 examples of fully random patterns were included in the experimental runs. The topological data analysis excluded the responses to the random patterns. For the experimental runs, a total of 2624 unique stimuli (64 × 41) each repeated 4 times (for a grand total of 10 496 stimuli) was shown. In 20 of the 28 analysed datasets, the 64 analysed responses to visual stimuli of the same type are part of a larger sequence of 128 responses to 128 unique stimuli; in these cases, we selected for our analysis the responses to the first 64 unique stimuli shown during the experiment. The patterns were presented in extended pseudorandom sequences at 100% contrast, surrounded by a mid-level grey background. Each pattern appeared for 0.32 s and was immediately followed by the next pattern in the sequence. A different pseudorandom sequence was used for each of the four times a sequence was run during a recording block. There was a pause of about 1 min between these repeats, during which time the grey background was displayed.

Analysed data
For each of the 28 datasets, the responses to all non-random stimulus types (each determined by a choice of an axis and a strength level for the correlations) for all repeats were then assembled to yield a total of 40 Â 4 ¼ 160 collections. To study the topology of the space that neural activity occupies, we selected the 80 collections in each dataset with the smallest number of empty responses (i.e. responses with 0 spikes). The majority of these collections had at least 60 non-empty responses; see electronic supplementary material, figure S1d. This choice for selecting collections is motivated by the fact that our analysis maps spike trains with at least one spike to distinct points in the topological space, while all empty spike trains are mapped to the same point. Thus, this selection criterion maximized the sampling of the space occupied by a collection.

Victor-Purpura distance
The (multi-neuron) Victor-Purpura distance [11] is a costbased distance between spike trains that has two parameters: a parameter q that controls the timescale used to quantify dissimilarity in spike timing, and a parameter k that determines the relevance of the neuron of origin of each spike. The distance between two neuronal responses is defined as the minimum cost of transforming one spike train into the other via a sequence of basic moves: addition or deletion of a single spike, with a cost of 1, shift of a single spike by a time interval Δt, with a cost q| Δt|, or change in the neuron of origin, with cost k. Thus, q = 0 corresponds to ignoring spike time but retaining spike count, q > 0 corresponds to considering temporal structure at a scale of 1/q. If a spike time in one train is within 1/q of a spike time in a second train, they are considered to correspond, as they contribute less than 1 unit to the dissimilarity. Conversely, if spikes are more than 2/q apart, they are considered to be unrelated, since shifting them into correspondence would incur a higher cost than deleting the spike from one train and inserting it into the other. Similarly, while for k = 0 changing the neuron of origin of a spike comes with no cost, for k = 1, the corresponding cost is equivalent to the insertion of a spike. Further background on the Victor-Purpura distance, as well as efficient dynamic-programming algorithms for calculating it, can be found in [7,47]. In our analysis, for each collection of neuronal responses, we determined the distance matrix D = (D ij ), where D ij is the Victor-Purpura distance with parameters (q, k) between the ith and the jth response in the collection. The parameters q and k ranged over a grid of values: q = 1, 2, 5, 10, 20, 50, 100, 200 (s −1 ), and k = 0, 1. The Victor-Purpura distance between spike trains is computed using the 'labdist_faster_qkpara_opt' Matlab function implemented by Thomas Kreuz, available at: www-users.med.cornell.edu/ jdvicto/labdist_faster_qkpara_opt.html.

TDA methods
We applied and extended the clique topology method introduced in [3], which corresponds to computing the Betti curves associated with our increasing filtration. Given a symmetric matrix M with zeros on its main diagonal, the first step of the method consists of considering the above-diagonal part of the matrix, transforming it by rank ordering its entries (thus replacing the original entry values by natural numbers 0, 1, …) and completing the below-diagonal part of the rank-ordered matrix by symmetry. The method then computes the persistent homology (see TDA software) of the rank-ordered matrix and finally determines the associated Betti curves. Since the input matrix is transformed by considering only the rank ordering of its entries, the output is invariant to monotonic transformation applied entry-wise to M. Given a dissimilarity or distance matrix, the clique topology method orders the entries increasingly. We observe that in [3] the method is also applied to a matrix C = (C ij ) of correlations by transforming it into a dissimilarity matrix D = (D ij ) via the application of any function that inverts the ordering between the absolute values |C ij | of the correlations and the entries D ij , e.g. D ij = −|C ij |. In this work, we apply the clique topology method by ordering the entries of our distance matrices both increasingly (increasing filtration) and decreasingly (decreasing filtration), the latter case corresponding to the weight rank clique filtration method introduced in [22]. A sequence of nested simple graphs G 0 , G 1 , · · · , G s , called a filtration, can be determined from a rank-ordered n × n symmetric matrix (typically representing pairwise distances) as follows. Each graph has the same set of n nodes, and the edges of a graph G k are determined by thresholding the rank-ordered matrix (setting to one all entries smaller than k, and setting to zero the remaining entries) and regarding it as an adjacency matrix. An edge is therefore present between nodes i and j of the graph G k if and only if the (i, j ) entry of the rank-ordered matrix is less than k. If the above-diagonal entries of the symmetric rank-ordered matrix are all different, the graph G 0 has no edges, the graph G 1 has one edge, corresponding to the smallest off-diagonal entry of the rank-ordered matrix, and so on. As the edges fill in, graphs can be enriched with higher-order connectivity information encoded by basic 'pieces' of different dimensions, called simplices. Specifically, a p-clique (subgraph of p all-to-all connected nodes) in a graph is regarded as a ( p − 1)-dimensional simplex: 2-cliques are the edges of the graph and are viewed as line segments (one-dimensional simplices), 3-cliques are viewed as triangles (two-dimensional simplices), 4-cliques are viewed as tetrahedra (three-dimensional simplices), etc. In this way, each graph becomes a simplicial complex, called the clique complex of the graph, in which arrangements of cliques can enclose 'holes' of different dimensions. The number of one-dimensional tunnels, two-dimensional voids and three-dimensional 'cavities', known as Betti numbers and denoted by β 1 , β 2 and β 3 respectively, are computed at each step of the filtration. The Betti numbers are viewed as functions of the edge density ρ, the number of filled-in edges divided by the number N = n(n − 1)/2 of potential edges. These functions (β 1 (ρ), β 2 (ρ) and β 3 (ρ)) are the Betti curves, which we computed over the range from ρ = 0 to 0.6 in our analysis.

TDA software
To compute clique topology and persistent homology, we use a faster and equivalent alternative to the original CliqueTop scripts (https://github.com/nebneuron/clique-top; see [3]), namely, Ripser [48] (https://github.com/Ripser/ripser) on a rank-royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220677 and EY07977 from the National Eye Institute of the NIH. K.P.P. acknowledges NIH grants no. EY034150 from the National Eye Institute and NS111019 from the National Institute of Neurological Disorders and Stroke. S.R. was supported by Ikerbasque (The Basque Foundation for Science). A.G. and S.R. were supported by the Basque Government through the BERC 2018-2021 programme and by the Ministry of Science, Innovation and Universities: BCAM Severo Ochoa accreditation SEV-2017-0718 and through project RTI2018-093860B-C21 funded by (AEI/FEDER, UE) with acronym 'MathNEURO'. A.G. acknowledges the support of the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation. M.D. and S.R. acknowledge suppport of Inria through the Associated Team 'NeuroTransfSF'.