Email updates

Keep up to date with the latest news and content from Neural Systems & Circuits and BioMed Central.

Open Access Research

A novel, jitter-based method for detecting and measuring spike synchrony and quantifying temporal firing precision

Ariel Agmon

Author Affiliations

Department of Neurobiology and Anatomy and the Sensory Neuroscience Research Center, West Virginia University, Morgantown, WV 26506-9303, USA

Neural Systems & Circuits 2012, 2:5  doi:10.1186/2042-1001-2-5

The electronic version of this article is the complete one and can be found online at: http://www.neuralsystemsandcircuits.com/content/2/1/5


Received:9 November 2011
Accepted:5 March 2012
Published:2 May 2012

© 2012 Agmon; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Precise spike synchrony, at the millisecond or even sub-millisecond time scale, has been reported in different brain areas, but its neurobiological meaning and its underlying mechanisms remain unknown or controversial. Studying these questions is complicated by the lack of a validated, well-normalized and robust index for quantifying synchrony. Previously used measures of synchrony are often improperly normalized and thereby are not comparable between different experimental conditions, are sensitive to variations in firing rate or to the firing rate differential between the two neurons, and/or rely on untenable assumptions of firing rate stationarity and Poisson statistics. I describe here a novel measure, the Jitter-Based Synchrony Index (JBSI), that overcomes these issues.

Results and discussion

The JBSI method is based on the introduction of virtual spike jitter. While previous implementations of the jitter method used it only to detect synchrony, the JBSI method also quantifies synchrony. Previous implementations of the jitter method used computationally intensive Monte Carlo simulations to generate surrogate spike trains, whereas the JBSI is computed analytically. The JBSI method does not assume any specific firing model, and does not require that the spike trains be locked to a repeating external stimulus. The JBSI can assume values from 1 (maximal possible synchrony) to −1 (minimal possible synchrony) and is therefore properly normalized. Using simulated Poisson spike trains with introduced controlled spike coincidences, I demonstrate that the JBSI is a linear measure of the spike coincidence rate, is independent of the mean firing frequency or the firing frequency differential between the two neurons, and is not sensitive to co-modulations in the firing rates of the two neurons. In contrast, several commonly used synchrony indices fail under one or more of these scenarios. I also demonstrate how the JBSI can be used to estimate the spike timing precision in the system.

Conclusions

The JBSI is a conceptually simple and computationally efficient method that can be used to compute the statistical significance of firing synchrony, to quantify synchrony as a well-normalized index, and to estimate the degree of temporal precision in the system.

Background

The most basic temporal relationship in the active nervous system is firing synchrony, that is, the coincidence of spikes fired by two or more neurons. Firing synchrony has been widely reported in the nervous system and is of much interest to neuroscientists both because of its potential role in encoding, transmitting and decoding of information in the brain, and because it may reveal something about the underlying synaptic (or other) interactions that generate it (for reviews, see [1,2]). Firing synchrony increases the reliability of firing in downstream neurons that receive the coincident synaptic inputs [3-5] and may therefore be an efficient way to transmit information through multiple layers of feed-forward excitatory networks [6-8]. Alternatively or in addition, synchrony may be used as a ‘binding’ signal [9] (but see [10]); may be the neural correlate of attention or expectation [11-13]; or may encode information about the stimulus, beyond that available from the firing rate alone [14,15]. Synchrony can arise when two or more neurons respond to a sudden onset or rapid modulation of a sensory stimulus, such as an auditory click [16], a visual saccade [17] or a high-velocity whisker deflection [18-20]. The focus here, however, will be on stimulus-independent synchrony that is observed (in some systems) after correcting for stimulus effects, or in the absence of a stimulus [5,21,22]. Whereas stimulus-dependent synchrony may encode information about the stimulus that drives it, stimulus-independent synchrony may teach the experimenter something about the connectivity of the underlying network. Two network motifs have most often been invoked to account for stimulus-independent synchrony: shared excitatory inputs from diverging presynaptic axons [23,24] and electrical coupling by gap junctions [25-27]. In a recent publication, we demonstrated that a third motif, reciprocal inhibitory chemical synapses, can drive stimulus-independent synchrony with sub-millisecond precision, in the absence of shared inputs or electrical coupling [28].

Since spike coincidence could also be a purely chance occurrence, an important first step in studying synchrony is to determine if the observed coincidence level is ‘real’, meaning statistically significant - in other words, to calculate its p-value or Z-score. We refer to this as detecting synchrony. However, if one wants to decipher the synaptic connectivity driving the synchrony, one needs to go beyond a p-value and to quantify the synchrony by a well-normalized ‘synchrony index’, which can be examined for correlation with synaptic parameters such as the strength of the electrical coupling or the rate of shared inputs. Various methods for achieving the twin goals of detecting and quantifying synchrony, some more practical than others, have been proposed and used in the past, and some of them are described and examined in detail in a recently published compendium [29]. Unfortunately, many previously proposed methods have serious limitations when used for detection and quantification of synchrony in real spike trains, as we will see below. Analysis techniques that fall under the general class of ‘jitter methods’ [30] seem to overcome many of the drawbacks of other methods when it comes to detecting synchrony; however, jitter methods have not been used previously to quantify synchrony. Here, I take jitter methods one step further and use them to compute a novel synchrony measure, the Jitter-Based Synchrony Index (JBSI).

In the remainder of this background section, I first show how synchrony can be detected when the two spike trains follow stationary Poisson statistics. I then demonstrate numerically how this detection method fails when firing rates are not stationary and explain how jitter methods overcome this hurdle. Finally, I define and describe two of the synchrony indices commonly employed by experimental neurophysiologists to quantify synchrony. In the Results, I first propose a set of five requirements that an ideal index of synchrony should fulfill, pointing out where existing measures fail to meet these requirements. I then explain how the JBSI is computed and follow with a comparison between the JBSI and previously used indices. Lastly, I demonstrate how the JBSI can be used to estimate the precision of spike timing in the system. Computational details are provided in three appendices.

Detecting synchrony: the Poisson case

Assume that we have recorded simultaneous spike trains from two neurons. I will refer to the slower-firing neuron as Neuron 1 or the reference neuron, and to the faster-firing neuron as Neuron 2 or the target neuron. Let us denote the number of spikes in the reference and target trains by n1 and n2 and the average firing rates by r1 and r2, respectively (so n1 ≤ n2 and r1 ≤ r2).

On the face of it, detecting and measuring synchrony seems quite simple. First, decide on a definition of ‘synchrony’: how far apart two spikes can occur and still be considered as synchronous. Throughout this paper, I will use τS to represent this synchrony span. Second, count how many of the spikes fired by the reference neuron occurred within the predefined synchrony span from any spike in the target train. I will denote this observed coincidence count as NC. Third, normalize the coincidence count so it can be used to compare paired spike trains of different lengths. A simple way to do this is to divide it by the total number of reference spikes, transforming the coincidence count to a coincidence rate RC ≡ NC/n1 (note that rate is used here in the sense of per spike rather than per unit time). The advantage of dividing by n1 rather than, say, n2 or the average of n1 and n2, is that RC can take values over the full range between 0 (no coincidences) to 1 (all spikes coincident).

On their own, however, NC and RC are not very informative, because for any two overlapping series of events there is some non-zero probability that occasionally they will coincide in time, purely by chance. In the case of two spike trains, we want to know whether the two neurons fired independently and the observed coincidences were therefore chance occurrences; this is the null hypothesis. The alternative is that the neurons were coupled, that is, interacted in some way with each other or with a third neuron, causing them to synchronize more (or perhaps less) than expected by chance. Our first task is therefore to determine if the observed coincidence was statistically significant. To do so, we need to know something about the distribution of all chance coincidence counts; if this distribution was known, we could determine the fractional area beyond NC under the tails of the distribution, in other words, the p-value, which would indicate how likely NC was to belong to this set. At the very least, we would like to know the mean of this distribution - that is, the expected coincidence count NC› - and its variance, which would then allow us to calculate a Z-score (distance of the observed value from the mean, in units of standard deviation). If the distribution is reasonably close to normal, the Z-score can be directly converted to a p-value (for example, p = 0.05 corresponds to Z = ~2 and p = 0.001 to Z = 3.3).

