overlay_bsep_spesvisual.m

%% This script overlays the BSEPs with the SPES+Visual traces
% 04/2026
%
%    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 = true;

switch sub

    case '2'
        ch = 175; % ROC1
        %ch = 176; % ROC2
        %ch = 177; % ROC3
        %ch = 90; % 88-90: RY10-12
        
        site = 'RMO3-RMO4';
        %site = 'ROC6-ROC7';

    case '1'
        ch = 205; % LOC1
        %ch = 206; % LOC2
        %ch = 215; % 214-215: LOC10-11
        %ch = 101; % LY11
        %ch = 153; % 151-153: LV13-15
        

        %site = 'LOC4-LOC5';
        site = 'LG6-LG7';
end
chName = channelsEphys.name{ch};

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

V = squeeze(M(ch, :, :));
VCcep = squeeze(dataCCEP(strcmp(stimSites, site)).MCcep(ch, :, :)); % assumes dataCCEP is in same order as stimSites (constructed this way)
ttCcep = dataCCEP(1).ttCcep; % these are commented out bc we not doin anything w them

% 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 array 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
badTrs = getBadTrsRound2(chName, site, 'erp'); % remove additional bad trials caught here
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
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
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

%% Load FIR model results

% Get the EVI models
loaddir = fullfile('output', 'mforge'); % to load results from mforge
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

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', sprintf('sub-%s', sub), chName, site, dsFactor)));
mdl = load(fullfile(loaddir, sprintf('firFullErp_%s_%s_%s_fitISI-%d_fitCoh-%d.mat', sprintf('sub-%s', 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

ttseg = mdl.ttseg;
ccep_fit = mdl.xMat(:, 1);


%% Plot observed BSEPs, observed 200 ms SPES-Visual, and FIR model fit together

% Get all the trials with 200 ms EVI
V_200ms = V(:, e2vs(4).idxSite);
e2vs_this = eventsSite.e2v(e2vs(4).idxSite);

% align around stim time
idx0 = find(tt == 0);
tt_range = [-0.4, 0.8];
tt_index1 = find(tt >= tt_range(1), 1); % first index without adjustment
tt_stimonset = tt(tt >= tt_range(1) & tt < tt_range(2));

V_200ms_stimonset = nan(length(tt_stimonset), length(e2vs_this));
for ii = 1:length(e2vs_this)
    V_200ms_stimonset(:, ii) = V_200ms(tt_index1-e2vs_this(ii) : tt_index1-e2vs_this(ii) + length(tt_stimonset) - 1, ii);
end

figure('Position', [200, 200, 400, 400]); hold on

% plot trials or CI
%plot(ttCcep, VCcep, '-', 'Color', [0.7, 0.7, 0.7]);
%plot(tt_stimonset, V_200ms_stimonset, '-', 'Color', [0.7, 0.7, 1]);
ieeg_plotCurvConf(ttCcep, VCcep', 'k', 0.3);
ieeg_plotCurvConf(tt_stimonset, V_200ms_stimonset', [0.12, 0.47, 1], 0.3);

% means
plot(ttCcep, mean(VCcep, 2), 'k-', 'LineWidth', 1.5);
plot(tt_stimonset, mean(V_200ms_stimonset, 2), '-', 'Color', [0.12, 0.47, 1], 'LineWidth', 1.5);

% model fit
plot(ttseg, ccep_fit, '-', 'Color', [1.0, 0.55, 0], 'LineWidth', 1.5);

xline(0.2, 'r', 'LineWidth', 1);
yline(0, 'Color', [0.5, 0.5, 0.5], 'LineWidth', 1);
xlim([-0.1 0.6]);
ylim([-400, 400]);

if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('overlay_bsep_sv_%s_%s', chName, site)), 'png');
    saveas(gcf, fullfile(outdir, sprintf('overlay_bsep_sv_%s_%s', chName, site)), 'svg');
end