analysisBBFIR1.m

%% Run this script for broadband analysis after loading preprocessed data with loadPreprocessedMef.m
% This version of the script fits electrical plus visual at different ISI. Simplest model
% Goal is to see if we can still extract a canonical "visual response" for each image condition
%
%    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) Extract data from ONE channel of interest and ONE stim site of interest, plot all trials

saveOutputs = false;

% use this block to configure subject and stim site
switch sub
    case '2'
        chName = 'ROC2';
        %chName = 'ROC1-ROC2';
        
        site = 'RMO3-RMO4'; ylimsOrig = [-2 15]; ylimsTrialsOrig = [-1, 50]; % ylims are for plotting averages, ylimTrials are for plotting individual trials 
        %site = 'ROC6-ROC7'; ylimsOrig = [-2 15]; ylimsTrialsOrig = [-1, 50]; 
        
    case '1'
        chName = 'LOC2';
        %chName = 'LOC1-LOC2';
  
        %site = 'LOC4-LOC5'; ylimsOrig = [-5 20]; ylimsTrialsOrig = [-1, 50]; 
        %site = 'LG6-LG7'; ylimsOrig = [-2 15]; ylimsTrialsOrig = [-1, 50]; 
end

if contains(chName, '-') % bipolar given
    fprintf('Extracting bipolar channel %s\n', chName);
    chNum = find(strcmpi(bipolarNames, chName));
    assert(length(chNum) == 1, 'Error: Did not find exactly one Bipolar channel');
    V = squeeze(Mbip(chNum, :, :));
    VCcep = squeeze(dataCCEP(strcmp(stimSites, site)).Mbip(chNum, :, :)); % CCEP-only data
else
    fprintf('Extracting CAR channel %s\n', chName);
    chNum = find(strcmpi(channelsEphys.name, chName));
    assert(length(chNum) == 1, 'Error: Did not find exactly one CAR channel');
    V = squeeze(M(chNum, :, :));
    VCcep = squeeze(dataCCEP(strcmp(stimSites, site)).MCcep(chNum, :, :));
end
ttCcep = dataCCEP(1).ttCcep;

% separate outdir for stimpair to recording electrode pair
outdir = fullfile('output', sprintf('sub-%s', sub), sprintf('%s_%s', chName, site), 'broadband');
mkdir(outdir);

% keep only trials corresponding to stim site of interest or sham stim (no e2v), and remove nan trials in V
idxSite = strcmp(events.electrical_stimulation_site, site) | isnan(events.e2v);
idxSite = idxSite & ~any(isnan(V), 1)';
V = V(:, idxSite);
eventsSite = events(idxSite, :);
for ii = 1:length(e2vs) % do the same with e2v indices
    e2vs(ii).idxSite = e2vs(ii).idx(idxSite); % e2v indices for stim site of interest only
end 
for ii = 1:length(imgs) % do the same with img indices
    imgs(ii).idxSite = imgs(ii).idx(idxSite);
end

% create a single numerical column in eventsSite to indicate which e2v type
eventsSite.e2vType = zeros(height(eventsSite), 1); % 0 means no stim
for ii = 2:length(e2vs)
    eventsSite.e2vType(e2vs(ii).idxSite) = ii - 1;
end

% remove trials that did not fit into an e2v category, as well as those that appear additionally too noisy with bipolar re-referencing.
% We remove these even when analyzing monopolar to keep the results consistent across reference formats
badTrs = getBadTrsRound2(chName, site, 'broadband');
badE2vs = find(~isnan(eventsSite.e2v) & eventsSite.e2vType == 0); % trials that were too far away to assign a type
if ~isempty(badE2vs) || ~isempty(badTrs)
    warning('%d trials are additionally annotated as bad, removing from analysis', length(badTrs)); 
    warning('%d trials are missing e2v type, removing from analysis', length(badE2vs));
    badTogether = union(badTrs, badE2vs);
    eventsSite(badTogether, :) = [];
    V(:, badTogether)= [];
    for ii = 1:length(e2vs)
        e2vs(ii).idxSite(badTogether) = [];
    end
    for ii = 1:length(imgs)
        imgs(ii).idxSite(badTogether) = [];
    end
end

