checkSummative.m

%% For each stim/recording pair, check if EPs/BBs are summative or multiplicative
%
%    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

load('broc10.mat');
imgCondColors = [0, 0, 0; broc10([3, 4, 5, 8, 7, 6], :)]; % for the 7 img conds
defaultCols = [0 0.4470 0.7410; 0.8500 0.3250 0.0980; 0.9290 0.6940 0.1250]; % first 3 cols of matlab cm (matching ISIs)

saveOutputs = true;

outdir = fullfile('output', 'checkSummative');
mkdir(outdir);

loaddir = fullfile('output', 'mforge');

% configure percentage threshold and total duration to analyze signal on
durEP = 0.5; % 0.5s matches the COD calculations for EP
durBB = 0.5; % 0.5s for broadband too. if responses occur they should be predominantly in the first 500ms
minFrac = 0.2; % minimum fraction of time points required. 0.2 means 100 ms out of first 500 ms

%% 1) Check EVOKED POTENTIAL summation for ch-site pairs

%sub = '1';
%ch = 'LOC2';
%ch = 'LY11';
%ch = 'LV15';
%site = 'LOC4-LOC5';
%site = 'LG6-LG7';

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

% for this pair, get if the best fit was simple vs image model. Note we can assume fitISI=0 because these are for non-interactive conditions
fitCohEP = getFitCohEP(ch, site);
fprintf('Loading EP fitCoh=%s bootstraps for %s, %s, %s\n', string(fitCohEP), sub, ch, site);

% load the bootstrapped EP fits
nBoots = 1000;
temp = load(fullfile(loaddir, sprintf('firBootstrap_%s_%s_%s_fitISI-0_fitCoh-%d', sub, ch, site, double(fitCohEP)), ...
            sprintf('firBootstrap_%s_%s_%s_fitISI-0_fitCoh-%d_seed-1.mat', sub, ch, site, double(fitCohEP))));
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-0_fitCoh-%d', sub, ch, site, double(fitCohEP)), ...
            sprintf('firBootstrap_%s_%s_%s_fitISI-0_fitCoh-%d_seed-%d.mat', sub, ch, site, double(fitCohEP), ii)));
    xMatBoot(:, :, ii) = mdl.xMat;
end
ttseg = mdl.ttseg;

% plot the bootstrap fits to make sure it matches our figures
xlims = [-0.2, 1];
ylims = [-400, 400];
figure('Position', [200, 200, 500, 300]);
subplot(1, 2, 1); hold on
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);
if ~fitCohEP % simple model
    subplot(1, 2, 2); hold on
    yline(0);
    plotCurvConfSample(ttseg, squeeze(xMatBoot(:, 2, :))', defaultCols(1, :), 0.3, 95); % CCEP
    plot(ttseg, mean(xMatBoot(:, 2, :), 3), 'Color', defaultCols(1, :), 'LineWidth', 0.75);
    hold off
    set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylims);
else % img model
    subplot(1, 2, 2); hold on
    hs = zeros(1, 7);
    yline(0);
    for ii = 1:7
        hs(ii) = plot(ttseg, mean(xMatBoot(:, ii+1, :), 3), 'Color', imgCondColors(ii, :), 'LineWidth', 0.75);
    end
    hold off
    set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylims);
end

% Quantify CCEP percent of tp significantly different from 0
conf95Ccep = prctile(squeeze(xMatBoot(:, 1, :)), [2.5, 97.5], 2);
conf95CcepWin = conf95Ccep(ttseg >= 0 & ttseg < durEP, :); % on window of interest for broadband
fracSigCcep = sum(conf95CcepWin(:, 1) > 0 | conf95CcepWin(:, 2) < 0) / size(conf95CcepWin, 1);

% Quantify Visual percent of tp significantly different from 0
if fitCohEP % if image model, we assume there was a visual response and just set fraction to 1 (all tps) by default
    fracSigVis = 1;
else
    conf95Vis = prctile(squeeze(xMatBoot(:, 2, :)), [2.5, 97.5], 2);
    conf95VisWin = conf95Vis(ttseg >= 0 & ttseg < durEP, :); % on window of interest for broadband
    fracSigVis = sum(conf95VisWin(:, 1) > 0 | conf95VisWin(:, 2) < 0) / size(conf95CcepWin, 1);