One way to estimate the distribution of chance coincidences is to assume a specific statistical firing model, and very often the assumption made is that the spike trains are a stationary Poisson process. Under this assumption, the expected number of spikes occurring within an interval ∆t is t, where r is the time-independent (stationary) average firing rate. Going back to our paired spike trains, we observe that, for a sufficiently small synchrony span τS, the probability that a reference spike will occur within ±τS from any given target spike is equal to the average number of reference spikes in a window of width 2τS, or 2τS·r1. Because the target neuron fired n2 spikes during the recording epoch, we should expect a total of <a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M1">View MathML</a> spike coincidences. For a recording epoch of length T, we have r1= n1/T and r2= n2/T, so we can re-write the expected number of coincidences either as:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M2">View MathML</a>

(1)

or as:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M3">View MathML</a>

(2)

In a Poisson process, the variance of the observed counts is equal to the expected counts: we can therefore test the null hypothesis of independent firing by calculating the likelihood that NC, the observed number of coincidences, belongs to a distribution with a mean and variance of NC.

Detecting synchrony: when Poisson cannot be assumed

The pitfall in the procedure described above is that, in actuality, we are testing a dual null hypothesis: that (a) the neurons fired independently and that (b) each spike train is a Poisson process with a stationary average firing rate. If (b) happens to be false, we stand the risk of rejecting the null hypothesis, even when (a) is true. We can illustrate this with a simple numerical example.

Assume that both neurons fired 10,000 spikes each during a 250 s recording epoch, so r1=r2 = 40 Hz. If τS = 0.5 ms, then the number of coincidences expected by chance is NC = 0.001 × 40 × 40 × 250 = 400, with a standard deviation of 20. Now say that we actually observed 500 coincidences. Because this is five standard deviations above the expected mean, we conclude that this degree of precise synchrony was highly unlikely to be a random deviation (less than one in a million probability) and that the two neurons must have been coupled. Unbeknownst to us, however, both neurons received bursts of inhibitory inputs from a common subset of inhibitory neurons during the recording epoch. These bursts occurred randomly, but were powerful enough to silence both neurons (concurrently) for 20% of the recording epoch. The actual firing epoch (that is, the epoch during which the neurons were ‘allowed’ to fire) was therefore only 200 s, so the firing rate of each neuron was actually 10,000/200 = 50 Hz, and the true expected number of spike coincidences was NC = 0.001 × 50 × 50 × 200 = 500. In other words, the observed synchrony was precisely at chance level. While one could argue that the two neurons were, in a sense, coupled by the common inhibitory input, and thus violated both assumptions (a) and (b), it would clearly be erroneous to conclude that this common inhibitory input caused the neurons to synchronize with a precision of ±0.5 ms.

The example above underscores the fact that it is the ‘local’ firing rate, rather than the global average firing rate, that needs to be taken into account when determining the expected coincidences, since global averaging ignores modulations in firing rates. If these bursts of inhibitory inputs were all identical or nearly so, if they repeated at regular intervals and if we knew their times of occurrence, we could estimate how many coincidences were introduced because of these co-modulations by dividing the full record into segments aligned on the beginning (or end) of each inhibitory burst and exchanging the spike trains of one neuron between segments. The spike coincidences remaining after this procedure could be taken as a good estimate of NC and can then be subtracted from the observed coincidence count NC to yield an estimate of the ‘excess coincidence’, that is, the coincidence above chance level. This is the same ‘shuffling’ procedure commonly used by electrophysiologists to account for spurious coincidences introduced by a sensory stimulus [31]. However, unlike co-modulations resulting from an experimentally administered stimulus, co-modulations resulting from uncontrolled environmental or systemic influences are rarely reproducible, do not occur at regular intervals and are not locked to the stimulus (if there is a stimulus). Indeed, they may not even be recognized by the experimenter. Therefore, these co-modulations cannot be corrected by shuffling.

Another potential solution is to divide the firing epoch into smaller segments and calculate the local firing rates, and from them the expected coincidences, separately for each segment. This would allow us to correct for any modulations in firing rates that happen on a time scale slower (longer) than the width of each segment. Further scrutiny, however, reveals multiple problems with this approach. First, where should the boundaries between segments be placed? Even a small segment may happen to straddle two epochs with very different firing rates, and may need to be further subdivided. Second, how local is ‘local’ - how small should each segment be? To account for as many co-modulations as possible, some of which may have fast time courses, we would need to make the segments as small as possible; but with very small segments, firing rate estimation would become very imprecise. Ideally, we would like to know the instantaneous firing rate in the immediate vicinity of each spike, but this is impossible to determine from a non-repeating spike train.

The jitter method offers a solution to this dilemma. It allows us to estimate the coincidence count expected from the local firing rate, without explicitly determining the rate itself. To do so, we replace one of the two spike trains with a virtual one (a ‘surrogate’), in which each spike is slightly shifted (or ‘jittered’) from its original time of occurrence by a random amount, within a predetermined ‘jitter window’ of ±τJ. This procedure destroys any spike coincidences that may have resulted from interactions on a time scale faster than 2τJ, but fully preserves local firing rates (and therefore any spike coincidences attributable to co-modulations in these rates) at all slower time scales. Again, we can take the number of jitter-resistant coincidences as an estimator of ‹NC›, the expected number of coincidences, and use it to estimate the excess synchrony. Moreover, if we can determine the distribution of all possible coincidence counts likely to be observed after a jitter, we can determine if the experimentally observed count NC is an exceptional value (that is, falls in the tails of this distribution) and therefore reflects true, statistically significant synchrony, or if it is likely to be a chance occurrence. Note that the jitter procedure relies on the assumption that co-modulations in firing rates occur on a relatively slow time scale compared with the time scale of the spike coincidences we are interested in. Thus, we need to set τJ judiciously: if we set τJ too large, we risk destroying some of the chance coincidences caused by co-modulations, and we will then underestimate ‹NC› and overestimate the true degree of synchrony, just as in the numerical example above. Conversely, it would be meaningless to make τJ smaller than the synchrony span τS, because such a small jitter would likely preserve all coincidences, both chance and real, and not help us in distinguishing between the two.

The rather intuitive idea of introducing virtual spike jitter was rigorously formulated and explored theoretically and in computer simulations by Geman and colleagues [30,32-34], and applied to experimental data by them [32,35] and others [36-40], in some cases critically [41]. The approach taken here differs in two important ways from these previous implementations. First, it differs in how it determines the distribution of all possible coincidence counts. Previous implementations used the Monte Carlo approach: to generate many (hundreds or thousands of) realizations of jittered spike trains and count the number of spike coincidences for each one, which is computationally intensive. The current approach is based on an analytical computation of the exact probability distribution, thus considerably reducing the computational effort. Even more importantly, previous implementations used the jitter method only to detect synchrony, not to quantify it. Here, I use the jitter method to define a novel synchrony index, the JBSI, which can be used to quantify synchrony and (as we will see) is robust under a wide variety of realistic conditions under which other synchrony indices fail.

Quantifying synchrony: how to compare different cell pairs

Detecting synchrony is only the first task facing us. Our second task is to quantify the synchrony to allow valid comparisons of synchrony strength between different cell pairs, even when recorded in different experiments or even by different investigators. At first glance, it may seem that the very measures of statistical significance could also be used to quantify the strength of the synchrony. However, measures of statistical significance can be made arbitrarily high (Z-scores) or arbitrarily close to zero (p-values) by increasing the sample size, that is, by using longer spike trains, even if the rate of synchrony (per spike and per time) remains unchanged. Measures of significance cannot, therefore, be used directly to quantify synchrony - to do so we need a synchrony index.

