Measures of spike train synchrony
Measures that estimate the degree of synchrony between spike trains are important tools for many applications. Among others, they can be used to quantify the reliability of neuronal responses upon repeated presentations of a stimulus (Mainen and Sejnowski, 1995) or to test the performance of neuronal models (Jolivet et al., 2008).
Time-scale dependent spike train distances
Arguably the most prominent task for spike train analysis is to address questions regarding the nature of the neuronal code (for an overview see Victor, 2005). This task is typically addressed by spike train distances which consider spike trains to be points in an abstract metric space and quantify their dissimilarity by non-negative values. Most of these spike train distances depend on a parameter which determines the temporal scale in the spike trains to which the distances are sensitive. While in one limit of the parameter range these distances are sensitive to the difference in spike number, they detect spike coincidences in the other limit. These limits reflect the characteristics of a rate code and a coincidence code, respectively.
In a typical experimental setup, single or multi-unit responses are recorded for repeated presentations of a set of stimuli. Neural coding can be assessed by applying a clustering analysis to the pairwise spike train distance matrices obtained for different sensitivities, e.g., different values of the time scale parameter. The time-scale \(\tau_d\) for which the responses to different stimuli are best distinguished (i.e., for which lowest distances are found between responses to the same stimulus and highest distances between responses to any two different stimuli) is assumed to be the discriminative precision of the neural code. The most widely used time-scale dependent measures are the Victor–Purpura distance (Victor and Purpura, 1996, 1997), the van Rossum distance (van Rossum, 2001) and the Schreiber et al. similarity measure (Schreiber et al., 2003). For comparison of these measures on simulated data, refer to Schrauwen and Campenhout, 2007, Paiva et al., 2010, and Chicharro et al., 2011.
Victor–Purpura distance
The spike train distance introduced in Victor and Purpura (1996, 1997) defines the distance between two spike trains in terms of the minimum cost of transforming one spike train into the other by means of just three basic operations ( Figure 1): spike insertion (cost \(1\)), spike deletion (cost \(1\)) and shifting a spike by some interval \(\Delta t\) (cost \(q|\Delta t|\)). The cost per time unit \(q\) sets the time scale of the analysis. For \(q = 0\) the distance is equal to the difference in spike counts, while for large \(q\) the distance approaches the number of non-coincident spikes, as it becomes more favorable to delete and reinsert all non-coincident spikes rather than shifting them. Thus, by increasing the cost, the distance is transformed from a rate distance to a temporal distance. The source code for the Victor-Purpura distance can be downloaded at [1].
In addition to this distance which is sensitive to the timing of individual spikes, two complementary cost-based distances have been proposed which are sensitive to interspike intervals and to temporal patterns of spikes (‘motifs’) (Victor and Purpura, 1997).
van Rossum distance
A second spike train distance was introduced in van Rossum (2001). Here the discrete spike trains \(X\) and \(Y\) are transformed into continuous functions by convolving each spike \(t_k\) with an exponential kernel
\[\tag{1} H(t)\exp(-\frac{t}{\tau_R})\ .\]
Here \(\tau_R\) is the time constant and \(H\) is the Heaviside step function with \(H(t)=0\) if \(t<0\) and \(H(t)=1\) if \(t \geq 0\ .\) This kernel shape is motivated by its causality and its resemblance to the shape of postsynaptic currents ( Figure 2). From the resulting waveforms \(\tilde{x}(t)\) and \(\tilde{y}(t)\ ,\) the van Rossum distance \(D_R\) can be calculated as
\[\tag{2} D_R (\tau_R) = \frac{1}{\tau_R} \int_0^\infty [ \tilde{x}(t) - \tilde{y}(t) ]^2 dt\ .\]
For this method, the time constant \(\tau_R\) of the exponential acts as the parameter that sets the time scale. It is inversely related to Victor and Purpura’s cost parameter, i.e., the temporal relationship between the spikes is evaluated for low \(\tau_R\ ,\) while for high \(\tau_R\) the distance is only sensitive to differences in rate.
For this latter case the interval containing relevant contributions to the integral generally extends beyond the end of the spike trains which can make the calculation time consuming. However, Paiva et al. (2009) have recently shown that the distance can be evaluated in terms of a computationally more efficient estimator which only scales with the product of the numbers of spikes in both spike trains.
Schreiber et al. similarity measure
In this approach (Schreiber et al., 2003) each spike \(\!t_k\) is convolved with a Gaussian filter
\[\tag{3} \frac{1}{\sqrt{2 \pi \sigma_S^2}} \exp(-\frac{t^2}{2 \sigma_S^2})\]
of width \(\sigma_S\) to form continuous signals \(\tilde{x}'(t)\) and \(\tilde{y}'(t)\) which are then normalized and cross correlated ( Figure 2):
\[\tag{4} C_S(\sigma_S) = \frac{\tilde{x}' \tilde{y}'}{|\tilde{x}'||\tilde{y}'|}\]
As for the van Rossum distance, a computationally efficient estimator has been proposed (Paiva et al., 2009).
In this case the width of the convolving Gaussian filter \(\sigma_S\) sets the time scale of interaction between the two spike trains. However, unlike \(D_V\)and \(D_R\ ,\) the inverted measure \(D_S = 1 - C_S\) cannot be used to distinguish between rate and temporal coding since it does not cover the range of time scales sensitivities from a coincidence detector to a rate code distance (Chicharro et al., 2011). Instead, in the limit \(\sigma_S = \infty\) it attains values close to zero regardless of the rate difference. Furthermore, even the applicability to estimate reliability is limited because of the individual normalization of the spike trains, which renders the measure sensitive only to similarities in the temporal modulation of the individual rate profiles but neglects differences in the absolute spike count.
Population extensions
In order to uncover if and how populations of neurons interact and cooperate to encode a sensory input, the Victor-Purpura and the van Rossum distances have been extended to measures that can estimate the dissimilarity between different responses recorded from a population of neurons (Aronov et al., 2003, and Houghton and Sen, 2008, respectively). Both of these extensions introduce a second parameter which describes the importance of distinguishing spikes fired in different cells by interpolating between the two extreme coding strategies for neuronal populations: the summed population (SP) code where for each response the spike trains from different neurons are superimposed before the distances between different responses are calculated, and the labeled line (LL) code where the distances between different responses are calculated separately for each neuron and then added.
Time-scale independent spike train distances
Complementary to the time-scale dependent approaches, in recent years spike train distances have been proposed which are parameter-free and time-scale-adaptive. While not allowing the functional characterization and precision analysis described above, single-valued methods give an objective and comparable estimate of neuronal variability. They can be preferable in applications to real data for which there is no validated knowledge about the relevant time scales. The computational cost is reduced since there is no need for parameter optimization. In fact, it is not at all guaranteed that there exists an optimal parameter. For example, spike trains that include different time-scales such as regular spiking and bursting might result in misleading conclusions, since any fixed parameter will misrepresent either one of these dynamics. Measures that do not rely on a time scale include event synchronization (Quian Quiroga et al., 2002), and the ISI- and SPIKE-distances (Kreuz et al., 2007a, 2011).
Event synchronization
Event synchronization (Quian Quiroga et al., 2002) acts as a coincidence detector and quantifies the level of synchrony from the number of quasi-simultaneous appearances of spikes. Similar to the previous methods, the temporal resolution can be adjusted by means of a coincidence window of fixed size \(\tau\ ,\) but there also exists a variant which is parameter- and scale-free since the maximum time lag \(\!\tau_{ij}\) up to which two spikes \(t_i^x\) and \(t_j^y\) are considered to be synchronous is adapted to the local spike rates according to
\[\tag{5} \tau_{ij} = \min \{t_{i+1}^x - t_i^x, t_i^x - t_{i-1}^x,t_{j+1}^y - t_j^y, t_j^y - t_{j-1}^y\}/2\ .\]
Denoting the respective number of spikes as \(S_x\) and \(S_y\ ,\) the number of appearances of a spike in \(X\) shortly after a spike in \(Y\) is given by
\[\tag{6} c (x|y) = \sum_{i=1}^{S_x} \sum_{j=1}^{S_y} J_{ij}\]
with
\[\tag{7} J_{ij} = \begin{cases} 1 & {\rm if} ~~ 0 < t_i^x - t_j^y \leq \tau_{ij} \\ 1/2 & {\rm if} ~~ t^x_i = t^y_j \\ 0 & {\rm else}. \end{cases}\]
With the opposite value \(c (y|x)\) defined accordingly, the normalized event synchronization is obtained as
\[\tag{8} Q = \frac {c (y|x) + c (x|y)} {\sqrt{S_x S_y}}\ .\]
Additionally, a directed variant
\[\tag{9} q = \frac {c (y|x) - c (x|y)} {\sqrt{S_x S_y}}\]
which is normalized between \(-1\) and \(1\) is able to characterize the relative delays between the events ( Figure 3). For a renormalization which accounts for simultaneous events expected just by chance refer to Kreuz et al. (2007b). The Matlab source code for calculating and visualizing event synchronization can be downloaded at [2].
ISI- and SPIKE-distance
The ISI- and SPIKE-distances \(D_I\) and \(D_S\) rely on instantaneous values in the sense that in a first step the sequences of discrete spike times are transformed into (quasi-)continuous temporal profiles \(X (t)\) with one value for each sampling point. For the ISI-distance these temporal profiles are derived from the interspike intervals, while in the case of the SPIKE-distance they are extracted from differences between the spike times of the two spike trains. Both distances are then defined as the temporal average of the respective time profile
\[\tag{10} D_X = \frac{1}{T} \int_{t=0}^T dt X (t)\;,\qquad X = I, S\]
The following equations show how to calculate the two different temporal profiles. For each neuron \(u = x,y\) one assigns to each time instant ( Figure 4) the time of the previous spike
\[\tag{11} t_{\mathrm {P}}^u (t) = \max(t_i^u | t_i^u \leq t) ,\]
and the time of the following spike
\[\tag{12} t_{\mathrm {F}}^u (t) = \min(t_i^u | t_i^u > t) ,\]
as well as the interspike interval
\[\tag{13} \nu_{\mathrm {ISI}}^u (t) = t_{\mathrm {F}}^u (t) - t_{\mathrm {P}}^u (t)\ .\]
ISI-distance
The ISI-distance (Kreuz et al., 2007a) is based on the instantaneous interspike intervals. To define a time-resolved, symmetric, and time-scale-adaptive measure of the relative firing rate pattern, the instantaneous ratio between \(\nu_{\mathrm {ISI}}^x\) and \(\nu_{\mathrm {ISI}}^y\) is calculated according to:
\[\tag{14} I (t) = \begin{cases} \nu_{\mathrm {ISI}}^x (t) / \nu_{\mathrm {ISI}}^y (t) - 1 & {\rm if} ~~ \nu_{\mathrm {ISI}}^x (t) \leq \nu_{\mathrm {ISI}}^y (t) \\ - (\nu_{\mathrm {ISI}}^y (t) / \nu_{\mathrm {ISI}}^x (t) -1) & {\rm otherwise}. \end{cases}\]
This quantity becomes \(0\) for identical ISI in the two spike trains, and approaches \(-1\) and \(1\ ,\) respectively, if the first or the second spike train is much faster than the other. Since all deviations from identical ISI count equally, the ISI-distance is calculated by temporal averaging over the absolute values \(|I (t)|\ .\)
SPIKE-distance
The ISI-distance is based on the relative length of simultaneous interspike intervals and is thus well-designed to quantify similarities in the neurons’ firing rate profiles. However, it is not optimally suited to track synchrony that is mediated by spike timing and in particular by changes in the fraction of coincident spikes. This particular kind of sensitivity is not only of theoretical importance but also of high practical relevance since coincidences of spikes have been proven to be of high prevalence in many different neuronal circuits. This issue is addressed by the SPIKE-distance which uniquely combines the properties of the ISI-distance with a specific focus on spike timing (Kreuz et al., 2011).
The instantaneous differences of previous and following spike times are denoted as
\[\tag{15} \Delta t_{\mathrm {P}} (t) = t_{\mathrm {P}}^x (t) - t_{\mathrm {P}}^y (t)\]
and
\[\tag{16} \Delta t_{\mathrm {F}} (t) = t_{\mathrm {F}}^x (t) - t_{\mathrm {F}}^y (t)\ ,\]
respectively.
An instantaneous spike time based measure of spike train distance is given by the regular average of the two absolute differences of previous and following spike times. In order to achieve time-scale invariance (i.e., all stretched and compressed spike trains should yield the same value) this average spike time difference is divided by the mean interspike interval. This function is piecewise constant and takes just one value within each interval of the pooled spike train. In order to be more local in time, a locally weighted average of the two differences \(|\Delta t_{\mathrm {j}} (t)|, \mathrm {j = P, F}\) is employed such that the difference of the spikes that are closer dominate. To this aim, for each neuron \(u = x, y\) the intervals from the time instant under consideration to the previous and the following spikes are denoted as
\[\tag{17} \nu_{\mathrm {P}}^u (t) = t - t_{\mathrm {P}}^u (t)\]
and
\[\tag{18} \nu_{\mathrm {F}}^u (t) = t_{\mathrm {F}}^u (t) - t\ ,\]
respectively. Inserting the inverse of the mean intervals \(\left \langle \nu_{\mathrm {j}}^u (t) \right \rangle_u, \mathrm {j = P, F}\) as weights \(f\) in the locally weighted average
\[\tag{19} \left \langle \Delta t_j (t) \right \rangle_u, \mathrm {j = P, F} = \frac{\sum_{\mathrm {j = P, F}} |\Delta t_j (t)| f (\nu_j^u (t))}{\sum_{\mathrm {j = P, F}} f (\nu_j^u (t))}\]
and making use of
\[\tag{20} \left \langle \nu_{\mathrm {P}}^u (t) \right \rangle_u + \left \langle \nu_{\mathrm {F}}^u (t) \right \rangle_u = \left \langle \nu_{\mathrm {ISI}}^u (t) \right \rangle_u\]
yields
\[\tag{21} S (t) = \frac{ |\Delta t_{\mathrm {P}} (t)| \left \langle \nu_{\mathrm {F}}^u (t) \right \rangle_u + |\Delta t_{\mathrm {F}} (t)| \left \langle \nu_{\mathrm {P}}^u (t) \right \rangle_u}{\left \langle \nu_{\mathrm {ISI}}^u (t) \right \rangle_u^2}\ .\]
Both the ISI- and the SPIKE-distance are bounded in the interval \([0, 1]\ .\) For the latter the limit value \(0\) is obtained only for perfectly identical spike trains while for the former it is also obtained for periodic spike trains with the same period.
The ISI- and the SPIKE-distance are conceptually simple, computationally efficient and easy to visualize in a time-resolved manner. By taking into account only the previous and the following spike in each spike train these distances rely on local information only. They are also time-scale adaptive since the information used is not contained within a window of fixed size but rather within a time frame whose size depends on the local rate of each spike train.
The Matlab source code for calculating and visualizing the ISI- and the SPIKE-distances can be downloaded at [3].
Multi-neuron-extensions
For both the ISI-and the SPIKE-distance, there exist two multivariate extensions which estimate the time-resolved level of dissimilarity within a group of spike trains (Kreuz et al., 2009, 2011). The first one is based on an average over all pairs of spike trains (which can, in principle, be carried out for all measures, but for the ISI-and the SPIKE-distance the averaging can be performed locally), the second one is based on the standard deviation and thus truly multivariate.
Comparison of measures
One of the main arguments for the use of time-scale-dependent measures of spike train (dis)similarity is their potential insight into the precision of the neuronal code (Victor and Purpura, 1996). This argument has recently been reevaluated in Chicharro et al. (2011). According to this study the optimal time-scale obtained from the cluster analysis is far from being conclusive. Rather it results in a non-trivial way from the interplay of many different factors such as the distribution of the information contained in different parts of the response and the degree of redundancy between them.
Despite these problems in the interpretation of the optimal timescale, the Victor–Purpura and the van Rossum distance are designed to test for neuronal codes ranging from a rate code to a coincidence detector. This is a level of generality somewhere in between two extremes. While some methods evaluate a very specific coding hypothesis (e.g., the classical correlation coefficient based on binning which focuses purely on coincidences), other methods are more general (e.g., the ISI- and the SPIKE-distances which are time-scale-adaptive and parameter-free). Measures on different ends of this scale are complementary in nature. If a particular coding scheme is assumed, specific measures are needed for a confirmatory analysis, otherwise more general measures are very well suited for an exploratory analysis (Kreuz et al., 2011).
References
- Aronov D, Reich DS, Mechler F, Victor JD (2003). Neural coding of spatial phase in V1 of the macaque monkey. J Neurophysiol 89:3304–3327.
- Chicharro D, Kreuz T, Andrzejak RG (2011). What can spike train distances tell us about the neural code? J Neurosci Methods 199, 146-65 (2011).
- Houghton C, Sen K (2008). A new multineuron spike train metric. Neural Comput 20:1495–1511.
- Jolivet R, Kobayashi R, Rauch A, Naud R, Shinomoto S, Gerstner W (2008). A benchmark test for a quantitative assessment of simple neuron models. J Neurosci Methods 169:417–424.
- Kreuz T, Haas JS, Morelli A, Abarbanel HDI, Politi A (2007a). Measuring spike train synchrony. J Neurosci Methods 165:151–161.
- Kreuz T, Kraskov A, Andrzejak RG, Mormann F, Lehnertz K, Grassberger P (2007b) Measuring synchronization in coupled model systems: a comparison of different approaches. Phys D 225:29–42.
- 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 (2009). Time-resolved and time-scale adaptive measures of spike train synchrony. J Neurosci Methods 195:92–106.
- Mainen Z, Sejnowski T (1995). Reliability of spike timing in neocortical neurons. Science 268:1503–1506.
- Paiva ARC, Park I, Principe JC (2009). A reproducing kernel hilbert space framework for spike train signal processing. Neural Computation 21:424–449.
- Paiva ARC, Park I, Principe JC (2010). A comparison of binless spike train measures. Neural Computing and Applications, 19:405–419.
- Quian Quiroga R, Kreuz T, Grassberger P (2002). Event synchronization: a simple and fast method to measure synchronicity and time delay patterns. Phys Rev E 66:041904.
- Schrauwen B, Campenhout JV (2007). Linking non-binned spike train kernels to several existing spike train metrics. Neurocomputing 70:1247–1253.
- Schreiber S, Fellous JM, Whitmer JH, Tiesinga PHE, Sejnowski TJ (2003). A new correlation based measure of spike timing reliability. Neurocomputing 52:925–931.
- van Rossum MCW (2001). A novel spike distance. Neural Comput 13:751–763.
- Victor JD, Purpura KP (1996). Nature and precision of temporal coding in visual cortex: a metric-space analysis. J Neurophysiol 76:1310–1326.
- Victor JD, Purpura KP (1997). Metric-space analysis of spike trains: theory, algorithms and application. Network 8:127–164.
- Victor JD (2005). Spike train metrics. Current Opinion in Neurobiology 15:585–592.
Internal references
- David Golomb (2007) Neuronal synchrony measures. Scholarpedia, 2(1):1347.
- James Meiss (2007) Dynamical systems. Scholarpedia, 2(2):1629.
- Arkady Pikovsky and Michael Rosenblum (2007) Synchronization. Scholarpedia, 2(12):1459.
- Thomas Kreuz (2011) Measures of neuronal signal synchrony. Scholarpedia, 6(12):11922.
