# SPIKE-distance

Post-publication activity

Curator: Thomas Kreuz

The SPIKE-distance is an estimator of the dissimilarity between two (or more) spike trains. In contrast to most other spike train distances (such as the Victor-Purpura distance) it is time-resolved and is able to track changes in instantaneous clustering, i.e., time-localized patterns of (dis)similarity among two or more spike trains. Additional features include selective and triggered temporal averaging as well as the instantaneous comparison of spike train groups. The SPIKE-distance can also be formulated as a causal measure which is defined such that the instantaneous values of dissimilarity rely on past information only so that time-resolved spike train synchrony can be estimated in real-time.

The final definition presented here is the one introduced in Kreuz et al., 2013, which improves considerably on the original proposal (Kreuz et al., 2011). Mathematical properties and expectation values for Poisson spike train can be found in Mulansky et al., 2015.

## Bivariate SPIKE-distance

The bivariate SPIKE-distance $$D_S$$ (Kreuz et al., 2013) relies on instantaneous values in the sense that in a first step the two sequences of discrete spike times are transformed into a (quasi-) continuous measure profile $$S (t)$$, that is, a temporal sequence of instantaneous dissimilarity values. The distance is then defined as the temporal average of the measure profile.

The instantaneous dissimilarity values $$S (t)$$ are derived from differences between the spike times of the two spike trains. First, for each neuron $$n = 1,2$$ one assigns to each time instant (Figure 1) three piecewise constant quantities. These are the time of the previous spike

$t_{\mathrm {P}}^{(n)} (t) = \max(t_i^{(n)} | t_i^{(n)} \leq t) ,$

and the time of the following spike

$t_{\mathrm {F}}^{(n)} (t) = \min(t_i^{(n)} | t_i^{(n)} > t) ,$

as well as the interspike interval

$x_{\mathrm {ISI}}^{(n)} (t) = t_{\mathrm {F}}^{(n)} (t) - t_{\mathrm {P}}^{(n)} (t)\ .$

Figure 1: Illustration of the local quantities (relative to time instant t) needed to calculate the instantaneous dissimilarity values on which the SPIKE-distance and its real-time variant are based. Modified from Kreuz et al., 2012.

The ambiguity regarding the definition of the very first and the very last interspike interval is resolved by adding to each spike train an auxiliary leading spike at time $$t = 0$$ (the beginning of the recording) and an auxiliary trailing spike at time $$t = T$$ (the end of the recording).

The instantaneous dissimilarity values are calculated in two steps: First for each spike the distance to the nearest spike in the other spike train is calculated, then for each time instant the local spike time differences are selected, weighted, and normalized. Each time instant is uniquely surrounded by four corner spikes: the preceding spike of the first spike train $$t_{\mathrm {P}}^{(1)}$$, the following spike of the first spike train $$t_{\mathrm {F}}^{(1)}$$, the preceding spike of the second spike train $$t_{\mathrm {P}}^{(2)}$$, and, finally, the following spike of the second spike train $$t_{\mathrm {F}}^{(2)}$$. Each of these corner spikes can be identified with a spike time difference to the nearest spike in the other spike train, for example, for the previous spike of the first spike train

$\Delta t_{\mathrm {P}}^{(1)} (t) = \min_i ( | t_{\mathrm {P}}^{(1)} (t) - t_i^{(2)} | ),$

and analogously for $$t_{\mathrm {F}}^{(1)}$$, $$t_{\mathrm {P}}^{(2)}$$, and $$t_{\mathrm {F}}^{(2)}$$.

Figure 2: Exemplary dissimilarity profiles for the SPIKE-distance. A. Bivariate example: Varying spike matches. B. Multivariate example with 50 spike trains: In the first half within the noisy background there are four regularly spaced spiking events with increasing jitter. The second half consists of ten spiking events with decreasing jitter but now without any noisy background. Modified from Kreuz et al., 2012.

For each spike train separately a locally weighted average is employed such that the differences for the closer spike dominate; the weighting factors depend on

$x_{\mathrm {P}}^{(n)} (t) = t - t_{\mathrm {P}}^{(n)} (t)$

and

$x_{\mathrm {F}}^{(n)} (t) = t_{\mathrm {F}}^{(n)} (t) - t\ ,$