Many different synchrony indices have been proposed in the literature but relatively few of these have gained popularity with the experimental neurophysiology community, and there seems to be no generally agreed upon ‘best’ index [5,12-14,39]. Neurophysiologists often depict the outcome of paired recordings graphically, by a cross-correlation histogram (cross-correlogram, CCG) [22,31], and most previously employed indices were calculated from the CCG. If the bin width of the raw CCG is chosen to be 2τS, then the height of the central bin will be equal to NC. To generate a synchrony index, the coincidence count is first corrected for chance coincidences by subtracting the expected from the observed count; this yields an estimate of excess coincidences. If the spike trains are locked to a repeating stimulus and if repeated spike trains are assumed to be stationary and reproducible, then this correction can be done by subtracting a shuffled cross-correlogram [31]. Otherwise, it is done by subtracting the average number of counts in a region of the CCG equal in width to but away from the central peak. For simplicity, we again assume that the central peak of the CCG is one-bin wide, so that the average count in a far-from-center bin is very close to the overall average count per bin, b·n1·n2/T, where b is the bin width. With a bin width of 2τS , the average count is therefore equal to our previously defined (see equation 2):

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M4">View MathML</a>

Thus, the excess coincidence count is equal to NCNC.

Since the excess coincidence count depends on the number of recorded spikes, it is typically normalized in some manner to yield a synchrony index that is not dependent on spike number. For comparison with the JBSI, I selected two of the more commonly used indices in the experimental literature. In the first index, which is usually notated in the literature by E but will be called here Excess Coincidence Index or ECI, the excess coincidence count is simply normalized by the number of spikes in the reference train [42-44]:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M5">View MathML</a>

(3)

In the second index, referred to in the literature as the cross-correlation coefficient (CCC) [4,45] or, somewhat loosely, as the correlation coefficient [19,46,47], the excess count is divided by a more complicated expression:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M6">View MathML</a>

(4)

This normalization factor is justified mathematically in [48]; for completeness, I derive the CCC from first principles in Appendix C.

Given the availability of these and other synchrony indices in the neurophysiological literature, the reader may wonder: how should an experimenter decide which is the best index to use, and why did the author bother to devise yet one more index? To answer these two questions, I first lay out a set of requirements that, I propose, should be met by an ideal synchrony index. I then show how both the ECI and the CCC fail to meet some of these requirements whereas a novel index, the JBSI, fulfills them all.

Results and discussion

What should we expect from an ideal synchrony index?

I propose the following as a minimum set of requirements to be met by an ideal measure of synchrony:

1. The measure should be applicable to both periodic and aperiodic data, and to spontaneous spike trains as well as to stimulated responses. This eliminates from the current discussion a large collection of methods for analyzing periodic signals and measures for quantifying synchrony within the frequency domain (for example, phase coherence) or for quantifying the degree of phase-locking to a precisely repeating stimulus (for example, vector strength) (reviewed in [49-51]).

2. As indicated above, the index should be properly normalized, to allow valid comparison between different experiments. This means that ‘perfect synchrony’ and ‘purely chance synchrony’ should assume the same finite values (say, 1 and 0, respectively), regardless of the experimental details. Several synchrony indices previously used in the literature are un-normalized or improperly so. For example, the synchrony index designated in the literature as k’ [23,52-54] and defined as the ratio NC/‹NC, is obviously unbounded, because it can assume very large values if NC (the expected coincidence count) is very low.

3. The synchrony measure should reflect the intrinsic strength of the network motif that drives it (common inputs, electrical coupling or mutual inhibition) but should not be sensitive to the mean firing rate, because the latter depends mostly on factors extrinsic to the two neurons such as background synaptic input or experimentally injected current. As previously demonstrated [53,55-57] and as is confirmed in Results, the ECI and several other previously used indices exhibit a strong negative dependence on the firing rate and therefore do not fulfill this requirement.

4. The ideal measure should be maximized whenever spikes fired by Neuron 1 are precisely synchronized with spikes of Neuron 2, even if Neuron 2 fires many more spikes than Neuron 1. As shown in Results and in Appendix C, the CCC is highly sensitive to the firing rate differential between the two trains. Similarly, so-called spike train metrics [58-61] are sensitive to differences between spike times as well as differences between spike numbers and therefore do not meet this requirement.

5. The method should be applicable to any arbitrary pair of concurrent spike trains, without assuming any specific firing statistics or firing rate stationarity. Many published methods for detecting and quantifying synchrony assume, explicitly or implicitly, that (in the case of spontaneous firing) spike trains are stationary Poisson processes (for example, [62]), or, in the case of evoked responses, that responses to repeating stimuli are reproducible (for example, [31,63]). However, real spike trains are, in general, not stationary Poisson processes, as the firing rate is continuously modulated by factors beyond the experimenter’s control or even knowledge. Similarly, repeating stimuli may not generate reproducible responses due to trial-to-trial variation in subject motivation, state of anesthesia, viability of the preparation and so on. These fluctuations will often manifest themselves as modulations in firing rate and/or in latencies, common to both neurons. As previously noted in the literature [64-66] and as illustrated numerically above (see Background), these co-modulations can lead to erroneous detection of correlations or synchrony. I show below that both the ECI and the CCC erroneously indicate synchrony between independently generated spike trains when co-modulations of firing rate are introduced.

Calculation of the Jitter-Based Synchrony Index

Let us assume that we have recorded two simultaneous spike trains. We express each spike train as a vector of spike occurrence times, so we have two vectors: one from the reference neuron, containing n1 spikes designated by their time of occurrence t10t1it1n1-1, and one from the target neuron, containing n2 spikes designated t20t2kt2n2-1. We will assume that the time of occurrence of each spike is known to any desirable precision (say 0.1 ms). In Figure 1, reference spikes are depicted in blue and target spikes are depicted in red. The coincidence counting process can be represented graphically by drawing a synchrony window WSk of width 2τS centered on each target spike t2k (Figure 1B) and counting how many blue spikes fall within a synchrony window. Formally, we assign to each reference spike t1i a binary value Syn(i), defined as:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M7">View MathML</a>

(5)

and then define:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M8">View MathML</a>

(6)

Obviously, NC will depend on our choice of the synchrony span τS. If τS is increased beyond ½ of the smallest interspike interval in the target train, synchrony windows will begin to overlap, and a given spike in the reference train may be synchronous with more than one spike in the target train. According to (5), it will still only be included once in the coincidence count, thus guaranteeing that NC ≤ n1 or that RC ≤ 1. Also, note that our definition of coincidence does not rely on time binning and thereby avoids the pitfall of two near-coincident spikes falling into two adjacent bins and not being counted as synchronous.

thumbnailFigure 1. Using virtual spike jitter to quantify synchrony. (A) A hypothetical segment from two simultaneously recorded spike trains, with spikes from the reference and target train represented as blue and red bars, respectively. (B) By placing a synchrony window (red) of ±τS centered on each target spike, we see that three of the five reference spikes are synchronous. (C) A jitter window (blue) of ±τJ is centered on each reference spike. (D) Areas of overlap between jitter and synchrony windows are shaded; spikes are omitted for clarity. The probability pi that any given reference spike will be synchronous after jittering is given by the shaded fraction of its jitter window, and is indicated below the window. The probabilities are then used to determine the expected number of spike coincidences (∑pi), which is then subtracted from the observed number to yield an estimate of excess coincidences. The latter is normalized by the number of spikes in the reference train and multiplied by a scaling factor to yield the JBSI.

Next, we select a jitter span τJ, τJ > τS, and shift (jitter) each reference spike t1i by up to ±τJ. (Note that we jitter only the shorter, reference spike train, thereby conserving computational effort.) As shown below, it is advantageous to choose the ratio τJS to be 2. The number of coincidences observed after the above jitter procedure, NCJ, is a random variable with some probabilistic distribution; according to our null hypothesis, NC should be from the same distribution. To test this hypothesis, we need to calculate the probability of observing any given number of total spike coincidences N, 0 ≤ N ≤ n1, after a jitter. We call this probability PJ(N):

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M9">View MathML</a>

(7)

