diff --git a/README.txt b/README.txt
index a2f4765..bc722ba 100644
--- a/README.txt
+++ b/README.txt
@@ -1,7 +1,7 @@
If this code is used in a publication, please cite the manuscript:
"CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, GA Worrell, KJ Miller, and D Hermes.
-A preprint is available currently at doi: ___
+A preprint is available currently at https://doi.org/10.48550/arXiv.2310.00185.
Correspondence:
H Huang: huang.harvey@mayo.edu; D Hermes: hermes.dora@mayo.edu
diff --git a/figS1a_simCCEPTotalChs.m b/figS1a_simCCEPTotalChs.m
new file mode 100644
index 0000000..8903b1f
--- /dev/null
+++ b/figS1a_simCCEPTotalChs.m
@@ -0,0 +1,171 @@
+%% This script is derived from simCCEPLoop, except an outer loop varies the number of total channels
+% We simulate 30 reps for each combination of total channels and responsiveness and test CARLA's accuracy
+%
+% 2024/02/12
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 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 .
+%
+%% Control parameters of simulated data
+
+srate = 4800;
+tt = -0.5:1/srate:1-1/srate;
+
+ntrs = 12; % number of trials
+reps = 30; % number of reps for each condition
+
+cmSens = [1, 165/255, 0]; % use orange as color for the sensitive optimum
+
+Aglobal = 0;
+
+% Total number of channels to simulate for each condition
+nChsToTest = [25, 20, 15, 10];
+
+
+%% Outer loop through number of total channels
+
+rng(nChsToTest(1)); % Set seed depending on what number of channels we want, for reproducibility
+
+for nn = 1:length(nChsToTest)
+
+ nchs = nChsToTest(nn);
+ fprintf('Number of channels = %d\n', nchs);
+
+ %% Inner loop through number of responsive channels
+ for nresp = 0:0.2*nchs:0.8*nchs % 0 to 80% responsiveness, by 20% intervals
+
+ %% Create data
+
+ fprintf('%d of %d channels responsive\n', nresp, nchs);
+
+ outdir = fullfile('output', 'simLoopNChs', sprintf('nchs%d-%d', nchs, nresp));
+ mkdir(outdir);
+
+ chsResp = 1:nresp; % first x channels responsive, for ease. order doesn't matter. will sort and color them red
+
+ %%
+ for rr = 1:reps % multiple repetitions at each significance
+ fprintf('.');
+
+ % i) artifact only
+ V0 = zeros(length(tt), nchs);
+ Aart = 50 + rand(nchs, 1)*5; % slightly different artifact amplitudes for each channel
+ artifact = sin(2*pi*600*tt)';
+ artifact(tt < 0 | tt > 0.002) = 0;
+ V0 = V0 + artifact*Aart';
+
+ % A) Add the evoked potentials
+ A = 100;
+ V1 = V0;
+ sig = genRandSig(tt, length(chsResp), A);
+ V1(:, chsResp) = V0(:, chsResp) + sig;
+
+ % B) Option to add a global noise
+ if Aglobal
+ sigCommon = genRandSig(tt, 1, Aglobal);
+ else
+ sigCommon = 0;
+ end
+ V2 = V1 + sigCommon;
+
+ % C) Add common noise to all channels at each trial
+ V3 = repmat(V2', 1, 1, ntrs); % ch x time points x trial
+ phLN = rand(ntrs, 3)*2*pi; % LN phases
+ LN = zeros(length(tt), ntrs);
+ for ii = 1:ntrs
+ LN(:, ii) = 8*sin(2*pi*60*tt - phLN(ii, 1)) + 2*sin(2*pi*120*tt - phLN(ii, 2)) + 1*sin(2*pi*180*tt - phLN(ii, 3));
+ end
+
+ % brown noise shared across channels
+ BN = cumsum(0.4*randn(2*length(tt), ntrs)); % variable brown noise coefficient now
+ BN = ieeg_highpass(BN, srate, true);
+ BN = BN((0.5*length(tt)+1) : 1.5*length(tt), :);
+
+ noiseCommon = LN + BN;
+ V3 = V3 + shiftdim(noiseCommon, -1);
+
+ % D) add random brown noise
+ noiseRand = cumsum(0.4*randn(nchs, 2*length(tt), ntrs), 2); % give double the number of time points so we can highpass it
+ for ii = 1:nchs
+ noiseRand(ii, :, :) = ieeg_highpass(squeeze(noiseRand(ii, :, :)), srate, true);
+ end
+ noiseRand = noiseRand(:, (0.5*length(tt)+1) : 1.5*length(tt), :);
+ V4 = V3 + noiseRand;
+
+ %% Apply CARLA and plot outputs
+
+ [Vout, CAR, stats] = CARLA(tt, V4, srate, true); % get the sensitive output
+
+ % number of channels used for the CAR
+ nCAR = length(stats.chsUsed);
+ [~, nCARglob] = max(mean(stats.zMinMean, 2)); % number of channels at global maximum
+
+ % Plot average zmin across trials
+ figure('Position', [200, 200, 400, 300], 'Visible', 'off'); hold on
+ errorbar(mean(stats.zMinMean, 2), std(stats.zMinMean, 0, 2), 'k-o');
+ plot(nCARglob, mean(stats.zMinMean(nCARglob, :), 2), 'b*'); % global max as blue
+ if nCARglob ~= nCAR; plot(nCAR, mean(stats.zMinMean(nCAR, :), 2), '*', 'Color', cmSens); end
+ yline(0, 'Color', 'k');
+ saveas(gcf, fullfile(outdir, sprintf('zmin_rep%d', rr)), 'png');
+ saveas(gcf, fullfile(outdir, sprintf('zmin_rep%d', rr)), 'svg');
+
+ % Sort and plot channels by increasing covariance, draw line at cutoff
+ V4MeanSorted = mean(V4(stats.order, :, :), 3);
+ respBool = antifind(chsResp, nchs);
+ respBool = respBool(stats.order); % logical array of where responsive channels are
+ cm = zeros(nchs, 3);
+ cm(respBool, 1) = 1; % make red
+ figure('Position', [200, 200, 250, 600], 'Visible', 'off');
+ yspace = 80;
+ ys = ieeg_plotTrials(tt, V4MeanSorted', yspace, [], cm, 'LineWidth', 1);
+ yline(ys(nCARglob)-yspace/2, 'Color', 'b', 'LineWidth', 1.5);
+ if nCARglob ~= nCAR; yline(ys(nCAR)-yspace/2, 'Color', cmSens, 'LineWidth', 1.5); end
+ xlim([-0.1, 0.5]); set(gca, 'xtick', [0, 0.5]);
+ xlabel('Time (s)'); ylabel('Channels');
+ saveas(gcf, fullfile(outdir, sprintf('chsSorted_rep%d', rr)), 'png');
+ saveas(gcf, fullfile(outdir, sprintf('chsSorted_rep%d', rr)), 'svg');
+
+ close all;
+
+ % Accuracy values. positive means responsive/excluded from CAR
+ % We keep these variables as named here, but note that FN and FP are now renamed RCM and NCM in the manuscript.
+ TP = sum(find(respBool) > nCAR); % responsive channels successfully excluded from CAR (above the cutoff)
+ TN = sum(find(~respBool) <= nCAR); % NR channels successfully below or at cutoff
+ FN = sum(find(respBool) <= nCAR); % responsive channels incorrectly included in CAR. *This matters most
+ FP = sum(find(~respBool) > nCAR); % NR channels incorrectly excluded from CAR
+
+ % same for the global threshold
+ TPglob = sum(find(respBool) > nCARglob);
+ TNglob = sum(find(~respBool) <= nCARglob);
+ FNglob = sum(find(respBool) <= nCARglob);
+ FPglob = sum(find(~respBool) > nCARglob);
+
+ fid = fopen(fullfile(outdir, sprintf('accuracy_rep%d.txt', rr)), 'w');
+ fprintf(fid, 'TP\t%d\nTN\t%d\nFN\t%d\nFP\t%d\n', TP, TN, FN, FP);
+ fprintf(fid, 'TPglob\t%d\nTNglob\t%d\nFNglob\t%d\nFPglob\t%d', TPglob, TNglob, FNglob, FPglob);
+ fclose(fid);
+
+ end
+
+ fprintf('\n');
+
+ end
+end
diff --git a/figS1b_summarizeSimCCEPTotalChs.m b/figS1b_summarizeSimCCEPTotalChs.m
new file mode 100644
index 0000000..ad0b4a7
--- /dev/null
+++ b/figS1b_summarizeSimCCEPTotalChs.m
@@ -0,0 +1,137 @@
+%% This script takes outputs from figS1a_simCCEPTotalChs, and summarizes the accuracy for each total channel size and level of responsiveness
+% Generates outputs for Figure S1A
+%
+% 2024/02/12
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 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 .
+%
+%% Configure parameters and directories
+
+parCm = brighten(parula(15), -0.5); % black is the 50 ch condition, yellow is lowest number of channels (10)
+parCm(1, :) = [0, 0, 0];
+parCm(4, :) = [0.4940 0.1840 0.5560];
+
+% Configure the name of the folders to load accuracy scores from.
+rootdir = 'simLoopNChs';
+rootdir50 = 'simLoopNoGlob'; % We do not regenerate the 50 Chs condition, just pull from previous outputs
+
+nChsToTest = [50, 25, 20, 15, 10]; % total number of channels
+fracResp = 0:0.2:0.8; % fraction of all channels responsive
+nReps = 30;
+
+%% Read and store accuracy values at each number of total channels, then save to file
+
+% each field is in order of forloop in simCCEP_testNChs: nchs x responsiveness x nReps
+accuracy = struct();
+accuracy.TP = nan(length(nChsToTest), length(fracResp), nReps);
+accuracy.TN = nan(length(nChsToTest), length(fracResp), nReps);
+accuracy.FN = nan(length(nChsToTest), length(fracResp), nReps);
+accuracy.FP = nan(length(nChsToTest), length(fracResp), nReps);
+accuracy.TPglob = nan(length(nChsToTest), length(fracResp), nReps); % for the global maxes
+accuracy.TNglob = nan(length(nChsToTest), length(fracResp), nReps);
+accuracy.FNglob = nan(length(nChsToTest), length(fracResp), nReps);
+accuracy.FPglob = nan(length(nChsToTest), length(fracResp), nReps);
+
+% First store the 50 ch condition by pulling from previous output folder
+for nR = 1:length(fracResp)
+ fprintf('.');
+ nResp = round(fracResp(nR)*50);
+ indir = fullfile('output', rootdir50, sprintf('nchs50-%d', nResp));
+ for rep = 1:nReps
+
+ T_acc = readtable(fullfile(indir, sprintf('accuracy_50-%d_rep%d.txt', nResp, rep)), 'FileType', 'text', 'Delimiter', '\t');
+ accuracy.TP(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TP'));
+ accuracy.TN(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TN'));
+ accuracy.FN(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FN'));
+ accuracy.FP(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FP'));
+ accuracy.TPglob(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TPglob'));
+ accuracy.TNglob(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TNglob'));
+ accuracy.FNglob(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FNglob'));
+ accuracy.FPglob(1, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FPglob'));
+
+ end
+end
+
+% now load for each total number of channels from rootdir
+for nC = 2:length(nChsToTest)
+ for nR = 1:length(fracResp)
+ fprintf('.');
+ nChs = nChsToTest(nC);
+ indir = fullfile('output', rootdir, sprintf('nchs%d-%d', nChs, round(fracResp(nR)*nChs)));
+ for rep = 1:nReps
+ T_acc = readtable(fullfile(indir, sprintf('accuracy_rep%d.txt', rep)), 'FileType', 'text', 'Delimiter', '\t');
+ accuracy.TP(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TP'));
+ accuracy.TN(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TN'));
+ accuracy.FN(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FN'));
+ accuracy.FP(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FP'));
+ accuracy.TPglob(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TPglob'));
+ accuracy.TNglob(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TNglob'));
+ accuracy.FNglob(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FNglob'));
+ accuracy.FPglob(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FPglob'));
+ end
+ end
+end
+fprintf('\n');
+save(fullfile('output', rootdir, 'accuracies_all.mat'), 'accuracy');
+
+%% Plot FN, FP vs resp traces for each noise level, for first-peak method, as bars
+
+load(fullfile('output', rootdir, 'accuracies_all.mat'), 'accuracy');
+
+% displacements for the bars
+disps = [-0.04, -0.02, 0, 0.02, 0.04];
+
+% Based on more sensitive, first-peak optimum
+figure('Position', [200, 200, 700, 800]);
+subplot(2, 1, 1); % FN = responsive channels included in the CAR
+hold on
+for nc = 1:size(accuracy.FN, 1) % starting at 2 excludes the 50
+ FN = squeeze(accuracy.FN(nc, :, :));
+ FN = FN ./ (fracResp'*nChsToTest(nc)); % normalize by number of responsive channels
+ FN(isnan(FN)) = 0; % when there are 0 responsive channels, change these nans to 0
+ iqr = prctile(FN, [25, 75], 2);
+ med = median(FN, 2);
+ bar(fracResp + disps(nc), med, 0.08, 'FaceColor', parCm(nc*3-2, :));
+ errorbar(fracResp + disps(nc), med, med - iqr(:, 1), iqr(:, 2) - med, 'Color', parCm(nc*3-2, :), 'LineStyle', 'none', 'LineWidth', 1, 'CapSize', 8);
+end
+hold off
+xlim([-0.05, 0.85]); ylim([-0.05, 1.05]);
+set(gca, 'xtick', 0:0.2:0.8, 'xticklabels', 0:20:80, 'ytick', 0:0.2:1, 'yticklabels', 0:20:100, 'box', 'off');
+ylabel('% Responsive Channels Missed');
+
+subplot(2, 1, 2);
+hold on
+for nc = 1:size(accuracy.FN, 1)
+ FP = squeeze(accuracy.FP(nc, :, :));
+ FP = FP ./ ((1-fracResp)'*nChsToTest(nc)); % normalize by number of non-responsive channels
+ iqr = prctile(FP, [25, 75], 2);
+ med = median(FP, 2);
+ bar(fracResp + disps(nc), med, 0.08, 'FaceColor', parCm(nc*3-2, :));
+ errorbar(fracResp + disps(nc), med, med - iqr(:, 1), iqr(:, 2) - med, 'Color', parCm(nc*3-2, :), 'LineStyle', 'none', 'LineWidth', 1, 'CapSize', 8);
+end
+hold off
+xlim([-0.05, 0.85]); ylim([-0.05, 1.05]);
+set(gca, 'xtick', 0:0.2:0.8, 'xticklabels', 0:20:80, 'ytick', 0:0.2:1, 'yticklabels', 0:20:100, 'box', 'off');
+xlabel('Fraction of Channels with Response');
+ylabel('% Non-responsive Channels Missed');
+saveas(gcf, fullfile('output', rootdir, 'FNFP_bars_acrossNChs'), 'png');
+saveas(gcf, fullfile('output', rootdir, 'FNFP_bars_acrossNChs'), 'svg');
diff --git a/figS2a_simCCEPSNRs.m b/figS2a_simCCEPSNRs.m
new file mode 100644
index 0000000..b361c6c
--- /dev/null
+++ b/figS2a_simCCEPSNRs.m
@@ -0,0 +1,284 @@
+%% This script is derived from simCCEPLoop, except an outer loop varies the amplitude of brown noise (mulitiplicative coefficient)
+% We simulate 30 reps for each combination of SNR and responsiveness and test CARLA's accuracy
+% Generate data with responses at 0, 10, 20, 30, 40, 50, 60, 70, 80 percent of all channels
+% The last section of this script saves examples of how the same signal looks at varying levels of SNR (Figure S2A)
+%
+% 2024/02/12
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 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 .
+%
+%% Control parameters of simulated data
+
+srate = 4800;
+tt = -0.5:1/srate:1-1/srate;
+
+nchs = 50; % number of simulated channels
+ntrs = 12; % number of trials
+reps = 30; % number of reps
+
+cmSens = [1, 165/255, 0]; % use orange as color for the sensitive optimum
+
+Aglobal = 0;
+
+% window to calculate SNR on. Use same window as for calculating responses
+snrWin = [0.01, 0.3];
+
+%% Outer loop through brown noise level
+
+% Set seed depending on where brown coef range is starting, for reproducibility when parallelizing.
+% This is configured in main
+rng(brownCoefRange(1));
+
+for brownCoef = brownCoefRange(1):0.1:brownCoefRange(2)
+
+ fprintf('Brown noise coefficient = %0.1f\n', brownCoef);
+
+ %% Inner loop through number of responsive channels
+ for nresp = 0:5:40 % 0 to 80% responsiveness
+
+ %% Create data
+
+ fprintf('%d of %d channels responsive\n', nresp, nchs);
+
+ outdir = fullfile('output', 'simLoopSNR', sprintf('noise%0.1f_nchs%d-%d', brownCoef, nchs, nresp));
+ mkdir(outdir);
+
+ chsResp = 1:nresp; % first x channels responsive, for ease. order doesn't matter. will sort and color them red
+
+ %%
+ for rr = 1:reps % multiple repetitions at each significance
+
+ % i) artifact only
+ V0 = zeros(length(tt), nchs);
+ Aart = 50 + rand(nchs, 1)*5; % slightly different artifact amplitudes for each channel
+ artifact = sin(2*pi*600*tt)';
+ artifact(tt < 0 | tt > 0.002) = 0;
+ V0 = V0 + artifact*Aart';
+
+ % A) Add the evoked potentials
+ A = 100;
+ V1 = V0;
+ sig = genRandSig(tt, length(chsResp), A);
+ V1(:, chsResp) = V0(:, chsResp) + sig;
+
+ % B) Option to add a global noise
+ if Aglobal
+ sigCommon = genRandSig(tt, 1, Aglobal);
+ else
+ sigCommon = 0;
+ end
+ V2 = V1 + sigCommon;
+
+ % C) Add common noise to all channels at each trial
+ V3 = repmat(V2', 1, 1, ntrs); % ch x time points x trial
+ phLN = rand(ntrs, 3)*2*pi; % LN phases
+ LN = zeros(length(tt), ntrs);
+ for ii = 1:ntrs
+ LN(:, ii) = 8*sin(2*pi*60*tt - phLN(ii, 1)) + 2*sin(2*pi*120*tt - phLN(ii, 2)) + 1*sin(2*pi*180*tt - phLN(ii, 3));
+ end
+
+ % brown noise shared across channels
+ BN = cumsum(brownCoef*randn(2*length(tt), ntrs)); % variable brown noise coefficient now
+ BN = ieeg_highpass(BN, srate, true);
+ BN = BN((0.5*length(tt)+1) : 1.5*length(tt), :);
+
+ noiseCommon = LN + BN;
+ V3 = V3 + shiftdim(noiseCommon, -1);
+
+ % D) add random brown noise
+ noiseRand = cumsum(brownCoef*randn(nchs, 2*length(tt), ntrs), 2); % give double the number of time points so we can highpass it
+ for ii = 1:nchs
+ noiseRand(ii, :, :) = ieeg_highpass(squeeze(noiseRand(ii, :, :)), srate, true);
+ end
+ noiseRand = noiseRand(:, (0.5*length(tt)+1) : 1.5*length(tt), :);
+ V4 = V3 + noiseRand;
+
+ %% Calculate SNR
+
+ if nresp == 0
+ Psig = nan;
+ PnoiseRand = squeeze(sum(noiseRand(chsResp, tt >= snrWin(1) & tt < snrWin(2), :).^2, 2) / diff(snrWin));
+ snr = nan;
+ else
+
+ % power of signal (for each responsive channel created)
+ Psig = sum(sig(tt >= snrWin(1) & tt < snrWin(2), :).^2)' / diff(snrWin); % express as per second
+
+ % power of the aperiodic noise (common brown noise + random noise), calculated separately for each trial
+ noiseSum = noiseRand + shiftdim(BN, -1);
+ PnoiseRand = squeeze(sum(noiseSum(chsResp, tt >= snrWin(1) & tt < snrWin(2), :).^2, 2) / diff(snrWin));
+
+ % calculate snr for each trial separately (same Psig for all trials at one channel)
+ snr = repmat(Psig, 1, ntrs) ./ PnoiseRand;
+
+ % for simplicity, we only vary and consider random noise in SNR, assuming that periodic noise can be mostly attenuated by filtering
+
+ end
+ fprintf('%0.2f ', geomean(snr, 'all')); % geometric mean since power is lognormal
+
+ %% Apply CARLA and plot outputs
+
+ [Vout, CAR, stats] = CARLA(tt, V4, srate, true); % get the sensitive output
+
+ % number of channels used for the CAR
+ nCAR = length(stats.chsUsed);
+ [~, nCARglob] = max(mean(stats.zMinMean, 2)); % number of channels at global maximum
+
+ % Plot average zmin across trials
+ figure('Position', [200, 200, 400, 300], 'Visible', 'off'); hold on
+ errorbar(mean(stats.zMinMean, 2), std(stats.zMinMean, 0, 2), 'k-o');
+ plot(nCARglob, mean(stats.zMinMean(nCARglob, :), 2), 'b*'); % global max as blue
+ if nCARglob ~= nCAR; plot(nCAR, mean(stats.zMinMean(nCAR, :), 2), '*', 'Color', cmSens); end
+ yline(0, 'Color', 'k');
+ saveas(gcf, fullfile(outdir, sprintf('zmin_rep%d', rr)), 'png');
+ saveas(gcf, fullfile(outdir, sprintf('zmin_rep%d', rr)), 'svg');
+
+ % Sort and plot channels by increasing covariance, draw line at cutoff
+ V4MeanSorted = mean(V4(stats.order, :, :), 3);
+ respBool = antifind(chsResp, nchs);
+ respBool = respBool(stats.order); % logical array of where responsive channels are
+ cm = zeros(nchs, 3);
+ cm(respBool, 1) = 1; % make red
+ figure('Position', [200, 200, 250, 600], 'Visible', 'off');
+ yspace = 80;
+ ys = ieeg_plotTrials(tt, V4MeanSorted', yspace, [], cm, 'LineWidth', 1);
+ yline(ys(nCARglob)-yspace/2, 'Color', 'b', 'LineWidth', 1.5);
+ if nCARglob ~= nCAR; yline(ys(nCAR)-yspace/2, 'Color', cmSens, 'LineWidth', 1.5); end
+ xlim([-0.1, 0.5]); set(gca, 'xtick', [0, 0.5]);
+ xlabel('Time (s)'); ylabel('Channels');
+ saveas(gcf, fullfile(outdir, sprintf('chsSorted_rep%d', rr)), 'png');
+ saveas(gcf, fullfile(outdir, sprintf('chsSorted_rep%d', rr)), 'svg');
+
+ close all;
+
+ % Accuracy values. positive means responsive/excluded from CAR
+ % We keep these variables as named here, but note that FN and FP are now renamed RCM and NCM in the manuscript.
+ TP = sum(find(respBool) > nCAR); % responsive channels successfully excluded from CAR (above the cutoff)
+ TN = sum(find(~respBool) <= nCAR); % NR channels successfully below or at cutoff
+ FN = sum(find(respBool) <= nCAR); % responsive channels incorrectly included in CAR. *This matters most
+ FP = sum(find(~respBool) > nCAR); % NR channels incorrectly excluded from CAR
+
+ % same for the global threshold
+ TPglob = sum(find(respBool) > nCARglob);
+ TNglob = sum(find(~respBool) <= nCARglob);
+ FNglob = sum(find(respBool) <= nCARglob);
+ FPglob = sum(find(~respBool) > nCARglob);
+
+ fid = fopen(fullfile(outdir, sprintf('accuracy_rep%d.txt', rr)), 'w');
+ fprintf(fid, 'TP\t%d\nTN\t%d\nFN\t%d\nFP\t%d\n', TP, TN, FN, FP);
+ fprintf(fid, 'TPglob\t%d\nTNglob\t%d\nFNglob\t%d\nFPglob\t%d', TPglob, TNglob, FNglob, FPglob);
+ fclose(fid);
+
+ % save snr info
+ save(fullfile(outdir, sprintf('snr_rep%d.mat', rr)), 'snr', 'Psig', 'PnoiseRand');
+
+ end
+
+ fprintf('\n');
+
+ end
+end
+
+%return
+
+%% Save a few examples of how the channels look with the same signal but variable SNR
+
+outdir = fullfile('output', 'simLoopSNR');
+
+rng(25); % a representative seed with the 3 channels showing comparable SNR to the population geomeans
+
+nchsEx = 4;
+chsResp = 1:3;
+
+% i) artifact only
+V0 = zeros(length(tt), nchsEx);
+Aart = 50 + rand(nchsEx, 1)*5; % slightly different artifact amplitudes for each channel
+artifact = sin(2*pi*600*tt)';
+artifact(tt < 0 | tt > 0.002) = 0;
+V0 = V0 + artifact*Aart';
+
+% A) Add the evoked potentials
+A = 100;
+V1 = V0;
+sig = genRandSig(tt, length(chsResp), A);
+V1(:, chsResp) = V0(:, chsResp) + sig;
+
+% no global noise, copy V1 over
+V2 = V1;
+
+% keep the line noise the same for all
+phLN = rand(ntrs, 3)*2*pi; % LN phases
+LN = zeros(length(tt), ntrs);
+for ii = 1:ntrs
+ LN(:, ii) = 8*sin(2*pi*60*tt - phLN(ii, 1)) + 2*sin(2*pi*120*tt - phLN(ii, 2)) + 1*sin(2*pi*180*tt - phLN(ii, 3));
+end
+
+% create the same series of random noise for all noise levels, just varying amplitude
+BNBase = randn(2*length(tt), ntrs);
+noiseRandBase = randn(nchsEx, 2*length(tt), ntrs);
+
+for brownCoef = 0.4:0.1:1.0
+
+ % C) Add common noise to all channels at each trial
+ V3 = repmat(V2', 1, 1, ntrs); % ch x time points x trial
+
+ % brown noise shared across channels
+ BN = cumsum(brownCoef*BNBase); % variable brown noise coefficient now
+ BN = ieeg_highpass(BN, srate, true);
+ BN = BN((0.5*length(tt)+1) : 1.5*length(tt), :);
+
+ noiseCommon = LN + BN;
+ V3 = V3 + shiftdim(noiseCommon, -1);
+
+ % D) add random brown noise
+ noiseRand = cumsum(brownCoef*noiseRandBase, 2); % give double the number of time points so we can highpass it
+ for ii = 1:nchsEx
+ noiseRand(ii, :, :) = ieeg_highpass(squeeze(noiseRand(ii, :, :)), srate, true);
+ end
+ noiseRand = noiseRand(:, (0.5*length(tt)+1) : 1.5*length(tt), :);
+ V4 = V3 + noiseRand;
+
+ % Calculate the SNRs in the example
+ Psig = sum(sig(tt >= snrWin(1) & tt < snrWin(2), :).^2)' / diff(snrWin); % express as per second
+ noiseSum = noiseRand + shiftdim(BN, -1);
+ PnoiseRand = squeeze(sum(noiseSum(chsResp, tt >= snrWin(1) & tt < snrWin(2), :).^2, 2) / diff(snrWin));
+ snr = repmat(Psig, 1, ntrs) ./ PnoiseRand;
+ snrByCh = geomean(snr, 2); % average across trials for each channel
+ fprintf('%0.2f, ', snrByCh(:)); fprintf('\n');
+
+ % Plot and save
+ yspace = 200;
+ figure('Position', [200, 200, 600, 600]);
+ ys = ieeg_plotTrials(tt, mean(V4, 3)', yspace, [], 'k', 'LineWidth', 1);
+ hold on
+ for jj = 1:4 % plot the individual trials
+ plot(tt, ys(jj) + squeeze(V4(jj, :, :)), 'Color', [0.5, 0.5, 0.5]);
+ end
+ hold off
+ ylim([ys(end)-yspace, yspace]);
+ kids = get(gca, 'Children');
+ set(gca, 'Children', [kids(end-1:-2:end-7); kids(1:end-8); kids(end:-2:end-6)]); % first (top) the means, then the trials, lastly the horzline
+ xlabel('Time from Stim. (s)');
+ saveas(gcf, fullfile(outdir, sprintf('exampleSNRs_noise%0.1f_nchs%d-%d.png', brownCoef, nchsEx, length(chsResp))));
+ saveas(gcf, fullfile(outdir, sprintf('exampleSNRs_noise%0.1f_nchs%d-%d.svg', brownCoef, nchsEx, length(chsResp))));
+ close(gcf);
+end
diff --git a/figS2b_summarizeSimCCEPSNRs.m b/figS2b_summarizeSimCCEPSNRs.m
new file mode 100644
index 0000000..34be583
--- /dev/null
+++ b/figS2b_summarizeSimCCEPSNRs.m
@@ -0,0 +1,134 @@
+%% This script takes outputs from figS2a_simCCEPSNRs, and summarizes the accuracy for each level of responsiveness and SNR
+% Generates outputs for Figure S2B
+%
+% 2024/02/12
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 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 .
+%
+%% Configure parameters and directories
+
+parCm = brighten(parula(21), -0.5);
+parCm(1, :) = [0, 0, 0]; % set first row to black for extra contrast
+parCm(4, :) = [0.4940 0.1840 0.5560]; % purple
+
+rootdir = 'simLoopSNR';
+nChs = 50; % total number of channels
+nReps = 30;
+noiseCoefs = 0.4:0.1:1.0; % which noise coefficient levels to load
+nResponsives = 0:5:40; % number of responsive channels, from 0 to 80% by 10% increments
+
+%% Read and store accuracy values at each subset of responsiveness, save to file
+
+% each field is in order of forloop in simCCEP_testSNRs: noise x responsiveness x nReps
+accuracy = struct();
+accuracy.TP = nan(length(noiseCoefs), length(nResponsives), nReps);
+accuracy.TN = nan(length(noiseCoefs), length(nResponsives), nReps);
+accuracy.FN = nan(length(noiseCoefs), length(nResponsives), nReps);
+accuracy.FP = nan(length(noiseCoefs), length(nResponsives), nReps);
+
+% similar structure for snrs but cell to accommodate differeing amounts
+snrs = cell(length(noiseCoefs), length(nResponsives), nReps);
+
+for nC = 1:length(noiseCoefs)
+ for nR = 1:length(nResponsives)
+
+ fprintf('.');
+
+ indir = fullfile('output', rootdir, sprintf('noise%0.1f_nchs%d-%d', noiseCoefs(nC), nChs, nResponsives(nR)));
+
+ for rep = 1:nReps
+
+ T_acc = readtable(fullfile(indir, sprintf('accuracy_rep%d.txt', rep)), 'FileType', 'text', 'Delimiter', '\t');
+ accuracy.TP(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TP'));
+ accuracy.TN(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'TN'));
+ accuracy.FN(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FN'));
+ accuracy.FP(nC, nR, rep) = T_acc.Var2(strcmp(T_acc.Var1, 'FP'));
+
+ temp = load(fullfile(indir, sprintf('snr_rep%d.mat', rep)));
+ snrs{nC, nR, rep} = temp.snr;
+
+ end
+
+ end
+end
+fprintf('\n');
+save(fullfile('output', rootdir, 'accuracies_snrs_all.mat'), 'accuracy', 'snrs');
+
+%% Plot FN, FP vs resp traces for each noise level
+
+load(fullfile('output', rootdir, 'accuracies_snrs_all.mat'), 'accuracy', 'snrs');
+
+cmSens = [1, 165/255, 0]; % orange = sensitive threshold
+
+% Based on more sensitive, first-peak optimum
+figure('Position', [200, 200, 600, 900]);
+subplot(2, 1, 1); % FN = responsive channels included in the CAR
+hold on
+for nc = 1:length(noiseCoefs)
+ FN = squeeze(accuracy.FN(nc, :, :));
+ iqr = prctile(FN, [25, 75], 2);
+ med = median(FN, 2);
+ plot(nResponsives, med, '-o', 'Color', parCm(nc*3-2, :), 'LineWidth', 1);
+ errorbar(nResponsives, med, med - iqr(:, 1), iqr(:, 2) - med, 'Color', parCm(nc*3-2, :), 'LineWidth', 1, 'CapSize', 8);
+end
+hold off
+xlim([-1, 41]); ylim([-1, 50]);
+set(gca, 'xtick', 0:5:50, 'xticklabels', 0:5:50, 'box', 'off');
+ylabel('Misses');
+title('Responsive Channels Missed');
+
+subplot(2, 1, 2);
+hold on
+for nc = 1:length(noiseCoefs)
+ FP = squeeze(accuracy.FP(nc, :, :));
+ iqr = prctile(FP, [25, 75], 2);
+ med = median(FP, 2);
+ plot(nResponsives, med, '-o', 'Color', parCm(nc*3-2, :), 'LineWidth', 1);
+ errorbar(nResponsives, med, med - iqr(:, 1), iqr(:, 2) - med, 'Color', parCm(nc*3-2, :), 'LineWidth', 1, 'CapSize', 8);
+end
+hold off
+xlim([-1, 41]); ylim([-1, 50]);
+set(gca, 'xtick', 0:5:50, 'xticklabels', 0:5:50, 'box', 'off');
+xlabel('Number of responsive channels / 50');
+ylabel('Misses');
+title('Nonresponsive Channels Missed');
+saveas(gcf, fullfile('output', rootdir, sprintf('nchs%d_FNFP', nChs)), 'png');
+saveas(gcf, fullfile('output', rootdir, sprintf('nchs%d_FNFP', nChs)), 'svg');
+
+% Calculate mean SNR for each noise level (geo bc lognormal and arith)
+snrGeomean = nan(length(noiseCoefs), 1);
+snrGeoSD = nan(length(noiseCoefs), 1);
+snrMean = nan(length(noiseCoefs), 1);
+snrSD = nan(length(noiseCoefs), 1);
+for ii = 1:length(noiseCoefs)
+ snrsCurrCell = reshape(snrs(ii, :, :), 1, []);
+ snrsCurr = [];
+ for jj = 1:length(snrsCurrCell)
+ snrsCurr = [snrsCurr; snrsCurrCell{jj}(:)];
+ end
+ snrsCurr(isnan(snrsCurr)) = [];
+ snrGeomean(ii) = geomean(snrsCurr);
+ snrGeoSD(ii) = 10^std(log10(snrsCurr));
+ snrMean(ii) = mean(snrsCurr);
+ snrSD(ii) = std(snrsCurr);
+end
+
diff --git a/figS1_compareCARvCARLA.m b/figS3_compareCARvCARLA.m
similarity index 100%
rename from figS1_compareCARvCARLA.m
rename to figS3_compareCARvCARLA.m
diff --git a/figS3_simCCEPComponentsGlobalSig.m b/figS5_simCCEPComponentsGlobalSig.m
similarity index 100%
rename from figS3_simCCEPComponentsGlobalSig.m
rename to figS5_simCCEPComponentsGlobalSig.m
diff --git a/functions/CARLA.m b/functions/CARLA.m
index 6fd36b2..63d3cc9 100644
--- a/functions/CARLA.m
+++ b/functions/CARLA.m
@@ -158,7 +158,7 @@
nMin = ceil(0.1*nChs); % minimum number of channels to start at, for stability. Hardcoded at 10 %
zMMxTrs = mean(stats.zMinMean(1, :, :), 3); % mean ZMinMean across bootstrapped samples for each n
- ii = nMin;
+ ii = max(nMin, 2); % if there are 10 or fewer channels, require ii to start at 2 (minimum CAR size)
while ii <= nChs
if ii == nChs % at the end. no more comparisons needed
diff --git a/functions/CAR_subset.m b/functions/CAR_subset.m
new file mode 100644
index 0000000..2fb9c0b
--- /dev/null
+++ b/functions/CAR_subset.m
@@ -0,0 +1,205 @@
+function [Vout, chsUsed] = CAR_subset(tt, V, srate, badChs, grp, optsIn)
+%
+% This is copied from ccep_CARVariance.m from mnl_ieegBasics toolbox, renamed to be included as part of CARLA package in publication
+%
+% Applies adjusted common average referencing to iEEG using a subset of channels with lowest mean trial variance or cross-trial covariance on a response interval.
+% This function was created as a bridge between the old ccep_CAR functions in mnl_ieegBasics and the newer CARLA. An important change over previous methodology
+% is that this function uses the same channels as common average across all trials of the same group (e.g., stimulation site). The common average is still calculated
+% on a per-trial basis. This function was created with CCEP data in mind, but can be extrapolated to cognitive task data.
+%
+% Vout = CARVariance(tt, V); % minimum inputs, NOT RECOMMENDED
+% Vout = CARVariance(tt, V, srate); % recommended minimum inputs
+% [Vout, chsUsed] = CARVariance(tt, V, srate, badChs, grp); % most common usage
+% [Vout, chsUsed] = CARVariance(tt, V, srate, badChs, grp, opts); % full usage
+% tt = 1 x t num. Time points matching V, in seconds
+% V = n x t x k num, Signal data to be re-referenced, with n channels by t time points by k trials.
+% If there are any singular bad trials that shouldn't be considered in CAR, set them to nan in V. Similarly, channels that shouldn't be
+% considered only for a group (e.g. stimulated electrodes) should be set to nan for all trials corresponding to that group.
+% Percentile calculation ignores channels that are set to nan this way. Any nans at single time points will be propagated to the entire trial(s) containing the nan time point.
+% srate = num. Sampling frequency of data, used for notch filter if opts.notchfirst==true.
+% If not given and opts.notchfirst==true, this will be estimated from the tt input (not recommended because of possible imprecision).
+% badChs = (optional) 1 x m num, m < n. List of all bad channels to not consider for CAR (e.g., [36, 78, 100, 212:230]).
+% grp = (optional) 1 x k num, string array, or cell array(char). Grouping variable for trials. All trials in the same group must have identical values (numerical or character).
+% If not given, assumes the entire input is one group
+% optsIn = (optional) struct to adjust parameters, fields described below:
+% pctThresh = num, between 0 and 100. Percentile threshold on how many channels to use for common average. Default = 25
+% winResp = 1 x 2 num. [start, stop] of responsive time period to calculate variance and correlations on, in seconds matching tt. Default = [0.01, 0.3].
+% varType = char, 'var' or 'cov'. What metric to use when ranking channels. 'var' ranks channels based on (geometric) average variance across trials.
+% 'cov' ranks channels based on average cross-trial covariance. Default = 'cov'.
+% notchFirst = bool. Whether notch filters are applied at 60, 120, 180 Hz to reduce line noise on the data before ranking/optimizing by CARLA.
+% This only affects determining which channels to include, as the rereferenced output is calculated from the unfiltered input. Default = true.
+%
+% Returns:
+% Vout = n x t x k num, or n x t num. Re-referenced signal data.
+% chsUsed = q x 1 cell array num. For each of the q unique groups in , the corresponding cell contains the list of channel numbers used to construct the
+% common average for that group, in order of increasing .
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package.
+% Copyright (C) 2023 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 .
+%
+%
+% HH, DD Multimodal Neuroimaging Lab, Mayo Clinic, 2023
+%
+ %% Input checks, data formatting
+
+ nTrs = size(V, 3); % number of trials
+ assert(size(V, 2) == length(tt), 'Error: second dimension does not match tt. Data should be channels x timepoints x trials');
+
+ % propagate any single nan timepoints in input to entire trial(s) affected
+ nanmask = repmat(any(isnan(V), 2), 1, length(tt), 1);
+ V(nanmask) = nan;
+
+ % assume all are the same stim condition if none given
+ if nargin < 5 || isempty(grp), grp = ones(nTrs, 1); end
+ if nargin < 4, badChs = []; end % no bad channels
+
+ assert(isnumeric(badChs), 'Error: badChs must be given as a numeric list');
+
+ if isstring(grp), grp = cellstr(grp); end % convert to cell array of chars if given as string array
+ assert(length(grp) == nTrs, 'Error: grp input does not match trials in length');
+
+ % Default opts
+ opts.pctthresh = 25;
+ opts.winresp = [0.01, 0.3];
+ opts.vartype = 'cov';
+ opts.notchfirst = true;
+
+ % Configure opts if input is given
+ if exist('optsIn', 'var') && isa(optsIn, 'struct')
+ inFields = fieldnames(optsIn);
+ for ii = 1:length(inFields) % make case insensitive
+
+ optsIn.(lower(inFields{ii})) = optsIn.(inFields{ii});
+ fieldCurr = lower(inFields{ii});
+
+ assert(any(strcmp(fieldCurr, {'pctthresh', 'winresp', 'vartype', 'notchfirst'})), 'Error: "%s" is not a valid opts field', fieldCurr);
+
+ if strcmp(fieldCurr, 'alpha'), assert(optsIn.(fieldCurr) > 0 & optsIn.(fieldCurr) < 100, 'Error: pctthresh must be between 0 and 100'); end
+ if strcmp(fieldCurr, 'winresp'), assert(length(optsIn.(fieldCurr)) == 2 & all(optsIn.(fieldCurr) > 0) & all(optsIn.(fieldCurr) < tt(end)), ...
+ 'Error: winresp must be given as [t1, t2] where both are greater than 0 and before tt(end)'); end
+ if strcmp(fieldCurr, 'vartype'), assert(any(strcmp(optsIn.(fieldCurr), {'var', 'cov'})), 'Error: vartype has to be "var" or "cov"'); end
+ if strcmp(fieldCurr, 'notchfirst'), assert(islogical(optsIn.(fieldCurr)), 'Error: notchfirst must be logical (true or false)'); end
+
+ opts.(fieldCurr) = optsIn.(fieldCurr);
+ end
+ end
+
+ % Estimate sampling frequency from time points if necessary (not recommended)
+ if nargin < 3 || isempty(srate) && opts.notchfirst
+ srate = (length(tt) - 1) / (max(tt) - min(tt));
+ warning('sampling frequency for notch filter was ESTIMATED to be %0.02f', srate);
+ end
+
+ % Remove bad channels
+ VgoodChs = V;
+ VgoodChs(badChs, :, :) = [];
+ goodChs = setdiff(1:size(V, 1), badChs); % keep track of indices of good chs to put back later
+
+ nChs = size(VgoodChs, 1);
+ fprintf('%d good channels out of %d total channels considered for CAR\n', nChs, size(V, 1));
+
+ fprintf('Opts: %d percent of channels, [%0.03f - %0.03f] ms window, using %s, notchFirst=%s\n', ...
+ opts.pctthresh, opts.winresp(1), opts.winresp(2), opts.vartype, string(opts.notchfirst));
+
+ %% Apply notch filter to reduce line noise (if requested)
+
+ Vclean = VgoodChs;
+ if opts.notchfirst
+ fprintf('Applying notch filters: ');
+ for ff = 60:60:180
+ fprintf('%d Hz, ', ff);
+ dNotch = designfilt('bandstopiir', 'FilterOrder', 4, ... % matches filter in ieeg_notch
+ 'DesignMethod', 'butter', ...
+ 'HalfPowerFrequency1', ff-2, ...
+ 'HalfPowerFrequency2', ff+2, ...
+ 'SampleRate', srate);
+ for ii = 1:nTrs % filter per trial
+ VcleanTr = Vclean(:, :, ii); % all channels at current trial
+ VcleanTr(~any(isnan(VcleanTr), 2), :) = filtfilt(dNotch, VcleanTr(~any(isnan(VcleanTr), 2), :)')'; % filter channels with non-nan time points
+ Vclean(:, :, ii) = VcleanTr;
+ end
+ end
+ fprintf('\n');
+ end
+
+ %% Loop through trial conditions to determine which channels to use for each condition
+
+ % extract data on response interval to perform analysis on
+ Vseg = Vclean(:, tt >= opts.winresp(1) & tt <= opts.winresp(end), :);
+
+ grpTypes = unique(grp);
+ chsUsed = cell(length(grpTypes), 1);
+
+ fprintf('Conditions: ');
+ for ii = 1:length(grpTypes)
+
+ if isnumeric(grp) % handle numeric and string (not char array) inputs
+ VsegCurr = Vseg(:, :, grp == grpTypes(ii));
+ condStr = num2str(grpTypes(ii));
+ elseif iscell(grp)
+ VsegCurr = Vseg(:, :, strcmp(grp, grpTypes{ii}));
+ condStr = grpTypes{ii};
+ end
+ fprintf('%s (%d), ', condStr, size(VsegCurr, 3));
+ if ~mod(ii, 5), fprintf('\n'); end
+
+ if strcmp(opts.vartype, 'var') % Geomean variance across trials
+ vars = geomean(var(VsegCurr, 0, 2), 3, 'omitnan'); % omitnan ignores bad trials set to nan
+ elseif strcmp(opts.vartype, 'cov') % Mean covariance
+ if size(VsegCurr, 3) < 6, warning('Mean covariance may be unreliable for few trials due to spurious anticorrelations between trials. Use ''var'' instead'); end
+ vars = nan(nChs, 1);
+ for jj = 1:nChs
+ covCurr = cov(squeeze(VsegCurr(jj, :, :)));
+ vars(jj) = mean(covCurr(logical(triu(ones(size(covCurr)), 1))), 'all', 'omitnan'); % average the covariances above the diagonal, ignoring bad trials (nan-coded)
+ end
+ end
+ nanChs = find(isnan(vars));
+ if ~isempty(nanChs)
+ warning('All trials are nan at channels %s, not considered for CAR at current condition', join(num2str(goodChs(nanChs))));
+ end
+
+ % variance threshold for current condition
+ thresh = prctile(vars, opts.pctthresh);
+
+ % sort good channels in increasing order of variance so easier to troubleshoot
+ [vars, ord] = sort(vars);
+ goodChsSorted = goodChs(ord);
+ chsUsed{ii} = goodChsSorted(vars < thresh); % which channels fall under threshold, by original ch number, sorted in order of increasing variance
+
+ end
+ fprintf('\n');
+
+ %% Perform CAR on each condition
+
+ Vout = V;
+
+ for ii = 1:length(grpTypes)
+
+ if isnumeric(grp) % handle numeric and string (not char array) inputs
+ inds = grp == grpTypes(ii);
+ elseif iscell(grp)
+ inds = strcmp(grp, grpTypes{ii});
+ end
+
+ Vout(:, :, inds) = Vout(:, :, inds) - mean(Vout(chsUsed{ii}, :, inds), 'omitnan'); % omitnan ignores bad trials when calculating common average
+ end
+
+end
\ No newline at end of file
diff --git a/main.asv b/main.asv
new file mode 100644
index 0000000..773e252
--- /dev/null
+++ b/main.asv
@@ -0,0 +1,187 @@
+%% Use each section of this main script to output all results and figures for the manuscript.
+%
+% 2023/09/27
+%
+% If this code is used in a publication, please cite the manuscript:
+% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+% GA Worrell, KJ Miller, and D Hermes.
+%
+% CARLA manuscript package: main script
+% Copyright (C) 2023 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 .
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Unless specified, each section of code below can be run directly from this script to generate all outputs.
+% The scripts listed may also be opened so that individual sections within can be separately examined and run.
+% Note that in output filenames, "zminmean" refers to zeta in the manuscript
+%
+% Subject CCEP and anatomical data will be available in BIDS format upon manuscript publication.
+% In the meantime, all results on simulated data can be reproduced with the appropriate section below.
+%
+%% Configure paths and dependencies
+
+% Make sure you are currently working in the root manuscript directory (CARLA)
+addpath('functions');
+
+% Data should be downloaded from OpenNeuro and copied to a folder here (./data)
+dataPath = 'data';
+
+% Add paths to other dependencies. Change the dummy paths listed here to where they are located on your system.
+addpath(genpath('path/to/mnl_ieegBasics'));
+addpath(genpath('path/to/matmef'));
+addpath('path/to/mnl_seegview');
+addpath(genpath('path/to/vistasoft'));
+addpath('path/to/spm12'); addpath(fullfile('path/to/spm12', 'config'));
+addpath(genpath('path/to/freesurfer/7.1.1/matlab'));
+addpath(genpath('path/to/gifti'));
+
+% To save figures as vectorized objects
+set(0, 'DefaultFigureRenderer', 'painters');
+
+%% i) Generate T1 MRI slices with stimulation site electrodes overlaid (Figures 4A, 6C)
+
+% This was done using sliceGUI.m in the mnl_seegview repository:
+% https://github.com/MultimodalNeuroimagingLab/mnl_seegview
+
+% Define paths to defaced T1 MRI and matching electrodes table
+niiPath = fullfile(dataPath, 'derivatives', 'freesurfer', 'sub-1', 'sub-1_ses-mri01_T1w_acpc_deFaced.nii');
+electrodes = fullfile(dataPath, 'sub-1', 'ses-ieeg01', 'ieeg', 'sub-1_ses-ieeg01_electrodes.tsv');
+
+% Run sliceGUI, and manually interact with it
+sliceGUI(niiPath, electrodes);
+
+%% Table 1: Calculate the number of measurement electrodes (channels) used in each subject
+
+% outputs are directly printed to command window
+tab1_dispNumElectrodes
+
+%% Figure 1: Save simulated CCEP channels and their theoretical components
+
+% outputs saved to .\output\simComponents
+fig1_simCCEPComponents
+
+%% Figure 2: Methods illustration of how CARLA works, on single-trial simulated CCEP data
+
+% outputs saved to .\output\simulationTrial
+fig2_CARLAonSingleTrial
+
+%% Figure 3: Save example of enveloped sinusoids used to simulate evoked potentials
+
+% outputs saved to .\output\EPConstruct
+fig3_EPConstruct
+
+%% Figure 4: Simulate 50 channels with responsiveness varying from 0 to 50, with 30 repetitions each, and apply CARLA
+
+% First run this to generate outputs for each number of responsive channels
+% outputs saved to .\output\simLoop
+Aglobal = 0; % no global signal desired (0 amplitude)
+simCCEPLoop;
+
+% Then, keeping Aglobal in the workspace, run this to make summary figures (4C, D), as well as identify median candidates for 4A, B
+% .\output\simLoop will be renamed .\output\simLoopNoGlob, where outputs are saved.
+summarizeSimCCEPLoop
+
+%% Figure 5: Apply CARLA to real data (sub 1, stim site 1)
+
+% Configure the subject, data run, and stim site
+sub = 'sub-1';
+run = 2;
+site = 'RMO8-RMO9';
+
+% which channels to plot in panel F
+chs2Plot = [100, 203];
+
+% outputs saved to .\output\realCCEPs\. These correspond to figures 5B-F. To generate 5A (stim site overlaid on T1 MRI), see section i), above.
+applyCARLARealCCEPs
+
+%% Figure S1: Simulate datasets with fewer total channels (25, 20, 15, 10) and calculate CARLA accuracy
+
+figS1_simCCEPTotalChs
+
+%% Figure S2: Simulate datasets with lower signal-to-noise ratio (SNR) and calculate CARLA accuracy
+
+brownCoefRange
+
+%% Figure S3: Average channel responses for sub 1, stim site 1 after applying CAR vs. CARLA
+
+% This supplemental figure is a direct addendum to Figure 5, showing the mean signal for the first 50 channels after rereferencing by CARLA or CAR
+% Run after the previous section, keeping all workspace variables
+figS3_compareCARvCARLA
+
+%% Figure S4: Apply CARLA to real data at another 4 stim sites in subjects 1-4
+
+% config parameters for each subject shown in S1
+config = struct('sub', {'sub-1', 'sub-2', 'sub-3', 'sub-4'}, ...
+ 'run', {1, 12, 3, 2}, ...
+ 'site', {'ROC11-ROC12', 'RZ4-RZ5', 'LA2-LA3', 'RO3-RO4'});
+chs2Plot = 1; % dummy value, no individual channels are plotted in Figure S1
+
+for ii = 1:length(config)
+ sub = config(ii).sub;
+ run = config(ii).run;
+ site = config(ii).site;
+
+ % outputs saved to .\output\realCCEPs\.
+ % Omission of the middle section of sorted channels (subjects 1-3) in Figure S1 was done manually in illustrator.
+ applyCARLARealCCEPs
+end
+
+%% Figure S5: Save simulated CCEP channels, this time with a global signal, and their theoretical components
+
+% outputs saved to .\output\simComponentsGlobal
+figS5_simCCEPComponentsGlobalSig
+
+%% Figure 6A, B Simulate 50 channels with responsiveness from 0 to 50, like Figure 4, but now with a global signal added
+
+% First run this to generate outputs for each number of responsive channels
+% outputs saved to .\output\simLoop
+Aglobal = 25; % global signal with mean amplitude 25
+simCCEPLoop;
+
+% Then, keeping Aglobal in the workspace, run this to make summary figure (6B), and to identify a trial candidate for 6A
+% .\output\simLoop will be renamed .\output\simLoopWithGlob, where outputs are saved.
+summarizeSimCCEPLoop
+
+% To generate figure 6C, see section i) on using sliceGUI, above.
+
+%% Figure 6D-F: Apply CARLA to real data containing global signal (sub 1, stim site 2)
+
+% Configure the subject, data run, and stim site
+sub = 'sub-1';
+run = 3;
+site = 'RK5-RK6';
+
+% which channels to plot in panel F
+chs2Plot = [23, 51, 205];
+
+% outputs saved to .\output\realCCEPs\
+applyCARLARealCCEPs
+
+%% Figure 7: Loop through stimulation sites in each subject and compute optimal CARLA cutoff
+% This also calculates the cross-channel R-squared values for all possible adjusted common average sizes, to be summarized in Figure 8
+
+% First run this script. Outputs are saved to .\output\realCCEPsLoop\.
+applyCARLARealCCEPsLoop
+
+% Then run this script to generate outputs for Figure 7
+fig7_summarizeCARLARealCCEPs
+
+%% Figure 8: Calculate cross-channel R-squared values for no CAR vs different types of fixed CARs vs CARLA
+
+% outputs saved to .\output\realCCEPsLoop\.
+fig8_summarizeInterChCorr
+
diff --git a/main.m b/main.m
index 5253c9a..b157cb1 100644
--- a/main.m
+++ b/main.m
@@ -86,7 +86,7 @@
%% Figure 4: Simulate 50 channels with responsiveness varying from 0 to 50, with 30 repetitions each, and apply CARLA
-% First run this to generate outputs for each number of responsive channels
+% First, run this to generate outputs for each number of responsive channels
% outputs saved to .\output\simLoop
Aglobal = 0; % no global signal desired (0 amplitude)
simCCEPLoop;
@@ -95,6 +95,36 @@
% .\output\simLoop will be renamed .\output\simLoopNoGlob, where outputs are saved.
summarizeSimCCEPLoop
+%% Figure S1: Simulate datasets like Fig 4, but with fewer total channels (25, 20, 15, 10) and calculate CARLA accuracy
+
+% First, run this to generate outputs for each total channel size and fraction of responsive channels
+% outputs are saved to .\output\simLoopNChs
+figS1a_simCCEPTotalChs
+
+% Then run this to make summary figure (S1A). S1B was constructed with examples most similar to the median of this summary output
+figS1b_summarizeSimCCEPTotalChs
+
+%% Figure S2: Simulate datasets with lower signal-to-noise ratio (SNR) and calculate CARLA accuracy
+
+% Brown noise coefficients to save outputs for. Split up into chunks for more efficient parallelization
+% SNR is varied by changing the Brown Noise coefficient. The population mean SNR is estimated at the end with figS2b_summarizeSimCCEPSNRs, below.
+brownCoefRanges = [0.4, 0.5;
+ 0.6, 0.7;
+ 0.8, 0.9;
+ 1.0, 1.0];
+
+% First, run this to generate outputs for each SNR level and number of responsive channels. This also saves the plots for Figure S2A
+% outputs are saved to .\output\simLoopSNR
+% this loop can be parallelized if desired
+for ii = 1:4
+ brownCoefRange = brownCoefRanges(ii, :);
+ figS2a_simCCEPSNRs
+end
+
+% Then run this to make summary figure, S2B.
+% This also calculates the geomean (sigma) SNR for each Brown Noise coefficient condition (variable snrGeomean in workspace)
+figS2b_summarizeSimCCEPSNRs
+
%% Figure 5: Apply CARLA to real data (sub 1, stim site 1)
% Configure the subject, data run, and stim site
@@ -108,13 +138,13 @@
% outputs saved to .\output\realCCEPs\. These correspond to figures 5B-F. To generate 5A (stim site overlaid on T1 MRI), see section i), above.
applyCARLARealCCEPs
-%% Figure S1: Average channel responses for sub 1, stim site 1 after applying CAR vs. CARLA
+%% Figure S3: Average channel responses for sub 1, stim site 1 after applying CAR vs. CARLA
% This supplemental figure is a direct addendum to Figure 5, showing the mean signal for the first 50 channels after rereferencing by CARLA or CAR
% Run after the previous section, keeping all workspace variables
-figS1_compareCARvCARLA
+figS3_compareCARvCARLA
-%% Figure S2: Apply CARLA to real data at another 4 stim sites in subjects 1-4
+%% Figure S4: Apply CARLA to real data at another 4 stim sites in subjects 1-4
% config parameters for each subject shown in S1
config = struct('sub', {'sub-1', 'sub-2', 'sub-3', 'sub-4'}, ...
@@ -132,10 +162,10 @@
applyCARLARealCCEPs
end
-%% Figure S3: Save simulated CCEP channels, this time with a global signal, and their theoretical components
+%% Figure S5: Save simulated CCEP channels, this time with a global signal, and their theoretical components
% outputs saved to .\output\simComponentsGlobal
-figS3_simCCEPComponentsGlobalSig
+figS5_simCCEPComponentsGlobalSig
%% Figure 6A, B Simulate 50 channels with responsiveness from 0 to 50, like Figure 4, but now with a global signal added
diff --git a/simCCEPLoop.m b/simCCEPLoop.m
index adffd52..d1322ff 100644
--- a/simCCEPLoop.m
+++ b/simCCEPLoop.m
@@ -151,6 +151,7 @@
close all;
% Accuracy values. positive means responsive/excluded from CAR
+ % We keep these variables as named here, but note that FN and FP are now renamed RCM and NCM in the manuscript.
TP = sum(find(respBool) > nCAR); % responsive channels successfully excluded from CAR (above the cutoff)
TN = sum(find(~respBool) <= nCAR); % NR channels successfully below or at cutoff
FN = sum(find(respBool) <= nCAR); % responsive channels incorrectly included in CAR. *This matters most