% remove line noise on CCEPVisual
V = removeLineNoise_SpectrumEstimation(V', srate, 'LF = 60, NH = 3, HW = 3', false)';
VCcep = removeLineNoise_SpectrumEstimation(VCcep', srate, 'LF = 60, NH = 3, HW = 3', false)';

% plot summary map for all conditions (for inspection, no need saving bc same as FP plot)
figure('Position', [200, 200, 2000, 800]); tiledlayout(4, 7, 'TileSpacing', 'Compact');
for ee = 1:4 % tile e2v conditions vertically
    for ii = 1:length(imgs) % image conditions horizontally
        Vcurr = V(:, e2vs(ee).idxSite & imgs(ii).idxSite); % current e2v, current img condition
        fprintf('%d trials for %s, e2v=%0.0fms\n', size(Vcurr, 2), imgs(ii).label, 1e3*e2vs(ee).mode);
        
        nexttile;
        
        plot(tt, Vcurr);
        xlim([-0.3, 1]); ylim([-600, 600]);
        title(sprintf('%s, e2v=%0.0fms', imgs(ii).label, 1e3*e2vs(ee).mode));
    end
end

%% 2) Remove artifact, hilbert transform to broadband, then downsample data (+/- optional smoothing first), and clip in prep for solving

downsampFactor = 24;
smoothWinSz = 0.05*srate; % How many samples to smooth before downsampling. 0 means no smoothing
twinTrial = [-0.3, 1]; % Total data window to solve on. -0.3 at beginning to allow for small diffs in ISI at -200ms stim, if 50 ms pre-stim rise is alotted for broadband
twinBaseline = [-0.25, -0.05]; % window to take baseline on for CCEPVisual, using non-stim trials only
twinBaselineCcep = [-0.5, -0.05]; % similar window to calculate baseline on for CCEP, but relative to stimulation onset
hilbertBands = [70, 90; 90, 110; 130, 150; 150, 170];

% V1: remove artifact by linear interpolation
idx0 = find(tt >= 0, 1); % where 0 is visual locked
nsampsArt = ceil(0.004*srate); % how many samples to remove artifact on
V1 = V;
for ii = 1:height(eventsSite)
    e2v = eventsSite.e2v(ii);
    if isnan(e2v), continue; end

    idxesRm = idx0-e2v:idx0-e2v+nsampsArt-1;
    valsInterp = interp1(idxesRm([1, end]), V(idxesRm([1, end]), ii), idxesRm);
    V1(idxesRm, ii) = valsInterp;
end

% V1 for CCEP
V1Ccep = VCcep;
idxesRmCcep = find(ttCcep >= 0, 1):find(ttCcep >= 0, 1)+nsampsArt-1;
for ii = 1:size(VCcep, 2)
    valsInterp = interp1(idxesRmCcep([1, end]), VCcep(idxesRmCcep([1, end]), ii), idxesRmCcep);
    V1Ccep(idxesRmCcep, ii) = valsInterp;
end

% subtract evoked responses (mean of each category). Assumption here is that differences in ISI within same category are negligible
V2 = V1;
for ee = 1:4
    for ii = 1:length(imgs)
        V2(:, e2vs(ee).idxSite & imgs(ii).idxSite) = V2(:, e2vs(ee).idxSite & imgs(ii).idxSite) - mean(V2(:, e2vs(ee).idxSite & imgs(ii).idxSite), 2);
    end
end
V2Ccep = V1Ccep - mean(V1Ccep, 2);
%V2Ccep = V1Ccep;

% VBB: calculate broadband in LOG POWER. We will do all the baseline corrects, downsampling in this space
VBB = ieeg_getHilbert(V2, hilbertBands, srate, 'logpower');
VBBCcep = ieeg_getHilbert(V2Ccep, hilbertBands, srate, 'logpower'); % for CCEP


% previous incorrect baseline normalization by sham of all runs
%BBbaseline = mean(VBB(tt >= twinBaseline(1) & tt < twinBaseline(2), isnan(eventsSite.e2v)), 'all'); % only use non-stim trials to estimate baseline.
%VBB = VBB - BBbaseline;

% Baseline correct separately by run, using the sham stim trials in each run
for rr = 1:max(eventsSite.run)
    BBbaseline = mean(VBB(tt >= twinBaseline(1) & tt < twinBaseline(2), isnan(eventsSite.e2v) & eventsSite.run == rr), 'all');
    VBB(:, eventsSite.run == rr) = VBB(:, eventsSite.run == rr) - BBbaseline;
end
BBbaselineCcep = mean(VBBCcep(ttCcep >= twinBaselineCcep(1) & ttCcep < twinBaselineCcep(2), :), 'all');
VBBCcep = VBBCcep - BBbaselineCcep;

% VBB2: optional smoothing before downsampling (in log space)
if smoothWinSz ~= 0
    VBB2 = movmean(VBB, smoothWinSz);
    VBB2Ccep = movmean(VBBCcep, smoothWinSz);
    fprintf('Log broadband data smoothed by %d samples\n', smoothWinSz);
else
    VBB2 = VBB;
    VBB2Ccep = VBBCcep;
end

% downsample data
VBBDs = downsample(VBB2, downsampFactor);
ttDs = downsample(tt, downsampFactor);
VBBDsCcep = downsample(VBB2Ccep, downsampFactor);
ttDsCcep = downsample(ttCcep, downsampFactor);
sratedown = srate/downsampFactor;
fprintf('Broadband data downsampled by %d\n', downsampFactor);
eventsSite.e2vDown = round(eventsSite.e2v/downsampFactor); % new e2v in samples after downsampling

% Extract the time window of ccepVisual data for analysis
VBBseg = VBBDs(ttDs >= twinTrial(1) & ttDs <= twinTrial(2), :);
ttseg = ttDs(ttDs >= twinTrial(1) & ttDs <= twinTrial(2));

% now lets save all 3 versions of data: logpower, power, and amplitude
% a) power
VBBsegPower = 10.^VBBseg - 1; % the -1 is so that 0 means no change. Units = fold increase
VBBDsPowerCcep = 10.^VBBDsCcep - 1;