To compute PJ(N), we first calculate the probability pi that spike t1i will be synchronous with spike t2k, for at least one k, after applying the jitter. The process can be represented graphically (Figure 1C) by drawing a jitter window WJi of width 2τJ, centered on each reference spike t1i (blue). It then becomes apparent that pi is equal to the fraction of WJi that intersects (overlaps) the union set of all synchrony windows (the intersections are shaded in Figure 1D).

Formally:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M10">View MathML</a>

(8)

In Appendix A1 we compute the union set of all synchrony windows, and in Appendix A2 we compute pi. PJ(N) can then be computed exactly from the vector pi using an efficient recursive algorithm, proposed in [67] and provided for the reader’s convenience in Appendix A3. This calculation demonstrates that the distribution PJ(N) converges rapidly to a normal distribution with the same mean and standard deviation (not shown). We therefore do not need to compute PJ(N) explicitly to test our null hypothesis - we can use the fact that the fractional area under the tail of a normal distribution, that is, the p-value, can be determined directly from the standardized distance of the tail from the mean, the Z-score. To calculate the Z-score, we only need to know the expected value and the variance of NCJ, and these can be derived directly from the pi’s:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M11">View MathML</a>

(9)

and

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M12">View MathML</a>

(10)

(In the case of equal probabilities, pi= p for all i’s, and we get the well-known formulas for the mean and variance of a binomial distribution, Np and Np·(1-p), respectively.) The Z-score of the experimentally observed synchrony NC is therefore:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M13">View MathML</a>

(11)

The Z-score tells us if the coincidence count NC is a statistically significant observation. In other words, it allows us to detect synchrony. However, the Z-score cannot be used to quantify synchrony because, like the p-value and other measures of statistical significance, it depends on the sample size (the length of the spike train) and therefore is not directly comparable between different experiments. To quantify synchrony in a manner that would allow comparison between experiments, we define a normalized synchrony index - the JBSI:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M14">View MathML</a>

(12)

where β = 2 if τJS ≤ 2 and β = τJ/(τJ − τS) if τJS > 2.

As shown in Appendix B, the scaling factor β assures that the JBSI will attain its maximal value of 1 for the case of perfect synchrony. In the case of perfect asynchrony, however, the JBSI will attain its minimal value of –1 only if τJS ≤ 2 (see Appendix B). In all the simulations below, we select τJS= 2 and therefore β = 2.

Comparison of the Jitter-Based Synchrony Index with cross-correlogram-based synchrony indices

To test the JBSI and compare it with previously used synchrony indices, I generated simulated paired spike trains that (initially) followed Poisson statistics, by allowing each neuron to fire at random at an average firing rate that could be made time-dependent. Controlled spike coincidences were inserted into each train pair by randomly selecting reference spikes, at a probability D, to be shifted to within ±C of the next nearest target spike. Thus, the parameter D determined the injected coincidence rate (per spike), and the parameter C determined the precision of the synchrony (see Methods for details). For the simulations shown in Figure 2, C was maintained at 1 ms. The first second of firing from representative simulations is illustrated in the two rightmost columns of Figure 2, with reference and target spikes shown blue and red, respectively. Below each pair of spike trains is the CCG computed from the full train (about 1,000 spikes per neuron) using 2 ms-wide bins. The counts in each bin are normalized by the number of spikes in the reference train, so the height of the central peak is numerically equal to the coincidence rate RC.

thumbnailFigure 2. Comparison of the Jitter-Based Synchrony Index with the Excess Count Index, the cross-correlation coefficient and the corrected Excess Count Index . In the left column of each panel, the indices and their linear regression lines are plotted for a range of simulations in which one parameter of the simulation was varied. (A) The rate of inserted coincidences, D, was varied; (B) the average firing rate, r (defined as the geometrical average of the mean firing rate of the two neurons) was varied; (C) the differential in firing rates, ∆r, was varied; (D) the amplitude of firing rate co-modulation, M, was varied. In the two rightmost columns, 1 s segments of simulated spike trains, representing two extreme values of the varied parameter, are shown above the cross-correlograms (CCGs) computed from the same trains. CCGs were computed in 2 ms bins and the counts normalized by the number of reference spikes. The black regression lines in B and C represent the corrected Excess Count Index (ECIcor); for clarity, the ECIcor data points are not plotted. In A and D, the ECIcor data points precisely overlapped the CCC. Note that the JBSI was the only index that was robust against variations in all tested parameters.

I first compared the JBSI with two commonly used synchrony indices, the ECI and the CCC (see Background), by generating simulated spike trains in which the average firing rate in both neurons was maintained constant at about 70 Hz but the rate of injected coincidences D was increased parametrically, from 0 to 0.6. For each value of D, five runs of the simulation were generated. The three indices and their linear trends are plotted in Figure 2A, left panel. Clearly, all three indices increased more or less linearly with the rate of inserted coincidences, albeit with somewhat different slopes. Thus, under routine conditions, all three indices performed comparably well.

Next, I tested the indices under three challenges. The first challenge was a series of paired spike trains in which D was kept constant (at 0.25) while the firing rate in both neurons was parametrically increased from approximately 10 to 140 Hz (Figure 2B). According to requirement three of the set of requirements outlined at the beginning of this section, an ideal synchrony index should be independent of the firing rate. As seen in the left panel of Figure 2B, both the JBSI and the CCC correctly maintained a nearly constant value; in contrast, the ECI trend line dropped steeply with the increased firing rate, incorrectly implying loss of synchrony at the higher firing rates.

Although several previous studies note the negative dependence of the ECI (and other synchrony indices) on the firing rate [53,55-57], they differ in their explanation for this dependency and no remedies are proposed. To see why the ECI is dependent on the firing rate and how this flaw can be corrected, we note that the ECI is meant to estimate the true spike coincidence rate, RCtrue, by subtracting the chance spike coincidence rate RCchance from the total spike coincidence rate RC. However, the set of chance coincidences and the set of true coincidences are not mutually exclusive sets - their intersection is a set of coincidences that occur at a rate that equals the product of their individual rates, RCtrueRCchance. Therefore:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M15">View MathML</a>

(13)

This can be solved for RCtrue:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M16">View MathML</a>

(14)

We should therefore correct the definition of the ECI by dividing it by (1- RCchance). Since RCchance = NC›/n1, the corrected ECI, or ECIcor, is:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M17">View MathML</a>

(15)

Note that in the case of the JBSI such a correction is not needed, because chance spike coincidences are defined as spike coincidences that survive the jitter process, while true spike coincidences are those that are destroyed by the jitter process, and these two sets are mutually exclusive. As shown in Figure 2B (black trend line), the ECIcor was numerically very close to the JBSI and was not dependent on the firing rate.

The second challenge tested the fourth of the requirements outlined at the beginning of this section, that the index be independent of any firing rate differential between the two neurons. This was tested by simulated spike trains in which the injected coincidence rate D was maintained at 0.2, the geometric mean of the two average firing rates (r) was kept at approximately 45 Hz, but the difference between the two firing rates (∆r) was parametrically increased from 2.5 to 110 Hz (Figure 2C). Both the JBSI and the ECI correctly maintained a nearly constant value, and so did the ECIcor (black trend line), but the CCC decreased steeply with ∆r.

The dependence of the CCC on the firing-rate differential is a direct result of its definition. As shown in Appendix C, for a fixed number of spikes per train, the value of the CCC for the case of perfect synchrony, CCCmax, will not attain 1 unless the number of spikes in the two trains is equal; otherwise, CCCmax will be <1 and will decrease with an increasing firing rate differential. This suggests a simple remedy - divide the CCC by CCCmax to yield a corrected CCC index, the CCCcor. Interestingly, as shown in Appendix C, the CCCcor turns out to be identical to the ECIcor above. Thus, the ECIcor (or CCCcor) satisfies both requirements three and four, and is therefore preferable to either the ECI or the CCC. Indeed, when firing rates are stationary, the ECIcor may be the index of choice because it has all the advantages of the JBSI but is much easier to compute.

