bachelor_thesis/DataAnalysisToolbox/ampAtStim.m

function varargout = ampAtStim(files, varargin)
%% Function to obtain the neurons amplitude during a stimulation
% The function calculates the amplitude for both cells during a stimulation
% out of the perspective of the stimulated cell. The stimuli for which the 
% calculation should be performed can be set, as well as which cell should 
% be the stimulated one. The amplitudes can either be plotted against the
% amplitude of the current pulses (with error bars if needed) or the values
% are returned. The function can handle multiple files at a time.
% -------------------------------------------------------------------------
% Input     [additional arguments are positional, 
% -----      optional arguments are name-value pairs]
%
%   files:      ID-number(s) of the '.mat' data-file(s), 
%                   Type: integer or vector
%   trials:     number(s) of the trials to be analysed, type: integer or
%   (additional)    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, it will be set to "no" (default).
%   StimCell:   set whether to use the T cell ('T') or the interneuron 
%   (optional)      ('INT') as the stimulated cell. Default: 'T'
%   ErrorBars:  set whether error bars (standard deviation) should be 
%   (optional)      plotted (1 or 0). Default: 1 (yes)
%   Stimuli:    set for which stimuli the data should be analyzed, 
%   (optional)      Type: integer or vector or 'all', default: 'all'
% ------
% Output
% ------
%   varargout{1}: structure(s), containing the calculated values wrapped in
%                   a cell array, amplitudes for the stimuli are sorted 
%                   from smallest stimulus to the biggest stimulus
%                    - MeanAmplitudeT/MeanAmplitudeINT: averaged amplitudes
%                      for the trials where the 'stimCell' was stimulated
%                    - AmplitudeT/AmplitudeINT: amplitudes for the cells
%                      for all trials (not averaged)
%                    - sdT/sdINT: standard deviation for all amplitudes   
%                      over all trials for each stimulus
%                    - stimulatedT: trials with a stimulated T cell
%
%   visually output of the plot(s), if no output variable is set, new 
%   figure window(s) is/are created
% -------------------------------------------------------------------------
% Program by: Bjarne Schultze [last modified 02.03.2022]
% -------------------------------------------------------------------------

%% Preparations
% parse the user input
[file_nums, trials, filter_choice, tol, ~, stimCell, error_bars, ...
    stimuli, ~] =  parse_input_peakfun(files, varargin{:});

% preallocate a cell array for the output of the calculated data
fun_output = cell(length(file_nums),1);

% get the necessary data with the extract_data function
recData = extract_data(file_nums, filter_choice);

%% Determine the stimuli onsets
% find the onsets for the positive stimuli
[stim_peaks, stim_locs] = findpeaks(recData{1,1}.stimulus, 10000);

% find the onsets for the negativ stimuli (therefor inverting the
% stimulus data)
[min_stim_peaks, min_stim_locs] = ...
    findpeaks(-recData{1,1}.stimulus, 10000, 'MinPeakHeight', 0);
% invert the peak values back
min_stim_peaks = -min_stim_peaks;

% concatenating positive and nagative stimlus peak data
all_stim_locs = [min_stim_locs; stim_locs];
all_stim_peaks = [min_stim_peaks; stim_peaks];
% sort the values for the current in ascending order
[all_stim_peaks_sort, sort_order] = sort(all_stim_peaks);
all_stim_locs_sort = all_stim_locs(sort_order);
% select single stimuli if requested
if ~isequal(stimuli, 'all')
    [~, stim_indices] = intersect(all_stim_peaks_sort, stimuli);
    all_stim_locs_sort = all_stim_locs_sort(stim_indices);
    all_stim_peaks_sort = all_stim_peaks_sort(stim_indices);
end
  
%% Calculation of the amplitude and its standard deviation
% repeat all calculations for every file requested by the user-input
for file = 1:length(file_nums)
    timeVector = recData{file,1}.timeVector;

    % length of stimulus (it is the same for all stimuli)
    stim_len = 0.5;                         % in ms
    stim_num = length(all_stim_locs_sort);  % number of stimuli
    
    trial_num = size(recData{file,1}.recordingT,2);

    % pre-define vectors for the results of the calculations
    ampT = zeros(stim_num,trial_num);
    ampINT = ampT;
    
    % iterate through all trials of an experiment
    for exp = 1:trial_num
        % iterate through the stimuli calculating amplitude and sd for both 
        % cells
        for stim = 1:length(all_stim_locs_sort)
            % select the time frame in which one cell was stimulated +
            % tolerance
            search_index_amp = timeVector >= all_stim_locs_sort(stim) ...
                + tol & timeVector <= all_stim_locs_sort(stim) + stim_len;
            % select a time frame before the stimulus onset with the same 
            % length as the stimulus for a reference
            search_index_preamp = timeVector <= all_stim_locs_sort(stim) ...
                & timeVector >= all_stim_locs_sort(stim) - stim_len + tol;
    
            % calculate the amplitude of the membrane potential in the T
            % cell in relation to the amplitude in the reference time frame
            ampT(stim,exp) = ...
                mean(recData{file,1}.recordingT(search_index_amp,...
                exp)) - mean(recData{file,1}.recordingT(...
                search_index_preamp,exp));
            % if the amplitude is smaller than 3 sd in the reference time
            % frame it is considered as zero
            if abs(ampT(stim,exp)) < ...
                    3 * std(recData{file,1}.recordingT...
                    (search_index_preamp,exp))
                ampT(stim,exp) = 0;
            end
    
            % same calculations for the interneuron
            ampINT(stim,exp) = ...
                mean(recData{file,1}.recordingINT(search_index_amp,...
                exp)) - mean(recData{file,1}.recordingINT(...
                search_index_preamp,exp));
            if abs(ampINT(stim,exp)) < ...
                    3 * std(recData{file,1}.recordingINT...
                    (search_index_preamp,exp))
                ampINT(stim,exp) = 0;
            end
        end
    end
    
