SPIKY: A graphical user interface for monitoring spike train synchrony

Techniques for recording large-scale neuronal spiking activity are developing very fast. This leads to an increasing demand for algorithms capable of analyzing large amounts of experimental spike train data. One of the most crucial and demanding tasks is the identification of similarity patterns with a very high temporal resolution and across different spatial scales. To address this task, in recent years three time-resolved measures of spike train synchrony have been proposed, the ISI-distance, the SPIKE-distance, and event synchronization. The Matlab source codes for calculating and visualizing these measures have been made publicly available. However, due to the many different possible representations of the results the use of these codes is rather complicated and their application requires some basic knowledge of Matlab. Thus it became desirable to provide a more user-friendly and interactive interface. Here we address this need and present SPIKY, a graphical user interface which facilitates the application of time-resolved measures of spike train synchrony to both simulated and real data. SPIKY includes implementations of the ISI-distance, the SPIKE-distance and SPIKE-synchronization (an improved and simplified extension of event synchronization) which have been optimized with respect to computation speed and memory demand. It also comprises a spike train generator and an event detector which makes it capable of analyzing continuous data. Finally, the SPIKY package includes additional complementary programs aimed at the analysis of large numbers of datasets and the estimation of significance levels.


Introduction
Spike train distances are measures of the degree of synchrony between spike trains which yield low values for very similar and high values for very dissimilar spike trains.They are applied in two major scenarios: simultaneous and successive recordings.
The first scenario is the simultaneous recording of a neuronal population, typically in a spatial multi-channel setup.If different neurons emit spikes at the same time, these spikes are truly 'synchronous' (Greek: 'occurring at the same time').Synchronization between individual neurons has been proven to be of high prevalence in many different neuronal circuits (Tiesinga et al., 2008;Shlens et al., 2008).As of now, many open questions remain regarding the spatial scale and the nature of interactions (pairwise or higher order, see Nirenberg and Victor, 2007) as well as their functional significance for neuronal coding and information processing (Kumar et al., 2010).
In the second scenario the neuronal spiking response is recorded in different time intervals.In order to allow a meaningful comparison there has to be a temporal reference point which is typically set by some kind of trigger (e.g., the onset of an external stimulation).There are two prominent applications for this successive trials scenario.Repeated Email address: thomas.kreuz@cnr.it(Thomas Kreuz).
presentation of the same stimulus addresses the reliability of individual neurons (Mainen and Sejnowski, 1995), while different stimuli are used to investigate neuronal coding and to find the features of the response that provide the optimal discrimination (e.g., Victor, 2005, for a more general introduction to neural coding cf.Quian Quiroga et al., 2002).These two applications are related since for a good clustering performance one needs both a pronounced discrimination between stimuli (high inter-stimulus spike train distances) and a high reliability for the same stimulus (low intra-stimulus spike train distances).
Electrophysiology and other modern recording techniques are developing fast.For both simultaneous population and successive trial recordings they often provide more data than available methods of spike train analysis can handle.There is a lack of algorithms able to identify multiple spike train patterns across different spatial scales and with a high temporal resolution.This is noticeable in both scenarios.In epilepsy, the analysis of the varying similarity patterns of simultaneously recorded ensembles of neurons can lead to a better understanding of the mechanisms of seizure generation, propagation, and termination (Truccolo et al., 2011;Bower et al., 2012).Similarly, the analysis of neuronal responses to successive presentations of time-dependent stimuli will help to understand the relevance of synchronous firing in neural coding (Miller and Wilson, 2008).Moreover, in population recordings it would be even more advantageous to be able to monitor spike train synchrony in realtime.This would be a necessary condition for a prospective epileptic seizure prediction algorithm (Mormann et al., 2007), but it could also be very useful for the rapid online decoding needed to control prosthetics (Hochberg et al., 2006;Sanchez, 2008).
In recent years, two such time-resolved measures have been proposed: The ISI-distance (Kreuz et al., 2007) and the SPIKE-distance (Kreuz et al., 2013) both rely on instantaneous estimates of spike train dissimilarity which makes it possible to track changes in instantaneous clustering, i.e., time-localized patterns of (dis)similarity among multiple spike trains.Additionally, both measures are parameter-free and time-scale independent.Furthermore, the SPIKE-distance also comes in a causal variant (Kreuz et al., 2013) which is defined such that the instantaneous values of dissimilarity are derived from past information only so that time-resolved spike train synchrony can be estimated in real-time.Both measures have already been widely used in various contexts (e.g., for the most recent measure, the SPIKE-distance: Papoutsi et al., 2013;Di Poppa and Gutkin, 2013;Sacré and Sepulchre, 2014).Another time-scale independent and time-resolved method is event synchronization (Quian Quiroga et al., 2002), a sophisticated coincidence detector which quantifies the level of synchrony from the number of quasi-simultaneous appearances of spikes.Originally, it was proposed and used in a bivariate context only.In this paper it is adapted to the time-resolved SPIKY-framework and extended to the multivariate case.
With all of these measures spike trains can be analyzed on different spatial and temporal scales, accordingly there are several levels of information extraction (Kreuz, 2012).In the most detailed representation one instantaneous value is obtained for each pair of spike train.The most condensed representation successive temporal and spatial averaging leads to one single distance value that describes the overall level of synchrony for a group of spike trains over a given time interval.In between these two extremes are spatial averages (dissimilarity profiles) and temporal averages (pairwise dissimilarity matrices).This variety of representations makes a straightforward implementation of the measures in one simple program/function unfeasible.Other important goals are high computational speed, efficient memory management, and applicability to large datasets.What is needed is an intuitive and interactive tool for analyzing spike train data which is able to overcome all of these challenges.
Here we address this need and present the graphical user interface SPIKY.Given a set of real or simulated spike train data (importable from many different formats), SPIKY calculates the measures of choice and allows the user to switch between many different visualizations such as measure profiles, pairwise dissimilarity matrices, or hierarchical cluster trees.SPIKY also includes the possibility of generating movies which are very useful in order to track the varying patterns of (dis)similarity.SPIKY has been optimized with respect to both computation speed (by using MEX-files, i.e.C-based Matlab executables) and memory demand (by taking advantage of the piecewise linear nature of the dissimilarity profiles).Finally, the SPIKY-package includes two complementary programs.The first program SPIKY loop is meant to be used for the analysis of a large number of datasets.The second program SPIKY loop surro is designed to evaluate the statistical significance of the results obtained for the original dataset by comparing them against the results of spike train surrogates generated from that dataset.
The remainder of this paper is organized as follows.In Section 2 we present the different measures available in SPIKY and provide some details about their implementation.These measures include the ISI-distance and the SPIKE-distance as well as the latter's realtime variant, and, introduced here, its forward variant (Section 2.1).For the SPIKE-distance we propose a correction of the edge-effect (spurious decrease to zero due to auxiliary spikes).In Section 2.2 we introduce the modified and extended variant of event synchronization.Some improvements realized in the new implementation of the measures are presented in Section 2.3.In the same Section we also compare the performance of this new implementation with the one of previously published source codes.In Section 3 an overview of the different ways to extract information is given.SPIKY, our graphical user interface for monitoring spike train synchrony, is presented in Section 4. In Section 4.1 we explain how to get access to the GUI and its documentation.Subsequently, in Section 4.2 we introduce the structure and the workflow of SPIKY, in particular we show how to input spike train data, how to change the layout of the figures and how to export results.The two complementary programs SPIKY loop and SPIKY loop surro are introduced in Sections 4.3 and 4.4, respectively.Finally, in Section 5 we summarize the methods and the program and present an outlook on future developments.

