#### Describe the new feature or enhancement
Hi all,
mne-tools/mne-python#79…37 indicated there was an error in the PLI implementation and I also had the same problem when working with resting-state epoched data. But after looking more into it I found out it is not really an error but instead two different ways to analyze connectivity: over time or over trials.
Mike X Cohen explains the difference quite nicely in his video here:
[Youtube Link](https://youtu.be/2h2X4TVUF5I?t=554)
I also commented in that issue, but are now creating a standalone feature request.
The current PLI and PLV implementation in spectral_connectivity() is connectivity over trials, which are suited for ERP data. But for resting-state data it does not make sense to analyze connectivity over trials but should instead be over time.
#### Describe your proposed implementation
I tried implementing my own version to calculate PLV and wPLI over time and you can see a simple working example here:
```python
# Libraries
import numpy as np
import mne
from scipy import special
from mne.connectivity import spectral_connectivity
# Generate data
np.random.seed(42)
n_epochs = 5
n_channels = 3
n_times = 1000
data = np.random.rand(n_epochs, n_channels, n_times)
# Set sampling freq
sfreq = 250 # A reasonable random choice
# Choose what kind of data should be used for evaluation
# 0 = same trial repeated in all epochs
# 1 = 10Hz sinus waves with random phase differences in each channel and epoch
data_option = 1
if data_option == 0:
# Make all 5 epochs the same trial to show difference between connectivity
# over time and over trials. Here we expect con over trials = 1
for i in range(n_epochs):
data[i] = data[0]
elif data_option == 1:
# Generate 10Hz sinus waves to show difference between connectivity
# over time and over trials. Here we expect con over time = 1
for i in range(n_epochs):
for c in range(n_channels):
wave_freq = 10
epoch_len = n_times/sfreq
# Introduce random phase for each channel
phase = np.random.rand(1)*10
# Generate sinus wave
x = np.linspace(-wave_freq*epoch_len*np.pi+phase,
wave_freq*epoch_len*np.pi+phase,n_times)
data[i,c] = np.squeeze(np.sin(x))
else:
print("Data_option value chosen is invalid")
# Define freq bands
Freq_Bands = {"delta": [1.25, 4.0],
"theta": [4.0, 8.0],
"alpha": [8.0, 13.0],
"beta": [13.0, 30.0],
"gamma": [30.0, 49.0]}
n_freq_bands = len(Freq_Bands)
# Convert to tuples for the mne 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))])
# Connectivity methods
connectivity_methods = ["plv","wpli"]
n_con_methods=len(connectivity_methods)
# Number of pairwise ch connections
n_ch_connections = scipy.special.comb(n_channels,2, exact=True, repetition=False)
# Pre-allocatate memory
con_data = np.zeros((n_con_methods,n_channels,n_channels,n_freq_bands))
con_data[con_data==0] = np.nan # nan matrix as 0 is meaningful
# Calculate PLV and wPLI - the MNE python implementation is over trials
con, freqs, times, n_epochs, n_tapers = spectral_connectivity(
data, method=connectivity_methods,
mode="multitaper", sfreq=sfreq, fmin=fmin, fmax=fmax,
faverage=True, verbose=0)
# Save the results in array
con_data[0,:,:,:] = con[0] # PLV
con_data[1,:,:,:] = con[1] # WPLI
print("Alpha PLV over trials")
print(con_data[0,:,:,2]) # all 1, since it is the same trial
# so PLV across trials should be 1
# PLV and wPLI across time
# Make linspace array for morlet waves
freq_centers = np.arange(fmin[0],fmax[-1]+0.25,0.25)
# Prepare Morlets
morlets = mne.time_frequency.tfr.morlet(sfreq,freq_centers,n_cycles=3)
# Make freqs array for indexing
freqs0 = [0]*n_freq_bands
for f in range(n_freq_bands):
freqs0[f] = freq_centers[(freq_centers>=fmin[f]) & (freq_centers<=fmax[f])]
def calculate_PLV_WPLI_across_time(cwt_data):
"""
Parameters
----------
cwt_data : array, shape(n_channels, n_freq, n_times)
The continuous wavelet transform of the data in one epoch
Returns
-------
con_array : array, shape(2, n_channels, n_channels, n_freq_bands)
The connectivity matrix. Only the lower diagonal is calculated.
First axis indicates whether it is PLV (0) or wPLI(1)
"""
n_ch, n_freq0, n_time0 = cwt_data.shape
# The wavelet transform coefficients are complex numbers
# The real part correspond to the amplitude and the imaginary part can be used to calculate the phase
angles_data = np.apply_along_axis(np.angle,2,cwt_data)
# Prepare array with phase differences between all combinations
phase_diff_arr = np.zeros((n_ch_connections,n_freq0,n_time0))
phase_diff_arr_ch_idx = [0]*n_ch_connections
counter = 0
for ch_c in range(n_ch):
for ch_r in range(ch_c+1,n_ch): # only calculate lower diagonal
phase_diff_arr[counter] = angles_data[ch_r]-angles_data[ch_c]
phase_diff_arr_ch_idx[counter] = [ch_r,ch_c]
counter += 1
del angles_data # free up some memory
# =========================================================================
# PLV over time correspond to mean across time of the absolute value of
# the circular length of the relative phases. So PLV will be 1 if
# the phases of 2 signals maintain a constant lag
# In equational form: PLV = 1/N * |sum(e^i(phase1-phase2))|
# =========================================================================
# Convert phase difference to complex part i(phase1-phase2)
phase_diff_arr_im = 0*phase_diff_arr+1j*phase_diff_arr
# Take the exponential
expPhase_arr = np.apply_along_axis(np.exp,2,phase_diff_arr_im)
# Take mean and then the absolute value
meanexpPhase_arr = np.apply_along_axis(np.mean,2,expPhase_arr)
PLV_arr = np.apply_along_axis(np.abs,1,meanexpPhase_arr)
del phase_diff_arr_im, expPhase_arr # free up some memory
# =========================================================================
# PLI over time correspond to the sign of the sine of relative phase
# differences. So PLI will be 1 if one signal is always leading or
# lagging behind the other signal. But it is insensitive to changes in
# relative phase, as long as it is the same signal that leads.
# If 2 signals are almost in phase, they might shift between lead/lag
# due to small fluctuations from noise. This would lead to unstable
# estimation of "phase" synchronisation.
# The wPLI tries to correct for this by weighting the PLI with the
# magnitude of the lag, to attenuate noise sources giving rise to
# near zero phase lag "synchronization"
# In equational form: WPLI = |E{|phase_diff|*sign(phase_diff)}| / E{|phase_diff|}
# =========================================================================
# Calculate the magnitude of phase differences
phase_diff_mag_arr = np.abs(np.sin(phase_diff_arr))
# Calculate the signed phase difference (PLI)
sign_phase_diff_arr = np.sign(np.sin(phase_diff_arr))
# Calculate the nominator (abs and average across time)
WPLI_nominator = np.abs(np.mean(phase_diff_mag_arr*sign_phase_diff_arr,axis=2))
# Calculate denominator for normalization
WPLI_denom = np.mean(phase_diff_mag_arr, axis=2)
# Calculate WPLI
WPLI_arr = WPLI_nominator/WPLI_denom
del phase_diff_mag_arr, sign_phase_diff_arr, phase_diff_arr # free up some memory
# Calculate mean for each freq band
con_array0 = np.zeros((2,n_ch_connections,n_freq_bands))
for f in range(n_freq_bands):
freq_of_interest = freqs0[f] # find freqs in the freq band of interest
freq_idx = [i in freq_of_interest for i in freq_centers] # get idx
con_array0[0,:,f] = np.apply_along_axis(np.mean,1,PLV_arr[:,freq_idx])
con_array0[1,:,f] = np.apply_along_axis(np.mean,1,WPLI_arr[:,freq_idx])
# Save to final array with ch-ch connectivity in matrix form
con_array = np.zeros((2,n_ch,n_ch,n_freq_bands))
for com in range(n_ch_connections):
ch_r = phase_diff_arr_ch_idx[com][0] # get row idx
ch_c = phase_diff_arr_ch_idx[com][1] # get col idx
con_array[0,ch_r,ch_c,:] = con_array0[0,com]
con_array[1,ch_r,ch_c,:] = con_array0[1,com]
return con_array
# Pre-allocate memory
con_data_time = np.zeros((n_con_methods,n_channels,n_channels,n_freq_bands))
con_data_time[con_data_time==0] = np.nan # nan matrix as 0 is meaningful
con_data1 = np.zeros((n_con_methods,n_epochs,n_channels,n_channels,n_freq_bands))
for epoch in range(n_epochs):
# First the data in each epoch is retrieved
temp_data = data[i]
# Then continuous wavelet transform is used to decompose in time frequencies
temp_data_cwt = mne.time_frequency.tfr.cwt(temp_data,morlets)
# PLV and WPLI value is calculated across timepoints in each freq band and averaged into the 5 defined freq bands
PLV_WPLI_con = calculate_PLV_WPLI_across_time(temp_data_cwt)
# Save results
con_data1[0,epoch,:,:,:] = PLV_WPLI_con[0] # phase locking value
con_data1[1,epoch,:,:,:] = PLV_WPLI_con[1] # weighted phase lag index
# Take average across epochs for PLV and wPLI
con_data2 = np.mean(con_data1,axis=1)
# Save to final array
con_data_time[0,:,:,:] = con_data2[0] # phase locking value
con_data_time[1,:,:,:] = con_data2[1] # weighted phase lag index
print("Alpha PLV over time")
print(con_data_time[0,:,:,2])
```
Notice I added the option to choose what kind of data should be used. If data_option = 0 then the expected output is that PLV across trials are 1 since I just repeated the first epoch in all other epochs, while if data_option = 1 the PLV across time should be 1 as I use 10Hz sine waves with different phase differences.
It is still quite rough and not optimized but should give an idea of the procedure.
I am not completely certain about my implementation but the wPLI is calculated by following equation 1 in [1] and weighing by the absolute magnitude of the phase differences as described by [2]. PLV was calculated according to the equation in [3]. According to mne-tools/mne-python#8305 this equation should be equivalent to what is being used in spectral_connectivity() already, but my implementation over time gives another PLV than spectral_connectivity().
It would be helpful if anyone could have a look. After optimizing it the feature should probably just be added as an option to spectral_connectivity() and does not need its own version.
#### Describe possible alternatives
I don't know any alternatives as this suggestion is just about a slight different way of calculating the connectivity over time instead of over trials.
[1] Hardmeier, Martin, Florian Hatz, Habib Bousleiman, Christian Schindler, Cornelis Jan Stam, and Peter Fuhr. 2014. “Reproducibility of Functional Connectivity and Graph Measures Based on the Phase Lag Index (PLI) and Weighted Phase Lag Index (WPLI) Derived from High Resolution EEG.” PLoS ONE 9 (10). https://doi.org/10.1371/journal.pone.0108648.
[2] Vinck, Martin, Robert Oostenveld, Marijn Van Wingerden, Franscesco Battaglia, and Cyriel M.A. Pennartz. 2011. “An Improved Index of Phase-Synchronization for Electrophysiological Data in the Presence of Volume-Conduction, Noise and Sample-Size Bias.” NeuroImage 55 (4): 1548–65. https://doi.org/10.1016/j.neuroimage.2011.01.055.
[3] Lachaux, Jean Philippe, Eugenio Rodriguez, Jacques Martinerie, and Francisco J. Varela. 1999. “Measuring Phase Synchrony in Brain Signals.” Human Brain Mapping 8 (4): 194–208. https://doi.org/10.1002/(SICI)1097-0193(1999)8:4<194::AID-HBM4>3.0.CO;2-C.
