bachelor_thesis/DataAnalysisToolbox/spikeTrigAvg_postsyn.m

function varargout = spikeTrigAvg_postsyn(files, varargin)
%% Calculate spike triggered average of the postsynaptic response
% For each spike in the stimulated cell the membrane potential of the not
% stimulated cell is gathered and summed up from 5 ms before to 30 ms after
% the spike. Devided by the number of spikes the averaged postsynaptic
% membrane potential is obtained. If there is a clear peak (due to EPSPs 
% for example), the time of this peak is determined as the latency of the 
% reaction to the spike in the other cell. The latency and the number of 
% used spikes are returned. Alternatively the result can be plotted. The 
% function can handle multiple files at a time.
% -------------------------------------------------------------------------
% Input     [additional arguments are positional, 
% -----      optional arguments are name-value pairs]
%
%   file:           number(s) of the '.mat' file(s) in the directory,
%                       type: integer or vector
%   trials:         number(s) of the trials to be analysed, type: integer
%   (additional)        of vector, default: 'all'
%   Filter:         select between three filter options: 
%   (optional)          'hn'    applys a hum-noise-remove filter
%                       'hnhfn' additonally applys a low pass filter to
%                              remove high frequency noise
%                       If not set, no filter will be applyed (default).
%   MinPeakPromT:   set the value for 'MinPeakProminence' which is used in
%   (optional)          the 'findpeaks' function to find T-cell peaks. 
%                       Default: 17
%   StimCell:       set whether the data for a stimulated T cell or a 
%   (optional)          stimulated interneuron should be used.
%                       'INT' for the interneuron
%                       'T'   for the T-cell (default)
% ------
% Output
% ------
%   varargout{1}:   latency from spike of the stimulated cell to peak in 
%                   average response of the not stimullated cell (no plot)
%                       type: float (single file) or vector (multiple
%                             files)
%   varargout{2}:   number of spikes used from the stimulated cell
%                       type: float (single file) or vector (multiple
%                             files)
%
%   visually output of the plot, if no output variable set
% -------------------------------------------------------------------------
% Program by: Bjarne Schultze [last modified 02.03.2022]
% -------------------------------------------------------------------------

%% Preparations
% parse the user input
[file_nums, trials, filter_choice, ~, MinPeakPromT, ...
    stimCell, ~,~,~] = parse_input_peakfun(files, varargin{:});
% extract the requested data
recData = extract_data(file_nums, filter_choice);

% pre-define vectors for the output (optional)
if nargout > 0
    latencies = zeros(length(file_nums),1);
    spikesNum = zeros(length(file_nums),1);
end