end

% BOTH frac sig CCEP and Visual need to be >= threshold for this to be summative
isSummative = fracSigCcep >= minFrac && fracSigVis >= minFrac;
fprintf('EP to %0.2f: Frac sig CCEP = %0.2f, Frac sig Visual = %0.2f, summative = %s\n', durEP, fracSigCcep, fracSigVis, upper(string(isSummative)));

% save a text output
if saveOutputs
    f = fopen(fullfile(outdir, sprintf('checkSummativeEP_%s_%s_%s.txt', sub, ch, site)), 'w');
    fprintf(f, 'durEP = %0.2f, minFrac = %0.2f\n', durEP, minFrac);
    fprintf(f, 'Frac sig CCEP = %0.2f\n', fracSigCcep);
    fprintf(f, 'Frac sig Visual = %0.2f\n', fracSigVis);
    fprintf(f, 'Summative = %s\n', string(isSummative));
    fclose(f);
end

%% 2) Check BROADBAND summation for ch-site pairs

%sub = '1';
%ch = 'LOC2';
%ch = 'LY11';
%ch = 'LV15';
%site = 'LOC4-LOC5';
%site = 'LG6-LG7';

sub = '2';
%ch = 'ROC2';
ch = 'RY12';
%site = 'ROC6-ROC7';
site = 'RMO3-RMO4';

% for this pair, get if the best fit was simple vs image model. Note we can assume fitISI=0 because these are for non-interactive conditions
fitCohBB = getFitCohBB(ch, site);
fprintf('Loading BB fitCoh=%s bootstraps for %s, %s, %s\n', string(fitCohBB), sub, ch, site);

% load the bootstrapped fits
nBoots = 200;
temp = load(fullfile(loaddir, sprintf('firBBBootstrap_%s_%s_%s_bbType-power_fitISI-0_fitCoh-%d_errType-abs', sub, ch, site, double(fitCohBB)), ...
            sprintf('firBBBootstrap_%s_%s_%s_bbType-power_fitISI-0_fitCoh-%d_errType-abs_seed-1.mat', sub, ch, site, double(fitCohBB))));
initSz = size(temp.xMat);
xMatBoot = nan(initSz(1), initSz(2), nBoots);
for ii = 1:nBoots
    mdl = load(fullfile(loaddir, sprintf('firBBBootstrap_%s_%s_%s_bbType-power_fitISI-0_fitCoh-%d_errType-abs', sub, ch, site, double(fitCohBB)), ...
            sprintf('firBBBootstrap_%s_%s_%s_bbType-power_fitISI-0_fitCoh-%d_errType-abs_seed-%d.mat', sub, ch, site, double(fitCohBB), ii)));
    xMatBoot(:, :, ii) = mdl.xMat;
end
ttseg = mdl.ttseg;

% plot the bootstrap fits to make sure it matches our figures
xlims = [-0.2, 1];
ylims = [-1, 10]; % for power
figure('Position', [200, 200, 500, 300]);
subplot(1, 2, 1); hold on
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);
if ~fitCohBB % simple model
    subplot(1, 2, 2); hold on
    yline(0);
    plotCurvConfSample(ttseg, squeeze(xMatBoot(:, 2, :))', defaultCols(1, :), 0.3, 95); % CCEP
    plot(ttseg, mean(xMatBoot(:, 2, :), 3), 'Color', defaultCols(1, :), 'LineWidth', 0.75);
    hold off
    set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylims);
else % img model
    subplot(1, 2, 2); hold on
    hs = zeros(1, 7);
    yline(0);
    for ii = 1:7
        hs(ii) = plot(ttseg, mean(xMatBoot(:, ii+1, :), 3), 'Color', imgCondColors(ii, :), 'LineWidth', 0.75);
    end
    hold off
    set(gca, 'xtick', 0:0.2:1); xlim(xlims); ylim(ylims);
end

