globalAnalysis2.m

%% Generates outputs for non-EVC visually responsive electrodes, after running globalAnalysis1 and data on the cluster
%
%    If this code is used in a publication, please cite the manuscript:
%    "H Huang, KN Kay, NM Gregg, G Ojeda Valencia, M In, C Kapeller, Y Shu, GA Worrell, KJ Miller, and D Hermes.
%    Single pulse electrical stimulation in white matter modulates iEEG visual responses in human early visual cortex. (Under Review)"
%
%    A preprint is available currently at doi: https://doi.org/10.1101/2025.05.05.652264.
%
%    The dataset corresponding to this code and manuscript is in BIDS format (version 1.10.0) on OpenNeuro (ds006485),
%    and it will be made publicly available upon manuscript acceptance.
%
%    Copyright (C) 2025 Harvey Huang
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <https://www.gnu.org/licenses/>.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% Configure paths

saveOutputs = true;

%sub = '1'; sites = {'LOC4-LOC5', 'LG6-LG7'};
sub = '2'; sites = {'ROC6-ROC7', 'RMO3-RMO4'};

outdir = fullfile('output', sub, 'global');
mkdir(outdir);

%% 1) Load the ERP all model fits (copied form analysisERPFIR2) and make bargraphs

loaddir = fullfile('output', 'mforge'); % to load results from mforge. Use outdir for local
ylimsRsq = [-60, 50];

rsqTable = readtable(fullfile(outdir, sprintf('BBRsq_%s.tsv', sub)), 'FileType', 'text', 'Delimiter', '\t');

rsqTableSig = rsqTable(rsqTable.responsive == 1, :);
rsqTableSig = sortrows(rsqTableSig, 'rsq', 'descend'); % sort by highest to lowest rsq
fprintf('%d significant channels found\n', height(rsqTableSig));

for cc = 1:height(rsqTableSig)
    chName = rsqTableSig.name{cc};
    
    for ss = 1:2
        site = sites{ss};

        fprintf('Loading odd-even ERP fits for %s %s %s\n', sub, chName, site);

        mdl = load(fullfile(loaddir, sprintf('firOddEvenAllModelsErp_%s_%s_%s.mat', sub, chName, site)));
        cods = mdl.sigExpsTest; % use unnormalized R^2. Rationale: zero is meaningful and since we're taking different windows we don't want to drastically shift baseline
        sampSz = size(cods, 1); % sample size of the testing r-squareds
        twinCod = mdl.config.twincod; % get twinCod from last one loaded
        sems = std(cods)/sqrt(sampSz); % SE of mean

        % bargraph with SD error bars
        figure('Position', [200, 200, 400, 600]); hold on
        bar(mean(cods), 'FaceColor', [0.5, 0.5, 0.5]);
        errorbar(mean(cods), sems, 'LineWidth', 1, 'LineStyle', 'none', 'CapSize', 10, 'Color', 'k'); % 95% conf intervals
        set(gca, 'xtick', 1:5, 'xticklabels', {'Simple', 'Img', 'ISI', 'Full', 'Mean'});
        xlabel('Model'); ylim(ylimsRsq); ylabel('COD');
        if saveOutputs
            saveas(gcf, fullfile(outdir, sprintf('%s_%s_compare_oddeven_cod', chName, site)), 'png');
            saveas(gcf, fullfile(outdir, sprintf('%s_%s_compare_oddeven_cod', chName, site)), 'svg');
        end
        
        fprintf('Odd-even COD by model, mean+/-SEM [--, %0.2fs]: [Simple]%0.01f%%+/-%0.01f%%, [Coh]%0.01f%%+/-%0.01f%%, [ISI]%0.01f%%+/-%0.01f%%, [Full]%0.01f%%+/-%0.01f%%\n', ...
            twinCod(2), mean(cods(:, 1)), sems(1), mean(cods(:, 2)), sems(2), mean(cods(:, 3)), sems(3), mean(cods(:, 4)), sems(4));
        
        % Stats: check if img, isi models are better than simple. Then check if full model is better than others
        [~, pImg] = ttest(cods(:, 2) - cods(:, 1), [], 'Tail', 'Right'); % img better than simple
        [~, pIsi] = ttest(cods(:, 3) - cods(:, 1), [], 'Tail', 'Right');
        [~, pFulls] = ttest(cods(:, 4) - cods(:, 1:3), [], 'Tail', 'Right');
        [~, modelBest] = max(mean(cods));
        if saveOutputs
            f = fopen(fullfile(outdir, sprintf('%s_%s_compare_oddeven_cod_p.txt', chName, site)), 'w');
            fprintf(f, 'Best model = %d\n', modelBest);
            fprintf(f, 'Means: '); fprintf(f, '%0.01f%%, ', mean(cods)); fprintf(f, '\n');
            fprintf(f, 'p, img > simple = %0.3e\n', pImg);
            fprintf(f, 'p, isi > simple = %0.3e\n', pIsi);
            fprintf(f, 'p, full > simple, img, isi: '); fprintf(f, '%0.3e, ', pFulls); fprintf(f, '\n');
            fclose(f);
        end
        
    end

