SEEG Time-Frequency Fingerprint Analysis for Epileptogenic Zone Localization

Authors: Yash Shashank Vakilna, Chinmay Chinara, Johnson Hampson, Takfarinas Medani, Raymundo Cassani, John Mosher, Sylvain Baillet, Richard Leahy

fc_3.png

This tutorial guides users through computing time-frequency decomposition maps to identify the epileptogenic zone (EZ) using both ictal and interictal SEEG data.

Make sure you complete the previous tutorials on CT to MRI co-registration and Contact localization and labeling before proceeding any further.

Access the recordings

Import the contacts positions

In order to generate epileptogenicity maps, we need accurate 3D positions for the contacts of the depth electrodes. Placing the contacts requires a good understanding of the implantation scheme reported by the neurosurgeon, and some skills in reading MRI scans.

Display the depth electrodes

3D figures

MRI Viewer

Panel iEEG

To know more about ways to display the SEEG recordings in Brainstorm refer to the Epileptogenicity tutorial.

Review recordings

Power spectrum

We recommend that you start your data analysis with a power spectral density estimation of the recordings to check the quality of sensor recording. This is described in more details in the Power spectrum tutorial.

Add events

We need to mark seizure onset event for the ictal and LVFA and wave recordings and spike event for interictal recording. There are events already available in recordings, that were marked for clinical use, to jump quickly to the page of interest. More details can be found in the tutorial Event Markers.

Import epochs of interest

At this point of the analysis, we are still looking at the original files, no SEEG data was copied to the database. The montages are saved in the Brainstorm preferences, the new events are saved in the links of the database.

We are now going to import the three segment of recordings i.e. LVFA and wave , ictal repetitive spike and interictal spike which are a subset of the Baseline recording.

Import in database

Bipolar montage

We will run the rest of the analysis using a bipolar montage (bipolar-2). The montage selected in the Record tab is for visualization only, most processes ignore this selection and work only on the original common-referential montage. To compute bipolar montage on time series, we need to explicitly apply the montage to the recordings. More details can be found in tutorials Montage editor and Epileptogenicity.

Head modeling

The forward models depend on the subject's anatomy, including head size and geometry, tissue conductivity, the computational method, and sensor characteristics. In this section, we will use the Boundary Element Method (BEM) approach available in Brainstorm for constructing the head model for sEEG.

Compute noise covariance matrix

Modeling interictal spikes

Compute inverse model

Display sensor time series

View sources

Atlases and scouts

Modeling ictal wave

Switch to the folder LVFA_and_wave (not the RAW folder) and repeat the steps to compute inverse model as per the section above and study the sensor time series and inverse modeling results.

69_disp_ts_wave.png

70_view_inv_model_wave.png

Modeling ictal onset with LVFA

Sensor space

Compute time-frequency decomposition

View time-frequency maps

Source space

Extract scout time series

Compute time-frequency decomposition

View time-frequency maps

Modeling ictal onset with repetitive spiking

Sensor space

Display time-series

Compute time-frequency decomposition

Same as above sections. Only one change, set sensor type: PIN5-PIN6 for Time-frequency (Morlet wavelets) process.

78_time_freq_decomp_ictal.png

View time-frequency maps

Source space

Compute inverse model

View sources

Additional Documentation

Forum discussions

Articles

Advanced

Scripting

The following script from the Brainstorm distribution reproduces the analysis presented in this tutorial page: brainstorm3/toolbox/script/tutorial_seizure_fingerprinting.m

