Example of neural data analysis

Edgar Bermudez, PhD edgar.bermudez@uleth.ca

March 2013

Neural data analysis gives a better idea about the response or information processing/computation in a group of neurons.

Contents

Experiment description

The experiment consisted in recording neuronal data in anesthetized rats receving tactile stimulation under two conditions:

Recording system

The recording system receives the electrical signal from the probes and the stimulation signal with a determined sampling rate.

Silicon probes

The silicon probes have recording sites (called tetrodes) with different configurations. Each probe has an array of shanks where the recording sites are. In this case, one probe was inserted in S1 and another in mPFC.

Important data structures in a electrophysiology recording

Stimulus triggered activity

In order to have an idea about how neurons in the different areas respond to stimulation, we can display the stimulus triggered response (LFP and PSTH).

Load the data

datadir = 'C:\Users\edgar.bermudez\Documents\';
load([datadir 'LK10' ]);

Neurons for S1 and PFC

pfc_nns = find( SpkEle<=8); % pfc neurons
s1_nns = find( SpkEle > 8); % s1 neurons
nr_pfc_nns = length( pfc_nns ); % number of pfc neurons

Extracting useful information In order to target our analysis, we divide and label the data to have a better control. For example, we need to know:

Differentiate the two conditions

stim_mark = find(stim==3 | stim== -1 );
stim_gaps = diff( stim_mark ); %

Assign a label for each stage in the experiment

exper_stages = zeros(size(stim));
stim_cond = find(stim_gaps> median(stim_gaps)*100); %
stg_xtr = sort([1; stim_mark(1); stim_mark(stim_cond); stim_mark(stim_cond+1); stim_mark(end); length(stim)]);
for i=1:length(stg_xtr)-1;
	exper_stages(stg_xtr(i):stg_xtr(i+1)) = i;
end
exper_stages = exper_stages - 2;

Plot experiment setup

figure;plot( (1:length(exper_stages))*3.2/1000/60, exper_stages,'--r'); hold on;
plot( (1:length(exper_stages))*3.2/1000/60, stim/5);
title('Experiment example');
snapnow;
% 1: SP ure, 2: stim ure, 3: SP ure, 5:SP Amph, 6:stim Amph, 7:SP Amph

Some useful info

Firing rate for S1 neurons during each stage

ev_vals = [1 2 3 5 6 7];
for i = 1:length(ev_vals)
	fr(i,:) = mean(SpkCnt(find(exper_stages == ev_vals(i)),s1_nns));
end

Selection of active S1 neurons

select neurons with FR above threshold s1 and PFC

fr_th = 0.002;
act_s1nns = find( fr(2,:) > fr_th & fr(3,:) > fr_th);

Slicing activity into trial locked time windows

% tactile stim trig activity before drug
f_tac1 = find(stim==3 & exper_stages == 2);
f_tac1( find(diff( f_tac1 ) < 10 | diff( f_tac1 ) > 700)+1) = [];
w1=100; w2=750;
spk1=zeros(length(f_tac1),w1+w2+1,size(SpkCnt,2));


for tr=1:length(f_tac1)
	spk1(tr,:,:)=SpkCnt(f_tac1(tr)-w1:f_tac1(tr)+w2,:);
end
% after drug: find the trial times (3 is for tact stim and 6 after
% drug)
f_tac2 = find(stim==3 & exper_stages == 6);
f_tac2( find(diff( f_tac2 ) < 10 | diff( f_tac2 ) > 700)+1) = [];
spk2=zeros(length(f_tac2),w1+w2+1,size(SpkCnt,2));


for tr=1:length(f_tac2)
	% take activity before and after trials (trial-w1:trial+w2)
	spk2(tr,:,:)=SpkCnt(f_tac2(tr)-w1:f_tac2(tr)+w2,:);
end

% range of the reponse plots in ms
range1 = (-100:750)*3.2;
range2 = (-100:750)*3.2;

Plotting the response

figure();
subplot(1,2,1); title(['ureth']);
imagesc(sq(mean(spk1))');
subplot(1,2,2); title(['amph']);
imagesc(sq(mean(spk2))');

figure(); imagesc(sq(mean(spk1(:,:,s1_nns(act_s1nns))))')

figure();
subplot(1,2,1); hold on; title(['S1 mean PSTH']);
plot(range1, 312.5*sq(mean(sum(spk1(:,:,s1_nns),3))));
plot(range2, 312.5*sq(mean(sum(spk2(:,:,s1_nns),3))),'r');
ylim([0 450]);
xlabel('Time (ms)'); ylabel(' Hz ');
subplot(1,2,2); hold on; title(['PFC mean PSTH']);
plot(range1, 312.5*sq(mean(sum(spk1(:,:,pfc_nns),3))));
plot(range2, 312.5*sq(mean(sum(spk2(:,:,pfc_nns),3))),'r');
ylim([0 450]);
xlabel('Time (ms)'); ylabel(' Hz ');

Template matching

Template construction

Template is defined as the (summed evoked act - mean) / std average stim triggered activity of the S1 neurons

t_len = 60;
t_window = 101:100+t_len;
X = squeeze(sum(spk1(:,t_window,s1_nns(act_s1nns)),1))';

zs = zscore(X')';

% In order to have a less specific template, we smooth the activity pattern
% for every neuron
g_sigma = 40;
gs = normpdf(-g_sigma*2:g_sigma*2,0, g_sigma/2 )';

for n=1:size(zs,1)
	temp = zs( n, :);
	tmplt(n ,:) = conv2(  temp', gs, 'same');
end
template = tmplt(:,:);

Matching

shift = 1;
norm_sc = [];  tgt = [];
num_windows = 50000;
analysis_per = find(exper_stages == 3);
for j=1:num_windows
	winanalysis = analysis_per(1)+((j-1)*shift) : analysis_per(1)+((j-1)*shift)+t_len-1;
	%end
	% Keep track of the spike activity within the window
	% analysis
	tgt = SpkCnt(winanalysis, s1_nns(act_s1nns))';
	tgt_z = zscore(tgt')';

	norm_aft(j) = dot(template(:), tgt_z(:));
end

analysis_per = find(exper_stages == 1);
for j=1:num_windows
	winanalysis = analysis_per(1)+((j-1)*shift) : analysis_per(1)+((j-1)*shift)+t_len-1;
	%end
	% Keep track of the spike activity within the window
	% analysis
	tgt = SpkCnt(winanalysis, s1_nns(act_s1nns))';
	tgt_z = zscore(tgt')';

	norm_bef(j) = dot(template(:), tgt_z(:));
end

Results

figure();
subplot(1,2,1); hist(norm_bef);
title('TM before stim');
xlim([-50 50]);
ylim([0 35000]);
subplot(1,2,2); hist(norm_aft);
title('TM after stim');
xlim([-50 50]);
ylim([0 35000]);

References for further reading

Published with MATLAB® 7.14