Burstanalysis.m

% 
function [Results, accommresult] = Burstanalysis(abffile,graph_on,start,stop,xminmax,channelidx,currentidx)

    if nargin ~= 7
        error('Not enough input arguments')
    end
    [~,filename,~] = fileparts(abffile);    
    % abffile = '12d06012.abf'
    [d,si,h]=abfload(abffile); 
    numsweeps = size(d,3); % number of sweeps
	accommodating = [];
%%    
    for sweep = 1:numsweeps
            Results(sweep).filename = abffile;
            data = [];
            data = d(:,channelidx,sweep);
            timemes = 1:length(data)*(si/1000);
            Results(sweep).data = data;
            Results(sweep).SI = si;
            Results(sweep).currentstimulus = d(:,currentidx,sweep);
            %% Get baseline data
            [Results(sweep).baseline_potential, Results(sweep).baseline_potentialstd] = baseline(data,start,stop);
            Results(sweep).baseline_timerange = [start*(si/1000); stop*(si/1000)]; % need to convert to time. 
            meantracebl(sweep,1) = Results(sweep).baseline_potential;
            Results(sweep).blcurrent = mean(d(start:stop,currentidx,sweep));
            %% Find peak relative to holding potential
            thresh = -10;
            [pks, loc] = findpeaks(data,'minpeakheight',thresh,'minpeakdistance',1/(si/1000));
            %% Check to make sure peaks aren't occuring within 1ms of each other
            IPI = diff(loc);%
            violates = [];
            violates = find(IPI < (1000/si));
            if ~isempty(violates)
                violates(:,2) = violates + 1;
                for i = 1:size(violates,1)
                    if pks(violates(i,1)) > pks(violates(i,2))
                        violates(i,1) = 0;
                    elseif pks(violates(i,1)) < pks(violates(i,2))
                        violates(i,2) = 0;
                    elseif pks(violates(i,1)) == pks(violates(i,2))
                        violates(i,2) = 0;
                    end
                end
                violates = violates(:,1) + violates(:,2);
                pks(violates) = [];
                loc(violates) = [];
            end
            %% Check if spikes occuring during baseline
            withinbl = ismember(loc,start:stop+1/(si/1000));
            if ismember(1,withinbl) % if there are spikes within baseline interval
                % calculate new baseline without spikes within it
                [Results(sweep).baseline_potential, Results(sweep).baseline_potentialstd, pointstouse]...
                    = nospikebaseline(data,si,start,stop+1/(si/1000),loc,withinbl);
                meantracebl(sweep,1) = Results(sweep).baseline_potential;
                Results(sweep).blcurrent = mean(d(pointstouse,currentidx,sweep));
                % removes peaks that are occuring within bl interval
                pks(withinbl) = [];
                loc(withinbl) = []; 
            end
            %% Calculate 1 second current (3.5sec_ds)
            start1sec = stop+1;
            stop1sec = stop+(3500/(si/1000));
%             stop150ms = stop+(50/(si/1000));
            Results(sweep).current1sec = mean(d(start1sec:stop1sec,currentidx,sweep));
            Results(sweep).current = Results(sweep).current1sec - Results(sweep).blcurrent;
            %%
            if ~isempty(pks)
                peak_times = loc*(si/1000);
                peak_amps = pks-Results(sweep).baseline_potential; 
                numspikes = length(peak_times);
                Results(sweep).num_spikes = numspikes;
                Results(sweep).peak_times = peak_times;
                Results(sweep).peak_idx = loc;
                Results(sweep).peak_amplitudes = pks;
                Results(sweep).peak_to_baseline = peak_amps;
                
                %% numspikes 1 sec and 150ms
                Results(sweep).spikes1sec_idx = loc(loc <= stop1sec);
                Results(sweep).spikes1sec_time = loc(loc <= stop1sec)*(si/1000);
                Results(sweep).spikes1sec_amp = pks(loc <= stop1sec);
                Results(sweep).num_spikes1sec = length(find(loc <=stop1sec));
                
