OSCR

Challenges in replay detection by TDLM in post-encoding resting state.

Code ↔ Paper

19 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 19 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] § Materials and methods › Material and methods study I: resting state analysis › MEG acquisition and pre-processing ↔ meg_tools.py, lines 328–476 · score 0.90 · find_bads_ecg, find_bads_eog, MNE, emg, muscle, picard
  2. [2] § Materials and methods › Material and methods study II: simulation › Replay density ↔ 2_run_study1.py, lines 98–144 · score 0.84 · 50–100 %, refractory period, linearly, reactivation, memory, 80 min
  3. [3] § Materials and methods › Material and methods study I: resting state analysis › Decoding framework and training ↔ 2_run_study1.py, lines 246–294 · score 0.81 · L1 regularization, 150–250 ms, scikit-learn, decoding accuracy, fold, 150 ms
  4. [4] § Materials and methods › Material and methods study II: simulation › Replay density ↔ 3_run_study2.py, lines 101–102 · score 0.75 · 50–100 %, linearly, events, Simulation
  5. [5] § Materials and methods › Material and methods study II: simulation › Synthetic simulation ↔ 7_run_discriminability_analysis.py, lines 307–348 · score 0.71 · logistic regression, classifier patterns, synthetic simulation, L1, noise, discriminability
  6. [6] § Materials and methods › Material and methods study I: resting state analysis › Sequential replay analysis ↔ tdlm/core.py, lines 136–220 · score 0.71 · auto correlations, transition matrices, rows, shifted, signal, GLM
  7. [7] § Materials and methods › Material and methods study II: simulation › Synthetic simulation ↔ 7_run_discriminability_analysis.py, lines 1–35 · score 0.69 · Wasserstein distance, pattern strength, classifier probabilities, separability, discriminability, ratio
  8. [8] § Materials and methods › Material and methods study II: simulation › Synthetic simulation ↔ 6_run_synthetic_simulation.py, lines 72–139 · score 0.68 · classifier patterns, synthetic simulation, replay events, inversely, noise, predicted
  9. [9] § Results › Previous synthetic simulation overestimates TDLM sensitivity ↔ 7_run_discriminability_analysis.py, lines 1–35 · score 0.65 · Wasserstein distance, synthetic simulation, Classifier probabilities, separability, quantify, discriminability
  10. [10] § Results › Previous synthetic simulation overestimates TDLM sensitivity ↔ 7_run_discriminability_analysis.py, lines 497–518 · score 0.61 · Wasserstein distance, pattern discriminability, pattern strength, mapping, synthetic, simulation
  11. [11] § Results › No sequential replay during resting state ↔ 2_run_study1.py, lines 625–655 · score 0.57 · cluster permutation, sign flip, post learning, backward, pre
  12. [12] § Results › No sequential replay during resting state ↔ 5_run_revision1.py, lines 212–232 · score 0.54 · sign flip permutation, fwd bkw, Post Learning, backwards sequenceness
  13. [13] § Results › No sequential replay during resting state ↔ 2_run_study1.py, lines 625–655 · score 0.54 · cluster permutation, sign flip, post learning, backward, pre
  14. [14] § Materials and methods › Material and methods study II: simulation › Simulating replay by inserting localizer data into the control resting state ↔ tdlm/utils.py, lines 328–396 · score 0.52 · reactivation events, replay event, chosen, zero, positions, onset
  15. [15] § Materials and methods › Material and methods study I: resting state analysis › MEG acquisition and pre-processing ↔ 3_run_study2.py, lines 134–160 · score 0.51 · eyes closed, decoding accuracy, ICA, Autoreject, onset, pre
  16. [16] § Materials and methods › Material and methods study I: resting state analysis › MEG acquisition and pre-processing ↔ 4_run_supplement.py, lines 129–156 · score 0.51 · eyes closed, decoding accuracy, ICA, Autoreject, onset, pre
  17. [17] § Materials and methods › Material and methods study I: resting state analysis › MEG acquisition and pre-processing ↔ meg_tools.py, lines 328–476 · score 0.51 · EOG, ECG, trigger, channels, filter, signal
  18. [18] § Materials and methods › Material and methods study I: resting state analysis › Sign-flip-permutation test ↔ MATLAB/uperms.m, the whole file · a weak match · score 0.51 · random subset, inconsistent, zero, Error, permutations
  19. [19] § Materials and methods › Material and methods study I: resting state analysis › Decoding framework and training ↔ 6_run_synthetic_simulation.py, lines 38–40 · score 0.50 · negative class, logistic regression, ratio, decoding

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

