createModelFigs.m

%% Create sample of how curves would look for each type of model, figure 4A
%
%    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/>.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
set(0, 'DefaultFigureRenderer', 'painters');

% Choose a set of fitted data to load so we can use its shapes as a starting point
sub = '1';
chName = 'LOC1';
site = 'LOC4-LOC5';
fitISI = true;
fitCoh = false;

saveOutputs = true;
outdir = fullfile('output', 'simulated');
loaddir = fullfile('output', 'mforge'); % to load results from mforge

% Configurations constant for all the plots
lw = 1; % line width of plots
ylims = [-250, 250];

%% Load data and pull out xMat

mdl = load(fullfile(loaddir, sprintf('firFullErp_%s_%s_%s_fitISI-%d_fitCoh-%d.mat', sub, chName, site, double(fitISI), double(fitCoh))));

ttseg = mdl.ttseg;
xMat = mdl.xMat;
figure; plot(ttseg, xMat);

% basic CCEP shape
SCcep0 = xMat(:, 1);


%% 1) Simple model

% basic visual response -- we'll use the 200 ms ISI prediction for this
SVisual0 = xMat(:, 4);

% process SCcep0 and use only part of it
SCcepTemp = processSig(ttseg, SCcep0);
SCcep1 = SCcep0;
SCcep1(ttseg > 0.1) = SCcepTemp(ttseg > 0.1);
SCcep1 = movmean(SCcep1, 20); % smooth out the kink
SCcep1 = [zeros(5, 1); SCcep1(1:end-5)]; % shift CCEP slightly to account for the anticausal rise caused by smoothing

% process visual
SVisual1 = processSig(ttseg, SVisual0);

% Plot components
figure('Position', [200, 200, 250, 300]);

% Plot CCEP 200 ms before visual
subplot(3, 2, 1); hold on
shiftSamps = sum(ttseg < (ttseg(1) + 0.2));
SCcep1Shift200 = circshift(SCcep1, -shiftSamps);
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(-0.2, 'Color', [1, 0.6, 0]); xline(0, 'Color', [0, 0, 1]);
plot(ttseg, SCcep1Shift200, 'k-', 'LineWidth', lw);
plot([0.1, 0.3], [ylims(2)-10, ylims(2)-10], 'k-', 'LineWidth', 3); % marker for length (200 ms)
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

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

% combined, with CCEP shifted 200 ms
subplot(3, 2, 5); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(-0.2, 'Color', [1, 0.6, 0]); xline(0, 'Color', [0, 0, 1]);
plot(ttseg, SVisual1 + SCcep1Shift200, 'k-', 'LineWidth', lw);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

% Now plot them at 0 ms
subplot(3, 2, 2); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(0, 'Color', [0, 0, 1]);
plot(ttseg, SCcep1, 'k-', 'LineWidth', lw);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

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

% combined, at 0 ms ISI
subplot(3, 2, 6); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(0, 'Color', [0, 0, 1]);
plot(ttseg, SVisual1 + SCcep1, 'k-', 'LineWidth', lw);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

if saveOutputs
    saveas(gcf, fullfile(outdir, 'erp_simulation_simple'), 'png');
    saveas(gcf, fullfile(outdir, 'erp_simulation_simple'), 'svg');
end

%% 2) Coherence model. For this we'll show the CCEP and the visual response, at a few different amplitudes corresponding to different images

% Plot components
figure('Position', [200, 200, 250, 300]);

SVisualLow = SVisual1*0.5;
SVisualHigh = SVisual1*1.5;

% Plot the CCEP again, now with 100 ms shifted so it doesn't convey a sense of ISI
subplot(3, 2, 3); hold on
shiftSamps = sum(ttseg < (ttseg(1) + 0.1));
SCcep1Shift100 = circshift(SCcep1, -shiftSamps);
yline(0, 'Color', [0.5, 0.5, 0.5]);
plot(ttseg, SCcep1Shift100, 'k-', 'LineWidth', lw);
xlim([-0.2, 0.5]); ylim(ylims);
axis off
hold off

% Low amplitude response
subplot(3, 2, 2); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
plot(ttseg, SVisualLow, 'k-', 'LineWidth', lw);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

% Normal amplitude response
subplot(3, 2, 4); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
plot(ttseg, SVisual1, 'k-', 'LineWidth', lw);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

% High amplitude response
subplot(3, 2, 6); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
plot(ttseg, SVisualHigh, 'k-', 'LineWidth', lw);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

if saveOutputs
    saveas(gcf, fullfile(outdir, 'erp_simulation_image'), 'png');
    saveas(gcf, fullfile(outdir, 'erp_simulation_image'), 'svg');
end

%% 3) ISI model. We'll make similar plots as for the simple model, but using a different response for 0 ms vs 200 ms

