analysisERPFIR2.m

%% Created to separate the parts of analysisERPFIR that require FITTING from then loading/summarizing
%   Created 03/03/2024
%
%    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/>.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% 1) Configure paths and which ERP fit to load

set(0, 'DefaultFigureRenderer', 'painters');

defaultCols = [0 0.4470 0.7410; 0.8500 0.3250 0.0980; 0.9290 0.6940 0.1250]; % for each ISI
defaultColsRev = defaultCols(end:-1:1, :);

% colormap for image coh condition
load('broc10.mat');
imgCondColors = [0, 0, 0; broc10([3, 4, 5, 8, 7, 6], :)];

sub = '1';
chName = 'LOC1';
site = 'LOC4-LOC5';
%site = 'LG6-LG7';

%sub = '2';
%chName = 'ROC1';
%chName = 'RY12';
%site = 'ROC6-ROC7';
%site = 'RMO3-RMO4';

outdir = fullfile('output', sprintf('sub-%s', sub), sprintf('%s_%s', chName, site));
mkdir(outdir);

saveOutputs = false;

%% 2a) Load one model, plot fitted components and the overall fit to the data

% Params for configuring which model to load
loaddir = fullfile('output', 'mforge'); % to load results from cluster
%loaddir = outdir; % to load from local
fitISI = true; % whether a separate independent predictor should be fit to each ISI
fitCoh = true; % whether a separate independent predictor should be fit to each img coherence condition

dsFactor = 8; % We need vsegData (data saved for model fitting) for the img, e2v info
try
    vsegData = load(fullfile(outdir, sprintf('Vseg_%s_%s_%s_ds-%d.mat', sub, chName, site, dsFactor)));
catch
    fprintf('Loading vsegData from global\n');
    vsegData = load(fullfile('output', sub, 'global', sprintf('Vseg_%s_%s_%s_ds-%d.mat', sub, chName, site, dsFactor)));
end
mdl = load(fullfile(loaddir, sprintf('firFullErp_%s_%s_%s_fitISI-%d_fitCoh-%d.mat', sub, chName, site, double(fitISI), double(fitCoh))));
assert(size(vsegData.Vseg, 2) == size(mdl.Vseg, 2), 'Error: mismatch in trials between vsegData and mdl'); % we just need trials to match cuz we wanna use img and e2v

% Pull out necessary fields
Vpred = mdl.Vpred;
Vseg = mdl.Vseg;
ttseg = mdl.ttseg;
imgs = vsegData.imgs;
e2vs = vsegData.e2vs;

ylims = [-450, 450];
ylimsFullModel = [-400, 200]; % what ylims to use when plotting the full model
%ylimsFullModel = ylims/2;
xlims = [-0.05, 0.6];
%xlims = [-0.2, 1];

figure('Position', [200, 200, 400, 300]);
subplot(1, 2, 1); hold on % plot CCEP
yline(0);
plot(mdl.ttseg, mdl.xMat(:, 1), 'k-', 'LineWidth', 1);
hold off
set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylims);
%xlabel('Time from Stim Onset'); ylabel('Voltage (\muV)');

if ~fitCoh % 1 visual or 3 visual conditions
    subplot(1, 2, 2); hold on
    yline(0);
    plot(mdl.ttseg, mdl.xMat(:, end:-1:2), 'LineWidth', 1);
    hold off
    set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylims);
    %xlabel('Time from Visual Onset');
elseif ~fitISI
    subplot(1, 2, 2); hold on
    hs = zeros(1, 7);
    yline(0);
    for ii = 1:7
        hs(ii) = plot(mdl.ttseg, mdl.xMat(:, ii+1), 'Color', imgCondColors(ii, :), 'LineWidth', 0.75);
    end
    hold off
    set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylims);
    legend(hs([1, 2, 5]), {'0 Coh', 'Img1', 'Img2'});
    %xlabel('Time from Visual Onset');
