% GJLFP, eLife 2026. Richa Sirmaur and Rishikesh Narayanan.
% Figure 8: Patterned inputs: spike-lfp phase differences: chem syn vs gap junctions

%% A) The following is for chemical synapses:
fs = 40000;
time = 0:0.025:2000;
T = time(1:80000).';
% 16 Hz band
fc1 = 10;
fc2 = 20;
[b,a] = butter(1, [fc1 fc2]/(fs/2), 'bandpass');
freq = 16;
period = 1000/freq;     % ms
win = 5 * period;       % 5 cycles
tStart = 500;
tEnd = tStart + win;
idxT = T >= tStart & T <= tEnd;
% 1)
cd('default_model_folder')
Mat_EFP=load('EFP.mat')   ;   
somav = load('Somavolt.txt');

efp = filtfilt(b,a,Mat_EFP(:,32)); % Mat2 electrode 32is 14.11 µm from soma
somav_f = filtfilt(b,a,somav);

figure; hold on;

plot(T(idxT), efp(idxT), 'm', 'LineWidth', 2);      % EFP (magenta)
plot(T(idxT), somav(idxT), 'k', 'LineWidth', 1.5); % soma (black)

xlabel('Time (ms)');
ylabel('Amplitude');
title('DEFAULT MODEL - 16 Hz (5 cycles)');
box off;

%2)
% cd('default_model_folder')
% filtering both
efp = filtfilt(b,a,efp_raw);
intra = filtfilt(b,a,somav);

figure; hold on;
plot(T(idxT), efp(idxT), 'm', 'LineWidth', 2);
plot(T(idxT), intra(idxT), 'k', 'LineWidth', 2);

xlabel('Time (ms)');
ylabel('Amplitude');
title('DEFAULT MODEL - BOTH FILTERED (16 Hz, 5 cycles)');
box off;

%3)
cd('no_sodium_folder')

Mat_EFP=load('EFP.mat')
somav = load('Somavolt.txt');
efp = filtfilt(b,a,Mat_EFP(:,32));

figure; hold on;
plot(T(idxT), efp(idxT), 'm', 'LineWidth', 2);   
plot(T(idxT), somav(idxT), 'k', 'LineWidth', 2); 
xlabel('Time (ms)');
ylabel('Amplitude');
title('NO SODIUM MODEL - 16 Hz (5 cycles)');
box off;

% do same as A) for gap junctions- load data from GJ folder.

%% B) Spike-LFP phase plot:
% clear; clc; close all;
fs = 40000;

freqs   = [1 2 4 8 16 32 64 128];
fc1_all = [0.5 1.5 3 6 10 30 60 120];
fc2_all = [1.5 3 5 10 20 40 70 132];

phi_all = cell(1, length(freqs));

load('MATS2.mat');              % Mat2
somav_full = load('Somavolt.txt');

