meaRtools: Basics of correlating spike trains

In meaRtools: Micro-Electro Array (MEA) Analysis

require(knitr)
options(width=60)
opts_chunk$set(cache=TRUE)

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))

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)

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))

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.

Computing STTCs per well

First read in some data (these 6 recordings are simulated, first three from P9, and the last three from P11 mouse retina - Demas et al. 2003).

demas_platelayout = list(n_well = 6,
                        wells = paste0("w", 1:6),
                        n_well_r = 2,
                        n_well_c = 3,
                        layout = c(3, 2),
                        n_elec_r = 8,
                        n_elec_c = 8,
                        xlim = c(-100, 7200),
                        ylim = c(0, 6000),
                        spacing = 200,
                        corr_breaks = 0
                        )
add_plateinfo("demas-6well", demas_platelayout)
times = system.file("extdata/textreader/demas.times", package="meaRtools")
pos = system.file("extdata/textreader/demas.pos", package="meaRtools")
s = read_spikelist_text(times, pos, array="demas-6well")

Compute the pairwise STTC. Results are then stored in a data frame:

sttc_results = compute_sttc_by_well(s)
head(sttc_results)

We can also make a lattice plot to show the pairwise correlations as a function of distance within each well.

require(lattice)
xyplot(STTC ~ Distance | Well, data = sttc_results,
       main = "STTC by well",
       pch=20, xlab = "Distance (um)")

Finally, to save the pairwise correlations, we can simply save this data frame:

sttc_file = tempfile(fileext=".csv")
write.csv(sttc_results, file=sttc_file, row.names=FALSE)

The top of the file looks like this:

cat(readLines(sttc_file, 10), sep='\n')

meaRtools documentation built on May 1, 2019, 7:32 p.m.