DentateSpikes/code/getDS.m

%detect putative dentate spikes, extract their features and voltage profiles
%2021 Dino Dvorak, New York University, neuraldino@gmail.com
%Dentate spikes and external control of hippocampal function,
%Cell Reports, Volume 36, Issue 5, 2021,109497,ISSN 2211-1247,
%https://doi.org/10.1016/j.celrep.2021.109497

clear; close all; clc;

eegFS = 2000; %sample rate

%input LFP file (channels x samples)
fnEEG = '2016_09_30_PP12_334_PRETRAIN_BASE_eeg.mat';

%DG channel (found using initial visual inspection)
chDG = 24;

%directories
drMAT = 'MAT_EEG/'; %LFPs
drOUT = 'DS/'; %OUTPUT

%create filter for detection
%pass-band 5-100Hz, 4Hz transition band, 40Hz attenuation
[b, gd] = getFIRbandpass( 5,100,4,40,eegFS );

chAround = 5; %+/- channels around for DS profile

winBA = [60 20]; %search window for maxima before/after in samples (-60..-20samples or +20..+60 samples)

thMin = []; %smallest threshold to accept (saves data)

%features
labelsFeat = {'sample','promZ','peakZ','peakRaw','promZbefore','promZafter',... %1..6
        'peakZbefore','peakZafter','minZbefore','minZafter',...%7..10
        'wDG','wBA','w50','tB','tA','tminB','tminA','acc','speed','max10_30before','max10_30after'}; %11..21

%load LFP
load([drMAT fnEEG]);

eegData = double(eegData);

Nchannels = size(eegData,1);

%DG electrode
eegDG = eegData(chDG,:);

%filter 2-100Hz
chST = chDG - chAround;
chED = chDG + chAround;
if chED > Nchannels; chED = Nchannels; end;
if chST < 1; chST = 1; end;
channelsProfiles = chST:chED;
chDGk = find(channelsProfiles == chDG);
eegF = filter(b,1,eegData(chST:chED,:),[],2);
eegF = cat(2,eegF(:,gd+1:end),zeros(size(eegF,1),gd));
%DG channel
eegFDG = eegF(chDGk,:);

%DG - maxima and minima (extract width,prominence)
[MaxesDG,kMaxesDG,wDG,pDG] = findpeaks(eegFDG);
[~,kMinsDG] = findpeaks(1.01*max(eegFDG) - eegFDG);

%for each maxima in DG
feat = nan(length(kMaxesDG),23);

for i = 10:length(kMaxesDG)-10 %skip first couple of minima for before/after min-max computations
    s = kMaxesDG(i); %sample
    pDGi = pDG(i); %prominence of peak in z-score signal
    wDGi = wDG(i)/eegFS*1000; %width of peak in z-score signal
    pDGB = pDG(i-1); %prominence of peak before in z-score signal
    pDGA = pDG(i+1); %prominence of peak after in z-score signal
    tDGB = (s - kMaxesDG(i-1))/eegFS*1000; %mseconds to maxima before
    tDGA = (kMaxesDG(i+1) - s)/eegFS*1000; %mseconds to maxima after
    mDG = MaxesDG(i); %max value in z-scored filtered signal
    acci = nan; %acceleration
    mDGB = MaxesDG(i-1); %preceding maxima in z-scored filtered signal
    mDGA = MaxesDG(i+1); %following maxima in z-scored filtered signal

    %indexes preceeding and following minima
    kB = kMinsDG(kMinsDG < s);
    kA = kMinsDG(kMinsDG > s);
    if isempty(kB) || isempty(kA); continue; end;
    kB = kB(end);
    kA = kA(1);
    %time to minima before
    tmB = (s - kB)/eegFS*1000; %mseconds to minima before
    tmA = (kA - s)/eegFS*1000; %mseconds to minima after

    %maxima in raw signal
    mr = max(eegDG(s-10:s+10));

    %peak to to minima before and after
    mB = eegFDG(kB);
    mA = eegFDG(kA);

    %speed
    spi = nan;

    %width between minima
    wBA = (kA-kB)/eegFS*1000;

    %width 50%
    w50 = nan;
    k50B = nan;
    k50A = nan;
    if mB > mA %before higher
        level = (mDG + mB)/2;
    else %after higher
        level = (mDG + mA)/2;
    end
    x = eegFDG(kB:kA);
    %cut before maxima
    k = find(x(1:end-1) <= level & x(2:end) > level);
    if ~isempty(k)
        k50B = k(1);
    end
    %cut after maxima
    k = find(x(1:end-1) >= level & x(2:end) < level);
    if ~isempty(k)
        k50A = k(end);
    end
    if ~isnan(k50B) && ~isnan(k50A)
        w50 = (k50A-k50B)/eegFS*1000;
    end

    %winBA = [60 20];
    %maxima before
    x = eegFDG(s-winBA(1):s-winBA(2));
    maxB = max(x);
    %maxima after
    x = eegFDG(s+winBA(2):s+winBA(1));
    maxA = max(x);

    %minima in raw before
    mrB = min(eegDG(kB-10:kB+10));

    %minima in raw after
    mrA = min(eegDG(kA-10:kA+10));

    feat(i,1) = s; %sample max sDG
    feat(i,2) = pDGi; %prominence of peak in z-score signal
    feat(i,3) = mDG;  %max value in z-scored filtered signal
    feat(i,4) = mr; %maxima in raw
    feat(i,5) = pDGB; %prominence of peak before in z-score signal
    feat(i,6) = pDGA; %prominence of peak after in z-score signal
    feat(i,7) = mDGB; %preceding maxima in z-scored filtered signal
    feat(i,8) = mDGA; %following maxima in z-scored filtered signal

    feat(i,9) = mB; %minima preceding peak in z-score signal
    feat(i,10) = mA; %minima following peak in z-score signal

    feat(i,11) = wDGi; %width of peak in z-score signal
    feat(i,12) = wBA; %width between minima
    feat(i,13) = w50; %width 50%

    feat(i,14) = tDGB; %mseconds to maxima before
    feat(i,15) = tDGA; %mseconds to maxima after
    feat(i,16) = tmB; %mseconds to minima before
    feat(i,17) = tmA; %mseconds to minima after

    feat(i,18) = acci; %acceleration
    feat(i,19) = spi; %speed

    feat(i,20) = maxB; %maxima 30-10ms before
    feat(i,21) = maxA; %maxima 10-30ms after

    feat(i,22) = mrB; %minima in raw before
    feat(i,23) = mrA; %minima in raw after
end

%remove nans
knan = ~isnan(feat(:,1));
feat = feat(knan,:);

%apply thMin
if ~isempty(thMin)
    k = feat(:,3) > thMin;
    %d = feat(:,3) - feat(:,9); %dist peak to minima before
    %k = d > 2;
    feat = feat(k,:);
end

%extract profiles around maxima
s = feat(:,1); %samples;
profilesDS = eegF(:,s);

%save putative DS events
save([drOUT fnEEG],'feat','chDG','profilesDS','channelsProfiles','labelsFeat');
     
%plot output
figure;
subplot(1,2,1);
plot(mean(profilesDS,2),1:size(profilesDS,1))
title('Volrage DS profile');

subplot(1,2,2);
histogram(feat(:,11));
title('Event width');