% Quantify CCEP percent of tp significantly different from 0
conf95Ccep = prctile(squeeze(xMatBoot(:, 1, :)), [2.5, 97.5], 2);
conf95CcepWin = conf95Ccep(ttseg >= 0 & ttseg < durBB, :); % on window of interest for broadband
fracSigCcep = sum(conf95CcepWin(:, 1) > 0 | conf95CcepWin(:, 2) < 0) / size(conf95CcepWin, 1);

% Quantify Visual percent of tp significantly different from 0
if fitCohBB % if image model, we assume there was a visual response and just set fraction to 1 (all tps) by default
    fracSigVis = 1;
else
    conf95Vis = prctile(squeeze(xMatBoot(:, 2, :)), [2.5, 97.5], 2);
    conf95VisWin = conf95Vis(ttseg >= 0 & ttseg < durBB, :); % on window of interest for broadband
    fracSigVis = sum(conf95VisWin(:, 1) > 0 | conf95VisWin(:, 2) < 0) / size(conf95CcepWin, 1);
end

% BOTH frac sig CCEP and Visual need to be >= threshold for this to be summative
isSummative = fracSigCcep >= minFrac && fracSigVis >= minFrac;
fprintf('BB to %0.2f: Frac sig CCEP = %0.2f, Frac sig Visual = %0.2f, summative = %s\n', durBB, fracSigCcep, fracSigVis, upper(string(isSummative)));

% save a text output
if saveOutputs
    f = fopen(fullfile(outdir, sprintf('checkSummativeBroadband_%s_%s_%s.txt', sub, ch, site)), 'w');
    fprintf(f, 'durBB = %0.2f, minFrac = %0.2f\n', durBB, minFrac);
    fprintf(f, 'Frac sig CCEP = %0.2f\n', fracSigCcep);
    fprintf(f, 'Frac sig Visual = %0.2f\n', fracSigVis);
    fprintf(f, 'Summative = %s\n', string(isSummative));
    fclose(f);
end

%% Functions

% evoked potentials
function fitCoh = getFitCohEP(ch, site)

    if strcmpi(ch, 'LOC1') && strcmpi(site, 'LG6-LG7')
        fitCoh = false;
    elseif strcmpi(ch, 'LOC2') && strcmpi(site, 'LG6-LG7')
        fitCoh = true;
    elseif strcmpi(ch, 'ROC1') && strcmpi(site, 'RMO3-RMO4')
        fitCoh = false;
    elseif strcmpi(ch, 'ROC2') && strcmpi(site, 'RMO3-RMO4')
        fitCoh = false;

    % global electrodes
    elseif strcmpi(ch, 'LOC10') && strcmpi(site, 'LG6-LG7') % re5
        fitCoh = true; % need to save, but visually obvious there is a response in both CCEP and visual
    elseif strcmpi(ch, 'LOC11') && strcmpi(site, 'LOC4-LOC5') % re6
        fitCoh = false;
    elseif strcmpi(ch, 'LY11') && strcmpi(site, 'LOC4-LOC5') % re 7
        fitCoh = false;
    elseif strcmpi(ch, 'LY11') && strcmpi(site, 'LG6-LG7')
        fitCoh = false;
    elseif strcmpi(ch, 'LV13') && strcmpi(site, 'LOC4-LOC5') % re 8
        fitCoh = false;
    elseif strcmpi(ch, 'LV13') && strcmpi(site, 'LG6-LG7')
        fitCoh = false;
    elseif strcmpi(ch, 'LV14') && strcmpi(site, 'LOC4-LOC5') % re 9
        fitCoh = false;
    elseif strcmpi(ch, 'LV14') && strcmpi(site, 'LG6-LG7')
        fitCoh = false;
    elseif strcmpi(ch, 'LV15') && strcmpi(site, 'LOC4-LOC5') % re 10
        fitCoh = false;
    elseif strcmpi(ch, 'LV15') && strcmpi(site, 'LG6-LG7')
        fitCoh = false;
    elseif strcmpi(ch, 'RY10') && strcmpi(site, 'ROC6-ROC7') % re11
        fitCoh = false;
    elseif strcmpi(ch, 'RY10') && strcmpi(site, 'RMO3-RMO4')
        fitCoh = false;
    elseif strcmpi(ch, 'RY11') && strcmpi(site, 'ROC6-ROC7') % re12
        fitCoh = false;
    elseif strcmpi(ch, 'RY11') && strcmpi(site, 'RMO3-RMO4')
        fitCoh = false;
    elseif strcmpi(ch, 'RY12') && strcmpi(site, 'ROC6-ROC7') % re13
        fitCoh = false;
    elseif strcmpi(ch, 'RY12') && strcmpi(site, 'RMO3-RMO4')
        fitCoh = false;
    else
        error('This is not a non-interactive ch/site pair');
    end

