% This code computes the Burstlet Rate & Fano factor 

% These analyses were performed for O'Neill et al., Time-dependent
% homeostatic mechanisms underlie BDNF action on neural circuitry. Comms
% Bio, 2023.

% This function was written by Erin D. Anderson and can be
% accessed at https://www.seas.upenn.edu/~molneuro/

% Last Updated: 11/14/2023

function spikeAnalysis = bdnfHomeostasisMEA_spikeAnalysis(folderName, overwrite)

% Collect files to run
fileNames = dir(folderName);
for jj = length(fileNames):-1:1 % get rid of all the irrelevant files
    if fileNames(jj).isdir
        fileNames(jj) = [];
    elseif ~strcmp(fileNames(jj).name(end-3:end),'.mat')
        fileNames(jj) = [];
    end
end

% Initialize table
columnNames = {'FileName','CellDensity','RecordingWidth','nNeurons','nRegions',...
    'BDNF1','BDNF2','InjuryFraction1','InjuryFraction2'...
    'BurstRate_Pre','BurstRate_During','BurstRate_Post','BurstRate_Norm','BurstRate_Delta',...
    'FF_Pre','FF_During1','FF_During2','FF_Post','FF_Norm','FF_Delta'};
columnTypes = [{'string'}; repmat({'double'},[8,1]); repmat({'cell'},[length(columnNames)-1-8,1])];

spikeTableMEA = table('Size',[length(fileNames) length(columnNames)],...
    'VariableNames',columnNames,'VariableTypes',columnTypes);

% Run analysis
for jj = 1:length(fileNames)
    mat = matfile(fullfile(folderName,fileNames(jj).name));
    if ismember('nRegions', who('-file',fullfile(folderName,fileNames(jj).name)))
        if ~ismember('spikeDataMEA',who('-file',fullfile(folderName,fileNames(jj).name))) || overwrite
            spikeDataMEA = struct;
            
            % Split into pre, during, and post segments
            [segmentedSpikeTimes,segmentedSpikeIndexes] = ...
                splitIntoSegments(mat.spikeTimesGrid,mat.spikeIndexesGrid,...
                mat.PreInjurySimTimeInSeconds,...
                [mat.SimTimeInSeconds,mat.InjurySimTimeInSeconds,mat.InjurySimTimeInSeconds,mat.PostInjurySimTimeInSeconds]);
            
            % Burst rate
            [~,~,~,spikeDataMEA.burstRate.pre,~,~,...
                spikeDataMEA.nSpikesPerBurst_std.pre] = INSIbyElectrode(segmentedSpikeIndexes{1},segmentedSpikeTimes{1},...
                mat.SimTimeInSeconds,0,mat.nRegions);
            [~,~,~,spikeDataMEA.burstRate.during1,~,~,...
                spikeDataMEA.nSpikesPerBurst_std.during1] = INSIbyElectrode(segmentedSpikeIndexes{2},segmentedSpikeTimes{2},...
                mat.InjurySimTimeInSeconds,0,mat.nRegions);
            [~,~,~,spikeDataMEA.burstRate.during2,~,~,....
                spikeDataMEA.nSpikesPerBurst_std.during2] = INSIbyElectrode(segmentedSpikeIndexes{3},segmentedSpikeTimes{3},...
                mat.InjurySimTimeInSeconds,0,mat.nRegions);
            [~,~,~,spikeDataMEA.burstRate.post,~,~,....
                spikeDataMEA.nSpikesPerBurst_std.post] = INSIbyElectrode(segmentedSpikeIndexes{4},segmentedSpikeTimes{4},...
                mat.PostInjurySimTimeInSeconds,0,mat.nRegions);
            spikeDataMEA.burstRate.norm = spikeDataMEA.burstRate.post./spikeDataMEA.burstRate.pre;
            spikeDataMEA.burstRate.delta = spikeDataMEA.burstRate.post - spikeDataMEA.burstRate.pre;
            
            
            % Fano factor
            spikeDataMEA.FF.pre = fanoFactorFn(segmentedSpikeIndexes{1}-1,segmentedSpikeTimes{1},mat.SimTimeInSeconds,...
                mat.PreInjurySimTimeInSeconds,mat.nRegions); % have to manually input pre-sim time for FF bc need to know exact start time
            spikeDataMEA.FF.during1 = fanoFactorFn(segmentedSpikeIndexes{2}-1,segmentedSpikeTimes{2},mat.InjurySimTimeInSeconds,...
                mat.PreInjurySimTimeInSeconds + mat.SimTimeInSeconds,mat.nRegions);
            spikeDataMEA.FF.during2 = fanoFactorFn(segmentedSpikeIndexes{3}-1,segmentedSpikeTimes{3},mat.InjurySimTimeInSeconds,...
                mat.PreInjurySimTimeInSeconds + mat.SimTimeInSeconds + mat.InjurySimTimeInSeconds,mat.nRegions);
            spikeDataMEA.FF.post = fanoFactorFn(segmentedSpikeIndexes{4}-1,segmentedSpikeTimes{4},mat.PostInjurySimTimeInSeconds,...
                mat.PreInjurySimTimeInSeconds + mat.SimTimeInSeconds + mat.InjurySimTimeInSeconds*2,mat.nRegions);
            spikeDataMEA.FF.norm = spikeDataMEA.FF.post./spikeDataMEA.FF.pre;
            spikeDataMEA.FF.delta = spikeDataMEA.FF.post - spikeDataMEA.FF.pre;
            