Measures and implementation
SPIKY implements four measures, one is multivariate and three are bivariate.The multivariate measure, included for comparison, is the standard Peri-Stimulus Time Histogram (PSTH) which measures the overall firing rate.In this Section we give an overview over the three bivariate measures.
The ISI-distance, the SPIKE-distance, and event synchronization all are time-resolved and rely on instantaneous values in the sense that in a first step the sequences of discrete spike times {t (1) i }, i = 1, ..., M 1 and {t (2) j }, j = 1, ..., M 2 (with M n denoting the number of spikes for spike train n with n = 1, 2) are transformed into dissimilarity profiles, I(t), S(t), and E(t), respectively.All three of these dissimilarity profiles are normalized between 0 and 1 and the same holds true for the distance values D I , D S and D E , which are defined as their temporal average.For all three of these bivariate measures there exists a straightforward extension to the case of more than two spike trains (spike train number N > 2), the averaged bivariate distance.This average over all pairs of spike trains commutes with the average over time, so it is possible to achieve the same kind of time-resolved visualization as in the bivariate case by first calculating the instantaneous average S a (t) (here for the SPIKE-distance) over all pairwise instantaneous values S mn (t), S mn (t). (1)

The ISI-and the SPIKE-distance
The dissimilarity profiles of the ISI-and the SPIKEdistance are based on three piecewise constant quantities which for each spike train n = 1, 2 are assigned to every time instant (see Fig. 1A).These are the time of the preceding spike the time of the following spike as well as the interspike interval The ambiguity regarding the definition of the very first and the very last interspike interval is resolved by placing for each spike train auxiliary leading spikes at time t = 0 and auxiliary trailing spikes at time t = T (but see also Section 2.1.2).

The ISI-distance
The ISI-distance, proposed as a bivariate measure in Kreuz et al. (2007) and extended to the multi spike train case in Kreuz et al. (2009) was the first spike train distance directly defined as the temporal average of an instantaneous dissimilarity profile.This profile is calculated as the absolute value of the instantaneous ratio between the interspike intervals x (1) ISI and x (2) ISI (see Fig. 1) according to: Since the ISI-values only change at the times of spikes, the dissimilarity profile is piecewise constant (with discontinuities at the spikes).The ISI-ratio equals 0 for identical interspike intervals in the two spike trains, and approaches 1 in intervals in which one spike train is much faster than the other.The ISI-distance is defined as the temporal average of this absolute ISI-ratio:

The SPIKE-distance
The SPIKE-distance (see Kreuz et al., 2011, for the original proposal and Kreuz et al., 2013;Kreuz, 2012 for the definite version presented here) is the centerpiece of SPIKY.In contrast to the ISI-distance, it considers the exact timing of the spikes.
The dissimilarity profile is calculated in two steps: First for each spike a spike time difference is calculated, then for each time instant the relevant spike time differences are selected, weighted, and normalized.Here 'relevant' means local; each time instant is surrounded by four corner spikes: the preceding spike from the first spike train t (1) P , the following spike from the first spike train t (1) F , the preceding spike from the second spike train t (2) P , and, finally, the following spike from the second spike train t (2) F .To each of these corner spikes one assigns the distance to the nearest spike in the other spike train, for example, for the previous spike of the first spike train: ∆t and analogously for t (1) P , and t (2) F (see Fig. 1).Subsequently, for each spike train separately, a locally weighted average is employed such that the differences for the closer spike dominate; for each spike train n = 1, 2 the weighting factors are the intervals to the previous and to the following spikes: x The local weighting for the spike time differences of the first spike train reads: S 1 (t) = ∆t (1) P (t)x (1) F (t)x (1) 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: This quantity sums the spike time differences for each spike train weighted 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.In order to get these ratios straight and to account for differences in firing rate, in a last step the two contributions from the two spike trains are locally weighted by their instantaneous interspike intervals.This leads to the definition of the dissimilarity profile: Again, the overall distance value is defined as the temporal average of the dissimilarity profile: Since the dissimilarity profile S(t) is obtained from a linear interpolation of piecewise constant quantities, it is piecewise linear (with potential discontinuities at the spikes).Both the dissimilarity profile S(t) and the SPIKEdistance D S are bounded in the interval [0, 1].The distance value D S = 0 is obtained for identical spike trains only.
Due to the finite recording time there is an ambiguity regarding the definitions of the initial distance to the preceding spike, the final distance to the following spike, as well as the very first and the very last interspike intervals.In previous implementations of the SPIKE-distance this ambiguity was resolved by adding to each spike train an auxiliary leading spike at time t = 0 and an auxiliary trailing spike at time t = T .This lead to spurious synchrony at the edges where by construction the dissimilarity profile reached the zero value.Here we partly correct this edge effect by incorporating all the information that is available.We describe the correction only for the beginning of the recording, an analogous procedure is applied at the end of the recording.
We count the auxiliary spikes as normal spikes which can be nearest neighbors to other spikes.But instead of calculating their spike time distance (which is always zero) we use the spike time difference of the first real spike.For the first interspike interval we know that it is at least the distance to the first spike t 1 −t 0 = t 1 but it could be longer.So to take the local firing rate into consideration we set where we use the length of the first known interspike interval t 2 − t 1 as a lower limit.In case t 1 is smaller than t 2 − t 1 we get at least a crude estimate of how much longer the first interspike interval could have been.

