bachelor_thesis/DataAnalysisToolbox/transJuncVolt.m

function transJuncVolt(files, varargin)
%% Analyse the effect of the transjunctional voltage on the electrical coupling
% This function plots the ratio of the amplitudes of two cells during a 
% stimulation in relation to the differnce in their resting membrane
% potentials. The latter is an estimation for the transjunctional voltage. 
% The plot helps to judge whether the transjunctional voltage has an effect 
% on the electrical coupling. If requested the trial numbers can be added 
% to check for other effects as well. The function can handle up to seven 
% stimuli at the same time. Alternatively up to 11 files can be analysed at 
% a time, but then only one stimulus is used.
% -------------------------------------------------------------------------
% Input     [optional arguments are name-value pairs]
% -----
%       files:      number of the dataset to examine
%       Filter:     select between no filter (default), hum noise ('hn') or
%       (optional)  hum and high frequency noise ('hnhfn')
%       StimCell:   set which cell should be the stimulated one ('T' or
%       (optional)  'INT')
% ------
% Output
% ------
%       graphic output, scatter plot, new figure window
% -------------------------------------------------------------------------
% Program by: Bjarne Schultze [last modified 17.12.2021]
% -------------------------------------------------------------------------

%% Preparations 
% select the stimuli to take for the calculations by their amperage (up to
% seven stimuli are possible at the same time)
stimAmps = [-2,-1]; 

% set whether the points in the scatter plot should be marked with the
% trial numbers
labelling = 1;

% set whether each file should get its own plot (0) or the results for all
% files are plottet in one plot (1)
onePlot = 1;
% with multiple files in one plot, only one stimulus is used
if onePlot
    stimAmps = stimAmps(1);
end

% parse the user input
[file_nums, ~, filter_choice, ~, ~, stimCell, ~,~,~] = ...
    parse_input_peakfun(files, varargin{:});
% extract the requested data
recordings = extract_data(file_nums, filter_choice);
% load a file containing the electrode resistance values
load('../AdditionalFiles/electrodeProps.mat','electrodeProps')

% preallocate a string array for the resistance output
d_resistance_output(1,1) = ...
    sprintf("\n\x0394 resistance values (StimCell %s)\n" + ...
    "------------------------\nFile\tT cell\t\tInt\n" + ...
    "------------------------\n", stimCell);

% list of all the amperages of the single stimuli in ascending order
amps_rows = [-2,-1,-0.5,0.03125,0.0625,0.125,0.25,0.5,0.75,1,1.25,1.5];

%% Main calculations
% getting the amplitudes for the both cells during stimulation
amplitudes = ampAtStim(files, varargin{:});