%                 Results(sweep).spikes150ms_idx = loc(loc <= stop150ms);
%                 Results(sweep).spikes150ms_time = loc(loc <= stop150ms)*(si/1000);
%                 Results(sweep).spikes150ms_amp = pks(loc <= stop150ms);
%                 Results(sweep).num_spikes150ms = length(find(loc <= stop150ms));
        
% % % % % %                 wind = 10/(si/1000); % look into 10ms past peak points
% % % % % %                 FastAHP = [];
% % % % % %                 for i = 1:numspikes
% % % % % %                     pklocation = loc(i);
% % % % % %                     if length(data) < (pklocation+wind)
% % % % % %                         stoppage = length(data);
% % % % % %                     else
% % % % % %                         stoppage = pklocation+wind;
% % % % % %                     end
% % % % % %                     [FastAHP(i,1) FastAHP(i,2)] = min(data(pklocation:stoppage));
% % % % % %                     FastAHP(i,2) = FastAHP(i,2) + pklocation;
% % % % % %                 end
% % % % % %                 Results(sweep).FastAHP_Voltage = FastAHP(:,1);
% % % % % %                 Results(sweep).FastAHP_Time = FastAHP(:,2)*(si/1000);
% % % % % %                 Results(sweep).FastAHP_Baseline = FastAHP(:,1) - Results(sweep).baseline_potential;
%                 %% Threshold
%                 [Results(sweep).threshold_time, Results(sweep).threshold_amplitude,Results(sweep).threshold_index,...
%                     Results(sweep).dvdtthreshold,Results(sweep).dvdt1,Results(sweep).dvdt2] = ...
%                     SpikeThreshold(data,dvdtthreshold,si,numspikes,xminmax,sweep,0,stop,loc);
%                 % plot the results with the fast AHP and peaks as well
%                 Results(sweep).threshold_baseline = Results(sweep).threshold_amplitude - Results(sweep).baseline_potential;
%                 
%                 %%%%%%%%%%%%%%%%%%%%add FastAHP here %%%%%%%%%%%%%%%%%%
%                 [FastAHP] = FastAHPfinder(data,si,Results(sweep).threshold_index);
%                 Results(sweep).FastAHP_Voltage = FastAHP(:,1);
%                 Results(sweep).FastAHP_Time = FastAHP(:,2)*(si/1000);
%                 Results(sweep).FastAHP_Baseline = FastAHP(:,1) - Results(sweep).baseline_potential;
%                 
                if isempty(loc(loc <= stop1sec))
                    error('Spikes in sweep but not in interval of interest. Check baseline interval so baseline offset corresponds to stimulus onset')
                end
                if graph_on == 1
                    fff = figure;
                    plot(1:length(data),data,'b',loc,pks,'r*')
                    hold on
                    plot(loc(loc <= stop1sec),pks(loc <= stop1sec),'c*')
%                     plot(loc(loc <= stop150ms),pks(loc <= stop150ms),'g*')
                    
                    
%                     plot(1:length(data),data,'b',Results(sweep).threshold_index,Results(sweep).threshold_amplitude,'m*')
                    xlim([xminmax])
%                     plot(FastAHP(:,2),FastAHP(:,1),'g.')
                    xlabel(sprintf('Points (multiply by %g for ms)',si/1000 ))
                    ylabel('Volts')
%                     movegui('northwest')
                    title(sprintf('Sweep %g File:%s',sweep,filename))
                end
