% GJLFP, eLife 2026. Richa Sirmaur and Rishikesh Narayanan.
% Figure 7: Patterned inputs: A)- B) Distance-dependence of Max spectral power for all frequencies, fourier of all electrodes for all frequencies - active vs passive

%%

fs = 40000;
freqs = [1 2 4 8 16 32 64 128];
fc1_all = [0.5 1.5 3 6 10 30 60 120]; % cutoff setting for 1, 2, 4, 8, 16, 32, 64, 128 Hz respectively
fc2_all = [1.5 3 5 10 20 40 70 132];  % cutoff setting for 1, 2, 4, 8, 16, 32, 64, 128 Hz respectively

nBins = 10;
distances = load('electrode_distances.txt');

edges = linspace(0,250,nBins+1);
binCenters = (edges(1:end-1)+edges(2:end))/2;
%% A) active 
mean_active = zeros(length(freqs), nBins);
sem_active  = zeros(length(freqs), nBins);

freqColors = [ % change to exact colors in illustrator
    1    190/255 30/255 45/255;      % red
    2    0       0       0.5;         % dark blue
    4    1       0       1;           % magenta
    8    0.3010  0.7450  0.9330;      % sky blue
    16   0.4660  0.6740  0.1880;      % green
    32   0.9290  0.6940  0.1250;      % dark yellow
    64   0.6350  0.0780  0.1840;      % maroon
    128  0.4940  0.1840  0.5560       % purple
];
L = size(MatA,1);
f = fs*(0:(L/2))/L;


for fidx = 1:length(freqs)
    f0  = freqs(fidx);
    fc1 = fc1_all(fidx);
    fc2 = fc2_all(fidx);

    [b,a] = butter(1,[fc1 fc2]/(fs/2),'bandpass');
    % load active efp data
    data = load(fullfile('active', sprintf('%dHz/allMATS.mat', f0))); % for specific freq folders
    Mat = cat(2, data.allMATS{:});  % combine 7 mats → 343 electrodes
    Mat(1:2000,:) = repmat(Mat(2001,:),2000,1); % remove transient
    peakPower = zeros(343,1);
    % L = size(Mat,1);
    % P_all = zeros(L/2+1, nElectrodes);

    for e = 1:343
        sig = filtfilt(b,a,Mat(:,e));
        P = abs(fft(sig)/length(sig));
        P = P(1:end/2+1);
        P(2:end-1) = 2*P(2:end-1);
        % P_all(:,e) = P;
        peakPower(e) = max(P);
    end
    % binning according to radial distance:
    for bidx = 1:nBins
        idx = distances >= edges(bidx) & distances < edges(bidx+1);
        vals = peakPower(idx);
        mean_active(fidx,bidx) = mean(vals);
        sem_active(fidx,bidx)  = std(vals)/sqrt(numel(vals));
    end
end

% figure row 1
figure; hold on;
for fidx = 1:length(freqs)
    c = freqColors(fidx,2:4);
    errorbar(binCenters, mean_active(fidx,:), sem_active(fidx,:), ...
        '-o', 'Color', c, ...
        'LineWidth', 2, ...
        'MarkerFaceColor', c);
end

xlabel('Radial distance (\mum)');
ylabel('Max LFP Power (a.u)');
title('ACTIVE');
xlim([0 250]);
box off;

% figure row 2: spectrum all frequencies:
for fidx = 1:length(freqs)
    for e = 1:343
        sig = filtfilt(b,a,Mat(:,e));
        Y = fft(sig)/L;
        P = abs(Y);
        P = P(1:L/2+1);
        P(2:end-1) = 2*P(2:end-1);
        figure; hold on; plot(f, P, 'Color', [c 0.08]); %  plot each electrode spectrum
        
    end
end
set(gca,'XScale','log');
xlabel('Frequency (Hz)');
ylabel('LFP magnitude (a.u)');
title('All electrodes Fourier spectra- all frequencies');
box off;

