In experiments measuring neural network activity, such as in studies using GCaMP, it’s often of interest to find neurons that have related activity patterns.

Finding neurons that are coactivated can easily be done using principal components analysis. However, you may also be interested in finding neurons that have a single upstream source driving their activity, but their activation timings are temporally offset.

Since these activation patterns don’t align, and are in a sea of other activity, there may not seem to be an easy way to identify these coupled neurons. Thankfully there is a method that can vastly reduce the complexity of this task. It’s a correlational method called cross-correlation, and it can be used to find, not just related signals, but also correct for offset in spike/signal timing.

Below, the script demonstrates how to use the cross-correlation function (here, in MATLAB using the xcorr() function) by loading three related but asynchronous signals into the workspace, computing their offset, and aligning them.

load relatedsig

figure(1)

ax(1) = subplot(3,1,1);
plot(s1)
ylabel('s_1')
axis tight

ax(2) = subplot(3,1,2);
plot(s2)
ylabel('s_2')
axis tight

ax(3) = subplot(3,1,3);
plot(s3)
ylabel('s_3')
axis tight
xlabel('Samples')

linkaxes(ax,'x')

Next, we use the xcorr() function to compute the cross-correlations between the three pairs of signals; then normalize them so their maximum value is one.

[C21,lag21] = xcorr(s2,s1);
C21 = C21/max(C21);

[C31,lag31] = xcorr(s3,s1);
C31 = C31/max(C31);

[C32,lag32] = xcorr(s3,s2);
C32 = C32/max(C32);

The locations of the maximum values of the cross-correlations indicate time leads or lags.

[M21,I21] = max(C21);
t21 = lag21(I21);

[M31,I31] = max(C31);
t31 = lag31(I31);

[M32,I32] = max(C32);
t32 = lag31(I32);

Plot the cross-correlations. In each plot display the location of the maximum.

figure

subplot(3,1,1)
plot(lag21,C21,[t21 t21],[-0.5 1],'r:')
text(t21+100,0.5,['Lag: ' int2str(t21)])
ylabel('C_{21}')
axis tight
title('Cross-Correlations')

subplot(3,1,2)
plot(lag31,C31,[t31 t31],[-0.5 1],'r:')
text(t31+100,0.5,['Lag: ' int2str(t31)])
ylabel('C_{31}')
axis tight

subplot(3,1,3)
plot(lag32,C32,[t32 t32],[-0.5 1],'r:')
text(t32+100,0.5,['Lag: ' int2str(t32)])
ylabel('C_{32}')
axis tight
xlabel('Samples')

We can see that s2 leads s1 by 350 samples; s3 lags s1 by 150 samples. Thus s2 leads s3 by 500 samples. Line up the signals by clipping the vectors with longer delays.

s1 = s1(-t21:end);
s3 = s3(t32:end);

ax(1) = subplot(3,1,1);
plot(s1)
ylabel('s_1')
axis tight

ax(2) = subplot(3,1,2);
plot(s2)
ylabel('s_2')
axis tight

ax(3) = subplot(3,1,3);
plot(s3)
ylabel('s_3')
axis tight
xlabel('Samples')

linkaxes(ax,'x')

The signals are now synchronized and ready for further processing.

---

Interested seeing more stuff like this? Head over to my GitHub page