% iterate through the single files as requested by the input
for file = 1:length(file_nums)
    % get the number or trials for the current file
    trial_count = size(recData{file,1}.recordingT,2);
    % set the length of the time frame to gather the triggered data
    trig_len = 0.03;                     % in seconds
    trig_len_dp = trig_len * 10000;      % in datapoints
    
    % distinguish between a stimulated T cell and a stimulated interneuron
    if isequal(stimCell, "INT")
        % get the columns with data for a stimulated interneuron
        stimulated = setdiff(1:trial_count, recData{file,1}.stimulatedT(:), ...
            'stable');
        % select the requested trials, ensure only valid trials
        if ~isequal(trials, 'all')
            stimulated = intersect(stimulated,trials);
        end
        % set an appropriate plot title
        plot_title = sprintf("Postsynaptic Spike Response - File: " + ...
            "%02.f - Interneuron stimulated", file_nums(file));
        % assign the data to working variables
        recDataNSTIM = recData{file,1}.recordingT(:,stimulated);
        recDataSTIM = recData{file,1}.recordingINT(:,stimulated);
    else
        % get the columns with data for a stimulated T-cell
        stimulated = recData{file,1}.stimulatedT(:);
        % select the requested trials, ensure only valid trials
        if ~isequal(trials, 'all')
            stimulated = intersect(stimulated,trials);
        end
        % set an appropriate plot title
        plot_title = sprintf("Postsynaptic Spike Response - File: " + ...
            "%02.f - T-cell stimulated", file_nums(file));
        % assign the data to working variables
        recDataSTIM = recData{file,1}.recordingT(:,stimulated);
        recDataNSTIM = recData{file,1}.recordingINT(:,stimulated);
        MinPeakPromSTIM = MinPeakPromT;
    end
    
    %% Main calculations
    % pre-define a variable to sum up the triggered membrane potentials
    triggered_sum = 0;                                         
    spike_counter = 0;
    % iterate through the trials
    for exp = 1:length(stimulated)
        % search for spikes in the stimulated cell
        if isequal(stimCell, 'INT')
            [~,spike_indexSTIM_all] = findIntSpikes(recDataSTIM(:,exp));
        else
            [~, spike_indexSTIM_all] = findpeaks(recDataSTIM(:,exp), ...
                'MinPeakProminence', MinPeakPromSTIM);
        end
        % exclude spikes at the very beginning and the end of one recording
        % to ensure a complete time frame for the triggered data
        spike_indexSTIM_sel = spike_indexSTIM_all > 50 & ...
            spike_indexSTIM_all < (size(recDataSTIM,1) - trig_len_dp);
        spike_indexSTIM = spike_indexSTIM_all(spike_indexSTIM_sel);
        
        % if no spikes are found the next trial is examined
        if isempty(spike_indexSTIM)
            continue
        else
            % select the membrane potentail of the not-stimulated cell in 
            % an index frame of the length 'trig_len_dp' datapoints for 
            % each spike of the stimulated cell
            for spike = 1:length(spike_indexSTIM)
                % select the data (triggered by a spike)
                triggered_data = ...
                    recDataNSTIM((spike_indexSTIM(spike)-50):...
                    (spike_indexSTIM(spike) + trig_len_dp), exp);
                % sum up the data
                triggered_sum = triggered_sum + triggered_data;
            end
            % count the spikes
            spike_counter = spike_counter + spike;
        end
    end
    % calculate the averaged postsynaptic membrane potential if more than 8
    % spikes were found in the T cell
    if spike_counter < 8
        afterSpkAvg = NaN;
        fprintf("File %g: Not enough T cell spikes (n < 8)!\n", ...
            file_nums(file))
    else
        afterSpkAvg = triggered_sum ./ spike_counter;
    
        % create a time vector for the selected data
        timeVector = [sort(-recData{file,1}.timeVector(1:50)); 0;...
            recData{file,1}.timeVector(1:trig_len_dp)];
        timeVector = timeVector * 1000;     % in ms
        
        % find the peaks in the average spike response if there are some
        [peaks, locs] = findpeaks(afterSpkAvg, timeVector, ...
            'MinPeakProminence',0.05);
        
        % exclude peaks with locations before zero
        sel_peaks = peaks(locs > 0); sel_locs = locs(locs > 0);
        
        % select peaks that are by at least 0.2 mV bigger than value at 0
        search_peaks = sel_peaks > (afterSpkAvg(51)+0.2);

        % subset the peak and corresponding location values
        sel_peaks = sel_peaks(search_peaks); 
        sel_locs = sel_locs(search_peaks);
        % select the peak with the location closest to zero
        % firstPeak_l is giving the latency (delay of the reaction to the
        % spike in the stimulated cell)
        [firstPeak_l,index] = min(sel_locs); 
        firstPeak_v = sel_peaks(index);
    end
    
    %% Plotting of the result and output
    % suitable output based on whether calculations were made (n >= 8) and
    % based on the number of output variables
    if nargout > 0 && spike_counter >=8
        % return the latency (firstPeak_l) if determined
        if isempty(firstPeak_l)
            latencies(file,1) = NaN;
        else
            latencies(file,1) = firstPeak_l;
        end
        % return the number of spikes used to the second output variable
        if nargout > 1
            spikesNum(file,1) = spike_counter;
        end
    elseif spike_counter >= 8
        % plot the averaged postsynaptic membrane potential
        figure();
        plot(timeVector, afterSpkAvg);
        title(plot_title);
        xlabel("time after spike [ms]"); 
        ylabel("average membrane potentail [mV]");
        % add a reference line at zero
        xline(0, 'Color', [0.6,0.6,0.6]);
        
        % add the distance from zero to the measured peak if determined
        if ~isempty(firstPeak_v)
            % set the coordinates for line and text
            distanceLine = [0,firstPeak_l]; 
            text_X = firstPeak_l/2; text_Y = firstPeak_v*0.997;
            % add the text
            text_T = sprintf("%g ms", firstPeak_l);
            text(text_X, text_Y, text_T, 'FontSize', 12, 'Color', ...
                [0.6,0.6,0.6]);
            % add the distance line
            hold on;
            plot(distanceLine, [firstPeak_v,firstPeak_v], '|-', ...
                'Color', [0.6,0.6,0.6])
            hold off;
            % adjust the limits of the y-axis
            max_data = max(afterSpkAvg);
            min_data = min(afterSpkAvg);
            range_data = abs(max_data-min_data);
            ylim([(min_data-range_data*0.1),(max_data+range_data*0.1)]);
        end
    else
        % give back NaN for no result
        latencies(file,1) = NaN;
    end
% end of outer for loop (iteration through files) 
end

% return the output values if requested
if nargout > 0
    varargout{1} = latencies;
    if nargout > 1
        varargout{2} = spikesNum;
    end
end


% end of function definition
end