end

%% 2) Load the Broadband all model fits (copied form analysisBBFIR2) and make bargraphs

loaddir = fullfile('output', 'mforge'); % to load results from mforge. Use outdir for local
ylimsRsq = [-30, 50];

dsFactor = 24;
bbType = 'power';
errType = 'abs';
twinErr = 1; % configure our own twinErr for calculation

rsqTable = readtable(fullfile(outdir, sprintf('BBRsq_%s.tsv', sub)), 'FileType', 'text', 'Delimiter', '\t');

rsqTableSig = rsqTable(rsqTable.responsive == 1, :);
rsqTableSig = sortrows(rsqTableSig, 'rsq', 'descend'); % sort by highest to lowest rsq
fprintf('%d significant channels found\n', height(rsqTableSig));

for cc = 1:height(rsqTableSig)
    chName = rsqTableSig.name{cc};
    if any(strcmpi(chName, {'LOC1', 'LOC2', 'ROC1', 'ROC2'})), continue; end % skip the main pairs, bc we delete these VBBsegs
    
    for ss = 1:2
        site = sites{ss};

        fprintf('Loading odd-even BB fits for %s %s %s\n', sub, chName, site);

        vbbsegData = load(fullfile(outdir, sprintf('VBBseg_%s_%s_%s_ds-%d.mat', sub, chName, site, dsFactor)));
        mdl = load(fullfile(loaddir, sprintf('firBBOddEvenAllModels_%s_%s_%s_bbType-%s_errType-%s.mat', sub, chName, site, bbType, errType)));
        assert(size(vbbsegData.VBBsegPower, 2) == size(mdl.VBBseg, 2), 'Error: mismatch in trials between vbbsegData and mdl'); % we just need trials to match cuz we wanna use img and e2v

        % Pull out necessary fields
        VBBseg = mdl.VBBseg;
        ttseg = mdl.ttseg;
        trialsTrain = mdl.trialsTrain;
        trialsTest = mdl.trialsTest;
        twinVisual = mdl.config.twinvisual;
        twinCcep = mdl.config.twinccep;
        imgs = vbbsegData.imgs;
        e2vs = vbbsegData.e2vs;

        % transform to logpower, get test trials
        VBBseg = transformBBType(VBBseg, bbType, 'logpower');
        VBBsegTest = VBBseg(:, trialsTest);
        eventsTest = mdl.eventsSite(trialsTest, :);

        % calculate cods for each model, with 5th column being prediction from condition mean
        codsTest = zeros(length(trialsTest), 5);
        t2 = find(ttseg < twinErr(end), 1, 'last'); % last time point index to calculate with
        for ii = 1:4
            Vpred = mdl.mdlFits(ii).Vpred;
            Vpred = transformBBType(Vpred, bbType, 'logpower'); % also transform to log power
        
            for jj = 1:height(eventsTest)
        
                e2vDown = eventsTest.e2vDown(jj); % downsampled e2v ISI
                firstVisualInd = find(ttseg >= twinVisual(1), 1, 'first'); % First visual-predicted sample; >= matches behavior in fitFIR
        
                if isnan(e2vDown) % sham stim, start window at first twinVisual sample
                    t1 = firstVisualInd;
                else
                    t1 = min(firstVisualInd, find(ttseg >= twinCcep(1), 1, 'first') - e2vDown); % stim, use that time point or visual onset, whichever earlier
                end
        
                % Pull out the necessary segment of data, and calculate the COD and sigExp
                Y = VBBsegTest(t1:t2, jj);
                Ypred = Vpred(t1:t2, jj);
                %codsTest(jj, ii) = 100*(1 - norm(Y - Ypred)^2 / norm(Y - mean(Y))^2);
                codsTest(jj, ii) = 100*(1 - norm(Y - Ypred)^2 / norm(Y)^2); % use unnormalized COD bc a) y=0 has meaning and b) window is always changing so will give constantly different baselines
            end
        end

        % also calculate COD using condition mean
        VBBsegTrain = VBBseg(:, trialsTrain);
        eventsTrain = mdl.eventsSite(trialsTrain, :);
        for jj = 1:height(eventsTest)
            e2vType = eventsTest.e2vType(jj);
            imgType = eventsTest.trial_type{jj};
            if contains(imgType, 'zeroCoh'), imgType = 'zeroCoh'; end % collapse across the 2 image zero-coh conditions
        
            % mean of all training trials of this condition is the "predicted response"
            meanSigCond = mean(VBBsegTrain(:, eventsTrain.e2vType == e2vType & contains(eventsTrain.trial_type, imgType)), 2);
        
            e2vDown = eventsTest.e2vDown(jj); % downsampled e2v ISI
            firstVisualInd = find(ttseg >= twinVisual(1), 1, 'first'); % First visual-predicted sample; >= matches behavior in fitFIR
            if isnan(e2vDown) % sham trial, use visual onset
                t1 = firstVisualInd;
            else % stimulated, use first ccep sample or visual sample, whichever earlier
                t1 = min(firstVisualInd, find(ttseg >= twinCcep(1), 1, 'first') - e2vDown);
            end
        
            Y = VBBsegTest(t1:t2, jj);
            Ypred = meanSigCond(t1:t2); % from the mean condition signal
            %codsTest(jj, 5) = 100*(1 - norm(Y - Ypred)^2 / norm(Y - mean(Y))^2); % use 5th col
            codsTest(jj, 5) = 100*(1 - norm(Y - Ypred)^2 / norm(Y)^2); % use 5th col
        end

        % Plot bargraphs and do stats
        sampSz = size(codsTest, 1);
        sems = std(codsTest)/sqrt(sampSz); % SE of mean
        
        % bargraph with SD error bars
        figure('Position', [200, 200, 400, 600]); hold on
        bar(mean(codsTest), 'FaceColor', [0.5, 0.5, 0.5]);
        errorbar(mean(codsTest), sems, 'LineWidth', 1, 'LineStyle', 'none', 'CapSize', 10, 'Color', 'k'); % 95% conf intervals
        set(gca, 'xtick', 1:5, 'xticklabels', {'Simple', 'Img', 'ISI', 'Full', 'Mean'});
        xlabel('Model'); ylim(ylimsRsq); ylabel('COD');
        if saveOutputs
            saveas(gcf, fullfile(outdir, sprintf('%s_%s_compare_oddeven_BBcod_bbType-%s_%dms', chName, site, bbType, 1e3*twinErr(end))), 'png');
            saveas(gcf, fullfile(outdir, sprintf('%s_%s_compare_oddeven_BBcod_bbType-%s_%dms', chName, site, bbType, 1e3*twinErr(end))), 'svg');
        end
        
        fprintf('Odd-even COD by model, mean+/-SEM [--, %0.2fs]: [Simple]%0.01f%%+/-%0.01f%%, [Coh]%0.01f%%+/-%0.01f%%, [ISI]%0.01f%%+/-%0.01f%%, [Full]%0.01f%%+/-%0.01f%%\n', ...
            twinErr(end), mean(codsTest(:, 1)), sems(1), mean(codsTest(:, 2)), sems(2), mean(codsTest(:, 3)), sems(3), mean(codsTest(:, 4)), sems(4));
        
        % Stats: check if img, isi models are better than simple. Then check if full model is better than others
        [~, pImg] = ttest(codsTest(:, 2) - codsTest(:, 1), [], 'Tail', 'Right'); % img better than simple
        [~, pIsi] = ttest(codsTest(:, 3) - codsTest(:, 1), [], 'Tail', 'Right');
        [~, pFulls] = ttest(codsTest(:, 4) - codsTest(:, 1:3), [], 'Tail', 'Right');
        [~, modelBest] = max(mean(codsTest));
        if saveOutputs
            f = fopen(fullfile(outdir, sprintf('%s_%s_compare_oddeven_BBcod_p_bbType-%s_%dms.txt', chName, site, bbType, 1e3*twinErr(end))), 'w');
            fprintf(f, 'Best model = %d\n', modelBest);
            fprintf(f, 'Means: '); fprintf(f, '%0.01f%%, ', mean(codsTest)); fprintf(f, '\n');
            fprintf(f, 'p, img > simple = %0.3e\n', pImg);
            fprintf(f, 'p, isi > simple = %0.3e\n', pIsi);
            fprintf(f, 'p, full > simple, img, isi: '); fprintf(f, '%0.3e, ', pFulls); fprintf(f, '\n');
            fclose(f);
        end

    end

