plotRecStimElecs.m

%% Use this script to plot recording and stimulation pair electrodes in the two subjects, with seegview and inflated gifti
% For Fig 2, S4A, S5A
%
%    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

addpath('path/to/freesurfer_code');

% path to freesurfer data
fsPath = fullfile('data', 'derivatives', 'freesurfer');

saveOutputs = true;

%sub = '1'; hemi = 'L'; recs = [132, 133]; % LOC1, LOC2
sub = '2'; hemi = 'R'; recs = [197, 198]; % ROC1, ROC2

% global electrodes
%sub = '1'; hemi = 'L'; recs = [132, 133, 141, 142, 38, 169, 170, 171]; % LOC1, LOC2, LOC10-11, LY11, LV13-15
%sub = '2'; hemi = 'R'; recs = [197, 198, 122, 123, 124]; % ROC1, ROC2, RY10, RY11, RY12

giiInfPath = fullfile(fsPath, sprintf('sub-%s', sub), sprintf('inflated.%s.surf.gii', upper(hemi)));
pialPath = fullfile(fsPath, sprintf('sub-%s', sub), sprintf('pial.%s.surf.gii', upper(hemi)));
whitePath = fullfile(fsPath, sprintf('sub-%s', sub), sprintf('white.%s.surf.gii', upper(hemi)));
sulcPath = fullfile(fsPath, sprintf('sub-%s', sub), 'sulc', sprintf('%sh.sulc', lower(hemi)));

if strcmpi(sub, '1') % use the no-contrast t1w image named differently
    niiPath = 'data/spesVisual_deidentified/data/derivatives/acpc/sub-1/sub-1_ses-mri01_T1w_deFaced.nii';
else
    niiPath = 'data/spesVisual_deidentified/data/derivatives/acpc/sub-2/sub-2_ses-mri01_T1w_acpc_deFaced.nii';
end
elecsPath = fullfile('data', sprintf('sub-%s', sub), 'ses-ieeg01', 'ieeg', sprintf('sub-%s_ses-ieeg01_electrodes.tsv', sub));

elecs = readtable(elecsPath, 'FileType', 'text', 'Delimiter', '\t');

%% Open sliceGUI to save slice at stim sites

opts = struct();
opts.elecDist = 4;
opts.elecNames = 'on';

sliceGUI(niiPath, elecsPath, opts);

%% Plot inflated brain with Benson or Wang Atlas and recording electrodes overlaid


% Read area labels and get cmap
namesAtlas = {'V1','V2','V3','hV4','VO1','VO2','LO1','LO2','TO1','TO2','V3b','V3a','FEF'};
cmAtlas = make_BensonColormap; cmAtlas = [cmAtlas; .2 .5 .2]; % add a row of color for FEF
surface_labels = MRIread(fullfile(fsPath, sprintf('sub-%s', sub), 'surf', sprintf('%sh.benson14_varea.mgz', lower(hemi)))); % use most labels from benson
wang_labels = MRIread(fullfile(fsPath, sprintf('sub-%s', sub), 'surf', sprintf('%sh.wang15_mplbl.mgz', lower(hemi))));
surface_labels.vol(wang_labels.vol == 25) = 13; % add wang label for FEF

vert_labels = surface_labels.vol(:);
assert(length(namesAtlas) == size(cmAtlas, 1));

% theta of medial view, determined from hemi
th = choose(strcmpi(hemi, 'r'), -90, 90);

% flip theta if you want lateral
%th = -th;

% load giftis
giiInf = gifti(giiInfPath);
white = gifti(whitePath); % for electrode inflating
pial = gifti(pialPath);
sulc = read_curv(sulcPath);

% inflate electrode and popout
xyzsRec = [elecs.x(recs), elecs.y(recs), elecs.z(recs)];
namesRec = elecs.name(recs);
xyzsInfRec = zeros(size(xyzsRec));

for ii = 1:length(recs)
    minDist2Pial = min(vecnorm(xyzsRec(ii, :) - pial.vertices, 2, 2));
    minDist2White = min(vecnorm(xyzsRec(ii, :) - white.vertices, 2, 2));
    fprintf('Rec #%d dist to pial = %0.2f, to white = %0.2f\n', ii, minDist2Pial, minDist2White);
    if minDist2Pial < minDist2White
        xyzsInf = getXyzsInf(elecs, hemi, pial, giiInf);
    else
        xyzsInf = getXyzsInf(elecs, hemi, white, giiInf);
    end
    xyzsInfPop = els_popout(xyzsInf, th, 0);
    xyzsInfRec(ii, :) = xyzsInfPop(recs(ii), :);
end

figure('Position', [200, 200, 1000, 800])
ieeg_RenderGiftiLabels(giiInf, vert_labels, cmAtlas, namesAtlas, sulc);
hold on
plot3(xyzsInfRec(:, 1), xyzsInfRec(:, 2), xyzsInfRec(:, 3), 'o', 'MarkerSize', 7, 'MarkerFaceColor', 'r', 'Color', [0.99, 0.99, 0.99]);
text(xyzsInfRec(:, 1), xyzsInfRec(:, 2), xyzsInfRec(:, 3), namesRec, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom');
hold off
ieeg_viewLight(th, 0);

recsString = join(elecs.name(recs), '-'); % join the recording electrode names for savin
if saveOutputs, saveas(gcf, fullfile('output', sub, sprintf('pialInf_%s_%s_th=%d.png', atlasChoice, recsString{1}, th))); end