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

TODO: Add flowchart

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

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








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


seeg/SeizureFingerprinting (last edited 2025-07-24 18:10:17 by ChinmayChinara)