%% Averaging of the amplitudes
    % pre-define vectors for the mean and standard deviation of the
    % amplitudes over all trials
    m_ampT = zeros(stim_num,1); m_ampINT = zeros(stim_num,1);
    sd_ampT = zeros(stim_num,1); sd_ampINT = zeros(stim_num,1); 
    % distinguish which cell should be the stimulated
    if isequal(stimCell, "INT")
        if isequal(trials,'all')
            stimulated = setdiff(1:trial_num, ...
                recData{file,1}.stimulatedT(:), 'stable');
        else
            stimulated = setdiff(trials, ...
                recData{file,1}.stimulatedT(:), 'stable');
        end
        plot_title = sprintf("File: %02.f - Interneuron stimulated", ...
            file_nums(file));
    else
        stimulated = recData{file,1}.stimulatedT;
        plot_title = sprintf("File: %02.f - T cell stimulated", ...
            file_nums(file));
        if ~isequal(trials,'all')
            stimulated = intersect(trials,1:trial_num);
        end
    end
    
    % calculate the average and standard deviation of the amplitude for 
    % each stimulus
    for stim = 1:stim_num
        % mean and sd of trials where the T-cell was stimulated
        m_ampT(stim) = mean(ampT(stim,stimulated));
        sd_ampT(stim) = std(ampT(stim,stimulated));
        
        % mean and sd for the interneurons
        m_ampINT(stim) = mean(ampINT(stim,stimulated));
        sd_ampINT(stim) = std(ampINT(stim,stimulated));
    end

    % define a threshold/minimum mean amplitude of 0.1 mV, everything 
    % beneath is set to zero
    under_thresholdT = abs(m_ampT) < 0.1;
    m_ampT(under_thresholdT) = 0; sd_ampT(under_thresholdT) = 0;
    under_thresholdINT = abs(m_ampINT) < 0.1;
    m_ampINT(under_thresholdINT) = 0; sd_ampINT(under_thresholdINT) = 0;
    
%% Return and/or plot the results 
    % if one output is requested the plot will not be shown
    if nargout > 0
        % gather the calculated data in a structure for the output
        output_data = struct('MeanAmplitudeT',m_ampT, ...
                            'MeanAmplitudeINT', m_ampINT, ...
                            'AmplitudeT', ampT, ...
                            'AmplitudeINT', ampINT, ...
                            'sdT', sd_ampT, 'sdINT', sd_ampINT, ...
                            'stimulatedT', recData{file,1}.stimulatedT);
        fun_output{file,1} = output_data;
        % assign the calculated data to output variable 1
        varargout{1} = fun_output;
    else
        % create a new figure window for a visible plot
        figure();
    end
    
    if isequal(nargout, 0) || nargout > 1
        % plot the results either with or without error bars
        if isequal(error_bars, 1)
            errorbar(all_stim_peaks_sort, m_ampT, sd_ampT);
            axis padded
            hold on;
            errorbar(all_stim_peaks_sort, m_ampINT, sd_ampINT);
            hold off;
            title(plot_title);
            xlabel("current [nA]"); ylabel("amplitude [mV]");
            legend("T-cell", "Interneuron", 'Location', 'best', ...
                'AutoUpdate', 'off');
            yline(0, 'Color', [0.6,0.6,0.6]);
            subtitle("- error bars representing standard deviation between trials -");
        else
            plot(all_stim_peaks_sort, m_ampT, 'o-',...
                all_stim_peaks_sort, m_ampINT, 'o-', 'MarkerSize', 4);
            axis padded
            title(plot_title);
            xlabel("stimulation current [nA]"); ylabel("amplitude [mV]");
            legend("T-cell", "Interneuron", 'Location', 'best',...
                'AutoUpdate', 'off');
            yline(0, 'Color', [0.6,0.6,0.6]);
        end
    end
    
% end of the iteration through the single files
end

%% end of function definition
end