%             save(fullfile(folderName,fileNames(jj).name),'spikeDataMEA','-append');
        else
            spikeDataMEA = mat.spikeDataMEA;
        end
        
        spikeTableMEA.FileName(jj) = fileNames(jj).name;
        spikeTableMEA.CellDensity(jj) = mat.CellDensity;
        spikeTableMEA.RecordingWidth(jj) = mat.RecordingWidth;
        spikeTableMEA.nNeurons(jj) = double(mat.N);
        spikeTableMEA.nRegions(jj) = mat.nRegions;
        spikeTableMEA.BDNF1(jj) = mat.BDNF(1,1);
        spikeTableMEA.BDNF2(jj) = mat.BDNF(1,2);
        spikeTableMEA.InjuryFraction1(jj) = mat.glutamate(1,1);
        spikeTableMEA.InjuryFraction2(jj) = mat.glutamate(1,2);
        
        spikeTableMEA.BurstRate_Pre{jj} = spikeDataMEA.burstRate.pre;
        spikeTableMEA.BurstRate_During1{jj} = spikeDataMEA.burstRate.during1;
        spikeTableMEA.BurstRate_During2{jj} = spikeDataMEA.burstRate.during2;
        spikeTableMEA.BurstRate_Post{jj} = spikeDataMEA.burstRate.post;
        spikeTableMEA.BurstRate_Norm{jj} = spikeDataMEA.burstRate.norm;
        spikeTableMEA.BurstRate_Delta{jj} = spikeDataMEA.burstRate.delta;
        
        spikeTableMEA.FF_Pre{jj} = spikeDataMEA.FF.pre;
        spikeTableMEA.FF_During1{jj} = spikeDataMEA.FF.during1;
        spikeTableMEA.FF_During2{jj} = spikeDataMEA.FF.during2;
        spikeTableMEA.FF_Post{jj} = spikeDataMEA.FF.post;
        spikeTableMEA.FF_Norm{jj} = spikeDataMEA.FF.norm;
        spikeTableMEA.FF_Delta{jj} = spikeDataMEA.FF.delta;
        
        clearvars spikeDataMEA
    else
        warndlg('Make sure to convert to MEA prior to running analysis!')
    end
end

% % save spikeTable for later
datestamp = clock;
datestamp = [num2str(datestamp(1)), sprintf('%02d',datestamp(2)),sprintf('%02d',datestamp(3))];
save(fullfile(folderName,['spikeTableMEA_',datestamp,'.mat']),'spikeTableMEA');

%% Split into different conditions
[combos,~,jj] = unique([spikeTableMEA.BDNF1,spikeTableMEA.BDNF2,...
    spikeTableMEA.InjuryFraction1,spikeTableMEA.InjuryFraction2],'rows');
RowsForEachCombo = accumarray(jj, find(jj), [], @(rows){rows});
TempLengths = [];

for i = 1:length(RowsForEachCombo)
    TempLengths = [TempLengths length(RowsForEachCombo{i})];
end

nCombos = size(combos,1);
nRegions = spikeTableMEA.nRegions(1); % they all have the same number of regions

burstRatePre = NaN(max(TempLengths)*nRegions,nCombos);
burstRatePost = NaN(max(TempLengths)*nRegions,nCombos);
fanoFactorPre = NaN(max(TempLengths)*nRegions,nCombos);
fanoFactorPost = NaN(max(TempLengths)*nRegions,nCombos);

for ii = 1:length(RowsForEachCombo)
    RowsForEachCombo{ii} = sort(RowsForEachCombo{ii}); % order matters because matching
    for jj = 1:numel(RowsForEachCombo{ii})
        burstRatePre(nRegions*(jj-1)+1:nRegions*jj,ii) = spikeTableMEA.BurstRate_Pre{RowsForEachCombo{ii}(jj)};
        burstRatePost(nRegions*(jj-1)+1:nRegions*jj,ii) = spikeTableMEA.BurstRate_Post{RowsForEachCombo{ii}(jj)};
        fanoFactorPre(nRegions*(jj-1)+1:nRegions*jj,ii) = spikeTableMEA.FF_Pre{RowsForEachCombo{ii}(jj)};
        fanoFactorPost(nRegions*(jj-1)+1:nRegions*jj,ii) = spikeTableMEA.FF_Post{RowsForEachCombo{ii}(jj)};
    end
end

%% Make Summary Table

comboNames = cell(nCombos,1);
for ii = 1:nCombos
    if combos(ii,1) == 0 && combos(ii,2) == 0 && ... % no BDNF
        combos(ii,3) == 0 && combos(ii,4) == 0 % or injury
            comboNames{ii} = 'Control';
    elseif combos(ii,3) ~= 0 % if injured
        if combos(ii,2) ~= 0 % and BDNF'd
            comboNames{ii} = 'Injury then BDNF';
        else % no BDNF
            comboNames{ii} = 'Injury';
        end
    end
end

columnNames = {'Name','BurstRatePre','BurstRatePost','FanoFactorPre','FanoFactorPost'};
columnTypes = {'string','double','double','double','double'};

spikeAnalysis = table('Size',[nCombos length(columnNames)],...
    'VariableNames',columnNames,'VariableTypes',columnTypes);

spikeAnalysis.Name = comboNames;
spikeAnalysis.BurstRatePre = burstRatePre';
spikeAnalysis.BurstRatePost = burstRatePost';
spikeAnalysis.FanoFactorPre = fanoFactorPre';
spikeAnalysis.FanoFactorPost = fanoFactorPost';