else
    subplot('Position', [0.52, 0.55, 0.35, 0.37]); hold on % plot on top-right, the by img averages
    yline(0);
    for ii = 1:7 % average each image across the ISIs
        plot(ttseg, mean(mdl.xMat(:, ii*3-1:ii*3+1), 2), 'Color', imgCondColors(ii, :), 'LineWidth', 0.75); % average across the 3 ISIs
    end
    hold off
    set(gca, 'xtick', []); xlim(xlims); ylim(ylimsFullModel);

    subplot('Position', [0.52, 0.11, 0.35, 0.37]); hold on % plot on bottom-right, the by ISI averages
    yline(0);
    for ii = 1:3 % average each ISI across the iamges
        plot(ttseg,  mean(mdl.xMat(:, ii+1:3:end), 2), 'Color', defaultColsRev(ii, :), 'LineWidth', 0.75); % average across the 7 images for each e2v. 0ms = yellow
    end
    hold off
    set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylimsFullModel);
end
if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_FIRResponses_fitISI=%s_fitCoh=%s', chName, site, string(fitISI), string(fitCoh))), 'png');
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_FIRResponses_fitISI=%s_fitCoh=%s', chName, site, string(fitISI), string(fitCoh))), 'svg');
end

% Plot for each condition with data and fitted response overlaid
figure('Position', [200, 200, 2000, 800]); tiledlayout(4, 7, 'TileSpacing', 'Compact');
for ee = 1:4 % tile e2v conditions vertically
    for ii = 1:7 % image conditions horizontally

        Vcurr = Vseg(:, e2vs(ee).idxSite & imgs(ii).idxSite);
        Vpredcurr = Vpred(:, e2vs(ee).idxSite & imgs(ii).idxSite);
        
        nexttile; hold on
        plot(ttseg, Vcurr, '-', 'Color', [0.5, 0.5, 0.5]);
        plot(ttseg, mean(Vcurr, 2), 'k-', 'LineWidth', 1.5);
        plot(ttseg, mean(Vpredcurr, 2), 'r-', 'LineWidth', 1.5);
        yline(0);
        xlim([-0.3, 1]); ylim([-600, 600]);
        title(sprintf('%s, e2v=%0.0fms', imgs(ii).label, 1e3*e2vs(ee).mode));
    end
end
if saveOutputs, saveas(gcf, fullfile(outdir, sprintf('FIRFitVsData_fitISI=%s_fitCoh=%s.png', string(fitISI), string(fitCoh)))); end

fprintf('COD/sigExp with fitISI=%d, fitCoh=%d model on [%0.2fs, %0.2fs] is %0.01f%%/%0.01f%%\n', ...
    double(fitISI), double(fitCoh), mdl.config.twincod(1), mdl.config.twincod(2), mdl.cod, mdl.sigExp);

%% 2b) Load a model, and plot a few conditions on the same axis (vertically)
% Used for Figure 4F, G

% Params for configuring which model to load
loaddir = fullfile('output', 'mforge'); % to load results from mforge
%loaddir = outdir; % to load from local
fitISI = true; % whether a separate independent predictor should be fit to each ISI
fitCoh = true; % whether a separate independent predictor should be fit to each img coherence condition

dsFactor = 8; % We need vsegData (data saved for model fitting) for the img, e2v info
vsegData = load(fullfile(outdir, sprintf('Vseg_%s_%s_%s_ds-%d.mat', sub, chName, site, dsFactor)));
mdl = load(fullfile(loaddir, sprintf('firFullErp_%s_%s_%s_fitISI-%d_fitCoh-%d.mat', sub, chName, site, double(fitISI), double(fitCoh))));
assert(size(vsegData.Vseg, 2) == size(mdl.Vseg, 2), 'Error: mismatch in trials between vsegData and mdl'); % we just need trials to match cuz we wanna use img and e2v

% Pull out necessary fields
Vpred = mdl.Vpred;
Vseg = mdl.Vseg;
ttseg = mdl.ttseg;
imgs = vsegData.imgs;
e2vs = vsegData.e2vs;

% Which e2v and img conditions to plot (match in length). E.g. 1,1 means zeroCoh, sham stim
e2vConds = [1, 4, 3, 2, 2]; % sham, 200ms, 100ms, 0ms, 0ms, feels more natural this way
imgConds = [1, 2, 6, 4, 7]; % 0coh, 25cohimg1, 50cohimg2, 100cohimg1, 100hcohimg2
condNames = arrayfun(@(x, y) sprintf('e2v-%d:img-%d', x, y), e2vConds, imgConds, 'UniformOutput', false);

% where to plot each condition on the axis
offsets = 800*(0:length(imgConds) - 1);
offsets = offsets(end:-1:1);