the intervals from the time instant under consideration to the previous and the following spikes for each neuron $$n = 1, 2$$. The local weighting for the spike time differences of the first spike train reads

$S_1 (t) = \frac{ \Delta t_{\mathrm {P}}^{(1)} (t) x_{\mathrm {F}}^{(1)} (t) + \Delta t_{\mathrm {F}}^{(1)} (t) x_{\mathrm {P}}^{(1)} (t)}{x_{\mathrm {ISI}}^{(1)} (t)}\ .$

and analogously $$S_2 (t)$$ is obtained for the second spike train. Averaging over the two spike train contributions and normalizing by the mean interspike interval yields

$S' (t) = \frac{ S_1 (t) + S_2 (t)}{2 \left \langle x_{\mathrm {ISI}}^{(n)} (t) \right \rangle_n}.$

This quantity weights the spike time differences for each spike train according to the relative distance of the corner spike from the time instant under investigation. This way relative distances within each spike train are taken care of, while relative distances between spike trains are not yet. In order to get these ratios straight, in a last step the two contributions from the two spike trains are locally weighted by their instantaneous interspike intervals. This yields the dissimilarity profile

$S (t) = \frac{ S_1 (t) x_{\mathrm {ISI}}^{(2)} (t) + S_2 (t) x_{\mathrm {ISI}}^{(1)} (t)}{2 \left \langle x_{\mathrm {ISI}}^{(n)} (t) \right \rangle_n^2}.$

Finally, integrating over time leads to the bivariate SPIKE-distance:

$D_S = \frac{1}{T} \int_{t=0}^T S (t) dt\;.\qquad$

The SPIKE-distance is bounded in the interval $$[0, 1]$$. The value zero is only obtained for perfectly identical spike trains.

An example plot of a dissimilarity profile obtained for two spike trains is shown in Figure 2A. The observed monotonous increase of the instantaneous values followed by a monotonous decrease mirrors exactly the actual change in the match between the spike times of the two spike trains.

## Averaged bivariate SPIKE-distance

In the case of multivariate datasets consisting of a larger number of spike trains ($$N > 2$$) it is convenient to average the bivariate SPIKE-distance over all pairs of spike trains to obtain the averaged bivariate SPIKE-distance. The same kind of time-resolved visualization as in the bivariate case is possible, because the two averages over time and over pairs of spike trains commute. One thus can first calculate the instantaneous average $$S^a (t)$$ over all pairwise instantaneous spike differences $$S^{mn} (t)$$

$S^a (t) = \frac{1}{N(N-1)/2}\sum_{n=1}^{N-1} \sum_{m=n+1}^N S^{mn} (t)$

and then average over time

$D_S^a = \frac{1}{T} \int_{t=0}^T S^a (t) dt\;.\qquad$

Like the bivariate SPIKE-distance, the averaged bivariate SPIKE-distance is restricted to the interval $$[0, 1]$$.

Figure 2B depicts an examplary dissimilarity profile for a multivariate example. In the first half high instantaneous values reflect the background noise. The presence of four spiking events with increasing jitter within this noisy background is indicated by less and less pronounced drops in the dissimilarity profile. In the second half in which there is no background noise, the evenly spaced firing events with increasing precision are correctly reflected by a rather monotonous decrease of the instantaneous synchrony level which reaches zero at the perfectly synchronous event in the end.

## Real-time SPIKE-distance

Figure 3: Real-time SPIKE-distance: Exemplary dissimilarity profiles including causal moving averages. The spike trains are the same as in Figure 2. Modified from Kreuz et al., 2012.

The real-time SPIKE-distance (Kreuz et al., 2012) is a modification of the SPIKE-distance with the key difference that the corresponding time profile $$S_r (t)$$ can be calculated online because it relies on past information only. From the perspective of an online measure, the information provided by the following spikes, both their position and the length of the interspike interval, is not yet available. Like the regular SPIKE-distance $$D_S$$, this causal variant is also based on local spike time differences but now only two corner spikes are available, and the spikes of comparison are restricted to past spikes, e.g., for the preceding spike of the first spike train

$\Delta t_{\mathrm {P}}^{(1)} (t) = \min_i ( | t_{\mathrm {P}}^{(1)} (t) - t_i^{(2)} | ), t_i < t$