Python · 1,024 lines · 40 KB · GPL-3.0 · 4 matches

  1. # -*- coding: utf-8 -*-
  2. """
  3. Created on Mon Apr 15 11:08:58 2024
  5. @author: simon.kern
  6. """
  7. import os
  8. import settings
  9. import tdlm
  10. import utils
  11. from tqdm import tqdm
  12. import time
  13. import random
  14. import numpy as np
  15. import matplotlib.pyplot as plt
  16. import seaborn as sns
  17. import pandas as pd
  18. import mne
  19. from meg_utils import plotting, decoding, misc
  21. from mne.stats import permutation_cluster_1samp_test
  22. import joblib
  23. from joblib import Parallel, delayed
  24. from settings import results_dir, cache_dir
  25. from load_funcs import load_localizers_seq12, load_neg_x_before_audio_onset
  26. from load_funcs import load_RS1, load_RS2
  27. from utils import get_best_timepoint, get_performance, load_pkl_pandas
  28. from utils import plot_correlation, zscore_multiaxis
  29. from scipy.stats import ttest_rel
  30. from scipy import stats
  31. from statsmodels.stats.multitest import multipletests
  32. import warnings
  33. np.random.seed(0) # just for safety
  35. # plotting for paper default settings
  36. plt.rc("font", size=14)
  37. meanprops = {
  38. "marker": "o",
  39. "markerfacecolor": "black",
  40. "markeredgecolor": "white",
  41. "markersize": "10",
  42. }
  44. warnings.filterwarnings(
  45. "ignore",
  46. message="Mean of empty slice",
  47. category=RuntimeWarning,
  48. )
  49. warnings.filterwarnings(
  50. "ignore",
  51. message="invalid value encountered in divide",
  52. category=RuntimeWarning,
  53. )
  55. utils.lowpriority() # make sure not to clog CPU on multi-user systems
  56. #%% Settings
  58. # this is the classifier that will be used
  59. C = 9.1 # determined previously via cross-validation, see below
  60. clf = decoding.LogisticRegressionOvaNegX(C=C, penalty="l1", neg_x_ratio=2, rng=0)
  62. # list the files
  63. files = utils.list_files(settings.data_dir, patterns=["*DSMR*"])
  64. subjects = [f"DSMR{subj}" for subj in sorted(set(map(utils.get_id, files)))]
  66. min_acc = 0.3 # minimum accuracy of decoders to include subjects
  67. min_perf = 0.5 # minimum memory performance to include subjects
  68. sfreq = 100 # downsample to this frequency. Changing is not supported.
  69. ms_per_point = 10 # ms per sample point
  71. bands = settings.bands_HP # only use HP filter
  73. baseline = None # no baseline correction of localizer
  74. uid = hex(random.getrandbits(128))[2:10] # get random UID for saving log files
  75. date = time.strftime("%Y-%m-%d")
  76. times = np.arange(-100, 510, 10)
  77. final_calculation = True # this can be set to use the leftout data of the RS
  79. proba_norm = 'lambda x: x/x.mean(0)'
  80. # proba_norm = 'lambda x:x'
  82. # static colors for plotting
  83. palette = sns.color_palette()
  84. c_fwd = palette[1]
  85. c_bkw = palette[2]
  87. #%% define TDLM parameters
  88. n_shuf = 1000 # do 1000 permutations
  89. max_lag = 30 # 500 ms time lag maximum
  90. alpha_freq = 10 # assume 10 Hz alpha freq
  91. # can be time or subject, either each permutation is zscored or all
  92. # permutations at the subject level, see norm_func
  93. zscore_axes = -1
  95. # create forward transition matrix
  96. tf = tdlm.seq2tf(settings.seq_12)
  98. # run the simulation twice, once scaling replace n_events from 50-100% once from 0-100% depending on performance
  99. names_perf_scale = ['best-case', 'linear']
  101. #%% simulation parameters
  103. mode = 'erp_diff_all' # take the ERP difference
  104. lag_sim = 8 # simulate replay at 80 milliseconds
  105. sequence = tdlm.utils.char2num(settings.seq_12)[:-1]
  106. tp = 31 # best decoding index precomputed, ~210ms after stim onset
  107. best_C = 9.1 # cross-validated regularization strength
  109. # refractory period to block before and after each event that we insert
  110. # before, twice the lag must be blocked, as other events can else be put right
  111. # beforehand and their second reactivation would reach into the next one
  112. # After, only one lag must be blocked. See tdlm.utils.insert_events for details
  113. refractory = [lag_sim*2, lag_sim]
  116. # turn off saving resting state simulation data, as the resulting array
  117. # is >64GB large. Mainly used for debugging purposes
  118. save_simulation_data_in_memory = False
  119. n_jobs = -1 # set this to a lower number (eg 1) if you run out of memory
  121. if not final_calculation:
  122. print('not running final calculation yet, only use this for final paper calculation')
  124. # best_tp will be calculated and replaced later, for debugging it is sometimes
  125. # easier to define it here, so you can run segments without computing everything
  126. best_tp = utils.load_pkl(f"{settings.cache_dir}/best_tp.pkl.zip", 31)
  128. palette = sns.color_palette("ch:start=.2,rot=-.3", n_colors=101)
  129. hues = {subj:palette[int(100*(get_performance(subj)-0.5)*2)] for subj in subjects}
  131. # gaussian window created by gaussian_filter1d(np.float_([0,0,1,0,0]), 1)
  132. gaussian_weighting = [0.05842299, 0.24210528, 1, 0.24210528, 0.05842299]
  133. #%% preload some data (e.g. localizer)
  134. localizer = {} # localizer training data
  135. rs1 = {}
  136. rs2 = {}
  137. neg_x = {} # pre-audio fixation cross neg_x of localizer
  138. seqs = {} # load sequences of that participant in this dict
  140. for subj in tqdm(subjects, desc="Loading data"):
  141. # data used for the localizer
  142. localizer[subj] = load_localizers_seq12(subj=subj, sfreq=sfreq, bands=bands, autoreject=settings.default_autoreject, ica=settings.default_ica_components)
  143. # negative examples from the fixation cross before audio cue onset
  144. neg_x[subj] = load_neg_x_before_audio_onset(subj=subj, sfreq=sfreq, bands=bands, autoreject=settings.default_autoreject, ica=settings.default_ica_components)
  146. # resting state data, both eyes open and eyes closed together
  147. rs1[subj] = load_RS1(subj=subj, sfreq=sfreq, bands=bands, final_calculation=final_calculation)
  148. rs2[subj] = load_RS2(subj=subj, sfreq=sfreq, bands=bands, final_calculation=final_calculation)
  150. # individual sequences for the trials (maybe not necessary)
  151. seqs[subj] = utils.get_sequences(subj)
  154. max_accuracy = [utils.get_decoding_accuracy(subj, clf=clf, n_splits=20)[0] for subj in subjects]
  155. test_performance = [utils.get_performance(subj, which='test') for subj in subjects]
  156. df_subjects = pd.DataFrame({'subject': subjects,
  157. 'accuracy': max_accuracy,
  158. 'performance': test_performance}).sort_values('accuracy')
  160. # define which participants are rejected due to preregistered criteria
  161. excluded_acc = [subj for subj in subjects if utils.get_decoding_accuracy(subj,
  162. clf=clf, n_splits=20)[0]<0.3]
  163. excluded_perf = [subj for subj in subjects if utils.get_performance(subj, which='test')<0.5]
  164. excluded_miss = [subj for subj in subjects if utils.get_responses_localizer(subj)['n_misses']>60*0.25]
  166. df_excluded = pd.DataFrame({'reason': ['low decoding']*len(excluded_acc) + \
  167. ['low performance']*len(excluded_perf) + \
  168. ['>25% misses']*len(excluded_miss)},
  169. index=excluded_acc + excluded_perf + excluded_miss)
  170. df_excluded.sort_index(inplace=True)
  172. subjects_incl = sorted(set(subjects).difference(set(df_excluded.index)))
  174. #%% %%% STUDY 1
  175. # %% general description
  177. # 0.1 - memory performance
  178. df_stats = pd.DataFrame()
  179. df_stats_blocks = pd.DataFrame()
  181. # # precomputation, this might take a whi
  182. for subj in tqdm(subjects, desc="cross validating decoding performance"):
  183. utils.get_decoding_accuracy(subj=subj, clf=clf)
  185. for subj in subjects:
  186. # retrieve precomputed cross validation results
  187. acc = utils.get_decoding_accuracy(subj=subj, clf=clf)[0]
  189. # get participant performance for learning and testing
  190. val_learn = get_performance(subj=subj, which="learning")
  191. val_test = get_performance(subj=subj, which="test")
  192. scatter = (np.random.rand() * -0.5) * 0.02 + 1
  193. df_tmp = pd.DataFrame(
  194. {
  195. "participant": subj,
  196. "performance": [val_learn[-1] * scatter, val_test * scatter],
  197. "type": ["learning\n(last block)", "retrieval"],
  198. }
  199. )
  200. df_tmp_blocks = pd.DataFrame(
  201. {
  202. "participant": subj,
  203. "performance": [val_learn[0], val_learn[-1], val_test],
  204. "block": ["first", "last", "test"],
  205. "n_blocks": len(val_learn),
  206. "acc": acc,
  207. }
  208. )
  209. df_stats = pd.concat([df_stats, df_tmp])
  210. df_stats_blocks = pd.concat([df_stats_blocks, df_tmp_blocks])
  212. # store some results for later usage
  213. df_stats_blocks.to_pickle(settings.cache_dir + "/df_stats_blocks.pkl")
  215. # start plotting
  216. fig = plt.figure(figsize=[4, 6])
  217. ax = fig.subplots(1, 1)
  218. sns.despine()
  219. df_stats_blocks = df_stats_blocks[df_stats_blocks.acc >= min_acc]
  220. sns.boxplot(
  221. data=df_stats_blocks,
  222. x="block",
  223. y="performance",
  224. color=sns.color_palette()[0],
  225. meanprops=meanprops,
  226. showmeans=True,
  227. ax=ax,
  228. )
  230. ax.set_title("Memory performance")
  231. sns.despine()
  232. plt.tight_layout()
  233. fig.savefig(results_dir + "/memory_performance.png")
  234. fig.savefig(results_dir + "/memory_performance.svg")
  235. fig.savefig(results_dir + "/memory_performance.eps")
  238. perf_first, perf_last, perf_test = df_stats_blocks.groupby("block")
  239. p = ttest_rel(perf_last[1].performance, perf_test[1].performance)
  240. mean_diff = perf_last[1].performance.mean() - perf_test[1].performance.mean()
  242. print(f"Learning performance first block {perf_first[1].performance.mean():.2f}")
  243. print(f"Learning performance last block {perf_last[1].performance.mean():.2f}")
  244. print(f"Learning increase last->test {p=}")
  246. #%% Localizer get best C value
  247. import sklearn
  248. utils.lowpriority()
  249. ex_per_fold = 4
  250. df_reg = pd.DataFrame()
  252. clf_x = sklearn.base.clone(clf)
  254. # we originally did this over a much larger selection of C, but to reduce
  255. # computation time, we select this subrange now. Most C outside the range were
  256. # useless.
  257. Cs = np.logspace(-1.2, 2.5, 25, dtype=np.float16)
  259. for C in tqdm(Cs):
  260. clf_x.set_params(C=C)
  261. resx = Parallel(4)(delayed(utils.get_best_timepoint)(*localizer[subj], subj=subj,
  262. clf=clf_x, ex_per_fold=8,
  263. add_null_data=True,
  264. verbose=False)
  265. for subj in subjects)
  266. df_reg = pd.concat([df_reg]+ resx, ignore_index=True)
  268. df_reg['C'] = np.repeat(Cs.astype(float), len(df_reg)/len(Cs))
  270. joblib.dump(df_reg, settings.results_dir + '/localizer_c.pkl.gz')
  273. # plot results
  274. fig, ax = plt.subplots(figsize=[8, 6])
  275. df_reg = joblib.load(settings.results_dir + '/localizer_c.pkl.gz')
  276. df_reg_mean = df_reg[(df_reg.timepoint>150) & (df_reg.timepoint<250)]
  277. df_reg_mean = df_reg_mean.groupby(['C', 'timepoint']).mean(True).reset_index()
  278. df_reg_peak = df_reg_mean.groupby(['C']).mean(True).reset_index()
  279. best_C = df_reg_peak.C[df_reg_peak.accuracy.argmax()]
  281. ax.set(xscale="log")
  282. sns.lineplot(df_reg_mean, x='C', y='accuracy', ax = ax)
  283. ax.vlines(best_C, 0.1, 0.3, color='black')
  284. best_C
  286. ax.set_xlabel('L1 regularization (C), logscale')
  287. ax.set_ylabel('Mean accuracy between 150-250 ms after stim onset')
  288. ax.set_title(f'Best decoding at L1 ~C={best_C:.1f}')
  289. plt.pause(0.1)
  290. plt.tight_layout()
  291. utils.savefig(fig, 'supplement/regularization_C_crossval.png')
  292. utils.savefig(fig, 'supplement/regularization_C_crossval.svg')
  293. utils.savefig(fig, 'supplement/regularization_C_crossval.eps')
  294. #%% localizer decoding accuracy
  296. # (basically a copy of previous publication code)
  297. # Procedure:
  298. # For each participant:
  299. # 1. calculate best timepoint
  300. # 2. save best tp in list of all participants
  301. #
  302. # Choose best classifier TP for each participants based on LOSO
  303. # - use data of all other participant for this participants best_tp
  304. # - in our case, this yields the same TP for all participants
  307. # preload data that has already been computed
  308. pkl_localizer = settings.cache_dir + "/1 localizer.pkl.zip"
  309. results_localizer = load_pkl_pandas(
  310. pkl_localizer,
  311. default=pd.DataFrame(columns=["subject"])
  312. )
  314. ### Calculation of best timepoint and leave-one-out cross-validation
  315. ex_per_fold = 1
  316. for i, subj in enumerate(tqdm(subjects, desc="subject")):
  318. # this is quite time consuming, so caching is implemented.
  319. # you can continue computation at a later stage if you interrupt.
  320. # if subj in results_localizer["subject"].values:
  321. # # skip subjects that are already computed
  322. # print(f"{subj} already computed")
  323. # continue
  325. data_x, data_y = localizer[subj]
  326. res = get_best_timepoint(
  327. data_x, data_y, subj=f"{subj}", n_jobs=-1, ex_per_fold=8, clf=clf
  328. )
  329. results_localizer = pd.concat([results_localizer, res], ignore_index=True)
  331. results_localizer.to_pickle(pkl_localizer) # store intermediate results
  334. max_acc_subj = [
  335. results_localizer.groupby("subject")
  336. .get_group(subj)
  337. .groupby("timepoint")
  338. .mean(True)["accuracy"]
  339. .max()
  340. for subj in subjects
  341. ]
  342. # only include participants with decoding accuracy higher than 0.3
  343. included_idx = np.array(max_acc_subj) >= min_acc
  344. included_subj = np.array(subjects)[included_idx]
  346. subj_accuracy_all = (
  347. results_localizer.groupby(["subject", "timepoint"]).mean(True).reset_index()
  348. )
  350. # compute best decoding time point
  351. best_tp_tmp = [
  352. subj_accuracy_all.groupby("subject").get_group(subj)["accuracy"].argmax()
  353. for subj in subjects
  354. ]
  356. best_tp_tmp = np.array(best_tp_tmp)
  357. best_acc_tmp = [
  358. subj_accuracy_all.groupby("subject").get_group(subj)["accuracy"].max()
  359. for subj in subjects
  360. ]
  361. best_acc_tmp = np.array(best_acc_tmp)
  363. # MEAN
  364. best_tp = int(np.round(np.mean(best_tp_tmp[included_idx])))
  365. best_acc = np.mean(best_acc_tmp[included_idx])
  367. # store best timepoint for later retrieval
  368. joblib.dump(best_tp, f"{settings.cache_dir}/best_tp.pkl.zip")
  370. results_localizer_included = pd.concat(
  371. [results_localizer.groupby("subject").get_group(subj) for subj in included_subj]
  372. )
  374. subj_acc_included = (
  375. results_localizer_included.groupby(["subject", "timepoint"])
  376. .mean(True)
  377. .reset_index()
  378. )
  380. # start plotting of mean decoding accuracy
  381. fig = plt.figure(figsize=[6, 6])
  382. ax = fig.subplots(1, 1)
  383. sns.despine()
  384. utils.plot_decoding_accuracy(
  385. results_localizer_included, x="timepoint", y="accuracy", ax=ax, color="tab:blue"
  386. )
  387. ax.set_ylabel("Accuracy")
  388. ax.set_xlabel("ms after stimulus onset")
  389. ax.set_title("Localizer: decoding accuracy")
  390. ax.set_ylim(0, 0.5)
  392. ax.axvspan(
  393. (best_tp * 10 - 100) - 5, (best_tp * 10 - 100) + 5, alpha=0.3, color="orange"
  394. )
  395. ax.legend(["decoding accuracy", "95% conf.", "chance", "peak"], loc="lower right")
  397. plt.tight_layout()
  398. plt.pause(0.1)
  399. fig.savefig(results_dir + "/localizer_decoding.svg", bbox_inches="tight")
  400. fig.savefig(results_dir + "/localizer_decoding.png", bbox_inches="tight")
  401. fig.savefig(results_dir + "/localizer_decoding.eps", bbox_inches="tight")
  403. #%% localizer heatmap transfer across time
  405. # store precomputed results
  406. pkl_loc = cache_dir + "/1b heatmap loc.pkl.zip"
  408. # get decoder accuracy per participants to exclude low decodable participants
  409. accuracies = [utils.get_decoding_accuracy(subj=subj, clf=clf)[0] for subj in subjects]
  410. subj_incl = [subj for subj, acc in zip(subjects, accuracies) if acc > 0.3]
  412. # calculate time by time decoding heatmap from localizer
  413. # basically: How well can a clf trained on t1 predict t2 of the localizer
  414. if os.path.isfile(pkl_loc):
  415. # either load precomputed results or compute and store
  416. maps_localizer = joblib.load(pkl_loc)
  417. else:
  418. # use a parallel pool, as this is computationally quite expensive
  419. pool = Parallel(n_jobs) # use parallel pool
  420. maps_localizer = pool(
  421. delayed(utils.get_decoding_heatmap)(clf, *localizer[subj], n_jobs=2)
  422. for subj in subj_incl
  423. )
  424. joblib.dump(maps_localizer, pkl_loc)
  427. # zscore maps for statistics
  428. maps_loc_norm = np.array(maps_localizer) - np.mean(maps_localizer)
  429. maps_loc_norm = maps_loc_norm / maps_loc_norm.std()
  431. # perform cluster permutation testing, basically same as when doing fMRI
  432. t_thresh = stats.distributions.t.ppf(1 - 0.05, df=len(maps_localizer) - 1)
  433. t_clust, clusters1, p_values1, H0 = permutation_cluster_1samp_test(
  434. maps_loc_norm,
  435. tail=1,
  436. n_jobs=None,
  437. threshold=t_thresh,
  438. adjacency=None,
  439. n_permutations=1000,
  440. out_type="mask",
  441. )
  444. # create mask from cluster for all clusters of p<0.05
  445. clusters_sum1 = (np.array(clusters1)[p_values1 < 0.05]).sum(0)
  447. # now plot the heatmaps with masking using MNE visualization functions
  448. fig = plt.figure(figsize=[8, 6])
  449. ax = fig.subplots(1, 1)
  450. x = mne.viz.utils._plot_masked_image(
  451. ax,
  452. np.mean(maps_localizer, 0),
  453. times=range(61),
  454. mask=clusters_sum1,
  455. cmap="viridis",
  456. mask_style="contour",
  457. vmin=0.1,
  458. vmax=0.4,
  459. )
  460. plt.colorbar(x[0])
  463. ax.set_title('Localizer generalization')
  464. times = np.arange(-100, 500, 50) # t
  465. ax.set_xticks(np.arange(0, 60, 5), times, rotation=40)
  466. ax.set_yticks(np.arange(0, 60, 5), times)
  467. ax.set_xlabel("localizer train time (ms after onset)")
  468. ax.set_ylabel(
  469. f'{"retrieval " if i==1 else "localizer "} test time (ms after onset)'
  470. )
  472. fig.tight_layout()
  473. fig.savefig(results_dir + "/classifier-transfer.svg")
  474. fig.savefig(results_dir + "/classifier-transfer.png")
  475. fig.savefig(results_dir + "/classifier-transfer.eps")
  477. #%% RS1 vs RS2
  479. tp = best_tp
  481. clf.set_params(C=best_C)
  482. # put forward and backward sequenceness per subject for RS1/RS2 in these arrays
  484. rs1_sf = np.full([len(subjects_incl), n_shuf, max_lag+1], np.nan)
  485. rs1_sb = np.full([len(subjects_incl), n_shuf, max_lag+1], np.nan)
  486. rs2_sf = np.full([len(subjects_incl), n_shuf, max_lag+1], np.nan)
  487. rs2_sb = np.full([len(subjects_incl), n_shuf, max_lag+1], np.nan)
  489. # # create axes to plot subject sequenceness into
  490. fig1, axs1 = utils.make_fig(n_axs=len(subjects_incl), bottom_plots=0, suptitle=f'RS control')
  491. fig2, axs2 = utils.make_fig(n_axs=len(subjects_incl), bottom_plots=0, suptitle=f'RS post-learning')
  493. for i, subj in enumerate(tqdm(subjects_incl, desc="subject")):
  494. # train our classifier using localizer data and negative data
  495. train_x, train_y = localizer[subj]
  496. clf.reset_random_state(misc.make_seed(i, "rs"))
  497. clf.fit(train_x[:, :, tp], train_y, neg_x=neg_x[subj], neg_x_ratio=2.0)
  499. # get probability estimates for item reactivation from the resting state
  500. proba_rs1 = clf.predict_proba(rs1[subj])
  501. proba_rs2 = clf.predict_proba(rs2[subj])
  503. # normalize probabilities
  504. proba_rs1 = eval(proba_norm)(proba_rs1)
  505. proba_rs2 = eval(proba_norm)(proba_rs2)
  507. # calculate sequenceness
  508. rs1_sf_subj, rs1_sb_subj = tdlm.compute_1step(proba_rs1, tf=tf, n_shuf=n_shuf,
  509. max_lag=max_lag, alpha_freq=alpha_freq,
  510. seed=misc.make_seed(subj, 'rs1'))
  511. rs2_sf_subj, rs2_sb_subj = tdlm.compute_1step(proba_rs2, tf=tf, n_shuf=n_shuf,
  512. max_lag=max_lag, alpha_freq=alpha_freq,
  513. seed=misc.make_seed(subj, 'rs2'))
  515. # zscore sequenceness scores
  516. rs1_sf[i, :] = zscore_multiaxis(rs1_sf_subj, axes=zscore_axes)
  517. rs1_sb[i, :] = zscore_multiaxis(rs1_sb_subj, axes=zscore_axes)
  518. rs2_sf[i, :] = zscore_multiaxis(rs2_sf_subj, axes=zscore_axes)
  519. rs2_sb[i, :] = zscore_multiaxis(rs2_sb_subj, axes=zscore_axes)
  521. # # # subject level plot
  522. tdlm.plot_sequenceness(rs1_sf[i, :], rs1_sb[i, :] , ax=axs1[i], which=['fwd', 'bkw'], rescale=False, plotsignflip=False)
  523. tdlm.plot_sequenceness(rs2_sf[i, :], rs2_sb[i, :] , ax=axs2[i], which=['fwd', 'bkw'], rescale=False, plotsignflip=False)
  524. axs1[i].text(0, 0, subj, alpha=0.5)
  525. axs2[i].text(0, 0, subj, alpha=0.5)
  526. plt.pause(0.01)
  528. utils.normalize_lims(axs1[:len(subjects_incl)])
  529. utils.normalize_lims(axs2[:len(subjects_incl)])
  531. # create sequenceness dictionary
  532. df_sequenceness = pd.DataFrame()
  533. pkl_sequencenes = f'{settings.cache_dir}/rs1rs2-sequenceness.pkl.zip'
  534. for s, subj in enumerate(subjects_incl):
  535. for c, condition in enumerate(['rs1', 'rs2']):
  536. sf = [rs1_sf, rs2_sf][c]
  537. sb = [rs1_sb, rs2_sb][c]
  539. for d, direction in enumerate(['fwd', 'bkw']):
  540. sx = np.abs([sf, sb][d])
  541. sx_subj = sx[s, 0]
  542. peak_time = np.nanargmax(np.nanmean(sx[:, 0, :], 0))
  543. peak_sequenceness = sx_subj[peak_time]
  544. mean_sequenencess = np.nanmean(sx_subj)
  545. df_tmp = pd.DataFrame({'subject': subj,
  546. 'peak': peak_sequenceness,
  547. 'mean': mean_sequenencess,
  548. 'peak_time': peak_time,
  549. 'condition': condition,
  550. 'direction': direction},
  551. index=[subj])
  552. df_sequenceness = pd.concat([df_sequenceness, df_tmp],
  553. ignore_index=True)
  554. joblib.dump(df_sequenceness, pkl_sequencenes)
  556. #%% RS1 and RS2 main plot
  557. # plotting of sequenceness curves
  558. fig, axs = plt.subplot_mosaic([['1', '1', 'D1', 'D1'],
  559. ['2', '2', 'D2', 'D2'],
  560. ['S1', 'S2', 'S3', 'S4'],
  561. ['A', 'B', 'C', 'D']], figsize=[14, 14])
  563. perf_test = {subj: get_performance(subj=subj, which="test") for subj in subjects_incl}
  564. perf_diff = {subj: perf_test[subj] - get_performance(subj=subj, which="learning")[-1] for subj in subjects_incl}
  566. sequence_names = ['Control', 'Post-Learn']
  568. c_fwd = [sns.color_palette("bright")[1], sns.color_palette('dark')[1]]
  569. c_bkw = [sns.color_palette("bright")[2], sns.color_palette('dark')[2]]
  571. #### plot forward and backward
  572. tdlm.plot_sequenceness(rs1_sf, rs1_sb, which='fwd', title=f'Forward Sequenceness ',
  573. ax=axs['1'], rescale=True, color=c_fwd[0], clear=True, label='Control')
  574. tdlm.plot_sequenceness(rs2_sf, rs2_sb, which='fwd', title=f'Forward Sequenceness',
  575. ax=axs['1'], rescale=True, color=c_fwd[1], clear=False, label='Post-Learn')
  577. tdlm.plot_sequenceness(rs1_sf, rs1_sb, which='bkw', title=f'Backward Sequenceness',
  578. ax=axs['2'], rescale=True, color=c_bkw[0], clear=True, label='Control')
  579. tdlm.plot_sequenceness(rs2_sf, rs2_sb, which='bkw', title=f'Backward Sequenceness',
  580. ax=axs['2'], rescale=True, color=c_bkw[1], clear=False, label='Post-Learn')
  582. for ax in (axs_seq:=[axs['1'], axs['2']]):
  583. ax.set_ylabel('sequenceness\n(u.u., zscored)')
  584. ax.set_xlabel('time lag (milliseconds)')
  585. ax.legend(handles=[ln for ln in ax.lines if ln.get_label()])
  586. ax.plot([0, 0.000001], [0, .000000001], color='gray', linestyle='--', linewidth=1) # just for legend creation
  587. ax.plot([0, 0.000001], [0, .000000001], color='gray', linestyle='--', linewidth=1.5) # just for legend creation
  588. hndls = ax.get_legend_handles_labels()
  589. ax.legend([hndls[0][i] for i in [0, 3, 4, 5]],
  590. [hndls[1][i] for i in [0, 3, 4, 5]],
  591. loc='lower right',
  592. ncols=2, fontsize=12)
  595. utils.normalize_lims(axs_seq)
  597. for i, direction in enumerate(['forward', 'backward']):
  598. ax = [axs['1'], axs['2']][i]
  599. c = [c_fwd, c_bkw][i][0]
  600. rs_pre = [rs1_sf, rs1_sb][i]
  601. peak_lag = np.nanargmax(abs(np.nanmean(rs_pre[:, 0, :], 0)))
  602. peak_vals = dict(zip(subjects_incl, rs_pre[:, 0, peak_lag]))
  603. ax.vlines(peak_lag*10, *ax.get_ylim(), color=c, alpha=0.5)
  604. ax = [axs['A'], axs['C']][i]
  605. r1, pval1 = plot_correlation(peak_vals, values=perf_test,
  606. color=c, ax=ax,title=f"Control")
  608. ax.text(ax.get_xlim()[1], .5, f'r={r1:.2f} p={pval1:.3f}', horizontalalignment='right')
  610. ax = [axs['1'], axs['2']][i]
  611. c = [c_fwd, c_bkw][i][1]
  612. rs_post = [rs2_sf, rs2_sb][i]
  613. peak_lag = np.nanargmax(abs(np.nanmean(rs_post[:, 0, :], 0)))
  614. peak_vals = dict(zip(subjects_incl, rs_post[:, 0, peak_lag]))
  615. ax.vlines(peak_lag*10, *ax.get_ylim(), color=c, alpha=0.5)
  616. ax = [axs['B'], axs['D']][i]
  617. r1, pval1 = plot_correlation(peak_vals, values=perf_test,
  618. color=c, ax=ax,title=f"Post-Learn")
  620. ax.text(ax.get_xlim()[1], .5, f'r={r1:.2f} p={pval1:.3f}', horizontalalignment='right')
  622. utils.normalize_lims([ax for desc, ax in axs.items() if not desc.isnumeric() and len(desc)==1])
  625. #### signflip test
  626. colors = [sns.color_palette("bright")[1:3], sns.color_palette("dark")[1:3]]
  628. for s, session in enumerate(['Control', 'Post-Learn']):
  629. for d, direction in enumerate(['Forward', 'Backward']):
  630. sx = [[rs1_sf, rs1_sb], [rs2_sf, rs2_sb]][s][d]
  631. pval, t_base, t_perm = tdlm.signflit_test(sx[:, 0, :], rng=misc.make_seed(s, d))
  632. print(f'{session=} {direction=} {pval=:.3f} {t_base=:.3f}')
  633. ax = axs[f'S{d*2+s+1}']
  634. ax.clear()
  635. color = colors[s][d]
  636. tdlm.plot_tval_distribution(t_base, t_perm, color=color, ax=ax)
  637. ax.set_xlabel('sign-flip t-value distribution')
  638. ax.set_title(f'{session}')
  640. leg = ax.get_legend()
  641. plt.setp(leg.get_texts(), fontsize='smaller')
  642. plt.setp(leg.get_title(), fontsize='smaller')
  646. #### RS2-RS1 post cluster permutation
  647. from mne.stats import permutation_cluster_1samp_test
  648. fontdict ={ 'fontsize': 18, 'horizontalalignment':'center'}
  650. rs_pre_sf = rs1_sf[:, 0, 1:]
  651. rs_pre_sb = rs1_sb[:, 0, 1:]
  652. rs_post_sf = rs2_sf[:, 0, 1:]
  653. rs_post_sb = rs2_sb[:, 0, 1:]
  655. n_subj = len(rs_pre_sf)
  657. # create differences between post and pre
  658. diff_sf = rs_post_sf - rs_pre_sf
  659. diff_sb = rs_post_sb - rs_pre_sb
  661. df_diff = pd.DataFrame()
  662. for lag in range( max_lag):
  663. df_diff = pd.concat([df_diff,
  664. pd.DataFrame({'participant': list(range(n_subj))*2,
  665. 'direction': ['Forward'] * n_subj + ['Backward']*n_subj,
  666. 'timelag': 10+lag*10,
  667. 'difference (u.u.)': list(diff_sf[:, lag]) + list(diff_sb[:, lag])
  668. })],
  669. ignore_index=True)
  671. for i, direction in enumerate(['Forward', 'Backward']):
  672. ax = axs[f'D{i+1}']
  673. ax.clear()
  674. ax.hlines(0, 0, max_lag*10+10, color='gray', linestyle='--', alpha=0.2, label='_nolegend_')
  675. sns.lineplot(df_diff[df_diff.direction==direction], x='timelag', y='difference (u.u.)', ax=ax,
  676. color=sns.color_palette()[i+1], err_style="bars", errorbar=("se", 2),
  677. err_kws={'fmt':'o-', 'capsize':5})
  678. ax.set_ylim([x*1.5 for x in ax.get_ylim()])
  679. ax.set_title(f'{direction} Difference Post minus Control')
  680. ax.legend(['diff', 'SE'], loc='lower right')
  682. utils.normalize_lims( [axs[f'D1'], axs[f'D2']])
  684. p_fwd = stats.ttest_rel(rs_pre_sf, rs_post_sf, axis=0, nan_policy='omit')
  685. p_bkw = stats.ttest_rel(rs_pre_sb, rs_post_sb, axis=0, nan_policy='omit')
  687. _, p_fwd_corr, _, _ = multipletests(p_fwd.pvalue, method='fdr_bh')
  688. _, p_bkw_corr, _, _ = multipletests(p_bkw.pvalue, method='fdr_bh')
  691. for idx in np.where(p_fwd.pvalue<0.05)[0]:
  692. # axs[f'D1'].text(idx*10+10, axs[f'D1'].get_ylim()[1]*0.75,'(*)', fontdict=fontdict)
  693. print((idx+1)*10, f'{p_fwd.pvalue[idx]=:.3f} {p_fwd_corr[idx]=:.3f}')
  695. for idx in np.where(p_bkw.pvalue<0.05)[0]:
  696. # axs[f'D2'].text(idx*10+10, axs[f'D2'].get_ylim()[1]*0.75,'(*)', fontdict=fontdict)
  697. print((idx+1)*10, f'{p_bkw.pvalue[idx]=:.3f} {p_bkw_corr[idx]=:.3f}')
  699. plt.pause(0.1)
  700. fig.tight_layout()
  702. for i, diff in enumerate([diff_sf, diff_sb]):
  703. t_thresh = stats.distributions.t.ppf(1 - 0.05, df=len(diff) - 1)
  704. t_clust, clusters1, p_values1, H0 = permutation_cluster_1samp_test(
  705. diff,
  706. tail=1,
  707. n_jobs=None,
  708. threshold=t_thresh,
  709. adjacency=None,
  710. n_permutations=10000,
  711. out_type="mask",
  712. )
  713. print(clusters1, p_values1)
  715. utils.savefig(fig, f'figure/sequenceness_rs.png')
  716. utils.savefig(fig, f'figure/sequenceness_rs.svg')
  717. utils.savefig(fig, f'figure/sequenceness_rs.eps')
  719. #%% correlation with behaviour
  720. fig2, axs2 = plt.subplots(2, 2, figsize=[12,8 ])
  722. for i, direction in enumerate(['forward', 'backward']):
  723. ax = axs2[i][0]
  724. ax.clear()
  725. c = [c_fwd, c_bkw][i][0]
  726. rs_pre = [rs1_sf, rs1_sb][i]
  727. peak_lag = np.nanargmax(abs(np.nanmean(rs_pre[:, 0, :], 0)))
  728. peak_vals = dict(zip(subjects_incl, rs_pre[:, 0, peak_lag]))
  729. r1, pval1 = plot_correlation(peak_vals, values=perf_diff,
  730. color=c, ax=ax, title=f"Control")
  731. ax.text(ax.get_xlim()[1], -0.2, f'r={r1:.2f} p={pval1:.3f}', horizontalalignment='right')
  732. ax.set_ylabel('perf. diff. post-pre')
  734. ax = axs2[i][1]
  735. c = [c_fwd, c_bkw][i][1]
  736. rs_post = [rs2_sf, rs2_sb][i]
  737. peak_lag = np.nanargmax(abs(np.nanmean(rs_post[:, 0, :], 0)))
  738. peak_vals = dict(zip(subjects_incl, rs_post[:, 0, peak_lag]))
  739. r1, pval1 = plot_correlation(peak_vals, values=perf_diff,
  740. color=c, ax=ax,title=f"Post-Learn")
  741. ax.text(ax.get_xlim()[1], -0.2, f'r={r1:.2f} p={pval1:.3f}', horizontalalignment='right')
  742. ax.set_ylabel('perf. diff. post-pre')
  743. utils.normalize_lims(axs2.flatten())
  745. fig2.suptitle('Correlation between performance deltas and peak sequenceness')
  746. plt.pause(0.1)
  747. fig2.tight_layout()
  748. utils.savefig(fig2, f'supplement/correlation_pre-post-difference.png')
  749. utils.savefig(fig2, f'supplement/correlation_pre-post-difference.svg')
  750. utils.savefig(fig2, f'supplement/correlation_pre-post-difference.eps')
  751. #%% RS1RS2 all sequenceness curves in one plot
  753. df_rs1rs2 = pd.DataFrame()
  754. fig, axs = plt.subplots(2, 2, figsize=[12, 12])
  755. axs = axs.flatten()
  756. for i, sx in enumerate([rs1_sf, rs2_sf, rs1_sb, rs2_sb]):
  757. condition = ['forward control', 'backward post-learn', 'backward control', 'backward post-learn'][i]
  758. sx = zscore_multiaxis(sx[:, :, :], axes=zscore_axes)
  759. seq = sx[:, 0, :].ravel()
  760. time_lags = list(range(0, sx.shape[-1]*10, 10)) * len(subjects_incl)
  761. df_tmp = pd.DataFrame({'subject': np.repeat(subjects_incl, sx.shape[-1]),
  762. 'sequenceness': seq,
  763. 'time lag': time_lags})
  764. ax = axs[i]
  765. sns.lineplot(df_tmp, x='time lag', y='sequenceness', hue='subject', ax=ax,
  766. palette=hues, legend=False)
  767. ax.set_title(f'{condition}')
  769. fig.suptitle('Individual Sequenceness Curves of Participants')
  770. utils.normalize_lims(axs)
  771. utils.savefig(fig, 'supplement/sequenceness-individualplots.png')
  772. utils.savefig(fig, 'supplement/sequenceness-individualplots.svg')
  773. utils.savefig(fig, 'supplement/sequenceness-individualplots.eps')
  775. #%% RS1 vs RS2 segments
  777. tp = 31#np.mean(list(best_tp.values())).astype(int)
  779. # # create forward transition matrix
  780. tf = tdlm.seq2tf(settings.seq_12)
  782. # put forward and backward sequenceness per segment for RS1/RS2 in these arrays
  783. rs1_sf_seg = np.zeros([len(subjects_incl), 8, n_shuf, max_lag+1])
  784. rs1_sb_seg = np.zeros([len(subjects_incl), 8, n_shuf, max_lag+1])
  785. rs2_sf_seg = np.zeros([len(subjects_incl), 8, n_shuf, max_lag+1])
  786. rs2_sb_seg = np.zeros([len(subjects_incl), 8, n_shuf, max_lag+1])
  788. # # create axes to plot subject sequenceness into
  789. fig1, axs1 = utils.make_fig(n_axs=len(subjects_incl), bottom_plots=0, suptitle='RS control')
  790. fig2, axs2 = utils.make_fig(n_axs=len(subjects_incl), bottom_plots=0, suptitle='RS post-learning')
  792. n_segments = 8
  794. for i, subj in enumerate(tqdm(subjects_incl, desc="subject")):
  795. # train our classifier using localizer data and negative data
  796. train_x, train_y = localizer[subj]
  797. clf.reset_random_state(misc.make_seed(i, "rs"))
  798. clf.fit(train_x[:, :, tp], train_y, neg_x=neg_x[subj], neg_x_ratio=2.0)
  800. # there are 8 segments that are approximately 30 seconds long
  801. segments_rs1 = np.split(rs1[subj], n_segments)
  802. segments_rs2 = np.split(rs2[subj], n_segments)
  804. for seg in range(n_segments):
  806. # get probability estimates for item reactivation from the resting state
  807. proba_rs1seg = clf.predict_proba(segments_rs1[seg])
  808. proba_rs2seg = clf.predict_proba(segments_rs2[seg])
  809. proba_rs1seg = eval(proba_norm)(proba_rs1seg)
  810. proba_rs2seg = eval(proba_norm)(proba_rs2seg)
  812. # calculate sequenceness
  813. rs1_sf_subj, rs1_sb_subj = tdlm.compute_1step(proba_rs1seg, tf=tf,
  814. n_shuf=n_shuf, max_lag=max_lag,
  815. alpha_freq=alpha_freq,
  816. seed=misc.make_seed(subj, seg))
  817. rs2_sf_subj, rs2_sb_subj = tdlm.compute_1step(proba_rs2seg, tf=tf,
  818. n_shuf=n_shuf, max_lag=max_lag,
  819. alpha_freq=alpha_freq,
  820. seed=misc.make_seed(subj, seg))
  822. rs1_sf_seg[i, seg, :] = zscore_multiaxis(rs1_sf_subj, axes=zscore_axes)
  823. rs1_sb_seg[i, seg, :] = zscore_multiaxis(rs1_sb_subj, axes=zscore_axes)
  824. rs2_sf_seg[i, seg, :] = zscore_multiaxis(rs2_sf_subj, axes=zscore_axes)
  825. rs2_sb_seg[i, seg, :] = zscore_multiaxis(rs2_sb_subj, axes=zscore_axes)
  828. # plot as heatmap across segments
  829. sequenceness = [(rs1_sf_seg, 'forward', 'control'),
  830. (rs2_sf_seg, 'forward', 'post-learn'),
  831. (rs1_sb_seg, 'backward', 'control'),
  832. (rs2_sb_seg, 'backward', 'post-learn')]
  834. fig, axs = plt.subplots(2, 2, figsize=[14, 8])
  835. axs = axs.flatten()
  837. vmins = [np.nanmin(np.nanmean(sx[:,:,1:,:],0)) for sx, *_ in sequenceness]
  838. vmaxs = [np.nanmax(np.nanmean(sx[:,:,1:,:],0)) for sx, *_ in sequenceness]
  841. for i, (sx, direction, condition) in enumerate(sequenceness):
  842. heatmap = np.nanmean(sx[:, :, 0, :], 0)
  843. ax = axs[i]
  844. im = ax.imshow(heatmap, aspect='auto', cmap='PiYG', vmin=min(vmins), vmax=max(vmaxs))
  845. ax.set_xlabel('time lag')
  846. ax.set_xticks(np.arange(0, 31, 5), np.arange(0, 310, 50))
  847. ax.set_ylabel('min. in resting state')
  848. ax.set_yticks(np.arange(8), np.arange(1, 9))
  849. ax.set_title(f'{condition} - {direction}')
  851. cbar_ax = fig.add_axes([0.92, 0.1, 0.02, 0.8])
  852. fig.colorbar(im, cax=cbar_ax)
  854. plt.pause(0.1)
  855. fig.tight_layout(rect = [0, 0, 0.88, 1])
  856. cbar_ax.set_ylabel('sequenencess\n(u.u. zscored)', labelpad=10)
  857. cbar_ax.yaxis.set_label_position("left")
  858. utils.savefig(fig, 'figure/sequencess_blocks_heatmap.png', tight=False)
  859. utils.savefig(fig, 'figure/sequencess_blocks_heatmap.svg', tight=False)
  860. utils.savefig(fig, 'figure/sequencess_blocks_heatmap.eps', tight=False)
  862. # plot separately, as otherwise it get's too crowded in the plot
  863. fig1, axs = plt.subplots(2, 4, figsize=[16, 8], sharex=True, sharey=True)
  864. axs = axs.flatten()
  865. for seg in range(8):
  866. # first plot RS sequenceness individually
  867. tdlm.plot_sequenceness(rs1_sf_seg[:,seg, ...], rs1_sb_seg[:,seg, ...], which=['fwd', 'bkw'],
  868. title=f'Minute {seg+1}', ax=axs[seg])
  869. fig1.suptitle('Sequenceness Control Resting State ')
  870. axs[seg].legend(loc='upper left', ncols=2)
  871. plt.pause(0.1)
  872. utils.savefig(fig1, 'supplement/segments-rs1.png')
  873. utils.savefig(fig1, 'supplement/segments-rs1.svg')
  874. utils.savefig(fig1, 'supplement/segments-rs1.eps')
  876. fig2, axs = plt.subplots(2, 4, figsize=[16, 8], sharex=True, sharey=True)
  877. axs = axs.flatten()
  878. for seg in range(8):
  879. # first plot RS sequenceness individually
  880. tdlm.plot_sequenceness(rs2_sf_seg[:,seg, ...], rs2_sb_seg[:,seg, ...], which=['fwd', 'bkw'],
  881. title=f'Minute {seg+1}', ax=axs[seg])
  882. fig2.suptitle('Sequenceness Post-Learning Resting State ')
  883. axs[seg].legend(loc='upper left', ncols=2)
  884. plt.pause(0.1)
  886. utils.savefig(fig2, 'supplement/segments-rs2.png')
  887. utils.savefig(fig2, 'supplement/segments-rs2.svg')
  888. utils.savefig(fig2, 'supplement/segments-rs2.eps')
  894. #%% #### SUPPLEMENT
  895. # here are calculations that I need for the supplement
  896. #%% SUPPL: sensor patterns
  898. patterns = {img:[] for img in ['berg',
  899. 'schreibtisch',
  900. 'pinsel',
  901. 'kuchen',
  902. 'apfel',
  903. 'zebra',
  904. 'clown',
  905. 'fahrrad',
  906. 'tasse',
  907. 'fuß']}
  909. for i, subj in enumerate(tqdm(subjects, desc="subject")):
  910. # train our classifier using localizer data and negative data
  911. train_x, train_y = localizer[subj]
  912. clf.fit(train_x[:, :, tp], train_y, neg_x=neg_x[subj], neg_x_ratio=2.0)
  914. # get data from localizer that is inserted into the resting state
  915. insert_data_tp = train_x[:, :, tp] # take the ERP component at peak decodability
  917. # ERP of all visual evoked activity
  918. erp = insert_data_tp.mean(0)
  920. # now for each class, subtract the mean EPR from the class ERP.
  921. # this way, we should isolate the pattern that is responsible
  922. # for the differences
  923. names_subj = utils.get_image_names(subj)
  925. for i, name in enumerate(names_subj):
  926. insert_class_pattern = [insert_data_tp[train_y==i, :].mean(0)-erp]
  927. patterns[name] += insert_class_pattern
  930. fig = plt.figure(figsize=[10, 8])
  931. axs = fig.subplots(3, 4)
  932. axs = axs.flatten()
  933. cmap = sns.color_palette("vlag", as_cmap=True)
  935. vmin = np.min([np.mean(values, 0) for values in patterns.values()])/1.2
  936. vmax = np.max([np.mean(values, 0).max() for values in patterns.values()])/1.2
  938. for i, name in enumerate(patterns):
  939. plotting.plot_sensors(np.mean(patterns[name], 0),
  940. ax=axs[i],
  941. vmin=vmin,
  942. vmax=vmax,
  943. mode="color",
  944. title=name, cmap=cmap)
  946. axs[-1].axis("off")
  947. axs[-2].axis("off")
  949. utils.savefig(fig, f'supplement/insertion_patterns.png')
  950. utils.savefig(fig, f'supplement/insertion_patterns.svg')
  951. utils.savefig(fig, f'supplement/insertion_patterns.eps')
  953. #%% SUPPL: visualize ERP
  954. from meg_utils import plotting
  955. from scipy.stats import zscore
  957. erps = np.zeros([len(localizer), 10, 306])
  958. for s, subj in enumerate(localizer):
  959. data_x, data_y = localizer[subj]
  960. names = utils.get_image_names(subj)
  961. names_base = sorted(names)
  962. for y in range(10):
  963. data_img = data_x[data_y==y]
  964. erps[s, names_base.index(names[y])] = data_img[:, :, 31].mean(0)
  966. erps_mean = erps.mean(0)
  968. fig = plt.figure(figsize=[10, 7])
  969. axs = fig.subplots(3, 4)
  970. axs = axs.flatten()
  971. plotting.plot_sensors(erps_mean.mean(0), ax=axs[0], mode="size",
  972. title='mean', cmap='RdYlBu')
  973. for i in range(10):
  974. plotting.plot_sensors(erps_mean[i], ax=axs[i+1], mode="size",
  975. title=names_base[i], cmap='RdYlBu')
  976. axs[-1].axis("off")
  977. axs[-2].axis("off")
  979. fig.suptitle("Sensor value of images")
  980. plt.pause(0.1)
  981. fig.tight_layout
  982. fig.savefig(results_dir + "/image-ERP.svg")
  983. fig.savefig(results_dir + "/image-ERP.png")
  984. fig.savefig(results_dir + "/image-ERP.eps")
  986. #%% SUPPL: visualize sensors betas
  987. from meg_utils import plotting
  989. def fit(clf, subj, *args, **kwargs):
  990. np.random.seed(misc.make_seed(subj))
  991. return clf.fit(*args, **kwargs)
  994. clfs = Parallel(len(localizer) - 1)(
  995. delayed(fit)(clf, subj, X=localizer[subj][0][:, :, 31], y=localizer[subj][1])
  996. for subj in localizer
  997. )
  999. # these are the orders of images as they were assigned to the classes
  1000. names = [utils.get_image_names(subj) for subj in subjects]
  1001. names_base = sorted(names[0])
  1003. # we need to sort the rows of the matrix accordingly
  1004. sorting = [[names_base.index(x) for x in name] for name in names]
  1005. betas = np.array([clf.coef_[idxs, :] for clf, idxs in zip(clfs, sorting)])
  1006. sensors_active = np.mean(betas>0, 0)
  1008. fig = plt.figure(figsize=[10, 7])
  1009. axs = fig.subplots(3, 4)
  1010. axs = axs.flatten()
  1011. plotting.plot_sensors(sensors_active.mean(0), ax=axs[0], mode="size",
  1012. title='mean')
  1013. for i in range(10):
  1014. plotting.plot_sensors(sensors_active[i], ax=axs[i+1], mode="size",
  1015. title=names_base[i])
  1016. axs[-1].axis("off")
  1017. axs[-2].axis("off")
  1019. fig.suptitle("Sensor distribution for different images")
  1020. plt.pause(0.1)
  1021. fig.tight_layout
  1022. fig.savefig(results_dir + "/S6 sensorlocation.svg")
  1023. fig.savefig(results_dir + "/S6 sensorlocation.png")
  1024. fig.savefig(results_dir + "/S6 sensorlocation.eps")

