OSCR

Alpha frequency shapes perceptual sensitivity by modulating optimal phase likelihood.

Code ↔ Paper

4 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 4 matches · 1 of them tie a paragraph to a whole file, not to given lines: a weak match, whose lines are not tinted
  1. [1] § Methods › Extract instantaneous alpha frequency and phase ↔ SCRIPT/restingIAF_natcomm.m, the whole file · a weak match · score 0.64 · alpha band, spectral, derivative, spectopo, filtered, noise
  2. [2] § Methods › Extract instantaneous alpha frequency and phase ↔ SCRIPT/Script_NC.m, lines 1004–1130 · score 0.64 · single trial IAF, instantaneous frequency, IAF accuracy, EEGLAB, FFT, pre
  3. [3] § Methods › Assessing the specificity of the relationship between IAF and perceptual accuracy relative to alpha power ↔ SCRIPT/Script_NC.m, lines 780–913 · score 0.62 · alpha power, trial fluctuation, perceptual accuracy, amplitude, EEG, incorrect
  4. [4] § Methods › IAF single-trial regression: accuracy ↔ SCRIPT/Script_NC.m, lines 1004–1130 · score 0.54 · single trial IAF, instantaneous frequency, row, FFT, pre, stimulus

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

MATLAB · 1,130 lines · 44 KB · no license · 3 matches

  1. %% ------------------------- EXTRACT PEAK & ELECTRODES ---------------------
  2. clc; clear; close all
  3. eegfolder = ''%path
  4. cd(EEGFOLDER);
  5. dEEG = dir(EEGFOLDER);
  6. EEG_files = struct2cell(dEEG)';
  7. EEG_files = EEG_files(3:end,1);
  8. start_from = 2;
  9. subji = start_from:2:length(EEG_files);
  10. dataEEG = cell(1, numel(subji));
  11. cd(EEGFOLDER);
  12. for p = 1:numel(subji)
  13. idxEEG = subji(p);
  14. filename = char(EEG_files{idxEEG});
  15. EEG = pop_loadset('filename', filename);
  16. EEG = eeg_checkset(EEG);
  17. dataEEG{p} = EEG;
  18. end
  19. electrodes_peak = cell(1, numel(subji));
  20. peak = nan(1, numel(subji));
  21. for soggetto = 1:numel(subji)
  22. EEG = dataEEG{soggetto};
  23. EEG.data = diff(EEG.data, 1, 2);
  24. data = double(EEG.data);
  25. chosedelectrodes = [18,49,50,17,48];
  26. time2anal = [-800 -100];
  27. time2idx = dsearchn(EEG.times', time2anal');
  28. datawin = data(chosedelectrodes, time2idx(1):time2idx(2), :);
  29. [a,b] = spectopo(datawin, 0, EEG.srate, 'nfft', EEG.srate*100, 'plot', 'off', 'verbose', 'off');
  30. alphalim = [7 13];
  31. alphaidx = dsearchn(b, alphalim');
  32. alphapow = mean(a(:, alphaidx(1):alphaidx(2)), 2);
  33. [~, which_elec] = sort(alphapow, 'descend');
  34. electrodes_peak{soggetto} = EEG.chanlocs(chosedelectrodes(which_elec(1))).labels;
  35. data1 = squeeze(EEG.data(chosedelectrodes(which_elec(1)), :, :));
  36. data1 = reshape(data1, 1, size(data1,1), size(data1,2));
  37. cmin=1; fRange=[1 40]; w=[7 13]; Fw=11; k=5; Fs=EEG.srate;
  38. data_colcoran = reshape(data1(:, time2idx(1):time2idx(2), :), 1, []);
  39. [~, colcalpha] = restingIAF(data_colcoran, 1, cmin, fRange, Fs, w, Fw, k);
  40. peak(soggetto) = colcalpha.peaks;
  41. for elec = 2:length(chosedelectrodes)
  42. if isnan(peak(soggetto))
  43. electrodes_peak{soggetto} = EEG.chanlocs(chosedelectrodes(which_elec(elec))).labels;
  44. data_alt = squeeze(EEG.data(chosedelectrodes(which_elec(elec)), time2idx(1):time2idx(2), :));
  45. data_alt = reshape(data_alt, 1, size(data_alt,1), size(data_alt,2));
  46. data_colcoran = reshape(data_alt, 1, []);
  47. [~, colcalpha] = restingIAF(data_colcoran, 1, cmin, fRange, Fs, w, Fw, k);
  48. peak(soggetto) = colcalpha.peaks;
  49. end
  50. end
  51. end
  52. %% --------------------- EXTRACT INST_IAF --------------------
  53. clc; clear; close all
  54. load('AllSubjects.mat')
  55. eegfolder = ''%path
  56. cd(EEGFOLDER);
  57. dEEG = dir(EEGFOLDER);
  58. EEG_files = struct2cell(dEEG)';
  59. EEG_files = EEG_files(3:end,1);
  60. start_from = 2;
  61. subji = start_from:2:length(EEG_files);
  62. dataEEG = cell(1, numel(subji));
  63. cd(EEGFOLDER);
  64. for p = 1:numel(subji)
  65. idxEEG = subji(p);
  66. filename = char(EEG_files{idxEEG});
  67. EEG = pop_loadset('filename', filename);
  68. EEG = eeg_checkset(EEG);
  69. dataEEG{p} = EEG;
  70. end
  71. electrodes_to_analyze_all = readcell('electrodes_to_analyze.xlsx');
  72. peak = readmatrix('peak_to_analyze.xlsx');
  73. for subject = 1:numel(databehx)
  74. EEG = dataEEG{subject};
  75. EEG.data = diff(EEG.data, 1, 2);
  76. electrodes_to_analyze = find(strcmpi({EEG.chanlocs.labels}, electrodes_to_analyze_all{subject}));
  77. alpha_lim = [peak(subject)-2, peak(subject)+2];
  78. [inst_freq, inst_ph] = instantaneous_NatComm(EEG.data(electrodes_to_analyze,:,:), EEG.srate, alpha_lim);
  79. inst_freq = squeeze(inst_freq);
  80. inst_ph = squeeze(inst_ph);
  81. end
  82. %% Pre-stimulus IAF is linked to variations in perceptual sensitivity: BIN analysis
  83. clc,clear
  84. load('IAF_allSubjects.mat');
  85. nS = numel(datax);
  86. d = zeros(nS,2);
  87. c = zeros(nS,2);
  88. accuracy = zeros(nS,2);
  89. iaf_mean_terciles = zeros(nS,2);
  90. for s = 1:nS
  91. data = datax{s};
  92. data(:,4) = mean(data(:,4:end), 2);
  93. data(:,5:end) = [];
  94. [~, h] = sort(data(:,4));
  95. data = data(h,:);
  96. N = size(data,1);
  97. iLow = 1:ceil(N/3);
  98. iHigh = (ceil(2*N/3)+1):N;
  99. % ========== FIRST TERCILE ===============
  100. sub = data(iLow,:);
  101. iaf_mean_terciles(s,1) = mean(sub(:,4));
  102. stim_pres = sub(sub(:,1)==1,:);
  103. stim_abs = sub(sub(:,1)==0,:);
  104. nPres = size(stim_pres,1);
  105. nAbs = size(stim_abs,1);
  106. hit = sum(stim_pres(:,1)==stim_pres(:,2)) / nPres;
  107. fa = sum(stim_abs(:,1)~=stim_abs(:,2)) / nAbs;
  108. if hit == 0, hit = 0.5/nPres; elseif hit == 1, hit = (nPres-0.5)/nPres; end
  109. if fa == 0, fa = 0.5/nAbs; elseif fa == 1, fa = (nAbs-0.5)/nAbs; end
  110. d(s,1) = norminv(hit) - norminv(fa);
  111. c(s,1) = -(norminv(hit) + norminv(fa)) / 2;
  112. accuracy(s,1) = mean(sub(:,1) == sub(:,2));
  113. % ========== THIRD TERCILE ==============
  114. sub = data(iHigh,:);
  115. iaf_mean_terciles(s,2) = mean(sub(:,4));
  116. stim_pres = sub(sub(:,1)==1,:);
  117. stim_abs = sub(sub(:,1)==0,:);
  118. nPres = size(stim_pres,1);
  119. nAbs = size(stim_abs,1);
  120. hit = sum(stim_pres(:,1)==stim_pres(:,2)) / nPres;
  121. fa = sum(stim_abs(:,1)~=stim_abs(:,2)) / nAbs;
  122. if hit == 0, hit = 0.5/nPres; elseif hit == 1, hit = (nPres-0.5)/nPres; end
  123. if fa == 0, fa = 0.5/nAbs; elseif fa == 1, fa = (nAbs-0.5)/nAbs; end
  124. d(s,2) = norminv(hit) - norminv(fa);
  125. c(s,2) = -(norminv(hit) + norminv(fa)) / 2;
  126. accuracy(s,2) = mean(sub(:,1) == sub(:,2));
  127. end
  128. %% Time-resolved binning analysis revealed that the IAF effect on sensitivity
  129. %% is broadly extended over the pre-stimulus period: Bin analysis timeresolved
  130. clc,clear
  131. load('IAF_allSubjects.mat');
  132. nperms = 1000;
  133. rng('Shuffle')
  134. ntpns = size(datax{1},2) - 3;
  135. nS = numel(datax);
  136. iaf_mean_timeresolved = zeros(nS,ntpns,2);
  137. d = zeros(nS,ntpns,2);
  138. c = zeros(nS,ntpns,2);
  139. accuracy = zeros(nS,ntpns,2);
  140. for timepointi = 1:ntpns
  141. for participant = 1:nS
  142. data = datax{participant};
  143. data(:,4) = data(:, 3+timepointi);
  144. data(:,5:end)= [];
  145. [~, h] = sort(data(:,4));
  146. data_real = data(h,:);
  147. N = size(data_real,1);
  148. iL = 1:ceil(N/3);
  149. data_low_alpha = data_real(iL,:);
  150. stim_pres = data_low_alpha(data_low_alpha(:,1)==1,:);
  151. hit = (sum(stim_pres(:,1) == stim_pres(:,2)))/(length(stim_pres));
  152. if hit == 1, hit = (length(stim_pres) - 0.5)/length(stim_pres); end
  153. if hit == 0, hit = 0.5/length(stim_pres); end
  154. stim_abs = data_low_alpha(data_low_alpha(:,1)==0,:);
  155. fa = (sum(stim_abs(:,1) ~= stim_abs(:,2)))/(length(stim_abs));
  156. if fa == 0, fa = 0.5/length(stim_abs); end
  157. d(participant,timepointi,1) = norminv(hit) - norminv(fa);
  158. c(participant,timepointi,1) = -(norminv(hit) + norminv(fa))/2;
  159. accuracy(participant,timepointi,1) = sum(data_low_alpha(:,1)==data_low_alpha(:,2))/length(data_low_alpha);
  160. iH = (ceil(2*N/3)+1):N;
  161. data_hig_alpha = data_real(iH,:);
  162. stim_pres = data_hig_alpha(data_hig_alpha(:,1)==1,:);
  163. hit = (sum(stim_pres(:,1) == stim_pres(:,2)))/(length(stim_pres));
  164. if hit == 1, hit = (length(stim_pres) - 0.5)/length(stim_pres); end
  165. if hit == 0, hit = 0.5/length(stim_pres); end
  166. stim_abs = data_hig_alpha(data_hig_alpha(:,1)==0,:);
  167. fa = (sum(stim_abs(:,1) ~= stim_abs(:,2)))/(length(stim_abs));
  168. if fa == 0, fa = 0.5/length(stim_abs); end
  169. d(participant,timepointi,2) = norminv(hit) - norminv(fa);
  170. c(participant,timepointi,2) = -(norminv(hit) + norminv(fa))/2;
  171. accuracy(participant,timepointi,2) = sum(data_hig_alpha(:,1)==data_hig_alpha(:,2))/length(data_hig_alpha);
  172. end
  173. end
  174. d = d; % d, accuracy, c
  175. d_t_t_lowalpha = mean(d(:,:,1));
  176. d_t_t_highalpha = mean(d(:,:,2));
  177. load('timeperiod_iaf.mat')
  178. x = timeperiod;
  179. sig_real = zeros(ntpns,1);
  180. t_real = zeros(ntpns,1);
  181. for timi = 1:ntpns
  182. lowalpha = d(:,timi,1);
  183. highalpha = d(:,timi,2);
  184. [sig_real(timi),~,~,e] = ttest(lowalpha,highalpha);
  185. t_real(timi) = e.tstat;
  186. end
  187. t_real(~sig_real) = 0;
  188. islands = bwconncomp(t_real);
  189. clustsizes = islands.PixelIdxList;
  190. sum_t = zeros(length(clustsizes),1);
  191. for clusti = 1:length(clustsizes)
  192. ii = clustsizes{clusti};
  193. sum_t(clusti) = sum(abs(t_real(ii(1):ii(end))));
  194. end
  195. sig_fake = zeros(ntpns,nperms);
  196. t_fake = zeros(ntpns,nperms);
  197. max_cluster_sizes = zeros(nperms,1);
  198. for perm = 1:nperms
  199. lowalpha = d(:,:,1);
  200. highalpha = d(:,:,2);
  201. randorder = randperm(nS);
  202. randorder_1 = rem(randorder,2);
  203. lowalpha_fake = zeros(size(lowalpha));
  204. highalpha_fake = zeros(size(highalpha));
  205. for part = 1:nS
  206. if randorder_1(part)==1
  207. lowalpha_fake(part,:) = lowalpha(part,:);
  208. highalpha_fake(part,:) = highalpha(part,:);
  209. else
  210. lowalpha_fake(part,:) = highalpha(part,:);
  211. highalpha_fake(part,:) = lowalpha(part,:);
  212. end
  213. end
  214. for timi = 1:ntpns
  215. [sig_fake(timi,perm),~,~,e] = ttest(lowalpha_fake(:,timi),highalpha_fake(:,timi));
  216. t_fake(timi,perm) = e.tstat;
  217. end
  218. t_fake(~sig_fake(:,perm),perm) = 0;
  219. islands = bwconncomp(t_fake(:,perm));
  220. tempclustsizes = islands.PixelIdxList;
  221. sum_t_fake = zeros(max([length(tempclustsizes),1]),1);
  222. if ~isempty(tempclustsizes)
  223. for clusti = 1:length(tempclustsizes)
  224. ii = tempclustsizes{clusti};
  225. sum_t_fake(clusti) = sum(abs(t_fake(ii(1):ii(end),perm)));
  226. end
  227. end
  228. max_cluster_sizes(perm) = max(sum_t_fake);
  229. end
  230. p_value_cluster = zeros(length(clustsizes),1);
  231. S = sort(max_cluster_sizes);
  232. for clusti = 1:length(clustsizes)
  233. p_value_cluster(clusti) = (nperms - dsearchn(S,sum_t(clusti)')) / nperms;
  234. end
  235. %% IAF is higher in correct vs. incorrect trials.
  236. clc,clear
  237. load('IAF_allSubjects.mat');
  238. nS = length(datax);
  239. iaf = zeros(nS,2);
  240. for participant = 1:nS
  241. data = datax{participant};
  242. data(:,4) = mean(data(:,4:end),2);
  243. data(:,5:end) = [];
  244. accuracy = (data(:,1) == data(:,2));
  245. iaf_acc = data(accuracy ,4);
  246. iaf_inc = data(~accuracy,4);
  247. iaf(participant,1) = mean(iaf_acc);
  248. iaf(participant,2) = mean(iaf_inc);
  249. end
  250. %% Trial-by-trial fluctuations in IAF predict the accuracy of perceptual report.
  251. clc,clear
  252. load('IAF_allSubjects.mat');
  253. rng shuffle
  254. nperms = 2000;
  255. nS = length(datax);
  256. b_perms = zeros(nperms,1);
  257. b_z = zeros(nS,1);
  258. b = zeros(nS,1);
  259. for participant = 1:nS
  260. participant
  261. data = datax{participant};
  262. data(:,4) = mean(data(:,4:end),2);
  263. data(:,5:end) = [];
  264. data(:,4) = zscore(data(:,4));
  265. accuracy = (data(:,1) == data(:,2));
  266. X = [ones(length(data),1), data(:,4)];
  267. t = (X'*X)\(X'*accuracy);
  268. b(participant) = t(2);
  269. datat = data;
  270. for permi = 1:nperms
  271. datat(:,4) = Shuffle(data(:,4));
  272. Xp = [ones(length(data),1), datat(:,4)];
  273. tp = (Xp'*Xp)\(Xp'*accuracy);
  274. b_perms(permi) = tp(2);
  275. end
  276. b_z(participant) = (b(participant) - mean(b_perms)) / std(b_perms);
  277. end
  278. %% Alpha phase effects on perceptual sensitivity are moderated by IAF.
  279. clc,clear
  280. rng('Shuffle')
  281. load('timeperiod_phase.mat');
  282. NT = numel(timeperiod);
  283. which_group = 'below'; % 'below' | 'above'
  284. switch which_group
  285. case 'below'
  286. load('PHASE_below.mat');
  287. datax = data_phase_below;
  288. group = 1:numel(datax);
  290. case 'above'
  291. load('PHASE_above.mat');
  292. datax = data_phase_above;
  293. group = 1:numel(datax);
  294. end
  295. acc_phase = zeros(length(group), NT, 2);
  296. for participant = 1:length(group)
  297. X = datax{group(participant) };
  298. acc = (X(:,1) == X(:,2));
  299. for t = 1:NT
  300. phase_degree = rad2deg( X(:, t+3) );
  301. numBins = 2;
  302. binEdges = linspace(-180, 180, numBins+1);
  303. binIndices = zeros(size(phase_degree));
  304. for i = 1:numBins
  305. binIndices( phase_degree >= binEdges(i) & phase_degree <= binEdges(i+1) ) = i;
  306. end
  307. peak = (binIndices==2);
  308. trough = (binIndices==1);
  309. acc_phase(participant,t,1) = sum(acc(peak)) / sum(peak);
  310. acc_phase(participant,t,2) = sum(acc(trough)) / sum(trough);
  311. end
  312. end
  313. sign = zeros(NT,1);
  314. tval = zeros(NT,1);
  315. for t = 1:NT
  316. [sign(t),~,~,st] = ttest(acc_phase(:,t,1), acc_phase(:,t,2));
  317. tval(t) = st.tstat;
  318. end
  319. tval(sign<1) = 0;
  320. CC = bwconncomp(tval);
  321. clustsizes = CC.PixelIdxList;
  322. sum_t = zeros(numel(clustsizes),1);
  323. for c = 1:numel(clustsizes)
  324. idx = clustsizes{c};
  325. sum_t(c) = sum(abs(tval(idx(1):idx(end))));
  326. end
  327. nperms = 1000;
  328. sig_fake = zeros(NT,nperms);
  329. t_fake = zeros(NT,nperms);
  330. max_cluster_sizes = zeros(nperms,1);
  332. for perm = 1:nperms
  333. lowalpha = acc_phase(:,:,1);
  334. highalpha = acc_phase(:,:,2);
  335. randorder = randperm(size(acc_phase,1));
  336. randorder_1 = rem(randorder,2);
  337. lowalpha_fake = zeros(size(lowalpha));
  338. highalpha_fake = zeros(size(highalpha));
  339. for s = 1:size(acc_phase,1)
  340. if randorder_1(s)==1
  341. lowalpha_fake(s,:) = lowalpha(s,:);
  342. highalpha_fake(s,:) = highalpha(s,:);
  343. else
  344. lowalpha_fake(s,:) = highalpha(s,:);
  345. highalpha_fake(s,:) = lowalpha(s,:);
  346. end
  347. end
  348. for t = 1:NT
  349. [sig_fake(t,perm),~,~,e] = ttest(lowalpha_fake(:,t), highalpha_fake(:,t));
  350. t_fake(t,perm) = e.tstat;
  351. end
  352. t_fake(~sig_fake(:,perm),perm) = 0;
  353. CCp = bwconncomp(t_fake(:,perm));
  354. tmpMass = zeros(max([CCp.NumObjects,1]),1);
  355. if CCp.NumObjects>0
  356. for c = 1:CCp.NumObjects
  357. idx = CCp.PixelIdxList{c};
  358. tmpMass(c) = sum(abs(t_fake(idx(1):idx(end),perm)));
  359. end
  360. else
  361. tmpMass(1) = 0;
  362. end
  363. max_cluster_sizes(perm) = max(tmpMass);
  364. end
  365. p_value_cluster = zeros(numel(clustsizes),1);
  366. mx = sort(max_cluster_sizes);
  367. for c = 1:numel(clustsizes)
  368. p_value_cluster(c) = (nperms - dsearchn(mx, sum_t(c))) / nperms;
  369. end
  370. %% Trial-by-trial fluctuations in alpha phase predict perceptual sensitivity only in lower IAF trials.
  371. clc,clear
  372. rng shuffle
  373. load('PHASE_all.mat');
  374. load('IAF_allSubjects.mat')
  375. load('timeperiod_phase.mat')
  376. phase2analyse = 5;
  377. phaseidx = dsearchn(timeperiod', phase2analyse');
  378. phaseidx = unique(phaseidx);
  379. nS = numel(datax);
  380. iaf_coeff = zeros(nS,1);
  381. pha_coeff = zeros(nS,1);
  382. iaf_pha_c = zeros(nS,1);
  383. for participant = 1:nS
  384. X_iaf = datax{participant};
  385. accuracy = (X_iaf(:,1) == X_iaf(:,2)); % 0/1
  386. iaf = zscore( mean(X_iaf(:,4:end),2) );
  387. X_phase = data_phase_all{participant};
  388. phase = X_phase(:, phaseidx(1)+3);
  389. phase_degree = rad2deg(phase);
  390. numBins = 2;
  391. binEdges = linspace(-180, 180, numBins + 1);
  392. binIndices = zeros(size(phase_degree));
  393. for i = 1:numBins
  394. binIndices( phase_degree >= binEdges(i) & phase_degree <= binEdges(i+1) ) = i;
  395. end
  396. peak = (binIndices == 2);
  397. X = [ones(size(X_phase,1),1), iaf, double(peak), iaf .* double(peak)];
  398. t = (X'*X)\(X'*accuracy);
  399. iaf_coeff(participant) = t(2);
  400. pha_coeff(participant) = t(3);
  401. iaf_pha_c(participant) = t(4);
  402. iaf_coeff_p = zeros(1000,1);
  403. pha_coeff_p = zeros(1000,1);
  404. iaf_pha_p = zeros(1000,1);
  405. for permi = 1:1000
  406. acc_perm = Shuffle(accuracy);
  407. tperm = (X'*X)\(X'*acc_perm);
  408. iaf_coeff_p(permi) = tperm(2);
  409. pha_coeff_p(permi) = tperm(3);
  410. iaf_pha_p(permi) = tperm(4);
  411. end
  412. iaf_coeff(participant) = (iaf_coeff(participant) - mean(iaf_coeff_p)) ./ std(iaf_coeff_p);
  413. pha_coeff(participant) = (pha_coeff(participant) - mean(pha_coeff_p)) ./ std(pha_coeff_p);
  414. iaf_pha_c(participant) = (iaf_pha_c(participant) - mean(iaf_pha_p)) ./ std(iaf_pha_p);
  415. end
  416. %% Correct vs. incorrect decisions are associated with a different phase angle only in the low IAF group
  417. clc,clear
  418. rng Shuffle
  419. load('timeperiod_phase.mat')
  420. T = numel(timeperiod);
  421. which_group = 'below'; % 'below' | 'above'
  422. switch which_group
  423. case 'below'
  424. load('PHASE_below.mat');
  425. datax = data_phase_below;
  426. group = 1:numel(datax);
  428. case 'above'
  429. load('PHASE_above.mat');
  430. datax = data_phase_above;
  431. group = 1:numel(datax);
  432. end
  433. circolar_mean = zeros(numel(datax), T, 2);
  435. for p = 1:numel(datax)
  436. data = datax{p};
  437. accuracy = (data(:,1) == data(:,2));
  438. for t = 1:T
  439. circolar_mean(p,t,1) = circ_mean(data( accuracy==1, t+3)); % Correct
  440. circolar_mean(p,t,2) = circ_mean(data( accuracy==0, t+3)); % Incorrect
  441. end
  442. end
  443. p_obs = zeros(T,1);
  444. F_obs = zeros(T,1);
  445. for t = 1:T
  446. [p_obs(t), tbl] = circ_wwtest(circolar_mean(:,t,1), circolar_mean(:,t,2));
  447. F_obs(t) = tbl{2,5};
  448. end
  449. is_sig = (p_obs < 0.05);
  450. CC = bwconncomp(is_sig);
  451. clusters_obs = CC.PixelIdxList;
  452. mass_obs = zeros(numel(clusters_obs),1);
  453. for c = 1:numel(clusters_obs)
  454. idx = clusters_obs{c};
  455. mass_obs(c) = sum(F_obs(idx));
  456. end
  457. nperms = 1000;
  458. T = size(circolar_mean,2);
  459. nSubj = size(circolar_mean,1);
  460. max_mass = zeros(nperms,1);
  461. for perm = 1:nperms
  462. cm_tmp = circolar_mean;
  463. swap_mask = false(nSubj,1);
  464. swap_mask(randperm(nSubj, floor(nSubj/2))) = true;
  465. tmpMass = cm_tmp(swap_mask,:,1);
  466. cm_tmp(swap_mask,:,1) = cm_tmp(swap_mask,:,2);
  467. cm_tmp(swap_mask,:,2) = tmpMass;
  468. pvec = zeros(T,1); Fvec = zeros(T,1);
  469. for t = 1:T
  470. [pvec(t), tbl] = circ_wwtest(cm_tmp(:,t,1), cm_tmp(:,t,2));
  471. Fvec(t) = tbl{2,5};
  472. end
  473. is_sig_p = (pvec < 0.05);
  474. CCp = bwconncomp(is_sig_p);
  475. if CCp.NumObjects==0
  476. max_mass(perm) = 0;
  477. else
  478. masses = zeros(CCp.NumObjects,1);
  479. for c = 1:CCp.NumObjects
  480. ii = CCp.PixelIdxList{c};
  481. masses(c) = sum(Fvec(ii));
  482. end
  483. max_mass(perm) = max(masses);
  484. end
  485. end
  486. p_value_cluster = nan(numel(mass_obs),1);
  487. sort_max = sort(max_mass);
  488. for c = 1:numel(mass_obs)
  489. p_value_cluster(c) = (nperms - dsearchn(sort_max, mass_obs(c))) / nperms;
  490. end
  491. %% Alpha-phase clustering increases for correct responses only in low IAF individuals (Time-resolved)
  492. clc, clear
  493. rng('Shuffle')
  494. load('timeperiod_phase.mat')
  495. NT = numel(timeperiod);
  496. which_group = 'below'; % 'below' | 'above'
  497. switch which_group
  498. case 'below'
  499. load('PHASE_below.mat');
  500. datax = data_phase_below;
  501. group = 1:numel(datax);
  502. case 'above'
  503. load('PHASE_above.mat');
  504. datax = data_phase_above;
  505. group = 1:numel(datax);
  506. end
  507. no_perm = 500;
  508. itpc = zeros(length(datax), NT, 2);
  509. for participant = 1:length(datax)
  510. participant
  511. data = datax{participant};
  512. accuracy = (data(:,1) == data(:,2));
  513. nCorr = sum(accuracy==1);
  514. nErr = sum(accuracy==0);
  515. nMin = min(nCorr,nErr);
  516. for t = 1:NT
  517. itpc_s = zeros(no_perm,1);
  518. for r = 1:no_perm
  519. temp = Shuffle(find(accuracy==1));
  520. temp = temp(1:nMin);
  521. itpc_s(r,1) = abs(mean(exp(1i*data(temp , t+3))));
  522. end
  523. itpc(participant,t,1) = mean(itpc_s);
  524. itpc(participant,t,2) = abs(mean(exp(1i*data(accuracy==0, t+3))));
  525. end
  526. end
  527. signP = zeros(NT,1);
  528. tstat = zeros(NT,1);
  529. for t = 1:NT
  530. [~,signP(t),~,st] = ttest(itpc(:,t,1), itpc(:,t,2));
  531. tstat(t) = st.tstat;
  532. end
  533. tmask = tstat; tmask(signP>0.05) = 0;
  534. CC = bwconncomp(tmask);
  535. clists = CC.PixelIdxList;
  536. sum_t = zeros(numel(clists),1);
  537. for c = 1:numel(clists)
  538. idx = clists{c};
  539. sum_t(c) = sum(abs(tmask(idx(1):idx(end))));
  540. end
  541. nperms = 1000;
  542. sig_fake = zeros(NT,nperms);
  543. t_fake = zeros(NT,nperms);
  544. max_cluster_sizes = zeros(nperms,1);
  545. for perm = 1:nperms
  546. lowalpha = itpc(:,:,1); % Correct(bal)
  547. highalpha = itpc(:,:,2); % Incorrect
  548. randorder = randperm(size(itpc,1));
  549. randorder_1 = rem(randorder,2);
  550. lowalpha_fake = zeros(size(lowalpha));
  551. highalpha_fake = zeros(size(highalpha));
  552. for s = 1:size(itpc,1)
  553. if randorder_1(s)==1
  554. lowalpha_fake(s,:) = lowalpha(s,:);
  555. highalpha_fake(s,:) = highalpha(s,:);
  556. else
  557. lowalpha_fake(s,:) = highalpha(s,:);
  558. highalpha_fake(s,:) = lowalpha(s,:);
  559. end
  560. end
  561. for t = 1:NT
  562. [sig_fake(t,perm),~,~,e] = ttest(lowalpha_fake(:,t), highalpha_fake(:,t));
  563. t_fake(t,perm) = e.tstat;
  564. end
  565. t_fake(~sig_fake(:,perm),perm) = 0;
  566. CCp = bwconncomp(t_fake(:,perm));
  567. tmpMass = zeros(max([CCp.NumObjects,1]),1);
  568. if CCp.NumObjects>0
  569. for c = 1:CCp.NumObjects
  570. idx = CCp.PixelIdxList{c};
  571. tmpMass(c) = sum(abs(t_fake(idx(1):idx(end),perm)));
  572. end
  573. else
  574. tmpMass(1) = 0;
  575. end
  576. max_cluster_sizes(perm) = max(tmpMass);
  577. end
  578. p_value_cluster = zeros(numel(clists),1);
  579. mx = sort(max_cluster_sizes);
  580. for c = 1:numel(clists)
  581. p_value_cluster(c) = (nperms - dsearchn(mx, sum_t(c))) / nperms;
  582. end
  583. %% Alpha-phase clustering increases for correct responses only in low IAF individuals (Time-Frequency)
  584. clc, clearvars
  585. rng('shuffle');
  586. load('timeperiod_phase.mat')
  587. load('idx_group.mat')
  588. load('PHASE_all.mat');
  589. below_median_indices = idx_group(:,1);
  590. above_median_indices = idx_group(:,2);
  591. eegfolder = ''%path
  592. group_mode = 'below';
  593. switch group_mode
  594. case 'below'
  595. subj_idx = below_median_indices;
  596. case 'above'
  597. subj_idx = above_median_indices;
  598. end
  599. EEG_struct = dir(fullfile(eegfolder, '*.set'));
  600. EEG_files = {EEG_struct.name}';
  601. EEG_files_sel = EEG_files(subj_idx);
  602. PHASE_files_sel = data_phase_all(subj_idx);
  603. elec_list_all = readcell('electrodes_to_analyze.xlsx');
  604. elec_list = elec_list_all(subj_idx);
  605. EEG = pop_loadset('filename', EEG_files_sel{1}, 'filepath', eegfolder);
  606. EEG = eeg_checkset(EEG);
  607. freqs2use = linspace(2,50,50);
  608. tvec = -1.5:1/EEG.srate:1.5;
  609. half_wavelet = (length(tvec)-1)/2;
  610. range_cycles = [3 11];
  611. s = logspace(log10(range_cycles(1)),log10(range_cycles(end)),numel(freqs2use)) ./ (2*pi*freqs2use);
  612. n_wavelet = numel(tvec);
  613. cmwX_template = cell(numel(freqs2use),1);
  614. for fi = 1:numel(freqs2use)
  615. wavelet = exp(1i*2*pi*freqs2use(fi).*tvec) .* exp(-tvec.^2./(2*s(fi)^2));
  616. cmwX_template{fi} = wavelet;
  617. end
  618. nSubj = numel(EEG_files_sel);
  619. NT = numel(timeperiod);
  620. itpc = zeros(nSubj, numel(freqs2use), NT, 2);
  621. for s = 1:nSubj
  622. s
  623. data = PHASE_files_sel{s};
  624. accuracy = (data(:,1) == data(:,2));
  625. eeg_file = fullfile(eegfolder, EEG_files_sel{s});
  626. EEG = pop_loadset('filename', EEG_files_sel{s}, 'filepath', eegfolder);
  627. EEG = eeg_checkset(EEG);
  628. elec_name = elec_list{s};
  629. chan_idx = find(strcmp({EEG.chanlocs.labels}, elec_name));
  630. n_data = size(EEG.data,2)*size(EEG.data,3);
  631. n_conv = n_wavelet + n_data - 1;
  632. TF_phase = zeros(numel(freqs2use), size(EEG.data,2), size(EEG.data,3));
  633. chandat = reshape(EEG.data(chan_idx,:,:), 1, []);
  634. dataX = fft(chandat, n_conv);
  635. for fi = 1:numel(freqs2use)
  636. wavelet = cmwX_template{fi};
  637. cmw = fft(wavelet, n_conv);
  638. cmw = cmw ./ max(cmw);
  639. as = ifft(cmw .* dataX);
  640. as = as(half_wavelet+1:end-half_wavelet);
  641. as = reshape(as, size(EEG.data,2), size(EEG.data,3));
  642. TF_phase(fi,:,:) = angle(as);
  643. end
  644. nCorr = sum(accuracy==1);
  645. nErr = sum(accuracy==0);
  646. nMin = min(nCorr,nErr);
  647. time2saveidx = dsearchn(EEG.times', [-800 200]');
  648. itpc_s = zeros(numel(freqs2use), size(TF_phase,2), 500);
  649. corr_idx = find(accuracy==1);
  650. for subs = 1:500
  651. temp = Shuffle(corr_idx);
  652. temp = temp(1:nMin);
  653. for fi = 1:numel(freqs2use)
  654. itpc_s(fi,:,subs) = abs(mean(exp(1i*TF_phase(fi,:,temp)),3));
  655. end
  656. end
  657. itpc(s,:,:,1) = mean(itpc_s(:, time2saveidx(1):time2saveidx(2), :), 3);
  658. err_idx = (accuracy==0);
  659. for fi = 1:numel(freqs2use)
  660. itpc(s,fi,:,2) = abs(mean(exp(1i*TF_phase(fi, time2saveidx(1):time2saveidx(2), err_idx)), 3));
  661. end
  662. end
  663. itpc_corr = itpc(:,:,:,1);
  664. itpc_incorr = itpc(:,:,:,2);
  665. itpc_diff = itpc_corr - itpc_incorr;
  666. npart = size(itpc_diff,1);
  667. pval = 0.05/2;
  668. zval = abs(norminv(pval));
  669. n_permutes = 1000;
  670. permmaps = zeros(n_permutes, numel(freqs2use), numel(timeperiod));
  671. tf3d = {itpc_corr, itpc_incorr};
  672. for permi = 1:n_permutes
  673. group_temp = zeros(numel(freqs2use), numel(timeperiod));
  674. randorder = randperm(npart);
  675. randorder_1 = rem(randorder,2)+1;
  676. randorder_2 = rem(randorder_1,2)+1;
  677. for subj = 1:length(randorder_1)
  678. A = tf3d{randorder_1(subj)}(subj,:,:);
  679. B = tf3d{randorder_2(subj)}(subj,:,:);
  680. group_temp = group_temp + squeeze(B(1,:,:)-A(1,:,:));
  681. end
  682. permmaps(permi,:,:) = group_temp./subj;
  683. end
  684. mean_h0 = squeeze(mean(permmaps,1));
  685. std_h0 = squeeze(std(permmaps,0,1));
  686. zmap = (squeeze(nanmean(itpc_diff,1)) - mean_h0) ./ std_h0;
  687. zmap(abs(zmap)<zval) = 0;
  688. max_cluster_sizes = zeros(1,n_permutes);
  689. for permi = 1:n_permutes
  690. threshimg = squeeze(permmaps(permi,:,:));
  691. threshimg = (threshimg-mean_h0) ./ std_h0;
  692. threshimg(abs(threshimg)<zval) = 0;
  693. islands = bwconncomp(threshimg);
  694. if numel(islands.PixelIdxList)>0
  695. tempclustsizes = cellfun(@length, islands.PixelIdxList);
  696. max_cluster_sizes(permi) = max(tempclustsizes);
  697. end
  698. end
  699. mx = sort(max_cluster_sizes);
  700. islands = bwconncomp(zmap);
  701. nClust = islands.NumObjects;
  702. p_cluster = zeros(nClust,1);
  703. for c = 1:nClust
  704. this_size = numel(islands.PixelIdxList{c});
  705. rank_idx = dsearchn(mx', this_size);
  706. p_cluster(c) = 1-(rank_idx/n_permutes);
  707. if p_cluster(c) >= 0.05
  708. zmap(islands.PixelIdxList{c}) = 0;
  709. end
  710. end
  711. %% SI: No significant differences observed in behavioural performance between first and third power terciles
  712. clc; clearvars;
  713. rng('shuffle');
  714. eegfolder = ''%path
  715. load('AllSubjects.mat')
  716. cd(eegfolder);
  717. EEG_struct = dir(fullfile(eegfolder, '*.set'));
  718. EEG_files = {EEG_struct.name}';
  719. elec_to_select = readcell('electrodes_to_analyze.xlsx');
  720. iaf_to_select = readcell('peak_to_analyze.xlsx');
  721. freqs2use = linspace(2,50,50);
  722. nSubj = numel(EEG_files);
  723. d = zeros(nSubj,2);
  724. c = zeros(nSubj,2);
  725. for s = 1:nSubj
  726. filename = EEG_files{s};
  727. EEG = pop_loadset('filename', filename, 'filepath', eegfolder);
  728. EEG = eeg_checkset(EEG);
  729. EEG.data = diff(EEG.data,1,2);
  730. elec2anal = elec_to_select{s};
  731. chan_idx = find(strcmpi({EEG.chanlocs.labels}, elec2anal));
  732. TF_abs = TF_IAF(EEG.data, EEG.srate, chan_idx); % [F x T x trials]
  733. mean_pow = mean(TF_abs,3);
  734. basewin = [-3200 -2800];
  735. baseidx = round(dsearchn(EEG.times', basewin'));
  736. bslpow = mean(mean_pow(:, baseidx(1):baseidx(2)), 2);
  737. BSLpow = repmat(bslpow, [1, size(TF_abs,2), size(TF_abs,3)]);
  738. ampl_all = 10*log10(TF_abs ./ BSLpow);
  739. iaf2select = iaf_to_select{s};
  740. iafrange2select = dsearchn(freqs2use', iaf2select');
  741. time2select = [-800 -100];
  742. timerange2select = dsearchn(EEG.times', [time2select(1) time2select(2)]');
  743. subpow = ampl_all(iafrange2select, timerange2select(1):timerange2select(2), :);
  744. meanampl = squeeze(mean(subpow, 2));
  745. meanampl = meanampl(:);
  746. [~, sorted_alpha_ampl] = sort(meanampl);
  747. nTrials = numel(sorted_alpha_ampl);
  748. n1 = ceil(nTrials/3);
  749. n3start = ceil(2*nTrials/3);
  750. first_tercile = sorted_alpha_ampl(1:n1);
  751. third_tercile = sorted_alpha_ampl(n3start:end);
  752. data = databehx{s};
  753. data = data(:,1:2);
  754. data_low_alpha = data(first_tercile,:);
  755. stim_pres = data_low_alpha(data_low_alpha(:,1)==1,:);
  756. hit = sum(stim_pres(:,1)==stim_pres(:,2)) / size(stim_pres,1);
  757. if hit==1, hit = (size(stim_pres,1)-0.5)/size(stim_pres,1); end
  758. if hit==0, hit = 0.5/size(stim_pres,1); end
  759. stim_abs = data_low_alpha(data_low_alpha(:,1)==0,:);
  760. fa = sum(stim_abs(:,1)~=stim_abs(:,2)) / size(stim_abs,1);
  761. if fa==0, fa = 0.5/size(stim_abs,1); end
  762. d(s,1) = norminv(hit) - norminv(fa);
  763. c(s,1) = -(norminv(hit) + norminv(fa))/2;
  764. data_high_alpha = data(third_tercile,:);
  765. stim_pres = data_high_alpha(data_high_alpha(:,1)==1,:);
  766. hit = sum(stim_pres(:,1)==stim_pres(:,2)) / size(stim_pres,1);
  767. if hit==1, hit = (size(stim_pres,1)-0.5)/size(stim_pres,1); end
  768. if hit==0, hit = 0.5/size(stim_pres,1); end
  769. stim_abs = data_high_alpha(data_high_alpha(:,1)==0,:);
  770. fa = sum(stim_abs(:,1)~=stim_abs(:,2)) / size(stim_abs,1);
  771. if fa==0, fa = 0.5/size(stim_abs,1); end
  772. d(s,2) = norminv(hit) - norminv(fa);
  773. c(s,2) = -(norminv(hit)+norminv(fa))/2;
  774. end
  775. [~,p_d,~,st_d] = ttest(d(:,1), d(:,2));
  776. bf_d = bf.ttest(d(:,1), d(:,2));
  777. [~,p_c,~,st_c] = ttest(c(:,1), c(:,2));
  778. bf_c = bf.ttest(c(:,1), c(:,2));
  779. %% SI: Correct and Incorrect trials are not characterized by difference in oscillatory amplitude.
  780. %% SI: Trial-by-trial fluctuations in alpha power do not account for perceptual accuracy.
  781. clc; clearvars;
  782. rng('shuffle');
  783. eegfolder = ''%path
  784. D = dir(eegfolder);
  785. folder_direeg = struct2cell(D)';
  786. folder_direeg = folder_direeg(3:end,1);
  787. folder_direeg = folder_direeg(2:2:end);
  788. all_idx = 1:116;
  789. group = all_idx;
  790. freqs2use = linspace(2,50,50);
  791. iaf_to_select = readmatrix('peak_to_analyze.xlsx');
  792. elec_to_select= readcell('electrodes_to_analyze.xlsx');
  793. load('IAF_allSubjects.mat');
  794. nSubj = numel(group);
  795. ampl = cell(nSubj,2);
  796. b_z_pow = zeros(nSubj,1);
  797. iaf_z_pow = zeros(nSubj,1);
  798. for pi = 1:nSubj
  799. subj_id = group(pi);
  800. data = datax{subj_id};
  801. accuracy = (data(:,1) == data(:,2));
  802. iaf_p = zscore(mean(data(:,4:end), 2));
  803. cd(eegfolder);
  804. filename = folder_direeg{subj_id};
  805. EEG = pop_loadset('filename', filename);
  806. EEG = eeg_checkset(EEG);
  807. EEG.data = diff(EEG.data, 1, 2);
  808. iaf_part = iaf_to_select(subj_id);
  809. elec_name = elec_to_select{subj_id};
  810. chan_idx = find(strcmpi({EEG.chanlocs.labels}, elec_name));
  811. iafidx = dsearchn(freqs2use', iaf_part);
  812. TF_amp = TF_IAF(EEG.data, EEG.srate, chan_idx);
  813. time2save = [-800 -100];
  814. time2saveidx = dsearchn(EEG.times', time2save');
  815. basewin = [-3200 -2800];
  816. baseidx = round(dsearchn(EEG.times', basewin'));
  817. mean_pow_corr = mean(TF_amp(:,:,accuracy), 3);
  818. bslpow_CORR = mean(mean_pow_corr(:, baseidx(1):baseidx(2)), 2);
  819. mean_pow_incorr = mean(TF_amp(:,:,accuracy==0), 3);
  820. bslpow_INCORR = mean(mean_pow_incorr(:, baseidx(1):baseidx(2)), 2);
  821. ampl{pi,1} = 10*log10( mean(TF_amp(:, time2saveidx(1):time2saveidx(2), accuracy), 3) ./ bslpow_CORR );
  822. ampl{pi,2} = 10*log10( mean(TF_amp(:, time2saveidx(1):time2saveidx(2), accuracy==0),3) ./ bslpow_INCORR );
  823. mean_pow_all = mean(TF_amp, 3);
  824. bslpow_all = mean(mean_pow_all(:, baseidx(1):baseidx(2)), 2);
  825. BSLpow = repmat(bslpow_all, [1, size(TF_amp,2), size(TF_amp,3)]);
  826. powtrialbytrial = TF_amp ./ BSLpow;
  827. ampl2use = squeeze( zscore( mean(powtrialbytrial(iafidx, time2saveidx(1):time2saveidx(2), :), 2) ) );
  828. ampl2use = ampl2use(:);
  829. X = [ones(length(ampl2use),1), ampl2use, iaf_p];
  830. beta = (X' * X) \ (X' * accuracy);
  831. b_pow = beta(2);
  832. b_iaf = beta(3);
  833. nPermReg = 2000;
  834. b_perms_ampl = zeros(nPermReg,1);
  835. b_perms_iaf = zeros(nPermReg,1);
  836. for permi = 1:nPermReg
  837. rp = randperm(length(ampl2use));
  838. pow_r = ampl2use(rp);
  839. iaf_r = iaf_p(rp);
  840. Xp = [ones(length(pow_r),1), pow_r, iaf_r];
  841. beta_p = (Xp' * Xp) \ (Xp' * accuracy);
  842. b_perms_ampl(permi) = beta_p(2);
  843. b_perms_iaf(permi) = beta_p(3);
  844. end
  845. b_z_pow(pi) = (b_pow - mean(b_perms_ampl)) ./ std(b_perms_ampl);
  846. iaf_z_pow(pi) = (b_iaf - mean(b_perms_iaf)) ./ std(b_perms_iaf);
  847. end
  848. [~,p_b,~,st_b] = ttest(b_z_pow);
  849. bf_b = bf.ttest(b_z_pow);
  850. [~,p_i,~,st_i] = ttest(iaf_z_pow);
  851. bf_i = bf.ttest(iaf_z_pow);
  852. [F,T] = size(ampl{1,1});
  853. correct_value = zeros(nSubj,F,T);
  854. incorrect_value = zeros(nSubj,F,T);
  855. for i = 1:nSubj
  856. correct_value(i,:,:) = ampl{i,1};
  857. incorrect_value(i,:,:) = ampl{i,2};
  858. end
  859. correct_responses = squeeze(mean(correct_value,1)); % [F x T]
  860. incorrect_responses = squeeze(mean(incorrect_value,1)); % [F x T]
  861. difference_matrix = correct_responses - incorrect_responses;
  862. pval = 0.05/2;
  863. zval = abs(norminv(pval));
  864. time2save = [-800 -100];
  865. time2saveidx = dsearchn(EEG.times', time2save');
  866. timeperiod = EEG.times(time2saveidx(1):time2saveidx(2));
  867. num_zeros = nSubj/2;
  868. num_ones = nSubj/2;
  869. vector = [zeros(1,num_zeros), ones(1,num_ones)] + 1;
  870. meanamplcorr = {correct_value, incorrect_value};
  871. nPermTF = 1000;
  872. difference_matrix_perm = zeros(nPermTF,F,T);
  873. for permi = 1:nPermTF
  874. v = vector(randperm(nSubj));
  875. sv = 3 - v;
  876. corrtemp = zeros(F,T);
  877. incorrtemp = zeros(F,T);
  878. for part = 1:nSubj
  879. corrtemp = corrtemp + squeeze(meanamplcorr{v(part)} (part,:,:));
  880. incorrtemp = incorrtemp + squeeze(meanamplcorr{sv(part)}(part,:,:));
  881. end
  882. difference_matrix_perm(permi,:,:) = (corrtemp - incorrtemp) ./ nSubj;
  883. end
  884. mean_h0 = squeeze(mean(difference_matrix_perm,1));
  885. std_h0 = squeeze(std(difference_matrix_perm,[],1));
  886. zmap = (difference_matrix - mean_h0) ./ std_h0;
  887. zmap(abs(zmap) < zval) = 0;
  888. max_cluster_sizes = zeros(1, nPermTF);
  889. for permi = 1:nPermTF
  890. thr = squeeze(difference_matrix_perm(permi,:,:));
  891. thr = (thr - mean_h0) ./ std_h0;
  892. thr(abs(thr) < zval) = 0;
  894. islands = bwconncomp(thr);
  895. if numel(islands.PixelIdxList) > 0
  896. tempclustsizes = cellfun(@length, islands.PixelIdxList);
  897. max_cluster_sizes(permi) = max(tempclustsizes);
  898. else
  899. max_cluster_sizes(permi) = 0;
  900. end
  901. end
  902. mx = sort(max_cluster_sizes);
  903. islands = bwconncomp(zmap);
  904. nClust = islands.NumObjects;
  905. p_cluster = zeros(nClust,1);
  906. for c = 1:nClust
  907. this_size = numel(islands.PixelIdxList{c});
  908. rank_idx = dsearchn(mx', this_size);
  909. p_cluster(c) = 1 - (rank_idx / nPermTF);
  910. if p_cluster(c) >= 0.05
  911. zmap(islands.PixelIdxList{c}) = 0;
  912. end
  913. end
  914. %% SI: Low and Fast IAF groups do not exhibit significant differences in oscillatory amplitude.
  915. clc; clearvars;
  916. rng('shuffle');
  917. load('idx_group.mat');
  918. group_low = idx_group(:,1);
  919. group_high = idx_group(:,2);
  920. nLow = numel(group_low);
  921. nHigh = numel(group_high);
  922. iaf_to_select = readmatrix('peak_to_analyze.xlsx');
  923. elec_to_select = readcell('electrodes_to_analyze.xlsx');
  924. eegfolder = '';%path
  925. cd(eegfolder);
  926. D = dir(eegfolder);
  927. folder_direeg = struct2cell(D)';
  928. folder_direeg = folder_direeg(3:end,1);
  929. folder_direeg = folder_direeg(2:2:end);
  930. nSubj = numel(folder_direeg);
  931. freqs2use = linspace(2,50,50);
  932. EEG = pop_loadset('filename', folder_direeg{1}, 'filepath', eegfolder);
  933. EEG = eeg_checkset(EEG);
  934. time_window = [-800 -100];
  935. time2idx = dsearchn(EEG.times', time_window');
  936. timeperiod = EEG.times(time2idx(1):time2idx(2));
  937. basewin = [-3200 -2800];
  938. baseidx = round(dsearchn(EEG.times', basewin'));
  939. nFreq = numel(freqs2use);
  940. nTimeWin = numel(timeperiod);
  941. tf_ampl_all = NaN(nSubj, nFreq, nTimeWin);
  942. for s = 1:nSubj
  943. part_elec = elec_to_select{s};
  944. filename = folder_direeg{s};
  945. EEG = pop_loadset('filename', filename, 'filepath', eegfolder);
  946. EEG = eeg_checkset(EEG);
  947. EEG.data = diff(EEG.data,1,2);
  948. elec_idx = find(strcmpi({EEG.chanlocs.labels}, part_elec));
  949. TF_abs = TF_IAF(EEG.data, EEG.srate, elec_idx);
  950. bslpow = mean(mean(TF_abs(:, baseidx(1):baseidx(2), :), 2), 3);
  951. meanampl = mean(TF_abs, 3);
  952. BSLpow = repmat(bslpow, [1, size(meanampl,2)]);
  953. ampl_all = 10*log10(meanampl ./ BSLpow);
  954. tf_ampl_all(s,:,:) = ampl_all(:, time2idx(1):time2idx(2));
  955. end
  956. low_tf = squeeze(mean(tf_ampl_all(group_low,:,:), 1));
  957. high_tf = squeeze(mean(tf_ampl_all(group_high,:,:), 1));
  958. diff_tf = low_tf - high_tf;
  959. pval = 0.05/2;
  960. zval = abs(norminv(pval));
  961. nPermTF = 1000;
  962. all_idx = [group_low(:); group_high(:)];
  963. nTot = numel(all_idx);
  964. perm_maps = zeros(nPermTF, nFreq, nTimeWin);
  965. for permi = 1:nPermTF
  966. rp = randperm(nTot);
  967. fake_low_idx = all_idx(rp(1:nLow));
  968. fake_high_idx = all_idx(rp(nLow+1:end));
  969. fake_low_tf = squeeze(mean(tf_ampl_all(fake_low_idx,:,:), 1));
  970. fake_high_tf = squeeze(mean(tf_ampl_all(fake_high_idx,:,:), 1));
  971. perm_maps(permi,:,:) = fake_low_tf - fake_high_tf;
  972. end
  973. mean_h0 = squeeze(mean(perm_maps,1));
  974. std_h0 = squeeze(std(perm_maps,0,1));
  975. zmap = (diff_tf - mean_h0) ./ std_h0;
  976. zmap(abs(zmap) < zval) = 0;
  977. max_cluster_sizes = zeros(1, nPermTF);
  978. for permi = 1:nPermTF
  979. thr = squeeze(perm_maps(permi,:,:));
  980. thr = (thr - mean_h0) ./ std_h0;
  981. thr(abs(thr) < zval) = 0;
  982. islands = bwconncomp(thr);
  983. if numel(islands.PixelIdxList) > 0
  984. temp_sizes = cellfun(@length, islands.PixelIdxList);
  985. max_cluster_sizes(permi) = max(temp_sizes);
  986. else
  987. max_cluster_sizes(permi) = 0;
  988. end
  989. end
  990. mx = sort(max_cluster_sizes(:));
  991. islands = bwconncomp(zmap);
  992. nClust = islands.NumObjects;
  993. p_cluster = zeros(nClust,1);
  994. for c = 1:nClust
  995. this_size = numel(islands.PixelIdxList{c});
  996. rank_idx = dsearchn(mx, this_size);
  997. p_cluster(c) = 1 - (rank_idx / nPermTF);
  998. if p_cluster(c) >= 0.05
  999. zmap(islands.PixelIdxList{c}) = 0;
  1000. end
  1001. end
  1002. %% SI : FFT-based method for estimating single-trial IAF confirmed the pattern
  1003. %% of results observed with the instantaneous frequency approach.
  1004. %% SI : The FFT-based and instantaneous frequency approaches result in comparable single-trial IAF values.
  1005. clc; clear; close all;
  1006. EEGFOLDER = '';%path
  1007. cd(EEGFOLDER);
  1008. dEEG = dir(EEGFOLDER);
  1009. EEG_files = struct2cell(dEEG)';
  1010. EEG_files = EEG_files(3:end,1);
  1011. start_from = 2;
  1012. subji = start_from:2:length(EEG_files);
  1013. dataEEG = cell(1, numel(subji));
  1014. cd(EEGFOLDER);
  1015. for p = 1:numel(subji)
  1016. idxEEG = subji(p);
  1017. filename = char(EEG_files{idxEEG});
  1018. EEG = pop_loadset('filename', filename);
  1019. EEG = eeg_checkset(EEG);
  1020. dataEEG{p} = EEG;
  1021. end
  1022. electrodes_to_analyze = readcell('electrodes_to_analyze.xlsx');
  1023. cmin = 1;
  1024. fRange = [1 40];
  1025. Fw = 11;
  1026. k = 5;
  1027. iaf_all = cell(1, size(subji,2));
  1028. god_all = cell(1, size(subji,2));
  1029. for soggetto = 1:size(subji,2)
  1030. EEG = dataEEG{soggetto};
  1031. EEG.data = diff(EEG.data,1,2);
  1032. time2anal = [-800 -100];
  1033. time2idx = dsearchn(EEG.times',time2anal');
  1034. Fs = EEG.srate;
  1035. idxElectrode = find(strcmpi({EEG.chanlocs.labels}, ...
  1036. electrodes_to_analyze{soggetto}));
  1037. w = [7 14];
  1038. ntrials = size(EEG.data,3);
  1039. iaf_p = zeros(ntrials,1);
  1040. god_p = zeros(ntrials,1);
  1041. for triali = 1:ntrials
  1042. data_colcoran = EEG.data(idxElectrode,time2idx(1):time2idx(2),triali);
  1043. [~,colcalpha] = restingIAF_natcomm(data_colcoran, 1, cmin, fRange, Fs, w, Fw, k,'nfft', EEG.srate*4);
  1044. iaf_p(triali) = colcalpha.peaks;
  1045. god_p(triali) = ~isnan(iaf_p(triali));
  1046. end
  1047. iaf_all{soggetto} = iaf_p;
  1048. god_all{soggetto} = god_p;
  1049. end
  1050. load('IAF_allSubjects.mat')
  1051. iaf = zeros(length(datax),2);
  1052. gof = zeros(length(datax),2);
  1053. nperms = 2000;
  1054. [d,c,accuracies] = deal(zeros(length(datax),2));
  1055. b = zeros(1,length(iaf_all));
  1056. b_z = zeros(1,length(iaf_all));
  1057. rho = NaN(numel(iaf_all),1);
  1058. for participant = 1:length(iaf_all)
  1059. databeh = datax{participant};
  1060. iaf_inst = mean(databeh(:,4:end),2);
  1061. databeh = databeh(:,1:3);
  1062. accuracy = databeh(:,1) == databeh(:,2);
  1063. data_iaf = iaf_all{participant};
  1064. rho(participant) = corr(data_iaf, iaf_inst,'Type','Spearman','Rows','complete');
  1065. iaf_acc = data_iaf(accuracy);
  1066. iaf_inc = data_iaf(~accuracy);
  1067. iaf(participant,1) = nanmean(iaf_acc);
  1068. iaf(participant,2) = nanmean(iaf_inc);
  1069. data_gof = god_all{participant};
  1070. gof_acc = data_gof(accuracy);
  1071. gof_inc = data_gof(~accuracy);
  1072. gof(participant,1) = mean(gof_acc);
  1073. gof(participant,2) = mean(gof_inc);
  1074. X = [ones(length(data_iaf),1), data_iaf];
  1075. todel = isnan(X(:,2));
  1076. X(todel,:) = [];
  1077. accuracy(todel) = [];
  1078. X(:,2) = zscore(X(:,2));
  1079. t = (X'*X)\(X'*accuracy);
  1080. b(participant) = t(2);
  1081. datat = X;
  1082. X_o = X;
  1083. b_perms = zeros(nperms,1);
  1084. for permi = 1:nperms
  1085. datat(:,2) = Shuffle(X_o(:,2));
  1086. Xt = [ones(length(datat),1), datat(:,2)];
  1087. tt = (Xt'*Xt)\(Xt'*accuracy);
  1088. b_perms(permi) = tt(2);
  1089. end
  1090. b_z(participant) = (b(participant)-mean(b_perms))./std(b_perms);
  1091. todel = isnan(data_iaf);
  1092. data_iaf(todel) = [];
  1093. databeh(todel,:) = [];
  1094. [~, h] = sort(data_iaf);
  1095. databeh = databeh(h,:);
  1096. data_low_alpha = databeh(1:ceil(length(data_iaf)/3),:);
  1097. data_hig_alpha = databeh(ceil(length(data_iaf)/3*2)+1:end,:);
  1098. stim_pres = data_low_alpha( data_low_alpha(:,1)==1 ,:);
  1099. hit = (sum(stim_pres(:,1) == stim_pres(:,2)))/(length(stim_pres));
  1100. stim_abs = data_low_alpha( data_low_alpha(:,1)==0 ,:);
  1101. fa = (sum(stim_abs(:,1) ~= stim_abs(:,2)))/(length(stim_abs));
  1102. if hit==0, hit=0.5/length(stim_pres); elseif hit==1, hit=(length(stim_pres)-0.5)/length(stim_pres); end
  1103. if fa==0, fa=0.5/length(stim_abs); elseif fa==1, fa=(nN-0.5)/length(stim_abs); end
  1104. d(participant,1) = norminv(hit) - norminv(fa);
  1105. c(participant,1) = -(norminv(hit) + norminv(fa))/2;
  1106. accuracies(participant,1) = sum(data_low_alpha(:,1) == data_low_alpha(:,2)) ...
  1107. / length(data_low_alpha(:,1));
  1108. stim_pres = data_hig_alpha( data_hig_alpha(:,1)==1 ,:);
  1109. hit = (sum(stim_pres(:,1) == stim_pres(:,2)))/(length(stim_pres));
  1110. stim_abs = data_hig_alpha( data_hig_alpha(:,1)==0 ,:);
  1111. fa = (sum(stim_abs(:,1) ~= stim_abs(:,2)))/(length(stim_abs));
  1112. if hit==0, hit=0.5/length(stim_pres); elseif hit==1, hit=(length(stim_pres)-0.5)/length(stim_pres); end
  1113. if fa==0, fa=0.5/length(stim_abs); elseif fa==1, fa=(nN-0.5)/length(stim_abs); end
  1114. d(participant,2) = norminv(hit) - norminv(fa);
  1115. c(participant,2) = -(norminv(hit) + norminv(fa))/2;
  1116. accuracies(participant,2) = sum(data_hig_alpha(:,1) == data_hig_alpha(:,2)) ...
  1117. / length(data_hig_alpha(:,1));
  1118. end
  1119. [a,b,cc,dd] = ttest(iaf(:,1),iaf(:,2))
  1120. [bf10a] = bf.ttest(iaf(:,1),iaf(:,2))
  1121. [a,b,cc,dd] = ttest(gof(:,1),gof(:,2))
  1122. [bf1] = bf.ttest(gof(:,1),gof(:,2))
  1123. [aa,bb,cc,ddd] = ttest(b_z)
  1124. [bf10b] = bf.ttest(b_z)
  1125. [a,b,cc,dd] = ttest(d(:,1),d(:,2))
  1126. [bf10c] = bf.ttest(d(:,1),d(:,2))
  1127. [a,b,cc,dd] = ttest(c(:,1),c(:,2))
  1128. [bf10] = bf.ttest(c(:,1),c(:,2))
  1129. [a,b,cc,dd] = ttest(accuracies(:,1),accuracies(:,2))
  1130. [bf10d] = bf.ttest(accuracies(:,1),accuracies(:,2))

Script_NC.m, no license · at the source

Overview

Authors: Vincenzo Romei1,2, Luca Tarasi1
ORCID iDs: Vincenzo Romei, Luca Tarasi
  1. Dipartimento di Psicologia, Università di Bologna and Centro studi e ricerche in Neuroscienze Cognitive, Università di Bologna, Cesena, Italy
  2. Universidad Antonio de Nebrija, Madrid, Spain
Institutions: Universidad Nebrija (Spain); University of Bologna (Italy)
Journal: Nature communications, volume 17, issue 1, article 3384
Dates: received 8 March 2024; accepted 13 February 2026; published online 3 March 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-70124-9 · PMID 41776179 · PMCID PMC13065768 · OpenAlex W7133300256
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: EEG (modality), human (organism), cognitive (subfield)
Methods: Spectral & time-frequency, Preprocessing, Connectivity, Statistics, Smoothing, state filtering, decompositions, Physiology & signal measures
Keywords: Human behaviour, Perception
MeSH: Alpha Rhythm*, Perception*, Visual Perception*, Bayes Theorem, Decision Making, Electroencephalography, Female, Humans, Male, Photic Stimulation (* major topic)
Topic: Visual perception and processing mechanisms (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Citations: cited by 6 papers (Europe PMC); 84 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

Its files are read in the Code ↔ Paper reader above, with 4 matches between paragraphs and lines of code.

OSF 6298n

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Languages: MATLAB (4)
Size: 20 files, 4 scripts
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: EEGLAB (2 files), CircStat (1 file), Image Processing Toolbox (1 file), Signal Processing Toolbox (1 file), Statistics and Machine Learning Toolbox (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers (HTTP 200)
  • 30 September 2026: the link answers (HTTP 200)
4 files
At the source:

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • it points to the authors' code: OSF 6298n

Read it in the paper: doi.org/10.1038/s41467-026-70124-9.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 4 scripts, each with its path and the digest of its content;
  • 4 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • no repository, dataset or request procedure was recognized in it

Read it in the paper: doi.org/10.1038/s41467-026-70124-9.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 2 authors, 2 keywords, 10 MeSH terms, 82 references.

Cite

This paper

Romei, V., & Tarasi, L. (2026). Alpha frequency shapes perceptual sensitivity by modulating optimal phase likelihood. Nature communications, 17(1), 3384. https://doi.org/10.1038/s41467-026-70124-9

BibTeX

@article{romei2026alpha,
author = {Romei, Vincenzo and Tarasi, Luca},
title = {{Alpha frequency shapes perceptual sensitivity by modulating optimal phase likelihood}},
journal = {Nature communications},
year = {2026},
month = mar,
volume = {17},
number = {1},
pages = {3384},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-70124-9},
url = {https://doi.org/10.1038/s41467-026-70124-9},
pmid = {41776179},
pmcid = {PMC13065768}
}

RIS

TY - JOUR
AU - Romei, Vincenzo
AU - Tarasi, Luca
TI - Alpha frequency shapes perceptual sensitivity by modulating optimal phase likelihood
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/03/03
VL - 17
IS - 1
SP - 3384
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-70124-9
UR - https://doi.org/10.1038/s41467-026-70124-9
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-70124-9",
"type": "article-journal",
"title": "Alpha frequency shapes perceptual sensitivity by modulating optimal phase likelihood",
"container-title": "Nature communications",
"author": [
{
"family": "Romei",
"given": "Vincenzo"
},
{
"family": "Tarasi",
"given": "Luca"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "3384",
"DOI": "10.1038/s41467-026-70124-9",
"PMID": "41776179",
"PMCID": "PMC13065768",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-70124-9",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
3
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1162/imag.a.1263
Individual alpha frequency predicts the sensitivity of time perception.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: cognitive, 13 references
[2] doi:10.7554/elife.110000 [code]
Alpha-band phase modulates perceptual sensitivity by changing internal noise and sensory tuning.
Journal: eLife
In common: CircStat, Image Processing Toolbox, Statistics and Machine Learning Toolbox, EEG, 10 references
[3] doi:10.1093/cercor/bhag117 [code]
Spontaneous changes in the prestimulus alpha power predict both objective and subjective aspects of perception-evidence from 4 paradigms and 471 experimental sessions.
Journal: Cerebral cortex (New York, N.Y. : 1991)
In common: EEG, cognitive, 8 references
[4] doi:10.7554/elife.108408 [code]
Frequency and laminar profile of feature-specific visual activity revealed by interleaved EEG-fMRI.
Journal: eLife
In common: Image Processing Toolbox, Signal Processing Toolbox, Statistics and Machine Learning Toolbox, EEG, 5 references
[5] doi:10.1523/jneurosci.0154-26.2026 [code]
Faster but less precise: expectation enhances response speed while reducing sensory fidelity.
Journal: The Journal of neuroscience : the official journal of the Society for Neuroscience
In common: CircStat, EEGLAB, Image Processing Toolbox, 2 other tools, EEG, cognitive, 2 references
[6] doi:10.1111/psyp.70271 [code]
Disentangling Respiratory Phase-Dependent and Phase-Independent Components of Anticipatory Cardiac Deceleration.
Journal: Psychophysiology
In common: CircStat, EEGLAB, Image Processing Toolbox, 2 other tools, EEG, cognitive, 1 reference
[7] doi:10.1016/j.isci.2025.113806 [code]
Beta-band frequency shifts signal decisions in human prefrontal cortex
Journal: n/a
In common: Statistics and Machine Learning Toolbox, 6 references
[8] doi:10.1162/imag.a.1229 [code]
40 Hz audiovisual stimulation improves sustained attention and related brain oscillations.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: CircStat, EEGLAB, Image Processing Toolbox, 2 other tools, EEG, 1 reference
[9] doi:10.1162/imag.a.1199 [code]
Sustained alpha oscillations serve attentional prioritization in working memory, not maintenance.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: CircStat, EEGLAB, Image Processing Toolbox, 2 other tools, EEG, cognitive, 1 reference
[10] doi:10.1371/journal.pbio.3003818 [code]
Human neuronal firing varies with the frequency of local field potential oscillations.
Journal: PLoS biology
In common: CircStat, EEGLAB, Image Processing Toolbox, 2 other tools, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.