Realtime SPIKE-distance
In contrast to the dissimilarity profile S(t) of the regular SPIKE-distance, the dissimilarity profile S r (t) of the realtime SPIKE-distance 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 profile of the regular SPIKE-distance, 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: Since there are no following spikes available, there is no local weighting.There is no interspike interval either, so the normalization is achieved by dividing the average corner spike difference by twice the average time interval to the preceding spikes (Eq.8).This yields a causal indicator of local spike train dissimilarity: S r (t) = ∆t (1) (16)

Forward SPIKE-distance
The dissimilarity profile S f (t) of the forward SPIKEdistance is 'inverse' to the profile of the realtime SPIKEdistance.Instead of relying on past information only it relies on forward information only.It can be used in triggered temporal averaging in order to evaluate the (causal) effect of certain spikes or of specific stimuli features on future spiking.Again, for each time instant there are just two corner spikes and the potential nearest spikes in the other spike train are future spikes only.Thus, the spike time difference for the following spike of the first spike train reads: and accordingly for the following spike of the second spike train.In analogy to Eq. 16, an indicator of local spike train dissimilarity is obtained as follows: (18)

Event synchronization
Event synchronization (Quian Quiroga et al., 2002; see also Kreuz et al., 2007 andKreuz, 2011) works as a coincidence detector which quantifies the degree of synchrony from the relative number of quasi-simultaneous appearances of spikes.In Quian Quiroga et al. (2002) event synchronization was proposed together with the delay asymmetry, an antisymmetric measure which distinguishes between the two possible relative orders among these quasisimultaneous spikes.Here we are only interested in the level of synchrony, so we can use a simplified notation with just one coincidence counter c(t).This cumulative function counts the appearances of past spikes (i.e.before time instant t) in the first spike train within a certain coincidence window around past spikes in the second spike train: where There are two variants which differ in the way a coincidence is defined.The main variant is parameter-and scalefree, since the coincidence window τ ij , i.e., the time lag up to which two spikes t (1) i and t (2) j are considered to be synchronous, is adapted to the local spike rates (see Fig. 1B): (21) In contrast, the second variant uses a fixed coincidence window τ instead of the adaptive τ ij .
In order to be able to normalize event synchronization properly we need two more counters which keep track of the number of past spikes in the two spike trains.The cumulative spike number function for the first spike train can be written as and analogously for the second spike train.
The cumulative event synchronization function can now be obtained as the coincidence counter normalized by the average spike number in the two spike trains.We here follow Quian Quiroga et al. (2002) and normalize to the geometric mean of the two spike numbers: Since both the individual spike counters and the coincidence counter can only be incremented at the times of the spikes, all contributing functions are piecewise constant and thus so is the cumulative event synchronization.This is the property which makes event synchronization suitable for efficient use in our time-resolved SPIKY-framework.With this aim, we now switch from a time-based representation of Q(t) to an interval-based representation Q(k) in which we represent each interval k of the pooled spike train by just the one value which Q(t) attains for the whole interval.The form of Eq. 23 stays the same, only the t-dependence is replaced by a k-dependence.
Then, in order to get a time-resolved similarity profile of event synchronization, we define a local version of the Q(k)-function which only considers the changes in the three cumulative functions from the beginning to the end of a window centered around each interval k: and ∆m 1 (k) and ∆m 2 (k) defined accordingly.
It is straightforward to transform the E(k) values back into a piecewise constant dissimilarity profile E(t).While this 'piecewise constant' property is shared with the dissimilarity profile I(t) of the ISI-distance, there are also some differences.The temporal averaging is not performed by means of the integration used in Eq. 6, rather it is carried out on the finite set of interval-based values E(k), k = 1, ..., K according to The reason is that for event synchronization it is only important if spikes coincide or not; the lengths of the intervals in the pooled spike train do not matter.Another difference is that E is not parameter-free.The temporal resolution is governed by the window size ∆k.For the second variant of event synchronization there is also a second parameter, the fixed coincidence window τ .In order to find reasonable results it is necessary to adapt these parameters to the spike train data.
Finally, in contrast to I(t) -and also S(t) -the timeresolved event synchronization is oriented as a measure of similarity, i.e., it yields E(t) = 1 if and only if all spikes within the window coincide (or if there are no spikes at all).This orientation is the same as the one for the Peri-Stimulus Time Histogram so it makes sense to compare these two measures (Fig. 2).Fig. 2A demonstrates the similarities and dissimilarities between the temporal profiles of the two measures on a rather general example.Fig. 2B-D shows that event synchronization, in contrast to the PSTH, is a true measure of spike train synchrony (see also Kreuz et al., 2011).Since the PSTH is invariant to shuffling spikes among the spike trains, it yields the same value regardless of how spikes are distributed among the different spike trains.