for fidx = 1:length(freqs)
    f0  = freqs(fidx);
    fc1 = fc1_all(fidx);
    fc2 = fc2_all(fidx);
    fprintf('Processing %d Hz...\n', f0);
    switch f0
        case 1
            somav = somav_full(1:320001);
        case 2
            somav = somav_full(1:240000);
        case 4
            somav = somav_full(1:180000);
        case 8
            somav = somav_full(1:160000);
        case {16,32}
            somav = somav_full(1:100000);
        case {64,128}
            somav = somav_full(1:60000);
    end

    N = length(somav);
    t = (0:N-1)'/fs;   
    Wn = [fc1 fc2] / (fs/2);
    [b,a] = butter(2, Wn, 'bandpass');
    FilteredMat2 = zeros(size(Mat2,1), size(Mat2,2));
    for ch = 1:size(Mat2,2)
        FilteredMat2(:,ch) = filtfilt(b,a,Mat2(:,ch));
    end

    LFP = FilteredMat2(:,32); % electrode closer to soma (14um)

    switch f0 % (1-4hz)15 cycles, (8-128hz)25 cycles : refer to suppli table 4
        case 1
            idx = 15001:598001;
        case 2
            idx = 20281:315601;
        case 4
            idx = 24121:171401;
        case 8
            idx = 26601:150801;
        case 16
            idx = 27401:89801;
        case 32
            idx = 28041:59001;
        case 64
            idx = 28561:43481;
        case 128
            idx = 28321:35761;
    end

    LFP = LFP(idx);
    somav = somav(idx);
    t = t(idx);
    % peak detection
    [~, peakVlocs] = findpeaks(somav, ...
        'MinPeakProminence', std(somav));
    peakVtimings = t(peakVlocs);

    [~, trough_locs] = findpeaks(-LFP); %LFP troughs
    troughtimings = t(trough_locs);

    %  first spike per cycle
    cycle_ranges = [troughtimings(1:end-1), troughtimings(2:end)];
    peakVtimings2 = peakVtimings;
    
    peakVtimings2 = peakVtimings2( ... %% keep only spikes within valid range
        peakVtimings2 >= min(troughtimings) & ...
        peakVtimings2 <= max(troughtimings));
    selected_timings_within_cycles = nan(size(cycle_ranges,1),1);

    for i = 1:size(cycle_ranges,1)
        cycle_start = cycle_ranges(i,1);
        cycle_end   = cycle_ranges(i,2);
        spikes_in_cycle = peakVtimings2( ...
            peakVtimings2 >= cycle_start & ...
            peakVtimings2 <= cycle_end);
        if ~isempty(spikes_in_cycle)
            selected_timings_within_cycles(i) = spikes_in_cycle(1); % first spike
        end
    end
    selected_timings_within_cycles = ...
        selected_timings_within_cycles(~isnan(selected_timings_within_cycles));
    % phase
    phi = [];
    for i = 1:length(selected_timings_within_cycles)
        Tspike = selected_timings_within_cycles(i);
        idx_left = find(troughtimings <= Tspike); %Find the indices of the points in vector awhere Tspike lies between
        if isempty(idx_left) || idx_left(end) == length(troughtimings)
            continue;
        end
        Tleft  = troughtimings(idx_left(end));
        Tright = troughtimings(idx_left(end)+1);
        phase = 360 * ((Tspike - Tleft) / (Tright - Tleft));
        phi(end+1) = phase;
    end

    phi = mod(phi, 360);
    phi_all{fidx} = phi(:);
end


figure; hold on;
lightMaroon = [0.7 0.2 0.2]; % change colors in illustrator
darkMaroon  = [0.4 0 0];
medians = nan(1,length(freqs));
for fidx = 1:length(freqs)
    y = phi_all{fidx};
    if isempty(y), continue; end
    jitter = (rand(size(y)) - 0.5) * 0.3;
    xvals = fidx + jitter;
    % spikelfp phase
    scatter(xvals, y, 12, lightMaroon, ...
        'filled', 'MarkerFaceAlpha', 0.4);
    % overlaying median
    medians(fidx) = median(y);
    scatter(fidx, medians(fidx), 80, darkMaroon, 'o', 'filled');
end

% connect medians
plot(1:length(freqs), medians, '-', ...
    'Color', darkMaroon, 'LineWidth', 2);

xticks(1:length(freqs));
xticklabels({'1','2','4','8','16','32','64','128'});
xlabel('Frequency (Hz)');
ylabel('Spike - LFP Phase (degrees)');
ylim([0 360]);
yticks(0:90:360);
title('Spike-LFP Phase');
box off;
set(gca, 'TickDir', 'out', 'LineWidth', 1);

%% Do the same as B) for gap junctions, and do hold on for plots (overlaying/ using shifted bin approach)..

%% C) Polar plots: PLOTS ARE MADE IN IGOR PRO, but medians were saved using this:


freqs = [1 2 4 8 16 32 64 128];
nTrials = 5;
all_medians = zeros(nTrials, length(freqs));
basePath = 'chemsyn/alltrials';  % same for gap junctions
for tr = 1:nTrials
    fprintf('Processing Trial %d...\n', tr);
    trialPath = fullfile(basePath, sprintf('Trial%d', tr));

    for fidx = 1:length(freqs)
        f0 = freqs(fidx);
        freqFolder = fullfile(trialPath, sprintf('%dHz', f0));
        dataFile = fullfile(freqFolder, 'phi.mat'); % all Spike-lfp phases (calculated as above are saved)
        if ~exist(dataFile, 'file')
            warning('Missing: %s', dataFile);
            all_medians(tr, fidx) = NaN;
            continue;
        end
        S = load(dataFile);   %  phi inside
        phi = S.phi;
        if isempty(phi)
            all_medians(tr, fidx) = NaN;
        else
            all_medians(tr, fidx) = median(mod(phi,360));
        end

    end
end

grand_median = nanmedian(all_medians, 1);
disp('Final 8 frequency medians:');
disp(grand_median);