end

%% 3) Plot EP simple model for all pairs

loaddir = fullfile('output', 'mforge'); % to load results from mforge. Use outdir for local
dsFactor = 8;

rsqTable = readtable(fullfile(outdir, sprintf('BBRsq_%s.tsv', sub)), 'FileType', 'text', 'Delimiter', '\t');

rsqTableSig = rsqTable(rsqTable.responsive == 1, :);
rsqTableSig = sortrows(rsqTableSig, 'rsq', 'descend'); % sort by highest to lowest rsq
fprintf('%d significant channels found\n', height(rsqTableSig));

xlims = [-0.2, 1];
ylims = [-300, 300];

for cc = 1:height(rsqTableSig)
    chName = rsqTableSig.name{cc};
    
    for ss = 1:2
        site = sites{ss};

        mdl = load(fullfile(loaddir, sprintf('firFullErp_%s_%s_%s_fitISI-0_fitCoh-0.mat', sub, chName, site)));
    
        xMat = mdl.xMat;
        ttseg = mdl.ttseg;

        figure('Position', [200, 200, 500, 300]);
        subplot(1, 2, 1); hold on % plot CCEP
        yline(0);
        plot(ttseg, xMat(:, 1), 'k-', 'LineWidth', 1);
        hold off
        xlim(xlims); xticks(0:0.2:1); ylim(ylims);
        
        subplot(1, 2, 2); hold on
        yline(0);
        plot(ttseg, xMat(:, 2), 'LineWidth', 1);
        hold off
        xlim(xlims); xticks(0:0.2:1); ylim(ylims);

        if saveOutputs
            saveas(gcf, fullfile(outdir, sprintf('%s_%s_FIRResponses_fitISI=0_fitCoh=0', chName, site)), 'png');
            saveas(gcf, fullfile(outdir, sprintf('%s_%s_FIRResponses_fitISI=0_fitCoh=0', chName, site)), 'svg');
        end

        close(gcf);
    end
