Which function should I use for resting-state EEG connectivity estimation?

Support & Discussions
, ,
1
  • MNE version: 1.90
  • operating system: Windows 10
  • MNE-Connectivity version: 0.70

Hello MNE community,

I’m currently working on a project analyzing EEG data from patients with a specific disorder and a healthy control group.

My plan is to use connectivity matrix heatmaps as features and feed them into a neural network model to classify the two groups. To this end, I referred to the following tutorial when writing my code:
https://mne.tools/mne-connectivity/dev/auto_examples/compare_connectivity_over_time_over_trial.html

Here’s my situation and question:

I am using eyes-closed resting-state EEG data. For each subject, I extracted 360 seconds of EEG and segmented it into 120 non-overlapping epochs of 3 seconds each.

In the tutorial above, it is mentioned that:

  • For repeated time-locked events such as ERP data, mne_connectivity.spectral_connectivity_epochs() is recommended.

  • For resting-state data, mne_connectivity.spectral_connectivity_time() might be more appropriate.

Also, I found the following comment in the Support & Discussion section:
https://mne.discourse.group/t/valueerror-occurs-while-computing-spectral-connectivity-time-over-epoched-data/8531/2

“estimating connectivity from a single epoch will give very noisy estimates.”

This left me somewhat confused. So my questions are:

  1. In my case, which function should I use to estimate connectivity? If I should use mne_connectivity.spectral_connectivity_epochs(), how can I estimate connectivity across the 120 individual epochs? (The tutorial above doesn’t seem to cover this.)

  2. If I use mne_connectivity.spectral_connectivity_time() to estimate connectivity from a single epoch, will that cause a serious problem?

I’d really appreciate guidance on what would be the most appropriate approach for my case.

Below is my code:

file_path = "E:\EEG\sub001_task-eyesclosed_eeg.set"
raw = mne.io.read_raw_eeglab(file_path, preload=True)

tmin = 60
tmax = 420
raw_cropped = raw.copy().crop(tmin=tmin, tmax=tmax)


epoch_duration = 3.0  
overlap = 0

epochs = mne.make_fixed_length_epochs(raw_cropped, duration=epoch_duration, overlap=overlap, preload = True)
epochs.load_data()

# Freq bands of interest
Freq_Bands = {"theta": [4.0, 8.0], "alpha": [8.0, 13.0], "beta": [13.0, 30.0]}
n_freq_bands = len(Freq_Bands)
min_freq = np.min(list(Freq_Bands.values()))
max_freq = np.max(list(Freq_Bands.values()))

# Provide the freq points
freqs = np.linspace(min_freq, max_freq, int((max_freq - min_freq) * 4 + 1))

# The dictionary with frequencies are converted to tuples for the function
fmin = tuple([list(Freq_Bands.values())[f][0] for f in range(len(Freq_Bands))])
fmax = tuple([list(Freq_Bands.values())[f][1] for f in range(len(Freq_Bands))])

# We will try two different connectivity measurements as an example
connectivity_methods = ["wpli", "plv",]
n_con_methods = len(connectivity_methods)

sfreq = epochs.info["sfreq"]

# Pre-allocatate memory for the connectivity matrices
con_time_array = np.zeros(
    (n_con_methods, 120, 19, 19, n_freq_bands)
)
con_time_array[con_time_array == 0] = np.nan  # nan matrix

# Compute connectivity over time
con_time = spectral_connectivity_time(
    epochs,
    freqs,
    method=connectivity_methods,
    sfreq=sfreq,
    mode="cwt_morlet",
    fmin=fmin,
    fmax=fmax,
    faverage=True,
)

# Get data as connectivity matrices
for c in range(n_con_methods):
    con_time_array[c] = con_time[c].get_data(output="dense")

ch_names = epochs.ch_names
foi = list(Freq_Bands.keys()).index("alpha")

def plot_con_matrix(con_data, n_con_methods, epoch_num):
    """Visualize the connectivity matrix for a specific epoch."""
    fig, ax = plt.subplots(1, n_con_methods, figsize=(6 * n_con_methods, 6))
    fig.suptitle(f"Epoch {epoch_num} Connectivity", fontsize=16)
    
    for c in range(n_con_methods):
        con_plot = ax[c].imshow(con_data[c, :, :, foi], cmap="jet", vmin=0, vmax=1)
        ax[c].set_title(connectivity_methods[c])
        fig.colorbar(con_plot, ax=ax[c], shrink=0.7, label="Connectivity")
        ax[c].set_xticks(range(len(ch_names)))
        ax[c].set_xticklabels(ch_names, rotation=45)
        ax[c].set_yticks(range(len(ch_names)))
        ax[c].set_yticklabels(ch_names)
        
    plt.tight_layout(rect=[0, 0, 1, 0.96])
    plt.show() 


    plot_con_matrix(con_data_epoch_i, n_con_methods, i)

I’d really appreciate your help. Thank you!

2

Hi @Ananta,

Just to clarify, you want to have connectivity values for each of your 120 epochs? In this case, you should use spectral_connectivity_time() and keep the average parameter as False (the default value). The 3rd code snippet you showed looks good.

Regarding connectivity estimates from single epochs being noisy, this will always be a problem to some extent. However, passing all of your epochs to spectral_connectivity_time() and getting epoch-wise connectivity estimates will be preferable to passing each epoch at a time to spectral_connectivity_epochs().

Hope this helps!

3

@tsbinns Thank you so much for your quick reply and helpful guidance!

I have one more question, if you don’t mind.

In EEG analysis we often deal with multiple frequency bands such as delta, theta, and alpha. When using methods like the Morlet wavelet, I’m wondering how the n_cycles parameter is typically handled:

  • Is it considered more rigorous to set different n_cycles values optimized for each frequency band?

  • Or is it acceptable to choose a single n_cycles value and apply it across all bands?

I’d be very grateful to hear your thoughts on this!

1 Like
4

Hi @Ananta,

This is one of those things that there probably isn’t a definitive answer to unfortunately (see this thread). It may be more appropriate to scale n_cycles with frequency (more smoothing may be desired at higher frequencies), but choosing a single n_cycles value for every frequency band isn’t wrong.

This example from FieldTrip discusses how n_cycles affects smoothing. That first thread also describes a strategy for identifying reasonable values, but it’s trial and error. You’ll need to find something that looks reasonable based on the characteristics you’re trying to study.

Sorry I can’t give a more concrete answer.

5

@tsbinns Your explanation was very clear and helpful — thank you so much for taking the time to guide me!

1 Like