% b) log power
VBBsegLogpower = VBBseg; % no change
VBBDsLogpowerCcep = VBBDsCcep;

% c) amplitude
VBBsegAmplitude = sqrt(10.^VBBseg) - 1; % units are fold change over the "Root-Geomean-Square" baseline
VBBDsAmplitudeCcep = sqrt(10.^VBBDsCcep) - 1;

% plot logpower tiled trials for all conditions (downsampled)
figure('Position', [200, 200, 2000, 800]); tiledlayout(4, 7, 'TileSpacing', 'Compact');
for ee = 1:4 % tile e2v conditions vertically
    for ii = 1:length(imgs) % image conditions horizontally
        
        nexttile;
        Vcurr = VBBDs(:, e2vs(ee).idxSite & imgs(ii).idxSite); % current e2v, current img condition
        hold on
        plot(ttDs, Vcurr);
        plot(ttDs, mean(Vcurr, 2), 'k-', 'LineWidth', 1.5);
        xlim([-0.5, 1]); ylim([-1, log10(ylimsTrialsOrig(2))]);
        yline(0, 'Color', 'k');
        hold off
        title(sprintf('%s, e2v=%0.0fms', imgs(ii).label, 1e3*e2vs(ee).mode));
    end
end
if saveOutputs, saveas(gcf, fullfile(outdir, sprintf('AllBBTrialsTiled_ds=%d_logpower.png', downsampFactor))); end

% plot CCEP broadband trials
figure; plot(ttDsCcep, VBBDsCcep); hold on
plot(ttDsCcep, mean(VBBDsCcep, 2), 'k-', 'LineWidth', 2);
hold off
xlim([-0.2, 1]); ylim([-1, log10(ylimsTrialsOrig(2))]);
yline(0, 'Color', 'k');
xlabel('Time from Stimulation (s)'); ylabel('Log power');
if saveOutputs, saveas(gcf, fullfile(outdir, sprintf('ccepBBTrials_ds=%d_logpower.png', downsampFactor))); end

% save the downsampled data so we can parallelize it (section 4)
if saveOutputs
    save(fullfile(outdir, sprintf('VBBseg_%s_%s_%s_ds-%d.mat', sub, chName, site, downsampFactor)), 'ttseg', 'VBBsegPower', 'VBBsegLogpower', 'VBBsegAmplitude', 'eventsSite', 'imgs', 'e2vs', 'sub', 'chName', 'site');
    fprintf('Saved downsampled VBBseg data power, logpower, and amplitude)\n');
end

return

%% Plot mean and trials of CCEP trials, in power and log power

