Atypical microstates topographic maps - Pycrostates

Hi,

I’m having issues obtaining topographic maps that look anything like what’s normally seen in the literature. See images. Basically, the maps are inconsistent and not reflecting what should be expected at the individual level, or the group level.

I suspect that my preprocessing may be wrong or incomplete, but it may be something else.

The files are individual subjects from Matlab. I’ve done preprocessing in a similar manner as Artoni & Michel (2025) How does Independent Component Analysis Preprocessing Affect EEG Microstates? with an additional 60-120Hz notch filter (see code below) and a different ICA method. I’ve included a screenshot of the psd and the full time scale of the channels in the post (please disregard channel 64 named ref). The full time scale of EEG looks a little rough but because the time period is large, it’s hard to tell. There are no artefacts, some drift at the end to the time segment. The psd looks fine and the filters seem to be doing their job.

I’m getting plenty of GFP peaks and the segmentation seems to be working. Until i get satisfying topographic maps, i’m not sure i should interpret any statistical test or even the quality of the segmentation.

I’ve tried many different n_cluster settings in ModKMeans(), different min_peak_distance settings in extract_gfp_peaks(). I’ve tried resampling from GFP or epochs but nothing changed.

See images and preprocessing code below.

GFP peaks:
Out[152]: < ChData | 12480 samples >

Thanks in advance for any help.

EEG:

  • Brain Products ActiCap/ActiChamp
  • montage: 64 Channels (+3 EOG); easycap-M43
  • Length of recording: 310secs (size77501)

Software:

  • Pycrostates 0.6.1
  • Spyder version: 5.4.5 (conda)
  • Python version: 3.11.5 64-bit
  • Qt version: 5.15.8
  • PyQt5 version: 5.15.9
  • Operating System: macOS-14.6.1-x86_64-i386-64bit
# Preprocessing step largely inspired by Artoni & Michel (2025) *How does Independent Component Analysis Preprocessing Affect EEG Microstates?*
import os
import numpy as np
import sys
import subprocess
import matplotlib
from pycrostates.preprocessing import extract_gfp_peaks, apply_spatial_filter, resample
from pycrostates.cluster import ModKMeans
from scipy.signal import cheby2, filtfilt
import scipy.io as sio
from scipy.stats import kurtosis
import mne
from mne_icalabel import label_components
from mne.channels import find_ch_adjacency
from mne.viz import plot_ch_adjacency, plot_topomap


# Verbose mne
mne.set_log_level('info')

# --- Load Groups from a Matlab File ---
x_file = "~/[...].mat" # <--- --- --- --- --- --- --- --- --- --- 
x_data = sio.loadmat(x_file)  
x_scores = x_data["x"].flatten()  
HS = np.where(x_scores > 8)[0] + 1  
LS = np.where(x_scores < 4)[0] + 1
MS = np.setdiff1d(np.arange(1, len(x_scores) + 1), np.concatenate((HS, LS)))

# --- Load Standard EasyCap Montage ---
montage = mne.channels.make_standard_montage("easycap-M43")

data_path = "/[...]folder_x" # <--- --- --- --- --- --- --- --- --- --- 
sfreq = 250  

# Storage structures
subject_data = {}  # Store individual subject data
group_data = {"HS": [], "LS": [], "MS": []}  # Store group-level data
all_data = {}  # Store everything

save_dir = "/[...]folder_x" # <--- --- --- --- --- --- --- --- --- --- 
os.makedirs(save_dir, exist_ok=True)

