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