% GJLFP, eLife 2026. Richa Sirmaur and Rishikesh Narayanan.
% Figure 1: waveforms (A-B),boxplots active vs passive (C-D), waveforms across entire neuron (E)

%% Load dataset and bandpass filter EFP signals

% loading the EFP data from all 7 arrays, all having 49 electrodes
dataFiles = {
    'S1_allLFP_E0-E48.xlsx'
    'S2_allLFP_E0-E48.xlsx'
    'S3_allLFP_E0-E48.xlsx'
    'S4_allLFP_E0-E48.xlsx'
    'S5_allLFP_E0-E48.xlsx'
    'S6_allLFP_E0-E48.xlsx'
    'S7_allLFP_E0-E48.xlsx'
};
dataRange = 'A1:AW80000'; 
nFiles = numel(dataFiles); % seven

% FILTER PARAMETERS ----------------
fs  = 40000;   % sampling frequency (Hz)
fc1 = 0.5;     % low cutoff (Hz)
fc2 = 300;     % high cutoff (Hz)
Wn = [fc1 fc2] / (fs/2);   % cutoff
n  = 2;                    % filter order
[b,a] = butter(n, Wn, 'bandpass');

% LOAD + FILTER ----------------
lfpData = [];
for f = 1:nFiles
    Mat = xlsread(dataFiles{f}, dataRange); % Load dataset
    lfpMat = zeros(size(Mat));     % filtered matrix
    % Filter each electrode trace
    for i = 1:size(Mat,2)
        lfpMat(:,i) = filtfilt(b, a, Mat(:,i)); %forward-backward filtering for no time shift
    end
    lfpData{f} = lfpMat;
end
filteredmat = cat(2, lfpData{:}); % concatenating all arrays column-wise (7 arrays each with 49 electrodes) in one matrix
dt = 1/fs; 
T = (0:size(Mat,1)-1)' * dt;
% T = time(22001:34001);   % specific plotting window
%% A-B plot example waveforms

% ####### A:

basePath = '/chemsyn/synchro/passive/EFP';   % change to '../active/EFP' for active models
% ===================== PROXIMAL (50–100 µm) =====================
cd(fullfile(basePath, 'electrode_numbers_50-100'));

% Load electrode indices
q = []; % if already loaded and saved call the mat file: load('electrode_numbers.mat'); q = electrode_numbers (n=15)
% combine proximal electrode data for traces
data1 = [];
for k = 1:numel(q)
    data1(:,k) = lfpMat1(:, q(k)); % arrays whose electrode numbers are proximal
end
data2 = [];
for k = 1:numel(q)
    data2(:,k) = lfpMat2(:, q(k)); % arrays whose electrode numbers are proximal
end
LFP_prox = [data1, data2]; % Combine
avg_prox = mean(LFP_prox, 2); % Average across proximal electrodes

figure; % Plot
hold on;
for j = 1:15
    plot(T, LFP_prox(:,j), 'b', 'LineWidth', 1.5); % 'r' for active
end

plot(T, avg_prox, 'k', 'LineWidth', 3); 
title('LFPs: 50–100 µm (proximal- passive)');
xlabel('Time');
ylabel('LFP (uV)');
hold off;

% ===================== distal (190-300 µm) =====================
cd(fullfile(basePath, 'electrode_numbers_190-300'));
% Load electrode indices
q = []; %  load('electrode_numbers.mat'); q = electrode_numbers (n=21)
% combine distal electrode data for traces
data1 = [];
for k = 1:numel(q)
    data1(:,k) = lfpMat6(:, q(k)); % arrays whose electrode numbers are picked as distal
end
data2 = [];
for k = 1:numel(q)
    data2(:,k) = lfpMat7(:, q(k)); % arrays whose electrode numbers are picked as distal
end
LFP_dist = [data1, data2]; % Combine
avg_dist = mean(LFP_dist, 2); % Average across distal electrodes

figure; % Plot
hold on;
for j = 1:21
    plot(T, LFP_dist(:,j), 'b', 'LineWidth', 1.5); % 'r' for active
end
plot(T, avg_dist, 'k', 'LineWidth', 3);
title('LFPs: 190-300 µm (distal- passive)'); % change to active
xlabel('Time (ms)');
ylabel('LFP (uV)');
hold off;

% ####### B: 
% SAME AS 'A' but for gap junctions: /gapj/synchro/passive/EFP, and /gapj/synchro/active/EFP

%% Boxplot- active vs passive

nBins = 10;
activePath  = '/sync/active/Negativeamps'; % files have negative amplitudes (negative peaks-baseline) from each electrodes divided into 10 bins based on radial distance from soma
passivePath = '/sync/passive/Negativeamps';
colors = [
    190/255, 30/255, 45/255;   % active
    27/255, 117/255, 187/255   % passive
];
activeColor  = colors(1,:);
passiveColor = colors(2,:);
activeData  = cell(nBins,1);
passiveData = cell(nBins,1);

for i = 1:nBins
    fname = sprintf('Bin%d.txt', i);
    activeData{i}  = load(fullfile(activePath, fname));
    passiveData{i} = load(fullfile(passivePath, fname));
end

meanActive  = cellfun(@mean, activeData);
meanPassive = cellfun(@mean, passiveData);

x = 1:nBins;
xActive  = x - 0.2;
xPassive = x + 0.2;
figure; hold on;
boxWidth = 0.35;
% boxplots: 
for i = 1:nBins

    % ---- ACTIVE ----
    boxplot(activeData{i}, 'Positions', xActive(i), ...
        'Widths', boxWidth, ...
        'Colors', activeColor, ...
        'Symbol', '');
    h = findobj(gca,'Tag','Box');
    patch(get(h(1),'XData'), get(h(1),'YData'), ...
        activeColor, 'FaceAlpha', 0.2, ...
        'EdgeColor', activeColor, 'LineWidth', 2);

    % ---- PASSIVE ----
    boxplot(passiveData{i}, 'Positions', xPassive(i), ...
        'Widths', boxWidth, ...
        'Colors', passiveColor, ...
        'Symbol', '');
    h = findobj(gca,'Tag','Box');
    patch(get(h(1),'XData'), get(h(1),'YData'), ...
        passiveColor, 'FaceAlpha', 0.2, ...
        'EdgeColor', passiveColor, 'LineWidth', 2);

end
% set(findobj(gca,'Type','line'),'LineStyle','-');
xticks(linspace(0,250,6));
xticklabels(string(linspace(0,250,6)));
xlabel('Radial distance (µm)');
ylabel('LFP');
title('Active vs Passive LFP');
set(gca, 'YScale', 'log');
box off;
hold off;

% Inset: line plot:  
fogure; hold on;
plot(x, meanActive, '-', ...
    'Color', activeColor, ...
    'LineWidth', 2);
plot(x, meanPassive, '-', ...
    'Color', passiveColor, ...
    'LineWidth', 2);
set(gca, 'XLim', [0 250]);
box off;
hold off;



%% waveforms spanning the entire neuron (E)

% basal waveforms from above mentioned filtered mat, load and filter EFP recorded from apical dendritic side
dataFiles_apical = {
    'S1_allLFP_E0-E48.xlsx'
    'S2_allLFP_E0-E48.xlsx'
    'S3_allLFP_E0-E48.xlsx'
    'S4_allLFP_E0-E48.xlsx'
    'S5_allLFP_E0-E48.xlsx'
    'S6_allLFP_E0-E48.xlsx'
    'S7_allLFP_E0-E48.xlsx'
};
dataRange = 'A1:AW80000'; 
nFiles = numel(dataFiles); % seven

lfpData_apical = [];
for f = 1:nFiles
    Mat_apical = xlsread(dataFiles_apical{f}, dataRange); % Load dataset
    lfpMat = zeros(size(Mat_apical));     % filtered matrix
    % Filter each electrode trace
    for i = 1:size(Mat_apical,2)
        lfpMat(:,i) = filtfilt(b, a, Mat_apical(:,i)); %forward-backward filtering for no time shift
    end
    lfpData_apical{f} = lfpMat_apical;
end
filteredmat_apical = cat(2, lfpData_apical{:}); % concatenating all arrays column-wise (7 arrays each with 49 electrodes) in one matrix
T = T(4001:48001,1);

% basal:
waveform_1 = FilteredMat2(4001:48001,32); waveforms(:,1) = waveform_1;
waveform_2 = FilteredMat3(4001:48001,32); waveforms(:,2) = waveform_2;
waveform_3 = FilteredMat3(4001:48001,39); waveforms(:,3) = waveform_3;
waveform_4 = FilteredMat1(4001:48001,34); waveforms(:,4) = waveform_4;
waveform_5 = FilteredMat3(4001:48001,45); waveforms(:,5) = waveform_5;
waveform_6 = FilteredMat1(4001:48001,36); waveforms(:,6) = waveform_6;
waveform_7 = FilteredMat6(4001:48001,47); waveforms(:,7) = waveform_7;
% waveform_8 = FilteredMat2(4001:48001,4); waveforms(:,8) = waveform_8;
% waveform_9 = FilteredMat3(4001:48001,1); waveforms(:,9) = waveform_9;
% waveform_10 = FilteredMat6(4001:48001,7); waveforms(:,10) = waveform_10;

% apical
Ap_waveforms_A = zeros(length(T), num_Ap_waveforms);
Ap_waveform_11 = filteredmat_apical2(4001:48001,32); Ap_waveforms_A(:,1) = Ap_waveform_11;
Ap_waveform_12 = filteredmat_apical3(4001:48001,32); Ap_waveforms_A(:,2) = Ap_waveform_12;
Ap_waveform_13 = filteredmat_apical3(4001:48001,39); Ap_waveforms_A(:,3) = Ap_waveform_13;
Ap_waveform_14 = filteredmat_apical1(4001:48001,34); Ap_waveforms_A(:,4) = Ap_waveform_14;
Ap_waveform_15 = filteredmat_apical3(4001:48001,45); Ap_waveforms_A(:,5) = Ap_waveform_15;
Ap_waveform_16 = filteredmat_apical1(4001:48001,36); Ap_waveforms_A(:,6) = Ap_waveform_16;
Ap_waveform_17 = filteredmat_apical6(4001:48001,47); Ap_waveforms_A(:,7) = Ap_waveform_17;
% Ap_waveform_18 = filteredmat_apical2(4001:48001,4); Ap_waveforms_A(:,8) = Ap_waveform_18;
% Ap_waveform_19 = filteredmat_apical3(4001:48001,1); Ap_waveforms_A(:,9) = Ap_waveform_19;
% Ap_waveform_20 = filteredmat_apical6(4001:48001,7); Ap_waveforms_A(:,10) = Ap_waveform_20;


figure;
num_waveforms = 7; % Number of waveforms in each set (taking 7 for plotting)
total_waveforms = 14; % Total number
for i = 1:total_waveforms
    subplot(total_waveforms, 1, i);
    if i <= num_waveforms
        % Plot the first 7 waveforms from waveforms (basal)
        plot(T, waveforms(:,i));
    else
        % Plot the next 7 waveforms from Ap_waveforms (apical)
        plot(T, Ap_waveforms_(:,i - num_waveforms));
    end
    xlabel('Time');
    ylabel('Amplitude');
    % title(['Waveform ', num2str(i)]); % Uncomment this line if you want to add titles
    % xlim([100 1200]);
    % ylim([-100 100]);
    box off;
end

% do the same for gap junctional EFPs for right side waveforms.


