OSCR

Hippocampal ripples initiate cortical dimensionality expansion for memory retrieval.

Code ↔ Paper

6 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 6 matches
  1. [1] § Methods › Theta-gamma phase-amplitude coupling ↔ ripple_stage_05b_PAC_across_time_shuffled_final.m, lines 650–709 · score 0.72 · theta frequency, 140 Hz, hippocampal channels, wavelet, Hanning, taper
  2. [2] § Methods › Theta-gamma phase-amplitude coupling ↔ ripple_stage_05a_PAC_final.m, lines 607–664 · score 0.71 · theta frequency, 140 Hz, hippocampal channels, wavelet, Hanning, taper
  3. [3] § Methods › Method details › Dimensionality transformation ↔ ripple_stage_04c_dimensionality_LME_final.m, lines 603–718 · score 0.61 · elbow point, explained variance, eigenvalues, threshold, PCA, component
  4. [4] § Methods › Method details › Dimensionality transformation ↔ ripple_stage_04a_dimensionality_final.m, lines 531–662 · score 0.59 · elbow point, explained variance, eigenvalues, PCA, component, transformation
  5. [5] § Results › Hippocampal ripple-induced dimensionality expansion increases the separability of cortical representations ↔ ripple_stage_04a_dimensionality_final.m, lines 531–662 · score 0.51 · sliding windows, latent, elbow, eigenvalue, variance, overlap
  6. [6] § Results › Hippocampal ripple-induced dimensionality expansion increases the separability of cortical representations ↔ ripple_stage_04a_dimensionality_shuffled_final.m, lines 530–647 · score 0.51 · sliding windows, latent, elbow, eigenvalue, variance, overlap

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,223 lines · 43 KB · no license · 2 matches

  1. %%
  2. % [ripple_stage_04a_dimensionality] find ripples in hippocampal channels,
  3. % extract and realign data based on ripple events.
  4. % Do PCA to estimate dimensionality of correct and incorrect.
  5. % Bernhard Staresina [[email hidden]]
  6. % Casper Kerren [[email hidden]]
  8. clear
  9. restoredefaultpath
  10. addpath('/Users/kerrenadmin/Desktop/Postdoc/Project_1/Analyses_matlab/general_scripts_matlab/fieldtrip-20230422')
  11. ft_defaults
  13. % [~,ftpath]=ft_version;
  15. %% path settings
  16. settings = [];
  17. settings.base_path_castle = '/Users/kerrenadmin/Desktop/Other_projects/Dimensionality_ripples_Casper_and_Bernhard/'; % '/castles/nr/projects/w/wimberm-ieeg-compute/';
  18. settings.subjects = char('CF', 'JM', 'SO', 'AH','FC', 'HW', 'AM', 'MH','FS', 'AS', 'CB', 'KK');
  19. settings.SubjectIDs = char('01_CF', '02_JM', '03_SO', '06_AH','07_FC', '08_HW', '09_AM', '10_MH','11_FS', '12_AS', '13_CB', '14_KK');
  21. load("colour_scheme.mat")
  22. settings.colour_scheme = colour_scheme;
  24. settings.data_dir = [settings.base_path_castle,'preprocessing/artifact_rejected_data/'];
  26. settings.save_dir = [settings.base_path_castle,'output_data/decoding/'];
  27. settings.data_dir_channels = [settings.base_path_castle,'ripple_project_publication_for_replication/templates'];
  28. settings.anatomy_dir = [settings.base_path_castle,'ripple_project_publication_for_replication/additional_analyses/visualisation/'];
  29. settings.AAL_dir = fullfile(settings.base_path_castle,'ripple_project_publication_for_replication/subfunctions/AAL3');
  30. settings.SPM_dir = fullfile(settings.base_path_castle,'/ripple_project_publication_for_replication/subfunctions/spm12');
  32. settings.scalp_channels = {'C3' 'C4' 'Cz' 'T3' 'T4' 'T5' 'T6' 'O1' 'O2' 'Oz' 'F3' 'F4' 'Fz' 'Cb1' 'Cb2'};
  34. subjects = {'CF', 'JM', 'SO', 'AH','FC', 'HW', 'AM', 'MH','FS', 'AS', 'CB', 'KK'};
  35. SubjectIDs = {'01_CF', '02_JM', '03_SO', '06_AH','07_FC', '08_HW', '09_AM', '10_MH','11_FS', '12_AS', '13_CB', '14_KK'};
  36. settings.nu_rep = [1, 1, 2, 1, 1, 2, 2, 1, 1, 1, 1, 1];
  39. settings.healthyhemi = {'R' 'LR' 'R' 'L' 'L' 'R' 'R' 'R' 'R' 'R' 'R' 'LR'};
  41. addpath(genpath('/Users/kerrenadmin/Desktop/Postdoc/Project_1/Analyses_matlab/general_scripts_matlab/MVPA-Light-master'))
  42. addpath(genpath([settings.base_path_castle,'ripple_project_publication_for_replication/main_analyses/Slythm']))
  43. addpath([settings.base_path_castle,'ripple_project_publication_for_replication/subfunctions'])
  44. addpath(genpath('/Users/kerrenadmin/Desktop/Postdoc/Project_1/Analyses_matlab/help_functions'))
  45. addpath(genpath('/Users/kerrenadmin/Desktop/Postdoc/Project_1/Analyses_matlab/general_scripts_matlab/plotting'))
  47. %% pre-decoding settings
  48. % ripple extraction
  50. settings.remove_falsepositives = 1; % decide whether or not to exclude ripples deemed false positives based on spectral peak detection
  51. settings.full_enc_trial = 1; % set to 0 if you want encoding trial to end with RT and to 1 if it should end at 3 sec
  52. settings.remove_ripple_duplicates = 1; % remove co-occuring ripples
  53. settings.time_to_excl_RT = .25; % exclude last 250 ms of trials, to make sure ripple event was in trial
  54. settings.solo_ripple = 1; % pick one (maxEnv) ripple per trial if multiple ripple events are found
  55. settings.ripple_latency = [.25 5]; % define time window at retrieval in which the ripple events need to occur
  56. settings.do_surrogates = 0; % switch time of ripples between trials, 1 == for all trials, 2 == for correct trials only
  58. %% decoding settings
  60. % data preprocessing
  61. settings.ori_fsample = 1000; % original sample frequency
  62. settings.do_resample = 100; % [] or sample frequency
  64. % baseline and zscoring settings
  65. settings.zscore_data4class = 1;
  66. settings.bs_correct = 1;
  67. settings.bs_period = [-.2 0]; % [-.5 -.1]
  68. settings.bs_trim = 0; % can be 0. amount of % to trim away when calculating the baseling
  71. % smoothing options
  72. settings.do_smoothdata = 1; % use matlabs smoothdata function for running average
  73. settings.smooth_win = .200; % .100, [] running average time window in seconds;
  75. % time of interest
  76. settings.TOI_train = [-.5 3];
  77. settings.timesteps_train = ((settings.ori_fsample/settings.do_resample)/settings.ori_fsample); % in s. if you want it to take less sample points multiply [e.g., ((settings.ori_fsample/settings.do_resample)/settings.ori_fsample)*2
  79. settings.TOI_test = [-1.2 1.2]; % time around ripple (I take this time window to get a proper estimate around the edges too. Only look at -1 to 1 later.
  80. settings.timesteps_test = ((settings.ori_fsample/settings.do_resample)/settings.ori_fsample); % in s
  82. settings.classifier = 'lda';
  83. settings.metric = 'auc';
  85. %% channel settings
  87. settings.channel_sel = 1; % 1 exclude hippo, 2 only hippo, 3 all channels
  89. %% settings pca
  91. settings.smooth_before_dim = 1; % smooth Nans before doing dim reduction
  92. settings.nu_time_points = 60; % time of sliding window in ms
  93. settings.prc_overlap = .9; % percentage overlap sliding window
  94. settings.decode_components = 1; % Decode the PCA components I picked.
  95. %% start for loop
  97. timeaxis = settings.TOI_test(1):settings.timesteps_test:settings.TOI_test(2);
  98. freqaxis = settings.TOI_train(1):settings.timesteps_train:settings.TOI_train(2);
  100. perf = cell(1,numel(subjects));
  101. channs = cell(1,numel(subjects));
  103. RT_all_subj_correct = cell(1,numel(subjects));
  104. RT_all_subj_incorrect = cell(1,numel(subjects));
  105. ripple_time_correct = cell(1,numel(subjects));
  106. ripple_time_incorrect = cell(1,numel(subjects));
  108. numWorkers = 8; %
  110. parpool('local', numWorkers);
  112. tic
  113. parfor isubject = 1:numel(subjects)
  115. fprintf('processing subject %01d/%02d\n',isubject,numel(subjects));
  117. %% LOAD data
  119. tmp = load([settings.data_dir,'eeg_session01_all_chan_nobadchan_cmntrim_artdet_',subjects{isubject}]);
  120. data = tmp.data;
  121. onsets_session = tmp.onsets_session;
  122. tmp = [];
  124. %% load channels and restrict to channels that are in hippocampus
  126. SubjectID = SubjectIDs{isubject};
  128. tmp = load([settings.base_path_castle,'ripple_project_publication_for_replication/templates/channels_hipp_ripples']);
  129. channels_hipp_ripples = tmp.channels_hipp_ripples;
  130. tmp = [];
  131. channels = channels_hipp_ripples(isubject,:);
  132. channels = channels(~cellfun(@isempty,channels));
  134. cfg = [];
  135. cfg.channel = channels;
  136. data = ft_selectdata(cfg, data);
  138. % If NaNs in recording, interpolate these.
  140. cfg = [];
  141. cfg.prewindow = 3;
  142. cfg.postwindow = 3;
  143. data = ft_interpolatenan(cfg,data);
  145. %% Find ripples
  147. [inData,alldat] = detect_ripples(data, settings);
  149. %% remove false positives from ripple data
  151. if settings.remove_falsepositives
  153. alldat = remove_false_positives(alldat);
  155. end
  157. %% delete co-occuring ripples
  159. if settings.remove_ripple_duplicates
  161. alldat = remove_ripple_duplicates(alldat);
  163. end
  166. %% Load data to realign based on cue onset encoding and based on ripples
  169. data_in = load([settings.data_dir,'eeg_session01_all_chan_nobadchan_cmntrim_artdet_',subjects{isubject}],'data','onsets_session');
  170. data = data_in.data;
  171. data_in = [];
  172. onsets = onsets_session - data.time{1}(1)*inData.fsample;
  173. data.sampleinfo = 1+data.sampleinfo - data.sampleinfo(1);
  174. data.time{1} = 1/data.fsample+data.time{1}-data.time{1}(1);
  176. %% channel selection for data
  178. % 1 exclude hippo, 2 only hippo, 3 all channels
  179. tmp = load([settings.data_dir_channels,'/channels_to_exclude_all_hipp_both_hem.mat']);
  180. channels_to_exclude_all_hipp = tmp.channels_to_exclude_all_hipp;
  181. tmp = load([settings.data_dir_channels,'/channels_hipp_ripples.mat']);
  182. channels_hipp_ripples = tmp.channels_hipp_ripples;
  184. cfg = [];
  185. cfg.channel = data.label;
  187. switch settings.channel_sel
  188. case 1
  189. cfg.channel = setdiff(setdiff([data.label],char(channels_to_exclude_all_hipp{isubject,:})),settings.scalp_channels);
  190. case 2
  191. cfg.channel = intersect(cellstr(setdiff([data.label],settings.scalp_channels)),char(channels_hipp_ripples{isubject,:}));
  192. case 3
  193. cfg.channel = setdiff(data.label,settings.scalp_channels);
  194. end
  196. %--- load channel info and coordinates
  197. [~,~,entries] = xlsread(fullfile(settings.anatomy_dir,'well01_ripples_ROIs_w_labels.xlsx'));
  199. subject_colums = entries(1,:);
  200. column_names = entries(2,:);
  201. these_labels = entries(3:end,strcmp(subject_colums,['s' SubjectIDs{isubject}]) & strcmp(column_names,'label'));
  202. these_x = cell2mat(entries(3:end,strcmp(subject_colums,['s' SubjectIDs{isubject}]) & strcmp(column_names,'x')));
  203. these_y = cell2mat(entries(3:end,strcmp(subject_colums,['s' SubjectIDs{isubject}]) & strcmp(column_names,'y')));
  204. these_z = cell2mat(entries(3:end,strcmp(subject_colums,['s' SubjectIDs{isubject}]) & strcmp(column_names,'z')));
  206. data = ft_selectdata(cfg, data);
  208. % write out info of retained channels
  210. for ichannel = 1:numel(data.label)
  212. idx = strcmp(data.label{ichannel},these_labels);
  214. channs{isubject}(ichannel).names = data.label{ichannel};
  215. channs{isubject}(ichannel).coords = [these_x(idx) these_y(idx) these_z(idx)];
  216. end
  218. %% Load subject file and change RT
  220. [numbers,strings] = xlsread([settings.base_path_castle,'well01_behavior_all.xls']);
  222. strings = strings(2:end,:);
  224. if isnan(numbers(1,1))
  225. numbers = numbers(2:end,:);
  226. end
  229. sel = find(strcmp(strings(:,2),SubjectID));
  230. sel = sel(1:numel(onsets));
  231. trls_enc = strcmp(strings(sel,4),'encoding');
  232. trls_ret = strcmp(strings(sel,4),'retrieval');
  234. Memory = cell(size(strings(sel,12))); % Convert Memory to a cell array of the same size
  235. Memory = strings(sel,12);
  236. Memory = Memory(trls_ret);
  237. idx_trial = find(trls_ret);
  238. Memory(:, 2) = num2cell(idx_trial); % Assign idx_trial to the second column
  241. RT = numbers(sel,11);
  242. RT(RT==-1 & trls_enc==1) = 3; % -1 no press in time - set to 3s at encoding
  244. if settings.full_enc_trial
  245. RT(trls_enc==1) = 3; % [optional] set all encoding to 3s
  246. end
  248. RT(RT==-1 & trls_ret==1) = 5; % -1 no press in time - set to 5s at encoding and 5s at retrieval
  250. trialinfo = [];
  252. for itrial = 1:numel(sel)
  254. trialinfo(itrial).SubjectID = strings(sel(itrial),2);
  255. trialinfo(itrial).RunNumber = numbers(sel(itrial),3);
  256. trialinfo(itrial).ExpPhase = strings(sel(itrial),4);
  257. trialinfo(itrial).TrialNumber = numbers(sel(itrial),5);
  258. trialinfo(itrial).EventNumber = itrial;
  259. trialinfo(itrial).BlockType = strings(sel(itrial),6);
  260. trialinfo(itrial).Subcat = strings(sel(itrial),7);
  261. trialinfo(itrial).Word = strings(sel(itrial),8);
  262. trialinfo(itrial).OldNew = strings(sel(itrial),9);
  263. trialinfo(itrial).Response = strings(sel(itrial),10);
  264. trialinfo(itrial).Memory = strings(sel(itrial),12);
  265. trialinfo(itrial).RT = RT(itrial);
  266. end
  268. sel = [];
  269. strings = [];
  271. %% create onset matrices (remove last 250ms to ensure ripple event in trial)
  273. onsetmat = [onsets; onsets+(RT'.*data.fsample)-(settings.time_to_excl_RT*data.fsample)]';
  275. %% Extract ripples (optional to select long and short duration ripples)
  277. trl_ripple = [];
  278. cnt = 0;
  280. for ichannel = 1:numel(alldat)
  282. evs = alldat{ichannel}.evtIndiv.maxTime;
  283. envSum = alldat{ichannel}.evtIndiv.envSum;
  285. ripple_dur = alldat{ichannel}.evtIndiv.duration;
  287. duration_sel = logical(ones(1,numel(ripple_dur)));
  289. evs = evs(duration_sel);
  290. envSum = envSum(duration_sel);
  292. %% for each detected ripple, find the corresponding trial
  294. for iripple = 1:numel(evs)
  296. this_event = evs(iripple) >= onsetmat(:,1) & evs(iripple) <= onsetmat(:,2);
  298. if any(this_event)
  299. cnt=cnt+1;
  301. trl_ripple(cnt,1) = find(this_event); % note down corresponding event number
  302. trl_ripple(cnt,2) = evs(iripple); % note down ripple sample
  303. trl_ripple(cnt,3) = (evs(iripple) - onsetmat(this_event,1))/data.fsample; % note down time of ripple in trial
  304. trl_ripple(cnt,4) = ichannel; % note down channel
  305. trl_ripple(cnt,5) = envSum(iripple);
  307. end
  308. end
  310. end
  312. trl_ripple = sortrows(trl_ripple,1);
  314. %% restrict ripple data to retrieval
  316. trl_ripple = trl_ripple(ismember(trl_ripple(:,1),find(trls_ret)),:);
  318. %% [optional] do surrogates by taking time of ripple from other trial (for all or only for correct trials)
  321. correct_mem = cell2mat(Memory(ismember(Memory(:,1),'SourceCorrect'),2));
  322. idx_correct = ismember(trl_ripple(:,1), correct_mem);
  323. trl_ripple_correct = trl_ripple(idx_correct,:);
  325. if settings.do_surrogates == 1 % 1 for all trials
  327. randtrials = circshift(1:size(trl_ripple,1),1);
  328. tmp = trl_ripple(:,2)-round(trl_ripple(:,3).*data.fsample); % find cue onset
  329. tmp = tmp+round(trl_ripple(randtrials,3).*data.fsample); % add another ripple's event time
  330. trl_ripple(:,2) = tmp;
  331. trl_ripple(:,3) = trl_ripple(randtrials,3);
  333. elseif settings.do_surrogates == 2 % 2 for only correct trials
  335. randtrials_correct = circshift(1:size(trl_ripple_correct,1),-1);
  336. tmp_correct = trl_ripple_correct(:,2) - round(trl_ripple_correct(:,3) .* data.fsample); % find cue onset
  337. tmp_correct = tmp_correct + round(trl_ripple_correct(randtrials_correct,3) .* data.fsample); % add shuffled ripple event time
  338. trl_ripple_correct(:,2) = tmp_correct;
  339. trl_ripple_correct(:,3) = trl_ripple_correct(randtrials_correct,3);
  341. trl_ripple(idx_correct,:) = trl_ripple_correct; % add to original structure with only correct trials swapped
  343. end
  345. %% Create trial structure around ripples
  347. pretrig = round(abs(settings.TOI_test(1)) * data.fsample); % enough time to baseline correct later
  348. posttrig = round(abs(settings.TOI_test(2)) * data.fsample);
  350. cfg = [];
  351. cfg.trl = [trl_ripple(:,2)-pretrig trl_ripple(:,2)+posttrig -pretrig*ones(size(trl_ripple,1),1)];
  353. data_ripples = ft_redefinetrial(cfg,data);
  355. % add trial info for each ripple trial, accounting for multiple ripples
  356. % per trial
  357. data_ripples.trialinfo = [];
  359. for itrial = 1:numel(data_ripples.trial)
  361. corresponding_trialinfo = find([trialinfo.EventNumber]==trl_ripple(itrial,1));
  363. trl_info = trialinfo(corresponding_trialinfo);
  364. trl_info.sample_ripple = trl_ripple(itrial,2);
  365. trl_info.time_ripple = trl_ripple(itrial,3);
  366. trl_info.channel = trl_ripple(itrial,4);
  367. trl_info.name_channel = {alldat{trl_ripple(itrial,4)}.evtIndiv.label};
  368. trl_info.envSum = trl_ripple(itrial,5);
  370. data_ripples.trialinfo{itrial,1} = trl_info;
  371. end
  373. %% [optional] if there are multiple ripples per trial - pick the one ripple with highest activity (captured in the summed envelope metric)
  375. if settings.solo_ripple
  377. tmp_trl_info = data_ripples.trialinfo;
  378. % find the trials in which there were more than one ripple
  379. trlinfo = cell2mat(data_ripples.trialinfo);
  380. ripple_trial = [trlinfo.EventNumber];
  381. envSum = [trlinfo.envSum];
  382. idx_unique = unique(ripple_trial);
  384. sel = [];
  385. counter = 1;
  387. for itrial = 1:numel(idx_unique)
  389. idx = find(ripple_trial==idx_unique(itrial));
  390. [~,max_effect] = max(envSum(idx));
  391. sel(counter) = idx(max_effect);
  393. counter = counter+1;
  394. end
  396. cfg = [];
  397. cfg.trials = sel;
  398. data_ripples = ft_selectdata(cfg, data_ripples);
  400. end
  402. %% [optional] pick ripples in a specific time window
  403. if any(settings.ripple_latency)
  405. trlinfo = cell2mat(data_ripples.trialinfo);
  406. sel = [trlinfo.time_ripple] > settings.ripple_latency(1) & [trlinfo.time_ripple] < settings.ripple_latency(2);
  408. cfg = [];
  409. cfg.trials = sel;
  410. data_ripples = ft_selectdata(cfg, data_ripples);
  412. end
  414. %% Realign trials based on stimulus onsets
  416. pretrig = round(abs(settings.TOI_train(1)) * data.fsample);
  417. posttrig = round(abs(settings.TOI_train(2)) * data.fsample);
  419. cfg = [];
  420. cfg.trl = [onsets'-pretrig onsets'+posttrig -pretrig*ones(numel(onsets),1)];
  422. data_stimuli = ft_redefinetrial(cfg,data);
  424. data_stimuli.trialinfo = {};
  425. for itrial = 1:numel(trialinfo)
  426. data_stimuli.trialinfo{itrial,1} = trialinfo(itrial);
  427. end
  429. %% [optional] resample
  431. if any(settings.do_resample)
  432. cfg = [];
  433. cfg.resamplefs = settings.do_resample;
  434. data_ripples = ft_resampledata(cfg, data_ripples);
  435. data_stimuli = ft_resampledata(cfg, data_stimuli);
  436. end
  438. %% Time-lock data
  440. cfg = [];
  441. cfg.keeptrials = 'yes';
  442. cfg.removemean = 'no';
  444. data_ripples = ft_timelockanalysis(cfg, data_ripples);
  445. data_stimuli = ft_timelockanalysis(cfg, data_stimuli);
  447. fsample = round(1/(data_ripples.time(2)-data_ripples.time(1)));
  450. %% [optional] Running mean filter
  452. if ~isempty(settings.smooth_win)
  453. data_ripples.trial = smoothdata(data_ripples.trial,3,'movmean',settings.smooth_win*fsample);
  454. data_stimuli.trial = smoothdata(data_stimuli.trial,3,'movmean',settings.smooth_win*fsample);
  455. end
  457. %% [optional] BL correct (us pre-stim baseline also for ripple-locked data)
  459. if settings.bs_correct == 1
  461. % training data
  462. bl_idx = nearest(data_stimuli.time,settings.bs_period(1)):nearest(data_stimuli.time,settings.bs_period(2));
  463. bldat = trimmean(data_stimuli.trial(:,:,bl_idx),settings.bs_trim,'round',3);
  464. blmat = repmat(bldat,[1 1 size(data_stimuli.trial,3)]);
  465. data_stimuli.trial = data_stimuli.trial - blmat;
  467. % testing data
  468. trlinfo = cell2mat(data_ripples.trialinfo);
  469. orig_events = [trlinfo.EventNumber];
  471. bl4ripples = nan(size(data_ripples.trial,1),size(bldat,2));
  473. for iripple = 1:size(data_ripples.trial,1)
  474. bl4ripples(iripple,:) = bldat(orig_events(iripple),:);
  475. end
  476. blmat = repmat(bl4ripples,[1 1 size(data_ripples.trial,3)]);
  478. data_ripples.trial = data_ripples.trial - blmat;
  480. end
  483. %% time-lock data
  485. cfg = [];
  486. cfg.keeptrials = 'yes';
  487. data_stimuli = ft_timelockanalysis(cfg,data_stimuli);
  488. data_ripples = ft_timelockanalysis(cfg,data_ripples);
  490. %% calculate time of ripple and RT for those trials
  492. ripple_tmp = [];
  493. RT_tmp = [];
  494. for itrial = 1:numel(data_ripples.trialinfo)
  496. ripple_tmp(itrial) = data_ripples.trialinfo{itrial, 1}.time_ripple;
  497. RT_tmp(itrial) = data_ripples.trialinfo{itrial, 1}.RT;
  499. end
  500. trlinfo = cell2mat(data_ripples.trialinfo);
  501. sel1 = ismember([trlinfo.ExpPhase],'retrieval') & ismember([trlinfo.Memory],'SourceCorrect');
  502. sel2 = ismember([trlinfo.ExpPhase],'retrieval') &(ismember([trlinfo.Memory],{'SourceIncorrect' ,'SourceDunno' 'dunno'}));
  504. RT_all_subj_correct{isubject} = RT_tmp(sel1);
  505. RT_all_subj_incorrect{isubject} = RT_tmp(sel2);
  506. ripple_time_correct{isubject} = ripple_tmp(sel1);
  507. ripple_time_incorrect{isubject} = ripple_tmp(sel2);
  510. %%
  512. %%%%%%%%%%%%%%%%%%%%%%
  513. %%%%%%%%% PCA %%%%%%%%
  514. %%%%%%%%%%%%%%%%%%%%%%
  516. trlinfo = cell2mat(data_stimuli.trialinfo);
  518. sel1 = ismember([trlinfo.ExpPhase],'encoding') & ismember([trlinfo.BlockType],'color');
  519. sel2 = ismember([trlinfo.ExpPhase],'encoding') & ismember([trlinfo.BlockType],'scene');
  521. dataToClassifyTraining = cat(1,data_stimuli.trial(sel1,:,:),data_stimuli.trial(sel2,:,:));
  522. clabelTraining = cat(1,1*ones(sum(sel1),1),2*ones(sum(sel2),1));
  523. samples_train = nearest(data_stimuli.time,settings.TOI_train(1)):settings.timesteps_train*fsample:nearest(data_stimuli.time,settings.TOI_train(2));
  525. dataToClassifyTraining = dataToClassifyTraining(:,:,samples_train);
  527. %% 1. category PCA [colours vs. scenes; coarse]
  529. trlinfo = cell2mat(data_ripples.trialinfo);
  531. %% PCA on the data to get the eigenvalues that explain 90% of the variance.
  532. % Do it with a sliding window of 60ms with 100Hz)
  536. for icond = 1:2
  538. if icond == 1
  540. sel1 = ismember([trlinfo.ExpPhase],'retrieval') & ismember([trlinfo.BlockType],'color') & ismember([trlinfo.Memory],'SourceCorrect');
  541. sel2 = ismember([trlinfo.ExpPhase],'retrieval') & ismember([trlinfo.BlockType],'scene') & ismember([trlinfo.Memory],'SourceCorrect') ;
  543. perf{isubject}.correct.trl_num_test = [sum(sel1) sum(sel2)];
  545. elseif icond == 2
  547. sel1 = ismember([trlinfo.ExpPhase],'retrieval') & ismember([trlinfo.BlockType],'color')& (ismember([trlinfo.Memory],{'SourceIncorrect' ,'SourceDunno' 'dunno'}));
  548. sel2 = ismember([trlinfo.ExpPhase],'retrieval') & ismember([trlinfo.BlockType],'scene')& (ismember([trlinfo.Memory],{'SourceIncorrect' ,'SourceDunno' 'dunno'}));
  550. perf{isubject} .incorrect.trl_num_test = [sum(sel1) sum(sel2)];
  552. end
  555. clabelTest = cat(1,1*ones(sum(sel1),1),2*ones(sum(sel2),1));
  556. dataToClassifyTest = cat(1,data_ripples.trial(sel1,:,:),data_ripples.trial(sel2,:,:));
  559. sumsel1 = sum(sel1);
  560. sumsel2 = sum(sel2);
  562. fs = 1/(data_ripples.time(2)-data_ripples.time(1));
  563. chunk_size = round(settings.nu_time_points/((1/fs)*1000)); % Size of each chunk
  564. overlap_size = round(chunk_size * settings.prc_overlap);
  565. num_chunks = floor((size(dataToClassifyTest, 3) - overlap_size) / (chunk_size - overlap_size));
  566. TOI_ripple = linspace(data_ripples.time(1), data_ripples.time(end),num_chunks);
  568. explained_variances = [];
  569. how_much_variance = [];
  570. accuracy = [];
  571. ripple_to_decode = [];
  572. for i = 1:num_chunks
  573. % Calculate the start and end indices of the current chunk
  574. start_idx = (i - 1) * (chunk_size - overlap_size) + 1;
  575. end_idx = start_idx + chunk_size - 1;
  577. % Extract data for the current chunk and reshape
  578. chunk_data = dataToClassifyTest(:, :, start_idx:end_idx);
  579. chunk_data = reshape(chunk_data, size(dataToClassifyTest, 1), []);
  581. if settings.smooth_before_dim == 1 % smooth NaNs through linear interpolation
  582. for ismooth = 1:size(chunk_data, 1)
  584. valid_indices = ~isnan(chunk_data(ismooth, :));
  585. chunk_data(ismooth, :) = interp1(find(valid_indices), chunk_data(ismooth, valid_indices), 1:size(chunk_data, 2), 'linear', 'extrap');
  586. end
  587. end
  589. % Perform PCA and compute explained variance
  590. [coefficients, ~, latent, ~, explained] = pca(chunk_data);
  592. % Compute explained variance
  593. explained_variance_pca = latent / sum(latent);
  595. % use a data-driven approach to get the first elbow point where
  596. % least variance is explained.
  597. curvature = diff(diff(explained_variance_pca));
  598. [~, elbow_index] = max(curvature);
  599. elbow_component = elbow_index + 1; % Add 1 because of diff operation
  602. explained_variances(i, :) = elbow_component;
  604. how_much_variance(i,:) = sum(explained_variance_pca(1:elbow_component));
  606. % Do PCA inverse transformation for later decoding
  607. if settings.decode_components == 1
  608. selected_components = coefficients(:,1:elbow_component);
  609. transformed_data = chunk_data * selected_components;
  610. reconstructed_data = transformed_data * selected_components';
  611. reconstructed_data = reshape(reconstructed_data,[size(dataToClassifyTest(:, :, start_idx:end_idx))]);
  612. ripple_to_decode(:,:,i) = nanmean(reconstructed_data,3); % take mean of those time points used in sliding window
  613. end
  614. end
  617. accuracy = explained_variances;
  620. if settings.decode_components == 1
  622. if settings.zscore_data4class
  623. dataToClassifyTraining = zscore(dataToClassifyTraining);
  624. ripple_to_decode = zscore(ripple_to_decode);
  625. end
  627. cfg = [];
  628. cfg.classifier = settings.classifier;
  629. cfg.metric = settings.metric;
  630. [accuracy_dec, ~] = mv_classify_timextime(cfg, dataToClassifyTraining, clabelTraining, ripple_to_decode, clabelTest);
  632. end
  634. if icond == 1
  635. perf{isubject}.correct.accuracy = accuracy;
  636. perf{isubject}.correct.exl_var = how_much_variance;
  637. if settings.decode_components == 1
  638. perf{isubject}.dec.correct.accuracy = accuracy_dec;
  639. end
  640. elseif icond == 2
  641. perf{isubject}.incorrect.accuracy = accuracy;
  642. perf{isubject}.incorrect.exl_var = how_much_variance;
  643. if settings.decode_components == 1
  644. perf{isubject}.dec.incorrect.accuracy = accuracy_dec;
  645. end
  646. end
  648. % add info
  649. perf{isubject}.channelcount_test = size(data_ripples.trial,2);
  650. perf{isubject}.channels_test = data_ripples.label;
  651. perf{isubject}.time_train = TOI_ripple;
  652. perf{isubject}.time_test = TOI_ripple;
  654. end
  656. % clearvars -except perf settings isubject subjects SubjectIDs RT_all_subj_correct RT_all_subj_incorrect
  657. end
  658. toc
  660. delete(gcp);
  662. return
  663. %% Stats and plots
  665. RT_correct = cellfun(@mean, RT_all_subj_correct);
  666. RT_max_correct = cellfun(@max, RT_all_subj_correct);
  667. RT_min_correct = cellfun(@min, RT_all_subj_correct);
  669. RT_incorrect = cellfun(@mean, RT_all_subj_incorrect);
  670. RT_max_incorrect = cellfun(@max, RT_all_subj_incorrect);
  671. RT_min_incorrect = cellfun(@min, RT_all_subj_incorrect);
  673. perf{1}.RT.correct = RT_all_subj_correct;
  674. perf{1}.RT.incorrect = RT_all_subj_incorrect;
  677. RT = RT_correct;
  679. ripple_time_mean = cellfun(@mean, ripple_time_correct);
  680. ripple_time_max = cellfun(@max, ripple_time_correct);
  681. ripple_time_min = cellfun(@min, ripple_time_correct);
  683. perf{1}.ripple_times.correct = ripple_time_correct;
  684. perf{1}.ripple_times.incorrect = ripple_time_incorrect;
  686. ripples_to_plot = [];
  687. rt_to_plot = [];
  688. explained_var_corr = [];
  689. explained_var_incorr = [];
  690. for participant = 1:numel(subjects)
  692. ripples_to_plot = [ripples_to_plot, ripple_time_correct{participant}];
  693. rt_to_plot = [rt_to_plot, RT_all_subj_correct{participant}];
  696. explained_var_corr(participant,:) = perf{1, participant}.correct.exl_var;
  697. explained_var_incorr(participant,:) = perf{1, participant}.incorrect.exl_var;
  698. end
  700. delay_ripple_rt = rt_to_plot-ripples_to_plot;
  702. figure;
  703. subplot(3,1,1)
  704. nhist(ripples_to_plot','proportion','color',settings.colour_scheme(8,:))
  705. title('Time of ripples')
  706. xlabel('time of ripples')
  707. ylabel('proportion')
  708. set(gca,'FontSize',14)
  709. set(gca,'TickDir','out')
  710. title(sprintf('Time of ripples, median = %.2fms',median(ripple_time_mean*1000)),'interpreter','none')
  711. subplot(3,1,2)
  712. nhist(rt_to_plot','proportion','color',settings.colour_scheme(8,:))
  713. title('Reaction time in trials of ripples')
  714. xlabel('Reaction time')
  715. ylabel('proportion')
  716. set(gca,'FontSize',14)
  717. set(gca,'TickDir','out')
  718. title(sprintf('Reaction time in trials of ripples, median = %.2fms',median(RT*1000)),'interpreter','none')
  719. subplot(3,1,3)
  720. nhist(delay_ripple_rt','proportion','color',settings.colour_scheme(8,:))
  721. xlabel('Delay ripple RT')
  722. ylabel('proportion')
  723. set(gca,'FontSize',14)
  724. set(gca,'TickDir','out')
  725. title(sprintf('Delay ripples RT, median = %.2fms',median(delay_ripple_rt*1000)),'interpreter','none')
  728. for participant = 1:numel(subjects)
  729. trl_num_test = perf{participant}.correct.trl_num_test;
  730. trl_correct(participant) = sum(trl_num_test);
  731. trl_num_test = perf{participant}.incorrect.trl_num_test;
  732. trl_incorrect(participant) = sum(trl_num_test);
  733. end
  735. data_nu_trl = {};
  736. data_nu_trl{1,1} = trl_correct;
  737. data_nu_trl{2,1} = trl_incorrect;
  739. [~,p_val,~,stats] = ttest(trl_correct,trl_incorrect)
  741. figure;
  742. h = rm_raincloud(data_nu_trl, [settings.colour_scheme(6,:)],0,'ks',[],settings.colour_scheme);
  744. h.p{1, 1}.FaceColor = settings.colour_scheme(1,:);
  745. h.s{1, 1}.MarkerFaceColor = settings.colour_scheme(1,:);
  746. h.m(1, 1).MarkerFaceColor = settings.colour_scheme(1,:);
  747. h.p{2, 1}.FaceColor = settings.colour_scheme(10,:);
  748. h.s{2, 1}.MarkerFaceColor = settings.colour_scheme(10,:);
  749. h.m(2, 1).MarkerFaceColor = settings.colour_scheme(10,:);
  751. hold on
  752. title(sprintf('number of trials for the two conditions, t-stat = %.2f', stats.tstat))
  753. set(gca,'TickDir','out')
  754. xlabel('number of trials')
  755. yticklabels({sprintf('incorrect %d',mean(trl_incorrect)), sprintf('correct %d',floor(mean(trl_correct)))})
  756. ylabel('condition')
  757. set(gca,'FontSize',20)
  758. axis tight
  760. [~,p_val,~,stats] = ttest(trl_correct,trl_incorrect)
  762. correct_incorrect = {};
  764. for isubject = 1:numel(subjects)
  766. correct_incorrect{1}.label = {'Channels'};
  767. correct_incorrect{1}.time = perf{1,1}.time_train;
  769. correct_incorrect{1}.individual(isubject,1,:) = perf{isubject}.correct.accuracy;
  770. % correct_incorrect{1}.individual(isubject,1,:) = explained_var_corr(isubject,:);
  771. correct_incorrect{1}.dimord = 'subj_chan_time';
  773. end
  775. correct_incorrect{1}.avg = squeeze(correct_incorrect{1}.individual);
  778. correct_incorrect{2} = correct_incorrect{1};
  781. for isubject = 1:numel(subjects)
  783. correct_incorrect{2}.individual(isubject,1,:) = perf{isubject}.incorrect.accuracy;
  784. % correct_incorrect{2}.individual(isubject,1,:) = explained_var_incorr(isubject,:);
  786. end
  788. correct_incorrect{2}.avg = squeeze(correct_incorrect{2}.individual);
  790. % decoding
  791. if settings.decode_components == 1
  794. correct = cell(1,numel(subjects));
  795. incorrect = cell(1,numel(subjects));
  796. correct_dec = [];
  797. incorrect_dec = [];
  799. for isubject = 1:numel(subjects)
  801. % correct
  802. correct{1,isubject} = struct;
  803. correct{1,isubject}.label = {'chan'};
  804. correct{1,isubject}.dimord = 'chan_freq_time';
  805. correct{1,isubject}.freq = freqaxis;
  806. correct{1,isubject}.time = perf{1, 1}.time_test;
  807. correct{1,isubject}.powspctrm(1,:,:) = perf{isubject}.dec.correct.accuracy;
  810. % incorrect
  811. incorrect{1,isubject} = correct{1,isubject};
  812. incorrect{1,isubject}.powspctrm(1,:,:) = perf{isubject}.dec.incorrect.accuracy;
  814. baseline_all{1,isubject} = correct{1,isubject};
  815. baseline_all{1,isubject}.powspctrm(1,:,:) = .5*ones(size(perf{isubject}.dec.correct.accuracy));
  817. correct_dec(isubject,:,:) = perf{isubject}.dec.correct.accuracy;
  818. incorrect_dec(isubject,:,:) = perf{isubject}.dec.incorrect.accuracy;
  819. end
  821. end
  824. %% FT stats (dimensionality)
  826. xlimits = nearest(correct_incorrect{1, 1}.time, -1):nearest(correct_incorrect{1, 1}.time, 1);
  827. xlimits = correct_incorrect{1, 1}.time(xlimits);
  829. cfg = [];
  830. cfg.latency = [xlimits(1) xlimits(end)];
  832. cfg.channel = 'all';
  833. cfg.statistic = 'depsamplesT';
  834. cfg.method = 'montecarlo'; % 'montecarlo' 'analytic';
  835. cfg.correctm = 'cluster'; % 'no', cluster;
  836. cfg.alpha = .05;
  837. cfg.clusteralpha = .05;
  838. cfg.tail = 0;
  839. cfg.correcttail = 'alpha'; % alpha prob no
  840. cfg.neighbours = [];
  841. cfg.minnbchan = 0;
  842. cfg.avgovertime = 'no'; % 'no' 'yes'
  843. cfg.avgoverchan = 'no';
  844. cfg.computecritval = 'yes';
  846. cfg.numrandomization = 'all';%'all';
  848. cfg.clusterstatistic = 'maxsum'; % 'maxsum', 'maxsize', 'wcm'
  849. cfg.clustertail = cfg.tail;
  850. cfg.parameter = 'individual';
  852. nSub = size(correct_incorrect{1, 1}.individual ,1);
  853. % set up design matrix
  854. design = zeros(2,2*nSub);
  855. for i = 1:nSub
  856. design(1,i) = i;
  857. end
  858. for i = 1:nSub
  859. design(1,nSub+i) = i;
  860. end
  861. design(2,1:nSub) = 1;
  862. design(2,nSub+1:2*nSub) = 2;
  864. cfg.design = design;
  865. cfg.uvar = 1;
  866. cfg.ivar = 2;
  868. % run stats
  869. [Fieldtripstats] = ft_timelockstatistics(cfg, correct_incorrect{:});
  870. length(find(Fieldtripstats.mask))
  873. %% plot significant vals
  875. d = squeeze(correct_incorrect{1}.individual);
  877. m = nanmean(d);
  878. s = nanstd(d)./sqrt(size(d,1));
  880. figure;
  881. boundedline(perf{1,1}.time_train,m,s,'b');
  882. plot(perf{1,1}.time_train,m,'k','linewidth',2);
  883. hold on
  885. d = squeeze(correct_incorrect{2}.individual);
  887. m = nanmean(d);
  888. s = nanstd(d)./sqrt(size(d,1));
  891. boundedline(perf{1,1}.time_train,m,s,'r');
  892. plot(perf{1,1}.time_train,m,'k','linewidth',2);
  893. hold on
  895. stats_time = nearest(correct_incorrect{1}.time,cfg.latency(1)):nearest(correct_incorrect{1}.time,cfg.latency(2));
  897. sigline = nan(1,numel(correct_incorrect{1}.time));
  898. % sigline(stats_time(Fieldtripstats.mask==1)) = m(stats_time(Fieldtripstats.mask==1));
  899. sigline(stats_time(Fieldtripstats.mask==1)) = 2;
  901. plot(correct_incorrect{1}.time,sigline,'r','linewidth',4);
  903. set(gca,'FontSize',16,'FontName','Arial')
  904. xlabel('ripple time (s)')
  905. ylabel('dimensionality difference')
  906. set(gca,'TickDir','out')
  908. axis tight
  909. vline(0)
  912. xlim([cfg.latency(1), cfg.latency(end)])
  917. m_exp_corr = mean(explained_var_corr(:,sigline==2),2);
  918. m_exp_incorr = mean(explained_var_incorr(:,sigline==2),2);
  920. mean(m_exp_corr)
  921. mean(m_exp_incorr)
  923. std(m_exp_corr)
  924. std(m_exp_incorr)
  927. [~,p_val_expl_var,~,stat] = ttest(m_exp_corr,m_exp_incorr)
  929. %% relate sigline to reaction time on a group level
  930. d = squeeze(correct_incorrect{1}.individual);
  932. RT = RT_correct;
  934. dimensionality_change = mean(d(:,sigline==2),2);
  936. figure;
  937. scatter(RT, dimensionality_change, 'filled', 'MarkerFaceColor', '#0072BD');
  938. xlabel('RT', 'FontSize', 12);
  939. ylabel('Dimensionality Change', 'FontSize', 12);
  940. title('Correlation', 'FontSize', 14);
  941. grid on;
  942. box on;
  943. hold on;
  945. % Fit a linear regression line
  946. p = polyfit(RT, dimensionality_change, 1);
  947. f = polyval(p, RT);
  948. plot(RT, f, 'r-', 'LineWidth', 1.5);
  951. % Add legend
  952. legend('Data', 'Linear Fit', 'Location', 'best');
  954. % Customize the plot appearance
  955. set(gca, 'FontSize', 10); % Set font size for axis labels
  956. set(gcf, 'Color', 'w'); % Set background color of the figure to white
  959. yfit = polyval(p, RT);
  960. yresid = dimensionality_change - yfit;
  961. SSresid = sum(yresid.^2);
  962. SStotal = (length(dimensionality_change)-1) * var(dimensionality_change);
  963. rsq = 1 - SSresid/SStotal;
  964. disp(['R-squared: ', num2str(rsq)]);
  966. [rho_rt_dim, p_rt_dim] = corr(RT',dimensionality_change, 'type', 'Spearman');
  968. %% correlate RT across time
  969. data_to_correlate = squeeze(correct_incorrect{1}.individual);
  970. RT = RT_correct;
  972. fs = 1/(correct_incorrect{1, 1}.time(2)-correct_incorrect{1, 1}.time(1));
  973. chunk_size = round(100/((1/fs)*1000)); % Size of each chunk
  974. overlap_size = round(chunk_size * settings.prc_overlap);
  975. num_chunks = floor((size(correct_incorrect{1}.individual, 3) - overlap_size) / (chunk_size - overlap_size));
  976. TOI_corr = linspace(correct_incorrect{1, 1}.time(1), correct_incorrect{1, 1}.time(end),num_chunks);
  978. rho_across_time = [];
  979. p_across_time = [];
  980. for itime = 1:num_chunks
  981. start_idx = (itime - 1) * (chunk_size - overlap_size) + 1;
  982. end_idx = start_idx + chunk_size - 1;
  985. [rho_tmp, p_tmp] = corr(mean(data_to_correlate(:,start_idx:end_idx),2), RT','type', 'spearman');
  986. rho_across_time(itime) = rho_tmp;
  987. p_across_time(itime) = p_tmp;
  988. end
  990. rho_across_time_perm = [];
  991. p_across_time_perm = [];
  992. for nu_perm = 1:500
  993. rand_rt = randperm(12);
  995. for itime = 1:num_chunks
  996. start_idx = (itime - 1) * (chunk_size - overlap_size) + 1;
  997. end_idx = start_idx + chunk_size - 1;
  999. [rho_tmp, p_tmp] = corr(mean(data_to_correlate(rand_rt,start_idx:end_idx),2), RT','type', 'spearman');
  1000. rho_across_time_perm(nu_perm,itime) = rho_tmp;
  1001. p_across_time_perm(nu_perm,itime) = p_tmp;
  1002. end
  1003. end
  1005. xlimits_ind = nearest(TOI_corr, -1):nearest(TOI_corr, 1);
  1008. alpha = 0.05;
  1009. z_threshold = norminv(1 - alpha/2);
  1011. zvalue_rho = (rho_across_time(xlimits_ind)-(mean(rho_across_time_perm(:,xlimits_ind))))./std(rho_across_time_perm(:,xlimits_ind));
  1013. xlimits = TOI_corr(xlimits_ind);
  1014. plot(xlimits, zvalue_rho,'linewidth', 3)
  1016. hold on
  1017. plot(xlimits, z_threshold * ones(size(xlimits)), 'r--'); % positive threshold
  1018. plot(xlimits, -z_threshold * ones(size(xlimits)), 'r--'); % negative threshold
  1019. below_threshold = zvalue_rho < -z_threshold;
  1021. scatter(xlimits(below_threshold), zvalue_rho(below_threshold), 'r', 'filled', 'MarkerFaceAlpha', 0.5)
  1023. ylim([-3 3])
  1024. xlim([-1 1])
  1025. set(gca,'FontSize',16)
  1026. set(gca,'TickDir','out')
  1028. hold off
  1030. xlabel('Ripple time (sec)')
  1031. ylabel('Z-transformed correlation')
  1032. title('Z-value of correlation Across Time')
  1034. legend('Z-value', 'Positive Threshold', 'Negative Threshold', 'Location', 'NorthEast')
  1037. %% plot t line
  1039. figure
  1040. plot(Fieldtripstats.time,Fieldtripstats.stat,'k','linewidth',2)
  1041. hold on
  1042. sigline05 = nan(1,numel(Fieldtripstats.prob));
  1043. sigline01 = nan(1,numel(Fieldtripstats.prob));
  1044. p05 = Fieldtripstats.prob < .05;
  1045. p01 = Fieldtripstats.prob < .01;
  1047. sigline05(p05) = Fieldtripstats.stat(p05);
  1048. sigline01(p01) = Fieldtripstats.stat(p01);
  1050. plot(Fieldtripstats.time,sigline05,'r','linewidth',5)
  1051. plot(Fieldtripstats.time,sigline01,'y','linewidth',2)
  1053. vline(0)
  1054. hline(0)
  1056. set(gca,'FontSize',8)
  1057. xlabel('time (sec)')
  1058. set(gca,'TickDir','out')
  1060. title('dimensionality');
  1061. xlim([cfg.latency(1), cfg.latency(end)])
  1065. %% correlate fine-grained decoding with dimensionality
  1067. % FT stats (dimensionality)
  1069. xlimits = nearest(correct_incorrect{1, 1}.time, -1):nearest(correct_incorrect{1, 1}.time, 1);
  1070. xlimits = correct_incorrect{1, 1}.time(xlimits);
  1072. cfg = [];
  1073. cfg.latency = [xlimits(1) xlimits(end)];
  1075. cfg.channel = 'all';
  1076. cfg.statistic = 'depsamplesT';
  1077. cfg.method = 'montecarlo'; % 'montecarlo' 'analytic';
  1078. cfg.correctm = 'cluster'; % 'no', cluster;
  1079. cfg.alpha = .05;
  1080. cfg.clusteralpha = .05;
  1081. cfg.tail = 0;
  1082. cfg.correcttail = 'alpha'; % alpha prob no
  1083. cfg.neighbours = [];
  1084. cfg.minnbchan = 0;
  1085. cfg.avgovertime = 'no'; % 'no' 'yes'
  1086. cfg.avgoverchan = 'no';
  1087. cfg.computecritval = 'yes';
  1089. cfg.numrandomization = 500;%'all';
  1091. cfg.clusterstatistic = 'maxsum'; % 'maxsum', 'maxsize', 'wcm'
  1092. cfg.clustertail = cfg.tail;
  1093. cfg.parameter = 'individual';
  1095. nSub = size(correct_incorrect{1, 1}.individual ,1);
  1096. % set up design matrix
  1097. design = zeros(2,2*nSub);
  1098. for i = 1:nSub
  1099. design(1,i) = i;
  1100. end
  1101. for i = 1:nSub
  1102. design(1,nSub+i) = i;
  1103. end
  1104. design(2,1:nSub) = 1;
  1105. design(2,nSub+1:2*nSub) = 2;
  1107. cfg.design = design;
  1108. cfg.uvar = 1;
  1109. cfg.ivar = 2;
  1111. % run stats
  1112. [Fieldtripstats] = ft_timelockstatistics(cfg, correct_incorrect{:});
  1113. length(find(Fieldtripstats.mask))
  1116. %% plot significant vals
  1118. d = squeeze(correct_incorrect{1}.individual-correct_incorrect{2}.individual);
  1120. m = nanmean(d);
  1121. s = nanstd(d)./sqrt(size(d,1));
  1123. figure;
  1124. boundedline(perf{1,1}.time_train,m,s,'k');
  1125. plot(perf{1,1}.time_train,m,'k','linewidth',2);
  1126. hold on
  1128. stats_time = nearest(correct_incorrect{1}.time,cfg.latency(1)):nearest(correct_incorrect{1}.time,cfg.latency(2));
  1130. sigline = nan(1,numel(correct_incorrect{1}.time));
  1131. % sigline(stats_time(Fieldtripstats.mask==1)) = m(stats_time(Fieldtripstats.mask==1));
  1132. sigline(stats_time(Fieldtripstats.mask==1)) = 0;
  1134. plot(correct_incorrect{1}.time,sigline,'r','linewidth',4);
  1136. set(gca,'FontSize',16,'FontName','Arial')
  1137. xlabel('ripple time (s)')
  1138. ylabel('dimensionality difference')
  1139. set(gca,'TickDir','out')
  1141. axis tight
  1142. vline(0)
  1143. hline(0)
  1145. xlim([cfg.latency(1), cfg.latency(end)])
  1147. % load original data
  1148. load('mask_t_vals_fine_enc_ripple.mat')
  1150. mask_t_vals_fine = squeeze(mask_t_vals_fine_enc_ripple);
  1152. % Load data from decoding analysis
  1153. load('correct_dec_fine_all.mat')
  1154. load('incorrect_dec_fine_all.mat')
  1156. correct_dec_all = correct_dec_fine_all;
  1157. incorrect_dec_all = incorrect_dec_fine_all;
  1159. time_dec = linspace(-.5,3,351);
  1160. idx_time = nearest(time_dec,-.2):nearest(time_dec,3);
  1162. correct_dec_all = correct_dec_all(:,idx_time,:);
  1164. incorrect_dec_all = incorrect_dec_all(:,idx_time,:);
  1166. for isubject = 1:size(settings.subjects,1)
  1168. tmp_sub = [];
  1169. tmp_sub = squeeze(correct_dec_all(isubject, :,:));
  1170. tmp_sub(~mask_t_vals_fine) = NaN;
  1171. correct_dec_all(isubject, :,:) = tmp_sub;
  1173. tmp_sub = [];
  1174. tmp_sub = squeeze(incorrect_dec_all(isubject, :,:));
  1175. tmp_sub(~mask_t_vals_fine) = NaN;
  1176. incorrect_dec_all(isubject, :,:) = tmp_sub;
  1177. end
  1179. correct_dec_mean = nanmean(nanmean(correct_dec_all,3),2);
  1180. incorrect_dec_mean = nanmean(nanmean(incorrect_dec_all,3),2);
  1182. dimensionality_change = mean(d(:,sigline==0),2);
  1184. [r, p] = corr(dimensionality_change,correct_dec_mean-incorrect_dec_mean,'tail','right');
  1186. % Number of data points
  1187. n = numel(dimensionality_change);
  1189. % Convert r to Fisher's z-score
  1190. z = atanh(r);
  1192. effect_size = z * sqrt(n - 3);
  1194. disp(['Effect size (Cohen''s d): ' num2str(effect_size)]);
  1196. %% visualise included channels
  1198. % elec_size = 15;
  1199. % transp = 0.25;
  1200. % extracolor = .5;
  1201. %
  1202. % epos_all = [];
  1203. %
  1204. % for isubject = 1:numel(subjects)
  1205. % epos_all = [epos_all;cell2mat({channs{isubject}.coords}')];
  1206. % end
  1207. %
  1208. % colvec = ones(1,size(epos_all,1));
  1209. %
  1210. % views = [-90 0;0 0;180 -90];
  1211. %
  1212. % for iview=1:size(views,1)
  1213. % figure
  1214. % plot_ecog(colvec, ...
  1215. % fullfile(ftpath,'template/anatomy/'),...
  1216. % epos_all,[-max(colvec) max(colvec)+extracolor], transp, views(iview,:), elec_size);
  1217. % colorbar off
  1218. % end
  1223. %%

ripple_stage_04a_dimensionality_final.m at commit 1d7307a, no license · at the source

Overview

Authors: Casper Kerrén1, Sebastian Michelmann2, Christian F. Doeller1,3
ORCID iDs: Casper Kerrén, Sebastian Michelmann, Christian F. Doeller
  1. Max Planck Institute for Human Cognitive and Brain Sciences,Leipzig, Germany
  2. Department of Psychology, New York University,New York, NY USA
  3. Kavli Institute for Systems Neuroscience, Centre for Neural Computation, Egil and Pauline Braathen and Fred Kavli Centre for Cortical Microcircuits, Jebsen Centre for Alzheimer’s Disease, NTNU Norwegian University of Science and Technology,Trondheim, Norway
Institutions: Max Planck Institute for Human Cognitive and Brain Sciences (Germany); New York University (United States); Kavli Institute for Systems Neuroscience (Norway)
Journal: Nature communications, volume 17, issue 1, article 6677
Dates: received 22 September 2025; accepted 26 June 2026; published online 20 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41467-026-75345-6 · PMID 42476972 · PMCID PMC13385388 · OpenAlex W4410097547
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.21203/rs.3.rs-6512178/v1
Status: code verified
Categories: EEG (modality), human (organism), epilepsy (population), cognitive (subfield)
Methods: Spectral & time-frequency, Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Evoked potentials, Single-unit activity, calcium imaging, Physiology & signal measures
Keywords: Learning and memory, Hippocampus
MeSH: Cerebral Cortex*, Hippocampus*, Mental Recall*, Drug Resistant Epilepsy, Electroencephalography, Female, Humans, Male, Memory, Episodic, Reaction Time (* major topic)
Topic: Memory and Neural Mechanisms (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Deutsche Forschungsgemeinschaft (437219953); Max-Planck-Gesellschaft
Citations: not cited yet (Europe PMC); 105 references in the paper

Abstract

How are past experiences reconstructed from memory? Learning is thought to compress external inputs into low-dimensional hippocampal representations, later expanded into high-dimensional cortical activity during recall. Hippocampal ripples, brief high-frequency bursts linked to retrieval, may initiate this expansion. Analysing intracranial EEG data from patients with pharmacoresistant epilepsy during an episodic memory task, we found that cortical dimensionality increased following ripple events during correct, but not incorrect, retrieval. This expansion correlated with faster reaction times and reinstatement of the target association. Crucially, hippocampal theta and cortical gamma phase-amplitude coupling emerged after ripples but before cortical expansion, suggesting a mechanism for ripple-driven communication. Ripple events also marked the separation of task-relevant variables in cortical state space, revealing how hippocampal output reshapes the geometry of memory representations to support successful recall.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

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

kerrencasper/Hippocampal-ripples-initiate-cortical-dimensionality-expansion-for-memory-retrieval

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 1d7307a6037aa27e50b7b489b04f79bb0aa77aa9, 16 April 2025
Languages: MATLAB (18)
Size: 60 files, 18 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: FieldTrip (17 files), Statistics and Machine Learning Toolbox (15 files), Parallel Computing Toolbox (10 files), boundedline (9 files), SPM (4 files), Image Processing Toolbox (2 files), Signal Processing Toolbox (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
19 files
At the source: github.com/kerrencasper/Hippocampal-ripples-initiate-cortical-dimensionality-expansion-for-memory-retrieval.git

Code availability

The code that supports the conclusions of this study is available at GitHub (https://github.com/kerrencasper/Hippocampal-ripples-initiate-cortical-dimensionality-expansion-for-memory-retrieval.git).

Reproduced under the paper's license (CC BY), from the paper cited above.

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;
  • 18 scripts, each with its path and the digest of its content;
  • 6 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

Datasets cited

Data availability

The data that support the conclusions of this study are available at Zenodo (10.5281/zenodo.18490239).

Reproduced under the paper's license (CC BY), from the paper cited above.

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 2, 28 September 2026

  • Funding: added Deutsche Forschungsgemeinschaft: 437219953; Max-Planck-Gesellschaft

Version 1, 27 September 2026: the first record

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

Cite

This paper

Kerrén, C., Michelmann, S., & Doeller, C. F. (2026). Hippocampal ripples initiate cortical dimensionality expansion for memory retrieval. Nature communications, 17(1), 6677. https://doi.org/10.1038/s41467-026-75345-6

BibTeX

@article{kerren2026hippocampal,
author = {Kerrén, Casper and Michelmann, Sebastian and Doeller, Christian F.},
title = {{Hippocampal ripples initiate cortical dimensionality expansion for memory retrieval}},
journal = {Nature communications},
year = {2026},
month = jul,
volume = {17},
number = {1},
pages = {6677},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-75345-6},
url = {https://doi.org/10.1038/s41467-026-75345-6},
pmid = {42476972},
pmcid = {PMC13385388}
}

RIS

TY - JOUR
AU - Kerrén, Casper
AU - Michelmann, Sebastian
AU - Doeller, Christian F.
TI - Hippocampal ripples initiate cortical dimensionality expansion for memory retrieval
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/07/20
VL - 17
IS - 1
SP - 6677
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-75345-6
UR - https://doi.org/10.1038/s41467-026-75345-6
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-75345-6",
"type": "article-journal",
"title": "Hippocampal ripples initiate cortical dimensionality expansion for memory retrieval",
"container-title": "Nature communications",
"author": [
{
"family": "Kerrén",
"given": "Casper"
},
{
"family": "Michelmann",
"given": "Sebastian"
},
{
"family": "Doeller",
"given": "Christian F."
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "6677",
"DOI": "10.1038/s41467-026-75345-6",
"PMID": "42476972",
"PMCID": "PMC13385388",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-75345-6",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
20
]
]
}
}

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.7554/elife.100642 [code]
Disrupted hippocampal theta-gamma coupling and spike-field coherence following experimental traumatic brain injury.
Journal: eLife
In common: Parallel Computing Toolbox, Image Processing Toolbox, Signal Processing Toolbox, 11 references
[2] doi:10.1038/s41593-026-02345-6 [code]
Human hippocampal ripples tune cortical responses based on predicted uncertainty.
Journal: Nature neuroscience
In common: boundedline, FieldTrip, Image Processing Toolbox, 2 other tools, cognitive, 9 references
[3] doi:10.1093/braincomms/fcag041
Exercise enhances hippocampal-cortical ripple interactions in the human brain.
Journal: Brain communications
In common: 13 references
[4] doi:10.7554/elife.110795 [code]
REM sleep prefrontal high-frequency oscillation chains mediate distinct cortical - hippocampal reactivation patterns compared to NREM sleep.
Journal: eLife
In common: boundedline, FieldTrip, Image Processing Toolbox, 2 other tools, 8 references
[5] doi:10.1038/s41467-026-70633-7 [code]
Global coincident bursts of high frequency oscillations across the human cortex coordinate large-scale memory processing.
Journal: Nature communications
In common: epilepsy, cognitive, 10 references
[6] doi:10.1162/imag.a.1218 [code]
Reliability and signal comparison of OPM-MEG, fMRI &amp; iEEG in a repeated movie viewing paradigm.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: FieldTrip, Image Processing Toolbox, Statistics and Machine Learning Toolbox, 5 references, author Sebastian Michelmann
[7] doi:10.1016/j.celrep.2026.117845 [code]
Hippocampal skill memory expansion drives online performance dynamics during skill learning.
Journal: Cell reports
In common: Image Processing Toolbox, cognitive, 9 references
[8] doi:10.7554/elife.108023 [code]
Challenges in replay detection by TDLM in post-encoding resting state.
Journal: eLife
In common: Statistics and Machine Learning Toolbox, 9 references
[9] doi:10.1038/s42003-025-08618-3 [code]
Physical activity simultaneously improves working memory and ripple-spindle coupling
Journal: n/a
In common: Statistics and Machine Learning Toolbox, EEG, cognitive, 8 references
[10] doi:10.1016/j.celrep.2026.117646 [code]
Medial entorhinal-hippocampal desynchronization parallels the emergence of memory impairment in a mouse model of Alzheimer's disease pathology.
Journal: Cell reports
In common: boundedline, FieldTrip, Parallel Computing Toolbox, 3 other tools, 3 references

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.