Unfortunately, the ECIcor fails, as do the ECI and the CCC, in regards to requirement five, in that these indices are sensitive to co-modulations in the firing rates. This is illustrated in Figure 2D, which represents simulations with independent spike trains (that is, D = 0) in which the total number of spikes was maintained nearly unchanged at about 1,000 per train but the mean firing rate of both neurons was co-modulated in time with a period of 0.5 s, with the amplitude of the modulating waveform varied parametrically. Even though no coincidences were inserted into these trains, the values of the ECI and the CCC were > 0 and increased with the modulation amplitude (note that the ECIcor is not plotted in Figure 2D, because in these simulations the firing rate differential was very low, so the ECIcor was virtually identical to the CCC). This is a fatal, non-remediable flaw in these commonly used synchrony indices, originating in the basic assumption that the two spike trains are a Poisson process, an assumption that is violated (as we saw in Background) when co-modulations are present. In contrast, the JBSI correctly reported no synchrony for these independent spike trains, even when co-modulations were present.

Comparison of the Jitter-Based Synchrony Index with the Jitter-Sensitive Synchrony Index

In a previous experimental publication [28] we introduced the Jitter-Sensitive Synchrony Index (JSSI) and used it to quantify sub-millisecond firing synchrony between inhibitory cortical interneurons. The JSSI was defined as the Z-score of the observed coincidence count (Equation 11) normalized by ((τJS-1)·n1)1/2. Like the JBSI, the JSSI is robust in regards to co-modulations in firing rates; however, the JSSI exhibits a negative dependency on the overall firing rate and therefore the JBSI is preferable to the JSSI.

Estimating the temporal precision of neuronal firing

How precise is neuronal firing in the brain, and how does one measure this precision? These questions have long occupied both experimental and theoretical neuroscientists, and do not yet have a satisfactory answer [7,68-71]. Jitter methods are particularly amenable for testing hypotheses regarding temporal precision. Jitter methods test a specific null hypothesis: that the temporal precision of firing is no better than ±τJ, for any desired τJ, and therefore a jitter of up to ±τJ should not reduce the observed coincidence count [30,32-34]. If the probability that the null hypothesis is true falls below any predetermined threshold, one can conclude that firing precision was better than ±τJ. One can then proceed to test increasingly smaller values of τJ until the null hypothesis can no longer be rejected; the smallest value of τJ that allows rejection of the null hypothesis can be considered an estimate of the firing precision in the system. This is illustrated in Figure 3A, in which three simulated paired spike trains were used to calculate Z-scores and JBSI for τJ values increasing from 1 to 16 ms in multiples of √2, while maintaining τJS = 2. The three superimposed plots correspond to simulations with the same rate of inserted coincidences (determined by the parameter D, which for these simulations was fixed at 0.2) but with different degrees of precision of synchrony (determined by the parameter C, which was varied between 1, 2 and 4 ms, as indicated in the figure legend). Each plot is an average of five different runs of the simulation. The black arrows in Figure 3A point to the intersection of the graphs with the line Z = 3.3, corresponding to p = 0.001. With this stringent threshold for significance, firing precision was better than about 1.5, 3 and 6 ms, respectively, for the three simulations, consistent with the corresponding precision parameter C used to generate each simulation.

thumbnailFigure 3. Estimating temporal precision of firing using the Jitter-Based Synchrony Index. The three superimposed plots represent different simulated paired spike trains; each plot is an average of five runs of the simulation (error bars indicate standard error of the mean). The average firing rate r (approximately 45 Hz) and the rate of inserted coincidences D (0.2) were identical for all pairs, but the precision of synchrony C was varied from 1 to 4 ms. (A) The Z-score is plotted against the jitter span τJ, which was varied from 1 to 16 in √2-fold increments while τS was maintained at τJ/2. The intersections of the three plots with the line Z = 3.3, representing a significance threshold of p = 0.001, correspond to τJ = 1.5, 3 and 6 ms (black arrows). (B) For the same simulations as in (A), the JBSI is plotted for increasing values of the synchrony span τS. Cutoff points, from which the JBSI fell steeply to the left, correspond to τS = 1, 2 and 4 ms (black arrows).

The problem with this approach is the arbitrariness of any selected significance threshold and the fact (already alluded to in this article) that one can increase statistical significance at will, simply by recording longer spike trains. To remove the dependence on the length of the train, one can use the JBSI instead of the Z-score, as is done in Figure 3B; however, the question then becomes what JBSI value one should use as a threshold. The solution is to determine firing precision based on the shape of the JBSI versus τS curve (Figure 3B) rather than on any particular threshold. As illustrated in Figure 3B, this curve is strongly asymmetric, sloping shallowly from its peak value rightwards, but sloping steeply to the left. It is instructive to examine the reason for these two different slopes. The shallow rightward decrease in JBSI values reflects the gradual increase in the widths of the synchrony windows WS: with the union of all synchrony windows occupying a larger fraction of the spike record, a jittered spike is more likely to fall into a synchrony window and thereby preserve or increase the expected coincidence count ‹NCJ›. The steep decline to the left results from a very different reason - it reflects the jitter parameter τJ falling to values that are too small to make any difference in the temporal structure of the spike train; in other words, τJ falling below the temporal precision of the system. This suggests that one should use the points at which the curves drop off steeply to the left as cutoff points for determining firing precision. The cutoff points indicated by arrows in Figure 3B correspond to τS values of 1, 2 and 4 ms, exactly the values of the corresponding precision parameter C.

Other applications of the Jitter-Based Synchrony Index method

In this manuscript I demonstrate the utility of the JBSI method for quantifying precise firing synchrony between two neurons; in principle, however, this method can be extended to other temporal relationships between spike trains. For example, instead of testing for synchrony (zero lag between the two spikes), one can test for a recurring lag LL ≠ 0, between the two spike trains, such as postulated by the ‘synfire chain’ hypothesis [46]. This is equivalent to shifting the red synchrony windows in Figure 1B by L. One can also use the JBSI to measure precision of firing between recurring trials recorded from a single neuron in response to a repetitive stimulus, by regarding the stimulus as the reference train. To generate a time-resolved index of synchrony, one can calculate the JBSI over a sliding window of any chosen width, provided that the firing rate is high enough to generate sufficient spikes within each window for a reliable computation. Finally, although a bivariate measure in its present form, the JBSI can be extended to synchrony between multiple simultaneously recorded spike trains, for example by averaging over all possible pairs. Future work may generalize the JBSI to a truly multivariate measure of synchrony.

Conclusions

I describe here a conceptually simple and computationally efficient method for determining the statistical significance of firing synchrony between two neurons and for quantifying synchrony as a normalized index, the JBSI. The method is based on the introduction of virtual spike jitter, but unlike previous implementations of this idea, it does not rely on computationally intensive generation of surrogate spike trains, and it uses jitter not only to test the statistical significance of spike coincidences but also to quantify synchrony. To evaluate the JBSI in comparison with previously used synchrony measures, I propose a set of five requirements from an ideal measure of synchrony. I show that the JBSI meets them all, unlike some commonly used synchrony indices such as the ECI and the CCC. First, the JBSI can be computed for any pair of spike trains, whether spontaneous or locked to a repeating stimulus and whether periodic or not. Second, the JBSI is well-normalized, in that it assumes values between 1 (highest possible synchrony for a given number of spikes in each train) and −1 (lowest possible synchrony), with 0 indicating chance-level synchrony. Third, the JBSI is independent of the firing rate, whereas the ECI is not. Fourth, the JBSI is independent of the firing rate differential between the two neurons, whereas the CCC is not. Finally, the JBSI is robust against co-modulations in firing rate, while both the ECI and the CCC show spurious synchrony when such co-modulations are present. I also show that a minor modification in the definitions of the ECI and the CCC results in an improved index, the ECIcor, which is robust under the third and fourth requirements and is therefore superior to both the ECI and the CCC. The ECIcor may indeed be the index of choice due to its computational simplicity, if firing rates are known to be stationary. However, the JBSI is the only index that meets all five requirements. By virtue of its robustness, the JBSI can be used to compare firing synchrony between experiments conducted under widely different experimental conditions and, as I demonstrate, it can also be used to estimate the temporal precision of firing in the system.