Figure 4: Real-time SPIKE-distance: Multivariate example with two spiking events which have identical first halves but different continuations. At the end of the first half some of the neurons have just spiked while others have been silent for some time, and for both events this large variability is correctly reflected by high values of the averaged bivariate dissimilarity profile $$S_r^a (t)$$. At these points begin the shaded areas which mark the spike information that is not yet available. Here for the first event the initial spiking in some of the spike trains turns into a reliable spiking event (one spike for each spike train) which results in a rapid decrease to very low values. For the second event this does not happen and accordingly the decrease is less pronounced. Modified from Kreuz et al., 2012.

Since there are no following spikes available, there is no local weighting, and since there is no interspike interval, the normalization is achieved by dividing the average corner spike difference by twice the average time interval to the preceding spikes. This yields a causal indicator of local spike train dissimilarity:

$S_r (t) = \frac{ \Delta t_{\mathrm {P}}^{(1)} (t) + \Delta t_{\mathrm {P}}^{(2)} (t)}{4 \left \langle x_{\mathrm {P}}^{(n)} (t) \right \rangle_n}\ .$

Again the averaged bivariate version is obtained by averaging over the dissimilarity profiles for all pairs of spike trains.

Example plots of the dissimilarity profiles obtained for the bivariate and the multivariate example already used in Figure 2 are shown in Figure 3. As can be seen for the bivariate example (Figure 3A), any spike time difference is considered most relevant right at the later of two spikes when $$S_r(t)$$ goes back to a local maximum value. In the case where the two preceding spikes are closest to each other, it goes back to its maximum value of one, since at these points the mean time interval to the two preceding spikes is exactly half their difference. Any successive period of common non-spiking leads to a decrease of the instantaneous distance values. This is a desired property since common non-spiking is as much a sign of synchrony as common spiking. The decrease is hyperbolic and its slope depends on the preceding spike time difference.

Similar short time-scale fluctuations in the dissimilarity profile of the real-time SPIKE-distance can also be observed in the multivariate example (Figure 3B). The dissimilarity peaks are due to the lack of knowledge about future spiking behaviour, they occur at time instants when it is not yet known whether there will in fact be a reliable spiking event or not. They can be eliminated by means of an appropriate (causal) moving average. This introduces a parameter, the size of the moving window, which is set depending on the time-scale of interest, e.g., short-term behavior or long-term trends.

The local uncertainty can only be resolved when more and more information becomes available. This is illustrated in Figure 4 with two spiking events which are identical (one spike per spike train) except for the omitted second half of the second event. While for both events the instantaneous values in the first part necessarily have to be identical (and very high since only some of the neurons have recently spiked), the differences in the second part are clearly reflected by different continuations of the dissimilarity profile.

## Representations

### Full matrix and cross sections

Figure 5: Instantaneous clustering for artificially generated spike trains whose clustering changes every 500 ms (dashed black lines). The matrices of pairwise instantaneous values for both the regular (top) and the real-time SPIKE-distance (bottom) reflect the changing cluster association as exemplified for the four time instants marked by the green lines. Modified from Kreuz et al., 2012.