figure('Position', [200, 200, 300, 600]); hold on
for ii = 1:length(imgConds)
    Vcurr = Vseg(:, e2vs(e2vConds(ii)).idxSite & imgs(imgConds(ii)).idxSite);
    Vpredcurr = Vpred(:, e2vs(e2vConds(ii)).idxSite & imgs(imgConds(ii)).idxSite);

    yline(offsets(ii), 'Color', [0.2, 0.2, 0.2]);
    plot(ttseg, Vcurr + offsets(ii), '-', 'Color', [0.5, 0.5, 0.5]);
    plot(ttseg, mean(Vcurr, 2) + offsets(ii), 'k-', 'LineWidth', 1);
    plot(ttseg, mean(Vpredcurr, 2) + offsets(ii), 'r-', 'LineWidth', 1);
    xlim([-0.25, 1]);
end
hold off
set(gca, 'ytick', offsets(end:-1:1), 'yticklabels', condNames(end:-1:1));
if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('FIRFitVsDataSubset_fitISI=%s_fitCoh=%s', string(fitISI), string(fitCoh))), 'png');
    saveas(gcf, fullfile(outdir, sprintf('FIRFitVsDataSubset_fitISI=%s_fitCoh=%s', string(fitISI), string(fitCoh))), 'svg');
end

%% 3c) Load odd-even fitting results, plot bargraphs
% Figures 4H, J, S2

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

mdl = load(fullfile(loaddir, sprintf('firOddEvenAllModelsErp_%s_%s_%s.mat', sub, chName, site)));
%cods = mdl.codsTest;
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
t95 = tinv(0.975, sampSz - 1); % n-1 degrees of freedom, 97.5 percentile (2-tailed 95 conf)
confs = sems*t95; % size of half-conf interval

ylimsRsq = [0, 80];

% 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

%% 4) Plot bootstrapped visual responses along with significantly different intervals for ISI=true case

% Configs
%loaddir = outdir;
loaddir = fullfile('output', 'mforge'); % to load results from mforge. Use outdir for local
fitISI = true; % whether a separate independent predictor should be fit to each ISI
fitCoh = false; % whether a separate independent predictor should be fit to each img coherence condition
nBoots = 1000; % how many bootstrapped outputs to load

% load the first one to initialize matrix
temp = load(fullfile(loaddir, sprintf('firBootstrap_%s_%s_%s_fitISI-%d_fitCoh-%d', sub, chName, site, double(fitISI), double(fitCoh)), ...
            sprintf('firBootstrap_%s_%s_%s_fitISI-%d_fitCoh-%d_seed-1.mat', sub, chName, site, double(fitISI), double(fitCoh))));
initSz = size(temp.xMat);
xMatBoot = nan(initSz(1), initSz(2), nBoots);
for ii = 1:nBoots
    mdl = load(fullfile(loaddir, sprintf('firBootstrap_%s_%s_%s_fitISI-%d_fitCoh-%d', sub, chName, site, double(fitISI), double(fitCoh)), ...
            sprintf('firBootstrap_%s_%s_%s_fitISI-%d_fitCoh-%d_seed-%d.mat', sub, chName, site, double(fitISI), double(fitCoh), ii)));
    xMatBoot(:, :, ii) = mdl.xMat;
end
ttseg = mdl.ttseg;

% 200ms, 100ms, 0ms
visualBoot = xMatBoot(:, 4:-1:2, :);
xlims = [-0.05, 1];
%ylimsWide = [ylims(1), ylims(2)+100]; % make a bit wider to accommodate confidence intervals

figure('Position', [200, 200, 500, 300]);
subplot(1, 2, 1); hold on % plot CCEP
yline(0);
plotCurvConfSample(ttseg, squeeze(xMatBoot(:, 1, :))', 'k', 0.3, 95); % CCEP
plot(ttseg, mean(xMatBoot(:, 1, :), 3), 'k-', 'LineWidth', 0.75);
hold off
set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylims);
xlabel('Time from Stim Onset'); ylabel('Voltage (\muV)');