1 function tutorial_seizure_fingerprinting(tutorial_dir, reports_dir) 2 % TUTORIAL_SEIZURE_FINGERPRINTING: Script that reproduces the results of the online tutorial "Seizure Fingerprinting". 3 % 4 % CORRESPONDING ONLINE TUTORIALS: 5 % https://neuroimage.usc.edu/brainstorm/Tutorials/SeizureFingerprinting 6 % 7 % INPUTS: 8 % - tutorial_dir : Directory where the tutorial_seizure_fingerprinting.zip file has been unzipped 9 % - reports_dir : Directory where to save the execution report (instead of displaying it) 10 11 % @============================================================================= 12 % This function is part of the Brainstorm software: 13 % https://neuroimage.usc.edu/brainstorm 14 % 15 % Copyright (c) University of Southern California & McGill University 16 % This software is distributed under the terms of the GNU General Public License 17 % as published by the Free Software Foundation. Further details on the GPLv3 18 % license can be found at http://www.gnu.org/copyleft/gpl.html. 19 % 20 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE 21 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY 22 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF 23 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY 24 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. 25 % 26 % For more information type "brainstorm license" at command prompt. 27 % =============================================================================@ 28 % 29 % Authors: Chinmay Chinara, 2025 30 % Yash Shashank Vakilna, 2025 31 % Raymundo Cassani, 2025 32 33 %% ===== PARSE INPUTS ===== 34 % Output folder for reports 35 if (nargin < 2) || isempty(reports_dir) || ~isfolder(reports_dir) 36 reports_dir = []; 37 end 38 % You have to specify the folder in which the tutorial dataset is unzipped 39 if (nargin == 0) || isempty(tutorial_dir) || ~file_exist(tutorial_dir) 40 error('The first argument must be the full path to the tutorial dataset folder.'); 41 end 42 43 %% ===== FILES TO IMPORT ===== 44 % Build the path of the files to import 45 tutorial_dir = bst_fullfile(tutorial_dir, 'tutorial_seizure_fingerprinting'); 46 MriFilePre = bst_fullfile(tutorial_dir, 'anatomy', 'pre_T1.nii.gz'); 47 CtFilePost = bst_fullfile(tutorial_dir, 'anatomy', 'post_CT.nii.gz'); 48 BaselineFile = bst_fullfile(tutorial_dir, 'recordings', 'Baseline.edf'); 49 IctalFile = bst_fullfile(tutorial_dir, 'recordings', 'ictal_repetitive_spike.edf'); 50 InterictalFile = bst_fullfile(tutorial_dir, 'recordings', 'interictal_spike.edf'); 51 LvfaFile = bst_fullfile(tutorial_dir, 'recordings', 'LVFA_and_wave.edf'); 52 ElecPosFile = bst_fullfile(tutorial_dir, 'recordings', 'Subject01_electrodes_mm.tsv'); 53 MriCat12Path = bst_fullfile(tutorial_dir, 'anatomy', 'cat12'); 54 CtProcessFile = bst_fullfile(tutorial_dir, 'anatomy', 'subjectimage_post_CT_spm_reslice_masked_spm.mat'); 55 % Check if the folder contains the required files 56 if ~file_exist(BaselineFile) 57 error(['The folder ' tutorial_dir ' does not contain the folder from the file tutorial_seizure_fingerprinting.zip.']); 58 end 59 % Check for already CAT12 segmented MRI 60 isMriSegmented = file_exist(bst_fullfile(MriCat12Path, 'Subject01.nii')); 61 % Check for coregistered resliced skull-stripped CT volume 62 isCtProcessed = file_exist(CtProcessFile); 63 % Subject name 64 SubjectName = 'Subject01'; 65 66 %% ===== CREATE PROTOCOL ===== 67 % The protocol name has to be a valid folder name (no spaces, no weird characters...) 68 ProtocolName = 'TutorialSeizureFingerprinting'; 69 % Start brainstorm without the GUI 70 if ~brainstorm('status') 71 brainstorm nogui 72 end 73 % Delete existing protocol 74 gui_brainstorm('DeleteProtocol', ProtocolName); 75 % Create new protocol 76 gui_brainstorm('CreateProtocol', ProtocolName, 0, 0); 77 % Start a new report 78 bst_report('Start'); 79 % Reset visualization filters 80 panel_filter('SetFilters', 0, [], 0, [], 0, [], 0, 0); 81 % Reset colormaps 82 bst_colormaps('RestoreDefaults', 'timefreq'); 83 bst_colormaps('RestoreDefaults', 'source'); 84 % Set the current time series display mode to 'column' 85 bst_set('TSDisplayMode', 'column'); 86 % Hide scouts 87 panel_scout('SetScoutShowSelection', 'none'); 88 89 %% ===== IMPORT MRI VOLUME ===== 90 % Import and segment MRI with CAT12 91 if ~isMriSegmented 92 % Process: Import MRI 93 bst_process('CallProcess', 'process_import_mri', [], [], ... 94 'subjectname', SubjectName, ... 95 'voltype', 'mri', ... % MRI 96 'comment', 'pre_T1', ... 97 'mrifile', {MriFilePre, 'ALL'}, ... 98 'nas', [104, 207, 85], ... 99 'lpa', [ 26, 113, 78], ... 100 'rpa', [176, 113, 78]); 101 % Process: Segment MRI with CAT12 102 bst_process('CallProcess', 'process_segment_cat12', [], [], ... 103 'subjectname', SubjectName, ... 104 'nvertices', 15000, ... 105 'tpmnii', {'', 'Nifti1'}, ... 106 'sphreg', 1, ... % Use spherical registration 107 'vol', 0, ... % No volume parcellations 108 'extramaps', 0, ... % No additional cortical maps 109 'cerebellum', 0); 110 111 % Or, import already CAT12 segmented MRI 112 else 113 % Process: Import anatomy folder 114 bst_process('CallProcess', 'process_import_anatomy', [], [], ... 115 'subjectname', SubjectName, ... 116 'mrifile', {MriCat12Path, 'CAT12'}, ... 117 'nvertices', 15000, ... 118 'nas', [104, 207, 85], ... 119 'lpa', [ 26, 113, 78], ... 120 'rpa', [176, 113, 78]); 121 end 122 123 %% ===== IMPORT MRI VOLUME ===== 124 % Import, coregister, reslice and skull-strip the raw CT volume 125 if ~isCtProcessed 126 % Process: Import CT 127 bst_process('CallProcess', 'process_import_mri', [], [], ... 128 'subjectname', SubjectName, ... 129 'voltype', 'ct', ... % CT 130 'comment', 'post_CT', ... 131 'mrifile', {CtFilePost, 'ALL'}); 132 % Get updated Subject 133 sSubject = bst_get('Subject', SubjectName); 134 % Reference MRI 135 DbMriFilePre = sSubject.Anatomy(sSubject.iAnatomy).FileName; 136 % Imported CT (last volume) 137 DbCtFilePost = sSubject.Anatomy(end).FileName; 138 % Register and reslice CT to reference MRI using 'SPM' 139 DbCtFilePostRegReslice = mri_coregister(DbCtFilePost, DbMriFilePre, 'spm', 1); 140 % Skull strip the CT volume using 'SPM' 141 DbCtFilePostSkullStrip = mri_skullstrip(DbCtFilePostRegReslice, DbMriFilePre, 'spm'); 142 143 % Or import coregistered resliced skull-stripped CT volume 144 else 145 % Process: Import CT processed 146 bst_process('CallProcess', 'process_import_mri', [], [], ... 147 'subjectname', SubjectName, ... 148 'voltype', 'ct', ... % CT 149 'comment', 'post_CT_spm_reslice_masked_spm', ... 150 'mrifile', {CtProcessFile, 'BST'}); 151 % Get updated Subject 152 sSubject = bst_get('Subject', SubjectName); 153 % Imported CT processed (last volume) 154 DbCtFilePostSkullStrip = sSubject.Anatomy(end).FileName; 155 end 156 157 %% ===== CREATE SEEG CONTACT IMPLANTATION ===== 158 iStudyImplantation = db_add_condition(SubjectName, 'Implantation'); 159 % Import locations and convert to subject coordinate system (SCS) 160 ImplantationChannelFile = import_channel(iStudyImplantation, ElecPosFile, 'BIDS-SCANRAS-MM', 1, 0, 1, 0, 2, DbCtFilePostSkullStrip); 161 % Snapshot: SEEG electrodes in MRI slices 162 hFigMri3d = view_channels_3d(ImplantationChannelFile, 'SEEG', 'anatomy', 1, 0); 163 bst_report('Snapshot', hFigMri3d, ImplantationChannelFile, 'SEEG electrodes in 3D MRI slices'); 164 close(hFigMri3d); 165 166 %% ===== ACCESS THE RECORDINGS ===== 167 % Process: Create link to raw file 168 sFilesRaw = bst_process('CallProcess', 'process_import_data_raw', [], [], ... 169 'subjectname', SubjectName, ... 170 'datafile', {{BaselineFile, LvfaFile, IctalFile, InterictalFile}, 'EEG-EDF'}, ... 171 'channelreplace', 0, ... 172 'channelalign', 0); 173 % Process: Add EEG positions 174 bst_process('CallProcess', 'process_channel_addloc', sFilesRaw, [], ... 175 'channelfile', {ImplantationChannelFile, 'BST'}, ... 176 'fixunits', 0, ... % No automatic fixing of distance units required 177 'vox2ras', 0); % Do not use the voxel=>subject transformation, already in SCS 178 179 %% ===== REVIEW RECORDINGS ===== 180 % Process: Power spectrum density (Welch) 181 sFilesPsd = bst_process('CallProcess', 'process_psd', sFilesRaw, [], ... 182 'timewindow', [], ... 183 'win_length', 5, ... 184 'win_overlap', 50, ... 185 'sensortypes', 'SEEG', ... 186 'edit', struct(... 187 'Comment', 'Power', ... 188 'TimeBands', [], ... 189 'Freqs', [], ... 190 'ClusterFuncTime', 'none', ... 191 'Measure', 'magnitude', ... 192 'Output', 'all', ... 193 'SaveKernel', 0)); 194 195 % Process: Snapshot: PSD with power line noise 196 panel_display('SetDisplayFunction', 'log'); 197 bst_process('CallProcess', 'process_snapshot', sFilesPsd, [], ... 198 'target', 10, ... % Frequency spectrum 199 'modality', 6, ... % SEEG 200 'Comment', 'Power spectrum density'); 201 202 % Process: Set channels type 203 % 'MPS16' channel needs to be excluded because for BEM head modeling it lies outside the inner skull 204 bst_process('CallProcess', 'process_channel_settype', sFilesRaw, [], ... 205 'sensortypes', 'MPS16', ... 206 'newtype', 'SEEG_NO_LOC'); 207 % Define event: LVFA & wave and ictal repetitive spike 208 sEvt1 = db_template('event'); 209 sEvt1.label = 'sEEG onset'; 210 sEvt1.epochs = 1; 211 sEvt1.times = 15; 212 % Define event: Interictal spike 213 sEvt2 = db_template('event'); 214 sEvt2.label = 'Interictal spike'; 215 sEvt2.epochs = 1; 216 sEvt2.times = 5; 217 % Process: Events: Import from file 218 bst_process('CallProcess', 'process_evt_import', sFilesRaw(2:3), [], ... 219 'evtfile', {sEvt1, 'struct'}, ... 220 'evtname', ''); 221 bst_process('CallProcess', 'process_evt_import', sFilesRaw(4), [], ... 222 'evtfile', {sEvt2, 'struct'}, ... 223 'evtname', ''); 224 225 %% ===== IMPORT RECORDINGS ===== 226 % Process: Import SEEG event for LVFA & wave and ictal repetitive spike to database 227 sFilesOnset = bst_process('CallProcess', 'process_import_data_event', sFilesRaw(2:3), [], ... 228 'subjectname', SubjectName, ... 229 'eventname', 'sEEG onset', ... 230 'epochtime', [-15, 15], ... 231 'createcond', 0, ... 232 'ignoreshort', 0, ... 233 'usessp', 0, ... 234 'baseline', 'all', ... % Remove DC offset: All recordings 235 'blsensortypes', 'SEEG'); % Sensor types to remove DC offset 236 % Process: Import SEEG event for interictal spike to database 237 sFileInterictalSpike = bst_process('CallProcess', 'process_import_data_event', sFilesRaw(4), [], ... 238 'subjectname', SubjectName, ... 239 'eventname', 'Interictal spike', ... 240 'epochtime', [-5, 5], ... 241 'createcond', 0, ... 242 'ignoreshort', 0, ... 243 'usessp', 0, ... 244 'baseline', 'all', ... % Remove DC offset: All recordings 245 'blsensortypes', 'SEEG'); % Sensor types to remove DC offset 246 % ===== Bipolar Montage ===== 247 MontageSeegBipName = [SubjectName, ': SEEG (bipolar 2)[tmp]']; 248 % Apply montage (create new folders) 249 sFilesOnsetBip = bst_process('CallProcess', 'process_montage_apply', sFilesOnset, [], ... 250 'montage', MontageSeegBipName, ... 251 'createchan', 1); 252 253 %% ===== HEAD MODELING ===== 254 % Process: Generate BEM surfaces 255 bst_process('CallProcess', 'process_generate_bem', [], [], ... 256 'subjectname', SubjectName, ... 257 'nscalp', 1922, ... 258 'nouter', 1922, ... 259 'ninner', 1922, ... 260 'thickness', 4, ... 261 'method', 'brainstorm'); 262 263 % Snapshot: BEM surfaces 264 sSubject = bst_get('Subject', SubjectName); 265 BemInnerSkullFile = sSubject.Surface(sSubject.iInnerSkull).FileName; 266 BemOuterSkullFile = sSubject.Surface(sSubject.iOuterSkull).FileName; 267 BemScalpFile = sSubject.Surface(sSubject.iScalp).FileName; 268 hFigSurf = view_surface(BemInnerSkullFile); 269 hFigSurf = view_surface(BemOuterSkullFile, [], [], hFigSurf); 270 hFigSurf = view_surface(BemScalpFile, [], [], hFigSurf); 271 figure_3d('SetStandardView', hFigSurf, 'left'); % Set orientation (left) 272 bst_report('Snapshot', hFigSurf, BemInnerSkullFile, 'BEM surfaces'); 273 close(hFigSurf); 274 275 % Process: Compute head model 276 bst_process('CallProcess', 'process_headmodel', sFileInterictalSpike, [], ... 277 'comment', '', ... 278 'sourcespace', 1, ... % Cortex surface 279 'meg', 1, ... % None 280 'eeg', 1, ... % None 281 'ecog', 1, ... % None 282 'seeg', 2, ... % OpenMEEG BEM 283 'openmeeg', struct(... 284 'BemSelect', [0, 0, 1], ... % Only compute on BEM inner skull 285 'BemCond', [1, 0.0125, 1], ... 286 'BemNames', {{'Scalp', 'Skull', 'Brain'}}, ... 287 'BemFiles', {{}}, ... 288 'isAdjoint', 0, ... 289 'isAdaptative', 1, ... 290 'isSplit', 0, ... 291 'SplitLength', 4000)); 292 % Copy head model to other folders 293 sHeadModel = bst_get('HeadModelForStudy', sFileInterictalSpike.iStudy); 294 db_set_headmodel(sHeadModel.FileName, 'AllConditions'); 295 % Process: Compute noise covariance in Baseline 296 bst_process('CallProcess', 'process_noisecov', sFilesRaw(1), [], ... 297 'baseline', [0, 300.9995], ... 298 'dcoffset', 1, ... % Block by block, to avoid effects of slow shifts in data 299 'identity', 0, ... 300 'copycond', 1, ... % Copy to other folders 301 'copysubj', 0); 302 303 % Process: Snapshot: Noise covariance 304 bst_process('CallProcess', 'process_snapshot', sFilesRaw(1), [], ... 305 'target', 3, ... % Noise covariance 306 'Comment', 'Noise covariance'); 307 308 %% ===== MODELING INTERICTAL SPIKES ===== 309 % Process: Compute sources [2018] (SEEG) 310 sFileInterictalSpikeSrc = bst_process('CallProcess', 'process_inverse_2018', sFileInterictalSpike, [], ... 311 'output', 1, ... % Kernel only: shared 312 'inverse', struct(... 313 'Comment', '', ... 314 'InverseMethod', 'minnorm', ... 315 'InverseMeasure', 'sloreta', ... 316 'SourceOrient', {{'fixed'}}, ... 317 'UseDepth', 0, ... 318 'NoiseMethod', 'diag', ... 319 'SnrMethod', 'fixed', ... 320 'SnrRms', 1e-06, ... 321 'SnrFixed', 3, ... 322 'ComputeKernel', 1, ... 323 'DataTypes', {{'SEEG'}})); 324 325 % Interictal: Snapshots 326 Time = 0.041; % First peak of SPS10-SPS11 at 41ms 327 TimeWindow = [-0.5 0.5]; % Time window: -500ms to 500ms 328 DataThresh = 0.26; % Source threshold (percentage) 329 GetSnapshotSensorTimeSeries(sFileInterictalSpike.FileName, [SubjectName, ': SPS (bipolar 2)[tmp]'], Time, TimeWindow); 330 GetSnapshotSensor2DLayout(sFileInterictalSpike.FileName, Time, TimeWindow); 331 GetSnapshotsSources(sFileInterictalSpikeSrc.FileName, 'srf3d', Time, DataThresh); 332 GetSnapshotsSources(sFileInterictalSpikeSrc.FileName, 'mri3d', Time, DataThresh); 333 GetSnapshotsSources(sFileInterictalSpikeSrc.FileName, 'mri2d', Time, DataThresh); 334 335 % ===== Create a Desikan-Killiany atlas with scouts only in the right hemisphere ==== 336 % Load the surface 337 [hFigSurf, iDS, iFig] = view_surface_data([], sFileInterictalSpikeSrc.FileName); 338 % Show the SEEG electrodes 339 figure_3d('PlotSensors3D', iDS, iFig); 340 % Set Desikan-Killiany as the current atlas 341 [~, ~, sSurf] = panel_scout('GetAtlas'); 342 iAtlas = find(strcmpi('Desikan-Killiany', {sSurf.Atlas.Name})); 343 panel_scout('SetCurrentAtlas', iAtlas); 344 % Set scout options to display all the scouts 345 panel_scout('SetScoutsOptions', 0, 0, 1, 'all', 0.7, 1, 0, 0); 346 % Select the scouts in the right hemisphere 347 iScoutsR = find(cellfun(@(s) ~isempty(regexp(s, 'R$', 'once')), {sSurf.Atlas(iAtlas).Scouts.Label})); 348 panel_scout('SetSelectedScouts', iScoutsR); 349 % Create a new atlas from selected scouts 350 panel_scout('CreateAtlasSelected', 0, 0); 351 % Rename the atlas 352 [sAtlas, iAtlas] = panel_scout('GetAtlas'); 353 sAtlas.Name = 'Desikan-Killiany_RH'; 354 panel_scout('SetAtlas', [], iAtlas, sAtlas); 355 % Subdivide the new atlas to get scouts with 5 cm sq. area each 356 panel_scout('SubdivideScouts', 1, 'area', 5); 357 % Get the scouts 358 sScouts = panel_scout('GetScouts'); 359 % Close the figure 360 close(hFigSurf); 361 362 %% ===== MODELING ICTAL WAVE ===== 363 % Process: Compute sources [2018] (SEEG) 364 sFileLvfaOnsetSrc = bst_process('CallProcess', 'process_inverse_2018', sFilesOnset(1), [], ... 365 'output', 1, ... % Kernel only: shared 366 'inverse', struct(... 367 'Comment', '', ... 368 'InverseMethod', 'minnorm', ... 369 'InverseMeasure', 'sloreta', ... 370 'SourceOrient', {{'fixed'}}, ... 371 'UseDepth', 0, ... 372 'NoiseMethod', 'diag', ... 373 'SnrMethod', 'fixed', ... 374 'SnrRms', 1e-06, ... 375 'SnrFixed', 3, ... 376 'ComputeKernel', 1, ... 377 'DataTypes', {{'SEEG'}})); 378 379 % Snapshots: Sensor and source time series 380 Time = 0.270; % % Wave activity at 270ms 381 TimeWindow = [-0.5 0.5]; % Time window: -500ms to 500ms 382 DataThresh = 0.45; % Source threshold (percentage) 383 GetSnapshotSensorTimeSeries(sFilesOnset(1).FileName, [SubjectName, ': SPS (bipolar 2)[tmp]'], Time, TimeWindow); 384 GetSnapshotSensor2DLayout(sFilesOnset(1).FileName, Time, TimeWindow); 385 GetSnapshotsSources(sFileLvfaOnsetSrc.FileName, 'srf3d', Time, DataThresh); 386 GetSnapshotsSources(sFileLvfaOnsetSrc.FileName, 'mri3d', Time, DataThresh); 387 GetSnapshotsSources(sFileLvfaOnsetSrc.FileName, 'mri2d', Time, DataThresh); 388 389 390 %% ===== MODELING ICTAL ONSET WITH LVFA (SENSOR SPACE) ===== 391 % Process: Time-frequency (Morlet wavelets) 392 sFilesOnsetBipTf = bst_process('CallProcess', 'process_timefreq', sFilesOnsetBip(1), [], ... 393 'sensortypes', 'SEEG', ... 394 'edit', struct(... 395 'Comment', 'Power,1-100Hz', ... 396 'TimeBands', [], ... 397 'Freqs', [1, 2.1, 3.3, 4.7, 6.1, 7.8, 9.6, 11.5, 13.7, 16.1, 18.7, 21.6, 24.8, ... 398 28.3, 32.1, 36.4, 41.1, 46.2, 51.9, 58.1, 64.9, 72.5, 80.8, 89.9, 100], ... 399 'MorletFc', 1, ... 400 'MorletFwhmTc', 6, ... 401 'ClusterFuncTime', 'none', ... 402 'Measure', 'power', ... 403 'Output', 'all', ... 404 'SaveKernel', 0, ... 405 'Method', 'morlet'), ... 406 'normalize2020', 'multiply2020'); % Spectral flattening: Multiply output power values by frequency 407 408 409 % Snapshots: timefrequency map for sensor 410 Brightness = 0.65; % Brightness -65% 411 Contrast = 0.49; % Contrast 49% 412 TFpoint = [0.6, 65]; % TimeFreq point [s, Hz] 413 GetSnapshotTimeFreq(sFilesOnsetBipTf.FileName, 'AllSensors', TFpoint, 0, Brightness, Contrast); 414 GetSnapshotTimeFreq(sFilesOnsetBipTf.FileName, 'SPS8-SPS9', TFpoint, 1, Brightness, Contrast); 415 416 %% ===== MODELING ICTAL ONSET WITH LVFA (SOURCE SPACE) ===== 417 % Process: Extract scout time series 418 sFileLvfaOnsetScoutTs = bst_process('CallProcess', 'process_extract_scout', sFileLvfaOnsetSrc, [], ... 419 'timewindow', [-15, 15], ... 420 'scouts', {'Desikan-Killiany_RH', {sScouts.Label}}, ... 421 'flatten', 0, ... 422 'scoutfunc', 'pca', ... % PCA 423 'pcaedit', struct(... 424 'Method', 'pcai', ... 425 'Baseline', [NaN, NaN], ... 426 'DataTimeWindow', [-15, 15], ... 427 'RemoveDcOffset', 'none'), ... 428 'isflip', 1, ... 429 'isnorm', 0, ... 430 'concatenate', 1, ... 431 'save', 1, ... 432 'addrowcomment', 1, ... 433 'addfilecomment', []); 434 435 % Snapshot: Scout time series 436 bst_process('CallProcess', 'process_snapshot', sFileLvfaOnsetScoutTs, [], ... 437 'target', 5, ... % Data 438 'modality', 6, ... % SEEG 439 'Comment', 'Scout time series (matrix)'); 440 441 %% Process: Time-frequency (Morlet wavelets) 442 sFileLvfaOnsetTf = bst_process('CallProcess', 'process_timefreq', sFileLvfaOnsetScoutTs, [], ... 443 'sensortypes', 'SEEG', ... 444 'edit', struct(... 445 'Comment', 'Power,1-100Hz', ... 446 'TimeBands', [], ... 447 'Freqs', [1, 2.1, 3.3, 4.7, 6.1, 7.8, 9.6, 11.5, 13.7, 16.1, 18.7, 21.6, 24.8, ... 448 28.3, 32.1, 36.4, 41.1, 46.2, 51.9, 58.1, 64.9, 72.5, 80.8, 89.9, 100], ... 449 'MorletFc', 1, ... 450 'MorletFwhmTc', 6, ... 451 'ClusterFuncTime', 'none', ... 452 'Measure', 'power', ... 453 'Output', 'all', ... 454 'SaveKernel', 0, ... 455 'Method', 'morlet'), ... 456 'normalize2020', 'multiply2020'); % Spectral flattening: Multiply output power values by frequency 457 458 % Snapshots: Time-frequency map for scouts 459 Brightness = 0.65; % Brightness -65% 460 Contrast = 0.49; % Contrast 49% 461 TFpoint = [0.6, 65]; % TimeFreq point [s, Hz] 462 GetSnapshotTimeFreq(sFileLvfaOnsetTf.FileName, 'AllSensors', TFpoint, 0, Brightness, Contrast); 463 GetSnapshotTimeFreq(sFileLvfaOnsetTf.FileName, 'postcentral R.3', TFpoint, 1, Brightness, Contrast); 464 465 %% ===== MODELING ICTAL ONSET WITH REPETITIVE SPIKING (SENSOR SPACE) ===== 466 % Process: Time-frequency (Morlet wavelets) 467 sFilesOnsetTf = bst_process('CallProcess', 'process_timefreq', sFilesOnsetBip(2), [], ... 468 'sensortypes', 'PIN5-PIN6', ... 469 'edit', struct(... 470 'Comment', 'Power,1-100Hz', ... 471 'TimeBands', [], ... 472 'Freqs', [1, 2.1, 3.3, 4.7, 6.1, 7.8, 9.6, 11.5, 13.7, 16.1, 18.7, 21.6, 24.8, ... 473 28.3, 32.1, 36.4, 41.1, 46.2, 51.9, 58.1, 64.9, 72.5, 80.8, 89.9, 100], ... 474 'MorletFc', 1, ... 475 'MorletFwhmTc', 6, ... 476 'ClusterFuncTime', 'none', ... 477 'Measure', 'power', ... 478 'Output', 'all', ... 479 'SaveKernel', 0, ... 480 'Method', 'morlet'), ... 481 'normalize2020', 'multiply2020'); % Spectral flattening: Multiply output power values by frequency 482 483 % Snapshot: Sensor time series (PIN bipolar montage) 484 Time = 7.7325; 485 GetSnapshotSensorTimeSeries(sFilesOnsetBip(2).FileName, [SubjectName ': PIN (orig)[tmp]'], Time); 486 487 % Snapshot: Time-frequency maps (one sensor) 488 Brightness = 0.60; % Brightness -60% 489 Contrast = 0.23; % Contrast 23% 490 TFpoint = [Time, 25]; % TimeFreq point [s, Hz] 491 GetSnapshotTimeFreq(sFilesOnsetTf.FileName, 'PIN5-PIN6', TFpoint, 1, Brightness, Contrast); 492 493 %% ===== MODELING ICTAL ONSET WITH REPETITIVE SPIKING (SOURCE SPACE) ===== 494 % Process: Compute sources [2018] (SEEG) 495 sFilesOnsetSrc = bst_process('CallProcess', 'process_inverse_2018', sFilesOnset(2), [], ... 496 'output', 1, ... % Kernel only: shared 497 'inverse', struct(... 498 'Comment', '', ... 499 'InverseMethod', 'minnorm', ... 500 'InverseMeasure', 'sloreta', ... 501 'SourceOrient', {{'fixed'}}, ... 502 'UseDepth', 0, ... 503 'NoiseMethod', 'diag', ... 504 'SnrMethod', 'fixed', ... 505 'SnrRms', 1e-06, ... 506 'SnrFixed', 3, ... 507 'ComputeKernel', 1, ... 508 'DataTypes', {{'SEEG'}})); 509 510 % Set freq filters 511 panel_filter('SetFilters', 1, 55, 1, 5, 0, [], 0, 0); 512 % Set colormap for sources 513 sColormapSrc = bst_colormaps('GetColormap', 'source'); 514 sColormapSrc.MaxMode = 'custom'; 515 sColormapSrc.MinValue = 0; 516 sColormapSrc.MaxValue = 2e-8; 517 bst_colormaps('SetColormap', 'source', sColormapSrc); % Save the changes in colormap 518 % Snapshot: Sensor time series (PIN) 519 Time = 9.719; 520 TimeWindow = [-0.5, 0.5] + Time; 521 DataThresh = 0.33; 522 GetSnapshotSensorTimeSeries(sFilesOnset(2).FileName, [SubjectName ': PIN (orig)[tmp]'], Time, TimeWindow); 523 % Snapshot: Sources (display on MRI Viewer) 524 GetSnapshotsSources(sFilesOnsetSrc.FileName, 'mri2d', Time, DataThresh); 525 % Reset freq filters and colormap for sources 526 panel_filter('SetFilters', 0, [], 0, [], 0, [], 0, 0); 527 bst_colormaps('RestoreDefaults', 'source'); 528 529 %% ===== SAVE AND DISPLAY REPORT ===== 530 ReportFile = bst_report('Save', []); 531 if ~isempty(reports_dir) && ~isempty(ReportFile) 532 bst_report('Export', ReportFile, reports_dir); 533 else 534 bst_report('Open', ReportFile); 535 end 536 537 disp([10 'DEMO> Seizure Fingerpriting tutorial completed' 10]); 538 539 % =================================================================% 540 % ===================== SNAPSHOTS FUNCTIONS =======================% 541 % =================================================================% 542 %% ===== SNAPSHOTS: SENSOR TIME SERIES ===== 543 function GetSnapshotSensorTimeSeries(SensorFile, MontageName, Time, TimeWindow) 544 if nargin < 4 545 TimeWindow = []; 546 end 547 % Figure: Sensor time series (set montage) 548 hFig = view_timeseries(SensorFile, 'SEEG'); 549 panel_montage('SetCurrentMontage', hFig, MontageName); 550 if ~isempty(TimeWindow) 551 h1 = findobj(hFig, 'Tag','AxesGraph','-or','Tag','AxesEventsBar'); 552 xlim(h1, TimeWindow); 553 end 554 panel_time('SetCurrentTime', Time); 555 bst_report('Snapshot', hFig, SensorFile, 'Sensor time series (SPS bipolar)', [200, 200, 400, 400]); 556 close(hFig); 557 end 558 559 %% ===== SNAPSHOTS: SENSOR 2D LAYOUT TIME SERIES ===== 560 function GetSnapshotSensor2DLayout(SensorFile, Time, TimeWindow) 561 % Snapshot: 2D layout sensor time series 562 hFig = view_topography(SensorFile, 'SEEG', '2DLayout'); 563 figure_topo('SetTopoLayoutOptions', 'TimeWindow', TimeWindow); 564 panel_time('SetCurrentTime', Time); 565 bst_report('Snapshot', hFig, SensorFile, '2D layout sensor time series', [200, 200, 400, 400]); 566 close(hFig); 567 end 568 569 %% ===== SNAPSHOTS: SOURCES TIME SLICE ===== 570 function GetSnapshotsSources(SourceFile, FigType, Time, DataThreshold) 571 [sStudy, iStudy] = bst_get('AnyFile', SourceFile); 572 % Get MRI reference 573 if ~isempty(regexp(FigType, '^mri', 'once')) 574 sSubjectFig = bst_get('Subject', sStudy.BrainStormSubject); 575 MriFile = sSubjectFig.Anatomy(sSubjectFig.iAnatomy).FileName; 576 end 577 % Get ChannelFile 578 if ~isempty(regexp(FigType, '3d$', 'once')) 579 ChanFile = bst_get('ChannelForStudy', iStudy); 580 end 581 582 switch FigType 583 % Display sources on its default surface (cortex) 584 case 'srf3d' 585 hFig = view_surface_data([], SourceFile); 586 panel_time('SetCurrentTime', Time); 587 hFig = view_channels(ChanFile.FileName, 'SEEG', 0, 0, hFig, 1); 588 panel_surface('SetDataThreshold', hFig, 1, DataThreshold); 589 displayStr = 'cortex'; 590 591 % Display sources on display on 3D MRI 592 case 'mri3d' 593 hFig = view_surface_data(MriFile, SourceFile); 594 panel_time('SetCurrentTime', Time); 595 hFig = view_channels(ChanFile.FileName, 'SEEG', 0, 0, hFig, 1); 596 figure_3d('JumpMaximum', hFig); 597 figure_3d('SetStandardView', hFig, 'right'); 598 displayStr = '3D MRI'; 599 600 % Display sources on MRI Viewer 601 case 'mri2d' 602 hFig = view_mri(MriFile, SourceFile); 603 panel_time('SetCurrentTime', Time); 604 bst_figures('GetFigureHandles', hFig); 605 figure_mri('JumpMaximum', hFig); 606 displayStr = 'MRI viewer'; 607 end 608 panel_surface('SetDataThreshold', hFig, 1, DataThreshold); 609 bst_report('Snapshot', hFig, SourceFile, ['Sources: Display on ' displayStr], [200, 200, 400, 400]); 610 close(hFig); 611 end 612 613 %% ===== SNAPSHOTS: TIME-FREQUENCY ===== 614 function GetSnapshotTimeFreq(TimefreqFile, RowName, TimeFreqPoint, doSlices, Brightness, Contrast) 615 WinPos = [200, 200, 600, 400]; 616 if nargin < 4 || isempty(doSlices) 617 doSlices = 0; 618 end 619 % TimeFreq display mode 620 DisplayMode = 'SingleSensor'; 621 if strcmpi(RowName, 'AllSensors') 622 DisplayMode = RowName; 623 end 624 % Snapshot: Time frequency maps (all sensors/sources) 625 hFigTF = view_timefreq(TimefreqFile, DisplayMode); 626 sOptions = panel_display('GetDisplayOptions'); 627 sOptions.Function = 'log'; % Log power 628 sOptions.HighResolution = 1; % Smooth display 629 if ~strcmpi(RowName, 'AllSensors') 630 sTimefreq = in_bst_timefreq(TimefreqFile, 0, 'RowNames'); 631 iRow = ~cellfun(@isempty, regexp(sTimefreq.RowNames, ['^' RowName])); 632 sOptions.RowName = sTimefreq.RowNames{iRow}; 633 end 634 panel_display('SetDisplayOptions', sOptions); 635 bst_colormaps('SetColormapAbsolute', 'timefreq', 0); % Turn off absolute value 636 sColormap = bst_colormaps('GetColormap', hFigTF); 637 sColormap.Contrast = Contrast; 638 sColormap.Brightness = Brightness; 639 % Apply modifiers (for brightness and contrast) 640 sColormap = bst_colormaps('ApplyColormapModifiers', sColormap); 641 % Save the changes in colormap 642 bst_colormaps('SetColormap', 'timefreq', sColormap); 643 % Update the colormap in figures 644 bst_colormaps('FireColormapChanged', 'timefreq'); 645 % Select Time-Frequency point 646 if length(TimeFreqPoint) == 2 647 panel_time('SetCurrentTime', TimeFreqPoint(1)); 648 panel_freq('SetCurrentFreq', TimeFreqPoint(2), 0); 649 end 650 bst_report('Snapshot', hFigTF, TimefreqFile, ['Time frequency map (' RowName ')'], WinPos); 651 % Power time series and power spectrum from TF representation 652 if doSlices && length(TimeFreqPoint) == 2 && ~strcmpi(RowName, 'AllSensors') 653 hFigT = view_spectrum(TimefreqFile, 'TimeSeries', sOptions.RowName); 654 hFigF = view_spectrum(TimefreqFile, 'Spectrum', sOptions.RowName); 655 bst_report('Snapshot', hFigF, TimefreqFile, ['Time series (' RowName ')'], WinPos); 656 bst_report('Snapshot', hFigT, TimefreqFile, ['Power spectrum (' RowName ')'], WinPos); 657 close([hFigT hFigF]); 658 end 659 close(hFigTF); 660 end 661 end








Feedback: Comments, bug reports, suggestions, questions
Email address (if you expect an answer):


seeg/SeizureFingerprinting (last edited 2025-10-01 17:21:41 by ChinmayChinara)