For both the regular and the real-time SPIKE-distance there are several levels of information reduction (Kreuz et al., 2012). The starting point is the most detailed representation in which one instantaneous value is obtained for each pair of spike trains. This results in a matrix of size ’number of sampled time instants’ × ’squared number of spike trains’ (i.e. $$\#(t_n)N^2$$).

From this matrix it is possible to extract any information desired. By selecting a pair of spike trains one obtains the bivariate dissimilarity profile for this pair of spike trains (as shown in Figure 2A and Figure 3A). Selecting a time instant yields an instantaneous matrix of pairwise spike train dissimilarities. This matrix can be used to divide the spike trains into instantaneous clusters, i.e., groups of spike trains with low intra-group and high intergroup dissimilarity.

Examples for four different time instants are shown in Figure 5 in which artificially generated spike trains exhibit a different clustering behavior every 500 ms (Figure 5A). This varying clustering structure is correctly reflected in the pairwise dissimilarity matrices (Figure 5B) of both the regular and the real-time SPIKE-distance. The main difference between the two measures is clearly visible in the first column where the regular SPIKE-distance averages over past and future behavior and thus superimposes the checkered pattern of the first interval with the more disordered clustering of the second interval. The real-time variant is not yet aware of the latter interval and reflects the checkered pattern of the past interval only. Another difference can be seen in the second column where the real-time SPIKE-distance exhibits the large instantaneous values obtained during uncompleted firing events.

Figure 6: Selective temporal averaging. Each column (from left to right) depicts the selective temporal average over the intervals marked by the horizontal lines in A (from top to bottom). Modified from Kreuz et al., 2012.

### Spatial and temporal averaging

Another way to reduce the information of the dissimilarity matrix is averaging. There are two possibilities which commute: the spatial average over spike train pairs and the temporal average. As could be seen in Figure 2B and Figure 3B, the local average over spike train pairs yields a dissimilarity profile for the whole population. Temporal averaging on the other hand leads to a bivariate distance matrix; examples are shown in Figure 6. Finally, in both cases application of the respective remaining average results in one distance value which describes the overall level of synchrony for a group of spike trains over a given time interval.

Figure 7: Internally triggered temporal averaging. A. Poisson spike trains with superimposed firing patterns: Five spike trains follow the spiking of the first spike train (with small amounts of jitter). Below the horizontal line denoting the overall average, triangles mark the spike times of the first spike train which are used as trigger instants. B. Dissimilarity matrices (for the regular SPIKE-distance only): Averaging over the whole interval (left) and averaging over the instantaneous matrices at all trigger instants only (right). C. Hierarchical cluster trees (dendrograms) obtained from the dissimilarity matrices in B. Modified from Kreuz et al., 2012.

### Triggered averaging

A further option is triggered temporal averaging (Kreuz et al., 2012). Here the matrices are averaged over certain trigger time instants only. The idea is to check whether this triggered temporal average is significantly different from the global average since this would indicate that something peculiar is happening at these trigger instants. These trigger times can either be obtained from internal conditions (such as the spike times of a certain spike train) or from external influences (such as the occurrence of certain features in a stimulus).

#### Internal triggering

An example of internal triggering can be seen in Figure 7. In this artificially generated setup there are $$20$$ simultaneously recorded neurons and almost all of them fire independently from each other following a Poisson statistic (Figure 7A). The exception is the first neuron which fires at a lower rate and is assumed to have a strong excitatory synaptic coupling to five of the other neurons (# 4, 8, 11, 16, 19). Correspondingly, the spike trains of these neurons were modified such that they contained slightly lagged and jittered copies of the spikes of the first spike train in addition to spikes generated independently. This represents a situation in which each spike in the first spike train causes (triggers) a spike in these spike trains but there are also other spikes (which might have been caused by different spike trains).

The task is to identify these five neurons. This is very difficult via visual inspection, and also the overall temporal average is unable to do so (Figure 7B, left column). However, these neurons can be identified by averaging over the the pairwise instantaneous values obtained for the spike times of the first spike train (the internal trigger instants) only. The resulting dissimilarity matrix shown in Figure 7B (right column) includes an irregular grid of very small distance values. In Figure 7C another representation of dissimilarity matrices is used, a hierarchical cluster tree also known as dendrogram. While the dendrogram obtained from the average dissimilarity matrix (left) does not fall into separate clusters, in the dendrogram of the triggered average (right) the five modified spike trains form one distinct cluster with the first spike train and can thus easily be identified.

Figure 8: Externally triggered temporal averaging. A. Artificially generated spike trains. While half of the spike trains are noisy, the synchrony of the other half is modulated via a chirp-like external stimulus (shown on top). These non-periodically varying levels of synchrony can be traced by triggering on time instants with common stimulus amplitude (marked by horizontal green lines and green triangles). B/C. Each column (from left to right) depicts the dissimilarity matrix (B) and the corresponding dendrogram (C) of the externally triggered temporal average for decreasing amplitudes of the chirp function (from top to bottom). As the emerging spike train cluster shows lower stimulus amplitudes leads to increased levels of synchrony. Modified from Kreuz et al., 2012.

#### External triggering

In contrast to internal triggering, externally triggered averaging allows certain (external) stimulus features to be related with spike train synchrony and might, thus, be a promising tool for the investigation of neuronal coding. While the example for internal triggering assumes a simultaneously recorded population of neurons, in the setup of the external triggering example (Figure 8) just one neuron is recorded for repeated stimulation with the same stimulus. This stimulus is a chirp-function, representative of a non-periodic time-varying stimulus. It is assumed that the neuron is sensitive to negative amplitudes and accordingly it exhibits higher (lower) reliability for local minima (maxima) of the chirp function (Figure 8A). However, in order to better illustrate the gradual increase in clustering half the trials were left unaffected.

As the amplitude of the chirp function varies so does the spike train synchrony of half of the trials, but due to the non-periodicity of the stimulus this is quite difficult to detect. Detection is facilitated by externally triggered averaging where the triggering is performed on common stimulus features, here the amplitude of the chirp function. As can be seen in the dissimilarity matrices (Figure 8B) and even better in the dendrogram (Figure 8C), a decrease of this stimulus amplitude leads to the emergence of a spike train cluster consisting of the modulated spike trains which indicates their increase in reliability.

Figure 9: Screenshot from a movie (see Supplementary Material). A. Artificially generated spike trains. B. Dissimilarity matrices obtained by averaging over two separate time intervals for both the regular and the real-time SPIKE-distance as well as their averages over subgroups of spike trains. C. Corresponding dendrograms. Modified from Kreuz et al., 2012.

### Comparison of spike train groups and movie

The capability to track spike train synchrony can be demonstrated best in a movie (see external links). This movie uses the artificially generated spike trains from Figure 5 and Figure 6 and includes instantaneous clustering, selective temporal averaging of individual or combined intervals, several examples of triggered averaging as well as the corresponding dendrograms. As can be seen in the screenshot of the movie shown in Figure 9, the movie includes one additional feature, the comparison of spike train groups, where the spike trains are manually assigned to subgroups, and a block matrix (and the corresponding dendrogram) is obtained by averaging over the respective submatrices of the original dissimilarity matrix.

## Potential applications in neuroscience

In some situations it is sufficient to evaluate spike train synchrony at a rather low temporal resolution, e.g., by means of a moving-window analysis where the level of synchronization within a certain interval is compressed into a single value and then compared for successive intervals. On the other hand, many applications require a high temporal resolution, e.g., in order to evaluate the information transfer between synaptically coupled neurons (Reyes, 2003), to detect replay of precisely timed sequential patterns of neural activity (Ji and Wilson, 2007), or to track spike train response variability within a neuronal population (Mitchell et al., 2007; Kreuz et al., 2009).

The fact that for the SPIKE-distance there are no limits to the temporal resolution allows further analyses such as selective and triggered temporal averaging. In the latter case the trigger times could be obtained from internal conditions (such as the spike times of a certain neuron which might help to detect converging or diverging patterns of firing propagation) or from external influences (such as the occurrence of certain features in a stimulus which allows to address questions of neuronal coding). Finally, spatial averaging over spike train groups becomes possible. In applications to real data, these groups could be different neuronal populations or responses to different stimuli, etc.

The capability to track the synchrony between two or more spike trains opens up new possibilities for several important applications. Synchronization has been hypothesized to play a pivotal role in neuronal coding (Miller and Wilson, 2008), and thus a real-time tracker of spike train synchrony could be an essential tool for a rapid online decoding with brain-machine interfaces in order to control prosthetic devices (Hochberg et al., 2006). In epilepsy patients, monitoring the spiking activity of large ensembles of single neurons could lead to a better understanding of the mechanisms of seizure generation, propagation, and termination (Bower et al., 2012; Schevon et al., 2012). Furthermore, if neuronal spiking patterns turn out to be specific predictors of seizure occurrence as reported by a recent study (Truccolo et al., 2011), the real-time SPIKE-distance could be a viable tool for the implementation of a prospective seizure prediction algorithm (Mormann et al., 2007).

In a similar manner the SPIKE-distance and its real-time variant can be applied to monitor synchrony in continuous data (which are first transformed into discrete data). Even the analysis of mixed continuous-discrete signals is possible (see also Andrzejak and Kreuz, 2011). A potential application could be the combined analysis of discrete spike trains and continuous neuronal oscillations. In particular, it would be interesting to investigate neuronal synchronization patterns in dependence of the phase of the local field potential, a scenario reminiscent of the ’neuronal communication through neuronal coherence’ scenario (Fries, 2005).

## Source codes

Main article: SPIKY

SPIKY, a graphical user interface written in Matlab that facilitates the application of time-resolved measures of spike train synchrony (including the SPIKE-distance) to both simulated and real data can be found on the SPIKY download page. The same source code can also be used to calculate and visualize the ISI-distance (Kreuz et al., 2007) and SPIKE synchronization (Kreuz et al., 2015), two measures that are complementary to the SPIKE-distance. On this webpage there are also two movies which demonstrate the SPIKE-distance best.

## References

• Andrzejak RG, Kreuz T (2011). Characterizing unidirectional couplings between point processes and flows. Europhysics Letters 96:50012.
• Bower MR, Stead M, Meyer FB, Marsh WR, Worrell GA (2012). Spatiotemporal neuronal correlates of seizure generation in focal epilepsy. Epilepsia 53:807-16.
• Fries P (2005). A mechanism for cognitive dynamics: neuronal communication through neuronal coherence. Trends in Cognitive Sciences 9:474-80.
• Hochberg LR, Serruya MD, Friehs GM, Mukand JA, Saleh M, Caplan AH, Branner A, Chen D, Penn RD, Donoghue JP (2006). Neuronal ensemble control of prosthetic devices by a human with tetraplegia. Nature 442:164-71.
• Ji D, Wilson MA (2007). Coordinated memory replay in the visual cortex and hippocampus during sleep. Nature Neuroscience 10:100-7.
• Kreuz T, Haas JS, Morelli A, Abarbanel HDI, Politi A (2007). Measuring spike train synchrony. J Neurosci Methods 165:151–161.
• Kreuz T, Chicharro D, Andrzejak RG, Haas JS, Abarbanel HDI (2009). Measuring multiple spike train synchrony. J Neurosci Methods 183:287–299.
• Kreuz T, Chicharro D, Greschner M, Andrzejak RG (2011). Time-resolved and time-scale adaptive measures of spike train synchrony. J Neurosci Methods 195:92–106.
• Kreuz T, Chicharro D, Houghton C, Andrzejak RG, Mormann F (2013). Monitoring spike train synchrony. JNeurophysiol 109:1457-72.
• Kreuz T, Mulansky M, Bozanic N (2015). SPIKY: A graphical user interface for monitoring spike train synchrony. JNeurophysiol 113, 3432.
• Miller EK, Wilson MA (2008). All My Circuits: Using Multiple Electrodes to Understand Functioning Neural Networks. Neuron 60:483-8.
• Mitchell JF, Sundberg KA, Reynolds JH (2007). Differential attention-dependent response modulation across cell classes in macaque visual area V4. Neuron 55:131-41.
• Mormann F, Andrzejak RG, Elger CE, Lehnertz K (2007). Seizure prediction: the long and winding road. Brain 130:314-33.
• Mulansky M, Bozanic N, Sburlea A, Kreuz T (2015). A guide to time-resolved and parameter-free measures of spike train synchrony. IEEE Proceeding on Event-based Control, Communication, and Signal Processing (EBCCSP), 1-8.
• Reyes AD (2003). Synchrony-dependent propagation of firing rate in iteratively constructed networks in vitro. Nature Neurosci 6:593–599.
• Schevon CA, Weiss SA, McKhann G Jr, Goodman RR, Yuste R, Emerson RG, Trevelyan AJ (2012). Evidence of an inhibitory restraint of seizure activity in humans. Nat Commun 3:1060.
• Truccolo W, Donoghue JP, Hochberg LR, Eskandar EN, Madsen JR, Anderson WS, Brown EN, Halgren E, Cash SS (2011). Single-neuron dynamics in human focal epilepsy. Nature Neurosci 14:635-41.

Internal references

• Arkady Pikovsky and Michael Rosenblum (2007) Synchronization. Scholarpedia, 2(12):1459.
• Nebojsa Bozanic, Mario Mulansky, Thomas Kreuz (2014) SPIKY. Scholarpedia, 9(12):32344.