fig1_FP_broadband_examples.m

%% This script creates traces (taken from the data) to illustrate the EP and broadband data, and to compare additive vs multiplicative models
% Run this with analysisERPFIR1, after pulling data from a channel.
%
%    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/>.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% Plot a few illustrative traces for figure 1A

% from sub-2
tr = 126; % sub 2, used for figure

V1 = V(:, tr);
V2 = V(:, tr+1);

% mark visual and stim samples with 1 and 2, respectively. rest = 3
markers1 = zeros(size(tt));
markers2 = zeros(size(tt));
markers1(tt == 0) = 1;
markers2(tt == 0) = 1;
markers1(find(tt == 0) - eventsSite.e2v(tr)) = 2;
markers2(find(tt == 0) - eventsSite.e2v(tr + 1)) = 2;
markers1(tt == 1) = 3;
markers2(tt == 1) = 3;

trEnd1 = eventsSite.onset(tr) + 1.5; % the abs time associated with 1.5s after tr1
idxFirst = find(tt > trEnd1 - eventsSite.onset(tr + 1), 1, 'first'); % first index in V2 that is after end of V1

% combine the 2 adjacent trials and their stim/visual markers, using first 1.5s of first run
Vcombined = [V1(tt <= 1.5); V2(idxFirst:end)];
markersCombined = [markers1(tt <= 1.5), markers2(idxFirst:end)];
ttCombined = [tt(tt <= 1.5), tt(idxFirst:end) + eventsSite.onset(tr+1) - eventsSite.onset(tr)];

% downsample V
VcombinedDown = downsample(Vcombined, 8);
ttCombinedDown = downsample(ttCombined, 8);

restTimes = ttCombined(markersCombined == 3);

figure('Position', [200, 200, 800, 200]); hold on % for markers we can use non-downsampled vectors
xline(ttCombined(markersCombined == 1), 'Color', [0, 0, 0]); % visual
xline(ttCombined(markersCombined == 2), 'Color', [1, 0, 0]); % electrical
xline(restTimes, 'Color', [0, 1, 0]); % rest onset
plot(ttCombinedDown, VcombinedDown);
hold off
xlim([-0.2, restTimes(2)+0.1]); ylim([-550, 550]);
saveas(gcf, fullfile('output', sprintf('ccepVisualExample_%s_%s_%s_trs=%d-%d.svg', sub, chName, site, tr, tr+1)));

% Transform the same data to broadband

% first interpolate through artifacts
% V1: remove artifact by linear interpolation
idx0s = find(markersCombined == 2); % where 0 is visual locked
nsampsArt = ceil(0.01*srate); % how many samples to remove artifact on
Vcombined1 = Vcombined;
for ii = 1:2
    idxesRm = idx0s(ii):idx0s(ii)+nsampsArt-1;
    valsInterp = interp1(idxesRm([1, end]), Vcombined(idxesRm([1, end])), idxesRm);
    Vcombined1(idxesRm) = valsInterp;
end

hilbertBands = [70, 90; 90, 110; 130, 150; 150, 170];
VcombinedBB = ieeg_getHilbert(Vcombined1, hilbertBands, srate, 'power');
VcombinedBB = VcombinedBB / geomean(VcombinedBB(ttCombined > -0.5 & ttCombined < -0.25)); % baseline normalize
%VcombinedBB = log10(VcombinedBB); % log
VcombinedBB = VcombinedBB - 1;

% smooth broadband
VcombinedBB = movmean(VcombinedBB, 0.05*srate);

% downsample
VcombinedBBDown = downsample(VcombinedBB, 8);

figure('Position', [200, 400, 800, 200]); hold on % for markers we can use non-downsampled vectors
xline(ttCombined(markersCombined == 1), 'Color', [0, 0, 0]); % visual
xline(ttCombined(markersCombined == 2), 'Color', [1, 0, 0]); % electrical
xline(restTimes, 'Color', [0, 1, 0]); % rest onset
plot(ttCombinedDown, VcombinedBBDown);
hold off
xlim([-0.2, restTimes(2)+0.1]);
saveas(gcf, fullfile('output', sprintf('ccepVisualBroadBandExample_%s_%s_%s_trs=%d-%d.svg', sub, chName, site, tr, tr+1)));