% figure row 3: total transmembrane current-
figure; hold on;
L = []; 
for fidx = 1:length(freqs)
    f0  = freqs(fidx);
    fc1 = fc1_all(fidx);
    fc2 = fc2_all(fidx);
    [b,a] = butter(1,[fc1 fc2]/(fs/2),'bandpass');
    data = load(fullfile('active', sprintf('%dHz/all_TM_mats.mat', f0))); % load total transmembrane currents- all intracellular sites - all frequencies
    Mat = all_TM_mats;
    L = size(Mat,1);
    f = fs*(0:(L/2))/L;
    c = freqColors(fidx,:);
    nE = size(Mat,2);
    %fft
    for e = 1:nE
        sig = filtfilt(b,a,Mat(:,e));
        Y = fft(sig)/L;
        P = abs(Y);
        P = P(1:L/2+1);
        P(2:end-1) = 2*P(2:end-1);
        plot(f, P, 'Color', [c 0.06]);
    end
end
set(gca,'XScale','log');
xlabel('Frequency (Hz)');
ylabel('Current magnitude (a.u.)');
title('FFT of TTM (all intra sites × all frequencies)');
box off;


%% B) passive:
mean_passive = zeros(length(freqs), nBins);
sem_passive  = zeros(length(freqs), nBins);

for fidx = 1:length(freqs)

    f0  = freqs(fidx);
    fc1 = fc1_all(fidx);
    fc2 = fc2_all(fidx);
    [b,a] = butter(1,[fc1 fc2]/(fs/2),'bandpass');

    data = load(fullfile('passive', sprintf('%dHz/allMATS.mat', f0)));
    Mat = cat(2, data.allMATS{:});
    Mat(1:2000,:) = repmat(Mat(2001,:),2000,1);
    peakPower = zeros(343,1);

    for e = 1:343
        sig = filtfilt(b,a,Mat(:,e));
        P = abs(fft(sig)/length(sig));
        P = P(1:end/2+1);
        P(2:end-1) = 2*P(2:end-1);
        peakPower(e) = max(P);
    end
    for bidx = 1:nBins
        idx = distances >= edges(bidx) & distances < edges(bidx+1);
        vals = peakPower(idx);
        mean_passive(fidx,bidx) = mean(vals);
        sem_passive(fidx,bidx)  = std(vals)/sqrt(numel(vals));
    end
end

% figure row 1:Max LFP Power
figure; hold on;
for fidx = 1:length(freqs)
    c = freqColors(fidx,2:4);
    errorbar(binCenters, mean_passive(fidx,:), sem_passive(fidx,:), ...
        '-o', 'Color', c, ...
        'LineWidth', 2, ...
        'MarkerFaceColor', c);
end
xlabel('Radial distance (\mum)');
ylabel('Max LFP Power (a.u)');
title('PASSIVE');
xlim([0 250]);
box off;

% % figure row 2: spectrum all frequencies:
for fidx = 1:length(freqs)
    for e = 1:343
        sig = filtfilt(b,a,Mat(:,e));
        Y = fft(sig)/L;
        P = abs(Y);
        P = P(1:L/2+1);
        P(2:end-1) = 2*P(2:end-1);
        figure; hold on; plot(f, P, 'Color', [c 0.08]); %  plot each electrode spectrum
        
    end
end
set(gca,'XScale','log');
xlabel('Frequency (Hz)');
ylabel('LFP magnitude (a.u)');
title('All electrodes Fourier spectra- all frequencies');
box off;

% figure row 3: total transmembrane current-
figure; hold on;
L = []; 
for fidx = 1:length(freqs)
    f0  = freqs(fidx);
    fc1 = fc1_all(fidx);
    fc2 = fc2_all(fidx);
    [b,a] = butter(1,[fc1 fc2]/(fs/2),'bandpass');
    data = load(fullfile('passive', sprintf('%dHz/all_TM_mats.mat', f0))); % load total transmembrane currents- all intracellular sites - all frequencies
    Mat = all_TM_mats;
    L = size(Mat,1);
    f = fs*(0:(L/2))/L;
    c = freqColors(fidx,:);
    nE = size(Mat,2);
    %fft
    for e = 1:nE
        sig = filtfilt(b,a,Mat(:,e));
        Y = fft(sig)/L;
        P = abs(Y);
        P = P(1:L/2+1);
        P(2:end-1) = 2*P(2:end-1);
        plot(f, P, 'Color', [c 0.06]);
    end
end
set(gca,'XScale','log');
xlabel('Frequency (Hz)');
ylabel('Current magnitude (a.u.)');
title('FFT of TTM (all intra sites × all frequencies)');
box off;

%% C) and D) same as A) and B) respectively, BUT FOR GAP JUNCTIONS. EFP loaded here will be from simulations of gap junctional inputs.

