OSCR

Structural Eigenmodes of the Brain to Improve the Source Localization of EEG: Application to Epileptiform Activity.

Code ↔ Paper

14 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 14 matches
  1. [1] § Results › Structural Eigenmodes Enable Source Localization of Seizure Spread › Testing, Validation, and Comparison to Other Approaches › Comparison Between Structural Eigenmode Approaches ↔ STC_Analysis.ipynb, lines 501 to 531 · score 0.70 · region localization error, squared error, variance explained, cosine score, correlated, SNR · Report a mismatch
  2. [2] § Methods › Simulated Evaluation Data › Mapping to Sensor Space with MNE Python ↔ gem_fns.py, lines 561 to 619 · score 0.66 · sensor space, inverse problem, efficiency, oriented, source space, cortex · Report a mismatch
  3. [3] § Results › Structural Eigenmodes Enable Source Localization of Seizure Spread › Performance Against Number of Modes — The Over‐ Versus Under‐Parameterized Regimes ↔ STC_Analysis.ipynb, lines 143 to 245 · score 0.63 · region localization error, EEG noise, EEG power, cosine score, SNR, metrics · Report a mismatch
  4. [4] § Methods › Source Localization ↔ gem_fns.py, lines 473 to 559 · score 0.62 · Small singular, EEG modes, component, SVD, rcond, truncate · Report a mismatch
  5. [5] § Methods › Evaluation › Simulation Test Conditions ↔ sample_sesl_pipeline.ipynb, lines 134 to 150 · score 0.62 · dSPM, eLORETA, sLORETA, MNE · Report a mismatch
  6. [6] § Results › Structural Eigenmodes Enable Source Localization of Seizure Spread › Testing, Validation, and Comparison to Other Approaches › Performance of Established Methods ↔ sample_sesl_pipeline.ipynb, lines 134 to 150 · score 0.60 · dSPM, eLORETA, sLORETA, loose, MNE, noise · Report a mismatch
  7. [7] § Results › Structural Eigenmodes Enable Source Localization of Seizure Spread › Testing, Validation, and Comparison to Other Approaches › Performance of Established Methods ↔ STC_Analysis.ipynb, lines 501 to 531 · score 0.58 · sLORETA, variance explained, cosine score, loose, MNE, metrics · Report a mismatch
  8. [8] § Methods › Structural Brain Data ↔ dynsim_fns.py, lines 101 to 196 · score 0.56 · streamline length, connectivity matrix, Spherical, parcellation, threshold · Report a mismatch
  9. [9] § Methods › Simulated Evaluation Data › Epileptor Equations ↔ dynsim_fns.py, lines 436 to 507 · score 0.55 · state variables, Epileptor model, spike, TVB, seizures · Report a mismatch
  10. [10] § Methods › Evaluation › Simulation Test Conditions ↔ STC_Analysis.ipynb, lines 28 to 88 · score 0.55 · inverse solutions, forward model, HCP, fsaverage, signal, EEG · Report a mismatch
  11. [11] § Methods › Evaluation › Simulation Metrics › Region Localization Error ↔ STC_Analysis.ipynb, lines 143 to 245 · score 0.54 · Region localization error, EEG power, Otsu, MNE, Threshold · Report a mismatch
  12. [12] § Methods › Eigenmode Creation › Connectome Eigenmodes ↔ STC_Analysis.ipynb, lines 306 to 384 · score 0.54 · subcortical regions, Connectome eigenmodes, geometric eigenmodes, discrete, eigenvalue · Report a mismatch
  13. [13] § Methods › Structural Brain Data ↔ dynsim_fns.py, lines 267 to 325 · score 0.53 · subcortical structures, cortical region, atlas, parcellation, vertex · Report a mismatch
  14. [14] § Methods › Source Localization ↔ gem_fns.py, lines 473 to 559 · score 0.51 · EEG modes, singular, vector, SVD, operator, squares · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 720 lines · 33 KB · BSD-3-Clause · 6 matches

  1. # %% [markdown]
  2. # # Paper STC Analysis Notebook
  3. #
  4. # This notebook provides an example of how the functions in `sesl` (structural eigenmode source localisation) can be used to perform source localisation on a simulated dataset.
  5. # It runs on one example simulation, uses fsaverage, and demonstrates how the figures in the original paper can be generated (if run on all simulations and with all the different input files).
  6. # Remember to put the 1000 region parcellation labels into the appropriate folder in your Freesurfer fsaverage directory.
  7. # The aim for this notebook is for reproducibility of the analysis in the paper, so if you are interested in running `sesl` on your own data, please refer to the other notebook.
  9. # %%
  10. import numpy as np
  11. import matplotlib.pyplot as plt
  12. from pathlib import Path
  14. from IPython.display import clear_output
  16. import mne
  17. from mne.minimum_norm import make_inverse_operator, apply_inverse
  18. from skimage.filters import threshold_otsu
  20. from useful_fns import *
  21. from gem_fns import *
  22. from dynsim_fns import *
  23. from plot_fns import plot_mean_and_range, plot_tria_mesh
  25. %load_ext autoreload
  26. %autoreload 2
  28. # %%
  29. ANALYSIS_TYPE = 'GEM_discrete' # 'GEM_discrete' or 'GEM_surface' or 'CEM_discrete' or 'SEM_discrete'
  31. SNR = None # Signal to noise ratio for the simulation and MN inverse solutions - None for no noise
  32. if SNR is not None:
  33. RCOND = 1/(10**(SNR/10))
  34. else:
  35. RCOND = None
  37. ELECTRODE_CONFIG = '1010'
  38. if ELECTRODE_CONFIG == '1020':
  39. ONLY_20ELECTRODES = True # Whether to only use the 20 electrodes in the 10-20 system
  40. MONTAGE = 'standard_1020' # EEG montage - 'standard_1020' or 'standard_1005'
  41. elif ELECTRODE_CONFIG == '1010':
  42. ONLY_20ELECTRODES = False
  43. MONTAGE = 'standard_1020' # EEG montage - 'standard_1020' or 'standard_1005'
  44. elif ELECTRODE_CONFIG == '1005':
  45. ONLY_20ELECTRODES = False
  46. MONTAGE = 'standard_1005' # EEG montage - 'standard_1020' or 'standard_1005'
  48. LOOSE = 0 # Fixed Orientation = 0, Loose = 1.
  49. DEPTH = 0 # Depth weighting - Best not to use for EMs
  50. DOWNSAMPLE = True # Whether to use the downsampled SRC or dense MRI mesh for calculating the GEM
  51. PICK_EEG = True # Whether to use EEG
  52. PICK_MEG = False # Whether to use MEG
  53. BINARISE_CEM = False # Whether to binarise the CEM before calculating the eigenmodes
  54. USE_FSAVERAGE_FWD = False # Whether to use the fsaverage forward model when solving the inverse problem
  55. USE_FSAVERAGE_GEM = False # Whether to use the fsaverage MRI when calculating the GEMs
  56. USE_GROUP_CON = False # Whether to use the group connectome when calculating the CEMs
  58. NUM_MODES = 1000
  59. STEP = 50
  60. num_modes = np.arange(start=STEP, stop=NUM_MODES+STEP, step=STEP)
  62. COMPARE_MNE = None # Whether to compare with MNE - Type the method or None/False
  63. SURROGATE_TEST = None # 'rotate' # 'roll' or 'rotate'
  64. SHUFFLE_ARRAY = np.arange(7, 700, 70) # Best to use an increment such that it's not a divisor of the number of sources in each hemisphere and don't start from 0 if using roll
  66. # Hyperparameters Unlikely To Change:
  67. ICO = 5 # Ideally this would be 5 for a higher resolution mesh, but 3 for this working example due to online memory limits.
  68. CON = 'con1019' # 'con1019' or 'con379'
  69. SFREQ = 78.125 # Hz - Default Sampling Frequency is 78.125Hz from 78.125Hz = 20kHz /256 = from 1/dt*1000 / scaling of 256 - Or it is /128 by the temporal averaging = 0.61 Hz
  71. # Path
  72. FS_PATH = Path("/Applications/freesurfer/7.3.2/")
  73. SUBJECTS_DIR = FS_PATH / "subjects"
  74. SUBJECT = "fsaverage"
  75. SIM_DIRECTORY = Path('Ignore/Simulations/')
  76. PARCELLATION = 'native.1000Parcels_Yeo2011_7Networks' # 'HCPMMP1'
  77. CON_PATH = f'External_Resources/HCP_data/HCP100206_{CON}'
  78. DISCRETE_FWD_PATH = Path(f'Ignore/{SUBJECT}_Con1000_twohemisrc_fwd.fif')
  80. # If None, will generate the GEMs - Path('Ignore/{SUBJECT}_ico5_gems.npz')
  81. # Better to recalculate to avoid mistakes, but using the ultra high res mesh is better and not realistic to recalculate.
  82. GEM_PATH = None
  84. SIM_TYPE = 'NMM'
  85. # '7networks_RH_Default_TempPole_1' # 'NMM' for TVB simualtion / a given time series or NAME_OF_REGION / INDEX_OF_REGION for a patch simulation
  86. EXTENT = 20 # Extent of the patch in mm
  87. SIM_FILE = Path(f'{SIM_DIRECTORY}/HCP100206_con1019_CIFTI_STRUCTURE_THALAMUS_LEFT_g_2.0.npz')
  88. print(f"Simulation File: {SIM_FILE}")
  90. # %%
  91. if SIM_TYPE == 'NMM':
  92. # Load the TVB Simulation
  93. raw_x = np.load(SIM_FILE)['raw_x'].astype(np.float32)
  94. raw_t = np.arange(raw_x.shape[0])/SFREQ
  96. # Reduce data to 60s
  97. raw_x = raw_x[:int(60*SFREQ),:,:,:]
  98. raw_t = raw_t[:int(60*SFREQ)]
  100. # Take only necessary state_variables
  101. raw_x = raw_x[:,0,:,0]
  102. # Transformations to the data:
  103. # Use the negative, so that the seizure-state has higher voltage
  104. raw_x = -raw_x
  105. # Add 1 to the data so neutral is 0V
  106. # raw_x = raw_x + 1
  107. # Lower the units
  108. raw_x = raw_x*10e-9
  110. # Subtracting by simple offset works best, such that resting is around 0
  111. offset = np.percentile(raw_x, 20)
  112. raw_x = raw_x - offset
  114. # %%
  115. if SIM_TYPE == 'NMM':
  116. plt.plot(raw_t, raw_x[:, 0])
  117. plt.xlabel('Time (s)', fontsize=14)
  118. plt.ylabel('Amplitude (Am)', fontsize=14)
  119. plt.xticks(fontsize=14)
  120. plt.yticks(fontsize=14)
  121. plt.tight_layout()
  122. plt.show()
  124. # %%
  125. if SIM_TYPE == 'NMM':
  126. # Plot a carpet plot (square)
  127. # Need to Change the Shape of the Simulation Output (2D and Transposed)
  128. plt.figure(figsize=(20, 10))
  129. cs=plt.imshow(raw_x.T, cmap="gray", aspect='equal',interpolation='none',origin='lower')
  130. plt.title('Raw Source Time Series', fontsize=20)
  131. axcb=plt.colorbar(cs)
  132. axcb.set_label('Amplitude (Am)', fontsize=20)
  133. axcb.ax.tick_params(labelsize=16)
  134. plt.ylabel('Region', fontsize=20)
  135. plt.xlabel('Time (s)', fontsize=20)
  136. xticks = np.arange(start=0,stop=raw_t.shape[0],step=int(np.ceil(raw_t.shape[0]/5)))
  137. plt.xticks(ticks=xticks,labels=np.round(raw_t[xticks],0).astype(int))
  138. plt.axis("tight")
  139. plt.xticks(fontsize=16)
  140. plt.yticks(fontsize=16)
  141. plt.show()
  143. # %% [markdown]
  144. # ## GEM Surface
  146. # %%
  147. if ANALYSIS_TYPE == 'GEM_surface':
  149. # Initialise the lists
  150. RLEs = []
  151. RLE_points = []
  152. coss = []
  153. ses = []
  154. cos_acrosstime = []
  155. expvars = []
  157. # Read fwd solution
  158. fwd = mne.read_forward_solution(f'Ignore/{SUBJECT}_ico{ICO}-fwd.fif') # mne.read_forward_solution(f'Ignore/{SUBJECT}_openmeeg_ico{ICO}-fwd.fif')
  159. # Construct the cortical mask
  160. mask = cortical_mask(SUBJECT, SUBJECTS_DIR, fwd['src'], downsample=True)
  161. info = generate_mne_info(sampling_freq=SFREQ, montage_name=MONTAGE, add_meg=PICK_MEG)
  162. if ONLY_20ELECTRODES:
  163. names_1020 = ['Fp1', 'Fp2', 'F7', 'F3', 'Fz', 'F4', 'F8', 'A1', 'T3', 'C3', 'Cz', 'C4', 'T4', 'T5', 'P3', 'Pz', 'P4', 'T6', 'O1', 'O2', 'A2']
  164. ch_inds = mne.pick_channels(ch_names=info['ch_names'], include=names_1020, )
  165. info = mne.pick_info(info, ch_inds)
  166. fwd = mne.pick_channels_forward(fwd, include=info['ch_names'], )
  167. fwd = mne.pick_types_forward(fwd, meg=PICK_MEG, eeg=PICK_EEG, exclude=info['bads'], )
  169. com_labels = None # If end up shifting the ordering of function calls around, com_labels is input to tvb_loadtomne
  170. # Load the TVB simulation
  171. stc, sim = tvb_loadtomne(raw_x, SUBJECT, fwd, info, PARCELLATION, selected_label=com_labels, fs_path= FS_PATH, )
  172. if SNR is not None:
  173. eeg_power = np.sum(sim.data**2)/sim.data.size
  174. eeg_noise = np.random.normal(0, np.sqrt(eeg_power)/np.sqrt(10**(SNR/10)), sim.data.shape)
  175. sim.data = sim.data + eeg_noise
  176. # Set EEG Reference if required - Doesn't matter if before or after adding noise
  177. if PICK_EEG:
  178. sim.set_eeg_reference('average',projection=True)
  180. # Compute covariance
  181. cov = mne.Covariance(data=np.identity(len(sim.ch_names)), names=sim.ch_names, bads=sim.info['bads'], projs=sim.info['projs'], nfree=1 )
  183. # Apply a highpass filter
  184. # sim = sim.filter(l_freq=1., h_freq=None )
  186. # Index of 50% of max EEG power
  187. hmp_i = np.where(np.sum(sim.data * sim.data, axis=0) > 0.5*np.sum(sim.data * sim.data, axis=0).max())[0][0]
  189. # GEM Solver
  190. # Graph for varying number of modes
  192. if USE_FSAVERAGE_FWD:
  193. fwd = mne.read_forward_solution(Path(f'Ignore/fsaverage_ico{ICO}-fwd.fif'))
  194. fwd = mne.pick_channels_forward(fwd, include=info['ch_names'], )
  195. fwd = mne.pick_types_forward(fwd, meg=PICK_MEG, eeg=PICK_EEG, exclude=info['bads'], )
  196. stc, _ = tvb_loadtomne(raw_x, 'fsaverage', fwd, info, PARCELLATION, selected_label=com_labels, fs_path= FS_PATH, )
  197. SUBJECT = 'fsaverage'
  199. if GEM_PATH:
  200. # Load the GEMs
  201. eval, emodes = np.load(GEM_PATH)['evals'], np.load(GEM_PATH)['emodes']
  202. else:
  203. # Calculate the eigenmodes
  204. if USE_FSAVERAGE_GEM:
  205. fwd_fsaverage = mne.read_forward_solution(Path(f'Ignore/fsaverage_ico{ICO}-fwd.fif'))
  206. mask = cortical_mask('fsaverage', SUBJECTS_DIR, fwd_fsaverage['src'], downsample=True)
  207. evals, emodes = both_hemi_calc_eig(fwd_fsaverage['src'], mask=mask, num_modes=NUM_MODES)
  208. del fwd_fsaverage
  209. else:
  210. evals, emodes = both_hemi_calc_eig(fwd['src'], mask=mask, num_modes=NUM_MODES)
  211. weight = 1 / (2*np.pi + np.sqrt(np.abs(evals)))
  212. emodes = emodes * weight[np.newaxis, :]
  214. for num_mode in num_modes:
  216. # 1. Obtain the eigenmodes
  217. emodes_subset = np.concatenate((emodes[:, :num_mode//2], emodes[:, NUM_MODES//2:NUM_MODES//2+num_mode//2]), axis=1)
  219. # 2. Apply the solver
  220. loose, depth = LOOSE, DEPTH
  221. solver_kwargs = dict(emodes=emodes_subset, rcond=RCOND)
  222. stc_est_region, var_exp, w = apply_solver(gem_solver, sim, fwd, cov, loose, depth, **solver_kwargs)
  224. # 3. Calculate threshold
  225. threshold = threshold_otsu(np.abs(stc_est_region.data,))/ np.max(np.abs(stc_est_region.data)) * 100
  226. # threshold = 40
  227. # print(f"Otsu's Threshold = {threshold}")
  229. # 4. Calculate metrics
  230. RLE = mne.simulation.metrics.region_localization_error(stc, stc_est_region, fwd['src'], threshold=f'{threshold}%', per_sample=False)
  231. RLE_point = mne.simulation.metrics.region_localization_error(stc.copy().crop(tmin=raw_t[hmp_i], tmax=raw_t[hmp_i]), stc_est_region.copy().crop(tmin=raw_t[hmp_i], tmax=raw_t[hmp_i]), fwd['src'], threshold=f'{threshold}%', per_sample=False)
  232. cos = mne.simulation.metrics.cosine_score(stc, stc_est_region, per_sample=False)
  233. aligned_stc_est = stc_aligner(stc, stc_est_region)
  234. nse = (aligned_stc_est.data/np.abs(aligned_stc_est.data).max() - stc.data/np.abs(stc.data).max())**2
  235. varexp_stc = log_exp_var(stc.data/np.abs(stc.data).max(), np.sqrt(nse))
  236. se = np.mean(nse)
  237. cos_t = mne.simulation.metrics.cosine_score(stc, stc_est_region, per_sample=True)
  238. varexp_stc = log_exp_var(stc.data/np.abs(stc.data).max(), np.sqrt(nse))
  240. expvars.append(varexp_stc)
  241. RLEs.append(RLE)
  242. RLE_points.append(RLE_point)
  243. coss.append(cos)
  244. ses.append(se)
  245. cos_acrosstime.append(cos_t)
  247. # %% [markdown]
  248. # ## CON1019 Discrete Source Space (1000 cortical regions)
  250. # %%
  251. if ANALYSIS_TYPE != 'GEM_surface':
  253. # Load Connectivity
  254. if SIM_TYPE == 'NMM':
  255. load_centres=False
  256. else:
  257. load_centres=True
  258. if USE_GROUP_CON or SUBJECT == 'fsaverage':
  259. con = con_loader(Path('External_Resources/HCP_data/SC_con1019/group_streamline_count_atlas.csv'), load_tracts=True, threshold=None, weights_only=True, load_centres=load_centres, subject=SUBJECT, subjects_dir=SUBJECTS_DIR)
  260. con.region_labels = np.loadtxt(CON_PATH / Path('region_names.txt'), dtype=str)
  261. else:
  262. con = con_loader(CON_PATH, load_tracts=True, threshold=None, weights_only=False, load_centres=load_centres, subject=SUBJECT, subjects_dir=SUBJECTS_DIR)
  263. selected_label = mne.read_labels_from_annot(SUBJECT, PARCELLATION, subjects_dir=SUBJECTS_DIR)
  264. # Remove the ??? regions
  265. selected_label = [label for label in selected_label if '???' not in label.name]
  267. # Load fwd, needs to be the discrete one
  268. fwd = mne.read_forward_solution(DISCRETE_FWD_PATH)
  269. info = generate_mne_info(sampling_freq=SFREQ, montage_name=MONTAGE, add_meg=PICK_MEG)
  270. if ONLY_20ELECTRODES:
  271. names_1020 = ['Fp1', 'Fp2', 'F7', 'F3', 'Fz', 'F4', 'F8', 'A1', 'T3', 'C3', 'Cz', 'C4', 'T4', 'T5', 'P3', 'Pz', 'P4', 'T6', 'O1', 'O2', 'A2']
  272. ch_inds = mne.pick_channels(ch_names=info['ch_names'], include=names_1020, )
  273. info = mne.pick_info(info, ch_inds)
  274. fwd = mne.pick_channels_forward(fwd, include=info['ch_names'], )
  275. fwd = mne.pick_types_forward(fwd, meg=PICK_MEG, eeg=PICK_EEG, exclude=info['bads'], )
  277. # 1. Generate the information from the dense source space Freesurfer labels
  278. region_centres = centre_loader(con.region_labels, coords='MRI_RAS', parcellation=PARCELLATION, subject=SUBJECT, subjects_dir=SUBJECTS_DIR, ) / 1000 # Convert to m
  279. label_hemi = np.array([label.hemi for label in selected_label])
  280. label_names = np.array([label.name[:-3] for label in selected_label])
  281. com_labels = generate_com_labels(region_centres, label_names, label_hemi, SUBJECT)
  283. if SIM_TYPE == 'NMM':
  284. # Load the TVB simulation and convert to EEG
  285. stc, sim = tvb_loadtomne(raw_x, SUBJECT, fwd, info, PARCELLATION, selected_label=com_labels, fs_path= FS_PATH, )
  286. else:
  287. # Create a single timepoint patch simulation with the centre region being the string in SIM_TYPE
  288. stc, sim = patch_simulator(con, SIM_TYPE, fwd, info, selected_label=com_labels, radius=EXTENT, restrict_hemi=True, full_src=True)
  289. raw_t = np.array([0])
  291. if SNR is not None:
  292. eeg_power = np.sum(sim.data**2)/sim.data.size
  293. eeg_noise = np.random.normal(0, np.sqrt(eeg_power)/np.sqrt(10**(SNR/10)), sim.data.shape)
  294. sim.data = sim.data + eeg_noise
  296. # Set EEG Reference if required - Doesn't matter if before or after adding noise
  297. if PICK_EEG:
  298. sim.set_eeg_reference('average',projection=True)
  299. # sim.apply_proj()
  301. # Compute covariance
  302. cov = mne.Covariance(data=np.identity(len(sim.ch_names)), names=sim.ch_names, bads=sim.info['bads'], projs=sim.info['projs'], nfree=1 )
  303. # Index of 50% of max EEG power
  304. hmp_i = np.where(np.sum(sim.data * sim.data, axis=0) > 0.5*np.sum(sim.data * sim.data, axis=0).max())[0][0]
  306. # %%
  307. # Load or generate the EMs
  308. if ANALYSIS_TYPE != 'GEM_surface':
  309. # In case we want to use a different fwd to the one used for the simulation (i.e. using group average fwd)
  310. if USE_FSAVERAGE_FWD:
  311. SUBJECT = 'fsaverage'
  312. fwd = mne.read_forward_solution(Path(f'Ignore/{SUBJECT}_Con1000_twohemisrc_fwd.fif'))
  314. # Read fwd solution
  315. fwd_surf = mne.read_forward_solution(Path(f'Ignore/{SUBJECT}_ico{ICO}-fwd.fif'))
  317. # Geometric eigenmodes
  318. if ANALYSIS_TYPE == 'GEM_discrete':
  319. # Load the pre-computed GEMs
  320. if GEM_PATH:
  321. evals, emodes = np.load(GEM_PATH)['evals'], np.load(GEM_PATH)['emodes']
  322. # Calculate the eigenmodes
  323. elif USE_FSAVERAGE_GEM:
  324. fwd_surf = mne.read_forward_solution(Path(f'Ignore/fsaverage_ico5-fwd.fif'))
  325. mask = cortical_mask('fsaverage', SUBJECTS_DIR, fwd_surf['src'], downsample=DOWNSAMPLE)
  326. evals, emodes = both_hemi_calc_eig(fwd_surf['src'], downsample=DOWNSAMPLE, mask=mask, num_modes=NUM_MODES, normalise=False )
  327. emodes = subset_eigenmodes(emodes, fwd_surf['src'], method='average', downsample=DOWNSAMPLE, parcellation=PARCELLATION, subjects_dir=SUBJECTS_DIR, region_labels=con.region_labels, subject='fsaverage', normalise=True)
  328. else:
  329. mask = cortical_mask(SUBJECT, SUBJECTS_DIR, fwd_surf['src'], downsample=DOWNSAMPLE)
  330. evals, emodes = both_hemi_calc_eig(fwd_surf['src'], downsample=DOWNSAMPLE, mask=mask, num_modes=NUM_MODES, normalise=False )
  331. # np.savez_compressed(Path(f'Ignore/{SUBJECT}_ico5_gems.npz'), evals=evals, emodes=emodes)
  332. emodes = subset_eigenmodes(emodes, fwd_surf['src'], method='average', downsample=DOWNSAMPLE, parcellation=PARCELLATION, subjects_dir=SUBJECTS_DIR, region_labels=con.region_labels, subject=SUBJECT, normalise=True)
  333. # np.savez_compressed(Path(f'Ignore/{SUBJECT}_con1000_gems.npz'), evals=evals, emodes=emodes)
  334. # Now weight the emodes by the inverse of the eigenvalues
  335. weight = 1 / (2*np.pi + np.sqrt(np.abs(evals)))
  336. emodes = emodes * weight[np.newaxis, :]
  338. # Connectivity eigenmodes
  339. elif ANALYSIS_TYPE == 'CEM_discrete':
  340. evals, emodes = calc_con_eigenmodes(con.weights, NUM_MODES, threshold=0, normed=True, binarise=BINARISE_CEM)
  341. # Now weight the emodes by the inverse of the eigenvalues
  342. weight = 1 / (0.1*np.abs(evals).mean() + np.abs(evals))
  343. emodes = emodes * weight[np.newaxis, :]
  344. # Truncate the eigenmodes to remove subcortical regions
  345. if CON == 'con1019':
  346. emodes = emodes[:-19,:]
  347. elif CON == 'con379':
  348. emodes = emodes[19:,:]
  350. # Both types of eigenmodes
  351. elif ANALYSIS_TYPE == 'SEM_discrete':
  352. mask = cortical_mask(SUBJECT, SUBJECTS_DIR, fwd_surf['src'], downsample=DOWNSAMPLE)
  353. g_evals, emodes = both_hemi_calc_eig(fwd_surf['src'], downsample=DOWNSAMPLE, mask=mask, num_modes=NUM_MODES, )
  354. g_emodes = subset_eigenmodes(emodes, fwd_surf['src'], method='average', downsample=DOWNSAMPLE, parcellation=PARCELLATION, subjects_dir=SUBJECTS_DIR, region_labels=con.region_labels, subject=SUBJECT)
  355. # Order g_emodes by eigenvalue
  356. # Create interleaved index array for ordering eigenmodes
  357. idx = np.array([[i, i+NUM_MODES//2] for i in range(1, NUM_MODES//2+1)]).flatten()
  358. idx = idx - 1 # Convert to 0-based indexing
  359. g_emodes = g_emodes[:, idx]
  360. g_evals = g_evals[idx]
  362. c_evals, c_emodes = calc_con_eigenmodes(con.weights, NUM_MODES, threshold=0, normed=True, binarise=BINARISE_CEM)
  363. c_emodes = c_emodes/ np.linalg.norm(c_emodes, axis=0)
  364. # Truncate the eigenmodes to remove subcortical regions
  365. if CON == 'con1019':
  366. c_emodes = c_emodes[:-19,:]
  367. elif CON == 'con379':
  368. c_emodes = c_emodes[19:,:]
  370. # # Normalise the evals
  371. g_evals = np.abs(g_evals)/np.abs(g_evals).mean()
  372. c_evals = np.abs(c_evals)/np.abs(c_evals).mean()
  374. g_weight = np.pi / (0.1 + np.sqrt(np.abs(g_evals)))
  375. c_weight = np.pi / (0.1 + c_evals)
  377. g_emodes = g_emodes * g_weight[np.newaxis, :]
  378. c_emodes = c_emodes * c_weight[np.newaxis, :]
  380. # Combine the emodes and evals
  381. evals = np.concatenate((g_evals, c_evals))
  382. idx = np.argsort(evals) # Get indices that would sort eigenvalues from low to high
  383. evals = evals[idx] # Sort eigenvalues
  384. emodes = np.concatenate((g_emodes, c_emodes), axis=1)[:, idx] # Combine and sort eigenmodes
  386. # %%
  387. if ANALYSIS_TYPE != 'GEM_surface':
  388. # Initialise the lists
  389. RLEs = []
  390. RLE_points = []
  391. coss = []
  392. ses = []
  393. cos_acrosstime = []
  394. expvars = []
  395. rle_acrosstime = []
  397. for num_mode in num_modes:
  398. # 1. Obtain the eigenmodes
  399. if ANALYSIS_TYPE == 'GEM_discrete':
  400. emodes_subset = np.concatenate((emodes[:, :num_mode//2], emodes[:, NUM_MODES//2:NUM_MODES//2+num_mode//2]), axis=1)
  401. evals_subset = np.concatenate((evals[:num_mode//2], evals[NUM_MODES//2:NUM_MODES//2+num_mode//2]), )
  402. # elif ANALYSIS_TYPE == 'SEM_discrete':
  403. # # This will make it so that there are an equal number of g and c modes
  404. # g_emodes_subset = g_emodes[:, :num_mode//2]
  405. # c_emodes_subset = c_emodes[:, :num_mode//2]
  406. # emodes_subset = np.zeros((g_emodes_subset.shape[0], num_mode))
  407. # emodes_subset[:, ::2] = g_emodes_subset # Even indices for g_emodes
  408. # emodes_subset[:, 1::2] = c_emodes_subset # Odd indices for c_emodes
  409. else:
  410. emodes_subset = emodes[:, :num_mode]
  411. evals_subset = evals[:num_mode]
  413. # 2. Apply the solver
  414. loose, depth = LOOSE, DEPTH
  415. solver_kwargs = dict(emodes=emodes_subset, rcond=RCOND,)
  416. stc_est_region, var_exp, w = apply_solver(gem_solver, sim, fwd, cov, loose, depth, pick_ori='normal', **solver_kwargs)
  418. # 3. Calculate threshold
  419. threshold = threshold_otsu(np.abs(stc_est_region.data,))/ np.max(np.abs(stc_est_region.data)) * 100
  420. # threshold = 40
  421. print(f"Otsu's Threshold = {threshold}")
  423. # 4. Calculate metrics
  424. RLE = mne.simulation.metrics.region_localization_error(stc, stc_est_region, fwd['src'], threshold=f'{threshold}%', per_sample=False)
  425. rle_t = mne.simulation.metrics.region_localization_error(stc, stc_est_region, fwd['src'], threshold=f'{threshold}%', per_sample=True)
  426. RLE_point = mne.simulation.metrics.region_localization_error(stc.copy().crop(tmin=raw_t[hmp_i], tmax=raw_t[hmp_i]), stc_est_region.copy().crop(tmin=raw_t[hmp_i], tmax=raw_t[hmp_i]), fwd['src'], threshold=f'{threshold}%', per_sample=False)
  427. cos = mne.simulation.metrics.cosine_score(stc, stc_est_region, per_sample=False)
  428. aligned_stc_est = stc_aligner(stc, stc_est_region) # Technically not needed for the discrete source case
  429. nse = (aligned_stc_est.data/np.abs(aligned_stc_est.data).max() - stc.data/np.abs(stc.data).max())**2
  430. varexp_stc = log_exp_var(stc.data/np.abs(stc.data).max(), np.sqrt(nse))
  431. se = np.mean(nse)
  432. cos_t = mne.simulation.metrics.cosine_score(stc, stc_est_region, per_sample=True)
  434. expvars.append(varexp_stc)
  435. RLEs.append(RLE)
  436. RLE_points.append(RLE_point)
  437. coss.append(cos)
  438. ses.append(se)
  439. cos_acrosstime.append(cos_t)
  440. rle_acrosstime.append(rle_t)
  442. # %% [markdown]
  443. # # Graphs
  445. # %%
  446. # Plot the results
  447. # Turn into arrays
  448. RLEs = np.array(RLEs)
  449. RLE_points = np.array(RLE_points)
  451. fig, axs = plt.subplots(3, 2, figsize=(20, 10))
  453. axs[0, 0].plot(num_modes[:len(RLEs)], RLEs*1e3, marker='o')
  454. axs[0, 0].axvline(fwd['nchan'], ymin=0, ymax=1, color='r', linestyle='--')
  455. axs[0, 0].set_title(f'RLE - Otsu\'s Threshold')
  456. axs[0, 0].set_xlabel('Number of Modes')
  457. axs[0, 0].set_ylabel(f'RLE (mm)')
  458. axs[0, 0].grid()
  460. axs[0, 1].plot(num_modes, RLE_points*1e3, marker='o')
  461. axs[0, 1].axvline(fwd['nchan'], ymin=0, ymax=1, color='r', linestyle='--')
  462. axs[0, 1].set_title(f'RLE 50% EEG at Otsu\'s Threshold - t = {raw_t[hmp_i]:.1f}s')
  463. axs[0, 1].set_xlabel('Number of Modes')
  464. axs[0, 1].set_ylabel('RLE 50% EEG (mm)')
  465. axs[0, 1].grid()
  467. # Plot the Explained Variance
  468. axs[1,0].plot(num_modes, expvars, marker='o')
  469. axs[1,0].axvline(fwd['nchan'], ymin=0, ymax=1, color='r', linestyle='--')
  470. axs[1,0].set_title('Explained Variance')
  471. axs[1,0].set_xlabel('Number of Modes')
  472. axs[1,0].set_ylabel('Explained Variance (%)')
  473. axs[1,0].grid()
  475. # Normalised Mean Squared Error
  476. axs[1, 1].plot(num_modes, ses, marker='o')
  477. axs[1, 1].axvline(fwd['nchan'], ymin=0, ymax=1, color='r', linestyle='--')
  478. # axs[1, 1].axvline(306, ymin=0, ymax=1, color='r', linestyle='--')
  479. # axs[1, 1].axvline(343, ymin=0, ymax=1, color='r', linestyle='--')
  480. axs[1, 1].set_title('Normalised Mean Squared Error')
  481. axs[1, 1].set_xlabel('Number of Modes')
  482. axs[1, 1].set_ylabel('Mean Squared Error')
  483. # axs[1, 1].set_yscale('log')
  484. axs[1, 1].grid()
  486. axs[2, 0].plot(num_modes, coss, marker='o')
  487. axs[2, 0].axvline(fwd['nchan'], ymin=0, ymax=1, color='r', linestyle='--')
  488. axs[2, 0].set_title('Cosine Score')
  489. axs[2, 0].set_xlabel('Number of Modes')
  490. axs[2, 0].set_ylabel('Cosine Score')
  491. axs[2, 0].grid()
  493. cos_acrosstime = np.array(cos_acrosstime)
  494. plot_mean_and_range(cos_acrosstime, raw_t, stat_type='max', show_minmax=True, ax=axs[2,1])
  495. axs[2,1].set_title(f'Cosine Score Across Time - {ANALYSIS_TYPE}')
  496. axs[2,1].set_xlabel('Time (s)')
  497. axs[2,1].set_ylabel('Cosine Score')
  499. plt.tight_layout()
  501. # %% [markdown]
  502. # ## MNE inverse operator for comparison
  504. # %%
  505. COMPARE_MNE = 'sLORETA'
  506. # MNE inverse operator for comparison
  507. if COMPARE_MNE:
  508. # Compute inverse operator
  509. inv = make_inverse_operator(sim.info, fwd, cov, loose=LOOSE, depth=DEPTH, )
  510. if SNR is not None:
  511. # Noting that MNE-Python takes in a linear SNR - (10**(SNR/10)) - 1 / lambda2**2
  512. stc_est_region_mne = apply_inverse(sim, inv, method=COMPARE_MNE, lambda2=1/(10**(SNR/10)), pick_ori=None)
  513. else:
  514. stc_est_region_mne = apply_inverse(sim, inv, method=COMPARE_MNE, lambda2=np.finfo(np.float64).eps, pick_ori=None)
  516. # 3. Calculate tresholds
  517. threshold = threshold_otsu(np.abs(stc_est_region_mne.data,))/ np.max(np.abs(stc_est_region_mne.data)) * 100
  518. # threshold = 40
  519. print(f"Otsu's Threshold = {threshold}")
  521. print('RLE timeseries = ', mne.simulation.metrics.region_localization_error(stc, stc_est_region_mne, fwd['src'], threshold=f'{threshold}%', per_sample=False))
  522. print('RLE 50% EEG = ', mne.simulation.metrics.region_localization_error(stc.copy().crop(tmin=raw_t[hmp_i], tmax=raw_t[hmp_i]), stc_est_region_mne.copy().crop(tmin=raw_t[hmp_i], tmax=raw_t[hmp_i]), fwd['src'],
  523. threshold=f'{threshold}%', per_sample=False))
  524. print('Cosine Score = ', mne.simulation.metrics.cosine_score(stc, stc_est_region_mne, per_sample=False))
  525. cos_t = mne.simulation.metrics.cosine_score(stc, stc_est_region_mne, per_sample=True)
  526. cos_t = cos_t.reshape(1,-1)
  527. aligned_stc_est = stc_aligner(stc, stc_est_region_mne)
  528. print('Squared Error = ', np.sum((aligned_stc_est.data - stc.data)**2))
  529. print('Standardised Square Error = ', np.sum((aligned_stc_est.data/np.abs(aligned_stc_est.data).max() - stc.data/np.abs(stc.data).max())**2))
  530. print('Variance Explained = ', log_exp_var(stc.data/np.abs(stc.data).max(), np.sqrt(np.sum((aligned_stc_est.data/np.abs(aligned_stc_est.data).max() - stc.data/np.abs(stc.data).max())**2))))
  531. # print('Correlation = ', np.corrcoef(aligned_stc_est.data, stc.data)[0, 1]) # Correlation takes too long to compute for ico5
  533. # %%
  534. if COMPARE_MNE:
  535. fig, axs = plt.subplots(3, 2, figsize=(20, 10))
  536. plot_mean_and_range(cos_t, raw_t, show_minmax=True, ax=axs[0,0])
  537. axs[0,0].set_title(f'Cosine Score Across Time - {COMPARE_MNE}', fontsize=12)
  538. axs[0,0].set_xlabel('Time (s)', fontsize=12)
  539. axs[0,0].set_ylabel('Cosine Score', fontsize=12)
  540. axs[0,0].set_title(f'Cosine Score Across Time - {COMPARE_MNE}', fontsize=12)
  541. axs[0,0].set_xlabel('Time (s)', fontsize=12)
  542. axs[0,0].set_ylabel('Cosine Score', fontsize=12)
  543. axs[0,0].get_legend().remove()
  544. plt.tight_layout()
  546. # %% [markdown]
  547. # ## Surrogate Test
  549. # %%
  550. if 'SHUFFLE_ARRAY' not in locals() and SURROGATE_TEST:
  551. # Create a shuffle array - Increments of 7 such that it's not a divisor of the number of sources in each hemisphere
  552. print('Shuffle Array not defined, using default with 99 surrogates')
  553. SHUFFLE_ARRAY = np.arange(7, 700, 7)
  555. # Perform the surrogate test for N interations
  556. if SURROGATE_TEST:
  558. RLEs = RLEs.tolist()
  559. RLE_points = RLE_points.tolist()
  561. for shuffle in SHUFFLE_ARRAY:
  562. # Generate Surrogate Modes using the Roll Method
  563. if SURROGATE_TEST == 'roll':
  564. if ANALYSIS_TYPE == 'GEM_discrete':
  565. emodes_subset = np.concatenate((emodes[:, :num_mode//2], emodes[:, NUM_MODES//2:NUM_MODES//2+num_mode//2]), axis=1)
  566. elif ANALYSIS_TYPE == 'CEM_discrete':
  567. emodes_subset = emodes[:, :num_mode]
  569. # Roll the eigenmodes
  570. emodes_subset[:emodes_subset.shape[0]//2,:emodes_subset.shape[1]//2] = np.roll(emodes_subset[:emodes_subset.shape[0]//2,:emodes_subset.shape[1]//2], shift=shuffle, axis=0)
  571. emodes_subset[emodes_subset.shape[0]//2:,emodes_subset.shape[1]//2:] = np.roll(emodes_subset[emodes_subset.shape[0]//2:,emodes_subset.shape[1]//2:], shift=shuffle, axis=0)
  573. # Generate Surrogate Modes using the Rotational Eigenstrapping Method
  574. if SURROGATE_TEST == 'rotate':
  575. # Normalise to be safe
  576. emodes = emodes / np.linalg.norm(emodes, axis=0)[np.newaxis, :]
  577. surrogate_emodes = both_hemi_surrogate_eig(emodes, evals, normalize=True, thresh=0, )
  578. surrogate_emodes = surrogate_emodes*weight[np.newaxis, :]
  579. emodes_subset = np.concatenate((surrogate_emodes[:, :num_mode//2], surrogate_emodes[:, NUM_MODES//2:NUM_MODES//2+num_mode//2]), axis=1)
  581. # 2. Apply the solver
  582. loose, depth = LOOSE, DEPTH
  583. solver_kwargs = dict(emodes=emodes_subset, rcond=RCOND)
  584. stc_est_region, var_exp, w = apply_solver(gem_solver, sim, fwd, cov, loose, depth, **solver_kwargs)
  586. # 3. Calculate threshold
  587. threshold = threshold_otsu(np.abs(stc_est_region.data), )/ np.max(stc_est_region.data) * 100
  589. # 4. Calculate metrics
  590. RLE = mne.simulation.metrics.region_localization_error(stc, stc_est_region, fwd['src'], threshold=f'{threshold}%', per_sample=False)
  591. RLE_point = mne.simulation.metrics.region_localization_error(stc.copy().crop(tmin=raw_t[hmp_i], tmax=raw_t[hmp_i]), stc_est_region.copy().crop(tmin=raw_t[hmp_i], tmax=raw_t[hmp_i]), fwd['src'], threshold=f'{threshold}%', per_sample=False)
  592. cos = mne.simulation.metrics.cosine_score(stc, stc_est_region, per_sample=False)
  593. aligned_stc_est = stc_aligner(stc, stc_est_region) # Technically not needed for the discrete source case
  594. nse = (aligned_stc_est.data/np.abs(aligned_stc_est.data).max() - stc.data/np.abs(stc.data).max())**2
  595. varexp_stc = log_exp_var(stc.data/np.abs(stc.data).max(), np.sqrt(nse))
  596. se = np.mean(nse)
  598. expvars.append(varexp_stc)
  599. RLEs.append(RLE)
  600. RLE_points.append(RLE_point)
  601. coss.append(cos)
  602. ses.append(se)
  604. # Clear Cell Output
  605. clear_output(wait=True)
  607. RLEs = np.array(RLEs)
  608. RLE_points = np.array(RLE_points)
  610. # %%
  611. # Stop it here for now if run all
  612. 1/0
  614. # %% [markdown]
  615. # ## Videos of Sources
  616. #
  617. # ### Source Estimate
  619. # %%
  620. if ANALYSIS_TYPE != 'GEM_surface':
  621. # To convert the discrete source estimate to a surface estimate just run these two functions - then allows us to use MNE's plotting functions
  622. fwd_surface = mne.read_forward_solution(f'Ignore/{SUBJECT}_ico3-fwd.fif')
  623. fwd_surface = mne.pick_channels_forward(fwd_surface, include=info['ch_names'], )
  624. fwd_surface = mne.pick_types_forward(fwd_surface, meg=PICK_MEG, eeg=PICK_EEG, exclude=info['bads'], )
  625. stc_surface, _ = tvb_loadtomne(stc_est_region.data.T, SUBJECT, fwd_surface, info, PARCELLATION, fs_path= FS_PATH, )
  626. # time_max here is the frame shown in the video. Which can be interpretted as a picture.
  627. surfer_kwargs = dict(
  628. hemi='both', subjects_dir=SUBJECTS_DIR,
  629. # colormap='mne',
  630. # clim=dict(kind='value', lims=[stc_est_region.data.min(), (stc_est_region.data.max()+stc_est_region.data.min())/2, stc_est_region.data.max()]),
  631. clim = dict(kind='value', lims=[0, stc_est_region.data.max()/2, stc_est_region.data.max()]),
  632. # clim=dict(kind='value', lims=[0, stc.data.max()/2, stc.data.max()]),
  633. views='rostral', initial_time=raw_t[hmp_i],
  634. time_unit='s', size=(800, 800), smoothing_steps=10,
  635. backend='pyvistaqt',
  636. )
  638. if ANALYSIS_TYPE == 'GEM_surface':
  639. brain = stc_est_region.plot(**surfer_kwargs)
  640. else:
  641. brain = stc_surface.plot(**surfer_kwargs)
  643. brain.set_playback_speed(speed=1)
  645. # %% [markdown]
  646. # ### True Sources
  648. # %%
  649. if ANALYSIS_TYPE != 'GEM_surface':
  650. # Looking at true stc
  651. stc_surface, _ = tvb_loadtomne(stc.data.T, SUBJECT, fwd_surface, info, PARCELLATION, fs_path= FS_PATH, )
  653. # time_max here is the frame shown in the video. Which can be interpretted as a picture.
  654. surfer_kwargs = dict(
  655. hemi='both', subjects_dir=SUBJECTS_DIR,
  656. # colormap='mne',
  657. # clim=dict(kind='value', lims=[stc.data.min(), (stc.data.max()+stc.data.min())/2, stc.data.max()]),
  658. clim = dict(kind='value', lims=[0, stc.data.max()/2, stc.data.max()]),
  659. # clim=dict(kind='value', lims=[0, stc.data[:,hmp_i].max()/2, stc.data[:,hmp_i].max()]),
  660. views='rostral',
  661. initial_time=raw_t[hmp_i], #
  662. # background='white',
  663. time_unit='s', size=(800, 800), smoothing_steps=10)
  665. if ANALYSIS_TYPE == 'GEM_surface':
  666. brain = stc.plot(**surfer_kwargs)
  667. else:
  668. brain = stc_surface.plot(**surfer_kwargs)
  669. brain.set_playback_speed(speed=1)
  671. # %% [markdown]
  672. # ## Visualising the Eigenmodes
  674. # %%
  675. # Plot the eigenmodes
  676. l_points, l_trias = get_tria_from_src(fwd_surf['src'], hemi=0, downsample=DOWNSAMPLE)
  677. l_tria = lapy.TriaMesh(l_points, l_trias)
  678. mask = cortical_mask(SUBJECT, SUBJECTS_DIR, fwd_surf['src'], downsample=DOWNSAMPLE)
  679. evals, s_emodes = both_hemi_calc_eig(fwd_surf['src'], downsample=DOWNSAMPLE, mask=mask, num_modes=NUM_MODES, normalise=False )
  681. mode_number = 20
  682. plot_tria_mesh(l_tria, vfunc=s_emodes[:s_emodes.shape[0]//2, mode_number],showcaxis=True)
  684. # %%
  685. from mpl_toolkits.mplot3d import Axes3D
  686. %matplotlib widget
  688. # Check discrete source space - for the eigenmodes
  690. # 3D plot of Discrete Source space where colour is dependent on subset_emodes
  692. mode_number = 20
  694. # Extract positions and colors
  695. positions = np.concatenate((fwd['src'][0]['rr'], fwd['src'][1]['rr']), axis=0)
  696. # colors = disc_emodes[:, mode_number]
  697. colors = emodes[:, mode_number]
  700. # Create a 3D scatter plot
  701. fig = plt.figure(figsize=(10, 10))
  702. ax = fig.add_subplot(111, projection='3d')
  703. sc = ax.scatter(positions[:, 0], positions[:, 1], positions[:, 2], c=colors, cmap='coolwarm', s=50)
  704. # Ensure axes follow the same scale
  705. ax.set_box_aspect([np.ptp(positions[:, 0]), np.ptp(positions[:, 1]), np.ptp(positions[:, 2])])
  707. # Add color bar
  708. cb = plt.colorbar(sc, ax=ax, shrink=0.5)
  709. cb.set_label(f'Subset Eigenmodes (Mode {mode_number})', fontsize=12)
  711. # Set labels
  712. ax.set_xlabel('X', fontsize=12)
  713. ax.set_ylabel('Y', fontsize=12)
  714. ax.set_zlabel('Z', fontsize=12)
  715. ax.set_title('3D Plot of Discrete Source Space', fontsize=14)
  717. plt.show()
  719. # %%

STC_Analysis.ipynb at commit 81268ed, under BSD-3-Clause · at the source

Overview

Authors: Pok Him Siu1,2, Philippa J. Karoly1,2, Sina Mansour L.1,3, Artemio Soto‐Breceda4, Levin Kuhlmann5, Mark J. Cook1,2,6, David B. Grayden1,2,6
ORCID iDs: Pok Him Siu, Philippa J. Karoly, Sina Mansour L., Artemio Soto‐Breceda, Levin Kuhlmann, Mark J. Cook, David B. Grayden
  1. Department of Biomedical Engineering The University of Melbourne Melbourne Australia
  2. Graeme Clark Institute The University of Melbourne Melbourne Australia
  3. National University of Singapore Singapore Singapore
  4. Grenoble Institut des Neurosciences INSERM Grenoble (La Tronche) France
  5. Department of Data Science & AI Monash University Melbourne Australia
  6. Department of Medicine, St Vincent's Hospital, Melbourne The University of Melbourne Melbourne Australia
Institutions: The University of Melbourne (Australia); National University of Singapore (Singapore); Inserm (France); Grenoble Institute of Neurosciences (France); Monash University (Australia); St Vincent's Hospital Melbourne (Australia)
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany), volume 13, issue 31, article e16802
Dates: received 1 September 2025; accepted 12 March 2026; published online 23 March 2026; in print June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1002/advs.202516802 · PMID 41869755 · PMCID PMC13588046 · OpenAlex W7140078093
Open access: gold, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: EEG (modality), human (organism), epilepsy (population)
Methods: Connectivity, Statistics, Physiology & signal measures, Machine learning, fMRI & imaging
Keywords: connectome, EEG, epilepsy, geometric eigenmodes, source localization
MeSH: Brain*, Brain Mapping*, Connectome*, Electroencephalography*, Epilepsy*, Humans, Magnetic Resonance Imaging (* major topic)
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Department of Education and Training | Australian Research Council (ARC) (DP200102600)
Citations: cited by 1 paper (Europe PMC); 82 references in the paper

Abstract

A fundamental view of neuroscience is that, in addition to neuronal activity, the structure of the brain constrains and explains brain function. An alluring formalism in computational neuroscience has been the generation of structural eigenmodes of neural activity from a matrix representing the anatomy of the brain. Traditionally, brain connectomics has been the gold standard for the coupling between structure and function. However, it has recently been suggested that simpler brain geometry can provide more explanatory power in fMRI. An adjacent modality is the source localization problem of EEG, which aims to identify the underlying generators of EEG recordings. The underdetermined nature of the problem requires sufficient constraints to produce realistic and unique solutions of source activity. In this work, a simple framework for incorporating different forms of structural brain eigenmodes to constrain the source localization problem in epilepsy is presented. Geometric eigenmodes were found to be able to reconstruct the spread of a seizure through the brain slightly better than connectome eigenmodes, and both types of structural modes significantly outperformed commonly used approaches.

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 14 matches between paragraphs and lines of code.

Spokhim/sesl

License: BSD-3-Clause
State: the link answers, verified on 31 October 2026
Evidence: files inventoried
Commit: 81268edac5a70ea091d06cdce6ec2539075fdec0, 23 March 2026
Languages: Python (4), Jupyter (2)
Size: 11 files, 6 scripts
Software Heritage: not archived
Found in: “Data Availability Statement”
Holds: README, license file, environment (requirements.txt, requirements_sesl.txt), 2 notebooks
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (6 files), MNE-Python (5 files), Matplotlib (3 files), SciPy (3 files), BrainSpace (1 file), Plotly (1 file), scikit-image (1 file), scikit-learn (1 file), The Virtual Brain (1 file)
Availability: 2 checks, the latest on 31 October 2026: the link answers
  • 31 October 2026: the link answers
  • 30 September 2026: the link answers
8 files
At the source: github.com/Spokhim/sesl

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:

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

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

Data

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

Data Availability Statement

All code for the source localization pipeline can be found at: https://github.com/Spokhim/sesl. All code should be easily compatible with current MNE‐Python pipelines whether it is EEG or MEG. If there is enough interest, the author may look into a proper implementation into MNE‐Python.

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

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 5 keywords, 7 MeSH terms, 1 funder, 73 references.

Cite

This paper

Siu, P. H., Karoly, P. J., Mansour L., S., Soto‐Breceda, A., Kuhlmann, L., Cook, M. J., & Grayden, D. B. (2026). Structural Eigenmodes of the Brain to Improve the Source Localization of EEG: Application to Epileptiform Activity. Advanced science (Weinheim, Baden-Wurttemberg, Germany), 13(31), e16802. https://doi.org/10.1002/advs.202516802

BibTeX

@article{siu2026structural,
author = {Siu, Pok Him and Karoly, Philippa J. and Mansour L., Sina and Soto‐Breceda, Artemio and Kuhlmann, Levin and Cook, Mark J. and Grayden, David B.},
title = {{Structural Eigenmodes of the Brain to Improve the Source Localization of EEG: Application to Epileptiform Activity}},
journal = {Advanced science (Weinheim, Baden-Wurttemberg, Germany)},
year = {2026},
month = mar,
volume = {13},
number = {31},
pages = {e16802},
publisher = {Wiley},
issn = {2198-3844},
doi = {10.1002/advs.202516802},
url = {https://doi.org/10.1002/advs.202516802},
pmid = {41869755},
pmcid = {PMC13588046}
}

RIS

TY - JOUR
AU - Siu, Pok Him
AU - Karoly, Philippa J.
AU - Mansour L., Sina
AU - Soto‐Breceda, Artemio
AU - Kuhlmann, Levin
AU - Cook, Mark J.
AU - Grayden, David B.
TI - Structural Eigenmodes of the Brain to Improve the Source Localization of EEG: Application to Epileptiform Activity
T2 - Advanced science (Weinheim, Baden-Wurttemberg, Germany)
J2 - Adv Sci (Weinh)
PY - 2026
DA - 2026/03/23
VL - 13
IS - 31
SP - e16802
SN - 2198-3844
PB - Wiley
DO - 10.1002/advs.202516802
UR - https://doi.org/10.1002/advs.202516802
LA - en
ER -

CSL-JSON

{
"id": "10.1002/advs.202516802",
"type": "article-journal",
"title": "Structural Eigenmodes of the Brain to Improve the Source Localization of EEG: Application to Epileptiform Activity",
"container-title": "Advanced science (Weinheim, Baden-Wurttemberg, Germany)",
"author": [
{
"family": "Siu",
"given": "Pok Him"
},
{
"family": "Karoly",
"given": "Philippa J."
},
{
"family": "Mansour L.",
"given": "Sina"
},
{
"family": "Soto‐Breceda",
"given": "Artemio"
},
{
"family": "Kuhlmann",
"given": "Levin"
},
{
"family": "Cook",
"given": "Mark J."
},
{
"family": "Grayden",
"given": "David B."
}
],
"container-title-short": "Adv Sci (Weinh)",
"volume": "13",
"issue": "31",
"page": "e16802",
"DOI": "10.1002/advs.202516802",
"PMID": "41869755",
"PMCID": "PMC13588046",
"ISSN": "2198-3844",
"publisher": "Wiley",
"URL": "https://doi.org/10.1002/advs.202516802",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
23
]
]
}
}

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

Similar papers

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

[1] doi:10.1162/netn.a.543 [code]
High-resolution Bayesian Virtual Epileptic Patient using neural field models.
Journal: Network neuroscience (Cambridge, Mass.)
In common: The Virtual Brain, MNE-Python, SciPy, 2 other tools, epilepsy, 10 references
[2] doi:10.1002/advs.202405246 [code]
Seizure Sources Can Be Imaged from Scalp EEG by Means of Biophysically Constrained Deep Neural Networks.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: The Virtual Brain, MNE-Python, scikit-image, 3 other tools, epilepsy, EEG, 9 references
[3] doi:10.1073/pnas.2201128119 [code]
Deep neural networks constrained by neural mass models improve electrophysiological source imaging of spatiotemporal brain dynamics.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: The Virtual Brain, MNE-Python, scikit-image, 3 other tools, epilepsy, EEG, 8 references
[4] doi:10.1371/journal.pcbi.1011642 [code]
Parameter estimation in a whole-brain network model of epilepsy: Comparison of parallel global optimization solvers.
Journal: PLoS computational biology
In common: scikit-learn, SciPy, Matplotlib, 1 other tool, epilepsy, EEG, 10 references
[5] doi:10.1038/s43588-025-00841-6 [code]
Virtual brain twins for stimulation in epilepsy.
Journal: Nature computational science
In common: The Virtual Brain, MNE-Python, scikit-learn, 3 other tools, epilepsy, EEG, 7 references
[6] doi:10.1371/journal.pcbi.1012911 [code]
Virtual epilepsy patient cohort: Generation and evaluation.
Journal: PLoS computational biology
In common: The Virtual Brain, MNE-Python, scikit-learn, 3 other tools, epilepsy, EEG, 6 references
[7] doi:10.1093/braincomms/fcag217 [code]
Investigation of stress hormones across multiday seizure cycles.
Journal: Brain communications
In common: epilepsy, 3 authors
[8] doi:10.7554/elife.106194 [code]
Virtual Brain Inference (VBI), a flexible and integrative toolkit for efficient probabilistic inference on whole-brain models.
Journal: eLife
In common: MNE-Python, scikit-learn, SciPy, 2 other tools, 8 references
[9] doi:10.1038/s41586-023-06098-1 [code]
Geometric constraints on human brain function.
Journal: Nature
In common: BrainSpace, SciPy, NumPy, 10 references
[10] doi:10.1093/braincomms/fcac173 [code]
Multiple mechanisms shape the relationship between pathway and duration of focal seizures.
Journal: Brain communications
In common: epilepsy, 3 references, 2 authors

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.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1002/advs.202516802.svg)](https://openscicode.org/badge/doi/10.1002/advs.202516802)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1002/advs.202516802.svg
   :target: https://openscicode.org/badge/doi/10.1002/advs.202516802
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1002/advs.202516802"><img src="https://openscicode.org/badge/doi/10.1002/advs.202516802.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge 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: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

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