for file = 1:length(files)
    % select the dataset for the current iteration
    recData = recordings{file,1};
    
    % preallocate matrices to collect the calculated data
    stimNum = size(stimAmps,2); trialNum = size(recData.recordingT,2);
    search_StimAmp = zeros(stimNum,1);
    rmpT = zeros(stimNum,trialNum);
    rmpINT = zeros(stimNum,trialNum);

    % select and average 500 ms (5000 data points) before each stimulus as RMP
    for i = 1:stimNum
        % search the list of amperages for the given stimulus value
        search_StimAmp(i) = find(abs(amps_rows-stimAmps(i)) < 0.001);
        % find the times (datapoints) of the beginning of the required
        % stimulus
        dp_stimOnset = find((abs(recData.stimulus - stimAmps(i)) < 0.001),1);

        % select the recorded data in this time frames for all trials
        sel_dataT = recData.recordingT((dp_stimOnset-5000):...
            dp_stimOnset,:);
        sel_dataINT = recData.recordingINT((dp_stimOnset-5000):...
            dp_stimOnset,:);
        % average the membrane potentials within each time frame
        rmpT(i,:) = mean(sel_dataT);
        rmpINT(i,:) = mean(sel_dataINT);
    end

    % distinguish between a stimulated T cell and a stimulated interneuron
    if isequal(stimCell, 'INT')
        stimulated = setdiff(1:trialNum, recData.stimulatedT(:), 'stable');
        % select the RMP according to the stimulation information
        rmpNStim = rmpT(:,stimulated); rmpStim = rmpINT(:,stimulated);
        % select the amplitudes respectively 
        amps_NStim = ...
            amplitudes{file,1}.AmplitudeT(search_StimAmp,stimulated);
        amps_Stim = ...
            amplitudes{file,1}.AmplitudeINT(search_StimAmp,stimulated);
        % define a plot title
        plot_subtitle = ("- interneuron stimulated -");
    else
        stimulated = recData.stimulatedT;
        % select the RMP according to the stimulation information
        rmpStim = rmpT(:,stimulated); rmpNStim = rmpINT(:,stimulated);
        % select the amplitudes respectively 
        amps_Stim = ...
            amplitudes{file,1}.AmplitudeT(search_StimAmp,stimulated);
        amps_NStim = ...
            amplitudes{file,1}.AmplitudeINT(search_StimAmp,stimulated);
        % define a plot title
        plot_subtitle = ("- T cell stimulated -");
    end

    % calculate the difference in the RMPs
    rmp_diff = rmpStim - rmpNStim;

    % calculate the ratio between the amplitudes of the two cells
    amp_ratio = amps_NStim ./ amps_Stim;
    
    % calculate the changes in electrode resistances
    c_fileNum = files(file);
    if isequal(electrodeProps.cell_left(c_fileNum), 'T')
        d_resistance_T = electrodeProps.R_left_post(c_fileNum)-...
            electrodeProps.R_left(c_fileNum);
        d_resistance_INT = electrodeProps.R_right_post(c_fileNum)-...
            electrodeProps.R_right(c_fileNum);
    elseif isequal(electrodeProps.cell_left(c_fileNum),'I')
        d_resistance_INT = electrodeProps.R_left_post(c_fileNum)-...
            electrodeProps.R_left(c_fileNum);
        d_resistance_T = electrodeProps.R_right_post(c_fileNum)-...
            electrodeProps.R_right(c_fileNum);
    end
    % create a command window output for the delta resistance values
    if isnan(d_resistance_T)
        d_resistance_output(file+1,1) = sprintf("%02.f \t\t%02.f \t\t%02.f\n",...
        c_fileNum, d_resistance_T, d_resistance_INT);
    else
    d_resistance_output(file+1,1) = sprintf("%02.f \t\t%02.f \t\t\t%02.f\n",...
        c_fileNum, d_resistance_T, d_resistance_INT);
    end

    %% Show results
    % plot the results as a scatter plot
    if onePlot
        if file == 1
            figure();
            % obtain the maximum and minimum values for axes scaling
            maxX = max(rmp_diff, [], 'all');
            minX = min(rmp_diff, [], 'all');
            maxY = max(amp_ratio, [], 'all');
            minY = min(amp_ratio, [], 'all');
            ax = gca;
            % add new colors if more than 7 files are requested
            if length(files) > 7
                ax.ColorOrder((end+1:end+4),:) = [0.67 0.53 0.22; ...
                                                  0    0.56 0   ; ...
                                                  0.69 0.81 0.15; ...
                                                  0.83 0.06 0.06];
            end
        end
        % use the same color for all datasets if more than 11 are requested
        if length(files) > 11
            current_color = ax.ColorOrder(1,:);
        else
            current_color = ax.ColorOrder(file,:);
        end
        % plot the results of the current file
        hold on;
        plot(rmp_diff', amp_ratio', 'o', 'Color', current_color,...
            'MarkerSize', 5, 'MarkerFaceColor', current_color);
        hold off;
        % gather max and min values for axes scaling
        if max(rmp_diff, [], 'all') > maxX
            maxX = max(rmp_diff, [], 'all');
        end
        if min(rmp_diff, [], 'all') < minX
            minX = min(rmp_diff, [], 'all');
        end
        if max(amp_ratio, [], 'all') > maxY
            maxY = max(amp_ratio, [], 'all');
        end
        if min(amp_ratio, [], 'all') < minY
            minY = min(amp_ratio, [], 'all');
        end
    else
        % plot the results in one plot per file
        figure();
        plot(rmp_diff', amp_ratio', 'o');
    end

    % label the points with their trial number, if requested and if the 
    % results for each file are plotted  in a seperate figure
    if labelling && ~onePlot
        % set the points to place the labels 
        textX = rmp_diff';
        max_amp_ratio = max(amp_ratio, [], 'all');
        min_amp_ratio = min(amp_ratio, [], 'all');
        range_amp_ratio = abs(max_amp_ratio - min_amp_ratio);
        textY = (amp_ratio + range_amp_ratio * 0.04)';
        
        % create the labels
        textNums = num2str(stimulated');
        % get the current axes object to access the color order
        ax = gca;
        % label the datapoints
        for i = 1:stimNum
            text(textX(:,i), textY(:,i), textNums, 'FontSize', 7, ...
                'Color', ax.ColorOrder(i,:));
        end
        % re-scale the axes to have all numbers inside the plot
        if maxY ~= 0 && minY ~= 0
            ylim([(min_amp_ratio - range_amp_ratio * 0.1), ...
                (max_amp_ratio + range_amp_ratio * 0.1)]);
        end
        max_rmp_diff = max(rmp_diff, [], 'all');
        min_rmp_diff = min(rmp_diff, [], 'all');
        range_rmp_diff = abs(max_rmp_diff - min_rmp_diff);
        xlim([(min_rmp_diff - range_rmp_diff * 0.05),...
            (max_rmp_diff + range_rmp_diff * 0.05)]);
    end
    % add title and axes labels
    ylabel("ratio of amplitudes (NStim/Stim)"); 
    xlabel("transjunctional voltage (Stim-NStim) [mV]");
    title(sprintf("Effect of transjunctional voltage - File: %0.f", ...
        files(file)));
    subtitle(plot_subtitle);
    
    % add a suitable legend if single plots were created
    if ~onePlot
        leg_entries = "";
        for i = 1:stimNum
            leg_entries(i) = sprintf("%g nA stimulus", stimAmps(i));
        end
        legend(leg_entries, 'Location','best');
    end
% end of file-iteration-loop
end

if onePlot
    % rescale axes if multiple files in one plot
    rangeX = abs(maxX - minX); rangeY = abs(maxY - minY);
    xlim([(minX - 0.05*rangeX), (maxX + 0.05*rangeX)]);
    if minY ~= 0 && maxY ~= 0
        ylim([(minY - 0.1*rangeY), (maxY + 0.1*rangeY)]);
    end
    
    % add a suitable legend
    if length(files) < 12
        leg_entries = "";
        for i = 1:length(files)
            leg_entries(i) = sprintf("File: %g // %g nA stim",...
                files(i),stimAmps(1));
        end
        legend(leg_entries, 'Location', 'best');
    end
    title("Effect of transjunctional voltage");
end

% give back the delta resistance values to the command window
fprintf("%s", d_resistance_output);

%% end of function definition
end