SVisualISI0 = xMat(:, 2); % we'll use the 0ms ISI visual response for this
SVisualISI1 = processSig(ttseg, SVisualISI0);

% Plot components
figure('Position', [200, 200, 250, 300]);

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

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

% combined, with CCEP shifted 200 ms
subplot(3, 2, 5); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(-0.2, 'Color', [1, 0.6, 0]); xline(0, 'Color', [0, 0, 1]);
plot(ttseg, SVisual1 + SCcep1Shift200, 'k-', 'LineWidth', lw);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

% Now plot them at 0 ms
subplot(3, 2, 2); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(0, 'Color', [0, 0, 1]);
plot(ttseg, SCcep1, 'k-', 'LineWidth', lw);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

% visual only. Plot both, with the previous one in gray (we'll make dashed in illustrator).
% Now the idea will be we point an arrow from gray to black
subplot(3, 2, 4); hold on
yline(0, 'Color', [0.3, 0.3, 0.3]);
xline(0, 'Color', [0, 0, 1]);
plot(ttseg, SVisual1, '-', 'Color', [0.5, 0.5, 0.5], 'LineWidth', lw);
plot(ttseg, SVisualISI1, 'k-', 'LineWidth', lw);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

% combined, at 0 ms ISI
subplot(3, 2, 6); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(0, 'Color', [0, 0, 1]);
plot(ttseg, SVisualISI1 + SCcep1, 'k-', 'LineWidth', lw);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

if saveOutputs
    saveas(gcf, fullfile(outdir, 'erp_simulation_ISI'), 'png');
    saveas(gcf, fullfile(outdir, 'erp_simulation_ISI'), 'svg');
end

%% 4) Full model. For this we'll show a 2x2 set of visual responses with ellipses in the middle

% Create a low and high response version of the 0ms ISI-modulated response
SVisualISILow = 0.5*SVisualISI1;
SVisualISIHigh = SVisualISI1;

% Plot components
figure('Position', [200, 200, 375, 200]); % make this 2x3 so the (1,1) plot is the CCEP

% The ccep, same as in the coherence model plot
subplot(2, 3, 1); hold on
shiftSamps = sum(ttseg < (ttseg(1) + 0.1));
SCcep1Shift100 = circshift(SCcep1, -shiftSamps);
yline(0, 'Color', [0.5, 0.5, 0.5]);
plot(ttseg, SCcep1Shift100, 'k-', 'LineWidth', lw);
xlim([-0.2, 0.5]); ylim(ylims);
axis off
hold off

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

% 0 ms ISI, modulated low coherence visual response
subplot(2, 3, 3); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(0, 'Color', [0, 0, 1]);
plot(ttseg, SVisualISILow, 'k-', 'LineWidth', lw);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

% 200 ms, high coherence visual response
subplot(2, 3, 5); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(-0.2, 'Color', [1, 0.6, 0]); xline(0, 'Color', [0, 0, 1]);
plot(ttseg, SVisualHigh, 'k-', 'LineWidth', lw);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

% 0 ms ISI, modulated high coherence visual response
subplot(2, 3, 6); hold on
yline(0, 'Color', [0.5, 0.5, 0.5]);
xline(0, 'Color', [0, 0, 1]);
plot(ttseg, SVisualISIHigh, 'k-', 'LineWidth', lw);
xlim([-0.2, 0.6]); ylim(ylims);
axis off
hold off

if saveOutputs
    saveas(gcf, fullfile(outdir, 'erp_simulation_full'), 'png');
    saveas(gcf, fullfile(outdir, 'erp_simulation_full'), 'svg');
end

%% 5) Stich a few responses together for the "b" matrix in FIR model part

% 3 different signal types concatenated
b = [SCcep1 + SVisual1; SVisualHigh; SCcep1Shift200 + SVisualISI1];
ttStitch = [ttseg, ttseg, ttseg]'; % stitched time points

figure('Position', [200, 200, 600, 100]);
plot(b, 'k-'); xlim([-inf, inf]);
xline([find(ttStitch == 0)], 'Color', 'r');
axis off
if saveOutputs
    saveas(gcf, fullfile(outdir, 'erp_simulation_3stitched'), 'png');
    saveas(gcf, fullfile(outdir, 'erp_simulation_3stitched'), 'svg');
end

%% Functions

function S = processSig(tt, S0)

    % lowpass
    d = designfilt('lowpassiir', ...        % Response type
       'PassbandFrequency',20, ...     % Frequency constraints
       'StopbandFrequency',40, ...
       'PassbandRipple',3, ...          % Magnitude constraints
       'StopbandAttenuation',30, ...
       'DesignMethod','butter', ...      % Design method
       'MatchExactly','passband', ...   % Design method options
       'SampleRate',600);               % Sample rate

    S = filtfilt(d, S0);

    S = movmean(S, 50);

    S = S .* exp(-2*tt(:));
end