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

This tutorial introduces concepts specific to managing intracranial stereoencephalography (SEEG recordings) within the Brainstorm environment. It guides users through computing time-frequency decomposition maps to identify the epileptogenic zone (EZ) using both ictal and interictal SEEG data. The tutorial is based on a clinical case from the McGovern Medical School at the University of Texas Health Science Center at Houston, Texas, USA.

Please note that this tutorial is intended for users alerady familiar with Brainstorm. It does nto provide detailed explanations of the software's interface or theoretical foundations. For comprehensive introductory material, refer to the Brainstorm introduction tutorials.

NOT FOR CLINICAL USE:
The methods and software implementations presented in this tutorial have not been certified as medical devices. They are intended for research purposes only and should not be used for clinical decision-making.

Dataset description

License

This EEG, MRI, and CT data provided in this tutorial remain the property of the McGovern Medical School at the University of Texas Health Science Center at Houston. Use or distribution of this dataset outside the scope of the Brainstorm tutorials - including for research purposes - is striclty prohibited without prior written consent. For inquiries regarding permissions, please use the Brainstorm user forum.

Clinical description

The dataset featured in this tutorial was recorded at the Epilepsy Monitoring Unit (EMU) of UTHealth Houston. It pertains to a 25-year-old right-handed woman with drug-resistant epilepsy since the age of six. At 15, she underwent a right parietal opercular corticectomy. Despite this intervention, she continued to experience weekly focal aware seizures characterized by a left-hand tingling aura, as well as focal impaired awareness seizures presenting with staring and pouting.

During her EMU admission, intermittent right parietal slowing was observed, and ten habitual seizures were recorded, originating from the C4–P4 region. Magnetic resonance imaging (MRI) revealed bilateral perisylvian polymicrogyria (PMG), pachygyria, right posterior temporal periventricular nodular heterotopia, and post-surgical changes. Magnetoencephalography (MEG) localized discharges to the right superior parietal region adjacent to her previous resection.

SEEG implantation identified two distinct seizure onset patterns: 1) Low-voltage fast activity in the right superior parietal PMG during focal aware seizures; 2) Repetitive spiking in the posterior insular PMG during impaired awareness seizures.

Following a multidisciplinary evaluation, the patient underwent MR-guided laser interstitial thermal therapy (LITT) targeting the right superior parietal and posterior insular PMG. The procedure was uncomplicated, and she remained seizure-free at a one-year follow-up.

SEEG recordings

https://neuroimage.usc.edu/brainstorm/Tutorials/SeizureFingerprinting?action=AttachFile&do=get&target=pmt.png

The recording electrodes used in this dataset are PMT SEEG Depth Electrodes, with the following specifications:

Files

tutorial_seizure_fingerprinting/

References

Further details for this study can be found under the following link: https://zenodo.org/records/14807262

Download and installation

Import the anatomy

Pre-implantation MRI

While this is not applicable to the current dataset, some NIfTI MRI files may include a transformation matrix in their header. In such cases, a pop-up window may appear asking: Do you want to apply the transformation to the MRI file? Selecting Yes will apply the transformation and reorient the MRI into Brainstorm's standard orientation, allowing you to view the coronal, sagittal, and axial planes correctly.

Generate default surfaces using CAT12

We recommend generating cortical surfaces with CAT12, especially if you are interested in a realistic representation of the patient's cortical folding in 3D. Follow the CAT12 tutorial to generate the surfaces as under.

14_seg_cat12.png

These surfaces will be used later, in the computation of the epileptogenicity maps. Read the section Importing realistic surfaces for information on how to use realistic surfaces from BrainVISA or FreeSurfer.

Segmentation using CAT12 can take around 1 hour depending on your system. To save time, we provide the precomputed CAT12 segmented surfaces generated using the MRI above as part of the tutorial dataset (tutorial_seizure_fingerprinting/cat12. To import them just right click on the Subject01 > Import anatomy folder (auto) and select the precomputed CAT12 folder above. More details can be found in the CAT12 tutorial.

Post-implantation CT

The pre-implantation MRI imported above will serve as the anatomical reference for this subject.

We will now import a second scan acquired after SEEG electrode implantation, in which the electrode contacts are visible. In this dataset, the post-implantation volume is a CT scan, where the SEEG contacts appear as bright (hypersignal) spots.

You can also perform skull stripping using BrainSuite's Brain Surface Extractor. Installation steps are detailed in the BrainSuite for Brainstorm tutorial.

Electrode labeling and contact localization

If you do not have any recordings in the database, Brainstorm allows to creation and annotation of intracranial electrodes and contacts. Users can also then export these as a text file with all the positions that can be used in Brainstorm or any other program.

Generate isosurface

This creates a thresholded mesh from the CT by separating the contacts out from rest of the CT. This aids the user towards localization of the electrodes and its contacts more accurately.

Start implantation

MRI: MRI viewer loads up with MRI volume only.
CT: MRI viewer loads up with CT volume only.
MRI+CT: MRI viewer loads up with the CT overlayed on the MRI.
MRI+CT+IsoSurf: MRI viewer loads up with the CT overlayed on the MRI. 3D figure loads up with the isosurface and 3D MRI slices.

Create electrodes and plot contacts manually

Before we start the implantation a prior knowledge of the implantation scheme is required in order to have the correct labels of the various electrodes used. One way here is to have a look at the recordings file and get a knowledge of that. Brainstorm matches the channel names to that of the recordings while importing the positions to them.

In some cases, additional correction of the contacts may be required. To edit the individual contacts refer to the Edit the contacts positions advanced section.

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

Advanced

Edit the contacts positions

The trajectory of electrode while implantation may not always follow a straight line as there could be bending introduced when the neurosurgeon inserts the electrode. In such cases we need to move these contacts to more appropriate positions.

Advanced

Export the contacts position

You can export the contacts created in Brainstorm as a text file to be used later in Brainstorm or in an external software.

Additional Documentation

Forum discussions

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

Tutorials/SeizureFingerprinting (last edited 2025-06-18 15:43:59 by ChinmayChinara)