Comparison with other implementations
The very first all-in-one implementation of the ISI-and the SPIKE-distance was published online at the time of the publication of Kreuz et al. (2013).In between this first In the first half within the noisy background there are 4 regularly spaced spiking events with increasing jitter.The second half consists of 10 spiking events with decreasing jitter but now without any noisy background.In the noisy first half PSTH and event synchronization exhibit very similar profiles.The fact that the firing events become more distinct in the second half is indicated by event synchronization as a gradual increase to synchronization.In the PSTH the peaks become more and more narrow.The dissimilarity profile of the SPIKE-distance for the same example can be found in Fig. 2B of Kreuz et al. (2013) stricted to the dissimilarity profile and its temporal average (the overall dissimilarity).In contrast, SPIKY also allows the user to interactively access all the other different ways to extract information that will be introduced in Section 3.All of these codes for calculating the ISI-and the SPIKEdistance rely on equidistantly sampled dissimilarity profiles, and the same holds true for older codes of event synchronization.Typically the precision is set to the sampling interval of the neuronal recording.Since the dissimilarity profiles have to be calculated and stored for each pair of spike trains, one obtains, for each measure, a matrix of order 'number of sampled time instants' × 'number of spike train pairs' (i.e., #(t s ) × N (N − 1)/2).For small sampling intervals and large numbers of spike trains this leads to memory problems.
In SPIKY we use an optimized and more memoryefficient way of storing the results.We make use of the fact that the dissimilarity profiles I(t) of the ISI-distance and E(t) of event synchronization are piecewise constant and the dissimilarity profile S(t) of the SPIKE-distance is piecewise linear.Each constant/linear interval runs from one spike of the pooled spike train to the next.Thus, for each such interval (and for each pair of spike trains) we have to store only one value for the ISI-distance and for event synchronization and two values for the SPIKEdistance, one at the beginning and one at the end of the interval (see Fig. 3).The memory gain is proportional to the number of sample points per interspike interval in the pooled spike train and is typically much larger than one.
The dissimilarity profiles exhibit instantaneous jumps at the times of the spikes since this is where the lengths of the interspike intervals and the identity of the previous and the following spikes change abruptly or where a new coincidence is counted.For sampled dissimilarity profiles one has to 'cut the corners' of these instantaneous jumps.This leads to an estimation error which increases with the sampling interval.In contrast, within the new implementation each spike marks both the end of the previous and the beginning of the next interval and it becomes possible to store two dissimilarity values for these points.This way the integration from one spike of the pooled spike train to the next can be performed over the full interval.Thus, besides being far more memory-efficient, the new implementation also computes the exact distance values without any spurious dependence on the sampling interval.
The third effect of the new implementation is a considerable speed-up.To show this (and to provide relative computational costs for the different measures) we here calculate the speed gain achieved by going from the old equidistantly sampled implementation to the new minimally sampled implementation of the ISI-and the SPIKE-distance.
As benchmark we use the comprehensive performance comparison carried out by Rusu and Florian (2014).In this test the authors compared the performance of their newly proposed modulus-metric with the performance of previously proposed spike train distances including the ISIand the SPIKE-distance.Like them we used two random spike trains with different numbers of spikes.However, since we were also interested in applications to larger datasets, we extended the maximum number of spikes from 500 to 10000 spikes.The firing rate is kept constant, hence the duration of the trial increases with the number of spikes.For a fair comparison we implemented all algorithms in C++ and ran them on an Intel i7-4700MQ CPU @ 2.4 GHz.All speed gains were averaged over 10, 000 trials (Fig. 4).
In a first step we replicated the results of Rusu and Florian (2014) who had calculated the ISI-and SPIKEdistances using dissimilarity profiles that were equidistantly sampled with a fixed time step of dt = 1 ms.Rerunning the same algorithms on our computer, we could reproduce their results (for less than 500 spikes) within the same order of magnitude.Subsequently, we measured the running times using minimally sampled dissimilarity profiles and found that for both the ISI-and the SPIKEdistance the new implementation was considerably faster than the old implementation.As expected, the speed gain depends critically on the sampling rate used for the old implementation and is particularly large for densely sampled data (Fig. 4).
Apart from minimal sampling, SPIKY includes another improvement which strongly increases the overall performance.The source codes published along with Kreuz et al. (2013) were entirely written in Matlab.In contrast, the new SPIKY-implementation uses C-based Matlab executables (MEX)-files for the most time-consuming parts and this leads to another enormous performance boost (typically by a factor between 3 and 30).

Levels of information extraction
The ISI-and the SPIKE-distance combine a variety of properties that make them well suited for the application to real data.In particular, they are conceptually simple, computationally efficient, and easy to visualize in a timeresolved manner.By taking into account only the preced- ing and the following spike in each spike train, these distances rely on local information only.They are also timescale-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 same holds true for event synchronization which is also both time-resolved and time-scale adaptive.Moreover, the sensitivity to spike timing and the instantaneous reliability achieved by the SPIKE-distance opens up many new possibilities in multi-neuron spike train analysis (Kreuz et al., 2013).These build upon the fact that there are several ways to extract information all of which we describe in the following.As an illustration we use the detailed analysis of an artificially generated spike train dataset with the SPIKE-distance (see upper part of Fig. 5A for the rasterplot).

Full matrix and cross sections
The starting point is the most detailed representation in which one instantaneous value is obtained for each pair of spike trains (see Eq. 12).This representation could be viewed as a movie of a symmetric pairwise dissimilarity matrix in which each frame corresponds to one time instant (an example can be found in the supplementary material of Kreuz et al., 2013).For a movie of finite length the time axis necessarily has to be sampled but in principle this most detailed representation consists of an infinite number of values.However, since all dissimilarity profiles are piecewise linear there is a lot of redundancy.Using the most compact and memory-efficient representation one can store all pairwise dissimilarity profiles in a matrix of size 'number of interspike intervals in the pooled spike train' × 'number of spike train pairs' (×2 for the SPIKE-distance, see Section 2.3).From this two-dimensional matrix it is possible to extract both kinds of cross sections.By selecting a pair of spike trains, one obtains the bivariate dissimilarity profile S(t) for this pair of spike trains.Selecting a time instant t s (and using linear interpolation for time instants in between spikes) yields an instantaneous matrix of pairwise spike train dissimilarities S mn (t s ) (see Fig. 5B).This matrix can be used to divide the spike trains into instantaneous clusters, that is, groups of spike trains with low intra-group and high inter-group dissimilarity.

Spatial and temporal averaging
Another way to reduce the information of the dissimilarity matrix is averaging.There are two possibilities that commute: the spatial average over spike train pairs and the temporal average.Since the spatial average over spike train pairs can be done locally it yields a dissimilarity profile for the whole population.Examples for averages over four different spike train groups as well as over all spike trains are shown in the lower subplot of Fig. 5A.Temporal averaging over certain intervals on the other hand leads to a bi-variate distance matrix (see Fig. 5C and D for examples of non-continuous and continuous intervals).In real data, these temporal intervals could be chosen to correspond to different external conditions such as normal vs. pathological, asleep vs. awake, target vs. non-target stimulus, or presence/absence of a certain channel blocker.
A combination of temporal and spatial averaging can be seen in Fig. 5F.This dissimilarity matrix is obtained from the overall temporal average shown in Fig. 5D by (spatially) averaging over the 16 submatrices and thus depicts the pairwise spike train group dissimilarity (4 × 4 instead of 20 × 20).Fig. 5G shows the respective dendrogram.In applications to real data, these groups could be different neuronal populations or responses to different stimuli, depending on whether the spike trains were recorded simultaneously or successively.Finally, successive application of spatial average over all spike train pairs and temporal average over the whole interval results in just one distance value that describes the overall level of synchrony for the entire dataset.In Fig. 5A this value is stated in the upper right of the lower subplot.

Triggered averaging
The fact that there are no limits to the temporal resolution allows further analyses such as internally or externally triggered temporal averaging.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.The trigger times can either be obtained from external influences (such as the occurrence of certain features in a stimulus) or from internal conditions (such as the spike times of a certain spike train).External triggering is a standard tool to address questions of neural coding, for example, it can be used to evaluate the influence of localized stimulus features on the reliability of neurons under repeated stimulation.In multi-neuron data, internal triggering might help to uncover the connectivity in neural networks or to detect converging or diverging patterns of firing propagation.An example is shown in Fig. 5E.Here, neurons 2, 9, 14, and 19 follow the 7-th neuron during the 2nd and the 4th 500-ms subinterval.This can be revealed by triggering on the spike times of neuron 7 during these subintervals.The difference between the dissimilarity matrix of the full interval in Fig. 5D and the triggered dissimilarity matrix of 5E can be best seen by comparing the respective dendrograms in Fig. 5H and 5I.

SPIKY
SPIKY is a graphical user interface for monitoring synchrony between artificially simulated or experimentally recorded neuronal spike trains.It contains implementations of the ISI-distance, the SPIKE-distance, and event synchronization.For all three of these measures it allows interactive access to the different representations described in Section 3.Moreover, SPIKY includes the forward variant of the SPIKE-distance as well as a 'simulation' of the realtime SPIKE-distance.It calculates the values the realtime version would give but it does not work in real time (see also Section 5.2).All source codes are written in Matlab (MathWorks Inc, Natick, MA, USA) with the most time-consuming loops coded in MEX-files3 .Consequently, SPIKY is not stand-alone but requires Matlab to run.

Access to SPIKY and how to get started
SPIKY is distributed under the BSD licence (Copyright (c) 2014, Thomas Kreuz, Nebojsa Bozanic.All rights reserved.).A zip-package containing all the necessary files can be accessed for free on the download page4 .This package also contains a folder with documentation (such as a FAQfile and an introduction to all individual elements and all individual files of SPIKY).Further information and many demonstrations (both images and movies) can be found on the download page and on the SPIKY Facebook-page5 .Both of these pages are used to announce updates and distribute the latest information about new features.They also provide an opportunity to give feedback and ask questions.
The Facebook-page includes various screen recordings with voice-over in which the user is guided step by step through some of the most important features of SPIKY.All of these movies can also be viewed on the SPIKY Youtubechannel6 .
After downloading SPIKY the user has to first extract the zip-package which leaves all files in one folder named 'SPIKY'.If the system has a suitable MEX-compiler installed, the MEX-files can be compiled from within this folder by running the m-file 'SPIKY compile MEX'.The program is started with the m-file 'SPIKY'.
When SPIKY is running, the user can find quick information about the individual elements of the graphical user interface by activating the 'Hints'-checkbox in the 'Options'-Menu.Hovering with the mouse cursor above elements of interest will then show short hints.An overview of all the information contained in the hints can be found in the documentation file 'SPIKY-Elements'.Furthermore, at each step the suggested element for the next user action is highlighted by a bold font.
To get the user started quickly, SPIKY provides a few example datasets from previous publications.The most useful example is the entry 'Clustering' in the 'Selection: Data'listbox.This dataset has already been used in several fig- A typical SPIKY-session begins at the top with 'Get data', then goes clockwise and ends on the left with 'Get results'.For clarity we omitted two possible which can be performed from multiple positions on the workflow circle: From any point it is possible to reset to the very beginning ('reset' leads to the symbol ⊲), and from any later point it is possible to jump back before the measure calculation ('reset with the same data' leads to the symbol ▽).
ures as well as in the supplementary movie of Kreuz et al. (2013).The best way to get acquainted with SPIKY is to advance from panel to panel by pressing the highlighted button.When the end is reached the user can reset and run the same example again while changing some parameters in order to see the consequences.Note that it is not necessary to set all the parameters each time when SPIKY is started.Rather, it is possible to use the file 'SPIKY f user interface' to set and modify the spike train data as well as the parameters (again, the dataset 'Clustering' provides an example).

Structure and workflow of SPIKY
Overall, SPIKY has a rather linear workflow, however, it is much more interactive than previous implementations of the measures and there are many potential shortcuts and loops along the way.As can be seen in the SPIKY-flowchart in Fig. 6, the general flow is clearly directed from the input of spike train data to the output of results.So the first step the user has to do is to give SPIKY spike train data (i.e.sequences of spike times) to work with.There are four possible ways to do this: one can make use of predefined examples, load data from a file or from the Matlab workspace, detect discrete events from continuous data, or employ the spike train generator (see Section 4.2.1 for more details on the SPIKY input).
Once the full dataset is available, modification is still possible.The user can restrict the analysis to a specific subset, e.g., select a smaller time window and/or a subset of spike trains.It is also possible to impose some external structure on the rasterplot (spike trains vs time).For that, SPIKY allows the definition of two types of time markers (e.g.specific events such as seizure onset and offset in epilepsy, trigger onset during stimulation etc.) and two types of spike train separators (e.g.neurons from the left vs. neurons from the right hemisphere).The user can also define spike train groups.Depending on the setup these could be spike trains recorded in different brain regions or upon presentation of different stimuli.Fig. 7 shows an example of a raster plot with annotations marking all these different elements.
After updating all of these data parameters, the next step is to select the measures to be calculated.Options include the Peri-Stimulus Time Histogram as well as all measures described in detail in Section 2. In the same step the user can select successive frames for a temporal analysis of spike train patterns.These can be individual time instants for cross sections, temporal intervals for selective averages and sequences of time instants for triggered averages (see Section 3).Now the actual calculation of the measures takes place.For reasonably sized datasets this should take at most a few seconds.Very large datasets (typically datasets containing hundreds of spike trains and/or hundred thousands of spikes) are divided in smaller subintervals and the calculation will be performed in a loop which might take longer.It is from this point on that SPIKY becomes truly interactive.Now the user can switch between different representations of the results (such as dissimilarity profiles, dissimilarity matrices and dendrograms).Different matrices and dendrograms can be compared in the same figure or they can be viewed in sequence.
The presentation can be restricted to smaller time windows and/or subsets of spike trains, and temporal and spatial averaging (for example moving average and average over spike train groups) can be performed.Spike trains or spike train groups can also be sorted according to the num-ber of spikes (either within the whole spike train or within a certain interval) or according to the spike latency with respect to a specific time instant.The order can even depend on the result of the clustering analysis, i.e., spike trains belonging to the same cluster will appear next to each other.
Furthermore, it is possible to add further figure elements such as spike number histograms, overall averages, or dissimilarity profiles for individual spike train groups.At this stage the user can also retrospectively change the appearance of all the individual elements of the figure (see Section 4.2.2 for more details).Finally, SPIKY allows to extract both data and results to the Matlab workspace for further analysis, and it is also possible to save individual figures as postscript-file or a sequence of images as an 'avi'-movie (see Section 4.2.3 for more details on the SPIKY output).

Input
There are four different ways to input spike train data into SPIKY.
The first option is to select one of the predefined examples which are generated using Matlab-code.Initially these are the examples used in Kreuz et al. (2013) but one can also define new examples.
The second option is to load spike train data either from the Matlab workspace or from a file.Two different file formats are allowed, '.mat' and '.txt' (ASCII) files.For the mat-files SPIKY currently allows three different kinds of input formats (further formats can be added on demand): -cell arrays (ca) with just the spike times.This is the preferred format used by SPIKY since it is most memory efficient.The two other formats will internally be converted into this format; -regular matrices with each row being a spike train and zero padding (zp) in case the spike numbers are different; -matrices representing time bins where each zero/one (01) indicates the absence/presence of a spike.
In case of a mat-file or of the workspace, SPIKY looks for a variable called 'spikes'; if it cannot find it, the user has the chance to select the variable name (or field name) which contains the spikes via an input mask which provides a hierarchical structure tree of all the variables and fields contained in the mat-file or in the workspace.
In the text format spike times should be written as a matrix with each row being one spike train.The SPIKYpackage contains one example file for all four formats ('testdata ca.mat ', 'testdata zp.mat', 'testdata 01.mat' and 'testdata.txt').
The third option is to use the event detector in order to detect discrete events in continuous data.There are many different possibilities of defining an event.A variety of standard events (such as local maxima and minima and threshold crossings) with a number of parameters are already included.
The fourth option is to create new spike train data via the spike train generator.After setting some defining variables (number of spike trains, start and end time, sampling in-terval) the user can build spike trains from predefined spike train patterns (such as periodic, splay, uniform or Poisson) and/or by manually adding, shifting and deleting individual spikes or groups of spikes.

Figure-Layout
SPIKY is designed so that it can directly generate figures suitable for publication.With this aim the user is given control over the appearance of every individual element (e.g.fonts, lines etc.) in each type of figure.There are two ways to determine essential properties such as color, font size or line width.It is possible to change elements in the active figure while the program is already running.To do so the user has to simply click the right mouse button on the element to be changed.A context menu will appear which lets the user edit the properties of individual elements or of all elements of a certain type.Conveniently, one can also use the file 'SPIKY f user interface' to define the standard values for all the parameters that describe the principal layout of the figure .If a figure contains more than one subplot (besides the subplot containing the spike rasterplot and/or the dissimilarity profiles, these are typically subplots with dissimilarity matrices and dendrograms), it is possible to change their position and size.To do so, the user can move the cursor to the respective axis (either just left or just below the subplot) and click the right mouse button.Now one can edit all position variables by hand or change the x-position, the y-position, the width and the height individually.In case there are several dissimilarity matrices/dendrograms this can be done either for an individual matrix/dendrogram or for all of them at once.

Output
From within SPIKY it is possible to extract the spike trains and the results of the analysis (measure profiles, matrices, dendrograms) to the Matlab workspace for further processing.After clicking the right mouse button on the element whose data one wishes to extract, results will be stored in variables such as 'SPIKY spikes', 'SPIKY profile X 1', 'SPIKY profile Y 1', 'SPIKY profilename 1', 'SPIKY matrix 1' and 'SPIKY matrix name 1'.In addition, the results obtained during an analysis will automatically be stored in the output structure 'SPIKY results' which will have one field for each measure selected.Depending on the parameter selection within SPIKY, for each measure the structure can contain the following subfields that correspond to the different representations identified in Section 3: -SPIKY results.<Measure>.name:Name of selected measures -SPIKY results.<Measure>.distance:Level of dissimilarity over all spike trains and the whole interval.This is just one value, obtained by averaging over both spike trains and time -SPIKY results.<Measure>.matrix:Pairwise distance matrices, obtained by averaging over time -SPIKY results.<Measure>.time:Time-values of overall dissimilarity profile -SPIKY results.<Measure>.profile:Overall dissimilarity profile obtained by averaging over spike train pairs Note that the dissimilarity profiles are not equidistantly sampled.Rather they are stored as memory-efficiently as possible which means just one value for each interval of the pooled spike train for the ISI-and two values for the SPIKE-distance.Since this format can be more difficult to process, SPIKY includes three functions: 'SPIKY f selective averaging' for computing the selective average over time intervals, 'SPIKY f triggered averaging' for calculating the triggered average over time instants, and 'SPIKY f average pi' for averaging over many dissimilarity profiles.Furthermore, for the ISI-distance the function 'SPIKY f pico' can be used to obtain the average value as well as the x-and y-vectors for plotting.
Besides the standard way to work with Matlab-figures SPIKY also offers the opportunity to save each figure as a postscript-file.Finally, it is possible to save a sequence of images as an 'avi'-movie.

GUI vs. loop
SPIKY was mainly designed to facilitate the detailed analysis of one dataset.It enables the user to switch between different representations (see Section 3) and to zoom in on both spatial and temporal features of interest.However, SPIKY is not very suitable for the analysis of many different datasets when e.g. the statistics of a certain quantity such as an average over a specific time interval should be evaluated over all available datasets (e.g. over all trials of a stimulus setup or for recordings of all subjects etc.) in some kind of loop.For these purposes the SPIKY-package contains a program called SPIKY loop which is complementary to SPIKY.It is not a graphical user interface but it should be simple enough (and plenty of examples are provided) to allow everyone to run the same kind of analysis for many different datasets and to evaluate and compare their 'SPIKY results'.SPIKY loop uses the full functionality of SPIKY such as access to time instants, selective and triggered averages as well as averages over spike train groups.
So by combining these two programs it is possible to first use SPIKY for a rather exploratory but detailed analysis of a limited number of individual datasets and then use SPIKY loop and its output structure 'SPIKY loop results' to verify whether any effect discovered on the example dataset is consistently present within all of the datasets.
Both SPIKY and SPIKY loop require the storage of matrices of the order 'number of interspike intervals in the pooled spike train' × 'number of spike train pairs'.For very large datasets with many spike trains and/or spikes this can lead to memory problems.We addressed this issue by making the calculation sequential, i.e., by cutting the recording interval into smaller segments, and performing the averaging over all pairs of spike trains for each segment separately.In the end the dissimilarity profiles for the different segments (already averaged over pairs of spike trains) are concatenated, and its temporal average yields the distance value for the whole recording interval.During this sequential calculation SPIKY uses a waitbar (which displays the % completed) to continuously inform the user about the progress.

Spike train surrogates and significance
An important question that has not yet been asked is the one of statistical significance.Given a certain value of the SPIKE-distance how can one judge whether it reflects a significant decrease or increase in spike train synchrony and does not just lie within the range of values obtained for random fluctuations.One way to address this question is the use of spike train surrogates (Kass et al., 2005;Gruen, 2009;Louis et al., 2010).The idea is to compare the results obtained for the original dataset versus the results obtained for spike train surrogates generated from that dataset.If the value obtained for the original lies outside the range of values for the surrogates this value can be assumed to be significant to a level defined by the number of surrogates used (e.g.α = 0.05 for 19 surrogates or α = 0.001 for 999 surrogates).
The SPIKY-package contains a program 'Spiky loopsurro' which was designed to look at significance.So far it includes four different types of spike train surrogates.They differ in the properties that are preserved and maintain either the individual spike numbers (obtained by shuffling the spikes), the individual interspike interval distribution (obtained by shuffling the interspike intervals), the pooled spike train (obtained by shuffling spikes among the spike trains) or the Peri-Stimulus Time Histogram (PSTH).In the last case the spike train surrogates are obtained by means of inverse transform sampling, i.e., by resampling from the PSTH using its cumulative distribution function (CDF) (Ross, 1997).Other ways to calculate rate functions (e.g. based on kernels with different bandwidth) will be added in future releases.

Summary
In this paper we presented SPIKY, a graphical user interface which facilitates the application of methods of spike train analysis to both simulated and real datasets.Apart from the standard Peri-Stimulus Time Histogram, SPIKY contains three parameter-free and time-resolved measures (see Section 3).These measures are complementary to each other since each one addresses a different specific aspect of spike train similarity.While the ISI-distance quantifies local dissimilarities based on covariances of the neurons fir-ing rate profiles, both the SPIKE-distance and event synchronization capture the relative timing of local spikes7 .However, whereas the SPIKE-distance weights and normalizes the differences between nearest neighbor spikes, event synchronization acts as a binary coincidence detector, i.e. there is a cutoff at the (adaptive or fixed) time lag relative to which two neighboring spikes are either considered coincident or not and all detailed information both within or outside this coincidence window is discarded.
All of these measures yield instantaneous values for each pair of spike trains, and thus there are many different possible representations of the results (see Section 3).Often the most informative representation might depend on the amount and type of spike train data and SPIKY can be used to reveal it via some explorative and interactive analysis.SPIKY also allows to alter a given dataset before or after the actual analysis, e.g., to interactively select subintervals or spike train subsets, to define time markers and spike train separators, and to divide the dataset into different spike train groups.
In addition to the main GUI designed for the detailed analysis of one dataset, the SPIKY-package also includes two complementary programs.While SPIKY loop aims at the grand average analysis of large numbers of datasets, SPIKY loop surro allows the estimation of significance levels.In all of these programs we use MEX-files, i.e.C-based Matlab executables for the more time-consuming parts and we exploit the piecewise linearity of the dissimilarity profiles, thereby guaranteeing high computation speed and memory efficiency.

Outlook
One of the measures included in SPIKY is the realtime SPIKE-distance.The present algorithm calculates its instantaneous dissimilarity values by making use of past information only but it does so in a speed-optimized and parallelized manner which would not be compatible with an actual realtime-implementation.However, such a realtimeimplementation should actually be simpler and less problematic even for large numbers of neurons since the only information that would have to be stored at each time instant are the time stamps of the latest spikes of each spike train and their nearest neighbors in the other spike trains.
SPIKY was primarily designed to analyze neuronal spike trains.But in principle it is also applicable to any other kind of discrete data which comes in the form of sequences of time stamps (such as times of bouncing basketballs or time-coded social interactions in a psychological test).Fur-thermore, in Kreuz et al. (2013) the SPIKE-distance was already applied not only to discrete data but also to an example of continuous data (EEG).To clear the way for such applications SPIKY includes an event detector which allows to bridge the gap between continuous data and our methods for discrete data.
We chose to write SPIKY in Matlab due to its popularity in the neuroscience community, ease of use, and the engaging visualization capabilities of its GUI design.Because of its high level parallelization, Matlab provides a powerful tool for processing vectorized data, but it also includes a well-developed MEX-interface for integrating C functions for performance enhancements.C functions do not only lead to an increase in performance, they can also easily be incorporated into other programming platforms.Indeed, we have already ported some of its functionality to Python as part of the PySpike 8 module, an open-source project hosted on Github.As in SPIKY, in PySpike the computation intensive functions are implemented as C backends.
Another potential direction would be the extension of the methods used here from the quantification of synchrony within one population of two or more neurons to the quantification of synchrony between two (or more) neuronal populations.Similar extensions have been done for two other spike train distances (Aronov, 2003;Houghton and Sen, 2008) but both of these methods depend on not one but two parameters.So one particular challenge for a potential extension of our methods would be to make them parameter-free while maintaining their high temporal resolution.
There are still many open questions regarding neuronal synchrony.Among these questions are its role in dynamical diseases like epilepsy and its relevance for neural coding.While a thorough discussion of these issues is beyond the scope of this paper, we argue that SPIKY will be a very useful tool for investigation.If epileptic seizures and/or time-dependent stimulations lead to changes in spike train synchrony or spike train clustering, SPIKY should be able to detect them.

Fig. 2 .
Fig.2.Comparison of PSTH and event synchronization (first, adaptive variant with window size ∆k = 1).A. multivariate example with 50 spike trains.In the first half within the noisy background there are 4 regularly spaced spiking events with increasing jitter.The second half consists of 10 spiking events with decreasing jitter but now without any noisy background.In the noisy first half PSTH and event synchronization exhibit very similar profiles.The fact that the firing events become more distinct in the second half is indicated by event synchronization as a gradual increase to synchronization.In the PSTH the peaks become more and more narrow.The dissimilarity profile of the SPIKE-distance for the same example can be found in Fig.2BofKreuz et al. (2013).B-D.By construction the pooled spike train of these examples is identical consisting of 10 evenly spaced bursts.The only difference is the distribution of the spikes among the individual spike trains which varies from high via intermediate to low synchrony.Whereas the PSTH is the same for all three examples, event synchronization correctly indicates the decrease in synchrony.B. High reliability.Each spike train contains one spike per firing event.C. Random distribution of spikes among spike trains.D. Synfire chain of bursts.

Fig. 4 .
Fig. 4. Speed gain achieved by the new (minimally sampled) implementation of A. ISI-distance and B. SPIKE-distance with respect to the old (equidistantly sampled) implementation in dependence on the number of spikes and on the time step used by the old implementation.

Fig. 5 .
Fig. 5.The different levels of information extraction for the SPIKE-distance.A. Top: Spike rasterplot of 20 artificially generated spike trains divided in 4 spike train groups of 5 spike trains each.The clustering behavior changes every 500 ms.Bottom: Dissimilarity profiles of the SPIKE-distance for the four spike train groups (thin color-coded lines) and for all spike trains (thick black line).The overall dissimilarity is defined as the temporal average of the dissimilarity profile of all spike trains (0.192 in this case) and is marked by a dashed horizontal line.The green lines and symbols on top of the rasterplot mark temporal instants and intervals results of which are detailed in the subplots below.B-E.Matrices of pairwise instantaneous dissimilarity values for a single time instant (mark 1, subplot B), for two selective averages (over nonconsecutive intervals in mark 2, subplot C, and over the whole interval in mark 3, subplot D) and for a triggered average (mark 4, subplot E).F-G.For the overall average (mark 3) we also show the matrices of overall pairwise instantaneous dissimilarity values for the 4 spike train groups (F) and the corresponding dendrogram (G).H-I.Dendrogram of spike train matrices in D and E. Note that in contrast to the overall average (mark 3, subplot H) the triggered average (mark 4, subplot I) captures the local similarity between 5 of the spike trains.

Fig. 6 .
Fig.6.SPIKY-logo (center) and flowchart describing the workflow of SPIKY from the input of spike train data to the output of results.A typical SPIKY-session begins at the top with 'Get data', then goes clockwise and ends on the left with 'Get results'.For clarity we omitted two possible which can be performed from multiple positions on the workflow circle: From any point it is possible to reset to the very beginning ('reset' leads to the symbol ⊲), and from any later point it is possible to jump back before the measure calculation ('reset with the same data' leads to the symbol ▽).

Fig. 7 .
Fig. 7. Annotated screenshot from a movie.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 (denoted by < • > G ). C. Corresponding dendrograms.
. B-D.By construction the pooled spike train of these examples is identical consisting of 10 evenly spaced bursts.The only difference is the distribution of the spikes among the individual spike trains which varies from high via intermediate to low synchrony.Whereas the PSTH is the same for all three examples, event synchronization correctly indicates the decrease in synchrony.B. High reliability.Each spike train contains one spike per firing event.C. Random distribution of spikes among spike trains.D. Synfire chain of bursts.
(Thibeault et al., 2014)ation of both ISI-and SPIKE-distance courtesy of Rǎzvan Florian and hosted on GitHub 2 .The SPIKE-distance was also implemented in the commercially distributed HRLAnalysis T M software suite(Thibeault et al., 2014)designed for the analysis of large-scale spiking neural data.Note that all of these implementations are re-1 http://jeremy.fix.free.fr/Softwares/spike.html 2 https://github.com/modulus-metric/spike-train-metricsaTimeFig.3. Ilustration of the memory efficient storage of the dissimilarity profile for both the ISI-and the SPIKE-distance.