Methods

All computations and simulations were implemented in MathCad (PTC); the MathCad code is available from the author upon request. A MatLab routine for calculating the JBSI will be made available on the journal website.

Spike train simulations

For most simulations, the epoch length and/or the firing rate were adjusted so the generated trains consisted of about 1,000 spikes each. Each simulation started by generating two independent Poisson spike trains, as follows: for each 1 ms bin of the time epoch, a spike was considered fired if a randomly generated number between 0 and 1 was smaller than the predetermined firing rate (in spikes/ms); the precise time of the spike within the 1 ms bin was then randomly determined. To enforce a refractory period of 2 ms, if a spike was fired then no spikes were allowed in the next two bins. Next, firing coincidences, at a rate D and a precision C, were inserted as follows: for each spike of the reference train, if a randomly generated number was smaller than D, the spike was shifted forward to within ±C of the next spike of the target train. Finally, any spike in the shifted reference train that violated the refractory period was removed; typically, this resulted in the final spike count of the two spike trains differing by about 5% to 10% (if the initial, predetermined firing rates were equal). To generate co-modulations in the firing rate (Figure 2D), the predetermined firing rate was multiplied by a rectified sinusoidal function with a period of 1 s, raised to the power M: M was parametrically increased to vary the depth of the modulation.

Appendix A

A1. Algorithm for calculating the union set of all synchrony windows

Given the vector t2 representing the target spike train and the selected synchrony span τS, this algorithm returns a d by 2 matrix U, with each of the d rows representing a contiguous time segment [Um,0, Um,1], 0 ≤ m ≤ (d-1).

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M18">View MathML</a>

(A1.1)

A2. Algorithm for computing the probability vector pi (see equation 8)

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M19">View MathML</a>

The union of all synchrony windows, <a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M20">View MathML</a> , is first expressed as the d by 2 matrix U (Appendix A1). Next we calculate Ii,m, the intersection of the jitter window WJi with the segment [Um,0, Um,1], as follows (t1 is the vector representing the reference spike train):

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M21">View MathML</a>

(A2.1)

Note that if τS is smaller than ½ the smallest target interspike interval and therefore synchrony windows do not overlap, we have Uk,0= t2k- τS and Uk,1= t2k+ τS, and (A2.1) simplifies to

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M22">View MathML</a>

(A2.2)

Finally, we sum over all m segments and divide by the width of the jitter window to yield pi:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M23">View MathML</a>

(A2.3)

A3. Algorithm for calculating the probability of N successes in n1 trials with non-homogeneous success probabilities (adapted from [67])

We regard the reference spike train after a jitter as a Bernoulli series of n1 trials indexed on i, each with its own success (that is, synchrony) probability pi (calculated in Appendix A2) and failure probability qi = 1− pi. We calculate the probability PJ(N), 0 ≤ N ≤ n1, that exactly N of the spikes will be synchronous, by constructing a triangular n1 by (n1+ 1) matrix P recursively. The first three rows are:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M24">View MathML</a>

We define Pk,j = 0 for j > k.

Inspection shows that each row can be expressed in terms of the previous row:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M25">View MathML</a>

The final row gives us the desired probability, PJ(N) = Pn1,N.

Appendix B

A proof that the Jitter-Based Synchrony Index is bounded between −1 and 1

In the discussion below, the ratio of the jitter to synchrony windows will be indicated by α, that is, α = τJ/τS. The JBSI will be maximized when ‹NCJ›, that is, the expected coincidence count following a jitter, is minimized. This will happen when each of the spikes of Neuron 1 is at least τJ+ τS away from any spike of Neuron 2 other than spikes it is synchronized with, and is infinitesimally less than τS away from the spike it is synchronized with. Note that if one constructs a cross-correlogram of the two spike trains with a bin width of 2τS and the central bin symmetric about 0, the first condition is equivalent to saying that the two off-center bins should have 0 counts. Under these optimal conditions, the probability that a synchronous spike will remain synchronous after a random jitter is ½ for α 2 or 1/α for α ≥ 2, and PJ(N) becomes a binomial distribution with the binomial parameters p = 1/α and N = NC. A binomial distribution has an expected value of N·p, so for α ≥ 2 the JBSI becomes:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M26">View MathML</a>

(B1)

It is easily verified that the same holds for α ≤ 2.

In other words, in this optimal scenario, JBSI is equal to the observed coincidence rate RC. Since by definition RC ≤1, we have JBSI ≤1.

The JBSI can also attain negative values, indicating less-than-expected synchrony (not to be confused with ’anti-synchrony’, which is usually meant to indicate precise out-of-phase relationships). For example, if none of the spikes are synchronous but some of the reference spikes are within τJ+ τS of target spikes (that is, the center bin in the CCG is 0, but the off-center bins are not), there is a non-zero probability that some spikes will become synchronous after a jitter. This will render both the Z-score and the JBSI negative. To calculate the lowest possible values of the JBSI, we look at the extreme case in which each reference spike is infinitesimally more than τS away from some target spike, and therefore has a probability of 1/α (1/2 for α2264 2) of becoming synchronous after a jitter. Again we have a binomial distribution, yielding JBSI = −1 for α ≤ 2, or JBSI = −1/(α-1) for α ≥ 2. It is therefore both numerically convenient and advantageous to select α = 2, as this ratio will allow the maximal dynamic range for the JBSI.

Appendix C

Derivation of the Cross-Correlation Coefficient

Assume that we have recorded two simultaneous spike trains, during a recording epoch of duration T, with n1 and n2 spikes each, respectively. Assume n1 ≤ n2. We would like to calculate the probability P(N) of observing exactly N coincidences, 0 ≤ N ≤ n1. The CCC is based on the assumption that the two spike trains are independent and that the spikes occur randomly in time. We bin the epoch T into K bins, with bin width small enough so no more than one spike of each neuron can occur per bin. We can recreate the two spike trains in the following manner. First, we generate train #2, by distributing n2 spikes at random into n2 bins. We then generate train #1 by distributing the n1 spikes into n1 bins in the same random manner. What is the probability that N of these n1 bins already contain a spike from train #2? This can be solved by basic combinatorics, as follows. There are <a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M27">View MathML</a> distinct configurations of N objects (spikes, in our case) in n2 bins, where <a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M28">View MathML</a> is the binomial coefficient. For each of these configurations, there are <a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M29">View MathML</a> ways to place the remaining (n1-N) objects in the remaining (K-n2) bins, so in total there are <a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M30">View MathML</a> distinct configurations with exactly N spikes from train #1 falling into bins that already contain a spike from train #2. To convert this number into probability, we need to divide by the total number of possible configurations of n1 objects in K bins, which is <a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M31">View MathML</a> . The requested probability is therefore:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M32">View MathML</a>

(C1)

This is the well-known hypergeometric distribution, which has a mean (expected coincidences count) and a variance given, respectively, by:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M33">View MathML</a>

(C2)

and

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M34">View MathML</a>

(C3)

The CCC is defined as the Z-score of the observed coincidence count under the assumption of a hypergeometric distribution, normalized by √(K-1):

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M35">View MathML</a>

(C4)

To verify that this expression is indeed normalized, we substitute for NC the highest possible coincidence count, which is the number of spikes in the shorter train, n1:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M36">View MathML</a>

(C5)

Clearly, for all n1  n2, we will have CCCmax1, so the CCC is normalized. However, CCC as defined will only reach 1 if n1 = n2, that is, if the two firing rates are equal. This is a disadvantage compared to the JBSI, which will be 1 for perfectly synchronized trains, even if the two spike trains have very different rates. This suggests a simple way to correct the CCC -divide it by the CCCmax:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M37">View MathML</a>

If we choose the bin width to be 2τS, then K = T/2τS, and therefore (using equation 2) n1·n2/K = ‹NC

So, using equation 14:

