############################################################################# # R codes for Figure 4 in Takahashi et al. (2019) eNeuro # ############################################################################# ############################################################################# # Example spike count data - Replace neuron1 and neuron2 with your own data! ############################################################################# # example1 for Fig4(a), nonstationary pair neuron1 = c(400,400,240,320,440,600,400,960,760,920,880,1240,1000,760,680,680,640,760,880,680,760,1040,920,920,1040,960,680,720,920,880,1000,760,1040,960,1000,1360,1160,1040,960,1280)/40/1.9 neuron2 = c(1040,1200,1520,1920,800,1160,1160,1360,1960,1600,1520,2080,1720,2000,2000,2040,2400,1880,1600,2040,1720,2120,1920,2000,1920,1960,1760,1960,1880,2000,2000,2040,2200,2120,2360,2480,2400,1840,2040,1720)/40/1.9 # example2 for Fig4(b), stationary pair neuron1 = c(0,1640,120,440,0,200,200,0,0,0,1240,160,40,1800,480,560,40,1520,800,200,1800,0,0,1240,80,0,1360,1360,200,400,1160,440,680,0,80,0,200,1760,680,40)/40/1.5 neuron2 = c(0,1560,320,480,160,360,360,240,400,40,1160,760,80,1760,400,760,280,1360,840,360,1640,160,80,720,320,40,2040,1560,320,280,480,680,680,280,200,0,160,800,320,400)/40/1.5 # plot firing rates for left figures par(cex=1.4, mar=c(4,4,0.3,0.5)) plot(neuron1, type="o", ylim=c(0, max(neuron1, neuron2)), xlab="trial", ylab="firing rate (Hz)", lwd=1, col="purple", pch=16, cex=0.7) points(neuron2, type="o", lwd=1 , col="chartreuse3", pch=17, cex=0.7) ############################################################################# # computation of correlograms for spike count noise correlations ############################################################################# cor_conventional={}; p_conventional={}; cor_proposed={}; # initialization Shift=10; # maximum shift in correlograms for(shift in -Shift:Shift){ # trial shift for correlograms if (shift>=0){ x = neuron1[(1+shift):length(neuron1)]; y = neuron2[1:(length(neuron2)-shift)]; }; if (shift<0){ y = neuron2[(1-shift):length(neuron2)]; x = neuron1[1:(length(neuron1)+shift)]; }; # conventional correlation coefficient cor_conventional[shift+Shift+1] = cor(x,y); # conventional correlation coefficient p_conventional[shift+Shift+1] = cor.test(x,y)$p.value; # p-value for conventional cor # proposed method x1 = x[seq(1,length(x)-1,by=2)]; x2 = x[seq(2,length(x),by=2)]; y1 = y[seq(1,length(y)-1,by=2)]; y2 = y[seq(2,length(y),by=2)]; x3 = x[seq(3,length(x),by=2)]; x4 = x[seq(2,length(x)-1,by=2)]; y3 = y[seq(3,length(y),by=2)]; y4 = y[seq(2,length(y)-1,by=2)]; A1 = mean((x1^2+x2^2)-(x1+x2)^2 /2); C1 = mean((y1^2+y2^2)-(y1+y2)^2 /2); B1 = mean((x1*y1+x2*y2)-(x1+x2)*(y1+y2)/2); A2 = mean((x3^2+x4^2)-(x3+x4)^2 /2); C2 = mean((y3^2+y4^2)-(y3+y4)^2 /2); B2 = mean((x3*y3+x4*y4)-(x3+x4)*(y3+y4)/2); A = (A1+A2)/2; C = (C1+C2)/2; B = (B1+B2)/2; # take average of two styles cor_proposed[1+Shift+shift] = B/sqrt(A*C); } # plot correlograms for right figures X = seq(-Shift,Shift,by=2); plot(X, cor_conventional[seq(1,2*Shift+1,by=2)], ty='l', col='black', ylim=c(-1,1), xlab='trial shift', ylab='cross correlation', axes=T, lwd=3); lines(X, cor_proposed[seq(1,2*Shift+1,by=2)], lwd=3, col="red"); lines(c(-Shift-5,Shift+5),c(0,0),lty=3); lines(c(0,0), c(-1.2, 1.2), lty=3); legend("topright",c('conventional cor ', 'proposed cor '), col=c('black', 'red'), lwd=2, cex=0.8) ############################################################################# # simulation for accompanying statistical test - can take a few minutes ############################################################################# N = length(neuron1); # match simulation data size to your own data size cc_sim = rep(0, 100000); # initialization for(i in 1:100000){ # generate random data from standard normal distributions x = rnorm(N); y = rnorm(N) # compute the proposed correlation coefficients x1 = x[seq(1,length(x)-1,by=2)]; x2 = x[seq(2,length(x),by=2)]; y1 = y[seq(1,length(y)-1,by=2)]; y2 = y[seq(2,length(y),by=2)]; x3 = x[seq(3,length(x),by=2)]; x4 = x[seq(2,length(x)-1,by=2)]; y3 = y[seq(3,length(y),by=2)]; y4 = y[seq(2,length(y)-1,by=2)]; A1 = mean((x1^2+x2^2)-(x1+x2)^2 /2); C1 = mean((y1^2+y2^2)-(y1+y2)^2 /2); B1 = mean((x1*y1+x2*y2)-(x1+x2)*(y1+y2)/2); A2 = mean((x3^2+x4^2)-(x3+x4)^2 /2); C2 = mean((y3^2+y4^2)-(y3+y4)^2 /2); B2 = mean((x3*y3+x4*y4)-(x3+x4)*(y3+y4)/2); A = (A1+A2)/2; C = (C1+C2)/2; B = (B1+B2)/2; cc_sim[i] = B/sqrt(A*C); } # compute the p-value as a percentile in the simulation results (cc_sim) p_one_sided = mean(cor_proposed[Shift+1]