% Calculate significantly different intervals between 100ms and 0ms to the 200 ms condition
minDur = 0.02; % set a minimum duration of 20 ms for significant clusters
clusters100 = ttseg(clusterTestSampleNoperm(squeeze(visualBoot(:, 2, :)), squeeze(visualBoot(:, 1, :)), 99));
clusters0 = ttseg(clusterTestSampleNoperm(squeeze(visualBoot(:, 3, :)), squeeze(visualBoot(:, 1, :)), 99));
clusters0v100 = ttseg(clusterTestSampleNoperm(squeeze(visualBoot(:, 2, :)), squeeze(visualBoot(:, 3, :)), 99)); % also compare 0 and 100 ms conditions
clusters100(diff(clusters100, 1, 2) < minDur, :) = [];
clusters0(diff(clusters0, 1, 2) < minDur, :) = [];
clusters0v100(diff(clusters0v100, 1, 2) < minDur, :) = [];

subplot(1, 2, 2); hold on
yline(0);
for ii = 1:3
    plotCurvConfSample(ttseg, squeeze(visualBoot(:, ii, :))', defaultCols(ii, :), 0.3, 95);
    plot(ttseg, mean(visualBoot(:, ii, :), 3), 'Color', defaultCols(ii, :), 'LineWidth', 0.75);
end

% significant differences between the 100ms and 200ms conditions
for jj = 1:size(clusters100, 1)
    plot(clusters100(jj, :), ones(1, 2)*(ylims(2) - 20), 'Color', defaultCols(2, :), 'LineWidth', 2);
end
% significant differences between the 0ms and 200ms conditions
for jj = 1:size(clusters0, 1)
    plot(clusters0(jj, :), ones(1, 2)*(ylims(2) - 40), 'Color', defaultCols(3, :), 'LineWidth', 2);
end
% significant differences between the 0ms and 100ms conditions
for jj = 1:size(clusters0v100, 1)
    plot(clusters0v100(jj, :), ones(1, 2)*(ylims(2) - 60), 'Color', [0, 0.6, 0], 'LineWidth', 2);
end
hold off
set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylims);
xlabel('Time from Visual Onset'); ylabel('Voltage (\muV)');

if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_%s_FIRResponsesConf_fitISI=%s_fitCoh=%s_nboots=%d_mindur=%0dms', sub, chName, site, string(fitISI), string(fitCoh), nBoots, round(1e3*minDur))), 'png');
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_%s_FIRResponsesConf_fitISI=%s_fitCoh=%s_nboots=%d_mindur=%0dms', sub, chName, site, string(fitISI), string(fitCoh), nBoots, round(1e3*minDur))), 'svg');
end

%% 5) Add button press as predictor, and check if split-half R^2 improves, like in 3c

% where to load split-half results from
loaddir = fullfile('output', 'mforge'); % to load results from mforge. Use outdir for local
fitISI = true; % whether a separate independent predictor should be fit to each ISI
fitCoh = false; % whether a separate independent predictor should be fit to each img coherence condition

mdl = load(fullfile(loaddir, sprintf('firOddEvenToButtonErp_%s_%s_%s_fitISI-%d_fitCoh-%d.mat', sub, chName, site, fitISI, fitCoh)));
%cods = mdl.codsTest;
cods = mdl.sigExpsTest; % use unnormalized COD
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
t95 = tinv(0.975, sampSz - 1); % n-1 degrees of freedom, 97.5 percentile (2-tailed 95 conf)
confs = sems*t95; % size of half-conf interval

ylimsRsq = [0, 80];

% 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:2, 'xticklabels', {'-Button', '+Button'});
xlabel('Model'); ylim(ylimsRsq); ylabel('COD');
if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_fitISI=%d_fitCoh=%d_compare_oddeven_tobutton_cod', chName, site, mdl.config.fitisi, mdl.config.fitcoh)), 'png');
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_fitISI=%d_fitCoh=%d_compare_oddeven_tobutton_cod', chName, site, mdl.config.fitisi, mdl.config.fitcoh)), 'svg');
end

fprintf('Odd-even COD, mean+/-SEM [--, %0.2fs]: [-Button]%0.01f%%+/-%0.01f%%, [+Button]%0.01f%%+/-%0.01f%%\n', ...
    twinCod(2), mean(cods(:, 1)), sems(1), mean(cods(:, 2)), sems(2));

% Stats, ttest
[~, p] = ttest(cods(:, 2) - cods(:, 1), [], 'Tail', 'Right'); % right-tailed ttest to test if adding button is BETTER
fprintf('p = %0.3e (right-tailed ttest +Button > -Button)\n', p);