<a onClick="popup('http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.neuralsystemsandcircuits.com/content/2/1/5/mathml/M38">View MathML</a>

(C6)

Abbreviations

CCC: cross-correlation coefficient; CCCcor: corrected CCC; CCG: cross-correlogram; ECIcor: corrected Excess Count Index; ECI: Excess Count Index; JBSI: Jitter-Based Synchrony Index; JSSI: Jitter-Sensitive Synchrony Index.

Competing interests

The author declares no competing financial interests.

Acknowledgments

The author thanks Drs. Peter Latham, Paul Brown and George Spirou for helpful comments on the manuscript. The author is supported by grant NS050437 from the National Institutes of Health. The author of this manuscript is fully and solely responsible for its content.

References

  1. Lestienne R: Spike timing, synchronization and information processing on the sensory side of the central nervous system.

    Prog Neurobiol 2001, 65:545-591. PubMed Abstract | Publisher Full Text OpenURL

  2. Usrey WM, Reid RC: Synchronous activity in the visual system.

    Annu Rev Physiol 1999, 61:435-456. PubMed Abstract | Publisher Full Text OpenURL

  3. Wang HP, Spencer D, Fellous JM, Sejnowski TJ: Synchrony of thalamocortical inputs maximizes cortical reliability.

    Science 2010, 328:106-109. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  4. Roy SA, Alloway KD: Coincidence detection or temporal integration? What the neurons in somatosensory cortex are doing.

    J Neurosci 2001, 21:2462-2473. PubMed Abstract | Publisher Full Text OpenURL

  5. Alonso JM, Usrey WM, Reid RC: Precisely correlated firing in cells of the lateral geniculate nucleus.

    Nature 1996, 383:815-819. PubMed Abstract | Publisher Full Text OpenURL

  6. Reyes AD: Synchrony-dependent propagation of firing rate in iteratively constructed networks in vitro.

    Nat Neurosci 2003, 6:593-599. PubMed Abstract | Publisher Full Text OpenURL

  7. Kumar A, Rotter S, Aertsen A: Spiking activity propagation in neuronal networks: reconciling different perspectives on neural coding.

    Nat Rev Neurosci 2010, 11:615-627. PubMed Abstract | Publisher Full Text OpenURL

  8. Usrey WM, Reppas JB, Reid RC: Paired-spike interactions and synaptic efficacy of retinal inputs to the thalamus.

    Nature 1998, 395:384-387. PubMed Abstract | Publisher Full Text OpenURL

  9. Uhlhaas PJ, Pipa G, Lima B, Melloni L, Neuenschwander S, Nikolic D, Singer W: Neural synchrony in cortical networks: history, concept and current status.

    Front Integr Neurosci 2009, 3:17. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  10. Lamme VA, Spekreijse H: Neuronal synchrony does not represent texture segregation.

    Nature 1998, 396:362-366. PubMed Abstract | Publisher Full Text OpenURL

  11. Roy A, Steinmetz PN, Hsiao SS, Johnson KO, Niebur E: Synchrony: a neural correlate of somatosensory attention.

    J Neurophysiol 2007, 98:1645-1661. PubMed Abstract | Publisher Full Text OpenURL

  12. de Oliveira SC, Thiele A, Hoffmann KP: Synchronization of neuronal activity during stimulus expectation in a direction discrimination task.

    J Neurosci 1997, 17:9248-9260. PubMed Abstract | Publisher Full Text OpenURL

  13. Riehle A, Grun S, Diesmann M, Aertsen A: Spike synchronization and rate modulation differentially involved in motor cortical function.

    Science 1997, 278:1950-1953. PubMed Abstract | Publisher Full Text OpenURL

  14. McClurkin JW, Optican LM, Richmond BJ, Gawne TJ: Concurrent processing and complexity of temporally encoded neuronal messages in visual perception.

    Science 1991, 253:675-677. PubMed Abstract | Publisher Full Text OpenURL

  15. Mechler F, Victor JD, Purpura KP, Shapley R: Robust temporal coding of contrast by V1 neurons for transient but not for steady-state stimuli.

    J Neurosci 1998, 18:6583-6598. PubMed Abstract | Publisher Full Text OpenURL

  16. Eggermont JJ: Neural interaction in cat primary auditory cortex II. Effects of sound stimulation.

    J Neurophysiol 1994, 71:246-270. PubMed Abstract | Publisher Full Text OpenURL

  17. Ito J, Maldonado P, Singer W, Grun S: Saccade-related modulations of neuronal excitability support synchrony of visually elicited spikes.

    Cereb Cortex 2011, 21:2482-2497. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  18. Jadhav SP, Wolfe J, Feldman DE: Sparse temporal coding of elementary tactile features during active whisker sensation.

    Nat Neurosci 2009, 12:792-800. PubMed Abstract | Publisher Full Text OpenURL

  19. Temereanca S, Brown EN, Simons DJ: Rapid changes in thalamic firing synchrony during repetitive whisker stimulation.

    J Neurosci 2008, 28:11153-11164. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  20. Wang Q, Webber RM, Stanley GB: Thalamic synchrony and the adaptive gating of information flow to cortex.

    Nat Neurosci 2010, 13:1534-1541. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  21. Swadlow HA, Beloozerova IN, Sirota MG: Sharp, local synchrony among putative feed-forward inhibitory interneurons of rabbit somatosensory cortex.

    J Neurophysiol 1998, 79:567-582. PubMed Abstract | Publisher Full Text OpenURL

  22. Eggermont JJ: Pair-correlation in the time and frequency domain. In Analysis of Parallel Spike Trains. Edited by Grun S, Rotter S. New York: Springer; 2010:77-102.

    [Destexhe A, Brette R (Series Editor): Springer Series in Computational Neuroscience]

    OpenURL

  23. Sears TA, Stagg D: Short-term synchronization of intercostal motoneurone activity.

    J Physiol 1976, 263:357-381. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  24. Yeh CI, Stoelzel CR, Weng C, Alonso JM: Functional consequences of neuronal divergence within the retinogeniculate pathway.

    J Neurophysiol 2009, 101:2166-2185. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  25. Fanselow EE, Richardson KA, Connors BW: Selective, state-dependent activation of somatostatin-expressing inhibitory interneurons in mouse neocortex.

    J Neurophysiol 2008, 100:2640-2652. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  26. Galarreta M, Hestrin S: A network of fast-spiking cells in the neocortex connected by electrical synapses.

    Nature 1999, 402:72-75. PubMed Abstract | Publisher Full Text OpenURL

  27. Hu EH, Bloomfield SA: Gap junctional coupling underlies the short-latency spike synchrony of retinal alpha ganglion cells.

    J Neurosci 2003, 23:6768-6777. PubMed Abstract | Publisher Full Text OpenURL

  28. Hu H, Ma Y, Agmon A: Submillisecond firing synchrony between different subtypes of cortical interneurons connected chemically but not electrically.

    J Neurosci 2011, 31:3351-3361. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  29. Grun S, Rotter S (Eds): Analysis of Parallel Spike Trains. New York: Springer; 2010. OpenURL

  30. Amarasingham A, Harrison MT, Hatsopoulos NG, Geman S: Conditional modeling and the jitter method of spike resampling.

    J Neurophysiol 2012, 107:517-531. PubMed Abstract | Publisher Full Text OpenURL

  31. Perkel DH, Gerstein GL, Moore GP: Neuronal spike trains and stochastic point processes. II. Simultaneous spike trains.

    Biophys J 1967, 7:419-440. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  32. Date A, Geman S, Bienenstock E: On the temporal resolution of neural activity.

    Technical Report, Brown University;

    1998

    PubMed Abstract | Publisher Full Text OpenURL

  33. Geman S, Amarasingham A, Harrison M, Hatsopoulos N: The statistical analysis of temporal resolution in the nervous system.

    Technical Report, Brown University;

    2008

    PubMed Abstract | Publisher Full Text OpenURL

  34. Amarasingham A: Statistical methods for the assessment of temporal structure in the activity of the nervous system.

    Doctoral dissertation, Brown University 2004. OpenURL

  35. Hatsopoulos N, Geman S, Amarasingham A, Bienenstock E: At what time scale does the nervous system operate?

    Neurocomputing 2003, 52–54:25-29. PubMed Abstract OpenURL

  36. Abeles M, Gat I: Detecting precise firing sequences in experimental data.

    J Neurosci Methods 2001, 107:141-154. PubMed Abstract | Publisher Full Text OpenURL

  37. Lestienne R, Tuckwell HC: The significance of precisely replicating patterns in mammalian CNS spike trains.

    Neuroscience 1998, 82:315-336. PubMed Abstract | Publisher Full Text OpenURL

  38. Pauluis Q, Baker SN, Olivier E: Precise burst synchrony in the superior colliculus of the awake cat during moving stimulus presentation.

    J Neurosci 2001, 21:615-627. PubMed Abstract | Publisher Full Text OpenURL

  39. Pazienti A, Maldonado PE, Diesmann M, Grun S: Effectiveness of systematic spike dithering depends on the precision of cortical synchronization.

    Brain Res 2008, 1225:39-46. PubMed Abstract | Publisher Full Text OpenURL

  40. Shmiel T, Drori R, Shmiel O, Ben-Shaul Y, Nadasdy Z, Shemesh M, Teicher M, Abeles M: Temporally precise cortical firing patterns are associated with distinct action segments.

    J Neurophysiol 2006, 96:2645-2652. PubMed Abstract | Publisher Full Text OpenURL

  41. Stark E, Abeles M: Unbiased estimation of precise temporal correlations between spike trains.

    J Neurosci Methods 2009, 179:90-100. PubMed Abstract | Publisher Full Text OpenURL

  42. Alonso JM, Yeh CI, Stoelzel CR: Visual stimuli modulate precise synchronous firing within the thalamus.

    Thalamus Relat Syst 2008, 4:21-34. PubMed Abstract | PubMed Central Full Text OpenURL

  43. Dan Y, Alonso JM, Usrey WM, Reid RC: Coding of visual information by precisely correlated spikes in the lateral geniculate nucleus.

    Nat Neurosci 1998, 1:501-507. PubMed Abstract | Publisher Full Text OpenURL

  44. Datta AK, Farmer SF, Stephens JA: Central nervous pathways underlying synchronization of human motor unit firing studied during voluntary contractions.

    J Physiol 1991, 432:401-425. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  45. Eggermont JJ: Properties of correlated neural activity clusters in cat auditory cortex resemble those of neural assemblies.

    J Neurophysiol 2006, 96:746-764. PubMed Abstract | Publisher Full Text OpenURL

  46. Abeles M: Local Cortical Circuits. Berlin: Springer; 1982. OpenURL

  47. Eggermont JJ: Stimulus induced and spontaneous rhythmic firing of single units in cat primary auditory cortex.

    Hear Res 1992, 61:1-11. PubMed Abstract | Publisher Full Text OpenURL

  48. Palm G, Aertsen AM, Gerstein GL: On the significance of correlations among neuronal spike trains.

    Biol Cybern 1988, 59:1-11. PubMed Abstract | Publisher Full Text OpenURL

  49. Halliday DM, Rosenberg JR: Time and frequency domain analysis of spike train and time series data. In Modern Techniques in Neuroscience Research. Edited by Windhorst U, Johansson H. Berlin: Springer-Verlag; 1999:503-543. OpenURL

  50. Quian Quiroga R, Kraskov A, Kreuz T, Grassberger P: Performance of different synchronization measures in real data: a case study on electroencephalographic signals.

    Phys Rev E Stat Nonlin Soft Matter Phys 2002, 65:041903. PubMed Abstract | Publisher Full Text OpenURL

  51. Ashida G, Wagner H, Carr C: Processing of phase-locked spikes and periodic signals. In Analysis of Parallel Spike Trains. Edited by Grun S, Rotter S. New-York: Springer; 2010:59-74.

    [Destexhe A, Brette R (Series Editor): Springer Series in Computational Neuroscience]

    OpenURL

  52. Meister M, Lagnado L, Baylor DA: Concerted signaling by retinal ganglion cells.

    Science 1995, 270:1207-1210. PubMed Abstract | Publisher Full Text OpenURL

  53. Ellaway PH, Murthy KS: The origins and characteristics of cross-correlated activity between gamma-motoneurones in the cat.

    Q J Exp Physiol 1985, 70:219-232. PubMed Abstract | Publisher Full Text OpenURL

  54. Datta AK, Stephens JA: Synchronization of motor unit activity during voluntary contraction in man.

    J Physiol 1990, 422:397-419. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  55. Halliday DM, Rosenberg JR, Breeze P, Conway BA: Neural spike train synchronization indices: definitions, interpretations, and applications.

    IEEE Trans Biomed Eng 2006, 53:1056-1066. PubMed Abstract | Publisher Full Text OpenURL

  56. Nordstrom MA, Fuglevand AJ, Enoka RM: Estimating the strength of common input to human motoneurons from the cross-correlogram.

    J Physiol 1992, 453:547-574. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  57. Turker KS, Powers RK: The effects of common input characteristics and discharge rate on synchronization in rat hypoglossal motoneurones.

    J Physiol 2002, 541:245-260. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  58. Aronov D: Fast algorithm for the metric-space analysis of simultaneous responses of multiple single neurons.

    J Neurosci Methods 2003, 124:175-179. PubMed Abstract | Publisher Full Text OpenURL

  59. Houghton C, Sen K: A new multineuron spike train metric.

    Neural Comput 2008, 20:1495-1511. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  60. Victor JD: Spike train metrics.

    Curr Opin Neurobiol 2005, 15:585-592. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  61. Kreuz T, Haas JS, Morelli A, Abarbanel HD, Politi A: Measuring spike train synchrony.

    J Neurosci Methods 2007, 165:151-161. PubMed Abstract | Publisher Full Text OpenURL

  62. Grun S, Diesmann M, Aertsen A: Unitary events in multiple single-neuron spiking activity: I. Detection and significance.

    Neural Comput 2002, 14:43-80. PubMed Abstract | Publisher Full Text OpenURL

  63. Aertsen AM, Gerstein GL, Habib MK, Palm G: Dynamics of neuronal firing correlation: modulation of “effective connectivity”.

    J Neurophysiol 1989, 61:900-917. PubMed Abstract | Publisher Full Text OpenURL

  64. Brody CD: Slow covariations in neuronal resting potentials can lead to artefactually fast cross-correlations in their spike trains.

    J Neurophysiol 1998, 80:3345-3351. PubMed Abstract | Publisher Full Text OpenURL

  65. Brody CD: Correlations without synchrony.

    Neural Comput 1999, 11:1537-1551. PubMed Abstract | Publisher Full Text OpenURL

  66. Ventura V, Cai C, Kass RE: Trial-to-trial variability and its effect on time-varying dependency between two neurons.

    J Neurophysiol 2005, 94:2928-2939. PubMed Abstract | Publisher Full Text OpenURL

  67. Thomas M, Taub A: Calculating binomial probabilities when the trial probabilities are unequal.

    J Statist Comput Simulation 1982, 14:125-131. Publisher Full Text OpenURL

  68. de Charms RC, Zador A: Neural representation and the cortical code.

    Annu Rev Neurosci 2000, 23:613-647. PubMed Abstract | Publisher Full Text OpenURL

  69. London M, Roth A, Beeren L, Hausser M, Latham PE: Sensitivity to perturbations implies high noise and suggests rate coding in cortex.

    Nature 2010, 466:123-127. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  70. Chi Z, Margoliash D: Temporal precision and temporal drift in brain and behavior of zebra finch song.

    Neuron 2001, 32:899-910. PubMed Abstract | Publisher Full Text OpenURL

  71. Desbordes G, Jin J, Weng C, Lesica NA, Stanley GB, Alonso JM: Timing precision in population coding of natural scenes in the early visual system.

    PLoS Biol 2008, 6:e324. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL