% GJLFP, eLife 2026. Richa Sirmaur and Rishikesh Narayanan.
% Figure 2: waveforms (A-B, F-G),lfp deflection plots (C-E, H-J)

%% waveforms - all models
% this is chemical synapses, do the same for data from gap junctions

fs  = 40000;
fc1 = 0.5;
fc2 = 300;

Wn = [fc1 fc2] / (fs/2);
n  = 2;
[b,a] = butter(n, Wn, 'bandpass');

dataRange = 'A1:AW80000';
models = {'default','nona','noleak','nonaleak'};
for m = 1:numel(models)
    modelName = models{m};
    basePath = fullfile('/folderdata', modelName);  % set path 
    dataFiles = {
        fullfile(basePath,'S1_allLFP_E0-E48.xlsx')
        fullfile(basePath,'S2_allLFP_E0-E48.xlsx')
        fullfile(basePath,'S3_allLFP_E0-E48.xlsx')
        fullfile(basePath,'S4_allLFP_E0-E48.xlsx')
        fullfile(basePath,'S5_allLFP_E0-E48.xlsx')
        fullfile(basePath,'S6_allLFP_E0-E48.xlsx')
        fullfile(basePath,'S7_allLFP_E0-E48.xlsx')
    };

    nFiles = numel(dataFiles);

    % filtering 
    for f = 1:nFiles 
        Mat = xlsread(dataFiles{f}, dataRange);
        FilteredMat = zeros(size(Mat));
        for i = 1:size(Mat,2)
            FilteredMat(:,i) = filtfilt(b,a,Mat(:,i));
        end
     
    end

end


baseRoots = {
    '/chemsyn/synchro/active/EFP_default' 
    '/chemsyn/synchro/active/EFP_nona'
    '/chemsyn/synchro/active/EFP_noleak'
    '/chemsyn/synchro/active/EFP_nonaleak'
};
dt = 1/fs; 
T = (0:size(Mat,1)-1)' * dt;
colors = lines(numel(models));
for m = 1:numel(models)

    basePath = baseRoots{m};

    % ===================== proximal ==========================
    cd(fullfile(basePath, 'electrode_numbers_50-100'));

    % load('electrode_numbers.mat');   % gives q (n=15)
    % q = electrode_numbers; OR
    q=15
 
    data1 = [];
    for k = 1:numel(q)
        data1(:,k) = lfpMat1(:, q(k)); % arrays whose electrode numbers are proximal: lfpMat1= FilteredMat(:,1)
    end

    data2 = [];
    for k = 1:numel(q)
        data2(:,k) = lfpMat2(:, q(k)); % arrays whose electrode numbers are proximal lfpMat2= FilteredMat(:,2)
    end

    LFP_prox = [data1, data2];
    avg_prox = mean(LFP_prox, 2);

    figure(100+m); hold on;
    for j = 1:size(LFP_prox,2)
        plot(T, LFP_prox(:,j),'b', 'LineWidth', 1.5);
    end

    plot(T, avg_prox, 'k', 'LineWidth', 3);

    title(['LFPs: 50–100 µm (proximal) - ' models{m}]);
    xlabel('Time(ms)');
    ylabel('LFP (\muV)');
    box off;
    hold off;

    % ===================== distal ============================
    cd(fullfile(basePath, 'electrode_numbers_190-300'));
    % load('electrode_numbers.mat');   % gives q (n=21)
    % q = electrode_numbers; OR
    q=21

    data1 = [];
    for k = 1:numel(q)
        data1(:,k) = lfpMat6(:, q(k)); %arrays whose electrode numbers are distal lfpMat6= FilteredMat(:,6)
    end

    data2 = [];
    for k = 1:numel(q)
        data2(:,k) = lfpMat7(:, q(k)); %arrays whose electrode numbers are distal lfpMat7= FilteredMat(:,7)
    end

    LFP_dist = [data1, data2];
    avg_dist = mean(LFP_dist, 2);

    figure(200+m); hold on;
    for j = 1:size(LFP_dist,2)
        plot(T, LFP_dist(:,j), 'b', 'LineWidth', 1.5);
    end

    plot(T, avg_dist, 'k', 'LineWidth', 3);

    title(['LFPs: 190–300 µm (distal) - ' models{m}]); % change the color for different models in illustrator
    xlabel('Time(ms)');
    ylabel('LFP (\muV)');
    box off;
    hold off;

end

%% lfp deflection plots (C-E, H-J)
% this is chemical synapses, do the same for data from gap junctions

fs  = 40000;
fc1 = 0.5;
fc2 = 300;

Wn = [fc1 fc2] / (fs/2);
n  = 2;
[b,a] = butter(n, Wn, 'bandpass');

dataRange = 'A1:AW80000';
baseline_idx = 30001:33401;

models = {'default','nona','noleak','nonaleak'};
% colors = lines(numel(models));

figure; hold on;
for m = 1:numel(models)

    modelName = models{m};
    basePath = fullfile('/folderdata', modelName);  % set path 
    dataFiles = {
        fullfile(basePath,'S1_allLFP_E0-E48.xlsx')
        fullfile(basePath,'S2_allLFP_E0-E48.xlsx')
        fullfile(basePath,'S3_allLFP_E0-E48.xlsx')
        fullfile(basePath,'S4_allLFP_E0-E48.xlsx')
        fullfile(basePath,'S5_allLFP_E0-E48.xlsx')
        fullfile(basePath,'S6_allLFP_E0-E48.xlsx')
        fullfile(basePath,'S7_allLFP_E0-E48.xlsx')
    };

    nFiles = numel(dataFiles);

    minAmp_all = [];
    maxAmp_all = [];
    ttpAmp_all = [];

    % filtering 
    for f = 1:nFiles

        Mat = xlsread(dataFiles{f}, dataRange);
        FilteredMat = zeros(size(Mat));
        for i = 1:size(Mat,2)
            FilteredMat(:,i) = filtfilt(b,a,Mat(:,i));
        end

        minVals = min(FilteredMat,[],1);
        maxVals = max(FilteredMat,[],1);

        baseline = mean(FilteredMat(baseline_idx,:),1); % indices before the start of simulation

        minAmp = (minVals - baseline)';
        maxAmp = (maxVals - baseline)';
        ttpAmp = (maxVals - minVals)';

        minAmp_all = [minAmp_all; minAmp];
        maxAmp_all = [maxAmp_all; maxAmp];
        ttpAmp_all = [ttpAmp_all; ttpAmp];

    end

    % save
    writematrix(minAmp_all, ['LFP_minAmp_' modelName '.txt']);
    writematrix(maxAmp_all, ['LFP_maxAmp_' modelName '.txt']);
    writematrix(ttpAmp_all, ['LFP_TTPAmp_' modelName '.txt']);

end

nBins = 10;
maxDist = 250;
% LOAD distances = (343 x 1 vector) %radial distance of each electrodes from soma- saved
distances = distances(:);
edges = linspace(0, maxDist, nBins+1);

%  assign electrodes to bin
binID = zeros(size(distances));
for b = 1:nBins
    idx = distances >= edges(b) & distances < edges(b+1);
    binID(idx) = b;
end
binID(distances == maxDist) = nBins; % include max value
binCenters = (edges(1:end-1) + edges(2:end)) / 2; % if want
save('electrode_bins.mat', 'binID', 'binCenters', 'edges');

% load('electrode_bins.mat');  % if reloading
models = {'default','nona','noleak','nonaleak'}; % data saved in different folders with these names
colors = lines(numel(models));
nBins = numel(binCenters);

figure(1); hold on; title('Negative amplitude vs distance');
figure(2); hold on; title('Positive amplitude vs distance');
figure(3); hold on; title('TTP amplitude vs distance');

for m = 1:numel(models)

    modelName = models{m};

    % load
    minAmp = load(['LFP_minAmp_' modelName '.txt']);
    maxAmp = load(['LFP_maxAmp_' modelName '.txt']);
    ttpAmp = load(['LFP_TTPAmp_' modelName '.txt']);

    minAmp = minAmp(:);
    maxAmp = maxAmp(:);
    ttpAmp = ttpAmp(:);

    meanMin = zeros(nBins,1);
    meanMax = zeros(nBins,1);
    meanTTP = zeros(nBins,1);

    for b = 1:nBins
        idx = (binID == b);
        meanMin(b) = mean(minAmp(idx));
        meanMax(b) = mean(maxAmp(idx));
        meanTTP(b) = mean(ttpAmp(idx));
    end

    figure(1);
    plot(binCenters, meanMin, '-o', ...
        'Color', colors(m,:), ...       % colors will be changed in illustrator..
        'LineWidth', 2, ...
        'MarkerFaceColor', colors(m,:));

    figure(2);
    plot(binCenters, meanMax, '-o', ...
        'Color', colors(m,:), ...
        'LineWidth', 2, ...
        'MarkerFaceColor', colors(m,:));

    figure(3);
    plot(binCenters, meanTTP, '-o', ...
        'Color', colors(m,:), ... 
        'LineWidth', 2, ...
        'MarkerFaceColor', colors(m,:));

end

figure(1);
xlabel('Radial distance (\mum)');
ylabel('Negative deflection(\muV)');
legend(models); box off;

figure(2);
xlabel('Radial distance (\mum)');
ylabel('Positive deflection(\muV)');
legend(models); box off;

figure(3);
xlabel('Radial distance (\mum)');
ylabel('Total amplitude(\muV)');
legend(models); box off;