for idx in range(len(x_scores)):
    subject_id = idx + 1  
    file_name = os.path.join(data_path, f"S_{subject_id}_epoched_data.mat") # <--- --- --- --- --- --- --- --- --- --- 
    
    if os.path.exists(file_name):
        mat_data = sio.loadmat(file_name)  
        
        if "F" in mat_data:
            eeg_data = mat_data["F"]  
            
            num_channels = eeg_data.shape[0]
            
            # Channel ID
            raw_channel_names = montage.ch_names[:64] + ["EOG1", "EOG2", "EOG3"]
            ch_types = ["eeg"] * 64 + ["eog"] * 3
            info = mne.create_info(ch_names=raw_channel_names, sfreq=sfreq, ch_types=ch_types)
            raw = mne.io.RawArray(eeg_data, info)
            raw.set_montage(montage, on_missing="ignore")

            # --- Preprocessing ---
            # Chebyshev type II filters
            raw.filter(l_freq=1, h_freq=None, method='iir', iir_params=dict(order=24, ftype='cheby2', rs=25))
            raw.filter(l_freq=None, h_freq=45, method='iir', iir_params=dict(order=70, ftype='cheby2', rs=30))
            raw.notch_filter(freqs=[60, 120], method="fir")
            # 250Hz check
            if raw.info['sfreq'] != 250:
                raw.resample(250)

            def detect_bad_channels_kurtosis(raw, threshold=6): # Kurtosis threshold from Artoni & Michel 2025
                data = raw.get_data(picks="eeg")  # Pick only EEG channels (exclude EOG)
                kurt_vals = kurtosis(data, axis=1, fisher=False)  # Compute kurtosis for each channel
                # Use IQR method for robust thresholding
                Q1, Q3 = np.percentile(kurt_vals, [25, 75])
                iqr = Q3 - Q1
                upper_bound = Q3 + (threshold * iqr)
                bad_channels = [raw.ch_names[i] for i in range(len(kurt_vals)) if kurt_vals[i] > upper_bound]
            
                return bad_channels

            # Detect bad channels and update info
            bad_channels = detect_bad_channels_kurtosis(raw)
            raw.info['bads'] = bad_channels
            # Interpolating bad channels using Splines
            if raw.info['bads']:
                print(f"---------- ----------> Identified bad channels: {bad_channels} <---------- ---------- ")
                raw.interpolate_bads()
            else:
                print(f"\nNo bad channels detected. Skipping interpolation.\n")
            
            # Reference to average
            raw.set_eeg_reference("average")
            
            # --- Picard ICA --- (Converges better than fastICA)
            ica = mne.preprocessing.ICA(n_components=40, max_iter=2000, method='picard', random_state=42)
            
            # Run ICA on cleaned data
            ica.fit(raw)
            
            # Detect EOG-related artifacts using the 3 EOG channels AFTER fitting ICA
            eog_indices, eog_scores = ica.find_bads_eog(raw, ch_name=["EOG1", "EOG2", "EOG3"])
            
            # Apply a temporary 1-100 Hz filter for ICLabel classification (Error pop up even though there is a 1-100Hz filter)
            raw_for_iclabel = raw.copy()
            raw_for_iclabel.filter(l_freq=1.0, h_freq=100.0, method="iir", iir_params=dict(order=8, ftype='cheby2', rs=30))
            
            # Run ICLabel
            ic_labels = label_components(raw_for_iclabel, ica, method="iclabel")
            labels = ic_labels["labels"]
            scores = ic_labels["y_pred_proba"]
            
            # Identify and exclude high-confidence artifacts
            artifact_indices = [
                i for i, label in enumerate(labels)
                if label in ["muscle", "heart", "line_noise", "channel_noise", "other"]
                and scores[i] > 0.80]
            # Print detected artifacts
            if artifact_indices:
                print(f"Detected Artifacts:")
                for i in artifact_indices:
                    print(f"   - Component {i}: {labels[i]} (Confidence: {scores[i]:.2f})")
            else:
                print(f"No Artifacts Detected")
           
            # Combine artefacts with EOG components
            ica.exclude = list(set(eog_indices + artifact_indices))
            
            # Apply ICA correction
            raw = ica.apply(raw)

            # --- Apply Spatial Filtering to Enhance Microstates at the individual level 
            raw = apply_spatial_filter(raw, exclude_bads=True)
            
            subject_data[subject_id] = {
                "raw": raw,
                "bad_channels": bad_channels,
                "ICA_components": ica,
                "labels": ic_labels["labels"],}
            
            # Save raw data
            raw_save_path = os.path.join(save_dir, f"subject_{subject_id}_raw.fif") <--- --- --- --- 
            raw.save(raw_save_path, overwrite=True)
            if subject_id in HS:
                group_data["HS"].append(subject_data[subject_id])
            elif subject_id in LS:
                group_data["LS"].append(subject_data[subject_id])
            else:
                group_data["MS"].append(subject_data[subject_id])
        else:
            print(f"Skipping {file_name}: 'F' key not found in .mat file.")
    else:
        print(f"Skipping {file_name}: File not found.")




More images:



Even more images:


Hello,

The signal traces looks wierd, not just because of the timescale. ICA if a fine denoising method to apply before any kind of other processing, e.g. microstates analysis.
I would suggest to start by using MNE’s default for the preprocessing steps:

raw.filter(1, 45)
raw.notch_filter([60, 120])

Then skip the downsampling, keep the native sampling rate for now.
For bad channels, try first to manually label them instead of relying on automation. Don’t interpolate them, just let them marked as bads for now.
Similarly for the ICA, start by manually excluding ICs once before relying on automation. For instance, you automation is removing ‘line_noise’ components despite the notch filters, that’s not necessaray and might actually label as ‘line_noise’ wrong components.

Finally for the spatial filter, try with and without.
Once you get something sensible on one recording (e.g. the one you shared) through manual operation/minimal processing; then you can try to add back the automation/different elements.

Mathieu

Hi,

I think you’re one of the developers of the package, so thanks a lot for your work! :saluting_face:

I’ve been trying to alter the preprocessing steps for the past few days, adding steps one at a time. It seems like spatial filtering is making a big difference on the data and the plots. The two EEG signal plots are before and after spatial filtering, followed by the before and after n= 4 cluster plots for a single individual.

The output is really strange, it seems like it’s clustering at the back of the head, on the back electrodes. Then, when i apply spatial filtering, it messes a lot with the data, but there are still vivid electrodes, on 2 out of 4 topomaps, in very specific locations that are not at the back of the head. I’m starting to doubt that the montage could be causing this mess. The 2D view of the sensor plot doesn’t make sens, but the 3D view show the cap going around the back of the head.

Have you seen such strange clustering before? Everystep of preprocessing i was doing wasn’t improving much. Only the actual spatial filtering. Should i change the montage or get rid of some sensors? What do you think?



more photos:



The traces on the first plot looks good to me; on the second not so much but my experience with spatial filtering is limited. I would honestly not include it for now, you should be able to get decent maps without.

Montage could be an issue. I’m not sure I followed what you wrote about this; but you should indeed confirm that the location of the individual electode is correct.

Mathieu

Hi Mathieu,

The montage was the problem. I mixed up between two very similar Brain Products EEG caps which caused the order of channels to be completely messed up.

Thanks a lot for your help!



1 Like

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.