Correlating spike trains

We use the method proposed by (Cutts and Eglen, J Neurosci 2014) to compute the correlation between a pair of spike trains. This measure is called the “Spike Time Tiling Coefficient” (sttc). To illustrate these ideas, we first generate three related synthetic spike trains. A is a Poisson train at 1 Hz; B simply shifts each spike in A uniformly by up to +/- 0.1 seconds. Finally C is simply B shifted rightwards by 0.8 seconds.

(These correlation routines assume that the spike trains are sorted, smallest time first.)

require(meaRtools)
poisson_train <- function(n = 1000, rate = 1, beg = 0) {
  ## Generate a Poisson spike train with N spikes and firing rate RATE.
  ## BEG is time of start of recording
  ## Check that the histogram looks exponentially distributed.
  ## hist( diff(poisson_train()))
  x <- runif(n)
  isi <- log(x) / rate
  spikes <- beg - cumsum(isi)
  spikes
}

a = poisson_train()
b = sort(a + runif(length(a), -0.1, 0.1))
c = b + 0.8
allspikes = list(a=a, b=b, c=c)

To compute the correlation between A and B we need three extra bits of information. First, dt is the coincidence window (typically between 10 ms and 1 s); we then provide the beginning and end time in seconds of the two recordings.

range = range(unlist(allspikes))
beg = range[1]
end = range[2]
sttc(a, b, dt=0.01, rec_time = c(beg, end))
## [1] 0.09317701

To compute all pairwise correlations we can put all the spikes into a list. A upper-diagonal matrix is then returned with all pairwise correlations; the leading diagonal should be 1.0 for autocorrelations.

sttc_allspikes1(allspikes, 0.01, beg, end)
##      [,1]       [,2]        [,3]
## [1,]    1 0.09317701 0.002979672
## [2,]   NA 1.00000000 0.004984390
## [3,]   NA         NA 1.000000000

Since the STTC method was published, we have extended it to produce a correlogram-like profile. For a given value of tau, we cross correlate the spike trains A and B by adding tau to the time of each spike in B. This allows us to detect correlations when there are temporal delays. For this calculation, we simply also need to state the maximum absolute tau value (tau_max), and the step size for tau (tau_step) in sttcp. In the following plot, we compare our spike trains when using a range of coincidence windows, dt:

par(mfrow=c(2,2), las=1)
dt = 1.0;  plot(sttcp(a, b, dt = dt, beg = beg, end = end), main=sprintf("a x b: dt = %.3f",dt))
dt = 0.5;  plot(sttcp(a, c, dt = dt, beg = beg, end = end), main=sprintf("a x c: dt = %.3f",dt))
dt = 0.1;  plot(sttcp(a, c, dt = dt, beg = beg, end = end), main=sprintf("a x c: dt = %.3f",dt))
dt = 0.05; plot(sttcp(a, c, dt = dt, beg = beg, end = end), main=sprintf("a x c: dt = %.3f",dt))

plot of chunk unnamed-chunk-5

It clearly shows that a and c are ofset by about 0.8 seconds. We also see that for small values of dt, the cross-correlation peak drops significantly below 1.0.