%                 %% Spikewidth from Baseline
%                 for i = 1:numspikes
%                     threshold = (peak_amps(i)/2) + Results(sweep).baseline_potential;
%                     pklocation = loc(i);
%                     [Results(sweep).SpikeWidth_Baseline(i,1), ~, ~] = halfwidth(threshold,pklocation,data,si);
%                 end     
%                 %% Spikewidth from Threshold
%                 for i = 1:numspikes
%                     pklocation = loc(i); % need to go back to points 
%                     if ~isnan(Results(sweep).threshold_amplitude(i))
%                         threshold = ((pks(i) - Results(sweep).threshold_amplitude(i))/2) + Results(sweep).threshold_amplitude(i);
%                         [Results(sweep).SpikeWidth_Threshold(i,1), ~,~] = halfwidth(threshold, pklocation,data,si);
%                     else
%                         Results(sweep).SpikeWidth_Threshold(i,1) = NaN;
%                     end
%                 end    
%                 %% Spikewidth from First Spikewidth-Baseline 
%                 firstthreshold = (peak_amps(1)/2) + Results(sweep).baseline_potential;
%                 for i = 1:numspikes
%                     pklocation = loc(i); % need to go back to points 
%                     [Results(sweep).SpikeWidth_FirstSpike(i,1), ~,~] = halfwidth(firstthreshold, pklocation,data,si);
%                 end    
%                 
                %% Calculate ISI
                if numspikes > 1
                    Results(sweep).ISI = diff(peak_times) ;
                else
                    Results(sweep).ISI = NaN;
                end
                
                if Results(sweep).num_spikes1sec > 1
                    Results(sweep).ISIspikes1sec = diff(Results(sweep).spikes1sec_time);
                else
                    Results(sweep).ISIspikes1sec = NaN;
                end

%                 if Results(sweep).num_spikes150ms > 1
%                     Results(sweep).ISIspikes150ms = diff(Results(sweep).spikes150ms_time);
%                 else
%                     Results(sweep).ISIspikes150ms = NaN;
%                 end
                

%                 if sweep > numsweeps-3
                    starthz = stop+150/(si/1000);
                    duration = (stop1sec - starthz)*(si/1000);
                    numspk = length(find(loc >= starthz+1 & loc <= stop1sec));
                    hz = numspk/(duration/1000);
                    if hz <= 5 % <---- 5 hz cutoff for accommodation
                        accommodating = [accommodating; 1];
                    else
                        accommodating = [accommodating; 0];
                    end
                    Results(sweep).sweepaccommodating = accommodating;
                    Results(sweep).hz_150 = hz;
%                 end
                
 %%
            else
                Results(sweep).num_spikes = NaN;
                Results(sweep).peak_times = NaN;
                Results(sweep).peak_idx = NaN;
                Results(sweep).peak_amplitudes = NaN;
                Results(sweep).peak_to_baseline = NaN;
%                 Results(sweep).FastAHP_Voltage = NaN;
%                 Results(sweep).FastAHP_Time = NaN;
%                 Results(sweep).FastAHP_Baseline = NaN;
%                 Results(sweep).threshold_time = NaN;
%                 Results(sweep).threshold_amplitude = NaN;
%                 Results(sweep).threshold_index = NaN;
%                 Results(sweep).threshold_baseline = NaN;
%                 Results(sweep).dvdtthreshold = NaN;
%                 Results(sweep).dvdt1 = NaN;
%                 Results(sweep).dvdt2= NaN;
%                 Results(sweep).SpikeWidth_Baseline = NaN;
%                 Results(sweep).SpikeWidth_Threshold = NaN;
%                 Results(sweep).SpikeWidth_FirstSpike = NaN;
                Results(sweep).ISI = NaN;
                
                Results(sweep).spikes1sec_idx = NaN;
                Results(sweep).spikes1sec_time = NaN;
                Results(sweep).spikes1sec_amp = NaN;
                Results(sweep).num_spikes1sec = NaN;
                
%                 Results(sweep).spikes150ms_idx = NaN;
%                 Results(sweep).spikes150ms_time = NaN;
%                 Results(sweep).spikes150ms_amp = NaN;
%                 Results(sweep).num_spikes150ms = NaN;
                
%                 if sweep > numsweeps-3
                    accommodating = [accommodating; NaN];
                    Results(sweep).sweepaccommodating = accommodating;
                    Results(sweep).hz_150 = NaN;
%                 end
            end

    end
    
    %% Is the file Accommodating?
    % look at the last 3 sweeps
    if accommodating(numsweeps-2:end) == 1
        accommresult = 1;
        sprintf('File %s is accommodating',filename)
    elseif accommodating(numsweeps-2:end) == 0
        accommresult = 0;
        sprintf('File %s is non-accommodating',filename)
    else
        accommresult = NaN;
        sprintf('Last 3 sweeps are inconsistent')
        
    end
    for sweep = 1:numsweeps
        Results(sweep).fileaccommodating = accommresult;
    end
    
end