DentateSpikes/code/getDSCSD.m

%get CSD of DSs
%requires CSD Plotter https://github.com/espenhgn/CSDplotter
%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;

addpath(genpath('/Users/dino/DATA/matlab/CSDplotter'));

%directories
drDS = 'DS/';
drMAT = 'MAT_EEG/';
drOUT = 'DS_CSD/';

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

%load DS and LFP
load([drDS fnEEG]);
load([drMAT fnEEG]);

%LFP channels
Nchannels = size(eegData,1);
    
%CSD properties
h = 50e-6; %inter-electrode distance meters
ex_cond = 0.3; %external conductivity S/m
top_cond = 0.3; %S/m
diam = 0.5e-3; %mm
gauss_sigma = 0.05e-3; %mm
filter_range = 5*gauss_sigma; % numeric filter must be finite in extent
dt = 0.5; %sampling time ms (for 2k);
scale_plot = 1; %focus
max_plot = 0;
CSDres = 200; %CSD resolution (fixed)
methodSpline = 1; %CSD method

if Nchannels == 32
    el_pos = 0.1:0.05:1.65; %our tetrode distances (array)
    el_pos = el_pos*1e-3; %mm
elseif Nchannels == 16
    el_pos = 0.1:0.05:0.85;
    el_pos = el_pos*1e-3; %mm
end

%create filter for CSD filtering
eegFS = 2000;
[b, gd] = getFIRbandpass( 5,200,4,40,eegFS );

winLen = 10; %window for CSD extraction
     
%filter
eegF = filter(b,1,eegData,[],2);
eegF = cat(2,eegF(:,gd+1:end),zeros(size(eegF,1),gd));

eegData = eegF;

%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

%get features
s = feat(:,1); %sample
wDG = feat(:,11); %width

dB = feat(:,3) - feat(:,9); %distance to min before
dA = feat(:,3) - feat(:,10); %distance to min before

dB = log(dB); %log (gamma)
dA = log(dA); %log (gamma)

mx = mean(dB);
sx = std(dB);

dB = (dB-mx)/sx; %zscore
dA = (dA-mx)/sx;

%criteria for DS
k = dB > 0.75 & dA > 0.75 & wDG > 5 & wDG < 25;

%select samples
sk = s(k);

resultsCSD = zeros(CSDres,length(sk));

%get CSDs
for sI = 1:length(sk)

    s = sk(sI);

    eeg = eegData(:,s-winLen:s+winLen);

    if methodSpline == 0
        %compute CSD - standard method
        %adding vaknin electrodes not to loose first and last channel
        %duplicate first and last channel
        eeg = cat(1,eeg(1,:),eeg,eeg(end,:));
        eeg = eeg + eeg;

        Nch = size(eeg,1);
        CSD = -ex_cond*D1(Nch,h)*eeeegAvergX;

        %filter iCSD
        b0 = 0.54; %center
        b1 = 0.23; %neighbor
        [n1,n2]=size(CSD);            
        CSD_add(1,:) = zeros(1,n2);   %add top and buttom row with zeros
        CSD_add(n1+2,:)=zeros(1,n2);
        CSD_add(2:n1+1,:)=CSD;        %CSD_add has n1+2 rows
        CSD = S_general(n1+2,b0,b1)*CSD_add; % CSD has n1 row

    else
        % compute spline iCSD:
        Fcs = F_cubic_spline(el_pos,diam,ex_cond,top_cond);
        [zs,CSD_cs] = make_cubic_splines(el_pos,eeg,Fcs);
        if gauss_sigma~=0 %filter iCSD
          [~,CSD_cs]=gaussian_filtering(zs,CSD_cs,gauss_sigma,filter_range);
        end
        CSD = CSD_cs;
    end

    CSDp = CSD(:,winLen+1);

    resultsCSD(:,sI) = CSDp;
end

%plot results

imagesc(resultsCSD,[-30000 30000])
colormap(jet)

text(5000,100,'<- outer mol.layer of superior DG blade (DSL)')
text(5000,119,'<- middle mol.layer of superior DG blade (DSM)')
text(5000,148,'<- hilus')
text(5000,175,'<- middle mol.layer of inferior DG blade (DSM)')
text(5000,192,'<- outer mol.layer of inferior DG blade (DSL)')

save([drOUT fnEEG],'resultsCSD','sk');

xlabel('Putative DS events');
ylabel('Depth')