% removed section on alpha here

%% Also create models of broadband signals summing or interacting (fig 1)

srate = 4800;
dsSim = 24; % downsample factor of data to load
ylims = [-1, 4];

% sub 1 deidentified here
mdl = load(fullfile('output', 'mforge', 'firBBFull_1_LOC1-LOC2_LOC4-LOC5_bbType-power_fitISI-0_fitCoh-1_errType-abs.mat'));

ttseg = mdl.ttseg;
sigCcepOrig = mdl.xMat(:, 1);
sigVisOrig = mean(mdl.xMat(:, 2:end), 2);

% smooth each signal a little bit
sigCcep = movmean(sigCcepOrig, 0.05*srate/dsSim);
sigVis = movmean(sigVisOrig, 0.05*srate/dsSim);

% Taper the response with sigmoidal curves
sigmoidfn = @(t, k, t0) 1./(1+exp(-k*(t-t0)));
env1 = 1 - sigmoidfn(ttseg, 2, 0.8);
env1 = env1 - env1(end); % so that the end is 0
env2 = 1 - sigmoidfn(ttseg, 4, 0.5); % less steep drop off that doesn't go to 0 for visual
sigCcep = env1(:).*sigCcep;
sigVis = env2(:).*sigVis;

% CCEP
figure('Position', [200, 200, 250, 300]);
subplot(3, 2, 1); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(0, 'Color', [0, 0, 1]);
plot(ttseg, sigCcep, 'k-', 'LineWidth', 1);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

% visual
subplot(3, 2, 3); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(0, 'Color', [0, 0, 1]);
plot(ttseg, sigVis, 'k-', 'LineWidth', 1);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

% summative model
subplot(3, 2, 5); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(0, 'Color', [0, 0, 1]);
plot(ttseg, sigCcep+sigVis, 'k-', 'LineWidth', 1);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

% interactive model
subplot(3, 2, 6); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(0, 'Color', [0, 0, 1]);
plot(ttseg, sigCcep+sigVis, '-', 'Color', [0.5, 0.5, 0.5], 'LineWidth', 1);
plot(ttseg, sigCcep+1.5*sigVis, 'k-', 'LineWidth', 1);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

%saveas(gcf, fullfile('output', 'fig1_ccepVisualBroadBandModels.svg'));

%% Plot comparison of summative vs multiplicative model

figure('Position', [200, 200, 200, 300]);
subplot(3, 2, 1); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(0, 'Color', [0, 0, 1]);
plot(ttseg, sigCcep, 'k-', 'LineWidth', 1);
xlim([-0.1, 0.5]); ylim(ylims);
axis off
hold off

% visual
subplot(3, 2, 3); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(0, 'Color', [0, 0, 1]);
plot(ttseg, sigVis, 'k-', 'LineWidth', 1);
xlim([-0.1, 0.5]); ylim(ylims);
axis off
hold off

% summative model
subplot(3, 2, 5); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(0, 'Color', [0, 0, 1]);
plot(ttseg, sigCcep+sigVis, 'k-', 'LineWidth', 1);
xlim([-0.1, 0.5]); ylim(ylims);
axis off
hold off

% multiplicative model
subplot(3, 2, 6); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(0, 'Color', [0, 0, 1]);
plot(ttseg, (sigCcep + 1).*(sigVis + 1) - 1, 'k-', 'LineWidth', 1); % make sure we are multiplying fold changes
xlim([-0.1, 0.5]); ylim(ylims);
axis off
hold off

saveas(gcf, fullfile('output', 'fig5_ccepVisualBroadBandAddVsMultiply.svg'));