% stats for log power tps different from 0
tt2Test = ttDsCcep >= -0.1 & ttDsCcep < 1;
[~, p, ~, stats] = ttest(VBBDsCcep'); % uncorrected test
h = zeros(1, length(p));
h(p < 0.05 & stats.tstat > 1) = 1;
h(p < 0.05 & stats.tstat < 1) = -1;

% plot CCEP broadband trials: log power
figure('Position', [200, 200, 250, 300]); hold on
%plot(ttDsCcep, VBBDsCcep, 'Color', [0.5, 0.5, 0.5]);
ieeg_plotCurvConf(ttDsCcep, VBBDsCcep', 'k', 0.3);
plot(ttDsCcep, mean(VBBDsCcep, 2), 'k-', 'LineWidth', 1);
plot(ttDsCcep(h == 1), 1, 'r.', 'MarkerSize', 10);
plot(ttDsCcep(h == -1), 1, 'b.', 'MarkerSize', 10);
hold off
xlim([-0.2, 1]);
ylim([-0.5, 1.5]);
%ylim([-1, log10(ylimsTrialsOrig(2))]);
yline(0, 'Color', 'k');
xlabel('Time from Stimulation (s)');
if saveOutputs, saveas(gcf, fullfile(outdir, sprintf('%s_%s_ccepBBTrials_ds=%d_logpower', chName, site, downsampFactor)), 'svg'); end



% % plot CCEP broadband trials
% figure('Position', [300, 200, 250, 300]); hold on
% %plot(ttDsCcep, VBBDsPowerCcep, 'Color', [0.5, 0.5, 0.5]);
% h = ieeg_plotCurvConf(ttDsCcep, log(VBBDsPowerCcep + 1)', 'k', 0.3); % plot log confidence interval
% set(h, 'YData', exp(h.YData) - 1);
% plot(ttDsCcep, geomean(VBBDsPowerCcep + 1, 2) - 1, 'k-', 'LineWidth', 1);
% hold off
% xlim([-0.2, 1]);
% ylim([-1, 5]);
% %ylim([-1, log10(ylimsTrialsOrig(2))]);
% yline(0, 'Color', 'k');
% xlabel('Time from Stimulation (s)');
% if saveOutputs, saveas(gcf, fullfile(outdir, sprintf('%s_%s_ccepBBTrials_ds=%d_power', chName, site, downsampFactor)), 'svg'); end



%% i) Just enough variables to test the model fitting parts without having to load all the data

sub = '1';
chName = 'LOC1';
site = 'LOC4-LOC5';
outdir = fullfile('output', sprintf('sub-%s', sub), sprintf('%s_%s', chName, site), 'broadband');

downsampFactor = 24;


%% 3) Configure and solve for predictors using abs or lsq. Can be done using compiled code

% Configure fit parameters
fitISI = false; % 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

dataPath = fullfile(outdir, sprintf('VBBseg_%s_%s_%s_ds-%d.mat', sub, chName, site, downsampFactor));
configPath = fullfile('compile_fitFiR', sprintf('configBB_fitISI-%d_fitCoh-%d.tsv', fitISI, fitCoh));

fitBBFull(dataPath, configPath, outdir);

%% 4b) Solve for predictors for all 3 bbTypes, splitting into odd/even for training/testing
% for both EVC electrodes in each subject, main and control stim sites
% This is to general the results for figure S6, comparing fitting on power vs logpower.

% We'll use this by default as the "simplest model"
fitISI = false;
fitCoh = true;

dataPath = fullfile(outdir, sprintf('VBBseg_%s_%s_%s_ds-%d.mat', sub, chName, site, downsampFactor));
configPath = fullfile('compile_fitFIR', sprintf('configBB_fitISI-%d_fitCoh-%d.tsv', fitISI, fitCoh));

fitBBOddEvenAllBBs(dataPath, configPath, outdir);

%% 5) Solve for predictors on odd-even split, using all 4 ISI/Coh models

dataPath = fullfile(outdir, sprintf('VBBseg_%s_%s_%s_ds-%d.mat', sub, chName, site, downsampFactor));
configPath = fullfile('compile_fitFIR', 'configBB_allModels.tsv');

fitBBOddEvenAllModels(dataPath, configPath, outdir);

%% 6) Bootstrap full model to estimate confidence intervals on solved time series. Can be done on mforge with compiled code

% Configure fit parameters (these locate a specific config file)
fitISI = false; % 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

% number of resamples to bootstrap
nBoots = 10;
xMatBoot = [];
tic;
for nn = 1:nBoots
    showdot(nn);
    dataPath = fullfile(outdir, sprintf('VBBseg_%s_%s_%s_ds-%d.mat', sub, chName, site, downsampFactor));
    configPath = fullfile('compile_fitFIR', sprintf('configBB_fitISI-%d_fitCoh-%d.tsv', fitISI, fitCoh));
    [x, xMat] = fitBBBootstrap(dataPath, configPath, outdir, num2str(nn));
end
fprintf('\n'); toc;

%% 7) Fit button press as an additional predictor and compare to without button press, for desired ISI/Coh combo

% Configure fit parameters, which one to compare with button press
fitISI = false;
fitCoh = false;

dataPath = fullfile(outdir, sprintf('VBBseg_%s_%s_%s_ds-%d.mat', sub, chName, site, downsampFactor));
configPath = fullfile('compile_fitFIR', sprintf('configBB_fitISI-%d_fitCoh-%d.tsv', fitISI, fitCoh));

fitBBOddEvenToButton(dataPath, configPath, outdir);