# -*- coding: utf-8 -*-
"""
Updated Aug 7 2024
@author: Qianliang Li (glia@dtu.dk)
This is the main Python file containing the code that support the findings of
https://doi.org/10.1101/2024.05.06.592342
The data used in this analysis was previously described and preprocessed
by Zimmermann, M., Lomoriello, A. S., and Konvalinka, I.
Intra-individual behavioural and neural signatures of audience effects and
interactions in a mirror-game paradigm. Royal Society Open Science, 9(2) 2022
"""
# %% Load libraries
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import mne
import pickle
import mat73
import time
import seaborn as sns
import nolds
from tqdm import tqdm # progress bar
# import Python script for microstates [von Wegner & Lauf, 2018]
# originally downloaded from https://github.com/Frederic-vW/eeg_microstates
# I modified the script for estimating two-brain microstates
# by defining kmeans_return_all and kmeans_dualmicro
from eeg_microstates3 import (kmeans_return_all, kmeans_dualmicro)
# import helper functions
from helper import (numpy_arr_to_pandas_df, time_now)
from dualmicro_functions import (load_epoch_from_fieldtrip, prepare_1P_micro_arr,
plot_microstates, reorder_microstate_results,
single_micro_fit_all_feature_computation,
interbrain_microstate_feature_computation,
prepare_2P_micro_arr_collapsed_events,
plot_dualmicro, sign_swap_microstates,
dualmicro_fit_all_feature_computation,
load_microstate_arrays,
get_synch_events_from_pseudo_pairs,
combine_two_person_microstate_arrays,
pseudo_pair_dualmicro_backfitting,
dualmicro_fit_all_pseudo_pair_feature_computation,
compute_dualmicro_DFA, compute_dualmicro_DFA_pseudo,
shifted_interbrain_microstate_feature_computation)
# Style for matplotlib/seaborn
plt.style.use('default')
# Root for project
os.chdir("C:/Users/glia/Documents/MirrorGame")
# Paths
data_path = "./data/external/EEG/"
mov_data_path = "./data/external/movement/"
fig_save_path = "./reports/figures/"
feat_save_path = "./data/features/"
microstate_save_path = "./data/features/microstates2/"
mov_save_path = "./data/features/movement/"
# %% Load preprocessed EEG data
# The data was originally preprocessed in Fieldtrip by Marius Zimmermann
# Get filenames for the EEG data
files = []
for r, d, f in os.walk(data_path):
for file in f:
if (".mat" in file) & ("ppn" in file):
files.append(os.path.join(r, file))
# Sort the filenames
files.sort()
n_subjects = len(files)
# Get Subject_id
Subject_id = [0]*n_subjects
for i in range(n_subjects):
id_number = files[i].split("/")[-1].split(".")[0].split("pair")[-1].split("pair")[-1].replace("_ppn","")
Subject_id[i] = int(id_number)+1000 # add 1000 to keep the first 0
# There are data from 23 pairs
# Pair 21 and 25 were excluded in the original analysis
# After looking at the data, it seems pair 21, participant 1 and pair 25
# participant 2 only had 1254 and 1440 epochs respectively.
# Their data also do not end with resting-state condition
# All the other EEG data have around 2400 1s epochs and start and ends with rest
bad_subjects = [1211, 1212, 1251, 1252] # the whole pair is dropped
good_subject_idx = [not i in bad_subjects for i in Subject_id]
# Update Subject_id and files
Subject_id = list(np.array(Subject_id)[good_subject_idx])
n_subjects = len(Subject_id)
files = list(np.array(files)[good_subject_idx])
Pair_id = [0]*(n_subjects//2)
for i in range(n_subjects//2):
Pair_id[i] = int(str(Subject_id[2*i])[1:-1])
# Add 100 to pair_id to fix sorting for 1 digit numbers, e.g. 03
Pair_id = [ele+100 for ele in Pair_id]
n_pairs = len(Pair_id)
# Save the IDs as environmental variables to be used in functions
# from dualmicro_functions.py
os.environ["Subject_id"] = Subject_id
os.environ["Pair_id"] = Pair_id
event_id = {"rest":1, "uncoupled":2, "coupled": 3, "observe, actor": 4,
"observe, observer": 6, "imitate, leader": 5, "imitate, follower": 7,
"control": 8}
# Clarification of the labels
# Cond4: ppn1 is observer, ppn2 is actor
# Cond6: ppn1 is actor, ppn2 is observer
# Cond5: ppn1 is follower, ppn2 is leader
# Cond5: ppn1 is leader, ppn2 is follower
event_id_inv = {v: k for k, v in event_id.items()}
# We collapsed condition 4 and 6 & 5 and 7 for two-brain microstates
# By swapping the EEG of ppn1 and ppn2 so ppn1 is always observer/follower and
# ppn2 actor/leader
collapsed_event_id = {"rest":1, "uncoupled":2, "coupled": 3,
"observer_actor": 4, "follower_leader": 5, "control": 8}
collapsed_event_id_inv = {v: k for k, v in collapsed_event_id.items()}
# Load the first EEG to get info about sfreq and n_channels
i = 0
epoch, trialinfo = load_epoch_from_fieldtrip(0, files, event_id)
n_channels = epoch.info["nchan"]
sfreq = int(epoch.info["sfreq"])
# Visualize the data
# epoch.plot(scalings=40e-6, n_channels=32)
# mne.viz.plot_events(epoch.events, sfreq = 1, event_id = event_id, first_samp=-3) # sfreq set to epoch length in s to reflect experiment time
# We compute microstates for the three frequency ranges
alpha_range = [8.0, 13.0]
beta_range = [13.0, 30.0]
broadband_range = None # Data is already 1 to 40 Hz broadband filtered
freq_names = ["alpha","beta","broadband"]
all_freq_ranges = [alpha_range, beta_range, broadband_range]
# %% Intrabrain microstates fit all data
# All subjects from all pairs are concatenated to find common microstates
single_brain_event_id = {"rest":1, "uncoupled":2, "coupled": 3, "observer": 4,
"actor": 6, "follower": 5, "leader": 7, "control": 8}
ppn2_correction = {6:4, 4:6, 7:5, 5:7}
# Loop over frequencies
for f in len(all_freq_ranges):
ff = freq_names[f]
freq_range0 = all_freq_ranges[f]
# =========================================================================
# First the microstate topographies are determined
# It might be an advantage to run the estimation of microstates on a HPC
# =========================================================================
# Get data from all pairs before performing kmeans
np.random.seed(1234)
n_clusters=[3, 4, 5, 6, 7, 8, 9, 10]
n_runs = 100 # increased to 100 runs!
# Get current time
c_time1 = time_now(); print(c_time1)
# Save RAM by appending directly to array instead of making list and then array
sub_arr_indices = [0]
trialinfo_list = []
for i in range(n_subjects):
tmp_data, trialinfo = prepare_1P_micro_arr(i, ppn2_correction, sfreq,
freq_range=freq_range0, standardize=True)
sub_arr_indices.append(len(tmp_data))
trialinfo_list.append([Subject_id[i],trialinfo])
if i == 0: # first run initiation
micro_data_all = tmp_data
else:
micro_data_all = np.append(micro_data_all, tmp_data, axis=0)
del tmp_data # clear up space
print(f"Finished preparing microstate data for pair {Subject_id[i]}")
# Use cumulative sum to determine indices for each subjects's data
subject_indices = np.cumsum(sub_arr_indices)
# Save the trialinfos from all subjects, for easier access in later steps
with open(f"{microstate_save_path}Intrabrain_microstate_fit_all_{ff}_trialinfos.pkl", "wb") as filehandle:
pickle.dump(trialinfo_list, filehandle)
# # with args parser in hpc
# n_maps = n_clusters[(args.map_idx-1)]
# print(f"Running analysis for maps: {n_maps}")
# print("Memory used by the micro data array (GB):",micro_data_all.nbytes*9.31e-10)
# Run Kmeans
for n_maps in n_clusters: # Don't use for loop on the HPC!
# Run the 100 runs in batches of 10 to save underway in case the job script terminates
best_cv_crit = 9999 # initialize unreasonably high value
for r in range(10):
microstate_results = list(kmeans_return_all(micro_data_all, n_maps,
n_runs=int(n_runs/10),maxiter=1000))
# Overwrite the maps if a lower CV criterion was found for the initiation
if microstate_results[4] < best_cv_crit:
microstate_results.append(subject_indices)
# Save results
with open(f"{microstate_save_path}Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "wb") as filehandle:
pickle.dump(microstate_results, filehandle) # [maps, L, gfp_peaks, gev, cv_min, Subject_id]
print(f"Updated the microstates. Previous best CV: {best_cv_crit}",
f"new best CV criterion : {microstate_results[4]}")
# Update best cv criterion value
best_cv_crit = microstate_results[4]
print(f"Finished sub-run {r+1} out of 10")
print(f"Finished microstate analysis for n_maps = {n_maps}")
print("Started", c_time1, "\nCurrent",time_now())
# =========================================================================
# Evaluate microstates fitted to all data
# =========================================================================
# Get summary results
microstate_summary_results = []
for n_maps in n_clusters:
with open(f"{microstate_save_path}Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
# Also save summary results across n_maps
microstate_summary_results.append([microstate_results[0],microstate_results[3],microstate_results[4]])
# Use CV criterion to estimate best number of microstates
cv_gev_arr = np.zeros((len(n_clusters),2))
for imap in range(len(n_clusters)):
gev = np.sum(microstate_summary_results[imap][1])
cv = microstate_summary_results[imap][2]
cv_gev_arr[imap,:] = [cv, gev]
# Convert to Pandas dataframe
col_names = ["n_Microstates", "Fit_Criteria", "Value"]
Fit_Criteria = ["CV Criterion", "Global Explained Variance"]
dtypes = [int,str,"float64"]
cv_gev_df = numpy_arr_to_pandas_df(cv_gev_arr, col_names = col_names,
col_values = [n_clusters,Fit_Criteria],
dtypes = dtypes)
# Evaluate optimal n_Microstates
h_order = Fit_Criteria
g = sns.FacetGrid(data=cv_gev_df,row=None,
margin_titles=True, height=8, aspect=1.2)
g = g.map(sns.pointplot,"n_Microstates", "Value", "Fit_Criteria",
dodge=0, capsize=0.18, errorbar=None, linestyles=["-", "-"],
markers=["o", "o"], hue_order=h_order, palette=sns.color_palette())
g.add_legend()
plt.subplots_adjust(top=0.9, right=0.85, left=0.1)
g.fig.suptitle("Mean CV Criterion and GEV", fontsize=18)
g.set_axis_labels(x_var="Number of Microstates",
y_var="GEV and CV",
fontsize=14)
# The lower CV the better. Measure of residual variance
# But the higher GEV the better.
# Save file
g.savefig(f"{fig_save_path}Microstates/Fit_all_{ff}/"+"Single_micro_fit_all_{ff}_CV_Criterion_GEV"+".png")
# Count which number of microstates have the lowest cv criterion for each subject
min_idx = np.argmin(cv_gev_df.loc[cv_gev_df["Fit_Criteria"]=="CV Criterion","Value"])
cv_gev_df.loc[cv_gev_df["Fit_Criteria"]=="CV Criterion"].iloc[min_idx]
# Visualize all microstates prior to re-ordering
for ii in range(len(n_clusters)):
plot_microstates(n_clusters[ii], microstate_summary_results[ii][0], microstate_summary_results[ii][1], epoch.info)
# =========================================================================
# # Re-order intrabrain microstates
# =========================================================================
# This is only run once, after microstates are created
# The optimal number of microstates were 5, with 56% GEV
n_maps = 5
ii = n_clusters.index(n_maps)
with open(f"{microstate_save_path}Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
maps, m_labels, gfp_peaks, gev, cv_min, sub_idx = microstate_results
plot_microstates(n_maps, maps, gev)
# Make dictionary with n_maps and new order
manual_reordering_template = {"5_alpha":[4,1,3,2,0],
"5_beta":[3,2,1,4,0],
"5_broadband":[3,2,4,1,0]}
new_order = manual_reordering_template[f"{n_maps}_{ff}"]
# Re-order the microstates
maps, gev, m_labels = reorder_microstate_results(new_order, maps, gev, m_labels)
# Plot again to check it worked
plot_microstates(n_maps, maps, gev, epoch.info)
# Since neuronal activity is often oscillating, this causes polarity inversions
# Microstates ignores the sign, and hence the polarity in the map is arbitrary
# It is only the relative difference within the plot that is interesting
# depending on initiation. We can thus freely change the sign for visualization
# For two-person microstates, each person's map is sign-changed separately
manual_sign_correction = {"5_alpha":[1,-1,1,1,1],
"5_beta":[1,1,1,-1,-1],
"5_broadband":[-1,1,-1,1,1]}
sign_swap = manual_sign_correction[f"{n_maps}_{ff}"]
for m in range(n_maps):
maps[m] *= sign_swap[m]
# Plot a final time for last confirmation
plot_microstates(n_maps, maps, gev, epoch.info)
# Close all figures
plt.close("all")
### Save reordered results
n_maps = 5
ii = n_clusters.index(n_maps)
with open(f"{microstate_save_path}Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
maps, m_labels, gfp_peaks, gev, cv, sub_indices = microstate_results
# Re-order
new_order = manual_reordering_template[str(n_maps)]
maps, gev, m_labels = reorder_microstate_results(new_order, maps, gev, m_labels)
# Sign swap
for m in range(n_maps):
maps[m] *= sign_swap[m]
# Overwrite variable
microstate_results = maps, m_labels, gfp_peaks, gev, cv, sub_indices
# Save to new file
with open(f"{microstate_save_path}Reordered/Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "wb") as filehandle:
pickle.dump(microstate_results, filehandle) # [maps, L, gfp_peaks, gev, cv_min, sub_idx]
# Save topomaps for the microstates
save_path = f"{fig_save_path}Microstates/Fit_all_{ff}/"
with open(f"{microstate_save_path}Reordered/Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
maps, m_labels, gfp_peaks, gev, cv_min, sub_idx = microstate_results
fig = plot_microstates(n_maps, maps, gev, epoch.info)
fig.savefig(save_path+f"Intrabrain_fit_all_{ff}_maps{n_maps}"+".png")
# Save svg for Paper
fig.savefig(save_path+f"Intrabrain_fit_all_{ff}_maps{n_maps}"+".svg")
# =========================================================================
# # Estimate one-person microstate metrics/features
# # There might be a small error introduced due to gaps in the time series from
# # dropped segments, e.g. when calculating the transition probability as
# # the time series is discontinuous due to the gaps. But with the high sampling rate
# # only a very small fraction of the samples have discontinuous neighbors
# =========================================================================
# The observer_actor and observer_observe conditions have been separated
# So there are observer and actor conditions.
# And the same for leader and follower.
"""
Overview of common (intrabrain) microstate features:
1. Average duration a given microstate remains stable (Dur)
2. Frequency occurrence, independent of individual duration (Occ)
Average number of times a microstate becomes dominant per second
3. Ratio of total Time Covered (TCo)
4. Transition probabilities (TMx)
5. Ratio of shannon entropy relative to theoretical max chaos (Ent)
"""
# Hard-coded the optimal number of microstates based on CV criterion and GEV for dualmicro
n_maps = 5
# Load all microstate results
with open(f"{microstate_save_path}Reordered/Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
# Load all trialinfos
with open(f"{microstate_save_path}Intrabrain_microstate_fit_all_{ff}_trialinfos.pkl", "rb") as file:
trialinfo_list = pickle.load(file)
Microstate_names = [chr(ele) for ele in range(65,65+n_maps)]
m_labels = [0]*n_subjects
events = [0]*n_subjects
m_feats = [0]*n_subjects
for i in range(n_subjects):
m_labels[i], events[i], m_feats[i] = single_micro_fit_all_feature_computation(i,
n_maps, microstate_results, trialinfo_list, sfreq, event_id, single_brain_event_id)
print(f"Finished computing microstate features for Subject {Subject_id[i]}")
# Save the raw microstate features
with open(f"{microstate_save_path}/raw_features_single_micro_fit_all_{ff}_maps{n_maps}.pkl", "wb") as filehandle:
pickle.dump(m_feats, filehandle) # [Subject][Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr] [Event, map*]
# * the feature is calculated for each map, where applicable.
# Transition matrix is calculated for each map -> map transition probability
# with open(f"{microstate_save_path}/raw_features_single_micro_fit_all_{ff}_maps{n_maps}.pkl", "rb") as file:
# m_feats = pickle.load(file) # [Subject][Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr] [Event, map*]
### Convert all features to dataframes for further processing
col_names = ["Subject_ID", "Event_ID", "Microstate", "Value"]
col_values = [Subject_id,list(single_brain_event_id.keys()),Microstate_names]
dtypes = ["int64",str,str,"float64"]
# Mean duration
Dur_arr = np.stack([ele[0] for ele in m_feats]) # [Subject, event, n_map]
Dur_df = numpy_arr_to_pandas_df(Dur_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Duration"]*len(Dur_df)
Dur_df.insert(2, "Measurement", measurement_id)
# Save df
Dur_df.to_pickle(os.path.join(microstate_save_path,f"Single_micro_fit_all_{ff}_maps{n_maps}_duration_df.pkl"))
# Frequency of occurrence per sec
Occ_arr = np.stack([ele[1] for ele in m_feats]) # [Subject, event, n_map]
Occ_df = numpy_arr_to_pandas_df(Occ_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Occurrence"]*len(Occ_df)
Occ_df.insert(2, "Measurement", measurement_id)
# Save df
Occ_df.to_pickle(os.path.join(microstate_save_path,f"Single_micro_fit_all_{ff}_maps{n_maps}_occurrence_df.pkl"))
# Ratio total Time Covered
TCo_arr = np.stack([ele[2] for ele in m_feats]) # [Subject, event, n_map]
TCo_df = numpy_arr_to_pandas_df(TCo_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Time_covered"]*len(TCo_df)
TCo_df.insert(2, "Measurement", measurement_id)
# Save df
TCo_df.to_pickle(os.path.join(microstate_save_path,f"Single_micro_fit_all_{ff}_maps{n_maps}_ratio_time_covered_df.pkl"))
# Transition matrix should be read as probability of row to column
xi, xj = np.meshgrid(Microstate_names,Microstate_names)
_, arrow = np.meshgrid(Microstate_names,["->"]*n_maps)
transition_info = np.char.add(np.char.add(xj,arrow),xi)
TMx_arr = np.stack([ele[3] for ele in m_feats]) # [Subject, event, n_map, n_map]
TMx_arr = TMx_arr.reshape((n_subjects,len(single_brain_event_id),n_maps*n_maps)) # Flatten the maps to 1D
col_names = ["Subject_ID", "Event_ID", "Transition", "Value"]
col_values = [Subject_id,list(single_brain_event_id.keys()),transition_info.flatten()]
TMx_df = numpy_arr_to_pandas_df(TMx_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Probability"]*len(TMx_df)
TMx_df.insert(2, "Measurement", measurement_id)
# Save df
TMx_df.to_pickle(os.path.join(microstate_save_path,f"Single_micro_fit_all_{ff}_maps{n_maps}_transition_df.pkl"))
# Entropy
Ent_arr = np.stack([ele[4] for ele in m_feats]) # [Subject, event]
col_names = ["Subject_ID", "Event_ID", "Value"]
col_values = [Subject_id,list(single_brain_event_id.keys())]
dtypes = ["int64",str,"float64"]
Ent_df = numpy_arr_to_pandas_df(Ent_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Entropy"]*len(Ent_df)
Ent_df.insert(2, "Measurement", measurement_id)
# Save df
Ent_df.to_pickle(os.path.join(microstate_save_path,f"Single_micro_fit_all_{ff}_maps{n_maps}_ratio_entropy_df.pkl"))
# =========================================================================
# We also did it for 8 alpha microstates to use the same number as
# the two-brain microstates
# =========================================================================
# This is only run once, after microstates are created
ff = "alpha"
n_maps = 8
ii = n_clusters.index(n_maps)
with open(f"{microstate_save_path}Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
maps, m_labels, gfp_peaks, gev, cv_min, sub_idx = microstate_results
plot_microstates(n_maps, maps, gev)
# Make dictionary with n_maps and new order
manual_reordering_template = {"8":[6,0,5,1,7,2,3,4]}
new_order = manual_reordering_template[str(n_maps)]
# Re-order the microstates
maps, gev, m_labels = reorder_microstate_results(new_order, maps, gev, m_labels)
# Plot again to check it worked
plot_microstates(n_maps, maps, gev, epoch.info)
# Since neuronal activity is often oscillating, this causes polarity inversions
# Microstates ignores the sign, and hence the polarity in the map is arbitrary
# It is only the relative difference within the plot that is interesting
# depending on initiation. We can thus freely change the sign for visualization
# For two-person microstates, each person's map is sign-changed separately
manual_sign_correction = {"8":[-1,1,-1,1,1,1,-1,-1]}
sign_swap = manual_sign_correction[str(n_maps)]
for m in range(n_maps):
maps[m] *= sign_swap[m]
# Plot a final time for last confirmation
plot_microstates(n_maps, maps, gev, epoch.info)
# Close all figures
plt.close("all")
### Save reordered results
n_maps = 8
ii = n_clusters.index(n_maps)
with open(f"{microstate_save_path}Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
maps, m_labels, gfp_peaks, gev, cv, sub_indices = microstate_results
# Re-order
new_order = manual_reordering_template[str(n_maps)]
maps, gev, m_labels = reorder_microstate_results(new_order, maps, gev, m_labels)
# Sign swap
for m in range(n_maps):
maps[m] *= sign_swap[m]
# Overwrite variable
microstate_results = maps, m_labels, gfp_peaks, gev, cv, sub_indices
# Save to new file
with open(f"{microstate_save_path}Reordered/Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "wb") as filehandle:
pickle.dump(microstate_results, filehandle) # [maps, L, gfp_peaks, gev, cv_min, sub_idx]
# Save topomaps for the microstates
save_path = f"{fig_save_path}Microstates/Fit_all_{ff}/"
with open(f"{microstate_save_path}Reordered/Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
maps, m_labels, gfp_peaks, gev, cv_min, sub_idx = microstate_results
fig = plot_microstates(n_maps, maps, gev, epoch.info)
fig.savefig(save_path+f"Intrabrain_fit_all_{ff}_maps{n_maps}"+".png")
# Save svg for Paper
fig.savefig(save_path+f"Intrabrain_fit_all_{ff}_maps{n_maps}"+".svg")
# =========================================================================
# # Estimate one-person microstate metrics/features
# # There might be a small error introduced due to gaps in the time series from
# # dropped segments, e.g. when calculating the transition probability as
# # the time series is discontinuous due to the gaps. But with the high sampling rate
# # only a very small fraction of the samples have discontinuous neighbors
# =========================================================================
# The observer_actor and observer_observe conditions have been separated
# So there are observer and actor conditions.
# And the same for leader and follower.
"""
Overview of common (intrabrain) microstate features:
1. Average duration a given microstate remains stable (Dur)
2. Frequency occurrence, independent of individual duration (Occ)
Average number of times a microstate becomes dominant per second
3. Ratio of total Time Covered (TCo)
4. Transition probabilities (TMx)
5. Ratio of shannon entropy relative to theoretical max chaos (Ent)
"""
# Hard-coded the optimal number of microstates based on CV criterion and GEV for dualmicro
n_maps = 8
# Load all microstate results
with open(f"{microstate_save_path}Reordered/Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
# Load all trialinfos
with open(f"{microstate_save_path}Intrabrain_microstate_fit_all_{ff}_trialinfos.pkl", "rb") as file:
trialinfo_list = pickle.load(file)
Microstate_names = [chr(ele) for ele in range(65,65+n_maps)]
m_labels = [0]*n_subjects
events = [0]*n_subjects
m_feats = [0]*n_subjects
for i in range(n_subjects):
m_labels[i], events[i], m_feats[i] = single_micro_fit_all_feature_computation(i,
n_maps, microstate_results, trialinfo_list, sfreq, event_id, single_brain_event_id)
print(f"Finished computing microstate features for Subject {Subject_id[i]}")
# Save the raw microstate features
with open(f"{microstate_save_path}/raw_features_single_micro_fit_all_{ff}_maps{n_maps}.pkl", "wb") as filehandle:
pickle.dump(m_feats, filehandle) # [Subject][Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr] [Event, map*]
# * the feature is calculated for each map, where applicable.
# Transition matrix is calculated for each map -> map transition probability
# with open(f"{microstate_save_path}/raw_features_single_micro_fit_all_{ff}_maps{n_maps}.pkl", "rb") as file:
# m_feats = pickle.load(file) # [Subject][Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr] [Event, map*]
### Convert all features to dataframes for further processing
col_names = ["Subject_ID", "Event_ID", "Microstate", "Value"]
col_values = [Subject_id,list(single_brain_event_id.keys()),Microstate_names]
dtypes = ["int64",str,str,"float64"]
# Mean duration
Dur_arr = np.stack([ele[0] for ele in m_feats]) # [Subject, event, n_map]
Dur_df = numpy_arr_to_pandas_df(Dur_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Duration"]*len(Dur_df)
Dur_df.insert(2, "Measurement", measurement_id)
# Save df
Dur_df.to_pickle(os.path.join(microstate_save_path,f"Single_micro_fit_all_{ff}_maps{n_maps}_duration_df.pkl"))
# Frequency of occurrence per sec
Occ_arr = np.stack([ele[1] for ele in m_feats]) # [Subject, event, n_map]
Occ_df = numpy_arr_to_pandas_df(Occ_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Occurrence"]*len(Occ_df)
Occ_df.insert(2, "Measurement", measurement_id)
# Save df
Occ_df.to_pickle(os.path.join(microstate_save_path,f"Single_micro_fit_all_{ff}_maps{n_maps}_occurrence_df.pkl"))
# Ratio total Time Covered
TCo_arr = np.stack([ele[2] for ele in m_feats]) # [Subject, event, n_map]
TCo_df = numpy_arr_to_pandas_df(TCo_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Time_covered"]*len(TCo_df)
TCo_df.insert(2, "Measurement", measurement_id)
# Save df
TCo_df.to_pickle(os.path.join(microstate_save_path,f"Single_micro_fit_all_{ff}_maps{n_maps}_ratio_time_covered_df.pkl"))
# Transition matrix should be read as probability of row to column
xi, xj = np.meshgrid(Microstate_names,Microstate_names)
_, arrow = np.meshgrid(Microstate_names,["->"]*n_maps)
transition_info = np.char.add(np.char.add(xj,arrow),xi)
TMx_arr = np.stack([ele[3] for ele in m_feats]) # [Subject, event, n_map, n_map]
TMx_arr = TMx_arr.reshape((n_subjects,len(single_brain_event_id),n_maps*n_maps)) # Flatten the maps to 1D
col_names = ["Subject_ID", "Event_ID", "Transition", "Value"]
col_values = [Subject_id,list(single_brain_event_id.keys()),transition_info.flatten()]
TMx_df = numpy_arr_to_pandas_df(TMx_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Probability"]*len(TMx_df)
TMx_df.insert(2, "Measurement", measurement_id)
# Save df
TMx_df.to_pickle(os.path.join(microstate_save_path,f"Single_micro_fit_all_{ff}_maps{n_maps}_transition_df.pkl"))
# Entropy
Ent_arr = np.stack([ele[4] for ele in m_feats]) # [Subject, event]
col_names = ["Subject_ID", "Event_ID", "Value"]
col_values = [Subject_id,list(single_brain_event_id.keys())]
dtypes = ["int64",str,"float64"]
Ent_df = numpy_arr_to_pandas_df(Ent_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Entropy"]*len(Ent_df)
Ent_df.insert(2, "Measurement", measurement_id)
# Save df
Ent_df.to_pickle(os.path.join(microstate_save_path,f"Single_micro_fit_all_{ff}_maps{n_maps}_ratio_entropy_df.pkl"))
# %% Inter-brain microstates fit all data
# Based on the microstate topographies estimated on single-brian data
"""
Interbrain features:
1. Average duration of common interbrain microstates (IBDur)
2. Frequency occurrence of common interbrain microstates in the pair (IBOcc)
3. Ratio of total time covered by interbrain common microstates in the pair (IBCov)
4. Transition probability towards common interbrain microstates in the pair (IBTMx)
5. Ratio of joint shannon entropy relative to theoretical max chaos (IBEnt)
"""
for f in len(all_freq_ranges):
ff = freq_names[f]
# Hard-coded the optimal number of microstates based on CV criterion and GEV
n_maps = 5
# Load all microstate results
with open(f"{microstate_save_path}Reordered/Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
# Load all trialinfos
with open(f"{microstate_save_path}Intrabrain_microstate_fit_all_{ff}_trialinfos.pkl", "rb") as file:
trialinfo_list = pickle.load(file)
Microstate_names = [chr(ele) for ele in range(65,65+n_maps)]
# Insert Z as the symbol for non common microstate
Microstate_names.insert(0,"Z")
m_labels = [0]*(n_subjects//2)
events = [0]*(n_subjects//2)
m_feats = [0]*(n_subjects//2)
Pair_id = [0]*(n_subjects//2)
for i in range(n_subjects//2):
m_labels[i], events[i], m_feats[i] = interbrain_microstate_feature_computation(i,
n_maps, microstate_results, trialinfo_list, sfreq, event_id, collapsed_event_id)
Pair_id[i] = int(str(Subject_id[2*i])[1:-1])
print(f"Finished computing interbrain microstate features for pair {Pair_id[i]}")
Pair_id = [ele+100 for ele in Pair_id]
# Save the raw microstate features
with open(f"{microstate_save_path}/raw_interbrain_single_micro_fit_all_{ff}_maps{n_maps}.pkl", "wb") as filehandle:
pickle.dump([Pair_id, m_feats], filehandle) # [Subject][Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr] [Event, map*]
# * the feature is calculated for each map, where applicable.
# Transition matrix is calculated for each map -> map transition probability
# The first row and column correspond to the non common microstate, i.e.
# there is a different microstate in the pair
# with open(f"{microstate_save_path}/raw_interbrain_single_micro_fit_all_{ff}_maps.pkl", "rb") as file:
# Pair_id, m_feats = pickle.load(file) # [Subject][Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr] [Event, map*]
n_pairs = len(Pair_id)
### Convert all features to dataframes for further processing
col_names = ["Pair_ID", "Event_ID", "Microstate", "Value"]
col_values = [Pair_id,list(collapsed_event_id.keys()),Microstate_names]
dtypes = [int,str,str,"float64"]
# Mean duration
Dur_arr = np.stack([ele[0] for ele in m_feats]) # [Subject, event, n_map]
Dur_df = numpy_arr_to_pandas_df(Dur_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Duration"]*len(Dur_df)
Dur_df.insert(2, "Measurement", measurement_id)
# Save df
Dur_df.to_pickle(os.path.join(microstate_save_path,f"IB_Single_micro_fit_all_{ff}_maps{n_maps}_duration_df.pkl"))
# Frequency of occurrence per sec
Occ_arr = np.stack([ele[1] for ele in m_feats]) # [Subject, event, n_map]
Occ_df = numpy_arr_to_pandas_df(Occ_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Occurrence"]*len(Occ_df)
Occ_df.insert(2, "Measurement", measurement_id)
# Save df
Occ_df.to_pickle(os.path.join(microstate_save_path,f"IB_Single_micro_fit_all_{ff}_maps{n_maps}_occurrence_df.pkl"))
# Ratio total Time Covered
TCo_arr = np.stack([ele[2] for ele in m_feats]) # [Subject, event, n_map]
TCo_df = numpy_arr_to_pandas_df(TCo_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Time_covered"]*len(TCo_df)
TCo_df.insert(2, "Measurement", measurement_id)
# Save df
TCo_df.to_pickle(os.path.join(microstate_save_path,f"IB_Single_micro_fit_all_{ff}_maps{n_maps}_ratio_time_covered_df.pkl"))
# Transition matrix should be read as probability of row to column
xi, xj = np.meshgrid(Microstate_names,Microstate_names)
_, arrow = np.meshgrid(Microstate_names,["->"]*(n_maps+1))
transition_info = np.char.add(np.char.add(xj,arrow),xi)
TMx_arr = np.stack([ele[3] for ele in m_feats]) # [Subject, event, n_map, n_map]
TMx_arr = TMx_arr.reshape((n_pairs,len(collapsed_event_id),(n_maps+1)*(n_maps+1))) # Flatten the maps to 1D
col_names = ["Pair_ID", "Event_ID", "Transition", "Value"]
col_values = [Pair_id,list(collapsed_event_id.keys()),transition_info.flatten()]
TMx_df = numpy_arr_to_pandas_df(TMx_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Probability"]*len(TMx_df)
TMx_df.insert(2, "Measurement", measurement_id)
# Save df
TMx_df.to_pickle(os.path.join(microstate_save_path,f"IB_Single_micro_fit_all_{ff}_maps{n_maps}_transition_df.pkl"))
# Entropy
Ent_arr = np.stack([ele[4] for ele in m_feats]) # [Subject, event]
col_names = ["Pair_ID", "Event_ID", "Value"]
col_values = [Pair_id,list(collapsed_event_id.keys())]
dtypes = [int, str, "float64"]
Ent_df = numpy_arr_to_pandas_df(Ent_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Entropy"]*len(Ent_df)
Ent_df.insert(2, "Measurement", measurement_id)
# Save df
Ent_df.to_pickle(os.path.join(microstate_save_path,f"IB_Single_micro_fit_all_{ff}_maps{n_maps}_ratio_joint_entropy_df.pkl"))
# =========================================================================
# Repeat for 8 alpha microstates
# =========================================================================
ff = "alpha"
n_maps = 8
# Load all microstate results
with open(f"{microstate_save_path}Reordered/Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
# Load all trialinfos
with open(f"{microstate_save_path}Intrabrain_microstate_fit_all_{ff}_trialinfos.pkl", "rb") as file:
trialinfo_list = pickle.load(file)
Microstate_names = [chr(ele) for ele in range(65,65+n_maps)]
# Insert Z as the symbol for non common microstate
Microstate_names.insert(0,"Z")
m_labels = [0]*(n_subjects//2)
events = [0]*(n_subjects//2)
m_feats = [0]*(n_subjects//2)
Pair_id = [0]*(n_subjects//2)
for i in range(n_subjects//2):
m_labels[i], events[i], m_feats[i] = interbrain_microstate_feature_computation(i,
n_maps, microstate_results, trialinfo_list, sfreq, event_id, collapsed_event_id)
Pair_id[i] = int(str(Subject_id[2*i])[1:-1])
print(f"Finished computing interbrain microstate features for pair {Pair_id[i]}")
Pair_id = [ele+100 for ele in Pair_id]
# Save the raw microstate features
with open(f"{microstate_save_path}/raw_interbrain_single_micro_fit_all_{ff}_maps{n_maps}.pkl", "wb") as filehandle:
pickle.dump([Pair_id, m_feats], filehandle) # [Subject][Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr] [Event, map*]
# * the feature is calculated for each map, where applicable.
# Transition matrix is calculated for each map -> map transition probability
# The first row and column correspond to the non common microstate, i.e.
# there is a different microstate in the pair
# with open(f"{microstate_save_path}/raw_interbrain_single_micro_fit_all_{ff}_maps.pkl", "rb") as file:
# Pair_id, m_feats = pickle.load(file) # [Subject][Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr] [Event, map*]
n_pairs = len(Pair_id)
### Convert all features to dataframes for further processing
col_names = ["Pair_ID", "Event_ID", "Microstate", "Value"]
col_values = [Pair_id,list(collapsed_event_id.keys()),Microstate_names]
dtypes = [int,str,str,"float64"]
# Mean duration
Dur_arr = np.stack([ele[0] for ele in m_feats]) # [Subject, event, n_map]
Dur_df = numpy_arr_to_pandas_df(Dur_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Duration"]*len(Dur_df)
Dur_df.insert(2, "Measurement", measurement_id)
# Save df
Dur_df.to_pickle(os.path.join(microstate_save_path,f"IB_Single_micro_fit_all_{ff}_maps{n_maps}_duration_df.pkl"))
# Frequency of occurrence per sec
Occ_arr = np.stack([ele[1] for ele in m_feats]) # [Subject, event, n_map]
Occ_df = numpy_arr_to_pandas_df(Occ_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Occurrence"]*len(Occ_df)
Occ_df.insert(2, "Measurement", measurement_id)
# Save df
Occ_df.to_pickle(os.path.join(microstate_save_path,f"IB_Single_micro_fit_all_{ff}_maps{n_maps}_occurrence_df.pkl"))
# Ratio total Time Covered
TCo_arr = np.stack([ele[2] for ele in m_feats]) # [Subject, event, n_map]
TCo_df = numpy_arr_to_pandas_df(TCo_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Time_covered"]*len(TCo_df)
TCo_df.insert(2, "Measurement", measurement_id)
# Save df
TCo_df.to_pickle(os.path.join(microstate_save_path,f"IB_Single_micro_fit_all_{ff}_maps{n_maps}_ratio_time_covered_df.pkl"))
# Transition matrix should be read as probability of row to column
xi, xj = np.meshgrid(Microstate_names,Microstate_names)
_, arrow = np.meshgrid(Microstate_names,["->"]*(n_maps+1))
transition_info = np.char.add(np.char.add(xj,arrow),xi)
TMx_arr = np.stack([ele[3] for ele in m_feats]) # [Subject, event, n_map, n_map]
TMx_arr = TMx_arr.reshape((n_pairs,len(collapsed_event_id),(n_maps+1)*(n_maps+1))) # Flatten the maps to 1D
col_names = ["Pair_ID", "Event_ID", "Transition", "Value"]
col_values = [Pair_id,list(collapsed_event_id.keys()),transition_info.flatten()]
TMx_df = numpy_arr_to_pandas_df(TMx_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Probability"]*len(TMx_df)
TMx_df.insert(2, "Measurement", measurement_id)
# Save df
TMx_df.to_pickle(os.path.join(microstate_save_path,f"IB_Single_micro_fit_all_{ff}_maps{n_maps}_transition_df.pkl"))
# Entropy
Ent_arr = np.stack([ele[4] for ele in m_feats]) # [Subject, event]
col_names = ["Pair_ID", "Event_ID", "Value"]
col_values = [Pair_id,list(collapsed_event_id.keys())]
dtypes = [int, str, "float64"]
Ent_df = numpy_arr_to_pandas_df(Ent_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Entropy"]*len(Ent_df)
Ent_df.insert(2, "Measurement", measurement_id)
# Save df
Ent_df.to_pickle(os.path.join(microstate_save_path,f"IB_Single_micro_fit_all_{ff}_maps{n_maps}_ratio_joint_entropy_df.pkl"))
# %% Two-brain microstates fit all data
"""
The two observe and imitate conditions are collapesed
Instead of having ppn1 being observer/follower in 8 trials and actor/leader
in 8 trials, we will fix the topomap from "ppn1, top row" to always be
observer and follower. This means for condition 6 and 7, ppn2 will be treated
as ppn1 so the first topomap is still being fitted to the observer/follower!
So the first microstate (top row) will always correspond to the Observer and Follower
And the 2nd paired microstate (bot row) will always correspond to Actor and Leader
Additionally we compute features for 8 trials and then take the average instead
of all 16. This is done in order to compute it for the asymmetrical trials
without flipping, as the flip itself can create artefacts.
And the same process is repeated for the symmetrical conditions to be consistent,,
although it shouldn't have a big impact for those trials
"""
# Compute two-person microstates for each pair, fitted for all data
# We will concatenate the pairs along the channel axis
# Loop over frequencies
for f in len(all_freq_ranges):
ff = freq_names[f]
freq_range0 = all_freq_ranges[f]
# =========================================================================
# First the microstate topographies are determined
# It might be an advantage to run the estimation of microstates on a HPC
# =========================================================================
# Get data from all pairs before performing kmeans
np.random.seed(1234)
n_clusters=[3, 4, 5, 6, 7, 8, 9, 10]
n_runs = 100 # increased to 100 runs!
# Get current time
c_time1 = time_now(); print(c_time1)
# Save RAM by appending directly to array instead of making list and then array
pair_arr_indices = [0]
trialinfo_list = []
events_list = []
for i in range(n_pairs):
tmp_data, tmp_trialinfo, tmp_events = prepare_2P_micro_arr_collapsed_events(i, sfreq, event_id, freq_range=freq_range0, standardize=True)
pair_arr_indices.append(len(tmp_data))
trialinfo_list.append(tmp_trialinfo)
events_list.append(tmp_events)
if i == 0: # first run initiation
micro_data_all = tmp_data
else:
micro_data_all = np.append(micro_data_all,tmp_data, axis=0)
del tmp_data # clear up space
print(f"Finished preparing microstate data for pair {Pair_id[i]}")
# Use cumulative sum to determine indices for each pair's data
pair_indices = np.cumsum(pair_arr_indices)
# Save the trialinfos and events from all pairs, for easier access in later steps
with open(f"{microstate_save_path}Dualmicro_fit_all_{ff}_trial_events_infos.pkl", "wb") as filehandle:
pickle.dump([Pair_id,trialinfo_list,events_list], filehandle) # [maps, L, gfp_peaks, gev, cv_min, pair_idx]
# # with args parser in hpc
# n_maps = n_clusters[(args.map_idx-1)]
# print(f"Running analysis for maps: {n_maps}")
# print("Memory used by the micro data array (GB):",micro_data_all.nbytes*9.31e-10)
for n_maps in n_clusters: # Don't use for loop on the HPC!
# Run the 100 runs in batches of 10 to save underway in case the job script terminates
best_cv_crit = 9999 # initialize unreasonably high value
for r in range(10):
microstate_results = list(kmeans_dualmicro(micro_data_all, n_maps,
n_runs=int(n_runs/10),maxiter=1000))
# Overwrite the maps if a lower CV criterion was found for the initiation
if microstate_results[4] < best_cv_crit:
microstate_results.append(pair_indices)
# Save results
with open(f"{microstate_save_path}Dualmicro_fit_all_{ff}_data_maps{n_maps}.pkl", "wb") as filehandle:
pickle.dump(microstate_results, filehandle) # [maps, L, gfp_peaks, gev, cv_min, pair_idx]
print(f"Updated the microstates. Previous best CV: {best_cv_crit}",
f"new best CV criterion : {microstate_results[4]}")
# Update best cv criterion value
best_cv_crit = microstate_results[4]
print(f"Finished sub-run {r+1} out of 10")
print(f"Finished microstate analysis for n_maps = {n_maps}")
print("Started", c_time1, "\nCurrent",time_now())
# =========================================================================
# # Evaluate microstates fitted to all data
# =========================================================================
# Get summary results
microstate_summary_results = []
for n_maps in n_clusters:
with open(f"{microstate_save_path}Dualmicro_fit_all_{ff}_data_maps{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
# Also save summary results across n_maps
microstate_summary_results.append([microstate_results[0],microstate_results[3],microstate_results[4]])
# Use CV criterion to estimate best number of microstates
cv_gev_arr = np.zeros((len(n_clusters),2))
for imap in range(len(n_clusters)):
gev = np.sum(microstate_summary_results[imap][1])
cv = microstate_summary_results[imap][2]
cv_gev_arr[imap,:] = [cv, gev]
# Convert to Pandas dataframe
col_names = ["n_Microstates", "Fit_Criteria", "Value"]
Fit_Criteria = ["CV Criterion", "Global Explained Variance"]
dtypes = [int,str,"float64"]
cv_gev_df = numpy_arr_to_pandas_df(cv_gev_arr, col_names = col_names, col_values = [n_clusters,Fit_Criteria],
dtypes = dtypes)
# Evaluate optimal n_Microstates
h_order = Fit_Criteria
g = sns.FacetGrid(data=cv_gev_df,row=None,
margin_titles=True, height=8, aspect=1.5)
g = g.map(sns.pointplot,"n_Microstates", "Value", "Fit_Criteria",
dodge=0, capsize=0.18, errorbar=None, linestyles=["-", "-"],
markers=["o", "o"], hue_order=h_order, palette=sns.color_palette())
g.add_legend()
plt.subplots_adjust(top=0.9, right=0.85, left=0.1)
g.fig.suptitle("Mean CV Criterion and GEV", fontsize=18)
g.set_axis_labels(x_var="Number of Microstates",
y_var="GEV and CV",
fontsize=14)
# The lower CV the better. Measure of residual variance
# But the higher GEV the better.
# Save file
g.savefig(f"{fig_save_path}Microstates/Fit_all_{ff}/"+"Dualmicro_fit_all_{ff}_CV_Criterion_GEV"+".png")
# Count which number of microstates have the lowest cv criterion for each subject
min_idx = np.argmin(cv_gev_df.loc[cv_gev_df["Fit_Criteria"]=="CV Criterion","Value"])
cv_gev_df.loc[cv_gev_df["Fit_Criteria"]=="CV Criterion"].iloc[min_idx]
# Visualize the microstates
# Prior to re-ordering
for ii in range(len(n_clusters)):
plot_dualmicro(n_clusters[ii], microstate_summary_results[ii][0], microstate_summary_results[ii][1], epoch.info)
# =========================================================================
# # Re-order two-person microstates
# # This is only run once, after microstates are created
# # We only do it for 8 microstates, which was the optimal number
# =========================================================================
n_maps = 8
ii = n_clusters.index(n_maps)
with open(f"{microstate_save_path}Dualmicro_fit_all_{ff}_data_maps{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
maps, m_labels, gfp_peaks, gev, cv_min, pair_idx = microstate_results
plot_dualmicro(n_maps, maps, gev, epoch.info)
# Make dictionary with n_maps and new order
# All 4 top row consecutively, followed by 4 bot row
manual_reordering_template = {"8_alpha":[5,2,7,0,1,4,6,3],
"8_beta":[4,1,3,6,7,5,0,2],
"8_broadband":[6,3,4,0,2,1,5,7]}
new_order = manual_reordering_template[f"{n_maps}_{ff}"]
maps, gev, m_labels = reorder_microstate_results(new_order, maps, gev, m_labels)
# Plot again to check it worked
plot_dualmicro(n_maps, maps, gev, epoch.info)
# Since neuronal activity is often oscillating, this causes polarity inversions
# Microstates ignores the sign, and hence the polarity in the map is arbitrary
# It is only the relative difference within the plot that is interesting
# depending on initiation. We can thus freely change the sign for visualization
# For two-person microstates, each person's map is sign-changed separately
manual_sign_correction = {"8_alpha":[[-1,-1,-1,1,-1,1,1,1],[1,1,1,-1,-1,1,1,1]],
"8_beta":[[-1,-1,-1,-1,-1,1,1,1],[-1,-1,1,-1,-1,1,1,-1]],
"8_broadband":[[1,1,-1,-1,1,-1,-1,-1],[-1,1,-1,-1,1,-1,-1,-1]]}
sign_swap = manual_sign_correction[f"{n_maps}_{ff}"]
maps = sign_swap_microstates(sign_swap, maps, n_maps, n_channels)
# Plot a final time for last confirmation
plot_dualmicro(n_maps, maps, gev, epoch.info)
# Close all figures and repeat by changing n_maps
plt.close("all")
### Save reordered results
n_maps = 8
ii = n_clusters.index(n_maps)
with open(f"{microstate_save_path}Dualmicro_fit_all_{ff}_data_maps{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
maps, m_labels, gfp_peaks, gev, cv_min, pair_idx = microstate_results
# Re-order
new_order = manual_reordering_template[str(n_maps)]
maps, gev, m_labels = reorder_microstate_results(new_order, maps, gev, m_labels)
# Sign alignment
maps = sign_swap_microstates(sign_swap, maps, n_maps, n_channels)
# Overwrite variable
microstate_results = maps, m_labels, gfp_peaks, gev, cv_min, pair_idx
# Save to new file
with open(f"{microstate_save_path}Reordered/Dualmicro_fit_all_{ff}_data_maps{n_maps}.pkl", "wb") as filehandle:
pickle.dump(microstate_results, filehandle) # [maps, L, gfp_peaks, gev, cv_min, pair_idx]
# Save topomaps for the microstates
save_path = f"{fig_save_path}Microstates/Fit_all_{ff}/"
with open(f"{microstate_save_path}Reordered/Dualmicro_fit_all_{ff}_data_maps{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
maps, m_labels, gfp_peaks, gev, cv_min, pair_idx = microstate_results
fig = plot_dualmicro(n_maps, maps, gev, epoch.info)
fig.savefig(save_path+f"Dualmicro_fit_all_{ff}_maps{n_maps}"+".png")
# Save svg for Paper
fig.savefig(save_path+f"Dualmicro_fit_all_{ff}_maps{n_maps}"+".svg")
### Save svg with fixed color scales across all microstates
vlims = (np.min(maps), np.max(maps))
fig = plot_dualmicro(n_maps, maps, gev, vlims, epoch.info, vlims)
fig.savefig(save_path+f"Dualmicro_fit_all_{ff}_fixed_colorscale_maps{n_maps}"+".png")
fig.savefig(save_path+f"Dualmicro_fit_all_{ff}_fixed_colorscale_maps{n_maps}"+".svg")
# =========================================================================
# # Estimate two-person microstate metrics/features
# # There might be a small error introduced due to gaps in the time series from
# # dropped segments, e.g. when calculating the transition probability as
# # the time series is discontinuous due to the gaps. But with the high sampling rate
# # only a very small fraction of the samples have discontinuous neighbors
# =========================================================================
"""
Overview of common microstate features:
1. Average duration a given microstate remains stable (Dur)
2. Frequency occurrence, independent of individual duration (Occ)
Average number of times a microstate becomes dominant per second
3. Ratio of total Time Covered (TCo)
4. Transition probabilities (TMx)
5. Ratio of shannon entropy relative to theoretical max chaos (Ent)
"""
# Hard-coded the optimal number of microstates based on CV criterion and GEV
n_maps = 8
# Load all microstate results
with open(f"{microstate_save_path}Reordered/Dualmicro_fit_all_{ff}_data_maps{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
# Load all trialinfos
with open(f"{microstate_save_path}Dualmicro_fit_all_{ff}_trial_events_infos.pkl", "rb") as file:
trialinfo_list = pickle.load(file)
Microstate_names = [chr(ele) for ele in range(65,65+n_maps)]
m_labels = [0]*n_pairs
events = [0]*n_pairs
m_feats = [0]*n_pairs
for i in range(n_pairs):
m_labels[i], events[i], m_feats[i] = dualmicro_fit_all_feature_computation(i)
print(f"Finished computing microstate features for pair {Pair_id[i]}")
# Save the raw microstate features
with open(f"{microstate_save_path}/raw_dualmicro_fit_all_{ff}_features_maps{n_maps}.pkl", "wb") as filehandle:
pickle.dump(m_feats, filehandle) # [Subject][Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr] [Event, map*]
# * the feature is calculated for each map, where applicable.
# Transition matrix is calculated for each map -> map transition probability
# with open(f"{microstate_save_path}/raw_computed_dualmicro_fit_all_{ff}_features.pkl", "rb") as file:
# m_feats = pickle.load(file) # [Subject][Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr] [Event, map*]
### Convert all features to dataframes for further processing
col_names = ["Pair_ID", "Event_ID", "Microstate", "Value"]
col_values = [Pair_id,list(collapsed_event_id.keys()),Microstate_names]
dtypes = [int,str,str,"float64"]
# Mean duration
Dur_arr = np.stack([ele[0] for ele in m_feats]) # [Subject, event, n_map]
Dur_df = numpy_arr_to_pandas_df(Dur_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Duration"]*len(Dur_df)
Dur_df.insert(2, "Measurement", measurement_id)
# Save df
Dur_df.to_pickle(os.path.join(microstate_save_path,"Dualmicro_fit_all_{ff}_duration_df.pkl"))
# Frequency of occurrence per sec
Occ_arr = np.stack([ele[1] for ele in m_feats]) # [Subject, event, n_map]
Occ_df = numpy_arr_to_pandas_df(Occ_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Occurrence"]*len(Occ_df)
Occ_df.insert(2, "Measurement", measurement_id)
# Save df
Occ_df.to_pickle(os.path.join(microstate_save_path,"Dualmicro_fit_all_{ff}_occurrence_df.pkl"))
# Ratio total Time Covered
TCo_arr = np.stack([ele[2] for ele in m_feats]) # [Subject, event, n_map]
TCo_df = numpy_arr_to_pandas_df(TCo_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Time_covered"]*len(TCo_df)
TCo_df.insert(2, "Measurement", measurement_id)
# Save df
TCo_df.to_pickle(os.path.join(microstate_save_path,"Dualmicro_fit_all_{ff}_ratio_time_covered_df.pkl"))
# Transition matrix should be read as probability of row to column
xi, xj = np.meshgrid(Microstate_names,Microstate_names)
_, arrow = np.meshgrid(Microstate_names,["->"]*n_maps)
transition_info = np.char.add(np.char.add(xj,arrow),xi)
TMx_arr = np.stack([ele[3] for ele in m_feats]) # [Subject, event, n_map, n_map]
TMx_arr = TMx_arr.reshape((n_pairs,len(collapsed_event_id),n_maps*n_maps)) # Flatten the maps to 1D
col_names = ["Pair_ID", "Event_ID", "Transition", "Value"]
col_values = [Pair_id,list(collapsed_event_id.keys()),transition_info.flatten()]
TMx_df = numpy_arr_to_pandas_df(TMx_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Probability"]*len(TMx_df)
TMx_df.insert(2, "Measurement", measurement_id)
# Save df
TMx_df.to_pickle(os.path.join(microstate_save_path,"Dualmicro_fit_all_{ff}_transition_df.pkl"))
# Entropy
Ent_arr = np.stack([ele[4] for ele in m_feats]) # [Subject, event]
col_names = ["Pair_ID", "Event_ID", "Value"]
col_values = [Pair_id,list(collapsed_event_id.keys())]
dtypes = [int, str, "float64"]
Ent_df = numpy_arr_to_pandas_df(Ent_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Entropy"]*len(Ent_df)
Ent_df.insert(2, "Measurement", measurement_id)
# Save df
Ent_df.to_pickle(os.path.join(microstate_save_path,"Dualmicro_fit_all_{ff}_ratio_entropy_df.pkl"))
# %% Backfit two-person microstates to pseudo-pairs
# The pseudo-pairs are created for all participants except the real pair.
# This is fine for symmetrical tasks, e.g. rest and coupled.
# But not for assymmetrical tasks like observation and leader.
# We might have a leader - leader pseudo-pair.
# Hence we only look at ppn1 with ppn2 from different pairs and exclude
# ppn1 with ppn1 or ppn2 with ppn2
for f in len(all_freq_ranges):
ff = freq_names[f]
freq_range0 = all_freq_ranges[f]
# =========================================================================
# It might be an advantage to run the backfitting of microstates on a HPC
# =========================================================================
# To save time and prevent reloading the same EEG over and over, I divided
# the prepare array function into a load and combine function
# By loading all into memory, I can skip loading for every combination
# but this requires a very high memory, which is fortunately not a problem on the hpc
# I am limiting the pseudo-pairs to be where ppn1 ends with 1 and ppn2 with 2
# Which means we have 21 * 20 options
n_pseudo_pairs = n_pairs*(n_pairs-1)
# To not load data 420 times for two participants, we preload all EEG data to ram
c_time1 = time_now(); print("Starting load",c_time1)
all_micro_data = [0]*n_subjects
all_trial_data = [0]*n_subjects
for i in range(n_subjects):
all_micro_data[i], all_trial_data[i] = load_microstate_arrays(i)
print("Load finished", time_now())
# Get the prototypical alpha maps
n_maps = 8
with open(f"{microstate_save_path}Reordered/Dualmicro_fit_all_{ff}_data_maps{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
prototype_map = microstate_results[0]
# Start the backfitting
m_labels = [0]*n_pseudo_pairs
events = [0]*n_pseudo_pairs
GEVs = [0]*n_pseudo_pairs
counter = 0
pseudo_pair_id = []
for i in range(n_subjects):
for j in range(n_subjects):
# Skip if the subject is the same
if np.abs(Subject_id[i]-Subject_id[j]) == 0:
continue
# Skip if the subject are from the same pair
if np.abs(Subject_id[i]-Subject_id[j]) == 1:
continue
# Skip if ppn1 is not ending on 1, and ppn2 not ending on 2
if not (str(Subject_id[i])[-1] == "1") & (str(Subject_id[j])[-1] == "2"):
continue
# A valid pseudo pair
else:
# Get the synchronized events
event0 = get_synch_events_from_pseudo_pairs(all_trial_data[i],all_trial_data[j])
# Get the preloaded micro data
micro_data1 = all_micro_data[i]
micro_data2 = all_micro_data[j]
# Get the synchronized and concatenated micro data in alpha
micro_data0 = combine_two_person_microstate_arrays(micro_data1, micro_data2, event0, sfreq, freq_range=freq_range0)
# Backfit and get the labels
L, GEV = pseudo_pair_dualmicro_backfitting(micro_data0, prototype_map, event0, n_maps, sfreq)
# Save the results
m_labels[counter], GEVs[counter], events[counter] = L, GEV, event0
pseudo_pair_id.append(f"{Subject_id[i]}-{Subject_id[j]}")
# Move counter
counter += 1
print(f"Finished backfitting for pseudo pair {pseudo_pair_id[-1]}")
print("Started", c_time1, "\nCurrent",time_now())
backfit_results = [pseudo_pair_id, m_labels, GEVs, events]
# Save the results from all pseudo pairs
with open(f"{microstate_save_path}Reordered/Backfitting/Dualmicro_fit_all_{ff}_data_maps{n_maps}.pkl", "wb") as filehandle:
pickle.dump(backfit_results, filehandle) # [pseudo_pair_id, L, GEV, events]
# =========================================================================
# Estimate two-person microstate metrics/features
# There might be a small error introduced due to gaps in the time series from
# dropped segments, e.g. when calculating the transition probability as
# the time series is discontinuous due to the gaps. But with the high sampling rate
# only a very small fraction of the samples have discontinuous neighbors
# =========================================================================
"""
Overview of common microstate features:
1. Average duration a given microstate remains stable (Dur)
2. Frequency occurrence, independent of individual duration (Occ)
Average number of times a microstate becomes dominant per second
3. Ratio of total Time Covered (TCo)
4. Transition probabilities (TMx)
5. Ratio of shannon entropy relative to theoretical max chaos (Ent)
"""
n_maps = 8
# Load all the backfit pseudo-pair results
with open(f"{microstate_save_path}Reordered/Backfitting/Dualmicro_fit_all_{ff}_data_maps{n_maps}.pkl", "rb") as file:
backfit_results = pickle.load(file) # [pseudo_pair_id, L, GEV, events]
# Hard-coded the optimal number of microstates based on CV criterion and GEV
n_maps = 8
Microstate_names = [chr(ele) for ele in range(65,65+n_maps)]
pseudo_pair_id = backfit_results[0]
n_pseudo_pairs = len(pseudo_pair_id)
m_labels = [0]*n_pseudo_pairs
events = [0]*n_pseudo_pairs
m_feats = [0]*n_pseudo_pairs
for i in range(n_pseudo_pairs):
m_labels[i], events[i], m_feats[i] = dualmicro_fit_all_pseudo_pair_feature_computation(i,\
n_maps, backfit_results, sfreq, event_id, collapsed_event_id)
print(f"Finished computing microstate features for psuedo pair {pseudo_pair_id[i]}")
# Save the raw microstate features
with open(f"{microstate_save_path}/raw_dualmicro_fit_all_{ff}_pseudo_pairs_features_maps{n_maps}.pkl", "wb") as filehandle:
pickle.dump(m_feats, filehandle) # [Subject][Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr] [Event, map*]
# * the feature is calculated for each map, where applicable.
# Transition matrix is calculated for each map -> map transition probability
# with open(f"{microstate_save_path}/raw_computed_dualmicro_fit_all_{ff}_features.pkl", "rb") as file:
# m_feats = pickle.load(file) # [Subject][Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr] [Event, map*]
### Convert all features to dataframes for further processing
col_names = ["Pseudo_Pair_ID", "Event_ID", "Microstate", "Value"]
col_values = [pseudo_pair_id,list(collapsed_event_id.keys()),Microstate_names]
dtypes = [str,str,str,"float64"]
# Mean duration
Dur_arr = np.stack([ele[0] for ele in m_feats]) # [Subject, event, n_map]
Dur_df = numpy_arr_to_pandas_df(Dur_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Duration"]*len(Dur_df)
Dur_df.insert(2, "Measurement", measurement_id)
# Save df
Dur_df.to_pickle(os.path.join(microstate_save_path,"Dualmicro_fit_all_{ff}_pseudo_pairs_duration_df.pkl"))
# Frequency of occurrence per sec
Occ_arr = np.stack([ele[1] for ele in m_feats]) # [Subject, event, n_map]
Occ_df = numpy_arr_to_pandas_df(Occ_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Occurrence"]*len(Occ_df)
Occ_df.insert(2, "Measurement", measurement_id)
# Save df
Occ_df.to_pickle(os.path.join(microstate_save_path,"Dualmicro_fit_all_{ff}_pseudo_pairs_occurrence_df.pkl"))
# Ratio total Time Covered
TCo_arr = np.stack([ele[2] for ele in m_feats]) # [Subject, event, n_map]
TCo_df = numpy_arr_to_pandas_df(TCo_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Time_covered"]*len(TCo_df)
TCo_df.insert(2, "Measurement", measurement_id)
# Save df
TCo_df.to_pickle(os.path.join(microstate_save_path,"Dualmicro_fit_all_{ff}_pseudo_pairs_ratio_time_covered_df.pkl"))
# Transition matrix should be read as probability of row to column
xi, xj = np.meshgrid(Microstate_names,Microstate_names)
_, arrow = np.meshgrid(Microstate_names,["->"]*n_maps)
transition_info = np.char.add(np.char.add(xj,arrow),xi)
TMx_arr = np.stack([ele[3] for ele in m_feats]) # [Subject, event, n_map, n_map]
TMx_arr = TMx_arr.reshape((n_pseudo_pairs,len(collapsed_event_id),n_maps*n_maps)) # Flatten the maps to 1D
col_names = ["Pseudo_Pair_ID", "Event_ID", "Transition", "Value"]
col_values = [pseudo_pair_id,list(collapsed_event_id.keys()),transition_info.flatten()]
TMx_df = numpy_arr_to_pandas_df(TMx_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Probability"]*len(TMx_df)
TMx_df.insert(2, "Measurement", measurement_id)
# Save df
TMx_df.to_pickle(os.path.join(microstate_save_path,"Dualmicro_fit_all_{ff}_pseudo_pairs_transition_df.pkl"))
# Entropy
Ent_arr = np.stack([ele[4] for ele in m_feats]) # [Subject, event]
col_names = ["Pseudo_Pair_ID", "Event_ID", "Value"]
col_values = [pseudo_pair_id,list(collapsed_event_id.keys())]
dtypes = [str,str,"float64"]
Ent_df = numpy_arr_to_pandas_df(Ent_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Entropy"]*len(Ent_df)
Ent_df.insert(2, "Measurement", measurement_id)
# Save df
Ent_df.to_pickle(os.path.join(microstate_save_path,"Dualmicro_fit_all_{ff}_pseudo_pairs_ratio_entropy_df.pkl"))
# %% eLORETA on Intrabrain microstates
### Make forward solutions
# Computed using the fsaverage template MRI
# # First time setup will need to download fsaverage templates
# mne.datasets.fetch_fsaverage()
fs_dir = "C:/Users/glia/mne_data/MNE-fsaverage-data/fsaverage"
subjects_dir = os.path.dirname(fs_dir)
trans = "fsaverage"
src = os.path.join(fs_dir, "bem", "fsaverage-ico-5-src.fif")
bem = os.path.join(fs_dir, "bem", "fsaverage-5120-5120-5120-bem-sol.fif")
# Read the template sourcespace
sourcespace = mne.read_source_spaces(src)
# Since I use a template, I only need to make the forward operator once
# As we assume the channel positions are fixed approximately the same
# for all subjects using the same caps
subject_eeg = epoch.copy()
subject_eeg.set_eeg_reference(projection=True) # needed for inverse modelling
# Make forward solution
fwd = mne.make_forward_solution(subject_eeg.info, trans=trans, src=src,
bem=bem, eeg=True, mindist=5.0, n_jobs=1)
# # Save forward operator
# fname_fwd = "./Source_fwd/fsaverage_{}-fwd.fif".format(study_order[i])
# mne.write_forward_solution(fname_fwd, fwd, overwrite=True)
# # Check the alignment looks correct between EEG sensors and the template
# mne.viz.plot_alignment(
# subject_eeg.info, trans, src=src, fwd=fwd, dig=True,
# meg=["helmet", "sensors"], subjects_dir=subjects_dir, surfaces="auto")
### Load Parcellation
# Desikan-Killiany atlas (34 ROI from both hemispheres = 68 ROIs)
# Named aparc.annot in MNE python fsaverage folder
labels = mne.read_labels_from_annot("fsaverage", parc="aparc",
subjects_dir=subjects_dir)
labels = labels[:-1] # remove unknowns
label_names = [label.name for label in labels]
n_roi = len(labels)
# Prepare brain lobe information
Frontal_rois = ['superiorfrontal-lh','superiorfrontal-rh',
'rostralmiddlefrontal-lh','rostralmiddlefrontal-rh',
'caudalmiddlefrontal-lh','caudalmiddlefrontal-rh',
'parsopercularis-lh','parsopercularis-rh',
'parstriangularis-lh','parstriangularis-rh',
'parsorbitalis-lh','parsorbitalis-rh',
'lateralorbitofrontal-lh','lateralorbitofrontal-rh',
'medialorbitofrontal-lh','medialorbitofrontal-rh',
'precentral-lh','precentral-rh',
'paracentral-lh','paracentral-rh',
'frontalpole-lh','frontalpole-rh']
Parietal_rois = ['superiorparietal-lh','superiorparietal-rh',
'inferiorparietal-lh','inferiorparietal-rh',
'supramarginal-lh','supramarginal-rh',
'postcentral-lh','postcentral-rh',
'precuneus-lh','precuneus-rh']
Temporal_rois = ['superiortemporal-lh','superiortemporal-rh',
'middletemporal-lh','middletemporal-rh',
'inferiortemporal-lh','inferiortemporal-rh',
'bankssts-lh','bankssts-rh',
'fusiform-lh','fusiform-rh',
'transversetemporal-lh','transversetemporal-rh',
'entorhinal-lh','entorhinal-rh',
'temporalpole-lh','temporalpole-rh',
'parahippocampal-lh','parahippocampal-rh']
Occipital_rois = ['lateraloccipital-lh','lateraloccipital-rh',
'lingual-lh','lingual-rh',
'cuneus-lh','cuneus-rh',
'pericalcarine-lh','pericalcarine-rh']
Cingulate_rois = ['rostralanteriorcingulate-lh','rostralanteriorcingulate-rh',
'caudalanteriorcingulate-lh','caudalanteriorcingulate-rh',
'posteriorcingulate-lh','posteriorcingulate-rh',
'isthmuscingulate-lh','isthmuscingulate-rh']
Insular_rois = ['insula-lh','insula-rh']
Lobes = [Frontal_rois,Parietal_rois,Temporal_rois,Occipital_rois,Cingulate_rois,Insular_rois]
Brain_region_labels = ["Frontal","Parietal","Temporal","Occipital","Cingulate","Insular"]
Brain_region_hemi_labels = np.repeat(Brain_region_labels,2).astype("<U12")
Brain_region_hemi_labels[::2] = [ele+"-lh" for ele in Brain_region_labels]
Brain_region_hemi_labels[1::2] = [ele+"-rh" for ele in Brain_region_labels]
Brain_region = np.array(label_names, dtype = "<U32")
for l in range(len(Lobes)):
Brain_region[np.array([i in Lobes[l] for i in Brain_region])] = Brain_region_labels[l]
### Concatenate the microstates into one Raw Object to apply inverse on it
n_maps = 8
Microstate_names = [chr(ele) for ele in range(65,65+n_maps)]
for f in len(all_freq_ranges):
ff = freq_names[f]
with open(f"{microstate_save_path}Reordered/Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
# Get the microstates and reshape to have channels in the first dim
maps = microstate_results[0]
maps = maps.transpose()
raw_maps = mne.io.RawArray(maps,subject_eeg.info)
raw_maps._filenames = [""] # Fix error with NoneType for "filename" for raw created with RawArray
raw_maps.set_eeg_reference(projection=True) # needed for inverse modelling
# Using assumption about equal variance and no correlations I make a diagonal matrix as cov
noise_cov = mne.make_ad_hoc_cov(subject_eeg.info, None)
# Make inverse operator
# Using default depth parameter = 0.8 and free orientation (loose = 1)
inverse_operator = mne.minimum_norm.make_inverse_operator(subject_eeg.info,
fwd, noise_cov,
loose = 1, depth = 0.8,
verbose = 0)
src_inv = inverse_operator["src"]
# Compute inverse solution and retrieve the source localized microstate activities for each label
# Define regularization
snr = 3 # Default setting
# Use eLORETA and only keep the activity normal to the cortical surface
stc = mne.minimum_norm.apply_inverse_raw(raw_maps,inverse_operator,
lambda2 = 1/(snr**2),
pick_ori = "normal",
method = "eLORETA", verbose = 2)
# Get the source activity in the ROIs
label_activity = mne.extract_label_time_course(stc, labels, src_inv, mode="mean_flip",
return_generator=False, verbose=0)
# Visualize the microstates in source space
# This way of plotting makes the color scale fixed across microstates
brain = stc.plot(
hemi="lh",
subjects_dir=subjects_dir,
smoothing_steps=1,
)
### Convert Label Activity to Pandas DataFrame
# With ROI names and then add Brain Region label
col_names = ["ROI", "Microstate", "Value"]
col_names = ["Microstate", "ROI", "Value"]
col_val = [Microstate_names, label_names]
# Create the source microstate activity dataframe
sMicro_df = numpy_arr_to_pandas_df(label_activity.T, col_names = col_names, col_values = col_val)
assert sMicro_df.loc[(sMicro_df["ROI"]==label_names[4])&
(sMicro_df["Microstate"]==Microstate_names[3]),
"Value"].iloc[0] == label_activity[4,3]
# Add brain region information
sMicro_df.insert(2, "Brain_region", np.tile(Brain_region,int(sMicro_df.shape[0]/n_roi)))
sMicro_df["Brain_region"] = sMicro_df["Brain_region"].astype("category").\
cat.reorder_categories(Brain_region_labels, ordered=True)
# Add hemisphere information
sMicro_df.insert(3, "Hemisphere", [ele[-2:] for ele in sMicro_df["ROI"]])
# Add a colum that combines brain region and hemisphere for plotting
sMicro_df.insert(4, "Brain_region_hemi", [b+"-"+h for b, h in zip(sMicro_df["Brain_region"],sMicro_df["Hemisphere"])])
sMicro_df["Brain_region_hemi"] = sMicro_df["Brain_region_hemi"].astype("category").\
cat.reorder_categories(Brain_region_hemi_labels, ordered=True)
# Save the dataframe
sMicro_df.to_pickle(os.path.join(microstate_save_path,f"Single_micro_{ff}_source_activity_df.pkl"))
# %% eLORETA on two-brain microstates
# Continued based on fwd operator and template loaded for intrabrain
n_maps = 8
Microstate_names = [chr(ele) for ele in range(65,65+n_maps)]
for f in len(all_freq_ranges):
ff = freq_names[f]
with open(f"{microstate_save_path}Reordered/Dualmicro_fit_all_{ff}_data_maps{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
# Get the microstates
maps = microstate_results[0]
maps = maps.reshape(2*n_maps,n_channels)
# # Check the maps were split properly
# plot_microstates(n_maps, maps[:8], microstate_results[3])
# plot_microstates(n_maps, maps[8:], microstate_results[3])
# Maps are ordered as: ppn1 A, ppn2 A, ppn1 B, ppn2 B etc
# Transpose to have channels in the first dim
maps = maps.transpose()
raw_maps = mne.io.RawArray(maps,subject_eeg.info)
raw_maps._filenames = [""] # Fix error with NoneType for "filename" for raw created with RawArray
raw_maps.set_eeg_reference(projection=True) # needed for inverse modelling
# Using assumption about equal variance and no correlations I make a diagonal matrix as cov
noise_cov = mne.make_ad_hoc_cov(subject_eeg.info, None)
# Make inverse operator
# Using default depth parameter = 0.8 and free orientation (loose = 1)
inverse_operator = mne.minimum_norm.make_inverse_operator(subject_eeg.info,
fwd, noise_cov,
loose = 1, depth = 0.8,
verbose = 0)
src_inv = inverse_operator["src"]
# Compute inverse solution and retrieve the source localized microstate activities for each label
# Define regularization
snr = 3 # Default setting
# Use eLORETA and only keep the activity normal to the cortical surface
stc = mne.minimum_norm.apply_inverse_raw(raw_maps,inverse_operator,
lambda2 = 1/(snr**2),
pick_ori = "normal",
method = "eLORETA", verbose = 2)
# Get the source activity in the ROIs
label_activity = mne.extract_label_time_course(stc, labels, src_inv, mode="mean_flip",
return_generator=False, verbose=0)
# Visualize the microstates in source space
# This way of plotting makes the color scale fixed across microstates
brain = stc.plot(
hemi="lh",
subjects_dir=subjects_dir,
smoothing_steps=1,
)
# Visualize with different color scales for each microstate
Microstate_names2 = np.repeat(Microstate_names,2).astype("<U2")
Microstate_names2[::2] = [ele+"1" for ele in Microstate_names]
Microstate_names2[1::2] = [ele+"2" for ele in Microstate_names]
# Save source activations for each microstate
# Lateral and medial for each hemisphere + dorsal + flatmaps
save_path = f"{fig_save_path}Microstates/SourceDualmicroPrototypes/"
hemis = ["lh","rh"]
views = ["lateral","medial"]
for i in range(len(Microstate_names2)):
times0 = np.linspace(0,1,sfreq+1)[:2*n_maps+1]
stc0 = stc.copy().crop(times0[i],times0[i+1],include_tmax=False)
# Color bar limits defined as max saturation of top 1% (yellow or teal)
# middle at 5%, which means they will have alpha = 1 and progressively be
# closer to yellow or teal
# Lower boundary at 10%, which means they will be red/blue but with decreased
# transparency
clim_max = -(np.sort(-np.abs(stc0.data),axis=0)[stc0.shape[0]//100])[0]
clim_mid = -(np.sort(-np.abs(stc0.data),axis=0)[stc0.shape[0]//20])[0]
clim_min = -(np.sort(-np.abs(stc0.data),axis=0)[stc0.shape[0]//10])[0]
clim0 = {"kind":"value","pos_lims":[clim_min,clim_mid,clim_max]}
# Lateral and medial
for h in range(len(hemis)):
hh = hemis[h]
brain = stc0.plot(
hemi=hh,
subjects_dir=subjects_dir,
smoothing_steps=10, # spatial smoothing
colorbar=False,
background="white",
cortex="classic",
size=800,
transparent=True,
views=views[0],
clim=clim0,
)
brain.save_image(os.path.join(save_path, f"Dualmicro_source_{Microstate_names2[i]}_{hh}_{views[0]}"+".png"))
brain.show_view(views[1])
brain.save_image(os.path.join(save_path, f"Dualmicro_source_{Microstate_names2[i]}_{hh}_{views[1]}"+".png"))
# Dorsal map
brain = stc0.plot(
hemi="both",
subjects_dir=subjects_dir,
smoothing_steps=10, # spatial smoothing
colorbar=True,
background="white",
cortex="classic",
size=1500,
transparent=True,
views="dorsal",
clim=clim0,
)
brain.save_image(os.path.join(save_path, f"Dualmicro_source_{Microstate_names2[i]}_dorsal"+".png"))
# Flat map
brain = stc0.plot(
hemi="both",
surface="flat",
subjects_dir=subjects_dir,
smoothing_steps=10, # spatial smoothing
colorbar=False,
background="white",
cortex="classic",
size=1500,
transparent=True,
views="flat",
clim=clim0,
)
brain.save_image(os.path.join(save_path, f"Dualmicro_source_{Microstate_names2[i]}_flat"+".png"))
# Close all figures
mne.viz.close_all_3d_figures()
# Mean
brain = stc.mean().plot(
hemi="lh",
subjects_dir=subjects_dir,
smoothing_steps=10,
)
### Convert Label Activity to Pandas DataFrame
# With ROI names and then add Brain Region label
col_names = ["ROI", "Microstate", "Value"]
col_names = ["Microstate", "ROI", "Value"]
col_val = [Microstate_names2, label_names]
dtypes = [str, str, "float64"]
# Create the source microstate activity dataframe
sMicro_df = numpy_arr_to_pandas_df(label_activity.T, col_names, col_val, dtypes)
assert sMicro_df.loc[(sMicro_df["ROI"]==label_names[4])&
(sMicro_df["Microstate"]==Microstate_names2[3]),
"Value"].iloc[0] == label_activity[4,3]
# Add brain region information
sMicro_df.insert(2, "Brain_region", np.tile(Brain_region,int(sMicro_df.shape[0]/n_roi)))
sMicro_df["Brain_region"] = sMicro_df["Brain_region"].astype("category").\
cat.reorder_categories(Brain_region_labels, ordered=True)
# Add hemisphere information
sMicro_df.insert(3, "Hemisphere", [ele[-2:] for ele in sMicro_df["ROI"]])
# Add a colum that combines brain region and hemisphere for plotting
sMicro_df.insert(4, "Brain_region_hemi", [b+"-"+h for b, h in zip(sMicro_df["Brain_region"],sMicro_df["Hemisphere"])])
sMicro_df["Brain_region_hemi"] = sMicro_df["Brain_region_hemi"].astype("category").\
cat.reorder_categories(Brain_region_hemi_labels, ordered=True)
# Save the dataframe
sMicro_df.to_pickle(os.path.join(microstate_save_path,"Dualmicro_{ff}_source_activity_df.pkl"))
# %% LRTC with DFA on Two-person microstate label time series
# Using Detrended Fluctuation Analysis (DFA)
# Adapted from Python Implementation by Arthur-Ervin Avramiea <a.e.avramiea@vu.nl>
# From NBT2 toolbox
"""
See Hardstone et al, 2012 for more info
Perform DFA
1 Compute cumulative sum of time series to create signal profile
2 Define set of window sizes (see below)
3 Remove the linear trend using least-squares for each window
4 Calculate standard deviation for each window and take the mean
5 Plot fluctuation function (Standard deviation) as function
for all window sizes, on double logarithmic scale
6 The DFA exponent alpha correspond to Hurst exponent
f(L) = sd = L^alpha (with alpha as linear coefficient in log plot)
If 0 < alpha < 0.5: The process exhibits anti-correlations
If 0.5 < alpha < 1: The process exhibits positive correlations
If alpha = 0.5: The process is indistinguishable from a random process
If 1.0 < alpha < 2.0: The process is non-stationary. H = alpha - 1
Window sizes should be equally spaced on a logarithmic scale
Sizes should be at least 4 samples and up to 10% of total signal length
### Specific for our microstate DFA analysis
We have 8 microstates, but to compute the random walk we will partition
the microstate sequence into two classes (see reference on microstate Hurst
https://pubmed.ncbi.nlm.nih.gov/20921381/)
A/B/C/D will be assigned the positive direction, while E/F/G/H will be
assigned the negative direction, corresponding to whether ppn1 or ppn2
are in one of the canonical microstates, while the other have a non-specific
(average) topography.
Each 25s trial is too short to estimate LRTC on, so I will concatenate all
the trials corresponding to each condition.
This should yield up to 25s * 16 trials = 400s of data for each condition,
except rest which is up to 120s * 2 trials = 240s
DFA is computed from 8 trials and then averaged, to avoid the
problem of flipping in the asymmetric trials. We change windows size to 5-20s
To ensure consistency the same procedure is applied to the symmetric trials
"""
# Window sizes
compute_interval = [5,20] # the window sizes should be between 5s and 30s
# Compute DFA window sizes for the given Interval
window_sizes = np.floor(np.logspace(-1,3,40) * sfreq).astype(int) # %logspace from 0.1 seccond (10^-1) to 1000 (10^3) seconds
window_sizes = window_sizes[(window_sizes >= compute_interval[0]*sfreq) & \
(window_sizes <= compute_interval[1]*sfreq)]
for f in len(all_freq_ranges):
ff = freq_names[f]
# Nolds are already using all cores so multiprocessing with concurrent makes it slower
n_maps = 8
with open(f"{microstate_save_path}Reordered/Dualmicro_fit_all_{ff}_data_maps{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
# Load all trialinfos
with open(f"{microstate_save_path}Dualmicro_fit_all_{ff}_trial_events_infos.pkl", "rb") as file:
trialinfo_list = pickle.load(file)
# Pre-allocate memory
DFA_arr = np.zeros((n_pairs,len(collapsed_event_id)))
Fluctuation_arr = np.zeros((n_pairs,len(collapsed_event_id),len(window_sizes)))
# Get start time
c_time1 = time.localtime()
c_time1 = time.strftime("%a %d %b %Y %H:%M:%S", c_time1)
print("Started {}".format(c_time1))
# Nolds are already using all cores so concurrent futures with make it slower
for i in range(n_pairs):
# Compute DFA
dfa_temp, fluc_temp = compute_dualmicro_DFA(i, microstate_results,
trialinfo_list, sfreq, window_sizes, event_id, collapsed_event_id, True)
# Save to array
DFA_arr[i] = dfa_temp
Fluctuation_arr[i] = fluc_temp
print("Finished {} out of {} pairs".format(i+1,n_pairs))
# Get ending time
c_time2 = time.localtime()
c_time2 = time.strftime("%a %d %b %Y %H:%M:%S", c_time2)
print(("Started {} \nEnded Time {}".format(c_time1,c_time2)))
# Save the raw DFA analysis data
np.save(microstate_save_path+"DFA_arr.npy", DFA_arr)
np.save(microstate_save_path+"Fluctuation_arr.npy", Fluctuation_arr)
# Convert to Pandas dataframe (DFA exponent)
col_names = ["Pair_ID", "Event_ID", "Value"]
col_values = [Pair_id,list(collapsed_event_id.keys())]
dtypes = ["int64",str,"float64"]
DFA_df = numpy_arr_to_pandas_df(DFA_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["DFA"]*len(DFA_df)
DFA_df.insert(2, "Measurement", measurement_id)
# Save df
DFA_df.to_pickle(os.path.join(microstate_save_path,f"Dualmicro_{ff}_DFA_exponent_df.pkl"))
# %% DFA in pseudo-pairs
for f in len(all_freq_ranges):
ff = freq_names[f]
# Nolds are already using all cores so multiprocessing with concurrent makes it slower
n_maps = 8
# Load all the backfit pseudo-pair results
with open(f"{microstate_save_path}Reordered/Backfitting/Dualmicro_fit_all_{ff}_data_maps{n_maps}.pkl", "rb") as file:
backfit_results = pickle.load(file) # [pseudo_pair_id, L, GEV, events]
# Pre-allocate memory
DFA_arr = np.zeros((n_pairs,len(collapsed_event_id)))
Fluctuation_arr = np.zeros((n_pairs,len(collapsed_event_id),len(window_sizes)))
# Get start time
c_time1 = time.localtime()
c_time1 = time.strftime("%a %d %b %Y %H:%M:%S", c_time1)
print("Started {}".format(c_time1))
# Nolds are already using all cores so concurrent futures with make it slower
for i in range(n_pairs):
# Compute DFA
dfa_temp, fluc_temp = compute_dualmicro_DFA_pseudo(i, backfit_results,
sfreq, window_sizes, event_id, collapsed_event_id, True)
# Save to array
DFA_arr[i] = dfa_temp
Fluctuation_arr[i] = fluc_temp
print("Finished {} out of {} pairs".format(i+1,n_pairs))
# Get ending time
c_time2 = time.localtime()
c_time2 = time.strftime("%a %d %b %Y %H:%M:%S", c_time2)
print(("Started {} \nEnded Time {}".format(c_time1,c_time2)))
# Save the raw DFA analysis data
np.save(microstate_save_path+"DFA_arr.npy", DFA_arr)
np.save(microstate_save_path+"Fluctuation_arr.npy", Fluctuation_arr)
# Convert to Pandas dataframe (DFA exponent)
col_names = ["Pair_ID", "Event_ID", "Value"]
col_values = [Pair_id,list(collapsed_event_id.keys())]
dtypes = ["int64",str,"float64"]
DFA_df = numpy_arr_to_pandas_df(DFA_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["DFA"]*len(DFA_df)
DFA_df.insert(2, "Measurement", measurement_id)
# Save df
DFA_df.to_pickle(os.path.join(microstate_save_path,f"Dualmicro_{ff}_DFA_exponent_df.pkl"))
# %% Time-lagged inter-brain microstate synchrony
# Hard-coded the optimal number of microstates based on CV criterion and GEV
n_maps = 5
# The lag (number of samples) we will iterate over to find greatest time-lagged interbrain microstate synchrony
lag_search_range = sfreq # 1 second in both directions
lag_interval = np.linspace(-lag_search_range,lag_search_range,lag_search_range*2+1).astype(int)
Microstate_names = [chr(ele) for ele in range(65,65+n_maps)]
# Insert Z as the symbol for non common microstate
Microstate_names.insert(0,"Z")
# Loop over frequencies
for f in len(all_freq_ranges):
ff = freq_names[f]
# Load all microstate results
with open(f"{microstate_save_path}Reordered/Intrabrain_microstate_fit_all_{ff}{n_maps}.pkl", "rb") as file:
microstate_results = pickle.load(file)
# Load all trialinfos
with open(f"{microstate_save_path}Intrabrain_microstate_fit_all_{ff}_trialinfos.pkl", "rb") as file:
trialinfo_list = pickle.load(file)
m_labels = [0]*(n_subjects//2)
events = [0]*(n_subjects//2)
m_feats = [0]*(n_subjects//2)
shift_info = [0]*(n_subjects//2)
Pair_id = [0]*(n_subjects//2)
for i in tqdm(range(n_subjects//2)):
m_labels[i], events[i], m_feats[i], shift_info[i] = shifted_interbrain_microstate_feature_computation(i,
n_maps, microstate_results, trialinfo_list, sfreq,
event_id, collapsed_event_id, lag_search_range, lag_interval)
Pair_id[i] = int(str(Subject_id[2*i])[1:-1])
print(f"Finished computing interbrain microstate features for pair {Pair_id[i]}")
Pair_id = [ele+100 for ele in Pair_id]
# Save the raw microstate features
with open(f"{microstate_save_path}/raw_shifted_interbrain_single_micro_fit_all_{ff}_maps{n_maps}.pkl", "wb") as filehandle:
pickle.dump([Pair_id, m_feats, shift_info], filehandle) # [Subject][Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr][Event, map*]
# * the feature is calculated for each map, where applicable.
# Transition matrix is calculated for each map -> map transition probability
# The first row and column correspond to the non common microstate, i.e.
# there is a different microstate in the pair
# with open(f"{microstate_save_path}/raw_shifted_interbrain_single_micro_fit_all_{ff}_maps{n_maps}.pkl", "rb") as file:
# Pair_id, m_feats, shift_info = pickle.load(file) # [Subject][Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr] [Event, map*]
n_pairs = len(Pair_id)
### Convert all features to dataframes for further processing
col_names = ["Pair_ID", "Event_ID", "Microstate", "Value"]
col_values = [Pair_id,list(collapsed_event_id.keys()),Microstate_names]
dtypes = [int,str,str,"float64"]
# Ratio total Time Covered
TCo_arr = np.stack([ele[2] for ele in m_feats]) # [Subject, event, n_map]
TCo_df = numpy_arr_to_pandas_df(TCo_arr, col_names, col_values, dtypes)
# Add dummy variable to enabling combining of dataframes
measurement_id = ["Time_covered"]*len(TCo_df)
TCo_df.insert(2, "Measurement", measurement_id)
# Save df
TCo_df.to_pickle(os.path.join(microstate_save_path,f"Shifted_IB_Single_micro_fit_all_{ff}_maps{n_maps}_ratio_time_covered_df.pkl"))