2_run_study1.py at commit eea0f83, under GPL-3.0 · at the source

Overview

Authors: Simon Kern1,2,3,4, Juliane Nagel1,2,3,4, Lennart Wittkuhn5, Steffen Gais6, Raymond J Dolan7,8, Gordon B Feld1,2,3,4
ORCID iDs: Simon Kern, Juliane Nagel, Lennart Wittkuhn, Raymond J Dolan, Gordon B Feld
  1. Clinical Psychology, Central Institute of Mental Health, Medical Faculty Mannheim, University of Heidelberg Mannheim Germany
  2. Psychiatry and Psychotherapy, Central Institute of Mental Health, Medical Faculty Mannheim, University of Heidelberg Mannheim Germany
  3. Addiction Behavior and Addiction Medicine, Central Institute of Mental Health, Medical Faculty Mannheim, University of Heidelberg Mannheim Germany
  4. Department of Psychology, Ruprecht Karl University of Heidelberg Heidelberg Germany
  5. Institute of Psychology, Universität Hamburg Hamburg Germany
  6. Institute of Medical Psychology and Behavioral Neurobiology, Eberhard-Karls-University Tübingen Tübingen Germany
  7. Max Planck UCL Centre for Computational Psychiatry and Ageing Research London United Kingdom
  8. Wellcome Centre for Human Neuroimaging, University College London London United Kingdom