end

%% 4) Plot broadband simple model for all pairs

loaddir = fullfile('output', 'mforge'); % to load results from mforge. Use outdir for local

dsFactor = 24;
bbType = 'power';
errType = 'abs';

rsqTable = readtable(fullfile(outdir, sprintf('BBRsq_%s.tsv', sub)), 'FileType', 'text', 'Delimiter', '\t');

rsqTableSig = rsqTable(rsqTable.responsive == 1, :);
rsqTableSig = sortrows(rsqTableSig, 'rsq', 'descend'); % sort by highest to lowest rsq
fprintf('%d significant channels found\n', height(rsqTableSig));

xlims = [-0.2, 1];
ylims = [-1, 6];

for cc = 1:height(rsqTableSig)
    chName = rsqTableSig.name{cc};
    
    for ss = 1:2
        site = sites{ss};

        mdl = load(fullfile(loaddir, sprintf('firBBFull_%s_%s_%s_bbType-%s_fitISI-0_fitCoh-0_errType-%s.mat', sub, chName, site, bbType, errType)));
    
        xMat = mdl.xMat;
        ttseg = mdl.ttseg;

        figure('Position', [200, 200, 500, 300]);
        subplot(1, 2, 1); hold on % plot CCEP
        yline(0);
        plot(ttseg, xMat(:, 1), 'k-', 'LineWidth', 1);
        hold off
        xlim(xlims); xticks(0:0.2:1); ylim(ylims); yticks(0:2:10);
        
        subplot(1, 2, 2); hold on
        yline(0);
        plot(ttseg, xMat(:, 2), 'LineWidth', 1);
        hold off
        xlim(xlims); xticks(0:0.2:1); ylim(ylims); yticks(0:2:10);

        if saveOutputs
            saveas(gcf, fullfile(outdir, sprintf('%s_%s_FIRBBResponses_bbType=%s_fitISI=0_fitCoh=0_errType=%s', chName, site, bbType, errType)), 'png');
            saveas(gcf, fullfile(outdir, sprintf('%s_%s_FIRBBResponses_bbType=%s_fitISI=0_fitCoh=0_errType=%s', chName, site, bbType, errType)), 'svg');
        end

        close(gcf);
    end
end