main.m

%% Main script
%  Walks through code used to generate all manuscript results and figures, by section.
%
%    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/>.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%% Dependencies:
% mnl_ieegBasics: https://github.com/MultimodalNeuroimagingLab/mnl_ieegBasics
% knkutils: https://github.com/cvnlab/knkutils
% scientific color maps 8.0: https://www.fabiocrameri.ch/colourmaps/ (used as temporary plotting placeholders)
% distributionplot: (included) https://www.mathworks.com/matlabcentral/fileexchange/23661-violin-plots-for-plotting-multiple-distributions-distributionplot-m
%
%%% Keeping track of de-identified scripts:
% plotRecStimElecs.m
% tractography/ccepVisual_RenderTracks01.m
% preprocessFromMef
% preprocessCCEPForCCEPVisual.m
% loadPreprocessedMef.m
% analysisERPFIR1.m
% analysisERPFIR2.m
% fig1_FP_broadband_examples.m
% createModelFigs.m
% analysisBBFIR1.m
% compareBBTypes.m
% analysisBBFIR2.m
% globalAnalysis1.m
% globalAnalysis2.m
% checkSummative.m
% analysisPsychometric.m
% writeXmlCcepVisualDetPhase.m
% figure4A.m
% evaluateContrastShrinkage.m
%
%% Setup/paths
%
% - Be in the working directory containing this script
% - Download the data files from OpenNeuro (ds006485) and copy them to a "data" folder in the working directory (./data)
% - Note: code in ./compile_fitFIR meant to be copied to and executed on HPC cluster (GCP)
%

% Paths to some of the common dependencies
addpath('path/to/mnl_ieegBasics');
addpath('path/to/knkutils/mri');
addpath('path/to/scientific/colormaps')
addpath('external/distributionPlot');

% Local paths
addpath('functions');
addpath('functions/picturePrep');
addpath('compile_fitFIR'); % to test and run locally (e.g., in analysisERPFIR1.m)
addpath(genpath('tractography'));

%% 1) Task preparation
% Code involved in preparing task stimuli and stimuli order

% Produces the stimuli (noise-masked) from the original NSD color images. The resulting stimuli are included in data/derivatives/visual_stimuli
% Also produces the xml files that determine trial order
writeXmlCcepVisualDetPhase.m

% This script was used to find optimal contrast shrinkage toward the mean (0.3) before phase scrambling
evaluateContrastShrinkage.m

%% 2) Anatomical outputs
% Creates panels in Figures 2, S4A, S5A (overlays of recording, stimulating sites on brains)

% Plots inflated brains and coronal/axial slices
plotRecStimElecs.m

% Tractography code from Dora in ./tractography (for Figure 2A, D)
ccepVisual_RenderTracks01.m

%% 3) iEEG preprocessing

% Preprocess the SPES+Visual data.
% This merges the visual and electrical events files (i.e., creates derivatives/event_annotations/sub-X/sub-X_eventsCcepVisualRun0X.tsv),
% and doing so requires manual annotating of trial status (found in the "status" and "status_description" columns).
% This also saves the preprocessed iEEG data as 'output/sub-X/ccepVisualPreprocRun0X.mat' for each subject, run.
preprocessFromMef.m % for preprocessing

% Processes the SPES-only experimental data
% This also saves the preprocessed iEEG data as 'output/sub-X/CCEP_[stim_site]_preproc.mat' for each subject and stim site
preprocessCCEPForCCEPVisual.m

%% 4) Analyze SPES+Visual evoked potential data for EVC electrodes
% Figures 4, S1, S2

% Load the preprocessed data generated in section 3 into the workspace
loadPreprocessedMef.m

% Saves BSEPs in Figures 2C, F, and S1
% Prepares and downsamples SPES+Visual data for use in cluster (different file for each sub/stimsite/recordchannel combo)
analysisERPFIR1.m

% Calculates and saves outputs from BSEP-VEP model fitting. This includes Fig 4 and panels in S2. Also checks influence of button press
analysisERPFIR2.m

% Plots Figure 4A
figure4A.m

%% 5) Schematics for model figures
% Figures 1, 3

% This plots demonstrative traces in figure 1A, C
% Requires first loading data with analysisERPFIR1.m in section 4
% Note for 1C this creates the broadband traces only.
% Also creates the traces in S6A (broadband additive in power vs. log power)
fig1_FP_broadband_examples.m

% This plots the models in figure 3
% Also used this for EP trace in figure 1C
% also requires saving the necessary models (sub-1, LOC1, LOC4-5) via analysisERPFIR1.m first
createModelFigs.m

%% 6) Analyze SPES+Visual broadband data for EVC electrodes
% Figures 5, S3, S6

% As in section 4, load preprocessed data first into workspace
loadPreprocessedMef.m

% Prepares and downsamples SPES+Visual data for use in cluster, analogously to analysisERPFIR1
% Contains blocks that test running the compiled functions on the cluster
analysisBBFIR1.m

% compares the different broadband types to evaluate best fit (power vs. log power, figure S6)
% requires running models on cluster for all EVC electrode pairs
compareBBTypes.m

% Creates outputs and plots of broadband model fitting analysis, including testing button press
% Plots for Figures 5 and S3
analysisBBFIR2.m

%% 7) Analysis for non-EVC electrodes
% Table 1, Figures S4, S5

% As in section 4, load preprocessed data first into workspace
loadPreprocessedMef.m

% determine which electrodes are visually responsive, save PSD for all electrodes, save EP and BB downsampled model data
% Saves figures S4B, S5B
globalAnalysis1.m

% like in figure 2, use this to save inflated pial surface and slice plots for each visually responsive electrode
plotRecStimElecs.m

% Loops through all visually responsive electrodes and saves the EP, BB outputs (bargraphs, traces in S4C-D, S5C-D, stats)
% This saves all the simple model traces only. In cases where other models fit best, use analysisERPFIR2 and analysisBBFIR2 to save those traces
globalAnalysis2.m

% Checks whether, for each visually responsive electrode, EPs and BBs were additive or interactive. Requires manual iteration
% Generates results for Table 1
% Note helper function at bottom records manual calculation of whether simple model vs image model fit best (based on globalAnalysis2 stats)
checkSummative.m

%% 8) Behavioral (psychometric) analysis
% Figures S7, S8; tables S9-S10

% Update from resubmission: mdl5 is the combined logistic regression model on accuracy (S10)
analysisPsychometric.m

%% 9) New analyses added for 2026/04 resubmission

% S11 Figure - directly compare indepenendent SPES data, 200 ms EVI SPES+Visual, and model fits on same plot
overlay_bsep_spesvisual.m

% S12 Figure - time-frequency analysis using continuous wavelet transform
spectral_analysis.m