Measures of spike train synchrony for data with multiple time-scales
Eero Satuvuori, Mario Mulansky, Nebojsa Bozanic, Irene Malvestio, Fleur Zeldenrust, Kerstin Lenk, Thomas Kreuz
MMeasures of spike train synchrony for data with multiple time scales
Eero Satuvuori , Mario Mulansky , Nebojsa Bozanic , Irene Malvestio , Fleur Zeldenrust , Kerstin Lenk , ThomasKreuz Abstract
Background:
Measures of spike train synchrony are widely used in both experimental and computational neuroscience. Time-scale independent and parameter-free measures, such as the ISI-distance, the SPIKE-distance and SPIKE-synchronization, arepreferable to time scale parametric measures, since by adapting to the local firing rate they take into account all the time scales of agiven dataset.
New Method:
In data containing multiple time scales (e.g. regular spiking and bursts) one is typically less interested in thesmallest time scales and a more adaptive approach is needed. Here we propose the A-ISI-distance, the A-SPIKE-distance andA-SPIKE-synchronization, which generalize the original measures by considering the local relative to the global time scales. Forthe A-SPIKE-distance we also introduce a rate-independent extension called the RIA-SPIKE-distance, which focuses specificallyon spike timing.
Results:
The adaptive generalizations A-ISI-distance and A-SPIKE-distance allow to disregard spike time di ff erences that are notrelevant on a more global scale. A-SPIKE-synchronization does not any longer demand an unreasonably high accuracy for spikedoublets and coinciding bursts. Finally, the RIA-SPIKE-distance proves to be independent of rate ratios between spike trains. Comparison with Existing Methods:
We find that compared to the original versions the A-ISI-distance and the A-SPIKE-distanceyield improvements for spike trains containing di ff erent time scales without exhibiting any unwanted side e ff ects in other examples.A-SPIKE-synchronization matches spikes more e ffi ciently than SPIKE-Synchronization. Conclusions:
With these proposals we have completed the picture, since we now provide adaptive generalized measures thatare sensitive to firing rate only (A-ISI-distance), to timing only (ARI-SPIKE-distance), and to both at the same time (A-SPIKE-distance).
Keywords:
Data Analysis, Spike train distances, ISI-distance, SPIKE-distance, SPIKE-synchronization, Burst, Time-Scale
1. Introduction
In neuroscience the neuronal action potential and its com-plex molecular behavior (Bear et al., 2007) is often reducedto time-discrete events called spikes . Due to the all-or-nothingparadigm of neurons together with the long silent periods, the
Email addresses: [email protected] (Eero Satuvuori), [email protected] (Mario Mulansky), [email protected] (Nebojsa Bozanic), [email protected] (Irene Malvestio), [email protected] (Fleur Zeldenrust), [email protected] (Kerstin Lenk), [email protected] (Thomas Kreuz) Institute for Complex Systems, CNR, Sesto Fiorentino, Italy Department of Physics and Astronomy, University of Florence, SestoFiorentino, Italy MOVE Research Institute, Department of Human Movement Sciences,Vrije Universiteit Amsterdam, The Netherlands Department of Information and Communication Technologies, UniversitatPompeu Fabra, Barcelona, Spain Donders Institute for Brain Cognition and Behavior, Radboud Universiteit,Nijmegen, The Netherlands BioMediTech, Tampere University of Technology, Tampere, Finland DFG-Center for Regenerative Therapies Dresden, Technische UniversittDresden, Dresden, Germany time stamps of the spike events are considered to be an accurateenough description of the neuronal membrane potential (QuianQuiroga and Panzeri, 2013). These sequences of consecutivespikes are called spike trains . While spike trains do not di-rectly provide information about the connections between neu-rons, some form of link between two neurons is often inferredby the similarity of their spike trains. A spike train distancedoes not take into account the specific type of linkage, but sim-ply quantifies how (dis)similar the two spike trains are. Thismakes spike train distances universal and as such they can beapplied to all systems that can be reduced to point processes.In addition to the obvious neuroscience applications, they havealready been used to study inter-personal coordination (Rabi-nowitch and Knafo-Noam, 2015) and social cognition (Zapata-Fonseca et al., 2016) among many other fields.Over the years many di ff erent measures have been devel-oped in order to quantify similarities between two or more spiketrains (see Victor (2015), Naud et al. (2011) and Kreuz (2011)for an overview). The two most known time scale parametricmeasures, the Victor-Purpura (Victor and Purpura, 1996) andthe van Rossum distance (van Rossum, 2001), describe spiketrain (dis)similarity based on user-defined time scales to whichthe measures are mainly sensitive. One drawback of these mea- Preprint submitted to Elsevier May 31, 2017 a r X i v : . [ phy s i c s . d a t a - a n ] M a y ures is the fixed time scale, since it sets a boundary betweenrate and time coding for the whole recording. However, for realdata which typically contain many time scales (such as regularspiking and bursts), this is di ffi cult to detect with a measure thatis mainly sensitive to only one of them (Chicharro et al., 2011).The problem of having to choose one time scale has beeneliminated in the three time-resolved and time scale indepen-dent measures ISI-distance (Kreuz et al., 2007, 2009), SPIKE-distance (Kreuz et al., 2011, 2013) and SPIKE-synchronization(Kreuz et al., 2015). The ISI-distance (Kreuz et al., 2007) isa measure of instantaneous rate dissimilarity. It uses the in-terspike intervals (ISIs) to estimate the local firing rate of spiketrains and quantifies their di ff erences in a time-resolved manner.The SPIKE-distance (Kreuz et al., 2011) compares the spiketime accuracy between spike trains and uses instantaneous fir-ing rates to adapt to the local time scale. Finally, SPIKE-syn-chronization (Kreuz et al., 2015) is a discrete time-resolvedmeasure of similarity based on ISI-derived coincidence win-dows that are used to determine if two spikes from di ff erentspike trains are coincident or not. These measures have al-ready been successfully applied in many di ff erent contexts; forexample they have been used to detect determinism in pointprocesses (Andrzejak et al., 2014), to find correlations betweenspike trains and behaviour in an inverse neurocontroller (Dura-Bernal et al., 2016) and to evaluate a bio-inspired locomotionsystem in robotics (Espinal et al., 2016).Since they always adapt to the local firing rate, all three ofthese measures are time scale free. While they correctly iden-tify the relative firing rate di ff erences, they have no concept ofactual time scales and treat all time scales as equally important.This has the consequence that for very small time scales evenminor deviations from perfect synchrony lead to very high val-ues of dissimilarity. However, for real data the smallest timescales are often not very relevant and any dissimilarities therecan mostly be disregarded. Thus in this case the measures’ fo-cus on the local time scales results in a (spurious) amplificationof dissimilarities which compared to the global time scales arerather negligible.Here we address this problem by proposing generaliza-tions to the three measures called adaptive ISI-distance (A-ISI-distance), adaptive SPIKE-distance (A-SPIKE-distance) andadaptive SPIKE-synchronization (A-SPIKE-synchronization).These generalized definitions add a notion of the relative im-portance of local di ff erences compared to the global time scales.In particular, they start to gradually ignore di ff erences betweenspike trains for ISIs that are smaller than a minimum relevanttime scale (MRTS). The MRTS is implemented by an additionalvariable T which can either be defined as a parameter or esti-mated directly from the data.In some neuroscience applications only the similarity ofspike timing is important and rate di ff erentiation is not a de-sired property. While the A-ISI-distance is sensitive to firingrate alone and the A-SPIKE-distance responds to di ff erences inboth rate and timing, there is currently no measure that focusesonly on spike timing. Therefore, in a second step we extend theA-SPIKE-distance into the rate-independent adaptive SPIKE-distance (RIA-SPIKE-distance) which still identifies spike time di ff erences but ignores any rate deviations between the spiketrains.The remainder of this paper is organized as follows. In Sec-tion 2 we describe the generalized definitions of the three mea-sures, the A-ISI-distance (Section 2.1), the A-SPIKE-distance(Section 2.2), and A-SPIKE-synchronization (Section 2.3). InSection 2.4 we introduce a way to estimate the threshold valuedirectly from the data. We then investigate using both simulatedand real data how both the original measures and the adaptivegeneralizations deal with multiple time scales (Section 2.5). InSection 3 we add a rate-independent extension to A-SPIKE-distance (Section 3.1) and afterwards study the e ff ects of theextension (Section 3.2). The implications of the extensions arediscussed in Section 4. Finally, in the Appendix A we coversome non-trivial subtleties of the definitions for all three mea-sures. First we provide the definitions for the periods before thefirst and after the last spike in a spike train (where the interspikeinterval is not defined), and then we deal with the two specialcases of empty spike trains and spike trains with only one spike.The two experimental datasets used in Section 2.5 are describedin Appendix B.
2. Adaptive Generalizations
In this Section we introduce the adaptive generalizationsof the established measures ISI-distance (Kreuz et al., 2007),SPIKE-distance (Kreuz et al., 2011) and SPIKE-synchroniza-tion (Kreuz et al., 2015), which we will call A-ISI-distance, A-SPIKE-distance, and A-SPIKE-synchronization. All three gen-eralizations are built on a minimum relevant time scale (MRTS)which is implemented via the threshold parameter T . Thisthreshold is used to determine if a di ff erence between the spiketrains should be assessed in a local context or in relation to theglobal time scales. This threshold is used for all three measures,but the way it is applied varies. The generalized measures fallback on the original definitions when T =
0. In the followingthis is what we refer to whenever we talk of the original mea-sures. In this case even the smallest time scales matter and alldi ff erences are assessed in relation to the local context only.Note that the upcoming definitions only apply to the intervalbetween the first and the last spike. In Appendix A.1 and Ap-pendix A.2 they will be completed to range from the start ofthe recording t s to the end of the recording t e . Equally, some ofthe following equations are ill-defined when there are less thantwo spikes in a spike train. These special cases will be handledin Appendix A.3 and Appendix A.4.Throughout the paper we denote the number of spike trainsby N , indices of spike trains by n and m , spike indices by i and j and the number of spikes in spike train n by M n . The spiketimes of spike train n are denoted by { t ( n ) i } with i = , . . . , M n . The A-ISI-distance measures the instantaneous rate di ff er-ence between spike trains (see Fig. 1A). It relies on a time-resolved profile, meaning that a dissimilarity value is defined2 n A S p i k e t r a i n s t t (n)i t (n)i+1 t (m)j t (m)j+1 x ISI(n) (t)x
ISI(m) (t) mn B S p i k e t r a i n s t x (n)P (t) x (n)F (t)x (m)P (t) x (m)F (t) ∆ t i(n) ∆ t i+1(n) ∆ t j(m) ∆ t j+1(m) mn C t i(n) τ i(n) t j(m) τ j(m) t i(n) τ i(n) t j(m) τ j(m) Time S p i k e t r a i n s Figure 1: Schematic drawing for all three measures. (A) Illustration of thevariables that define the ISI-distance. The instantaneous interspike intervals x ( n ) IS I ( t ) are used as estimates of the local firing rate. (B) Additional variablesemployed in the definition of the SPIKE-distance. (C) Coincidence criterionfor SPIKE-synchronization. The coincidence window of each spike is derivedfrom its two surrounding interspike intervals. Here we illustrate two di ff erentexamples. The two spikes on the left side are considered coincident since bothlie in each other’s coincidence windows. On the right there is no coincidencesince the spike from the second spike train is outside of the coincidence windowfrom the spike of the first spike train. for each time instant. To obtain the profile, we assign to eachtime instant t the time of the previous spike t ( n )P ( t ) = max { t ( n ) i | t ( n ) i ≤ t } for t ( n )1 (cid:54) t (cid:54) t ( n ) M n (1)and the time of the following spike t ( n )F ( t ) = min { t ( n ) i | t ( n ) i > t } for t ( n )1 (cid:54) t (cid:54) t ( n ) M n . (2)From this for each spike train n an instantaneous ISI can becalculated as x ( n )ISI ( t ) = t ( n )F ( t ) − t ( n )P ( t ) . (3)For the A-ISI-distance we define the MRTS such that whenthe ISIs of both spike trains are smaller than a threshold value T , this value is used instead. The pairwise A-ISI-profile is thendefined as I An , m ( t ) = | x ( n )ISI ( t ) − x ( m )ISI ( t ) | max { x ( n )ISI ( t ) , x ( m )ISI ( t ) , T } . (4)The multivariate A-ISI-profile is obtained by averaging over all pairwise A-ISI-profiles I A ( t ) = N ( N − N − (cid:88) n = N (cid:88) m = n + I An , m ( t ) . (5)This is a non-continuous piecewise constant profile and a finalintegration over time gives the A-ISI-distance D AI = t e − t s (cid:90) t e t s I A ( t ) dt . (6)If the threshold T is set to zero, the generalized ISI-distance D AI falls back to the original ISI-distance D I . The A-SPIKE-distance measures the accuracy of spike timesbetween spike trains relative to local firing rates (see Fig. 1B).In order to assess the accuracy of spike events, each spike isassigned the distance to its nearest neighbor in the other spiketrain ∆ t ( n ) i = min j ( | t ( n ) i − t ( m ) j | ) . (7)These distances are then interpolated between spikes using forall times t the time di ff erences to the previous spike x ( n ) P ( t ) = t − t ( n ) i for t ( n ) i (cid:54) t (cid:54) t ( n ) i + , (8)and to the following spike x ( n ) F ( t ) = t ( n ) i + − t for t ( n ) i (cid:54) t (cid:54) t ( n ) i + . (9)These two quantities define a time-resolved dissimilarity profilefrom discrete values the same way as Eqs. 1 and 2 did for the A-ISI-distance. The instantaneous weighted spike time di ff erencefor a spike train can then be calculated as the interpolation fromone di ff erence to the next S n ( t ) = ∆ t ( n ) i ( t ) x ( n ) F ( t ) + ∆ t ( n ) i + ( t ) x ( n ) P ( t ) x ( n )ISI ( t ) , t ( n ) i (cid:54) t (cid:54) t ( n ) i + . (10)This function is analogous to the term x ( n )ISI for the ISI-distance,with the only di ff erence that it is piecewise linear instead ofpiecewise constant. It is also continuous.The pairwise A-SPIKE-distance profile is obtained by aver-aging the weighted spike time di ff erences, normalizing to thelocal firing rate average and, finally, weighting each profile bythe instantaneous firing rates of the two spike trains S Am , n ( t ) = S n x m ISI ( t ) + S m x n ISI ( t )2 (cid:104) x n , m ISI ( t ) (cid:105) max {(cid:104) x n , m ISI ( t ) (cid:105) , T } . (11)We define the MRTS by using a threshold, that replaces the de-nominator of weighting to spike time di ff erences if the meanis smaller than the threshold T . This profile is analogous tothe pairwise A-ISI-profile I An , m ( t ), but again it is piecewise lin-ear, not piecewise constant. Unlike S n ( t ) it is not continuous,since typically it exhibits instantaneous jumps at the times ofthe spikes. The multivariate A-SPIKE-profile is obtained the3 AA S ( t ) B C Time S A ( t ) Figure 2: An example spike train pair and its SPIKE-distance and A-SPIKE-distance profiles. (A) Two spike trains consisting of four events with five spikeseach. The sequence is the same for all four events, only the time scale is get-ting shorter and shorter. From a global perspective the first event consists ofnon-synchronous individual spikes, while the last event consists of coincidentbursts. The two events in the middle are intermediates. (B) The SPIKE-dis-tance considers only the local context and thus the profile shape is the same forall four events. (C) The A-SPIKE-distance takes into account also the globaltime scales. Like the SPIKE-distance it judges the first event as very dissimilar,but in contrast to the the SPIKE-distance it scales down the small spike timedi ff erences in the bursts and thus considers the coincident burst in the last eventas very similar. same way as the multivariate A-ISI-profile, by averaging overall pairwise profiles S A ( t ) = N ( N − N − (cid:88) n = N (cid:88) m = n + S Am , n ( t ) . (12)Finally, also the A-SPIKE-distance is calculated as the time in-tegral over the multivariate profile D AS = t e − t s (cid:90) t e t s S A ( t ) dt . (13)For T = ff ect of applying the threshold can be seen in Fig. 2.The first event of five spikes is compressed more and more un-til it becomes a single burst in the fourth event. The originalSPIKE-distance profile S ( T ) has the same proportions of dis-similarity for all events, since it uses local context only and thusconsiders all time scales as equal, while the A-SPIKE-distanceprofile S A ( t ) is scaled down when the di ff erences become smallcompared to the global time scales. A-SPIKE-synchronization quantifies how many of the possi-ble spike coincidences in a dataset are actually occurring (Fig.1C). While the A-ISI-distance and the A-SPIKE-distance aremeasures of dissimilarity which obtain low values for similarspike trains, A-SPIKE-synchronization measures similarity. Ifall the spikes are coincident with a spike in all the other spike trains, its value will be one. In contrast, if none of the spikesare coincident, it will be zero.The original SPIKE-synchronization (Kreuz et al., 2015) isparameter- and time scale-free, since it uses the adaptive coin-cidence detection first proposed for the measure event synchro-nization (Quian Quiroga et al., 2002). The coincidence win-dow, i.e., the time lag below which two spikes from two di ff er-ent spike trains, t ( n ) i and t ( m ) j , are considered to be coincident,is adapted to the local firing rate. Spikes are coincident only ifthey both lie in each other’s coincidence windows.For A-SPIKE-synchronization we generalize the definitionby introducing a threshold, which decides if the window is de-termined locally or if the global time scales should be taken intoaccount. As a first step, we define the ISI before the spike as x ( n ) iP = lim t → t i − x ( n )ISI ( t ) (14)and the ISI after the spike as x ( n ) iF = lim t → t i + x ( n )ISI ( t ) . (15)The coincidence window for spike i of spike train n is definedby determining the minimum coincidence window size for aspike as half the length of the two ISIs adjacent to the spike τ ( n ) i =
12 min { x ( n ) iP , x ( n ) iF } , (16)and allowing asymmetric coincidence windows based onMRTS. This is done by replacing τ ( n ) i with the threshold value T , if it is the smaller of the two. Since the threshold value isderived from ISIs and the coincidence window spans both sidesof the spike, only half of the threshold spans each side. For theA-ISI- and the A-SPIKE-distance the changes induced by thethreshold appear gradually, but for A-SPIKE-synchronizationthey occur as an abrupt jump from 0 to 1. Therefore, to com-pensate for the binary nature of A-SPIKE-synchronization, thethreshold is divided by two, resulting in an overall factor of1 /
4. The coincidence windows of neighboring spikes are notallowed to overlap, and thus each side is limited to half the ISIeven if the threshold is larger. Thus, the coincidence windowbefore the spike is determined as τ ( n ) iP = min { max( 14 T , τ ( n ) i ) , x ( n ) iP } (17)and the coincidence window after the spike as τ ( n ) iF = min { max( 14 T , τ ( n ) i ) , x ( n ) iF } . (18)The combined coincidence window for spikes i and j is thendefined as τ ( n , m ) i j = min { τ ( n ) iF , τ ( m ) jP } if t i (cid:54) t j min { τ ( n ) iP , τ ( m ) jF } otherwise . (19)The coincidence criterion can be quantified by means of acoincidence indicator C ( n , m ) i = j {| t ( n ) i − t ( m ) j |} < τ ( n , m ) i j . (20)4 S p i k e t r a i n S C AAAAAAAAAAAAAAAAAAAAAAAAAAA S p i k e t r a i n S CA BBBBBBBBBBBBBBBBBBBBBBBBBBB
Time S p i k e t r a i n Diff
CCCCCCCCCCCCCCCCCCCCCCCCCCC
Figure 3: SPIKE-synchronization (A), A-SPIKE-synchronization (B) and theirdi ff erence (C) illustrated using five spike trains with four simple events. For theoriginal measure (A) the small interspike intervals of spike doublets (first andsecond event) or bursts (third event) result in an unreasonably high demandfor spike timing accuracy. With the adaptive generalization (B) for all thesecases the likelihood increases that at least one of the spikes is part of a coinci-dence. On the other hand, if there are no doublets or bursts (last event), nothingchanges (best seen in C). Note that the color scales di ff er, for better visibilitywe use grey-black in A and B but white-black in C. This definition ensures that each spike can only be coincidentwith at most one spike in the other spike train. The coincidencecriterion assigns either a one or a zero to each spike dependingon whether it is part of a coincidence or not. For each spike ofevery spike train, a normalized coincidence counter C ( n ) i = N − (cid:88) m (cid:44) n C ( n , m ) i (21)is obtained by averaging over all N − i in spike train n .This way we have defined a coincidence indicator for eachindividual spike in the spike trains. In order to obtain one com-bined similarity profile, we pool the spikes of the spike trains aswell as their coincidence indicators by introducing one overallspike index k . This yields one pooled set of coincidence indica-tors { C k } = (cid:91) n { C ( n ) i } (22)from which the A-SPIKE-synchronization profile C A ( t k ) can beobtained via C A ( t k ) = C ( k ). Finally, A-SPIKE-synchronizationis defined as the average value of this discrete profile S AC = M M (cid:88) k = C A ( t k ) , (23)where M is the overall number of spikes. In Fig. 3we illustrate how the asymmetric coincidence windows of A-SPIKE-synchronization allow for a better coverage ofburst events which makes it easier to match spikes whencompared to the original SPIKE-synchronization (A-SPIKE-synchronization with T = ff erences below threshold adds coincidences and thus,since it is a measure of similarity, A-SPIKE-synchronizationcan only increase. In neuroscience typical time scales are in the range of mil-liseconds or sometimes seconds and any time scales below thiswill not be considered relevant. In fields such as meteorol-ogy the respective time scales could be hours and days or evenmonths and years. The relevant time scales clearly depend onthe system under consideration. Setting the minimum relevanttime scale (MRTS) for a given dataset might not be a simpletask. To address this, we propose a method to extract a thresh-old value from the spike trains, that is based on the proportionsof the di ff erent time scales present in the data.It is important to note that the selected MRTS is not an indi-cator of a time scale of the system; it just determines the out-come of the adaptive generalizations. It is also not a hard setlimiter neglecting everything below the threshold, but rather itmarks the time scale from which on di ff erences are consideredin the global instead of the local context. Thus from this timescale on deviations from synchrony are treated as less and lessrelevant the smaller they get, even if they are large in relation tothe local time scales.The purpose of the threshold is to act as an indicator of whatglobally is a high rate or inversely a small ISI. The originalnormalizations are based on the ISIs, so it is reasonable to de-termine the threshold from the pooled ISI-distribution. We usethe ISIs after the edge e ff ect has been corrected (see AppendixA.1). The threshold should fulfill two main criteria. First, itneeds to decrease proportionally to the spike count, so that in-creasing rates (or longer recordings with the same rate) do notchange the threshold. Second, the threshold should respond tochanges in the ISI-distribution so that it is able to adapt betweensingle and multiple time scale data sets. In Fig. 4 we use a sim-ple spike train motive of just four spikes to illustrate these twocriteria.The most straightforward threshold would be the meanlength of the ISIs (cid:104) L ISI (cid:105) = (cid:80) Gg = L g ISI M IS I = LM IS I . (24)Here L g ISI denotes the ISI-length and M IS I is the total number ofISIs in the pooled ISI-distribution. In the numerator the sum ofthe lengths of all ISI equals the overall length L of the pooledISIs. Apart from edge e ff ect corrections this is equal to theproduct of recording length and number of spike trains whichis a constant. Thus while the mean of ISIs depends on the num-ber of spikes (Fig. 4A), for a given number of spikes (numberof ISIs) it is completely independent of how the ISIs are dis-tributed around the mean (Fig. 4B). It adapts to the spike count5 A 1 Iteration T i m e [ a . u .] ThresholdMean
B 1 Iteration T i m e [ a . u .] Figure 4: Threshold value vs. the mean of the ISI-distribution. (A) Dependenceon the number of spikes (first criterion). In each iteration the number of spikesis increased by concatenating two half-length copies of the previous iteration.Both the mean and the threshold decrease with spike count. (B) Dependenceon the ISI-distribution (second criterion). From iteration to iteration the ISI-distribution is changed by halving the three short ISIs and prolonging the longISI accordingly. Since the spike count (and thus the number of ISIs) is keptconstant, the mean does not respond to this change. However, the thresholdcorrectly increases with the heightened importance of the long ISI. but not to the proportions in which the ISIs appear in the datathus fulfilling the first but not the second criterion.To fulfill both criteria one needs to not just count the inter-spike intervals but weight them by their length. This reduces theimportance of short ISIs and allows the long ISIs to influencethe threshold according to their contribution and not just num-ber. It is equivalent to taking the mean of the second momentsof the ISIs T = (cid:112) (cid:104) ( L ISI ) (cid:105) = (cid:115) (cid:80) Gg = L g ISI2 M IS I . (25)Note that in order to obtain a value with the right dimension thesquare root of the average must be taken. This threshold valuehas roughly the same dependence on the number of spikes asthe mean value (Fig. 4A), however, in contrast to the mean itis also sensitive to changes in the ISI-distribution. In summary,using T as the MRTS fulfills both criteria set for the threshold. In this Section we investigate how both the adaptive gener-alizations (with automated thresholding) and the original mea-sures deal with multiple time scales. For the A-ISI-distance andthe A-SPIKE-distance we use a test spike train set consisting ofsimulated and real spike trains to study the e ff ect of the gener-alized versions (Section 2.5.1). After that, in Section 2.5.2, westudy on real MEA recordings how A-SPIKE-synchronizationdi ff ers from SPIKE-synchronization. In Section 2.5.3 we sys-tematically test the influence of the amount of bursts on the S p i k e t r a i n s Figure 5: Spike train test set used to compare the generalized versus the orig-inal measures of spike train synchrony. Spike trains 1-25 are artificially con-structed examples which cover a range of archetypical spiking patterns, whereasspike trains 26-30 are selected examples of neuronal spiking data from a neu-ronal culture recorded on a micro electrode array (see Appendix B.1). All spiketrains are normalized by their total length. di ff erence between adaptive and original measures. Finally, inSection 2.5.4 we investigate how the adaptive versions changethe analysis of neuronal reliability in an experimental dataset. We address two points. First, we look at the sensitivity of theadaptive generalizations and verify that in the presence of burststhey perform better than the original measures. Second, we alsomake sure that the changes are specific, e.g., we confirm that inall other cases and especially if there are no bursts, the adaptivegeneralizations do not exhibit unwanted side e ff ects.To this aim, we use a test set composed of both artificial andreal spike trains (Fig. 5) to compare A-ISI-distance to ISI-dis-tance and A-SPIKE-distance to SPIKE-distance. We use twomodels to generate our samples. For the spike trains with per-fect periodicity we use a time varying steady rate (fixed ISI)model. For samples with more variability in spike timing weused a Poisson spiking model, where the rate is fixed for a cer-tain window at a time. In some cases we add small jitter noiseto both models. The artificial spike trains 1-25 are designed toexhibit a variety of stereotypical spiking behaviours includingboth single and multiple time scales. The experimental spiketrains 26-30 consist of short recordings from neuronal cultureson microelectrode arrays (see Appendix B.1 for details). Forthe adaptive versions the threshold is estimated from the data(see Section 2.4) for each pair separately.In the analysis every spike train is paired with all the others.Because for both the A-ISI-distance and the A-SPIKE-distancethe MRTS T can only reduce but never increase the dissimi-larity value, all pairs are found in the lower half of the scatterplot (Fig. 6). Furthermore, all values between pairs of spiketrains are close to or on the diagonal, which means that bothversions attain very similar values or even the same value. The6i ff erences between the two SPIKE-distances are slightly morepronounced than the di ff erences between the two ISI-distances.Such seemingly small di ff erences can still be of high signifi-cance since in a typical experimental setup it is rarely the abso-lute value of similarity that matters but rather the relative orderof similarity for di ff erent conditions. Moreover, in real data therange of similarity values obtained is usually quite small whichfurther increases the relative importance of small changes insimilarity.For one spike train at a time we then look at all its pairingsand sort the results in ascending order according to the originalversions. The results from the adaptive versions are arrangedin the same order. If the order of the spike train pairs does notmatch, there is a clear di ff erence in the way the two measuresconsider spike train similarity.We now investigate in more detail not only the largest abso-lute, but also the largest relative changes observed in Fig. 7.First, the largest absolute changes are identified by calculatingthe Euclidean distances between the results for the two spiketrain pairs. They typically take place for pairs of spike trainswith large distances. Next, since deviations from near perfectsynchrony are more prominent and easier to detect than di ff er-ences between various levels of high dissimilarity, we also lookat relative changes. These can be found by dividing each dis-tance by its corresponding ISI- and A-ISI-distance or SPIKE-and A-SPIKE-distance average. For both distances they mostlyoccur for pairs of very similar spike trains.For the A-ISI-distance, the spike train pairs showing thelargest absolute change compared to the ISI-distance can beseen in Fig. 7A. The two measures show a di ff erent order ofsimilarity; while the ISI-distance increases, the A-ISI-distancedecreases from the first to the second pair. Spike trains 27 and30 in the first pair are seen very similarly (deviation < < ff erence is found for the second spike train pair. Here the ISI-distance looks at the detailed structure and judges the interspikeintervals within the bursts as very dissimilar, whereas the A-ISI-distance sees simply matching bursts and attains a considerablylower distance value than the ISI-distance (0.100 vs 0.129).For the A-SPIKE-distance, the pairs of spike trains showingthe largest absolute change compared to the SPIKE-distance aredepicted in Fig. 7C. As there are no bursts in either of the twospike trains, both measures attain exactly the same value for thefirst spike train pair. This is a very good example for the speci-ficity of the generalized version. On the other hand, the originaldistance considers the second spike train pair (periodic spik-ing versus periodic bursts) as much more dissimilar (increase > ff erences in the middle of two bursts as largerthan the di ff erences in the middle of the burst.Finally, the largest relative change between the two SPIKE-distances is shown in Fig. 7D. Again, there is not much dif-ference between the two distances for the first spike train pair.However, the SPIKE-distance considers the second spike trainpair much more dissimilar ( >
62% higher) due to the large rel-ative deviations in spike timing within their coinciding bursts.In contrast, the A-SPIKE-distance puts much less weight on thedi ff erences within bursts, but still reacts to the spikes outside ofthe bursts. This is an example of the sensitivity of A-SPIKE-distance.All these results show that the e ff ect of both generalized ver-sions is strongest in situations with multiple time scales in thespike trains. A prominent example are bursts embedded in longsilent periods. In this case the long ISIs (of the silent periods)strongly influence the global time scales such that deviationsof synchrony on the smallest time scales (within the bursts) areweighted less. A-SPIKE-synchronization can not be meaningfully tested byusing the spike train set of Fig. 5. The perfect periodicity inmany spike trains makes analysis of the A-ISI-distance and theA-SPIKE-distance simple, but causes very abrupt changes inA-SPIKE-synchronization due to its binary nature. The valuescan be computed, but the largest di ff erences are not meaningfulwith this data set, since many spike trains with bursts jump fromzero to a large value and there is no way of ordering di ff erentpairs having zeros in the original measure. Thus, we here use aqualitative approach together with insights from the analysis ofA-ISI-distance and A-SPIKE-distance.As a side e ff ect of being time scale adaptive, SPIKE-syn-chronization demands very high spike timing accuracy duringfast firing. This leads to situations such as the one shown inFig. 3 (Section 2.3). For spike trains 3 and 4 the spikes inthe first event are considered coincident, but the doublet in be-tween them in spike train 2 is not judged as coincident witheither of them. In contrast, for A-SPIKE-synchronization thecoincidence windows are adapted to the distribution of all ISIsin the data set and the two sides of the coincidence window areallowed to be of di ff erent length. With this change each of thespikes in the doublet becomes coincident with one of the spikes(the respective closest one) in spike trains 3 and 4.As a by-product of the adaptation, A-SPIKE-synchronizationalso increases the coincidence window coverage within and atthe edges of a burst and thus matches as many spikes as pos-sible. For SPIKE-synchronization many of these spikes wouldbe ignored due to the unreasonably small coincidence windows.This phenomenon occurs very often with real data. An exam-ple containing two small and one large burst event is shown inFig. 8. In the first two events A-SPIKE-synchronization is ableto detect a few additional coincidences compared to SPIKE-synchronization. The di ff erence is much more pronouncedfor the third and largest event. Here for SPIKE-synchroniza-tion many potential matches are left out and this leads to a7 .2 0.4 0.6 0.8 1ISI-distance0.20.40.60.81 A - I S I - d i s t an c e A A - SP I KE - d i s t an c e B Figure 6: Scatter plots showing the A-ISI-distance between all pairs of spike trains versus the original ISI-distance (A) and the A-SPIKE-distance versus theSPIKE-distance (B). The diagonal line marks where the measures would show equal distance. The pairs were sorted according to rising order of the originaldistances- Thus, if the order changed for adaptive extension, there is a negative slope in a line connecting all the pairs. For each line with a negative slope wecalculated its length using the Euclidean distance. The five largest absolute changes are indicated in magenta, the five largest relative changes in cyan. In addition, amagenta (cyan) arrow points to the very largest absolute (relative) change. Overall, while the changes are seemingly small on an absolute scale, the relative changescan be very significant. rather low overall value of 0.238. Instead, when A-SPIKE-synchronization is used, there are almost 45% more matchedspikes within the burst and this strongly increases the overallsynchronization value to 0.345.Fig. 8C clearly shows that the additional spike matching ofA-SPIKE-synchronization only occurs in the high frequencyevents for which small di ff erences in the ISIs cause gaps be-tween the coincidence windows of adjacent spikes. Coinci-dences outside of these high frequency events are not a ff ected. Next we test how the e ff ect of the automated thresholdchanges when spikes are forming tighter bursts (Fig. 9). To dothis we first create two Poisson spike trains which are dividedinto four equally long segments. These segments are then in-creasingly compressed which prolongs the ISIs between themsuch that the total length remains constant (Fig. 9A). We usethe relative length of the interburst intervals R as a parameterand track the di ff erence between the adaptive and the origi-nal versions. The results for the ISI-distance and the SPIKE-distance are very similar and we only show the latter. FromFig. 9B we can see that the SPIKE-distance decreases almostlinearly with R since the relative importance of the common si-lence in the interburst intervals increases. The adaptive versiondecreases sub-linearly with the largest absolute di ff erence be-tween the two measures occurring around R = .
4. For higher R -values the reduction of the burst length overshadows the in-creases in similarity at burst times and the di ff erence increasesup to a point and then starts to decrease. The relative di ff erence increases over the whole interval (data not shown, but can be ap-preciated by observing the di ff erence approaching the SPIKE-distance value towards R = ff ect on the overall value, SPIKE-synchronization is sensitiveto the matching of spikes only and is based on one coincidenceindicator value per spike. Thus, the e ff ect is increasing onlyuntil all possible spike pairs within the bursts are matched. Forour example the increase saturates at R = . In order to demonstrate the e ff ects of the adaptive generaliza-tion in a more realistic scenario, we re-analyze data previouslyused to study the e ff ect of membrane potential resting state onneuronal reliability (Zeldenrust et al. (2013), see Appendix B.2for details on the recordings). When in the original study frozennoise was injected into thalamocortical relay cells of rats, it wasfound that the reliability of the cell response increases with de-polarization (Zeldenrust et al., 2013).Here we use both the original versions and the adaptive gen-eralizations of all three measures to assess the reliability of theresponses from the two neurons for which all four levels ofmembrane potential were recorded. The adaptive versions usea threshold obtained from the data by pooling all spike trains8 B S p i k e t r a i n s Time2730 A S p i k e t r a i n s Time 0.51 ISI−distance 0.126A−ISI−distance 0.1200 2 4 6 8 100.51 ISI−distance 0.129A−ISI−distance 0.100 I ( t ) ,I A ( t ) Time0.51 ISI−distance 0.852A−ISI−distance 0.8480 2 4 6 8 100.51 ISI−distance 0.883A−ISI−distance 0.847 I ( t ) ,I A ( t ) Time 515 D S p i k e t r a i n s Time1715 C S p i k e t r a i n s Time 0.51 SPIKE−distance 0.029A−SPIKE−distance 0.0280 2 4 6 8 100.51 SPIKE−distance 0.047A−SPIKE−distance 0.027 S ( t ) , S A ( t ) Time0.51 SPIKE−distance 0.333A−SPIKE−distance 0.3330 2 4 6 8 100.51 SPIKE−distance 0.368A−SPIKE−distance 0.331 S ( t ) , S A ( t ) Time
Figure 7: Largest absolute and largest relative change between ISI-distance and A-ISI-distance (A,B) as well as SPIKE distance and A-SPIKE-distance (C,D)for the spike train set shown in Fig. 5. In both cases the two measures show a di ff erent order of similarity. The original distances increase from the first to thesecond pair, while the adaptive extensions decrease. The first pairs attains distance values in between the second pairs, which results in a di ff erent order of similaritybetween the measures (see Fig. 6). Upper subplots: The spike train pairs with the largest changes. Lower subplots: Respective original distance vs. adaptive versionprofiles with the di ff erence between the two profiles emphasized. The distance values for the first (second) pair are shown on top (at the bottom). They are alsomarked by a dashed line for the original distance and by a solid line for the adaptive distance. S p i k e t r a i n SPIKE-synchronization 0.238 A S C S p i k e t r a i n A-SPIKE-synchronization 0.345 B S CA
91 92 93 94 95 96 97
Time [s] S p i k e t r a i n Difference 0.107 C Diff
Figure 8: Real data example from MEA recordings (see Appendix B.1 formore details). Ten spike trains from the data set are plotted and their coinci-dence windows are drawn as obtained by SPIKE-synchronization (A) and A-SPIKE-synchronization (B). The di ff erence is plotted in C. Due to the adaptivecoincidence windows, A-SPIKE-synchronization is able to match around 45%more spikes between bursting spike trains than SPIKE-synchronization. As inFig. 3, the color scale is grey-black in A and B and white-black in C. of each level and trial together. In Fig. 10 we show the re-sults of the cell with the more prominent e ff ect but we get sim-ilar results for the other cell as well. The cells analyzed wererecorded three times for each holding membrane potential andreliability was assessed by trial to trial variations. For the high-est hyperpolarization (Fig. 10A) the original SPIKE-distanceyields spuriously high values for the local dissimilarity dur-ing the bursts, since it only evaluates the local context. Evenwhen the A-SPIKE-distance takes the global context into ac-count, both measures agree that there are large dissimilaritiesin the spike trains.For the most depolarized state (Fig. 10B) the membrane po-tential is considerably closer to the action potential threshold.The patterns are closely matching the burst positions of Fig.10A, but also additional events appear. The neuron no longerresponds in clearly distinguished bursts and it is considerablymore di ffi cult to determine where a burst begins or ends. Sincethe generalized version adapts to time scales found in all thespike trains, it is able to distinguish when a burst-like patternemerges and considers them as more similar.As can be seen in Figs. 10C and 10D, the original versions,without adaptation and only using the local context, attain ahigher level of similarity for -60mV than for -50mV, whichcontradicts both the results in the original study and the re-sults for SPIKE-synchronization (Fig. 10E). Since the adap-tive versions are able to make use of the global context of allthe spike trains, they attain results without this spurious dissim- ilarity and thus for higher membrane potentials the similarityincreases monotonously.A-SPIKE-synchronization works slightly di ff erently (Fig.10E). Due to the tight bursts that cause excessively small co-incidence windows, the largest di ff erence occurs for the hyper-polarized states. However, both versions agree that the relia-bility as quantified by spike to spike matching in the responsepatterns clearly show a monotonous increase over the baselinemembrane potential. In summary, the results obtained by theA-SPIKE-distance and A-ISI-distance seem to be appropriateand more in line with the original results.
3. Rate-independent extension
Sometimes in neuroscience one is interested in the pure sim-ilarity of spike timing, independent of any di ff erences in spikerates. Thus there is the need for a measure which can identifydi ff erences in spike timing but is able to ignore any di ff erencesin rate between the spike trains. Here we propose such a rate-independent extension for the A-SPIKE-distance. In order to understand how rate-independence for A-SPIKE-distance is achieved, we need to separate Eq. 11 (Section 2.2)for the pairwise A-SPIKE-distance profile into its three compo-nents.The first two components are the mean of spike time dissim-ilarity and the normalization to firing rate S m , n ( t ) = S n ( t ) + S m ( t )2 · max {(cid:104) x n , m ISI ( t ) (cid:105) , T } , (26)where S n ( t ) and S m ( t ) are the weighted spike time di ff erencesfor spike trains n and m defined by Eq. 10. The third componentis a weighting of the spike time dissimilarity according to thefiring rate di ff erence that is applied to the first component S n ( t ) x m ISI ( t ) + S m ( t ) x n ISI ( t ) (cid:104) x n , m ISI ( t ) (cid:105) . (27)The rate-independent adaptive SPIKE-distance (RIA-SPIKE-distance) simply leaves out this last weighting and can thus bewritten as S RIAm , n ( t ) = S n ( t ) + S m ( t )2 max( (cid:104) x n , m ISI ( t ) (cid:105) , T ) . (28)The RIA-SPIKE-distance shares all the properties of A-SPIKE-distance, but it only evaluates normalized spike timing di ff er-ences, whereas the A-SPIKE-distance additionally uses di ff er-ences in rate to determine similarity. The A-ISI-distance is a measure of instantaneous rate di ff erence and arate-independent measure of rate di ff erence makes little sense. A-SPIKE-synchronization is rate-dependent by definition, since it is calculated as theaverage value of spike-based coincidence indicators (Eqs. 20 and 23). igure 9: E ff ect of bursts on the adaptive versions evaluated by using the relative length R of interburst intervals. The values are averages over 10realizations. (A) Five spike train pairs with increasing levels of burstiness for one example realization. (B) E ff ect of burstiness on the di ff erence betweenA-SPIKE-distance and SPIKE-distance. The graph for the ISI-distance looks very similar and is thus omitted. (C) Equivalent results for A-SPIKE-synchronization. The R-values of the examples in (A) are marked in (B) and (C) as dotted vertical lines. In this Section we compare the RIA-SPIKE-distance to theregular A-SPIKE-distance regarding their response to di ff er-ences in rate. First, in Fig. 11A we look at Poisson spiketrains with di ff erent rate ratios. The regular A-SPIKE-distanceexhibits a clear rate dependency obtaining its lowest value forspike trains with identical rates and increasing for higher ratedi ff erences. The RIA-SPIKE-distance on the other hand startsnear 0 .
25 and remains relatively constant for all rate ratios.These deviations from perfect rate-independence occur becauseof the irregularities of the Poisson spike trains. When we re-peat the same analysis with steady rate instead of Poisson spiketrains (Fig. 11B), thereby removing the e ff ects of the Poissonstatistics, the RIA-SPIKE-distance exhibits indeed perfect rate-independence.Regarding the original distances, in Fig. 11A they wouldshow very similar behavior to the adaptive generalizations.Only for rate ratios close to 1 there would be a small increasedue to coincident burst-like events within the Poisson spiketrains. In Fig. 11B the curves would overlap perfectly sincethere is only one time scale in steady rate spike trains (both re-sults not shown).
4. Discussion
In this manuscript we introduce adaptive generalizations tothe three existing measures ISI-distance, SPIKE-distance andSPIKE-synchronization as well as a rate-independent extensionto the generalized SPIKE-distance. These new measures ad-dress two distinct problems. The adaptive generalizations allow to disregard spike timedi ff erences that are not relevant on a more global scale. Bymeans of a specifically constructed library of both stereotyp-ical and real data spike trains, we can show that both A-ISI-distance and A-SPIKE-distance indeed yield improvements forpairs of spike trains containing di ff erent time scales withoutexhibiting any unwanted side e ff ects in other examples. Thusthe changes are both sensitive and specific. Regarding the sizeof the changes, even if they are seemingly small on an abso-lute scale, the relative changes can be very significant. Forour test set the largest relative change reaches 29% for the A-ISI-distance and even up to 62% for the A-SPIKE-distance.With a more qualitative approach we then show that A-SPIKE-synchronization fixes the problem of SPIKE-synchronizationwhich demands an unreasonably high accuracy for spike dou-blets and coinciding bursts. By introducing a global referenceframe, it manages to match spikes more e ffi ciently (for our testdata we found an increase of 45%).In order to test the adaptive measures methodologically wetested them in a controlled environment where two Poissonspike trains were split into bursts using increasingly large in-terburst intervals. We designed the adaptive extension to besensitive to bursting structure, therefore for increasing relativelength of interburst intervals we expect a larger di ff erence be-tween the original and the adaptive versions. We show that therelative di ff erence indeed increases monotonously with an in-crease in the ratio between interbursts interval and bursts.The absolute di ff erence obtains its maximal value when thedi ff erences ignored in the bursts are large and the bursts are longenough in comparison to the total length of the recording. Whenvery similar spike trains are compared their relative di ff erence11 S ( t ) , S A ( t ) s p i k e t r a i n s
20 21 22 23 24 25 26 27 28 29 3000.20.4 Time [s] S ( t ) , S A ( t ) −80 −70 −60 −5000.020.040.060.080.1 Voltage [mV] D i s t an c e BC ISIA−ISI −80 −70 −60 −5000.020.040.060.080.1 Voltage [mV] D i s t an c e D SPIKEA−SPIKE −80 −70 −60 −5000.20.40.60.81 Voltage [mV] D i s t an c e E SPIKE sync.A−SPIKE sync.
Figure 10:
Analysis of neuronal responses to multiple presentations of frozen noise for four di ff erent levels of the membrane potential. (A) Spike trainresponses (top) to two of the three noise presentations for a membrane potential of − mV and corresponding profiles for both A-SPIKE-distance andSPIKE-distance (bottom). The di ff erence between the two profiles is marked in red. (B) Same as in (A) but for a membrane potential of − mV . Resultsof the original and the adaptive measures for spike train sets of all three trials at four di ff erent voltage levels for the ISI-distance (C), the SPIKE-distance(D) and SPIKE-synchronization (E). becomes dominant and internal structures of coinciding burstsbecome less relevant.Additionally, we apply the measures to a dataset previouslyanalyzed for reliability and find that the adaptive methodsagree with the previous results better than the original ver-sions. The A-ISI-distance and the A-SPIKE-distance seemto yield more reasonable results than the ISI-distance and theSPIKE-distance. On the other hand when the coincidence win-dows of the original version get spuriously small, A-SPIKE-synchronization can match spikes much more e ffi ciently. Thee ff ect can be especially meaningful in applications in whichleader-follower relationships based on the temporal order ofspikes are determined (Kreuz et al., 2017).The rate-independent extension on the other hand focuseson spike time accuracy while disregarding rate di ff erences inthe two spike trains. The original SPIKE-distance considersspike time di ff erences but also has a feature that takes into ac-count the firing rate di ff erence between the spike trains. How-ever, sometimes only the spike time accuracy is of interestand for that purpose the RIA-SPIKE-distance disregards anydeviations in firing rate. We can show that the RIA-SPIKE-distance is approximately rate-independent for Poisson spike trains (apart from minor statistical e ff ects) and perfectly rate-independent for strictly periodic spike trains. With this finaladdition we have completed the picture, since we now havemeasures that are sensitive to rate only (A-ISI-distance), to tim-ing only (ARI-SPIKE-distance), and to both at the same time(A-SPIKE-distance).The adaptive generalizations are implemented for caseswhere we have prior knowledge of the system or where wewant to reduce the importance of very small details. However,one has to be careful with this method. If the threshold param-eter that defines the minimum relevant time scale (MRTS) ischosen too high, this can introduce spurious synchrony. To fa-cilitate the selection, we introduce a method for automaticallyextracting the threshold from the spike train data. This is doneby using the second moment the ISI-distribution of the wholedataset, thereby giving more weight to longer ISIs.Here it is important to note that while this automated esti-mation of MRTS gives us a threshold value for each dataset,one has to be very careful when comparing results obtainedwith di ff erent threshold values. Thus, one cannot use the adap-tive version for two recordings from the same source with-out using the same threshold for both recordings, even if the12 D i s t an c e A -3 -2 -1 Rate ratio00.20.40.60.81 D i s t an c e B A-ISIA-SPIKERIA-SPIKE
Figure 11: Rate-independent RIA-SPIKE-distance vs. A-SPIKE-distance andA-ISI- distance for Poisson (A) and steady-rate spike trains (B). (A) Distancesfor two Poisson spike trains with varying rate ratios. The overall number ofspikes in the two spike trains is kept constant. Each data point is an average over100 trials. In contrast to the clearly rate-dependent A-ISI- and A-SPIKE-dis-tances, the rate-independent RIA-SPIKE-distance exhibits an almost constantcurve. (B) For the steady-rate spike train curves each data point is an averageover 100 trials with random phase shifts between the two spike trains. In thiscase the line for the RIA-SPIKE-distance is indeed perfectly constant.
ISI-distributions di ff er. In such cases, the preferable optionwould be to combine the ISI-distributions before calculatingthe threshold and to use the resulting value for both record-ings. However, this might not work in all cases. For example,recordings before and during an epileptic seizure can have verydi ff erent ISI-distributions. This means that a globally mean-ingful threshold can not be extracted due to a very bi-modaldistribution of all the ISIs from the whole recording. The re-sulting threshold would be in between the two modes whichwould cause the adaptive measures to basically consider one ofthe recordings as a long burst and the other as an almost silentperiod. Thus, in cases where a suitable threshold can not befound, it is preferable to just set it to zero and consider onlylocal information. This is equivalent to using the original ver-sions.Many time scale parametric measures like the Victor-Purpuraand the van Rossum distance use a parameter to define the timescale of the system. The threshold set for the adaptive versionsis philosophically di ff erent in the sense that it does not define a single time scale, but sets a line below which the e ff ects ofthe smaller time scales are being toned down. All di ff erent timescales are still considered at the same time, but weighted di ff er-ently depending on how they compare to the threshold.Other measures that deal with multiple time scales exist. Forexample, Lyttle and Fellous have proposed a metric to specifi-cally assess the similarity of spike trains with bursts or commonsilent periods (Lyttle and Fellous, 2011). While in the proposedadaptive measures the time scale parameter is limiting full timescale independence of the original measures, in many measuresthe time scale is a fixed value. With the method proposed byLyttle and Fellous they can detect bursts as well as silent pe-riods. However, this comes with a cost, since the method re-quires two time scale parameters and three additional parame-ters; length of minimum silent period, length of burst ISI, min-imum number of spikes in a burst, scaling factor to decide howimportant bursts are in comparison to single spikes, and an-other factor to decide between importance of burst and silentperiod detection. While the large array of options gives the ex-perimenter a powerful tool and provides more control over theanalysis, it also increases the complexity of the overall experi-ment. This may cause problems, in particular when the data hasmany dimensions. Similarly, Rusu and Florian have introduceda new class of metrics (Rusu and Florian, 2014). The max-metric and the modulus-metric are well suited for measuringdistances between spike trains where information is encodedin bursts but single spike accuracy within burst is not relevant.The max-metric depends on the kernel chosen and a time scaleparameter deciding its size. The modulus-metric is parameterfree like the ISI-distance, the SPIKE-distance, and the SPIKE-synchronization. This is achieved by using a very simplifiedkernel. However, the results obtained with both methods arenot normalized. Thus based on the dissimilarity value alone itis not possible to say anything about the similarity of the twospike trains, but only about the order of di ff erent pairs.Another often used alternative to spike train distances arecorrelation measures (see e.g. Cutts and Eglen, 2014). How-ever, these measures traditionally require windowing or binningand this creates the problem that their performance can dependcrucially on the window length or bin size and also on the start-ing points and the overlap of the windows which clearly reducesthe objectivity of the results.The results confirmed our initial expectation that the maindi ff erences between the adaptive generalizations and the orig-inal measures is in their assessment of the similarity of burstydata. Since bursts are ubiquitous and have been identified asan important area of neuroscience research (see e.g. Izhikevichet al., 2003; Sherman, 2001), there is a strong need for this kindof similarity measurement. For the ISI-distance, a method hasbeen proposed for evaluating the similarity of bursty data byidentifying bursts and assigning spikes at the beginning of thebursts (Qu et al., 2016). However, burst detection is a notori-ously di ffi cult problem for which rather complicated methodshave been developed (see for example Kapucu et al., 2012).Thus, a measure based on assigning spikes to bursts inherits theproblems of burst detection. Another problem with the mea-sure proposed in Qu et al. (2016) is that it disregards di ff erences13n spiking behavior within the bursts. In contrast, our adaptiveversions do not detect bursts at all, but automatically adapt theirbehavior whenever there are burst-like features in the data.All the measures presented here are symmetric and thus in-variant to the order of the spike trains. Recently we have de-veloped a complementary directional approach consisting oftwo new measures called SPIKE-order and Spike Train Order(Kreuz et al., 2017). This approach utilizes the adaptive coinci-dence detection of SPIKE-synchronization to first sort multiplespike trains from leader to follower and then to quantify theconsistency of the spatio-temporal propagation patterns. A nat-ural continuation of the work presented in this article would beto use the adaptive measures for this new approach as well.We would like to finish by pointing out that the implemen-tations of both the original and the extended measures areprovided online in three separate free code packages calledSPIKY (Kreuz et al., 2015) (Matlab GUI), PySpike (Python)(Mulansky and Kreuz, 2016) and, most recently, cSPIKE (Matlab command line with MEX-files). Acknowledgements
E.S., I.M., and T.K. acknowledge support from the European Union’s Hori-zon 2020 research and innovation program under the Marie Sklodowska-CurieGrant Agreement No. 642563 ’Complex Oscillatory Systems: Modeling andAnalysis’ (COSMOS). M.M., N.B., and T.K. acknowledge support from theEuropean Commission through the Marie Curie Initial Training Network ’Neu-ral Engineering Transformative Technologies’ (NETT), project No. 289146.K.L. was supported by the 3DNeuroN project in the European Union’s Sev-enth Framework Programme, Future and Emerging Technologies [Grant Agree-ment No. 296590] and by the Tekes funded Human Spare Part project.We thank Emre Kapucu, Inkeri V¨alkki and Jari Hyttinen from the Hyttinenlab, BioMediTech, Tampere, Finland for the data from the MEA-cultures,andWytse Wadman and Pascal Chameau from Cellular and Systems Neurobiology,Swammerdam Institute for Life Sciences, University of Amsterdam, Amster-dam, the Netherlands for the current-clamp data. Finally, we thank Ralph G.Andrzejak for useful discussions and for carefully reading the manuscript.
Appendix A. Edge e ff ect correction and treatment of spe-cial cases Here, we deal with some subtle details in the definitions ofall three measures A-ISI-distance, A-SPIKE-distance and A-SPIKE-synchronization. First, in Appendix A.1 and AppendixA.2, we correct the edge e ff ect by providing definitions for theperiods before the first and after the last spike in a spike train(for which the interspike interval is not defined). This is nec-essary to guarantee that all measures are well-defined for thewhole recording interval. Subsequently, in Appendix A.3 andAppendix A.4, we deal with the two special cases of emptyspike trains and spike trains with only one spike. Even if somespike trains are empty or very sparse, all measures should stillbe defined in a way which is consistent with the regular defini-tions. http: // / users / thomas.kreuz / Source-Code / SPIKY.html http: // mariomulansky.github.io / PySpike / http: // / users / thomas.kreuz / Source-Code / cSPIKE.html Appendix A.1. Edge e ff ect correction for A-ISI- and A-SPIKE-distance Since the A-ISI- and the A-SPIKE-distance are time-resolvedand are based on ISIs defined by Eq. 3, there is ambiguity at theedges before the first spike and after the last spike. To resolvethis ambiguity we need to add auxiliary spikes. For the be-ginning of the spike train, we assign an auxiliary spike at themaximum of the distance between the start of the observationinterval and the first spike, and the first known ISI t ( n ) s aux = t ( n )1 − max { t ( n )1 − t s , t ( n )2 − t ( n )1 } . (A.1)This definition assumes that the rate stays the same at both sidesof the spike unless the edge is too far away for this to be true,in which case the auxiliary spike is assigned at the edge. Anal-ogously, the time of the auxiliary spike at the end is t ( n ) e aux = t ( n ) M + max { t e − t ( n ) M , t ( n ) M − t ( n ) M − } . (A.2)If the first or last spike is at the edge, no edge correction isnecessary at that end. This defines the ISI which is then usednot only for the ISI-distance but also for the A-SPIKE-distanceand A-SPIKE-synchronization.An auxiliary spike used for the edge e ff ect correction is ba-sically treated as any other spike, for example they can be thenearest neighbor to a real spike. But there is one exception: Inorder to avoid artificial synchrony at the edges in the A-SPIKE-distance, they use the distance to the nearest neighbor from thefirst / last real spike ∆ t ( n ) s aux = ∆ t ( n )1 and ∆ t ( n ) e aux = ∆ t ( n ) M n . (A.3) Appendix A.2. Edge e ff ect correction for A-SPIKE-synchronization For the A-SPIKE-synchronization profile we first apply theedge e ff ect correction described above and then calculate thecoincidence windows following Eqs. 14 and 15.For cases when there is a spike right at the edge, we use theone ISI that exists for setting the coincidence window of thespike to τ ( n )1 = x ( n )1 F and τ ( n ) M = x ( n ) MP . (A.4)We also determine that an auxiliary spike can under no circum-stance be part of a coincidence nor can it have a coincidencecounter. Finally, an auxiliary spike does not count as a spike inthe normalization. Appendix A.3. Special cases for A-ISI- and A-SPIKE-distance
Empty spike trains and spike trains with only one spike donot provide the ingredients needed to apply Eq. A.1 and A.2.In order to define the ISI of an empty spike train without anyspikes, we assign auxiliary spikes to its edges, the beginningand the end of the recording interval. This is the only intervalfor which we can guarantee that there were no spikes.However, while we can now use Eq. 3, Eq. A.3 for the dis-tance to the nearest neighbour of the auxiliary spikes is still14ll-defined, since there are no real spikes. In this case a value isassigned exactly as in Eq. 7 and the nearest neighbor can eitherbe a real or another auxiliary spike. A very reasonable impli-cation of this definition is that two empty spike trains will beconsidered equal by both measures.Similarly, it is not possible to assess the rate at either sideof a single spike. The most reasonable auxiliary spike locationis again at the edge of the recording. Thus for both cases, theauxiliary spikes are assigned at the edges as t ( n ) s aux = t s and t ( n ) e aux = t e (A.5)and this completes the definitions for the A-ISI- and the A-SPIKE-distance. Appendix A.4. Special cases for A-SPIKE-synchronization
For A-SPIKE-synchronization the situation is slightly di ff er-ent, since it is not continuous but only defined at the times ofthe spikes. This means that by definition an empty spike traincannot have synchronous spikes and thus has no value. In caseall spike trains are empty, we set A-SPIKE-synchronization to S AC =
1, i.e. empty spike trains are considered to be perfectlysynchronous. If a spike train contains only a single spike, weuse half the spike train length to define the coincidence windowfor the spike as τ ( n )1 =
12 ( t e − t s ) . (A.6)These special cases complete the definition of A-SPIKE-synchronization. Appendix B. Experimental recordings
Appendix B.1. Microelectrode array recordings from mousecortical cells
The electrophysiological data analyzed in Sections 2.5.1 and2.5.2 were recorded in the group of Prof. Jari Hyttinen at Tam-pere University of Technology / BioMediTech, Tampere, Fin-land. These recordings were performed prior to and indepen-dently from the design of this study.Between 5,000 and 25,0000 commercially available pri-mary mouse cortical cells (A15586, Gibco, Thermo Fisher)were plated on five microelectrode arrays (MEAs; four60MEA200 / / + ◦ C,5% CO2, 95% air) prior to and between recordings. Data wererecorded three times a week between the 4th and the 35th dayin vitro. Every recording lasted five minutes and was performedwith 25 kHz sampling rate. Spike detection was carried out bysetting an amplitude threshold at five times the standard devi-ation of the signal-noise level and the spike time stamps were stored with the Neuroshare Library for MATLAB (Multi Chan-nel Systems). We used two recordings for our examples andtest sets.The five real data spike trains used in the test set (spike trains26 to 30 in Fig. 5) were selected from these data by hand to rep-resent di ff erent time scales but chosen such that spike numberswere quite constant and comparable to the artificial examples. Appendix B.2. Patch clamp recordings of rat thalamocorticalrelay cells
The electrophysiological data analyzed in Section 2.5.4 wererecorded at the Swammerdam Institute for Life Sciences, Uni-versity of Amsterdam, the Netherlands. Again, these recordingswere performed prior to and independently from the design ofthis study. The experiments carried out on brain slices fromWistar rats (Harlan, Netherlands; postnatal days 12-16) wereapproved by the animal welfare committee of the University ofAmsterdam.For details on the animals, slice preparation and electrophysi-ological recordings, see Zeldenrust et al. (2013). In the current-clamp measurements the cell was injected with current that con-sisted of a DC component with superimposed noise: a computergenerated (MATLAB) time series of Gaussian distributed ran-dom numbers of a length of 300 s, filtered by an exponentialfilter with a time constant τ =
10 ms and a standard deviationof σ =
100 pA. A slow feedback system controlled the back-ground DC current to stabilize the membrane voltage at oneof the specified values (-80 mV, -70 mV, -60 mV or -50 mV)before the actual recording started; after the start this DC cur-rent component was fixed. The same frozen ( = an exactly re-produced computer generated) noise train was injected into thesoma of the TCR neuron for every repetition of the experiment.Signals were filtered at 510 kHz and sampled at 1020 kHz.The recordings consisted of trials from five di ff erent cells ofwhich only two included trials for all four levels of membranepotential. The cells analyzed were recorded three times andreliability was assessed by trial to trial variations. References
Andrzejak, R. G., Mormann, F., Kreuz, T., 2014. Detecting determinism frompoint processes. Physical Review E 90, 062906.Bear, M. F., Connors, B. W., Paradiso, M. A., 2007. Neuroscience: Exploringthe Brain. Lippincott Williams & Wilkins, Philadelphia, PA, USA.Chicharro, D., Kreuz, T., Andrzejak, R. G., 2011. What can spike train dis-tances tell us about the neural code? J Neurosci Methods 199, 146–165.Cutts, C. S., Eglen, S. J., 2014. Detecting pairwise correlations in spike trains:an objective comparison of methods and application to the study of retinalwaves. Journal of Neuroscience 34 (43), 14288–14303.Dura-Bernal, S., Li, K., Neymotin, S. A., Francis, J. T., Principe, J. C., Lytton,W. W., 2016. Restoring behavior via inverse neurocontroller in a lesionedcortical spiking model driving a virtual arm. Frontiers in Neuroscience 10.Espinal, A., Rostro-Gonzalez, H., Carpio, M., Guerra-Hernandez, E. I.,Ornelas-Rodriguez, M., Puga-Soberanes, H. J., Sotelo-Figuero, M. A.,Melin, P., 2016. Quadrupedal robot locomotion: A biologically inspired ap-proach and its hardware implementation. Computational Intelligence andNeuroscience, 5615618.Hales, C. M., Rolston, J. D., Potter, S. M., 2010. How to culture, record andstimulate neuronal networks on micro-electrode arrays (MEAs). Journal ofVisualized Experiments 39, e2056. zhikevich, E. M., Desai, N. S., Walcott, E. C., Hoppensteadt, F. C., 2003.Bursts as a unit of neural information: selective communication via reso-nance. Trends in Neurosciences 26, 161–167.Kapucu, F. E., Tanskanen, J. M. A., Mikkonen, J. E., Yl¨a-Outinen, L., Nark-ilahti, S., Hyttinen, J. A. K., 2012. Burst analysis tool for developing neu-ronal networks exhibiting highly varying action potential dynamics. Fron-tiers in Computational Neuroscience 6, 38.Kreuz, T., 2011. Measures of spike train synchrony. Scholarpedia 6 (10), 11934.Kreuz, T., Chicharro, D., Andrzejak, R. G., Haas, J. S., Abarbanel, H. D. I.,2009. Measuring multiple spike train synchrony. J Neurosci Methods 183,287.Kreuz, T., Chicharro, D., Greschner, M., Andrzejak, R. G., 2011. Time-resolvedand time-scale adaptive measures of spike train synchrony. J Neurosci Meth-ods 195, 92.Kreuz, T., Chicharro, D., Houghton, C., Andrzejak, R. G., Mormann, F., 2013.Monitoring spike train synchrony. J Neurophysiology 109, 1457.Kreuz, T., Haas, J. S., Morelli, A., Abarbanel, H. D. I., Politi, A., 2007. Mea-suring spike train synchrony. J Neurosci Methods 165, 151.Kreuz, T., Mulansky, M., Bozanic, N., 2015. SPIKY: A graphical user interfacefor monitoring spike train synchrony. J Neurophysiol 113, 3432.Kreuz, T., Satuvuori, E., Pofahl, M., Mulansky, M., 2017. Leaders and follow-ers: Quantifying consistency in spatio-temporal propagation patterns. NewJournal of Physics 19, 043028.Lyttle, D., Fellous, J. M., 2011. A new similarity measure for spike trains:Sensitivity to bursts and periods of inhibition. J Neurosci Methods 199, 296.Mulansky, M., Kreuz, T., 2016. Pyspike - A python library for analyzing spiketrain synchrony. Software X 5, 183.Naud, R., Gerhard, F., Mensi, S., Gerstner, W., 2011. Improved similarity mea-sures for small sets of spike trains. Neural Comput 23, 3016.Qu, J., Wang, R., Du, Y., Yan, C., 2016. Advances in Cognitive Neurodynamics(V). Springer Singapore, Ch. An Improved Method of Measuring MultipleSpike Train Synchrony, p. 777.Quian Quiroga, R., Kreuz, T., Grassberger, P., 2002. Event synchronization: Asimple and fast method to measure synchronicity and time delay patterns.Phys. Rev. E 66, 041904.Quian Quiroga, R., Panzeri, S., 2013. Principles of neural coding. CRC Taylorand Francis, Boca Raton, FL, USA.Rabinowitch, T. C., Knafo-Noam, A., 2015. Synchronous rhythmic interactionenhances children´s perceived similarity and closeness towards each other.PLoS ONE 10, e0120878.Rusu, C. V., Florian, R. V., 2014. A new class of metrics for spike trains. NeuralComputation 26 (2), 306–348.Sherman, S. M., 2001. Tonic and burst firing: Dual modes of thalamocorticalrelay. Trends in Neurosciences 24, 122–126.van Rossum, M. C. W., 2001. A novel spike distance. Neural Computation 13,751.Victor, J. D., 2015. Spike train distance. Encyclopedia of Computational Neu-roscience, 2808–2814.Victor, J. D., Purpura, K. P., 1996. Nature and precision of temporal coding invisual cortex: A metric-space analysis. J Neurophysiol 76, 1310.Zapata-Fonseca, L., Dotov, D., Fossion, R., Froese, T., 2016. Time-series analy-sis of embodied interaction: Movement variability and complexity matchingas dyadic properties. Frontiers in Psychology 7.Zeldenrust, F., Chameau, P. J. P., Wadman, W. J., 2013. Reliability of spike andburst firing in thalamocortical relay cells. Journal of Computational Neuro-science 5 (3), 317–334.URL http://dx.doi.org/10.1007/s10827-013-0454-8http://dx.doi.org/10.1007/s10827-013-0454-8