end

% broadband
function fitCoh = getFitCohBB(ch, site)

    if strcmpi(ch, 'LOC1') && strcmpi(site, 'LOC4-LOC5')
        fitCoh = false;
    elseif strcmpi(ch, 'LOC2') && strcmpi(site, 'LOC4-LOC5')
        fitCoh = false;
    elseif strcmpi(ch, 'LOC1') && strcmpi(site, 'LG6-LG7')
        fitCoh = false;
    elseif strcmpi(ch, 'LOC2') && strcmpi(site, 'LG6-LG7')
        fitCoh = true;
    elseif strcmpi(ch, 'ROC1') && strcmpi(site, 'ROC6-ROC7')
        fitCoh = true;
    elseif strcmpi(ch, 'ROC2') && strcmpi(site, 'ROC6-ROC7')
        fitCoh = false;
    elseif strcmpi(ch, 'ROC1') && strcmpi(site, 'RMO3-RMO4')
        fitCoh = true;
    elseif strcmpi(ch, 'ROC2') && strcmpi(site, 'RMO3-RMO4')
        fitCoh = true;

    % global channels
    elseif strcmpi(ch, 'LOC10') && strcmpi(site, 'LOC4-LOC5') % re5
        fitCoh = false;
    elseif strcmpi(ch, 'LOC10') && strcmpi(site, 'LG6-LG7')
        fitCoh = false;
    elseif strcmpi(ch, 'LOC11') && strcmpi(site, 'LOC4-LOC5') % re6
        fitCoh = false;
    elseif strcmpi(ch, 'LOC11') && strcmpi(site, 'LG6-LG7')
        fitCoh = false;
    elseif strcmpi(ch, 'LY11') && strcmpi(site, 'LOC4-LOC5') % re7
        fitCoh = false;
    elseif strcmpi(ch, 'LY11') && strcmpi(site, 'LG6-LG7')
        fitCoh = false;
    elseif strcmpi(ch, 'LV13') && strcmpi(site, 'LOC4-LOC5') % re8
        fitCoh = false;
    elseif strcmpi(ch, 'LV13') && strcmpi(site, 'LG6-LG7')
        fitCoh = false;
    elseif strcmpi(ch, 'LV14') && strcmpi(site, 'LOC4-LOC5') % re9
        fitCoh = false;
    elseif strcmpi(ch, 'LV14') && strcmpi(site, 'LG6-LG7')
        fitCoh = false;
    elseif strcmpi(ch, 'LV15') && strcmpi(site, 'LOC4-LOC5') % re10
        fitCoh = false;
    elseif strcmpi(ch, 'LV15') && strcmpi(site, 'LG6-LG7')
        fitCoh = false;
    elseif strcmpi(ch, 'RY10') && strcmpi(site, 'ROC6-ROC7') % re11
        fitCoh = false;
    elseif strcmpi(ch, 'RY10') && strcmpi(site, 'RMO3-RMO4')
        fitCoh = false;
    elseif strcmpi(ch, 'RY11') && strcmpi(site, 'ROC6-ROC7') % re12
        fitCoh = false;
    elseif strcmpi(ch, 'RY11') && strcmpi(site, 'RMO3-RMO4')
        fitCoh = true;
    elseif strcmpi(ch, 'RY12') && strcmpi(site, 'ROC6-ROC7') % re13
        fitCoh = false;
    elseif strcmpi(ch, 'RY12') && strcmpi(site, 'RMO3-RMO4')
        fitCoh = false;
    else
        error('This is not a non-interactive ch/site pair');
    end

end