Institutions: Central Institute of Mental Health (Germany); Heidelberg University (Germany); Universität Hamburg (Germany); University of Tübingen (Germany); Max Planck Society (Germany); University College London (United Kingdom)
Journal: eLife, volume 14, article RP108023
Dates: published online 9 April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.7554/elife.108023 · PMID 41954994 · PMCID PMC13065329 · OpenAlex W4414124402
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: MEG (modality), human (organism)
Methods: Spectral & time-frequency, Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Evoked potentials, fMRI & imaging
Keywords: replay, memory reactivation, decoding, MEG, TDLM, Human
MeSH: Brain*, Magnetoencephalography*, Memory*, Rest*, Computer Simulation, Humans, Linear Models (* major topic)
Journal subjects: Neuroscience
Topic: Memory and Neural Mechanisms (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Deutsche Forschungsgemeinschaft (10.3030/101170886, FE1617/2-1); Deutsche Gesellschaft für Schlafmedizin; German National Academic Foundation
Citations: cited by 1 paper (Europe PMC); 120 references in the paper
Research resources: RRID:SCR_002577, RRID:SCR_002630, RRID:SCR_004129, RRID:SCR_005972, using a Python RRID:SCR_008394, RRID:SCR_022515

Abstract

Using temporally delayed linear modeling (TDLM) and magnetoencephalography (MEG), we investigated whether items associated with an underlying graph structure are replayed during a post-learning resting state. In these same data, we previously provided evidence for replay during online (non-rest) memory retrieval. Despite successful decoding of brain activity during a localizer task, and contrary to predictions, we found no evidence for replay during a post-learning resting state. To better understand this, we performed a hybrid simulation analysis in which we inserted synthetic replay events into a control resting state recorded prior to the actual experiment. This simulation revealed that replay detection using our current pipeline requires an extremely high replay density to reach significance (>1 replay sequence per second, with ‘replay’ defined as a sequence of reactivations within a certain time lag). Furthermore, when scaling the number of replay events with a behavioral measure, we were unable to induce a strong correlation between sequenceness and this measure. We infer that even if replay was present at plausible rates in our resting state dataset, we would lack statistical power to detect it with TDLM. Finally, contrasting our novel hybrid simulation to existing purely synthetic simulations indicated that the latter approaches overestimate the sensitivity of TDLM. We discuss approaches that might optimize the analytic methodology, including identifying boundary conditions under which TDLM can be expected to detect replay. We conclude that solving these methodological constraints will be crucial for optimizing the non-invasive measurement of human replay using MEG.

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

Repositories

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

identifiers.org/rrid:scr_002630

License: none: the authors keep all their rights
State: unreachable at the last attempt, verified on 29 September 2026
Evidence: found in the paper
Software Heritage: not checked
Found in: the text, “MEG acquisition and pre-processing”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 2 checks, the latest on 29 September 2026: unreachable at the last attempt (HTTP 403)
  • 29 September 2026: unreachable at the last attempt (HTTP 403)
  • 29 September 2026: unreachable at the last attempt (HTTP 403)
At the source: identifiers.org/RRID:SCR_002630

cimh-clinical-psychology/desmrrest-tdlm-simulation

License: GPL-3.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: eea0f83745d8eca2fd046b83563d0d7e7377fd10, 10 April 2026
Languages: Python (12), MATLAB (7), JavaScript (6)
Size: 393 files, 25 scripts
Software Heritage: not archived
Found in: the references
Holds: README, license file, environment (requirements.txt), documentation
Not found: CITATION.cff, tests, continuous integration
Tools: NumPy (11 files), Matplotlib (8 files), pandas (8 files), seaborn (8 files), SciPy (6 files), scikit-learn (4 files), Statistics and Machine Learning Toolbox (3 files), MNE-Python (3 files), shadedErrorBar (2 files), autoreject (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
27 files
At the source: archive.softwareheritage.org/swh:1:dir:9fccd6b56a8af789ae3c7f4dfa5c2a52c0470b04;origin=https://github.com/CIMH-Clinical-Psychology/DeSMRRest-TDLM-Simulation;visit=swh:1:snp:3c3a5e25b762e5231534f917bca6832c16696c53;anchor=swh:1:rev:67119521f4348b93e7f362b38028aca66c010169

Zenodo 12623445

License: CC-BY-4.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (7 files), Statistics and Machine Learning Toolbox (6 files), SciPy (3 files), shadedErrorBar (3 files), Numba (2 files), scikit-learn (2 files), Matplotlib (1 file), pandas (1 file), seaborn (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
  • 29 September 2026: the link answers (HTTP 200)
30 files
At the source:

The paper's code and data availability statement is in the Data section.

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:

  • 3 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 53 scripts, each with its path and the digest of its content;
  • 19 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

Data availability MaxFiltered and anonymized MEG raw data are available at Zenodo in two parts: Localizer data is available as v1 of https://doi.org/10.5281/zenodo.8001755. Resting state data is available as v2 of https://doi.org/10.5281/zenodo.15629081. Code availability The code of the analysis as well as the experiment paradigm and the stimulus material are available at https://github.com/CIMH-Clinical-Psychology/DeSMRRest-TDLM-Simulation, copy archived at Kern, 2026.

The following previously published datasets were used:

KernS 2025Challenges in Replay Detection by TDLM in Post-Encoding Resting StateZenodo10.5281/zenodo.15629081PMC1306532941954994

KernS 2023Reactivation strength during cued recall is modulated by graph distance within cognitive mapsZenodo10.5281/zenodo.8001755PMC1113649338810249

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 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 6 authors, 6 keywords, 7 MeSH terms, 3 funders, 117 references, 6 RRIDs.

Cite

This paper

Kern, S., Nagel, J., Wittkuhn, L., Gais, S., Dolan, R. J., & Feld, G. B. (2026). Challenges in replay detection by TDLM in post-encoding resting state. eLife, 14, RP108023. https://doi.org/10.7554/elife.108023

BibTeX

@article{kern2026challenges,
author = {Kern, Simon and Nagel, Juliane and Wittkuhn, Lennart and Gais, Steffen and Dolan, Raymond J and Feld, Gordon B},
title = {{Challenges in replay detection by TDLM in post-encoding resting state}},
journal = {eLife},
year = {2026},
month = apr,
volume = {14},
pages = {RP108023},
publisher = {eLife Sciences Publications, Ltd},
issn = {2050-084X},
doi = {10.7554/elife.108023},
url = {https://doi.org/10.7554/elife.108023},
pmid = {41954994},
pmcid = {PMC13065329}
}

RIS

TY - JOUR
AU - Kern, Simon
AU - Nagel, Juliane
AU - Wittkuhn, Lennart
AU - Gais, Steffen
AU - Dolan, Raymond J
AU - Feld, Gordon B
TI - Challenges in replay detection by TDLM in post-encoding resting state
T2 - eLife
J2 - eLife
PY - 2026
DA - 2026/04/09
VL - 14
SP - RP108023
SN - 2050-084X
PB - eLife Sciences Publications, Ltd
DO - 10.7554/elife.108023
UR - https://doi.org/10.7554/elife.108023
LA - en
ER -

CSL-JSON

{
"id": "10.7554/elife.108023",
"type": "article-journal",
"title": "Challenges in replay detection by TDLM in post-encoding resting state",
"container-title": "eLife",
"author": [
{
"family": "Kern",
"given": "Simon"
},
{
"family": "Nagel",
"given": "Juliane"
},
{
"family": "Wittkuhn",
"given": "Lennart"
},
{
"family": "Gais",
"given": "Steffen"
},
{
"family": "Dolan",
"given": "Raymond J"
},
{
"family": "Feld",
"given": "Gordon B"
}
],
"container-title-short": "eLife",
"volume": "14",
"page": "RP108023",
"DOI": "10.7554/elife.108023",
"PMID": "41954994",
"PMCID": "PMC13065329",
"ISSN": "2050-084X",
"publisher": "eLife Sciences Publications, Ltd",
"URL": "https://doi.org/10.7554/elife.108023",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
9
]
]
}
}

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.1038/s41467-025-63281-w [code]
Trait anxiety is associated with reduced reward-related replay at rest
Journal: n/a
In common: shadedErrorBar, Statistics and Machine Learning Toolbox, 21 references, author Raymond J Dolan
[2] doi:10.1038/s41593-026-02357-2 [code]
Experience reorganizes content-specific memory traces in macaques.
Journal: Nature neuroscience
In common: Statistics and Machine Learning Toolbox, seaborn, scikit-learn, 4 other tools, 14 references
[3] doi:10.1038/s42003-025-08618-3 [code]
Physical activity simultaneously improves working memory and ripple-spindle coupling
Journal: n/a
In common: shadedErrorBar, Statistics and Machine Learning Toolbox, MEG, 8 references
[4] doi:10.1093/nc/niag029 [code]
A data-driven approach to identifying and evaluating connectivity-based neural correlates of conscious visual perception.
Journal: Neuroscience of consciousness
In common: autoreject, MNE-Python, statsmodels, 7 other tools, MEG, 3 references
[5] doi:10.1038/s41467-026-74822-2 [code]
A likelihood-based method for identifying replay from spike sequences.
Journal: Nature communications
In common: Statistics and Machine Learning Toolbox, 10 references
[6] doi:10.7554/elife.103956 [code]
Human brain-wide activation of sleep rhythms.
Journal: eLife
In common: pandas, NumPy, 10 references
[7] doi:10.1016/j.cub.2026.05.068 [code]
An abstract relational map emerges in the human medial prefrontal cortex with consolidation.
Journal: Current biology : CB
In common: statsmodels, Statistics and Machine Learning Toolbox, seaborn, 4 other tools, 6 references
[8] doi:10.1371/journal.pbio.3003938 [code]
Theta oscillations tag episodic memories for sleep-dependent consolidation.
Journal: PLoS biology
In common: shadedErrorBar, Statistics and Machine Learning Toolbox, 7 references
[9] doi:10.1038/s41467-026-75345-6 [code]
Hippocampal ripples initiate cortical dimensionality expansion for memory retrieval.
Journal: Nature communications
In common: Statistics and Machine Learning Toolbox, 9 references
[10] doi:10.1162/imag.a.1321 [code]
Phase similarity between similar objects indicates representational merging across retrieval training but not sleep.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: autoreject, MNE-Python, seaborn, 5 other tools, 4 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.