Brainwaves under medication: revealing class-specific neural signatures of psychotropic medication from 24,000 EEGs.

Code ↔ Paper

8 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 8 matches

[1] § Methods › Medicine groups ↔ Brainwaves_Under_Medication_Visualisation_and_PCA_analysis_code_clean.ipynb, lines 172–219 · score 0.83 · sedative hypnotic, NaSSA, AChE, Anticholinergic, opioid, SNRIs
[2] § Methods › Preprocessing ↔ preprocessing_REST_ASR_flexible_commented.m, lines 64–200 · score 0.80 · high pass filter, window criterion, channel interpolation, ASR, pipeline, FASTER
[3] § Methods › Preprocessing ↔ preprocessing_REST_ASR_flexible_commented.m, lines 1–62 · score 0.71 · CleanLine, infinity, kurtosis, EEGLAB, plugin, noise
[4] § Results › Data ↔ Brainwaves_Under_Medication_Visualisation_and_PCA_analysis_code_clean.ipynb, lines 172–219 · score 0.67 · AED Ca, AED Na, NaSSA, SARI, AP, atypical
[5] § Methods › Dimensionality reduction ↔ Brainwaves_Under_Medication_code_commented.m, lines 253–282 · score 0.62 · confidence interval, Principal Component, uncorrelated, coefficients, variance, PCA
[6] § Methods › EEG signal features ↔ Brainwaves_Under_Medication_code_commented.m, lines 37–58 · score 0.62 · feature exceeded, standard deviations, VAR, outliers, scored, patients
[7] § Methods › Medicine groups ↔ balance_groups_meds_DN.m, lines 1–55 · score 0.60 · Chi squared, classified, medication classes, binary, psychotropic, diagnosis
[8] § Results › Data ↔ Brainwaves_Under_Medication_Visualisation_and_PCA_analysis_code_clean.ipynb, lines 938–1070 · score 0.56 · hierarchical regression model, Holm, mixed, dimensionality, PCA, matched

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Jupyter notebook · 1,338 lines · 52 KB · no license · 3 matches

# %% [markdown]
# # Brainwaves Under Medication: visualisation and PCA analysis
#
# This notebook reads the outputs generated by the accompanying MATLAB analysis script and produces:
# - summary tables for significant PCA components,
# - overview plots showing the variance explained by significant PCs,
# - example PCA-loading plots,
# - visualisations of high-loading features for medication-related PCs.
#
# The code is organised so that project-specific inputs are concentrated in the setup cells below. Update the paths and, if needed, the feature-column mapping before running the notebook.
# %% [markdown]
# ## 1. Analysis settings
# %%
# Number of independent matched comparisons used in the MATLAB script.
n_matching_repetitions = 10
# A PCA component is treated as robustly significant only if it is significant
# in at least this many matched comparisons.
min_significant_matches = 7
# PCA variance threshold used in the MATLAB script. The notebook uses this
# threshold when summarising/loadings PCs from the PCA decomposition.
min_variance_explained_percent = 90
# Significance threshold after Holm-Bonferroni correction.
alpha = 0.05
# %% [markdown]
# ## 2. Imports and helper functions
# %%
import os
import warnings
import h5py
import mne
import numpy as np
import pandas as pd
import scipy.io
import seaborn as sns
import statsmodels.formula.api as smf
from statsmodels.stats.multitest import multipletests
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib import cm
from matplotlib.colors import Normalize
from mne.channels.layout import _find_topomap_coords
%matplotlib inline
warnings.filterwarnings("ignore")
# %% [markdown]
# ## 3. Define paths
# %%
# Folder containing MATLAB outputs from the published MATLAB script, EEG features and patient data.
# Expected files include:
# cleaned_data/data_new.h5
# cleaned_data/PCA_data.mat
# matched_indices.mat
# Effects_medicines_PCA_OvR.mat
# Effects_medicines_PCA_DN.mat
matlab_output_directory = r'...'
# CSV file describing each EEG feature/column in the EEG feature matrix.
feature_metadata_file = r'...'
# An example EEGLAB .set file. It is used only to read channel locations for
# topographic plots. The participant does not need to be part of the analysis.
example_eeg_file = r'...'
# Output folder for figures generated by this notebook.
figure_directory = os.path.join(matlab_output_directory, "graphs")
if not os.path.isdir(figure_directory):
os.mkdir(figure_directory)
# %% [markdown]
# ## 4. Load cleaned EEG data and metadata
# %%
cleaned_data_file = os.path.join(matlab_output_directory, "cleaned_data", "data_new.h5")
with h5py.File(cleaned_data_file, "r") as h5_file:
data_meds = np.transpose(np.array(h5_file["data_meds"]))
zdata = np.transpose(np.array(h5_file["zdata"]))
data = np.transpose(np.array(h5_file["data"]))
metadata = np.transpose(np.array(h5_file["metadata"]))
print("Medication matrix:", data_meds.shape) #N participants x N medicines
print("Z-scored EEG feature matrix:", zdata.shape) #N participants x N features
print("Participant metadata:", metadata.shape) #N participants x 5
# metadata columns from the MATLAB script:
# column 0 = sex
# column 1 = age
# column 2 = diagnosis
# column 3 = centred recording year
# column 4 = recording site
# %% [markdown]
# ## 5. Load and standardise feature descriptions
# %%
feature_metadata = pd.read_csv(feature_metadata_file)
# Standard column names expected by the rest of the notebook.
standard_feature_columns = [
"measures",
"events",
"channel1",
"channel2",
"freq",
"feat_group",
]
# Map each standard name to the corresponding column in your CSV file.
# If your file already uses the standard names, leave this dictionary as is.
# Example for a different CSV:
# feature_column_map = {
# "measures": "measure_name",
# "events": "condition",
# "channel1": "electrode_a",
# "channel2": "electrode_b",
# "freq": "frequency_hz",
# "feat_group": "feature_type",
# }
feature_column_map = {
"measures": "measures",
"events": "events",
"channel1": "channel1",
"channel2": "channel2",
"freq": "freq",
"feat_group": "feat_group",
}
missing_columns = [
source_column
for source_column in feature_column_map.values()
if source_column not in feature_metadata.columns
]
if missing_columns:
raise ValueError(
"The following columns are missing from feature_metadata_file: "
+ ", ".join(missing_columns)
)
# Rename user-specified columns to the standard names used below.
rename_to_standard = {
source_column: standard_column
for standard_column, source_column in feature_column_map.items()
}
feature_metadata = feature_metadata.rename(columns=rename_to_standard)
# Keep the standard columns first and preserve any additional columns after them.
extra_columns = [
column for column in feature_metadata.columns
if column not in standard_feature_columns
]
feature_metadata = feature_metadata[standard_feature_columns + extra_columns].copy()
# Make sure frequency values are numeric where possible.
feature_metadata["freq"] = pd.to_numeric(feature_metadata["freq"], errors="coerce")
feature_metadata.head()
# %%
# Basic labels used throughout the notebook.
event_names = np.unique(feature_metadata["events"].values)
channel_names = ['Fp1', 'Fp2', 'F7', 'F3', 'Fz', 'F4', 'F8', 'T3','C3', 'Cz', 'C4', 'T4', 'T5', 'P3', 'Pz', 'P4', 'T6', 'O1', 'O2']
measure_names = np.unique(feature_metadata["measures"].values)
frequency_grid = np.arange(4, 60.5, 0.5)
# Medication names must be in the same order as the columns of data_meds.
medication_names = [
"Anticholinergic",
"BDZ",
"NaSSA",
"SARI",
"SNRIs",
"SSRIs",
"AChEi",
"AP atypical",
"opioid",
"sedative-hypnotic",
"TCA",
"AP typical",
"AED_Na",
"AED_Ca"
]
# Short labels used in figure axes where long names would be difficult to read.
medication_short_names = [
"A-Ach",
"BDZ",
"NaSSA",
"SARI",
"SNRIs",
"SSRIs",
"AChEi",
"AP atyp",
"opioid",
"hypno",
"TCA",
"AP typ",
"AED (Na)",
"AED (Ca)"
]
if data_meds.shape[1] != len(medication_names):
raise ValueError(
"The number of medication names does not match the number of columns "
"in data_meds."
)
# %% [markdown]
# ## 6. Load MATLAB outputs
# %%
# Matched participant indices generated by the MATLAB script.
matched_indices_file = scipy.io.loadmat(
os.path.join(matlab_output_directory, "matched_indices.mat"),
squeeze_me=True,
)
matched_indices_ovr = matched_indices_file["matched_indices_OvR"]
matched_indices_drug_naive = matched_indices_file["matched_indices_DN"]
# Effect-size and p-value statistics generated by the MATLAB script.
effects_ovr = scipy.io.loadmat(
os.path.join(matlab_output_directory, "Effects_medicines_PCA_OvR.mat"),
squeeze_me=True,
)
effects_drug_naive = scipy.io.loadmat(
os.path.join(matlab_output_directory, "Effects_medicines_PCA_DN.mat"),
squeeze_me=True,
)
# PCA decomposition generated by the MATLAB script.
pca_file = os.path.join(matlab_output_directory, "cleaned_data", "PCA_data.mat")
with h5py.File(pca_file, "r") as h5_file:
pca_coefficients = np.array(h5_file["PCA_coeff"])
pca_scores = np.array(h5_file["PCA_score"])
pca_explained_variance = np.asarray(h5_file["explained"]).ravel()
print("PCA scores:", pca_scores.shape)
print("PCA coefficients:", pca_coefficients.shape)
print("Explained variance vector:", pca_explained_variance.shape)
# %%
np.shape(matched_indices_drug_naive[8][7])
# %%
# Optional check: print how many matched pairs are available for each medication
# in the first matching repetition of the OvR comparison.
for medication_index, medication_name in enumerate(medication_names):
n_matched_pairs = len(matched_indices_ovr[medication_index][0])
print(medication_name, n_matched_pairs)
# %% [markdown]
# ## 7. Load channel locations for topographic plots
# %%
# MNE uses channel locations stored in example_raw_eeg.info to draw scalp maps.
# The standard_1020 montage is applied so the plots use standard EEG positions.
example_raw_eeg = mne.io.read_raw_eeglab(example_eeg_file, preload=False)
montage_1020 = mne.channels.make_standard_montage("standard_1020")
example_raw_eeg.set_montage(montage_1020)
# Keep only channels that are present both in the feature metadata and in the
# example EEG file. This avoids plotting errors if the feature table contains
# channels absent from the example file.
channel_names = np.array([
channel for channel in channel_names
if channel in example_raw_eeg.ch_names
])
sensor_adjacency, sensor_names = mne.channels.find_ch_adjacency(
example_raw_eeg.info,
ch_type="eeg",
)
print("Channels available for plotting:", len(channel_names))
# %% [markdown]
# ## 8. Calculate mean effect sizes and robustly significant PCs
# %%
mean_hedges_g_ovr = []
mean_hedges_g_drug_naive = []
significant_pc_indices_ovr = []
significant_pc_indices_drug_naive = []
for medication_index in range(len(medication_names)):
# Average Hedges' g across repeated matched samples.
# Rows = PCA components; columns = matching repetitions.
hedges_g_ovr = np.asarray(effects_ovr["Hedges_g_all"][medication_index], dtype=float)
hedges_g_drug_naive = np.asarray(effects_drug_naive["Hedges_g_all"][medication_index], dtype=float)
mean_hedges_g_ovr.append(np.mean(hedges_g_ovr, axis=1))
mean_hedges_g_drug_naive.append(np.mean(hedges_g_drug_naive, axis=1))
# Identify PCs that are significant in at least min_significant_matches independent matched comparisons.
corrected_p_ovr = np.asarray(np.stack(effects_ovr["P_corr"][medication_index]), dtype=float)
corrected_p_drug_naive = np.asarray(np.stack(effects_drug_naive["P_corr"][medication_index]), dtype=float)
significant_pc_indices_ovr.append(
np.where(np.sum(corrected_p_ovr < alpha, axis=0) >= min_significant_matches)[0]
)
significant_pc_indices_drug_naive.append(
np.where(np.sum(corrected_p_drug_naive < alpha, axis=0) >= min_significant_matches)[0]
)
significant_pc_indices_ovr
# %%
np.shape(mean_hedges_g_ovr), np.shape(significant_pc_indices_ovr)
# %% [markdown]
# ### Show effect size and p-value for significant components (supplementary tables)
# %%
significant_effect_tables = {}
for medication_index, medication_name in enumerate(medication_names):
significant_pcs = significant_pc_indices_ovr[medication_index]
best_pcs = pd.DataFrame([], columns=['PC', 'var explained', 'p_OvR', 'g_OvR', 'g_CI_OvR','nmatch_OvR',
'p_DN', 'g_DN', 'Hg_CI_DN', 'nmatch_DN'])
#PC number, explained variance,
#corrected p, Hedges' g, 95% confidence intervals for Hedges' g, number of matches where component was significant
#in OvR comparison and in drug-naive comparison
if len(significant_pcs) > 0:
corrected_p_ovr = np.asarray(np.stack(effects_ovr["P_corr"][medication_index]), dtype=float)
corrected_p_dn = np.asarray(np.stack(effects_drug_naive["P_corr"][medication_index]), dtype=float)
for pc in significant_pcs:
#calculate median p-value
p1 = np.median(corrected_p_ovr[:, pc])
p2 = np.median(corrected_p_dn[:, pc])
#calculate in how many comparisons the effect was significant
nm1 = np.sum(corrected_p_ovr[:, pc] < 0.05, 0)
nm2 = np.sum(corrected_p_dn[:, pc] < 0.05, 0)
#mean effects and their CIs
hg1 = abs(mean_hedges_g_ovr[medication_index][pc])
ci_u = abs(np.mean(np.stack(effects_ovr['CI_upp_all'][medication_index]), 1)[pc])
ci_l = abs(np.mean(np.stack(effects_ovr['CI_low_all'][medication_index]), 1)[pc])
cis_ovr = '[' + str(round(min(ci_l, ci_u), 2)) + ',' + str(round(max(ci_l, ci_u), 2)) + ']'
hg2 = abs(mean_hedges_g_drug_naive[medication_index][pc])
ci_u = abs(np.mean(np.stack(effects_drug_naive['CI_upp_all'][medication_index]), 1)[pc])
ci_l = abs(np.mean(np.stack(effects_drug_naive['CI_low_all'][medication_index]), 1)[pc])
cis_dn = '[' + str(round(min(ci_l, ci_u), 2)) + ',' + str(round(max(ci_l, ci_u), 2)) + ']'
best_pcs.loc[len(best_pcs)] = [pc+1, round(pca_explained_variance[pc], 2),
round(p1, 3), round(hg1, 2), cis_ovr, nm1,
round(p2, 3), round(hg2, 2), cis_dn, nm2]
#sort by OvR effect size
best_pcs.sort_values('g_OvR', ascending = False, inplace=True)
significant_effect_tables[medication_name] = best_pcs
# %%
# Example table. Change "SSRIs" to any medication name from medication_names.
significant_effect_tables["SSRIs"]
# %% [markdown]
# ## 9. Visualization
# %%
# Choose which comparison to summarise in subsequent figures.
# Options:
# "OvR" = medication group versus other medicated participants
# "DN" = medication group versus drug-naive participants
comparison = "OvR"
# %% [markdown]
# ### 9.1 total variance explained by significant PCs
# %%
if comparison == "OvR":
significant_pcs = significant_pc_indices_ovr
elif comparison == "DN":
significant_pcs = significant_pc_indices_drug_naive
else:
raise ValueError("comparison_for_variance_plot must be 'OvR' or 'DN'.")
expl_var_meds = np.zeros(len(significant_pcs))
for medication_index in range(len(significant_pcs)):
if len(significant_pcs[medication_index]) > 0:
expl_var_meds[medication_index] = np.sum(pca_explained_variance[significant_pcs[medication_index]])
#expl_var_meds = expl_var_meds.flatten()
meds_signif = np.array(medication_names)[expl_var_meds > 0]
expl_signif = expl_var_meds[expl_var_meds > 0]
print(meds_signif)
print(expl_signif)
# %%
fig, ax = plt.subplots(figsize=(8, 6))
sns.barplot(
x=meds_signif,
y=expl_signif,
palette="tab10",
ax=ax,
)
ax.set_ylabel("Variance explained by \n significant PCs (%)", size=18)
ax.set_xlabel("Medication group", size=18)
ax.tick_params(axis="y", labelsize=16)
ax.set_xticklabels(
meds_signif,
size=16,
rotation=30,
ha="right",
)
plt.tight_layout()
variance_plot_png = os.path.join(
figure_directory,
"explained_variance_significant_pcs.png",
)
variance_plot_svg = os.path.join(
figure_directory,
"explained_variance_significant_pcs.svg",
)
plt.savefig(variance_plot_png, dpi=900)
plt.savefig(variance_plot_svg, dpi=900, format="svg")
plt.show()
# %% [markdown]
# ### 9.2 Identify high-loading features for each PCA component
# %%
# Count how many features belong to each feature group in the full dataset.
fgs, fgcnts = np.unique(feature_metadata['feat_group'].values, return_counts=True)
fgs, fgcnts
# %%
# Feature groups shown in summary bar plots.
# The first entry combines frequency-binned and frequency-resolved coherence
# features into one broader "coherence" category.
best_feats_allPC = {}
props_allPC = []
nbpca = np.where(np.cumsum(pca_explained_variance) > min_variance_explained_percent)[0][0]
for pc in range(0,nbpca):
best_feats_allPC[pc] = {}
feats_comp = pca_coefficients[pc, :]
#ADD feature components to each PC
feature_metadata_by_pc = feature_metadata.copy()
feature_metadata_by_pc['feature_eigenvalues'] = feats_comp
#calculate mean and standard deviation for absolsute PC coefficients and find features that have coefficients higher than 2sd above the mean
mm = np.mean(abs(feats_comp))
sd = np.std(abs(feats_comp))
best_feats = abs(feats_comp)>mm+2*sd
#change indices to feature descriptions
best_feats_allPC[pc] = feature_metadata_by_pc.iloc[best_feats, :]
#calculate the contribution of each feature group in this PCA relative to the contribution of this group in whole data
pca_fgs, pca_fgcnts = np.unique(best_feats_allPC[pc]['feat_group'].values, return_counts=True)
props_l_pca = [0,0,0,0,0]
for fi in range(0, len(pca_fgcnts)):
fi2 = list(fgs).index(pca_fgs[fi])
prop = pca_fgcnts[fi]/fgcnts[fi2]
props_l_pca[fi2] = prop
props_allPC.append(props_l_pca)
# %% [markdown]
# ### 9.3 Plot example PCA components
# %%
comps_to_plot = 10 #how many components to plot (integer between 1 and nbpca)
example_pca_directory = os.path.join(figure_directory, "example_PCA_components")
if not os.path.isdir(example_pca_directory):
os.makedirs(example_pca_directory)
# Colorbar placement within topomap axes.
colorbar_x_start = 1.03
colorbar_x_width = 0.04
colorbar_y_start = 0.05
colorbar_y_height = 0.90
for pc in range(comps_to_plot):
feats_pca = best_feats_allPC[pc]
fig, ax = plt.subplots(1,3, figsize=(18, 5), width_ratios=[1.1,1.1,0.8])
#FREQUENCY
freq_in_range_oi = []
for fqi, fq in zip(feats_pca.index.values, feats_pca['freq'].values):
if fq in frequency_grid:
freq_in_range_oi.append(fqi)
feats_pca1 = feats_pca.copy()
feats_pca1 = feats_pca1.loc[freq_in_range_oi, :]
freq_values, freq_counts = np.unique(feats_pca1['freq'].astype(float).values, return_counts = True)
#frequency plot
ax1 = ax[1]
ax1.set_xlim(0, 61)
ax1.bar(freq_values, freq_counts, width=0.5, color='navy')
ax1.set_xticks(range(0,61,10))
ax1.tick_params(axis='both', labelsize=18)
ax1.set_xlabel("Frequency [Hz]", fontsize=20)
ax1.set_ylabel('N features \n (eigenvalue > M+2SD)', fontsize=20)
ax1.text(0.5, 1.02, "Frequencies", horizontalalignment='center', verticalalignment='bottom', size=20, transform=ax1.transAxes)
#CHANNELS
feats_pca2 = feats_pca.copy()
feats_pca2['channel1'] = feats_pca['channel2']
feats_pca2 = feats_pca2.append(feats_pca)
names_ch, counts_ch = np.unique(feats_pca2['channel1'], return_counts = True)
counts_ch2 = np.zeros((np.shape(channel_names)))
for chi, ch in enumerate(channel_names):
if ch in names_ch:
counts_ch2[chi] = counts_ch[np.where(names_ch==ch)]
#channels topomap plot
ax2 = ax[2]
im, cm2 = mne.viz.plot_topomap(counts_ch2, example_raw_eeg.info, names=channel_names,
cmap = 'inferno', axes=ax2, show=False)
#increase fontsize for channel names
for tt in plt.findobj(fig, matplotlib.text.Text):
if tt.get_text() in example_raw_eeg.ch_names:
tt.set_fontsize(16)
#colorbar
cbar_ax = ax2.inset_axes([colorbar_x_start, colorbar_y_start, colorbar_x_width, colorbar_y_height])
clb = fig.colorbar(im, cax=cbar_ax)
clb.set_label("N features \n (eigenvalue > M+2SD)", size=18)
clb.ax.tick_params(labelsize=16)
ax2.text(0.5, 1.02, "Channels", horizontalalignment='center', verticalalignment='bottom', size=20, transform=ax2.transAxes)
#feature groups
ax3 = ax[0]
fg_names = ['coherence', 'phase \n synchrony', 'spectral', 'nonlinear']
fg_counts = [np.mean([props_allPC[pc][0]*100, props_allPC[pc][1]*100]), props_allPC[pc][2]*100, props_allPC[pc][3]*100, props_allPC[pc][4]*100]
ax3.barh(fg_names, fg_counts, color='darkgreen', edgecolor = 'white')
ax3.tick_params(axis='both', labelsize=20)
ax3.set_xlabel("% of all features \n (eigenvalue > M+2SD)", fontsize=20)
ax3.text(0.5, 1.02, "Feature types", horizontalalignment='center', verticalalignment='bottom', size=20, transform=ax3.transAxes)
evar = round(pca_explained_variance[pc], 2)
fig.suptitle('PCA' + str(pc+1) + ' (explained variance ' + str(evar) + '%)', fontsize=22, fontweight='bold')
plt.tight_layout()
plt.savefig(os.path.join(example_pca_directory, 'PCA' + str(pc+1) + '.svg'), dpi=900, format="svg")
plt.savefig(os.path.join(example_pca_directory, 'PCA' + str(pc+1) + '.png'), dpi=900)
plt.close()
#Note: in graph titles and filenames, PC's are labeled from 1, not from 0
# %% [markdown]
# ### 9.4 Plot high-loading features for significant medication-related PCs
# %%
#Function for connectivity plots
from matplotlib import cm
from matplotlib.colors import Normalize
from mne.channels.layout import _find_topomap_coords
def plot_sensor_connectivity_topomap(
info, connectivity_matrix, top_n=30,
v_min=0.015, v_max=0.050, cmap='viridis',
ax=None, save_path=None, cbar_label="z"):
"""
Plot 2D EEG sensor connectivity on top of sensor topomap.
Parameters
----------
info : instance of mne.Info
EEG channel info.
connectivity_matrix : ndarray, shape (n_channels, n_channels)
Symmetric matrix of connectivity values between EEG channels.
top_n : int
Number of top connections to display.
v_min, v_max : float
Color scale limits for connectivity strength.
cmap : str
Colormap to use.
ax : matplotlib Axes | None
If provided, plot into this axis (useful for subplots).
If None, a new figure is created.
save_path : str or None
If provided, path to save the figure.
"""
# Pick EEG sensors and get projected 2D coords
picks = mne.pick_types(info, meg=False, eeg=True)
n_channels = len(picks)
assert connectivity_matrix.shape == (n_channels, n_channels), \
f"Expected matrix of shape ({n_channels}, {n_channels})"
coords = _find_topomap_coords(info, picks)
# Extract top-N strongest connections
triu_inds = np.triu_indices(n_channels, k=1)
conn_vals = connectivity_matrix[triu_inds]
top_idx = np.argsort(abs(conn_vals))[-top_n:]
i1, i2 = triu_inds[0][top_idx], triu_inds[1][top_idx]
top_vals = conn_vals[top_idx]
# Create axis if not provided
own_fig = False
if ax is None:
fig, ax = plt.subplots(figsize=(5.5, 5))
own_fig = True
else:
fig = ax.figure
# Plot EEG sensor positions
mne.viz.plot_sensors(info, kind='topomap', show=False, axes=ax)
# Normalize colors
norm = Normalize(vmin=v_min, vmax=v_max)
cmap_func = cm.ScalarMappable(norm=norm, cmap=cmap)
# Draw connectivity lines
for idx_a, idx_b, val in zip(i1, i2, top_vals):
x1, y1 = coords[idx_a]
x2, y2 = coords[idx_b]
ax.plot([x1, x2], [y1, y2], color=cmap_func.to_rgba(val), linewidth=2)
# Colorbar
#ticksl = np.round(np.arange(v_min, v_max, 0.001), 3)
#ticksl = np.arange(v_min, v_max)
#cbar = fig.colorbar(cmap_func, ax=ax, fraction=0.046, pad=0.03, ticks=ticksl)
cbar = fig.colorbar(cmap_func, ax=ax, fraction=0.047, pad=0.02, shrink=0.9) #, ticks=ticksl)
#cbar.set_label('Mean eigenvalue', size=18)
cbar.set_label(cbar_label, size=20, labelpad=-50, y=1.05) #, rotation=90) #rotation??
cbar.ax.tick_params(labelsize=16)
# Save or show
if save_path:
if save_path.endswith('.svg'):
fig.savefig(save_path, dpi=900, bbox_inches='tight', format='svg')
else:
fig.savefig(save_path, dpi=900, bbox_inches='tight')
print(f"Saved plot to {save_path}")
elif own_fig:
plt.show()
if own_fig:
plt.close(fig)
# %%
if comparison == "OvR":
signif_pcs = significant_pc_indices_ovr
mean_effsizes = mean_hedges_g_ovr
comparison_label = "others"
elif comparison == "DN":
signif_pcs = significant_pc_indices_drug_naive
mean_effsizes = mean_hedges_g_drug_naive
comparison_label = "drug-naive"
else:
raise ValueError("comparison_for_pc_plots must be 'OvR' or 'DN'.")
significant_pc_figure_directory = os.path.join(
figure_directory,
"significant_PCs_{0}".format(comparison),
)
if not os.path.isdir(significant_pc_figure_directory):
os.makedirs(significant_pc_figure_directory)
plot_panel_titles = [
"Types of features",
"Coherence frequencies",
"PSD frequencies",
"Nonlinear channels",
"Coherence connections",
"PSD channels",
]
subplot_labels = ["A", "B", "C", "D", "E", "F"]
#feature_group_plot_labels = [label for label, _ in plot_feature_group_order]
top_channels_psd_by_medication = {}
top_channels_nonlinear_by_medication = {}
# Measures whose sign is reversed to match the interpretation used in the manuscript (higher values mean greater complexity).
#Edit this list if the feature definitions change.
nonlinear_measures_with_reversed_direction = ["Hjorth_complexity", "DFA_"]
# %%
for medication_index, medication_name in enumerate(medication_names):
top_channels_psd_by_medication[medication_index] = []
top_channels_nonlinear_by_medication[medication_index] = []
if not os.path.isdir(os.path.join(significant_pc_figure_directory, medication_name)):
os.mkdir(os.path.join(significant_pc_figure_directory, medication_name))
for pc in signif_pcs[medication_index]:
high_loading_features = best_feats_allPC[pc]
# Orient PCA loadings so that positive values correspond to larger
# values in the medication group relative to the comparison group.
pca_sign = 1
if mean_effsizes[medication_index][pc] < 0:
pca_sign = -1
high_loading_features['feature_eigenvalues'] = high_loading_features['feature_eigenvalues']*pca_sign
#reverse the sign for complexity features for which higher values mean less complexity
reverse_measure_mask = high_loading_features["measures"].isin(nonlinear_measures_with_reversed_direction)
high_loading_features.loc[reverse_measure_mask, 'feature_eigenvalues'] = -high_loading_features.loc[reverse_measure_mask, 'feature_eigenvalues'].values
#FIGURE
fig = plt.figure(figsize=(20, 12)) #, constrained_layout=True)
outer = gridspec.GridSpec(2, 3, wspace=0.3, hspace=0.3, height_ratios = [3,2])
conn_temptate = np.zeros((19,19))
for mg, sbi in zip(['coh_fq', 'linear_psd'], [1,2]):
df_mg = high_loading_features.loc[high_loading_features['feat_group']==mg, :]
df_mg['abs_pcacoef'] = abs(df_mg['feature_eigenvalues'])
df_mg_p = df_mg.loc[df_mg['feature_eigenvalues'] > 0, :]
df_mg_m = df_mg.loc[df_mg['feature_eigenvalues'] < 0, :]
#best channels
channels_mean_coeff = df_mg.groupby(['channel1']).mean('abs_pcacoef').sort_values(by = ['abs_pcacoef'], ascending=False)
top_channels_psd_by_medication[medication_index].append(channels_mean_coeff.index.values[0:5])
# ---------------Frequency plots (for psd and connectivity)------------------------------
nmf_p, cntf_p = np.unique(df_mg_p['freq'].astype(float).values, return_counts = True)
nmf_m, cntf_m = np.unique(df_mg_m['freq'].astype(float).values, return_counts = True)
inner = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=outer[sbi], wspace=0, hspace=0.2)
ax1 = plt.Subplot(fig, inner[0])
ax2 = plt.Subplot(fig, inner[1])
#plot frequencies where medicine group has higher values than the comparison group
#in the positive axis and red color and where it has lower values in negative axis and blue color
ax1.bar(nmf_p, cntf_p, width=0.5, color='indianred')
ax2.bar(nmf_m, cntf_m, width=0.5, color='royalblue')
ax1.set_xlim(0, 61)
ax1.set_xticks(range(0,61,10))
ax1.set_xticklabels([])
ax1.tick_params(axis='y', labelsize=18, direction='out', length=5)
ax1.set_ylabel(medication_short_names[medication_index] + ' > ' + comparison_label, fontsize=19, loc='center')
ax2.set_xlim(0, 61)
ax2.set_xticks(range(0,61,10))
ax2.set_xlabel("Frequency [Hz]", fontsize=20, loc='center')
ax2.set_ylabel(medication_short_names[medication_index] + ' < ' + comparison_label, fontsize=19, loc='center')
#inverted joined y axis
ax2.invert_yaxis()
ax2.tick_params(axis='both', labelsize=18, direction='out', length=5)
ax2.xaxis.tick_top()
#ensure the same limits for both positive and negative axes
lim = max(max(ax1.get_ylim()), max(ax2.get_ylim()))
ax1.set_ylim(0, lim)
ax2.set_ylim(lim, 0)
ax1.text(0.5, 1.02, plot_panel_titles[sbi], horizontalalignment='center', verticalalignment='bottom', size=20, transform=ax1.transAxes)
ax1.text(-0.1, 1.02, subplot_labels[sbi], weight = 'bold', horizontalalignment='right', verticalalignment='bottom', size=22, transform=ax1.transAxes)
ax1.text(-0.01, 1.00, 'N', horizontalalignment='right', verticalalignment='bottom', size=18, transform=ax1.transAxes)
#move the x-axis label to the center and remove the double zero label
yticks2 = ax2.yaxis.get_major_ticks()
yticks2[0].set_visible(False)
fig.add_subplot(ax1)
fig.add_subplot(ax2)
# ---------------Sensor topomap plots (for psd and connectivity)----------------------------
inner = gridspec.GridSpecFromSubplotSpec(1, 1, subplot_spec=outer[sbi+3], wspace=0, hspace=0) #, width_ratios=[5, 1])
ax = plt.Subplot(fig, inner[0])
cnt_ch = np.zeros((len(channel_names), ))
nm_ch1, cnt_ch1 = np.unique(df_mg['channel1'], return_counts = True)
nm_ch2, cnt_ch2 = np.unique(df_mg['channel2'], return_counts = True)
nm_ch1 = list(nm_ch1)
nm_ch2 = list(nm_ch2)
cnt_ch1 = list(cnt_ch1)
cnt_ch2 = list(cnt_ch2)
for chi in range(0, len(channel_names)):
if channel_names[chi] not in nm_ch1:
nm_ch1.append(channel_names[chi])
cnt_ch1.append(0)
if channel_names[chi] not in nm_ch2:
nm_ch2.append(channel_names[chi])
cnt_ch2.append(0)
cnt_ch[chi] = cnt_ch1[nm_ch1.index(channel_names[chi])] + cnt_ch2[nm_ch2.index(channel_names[chi])]
if mg == 'linear_psd':
im, cm2 = mne.viz.plot_topomap(cnt_ch, example_raw_eeg.info, names=channel_names,
cmap = 'Blues', axes=ax, show=False)
#colorbar
cbar_ax = ax.inset_axes([colorbar_x_start, colorbar_y_start, colorbar_x_width, colorbar_y_height])
clb = fig.colorbar(im, cax=cbar_ax)
clb.ax.tick_params(labelsize=16)
cbar_ax.text(0, 1.02, 'N', horizontalalignment='left', verticalalignment='bottom', size=18, transform=cbar_ax.transAxes)
elif mg == 'coh_fq':
for chi1 in range(len(channel_names)):
for chi2 in range(chi1, len(channel_names)):
df_ch = df_mg.loc[df_mg['channel1'] == channel_names[chi1]].loc[df_mg['channel2'] == channel_names[chi2]]
conn_temptate[chi1, chi2] = len(df_ch)
conn_temptate[chi2, chi1] = len(df_ch)
plot_sensor_connectivity_topomap(example_raw_eeg.info, conn_temptate, top_n=30,
v_min = np.min(conn_temptate[conn_temptate > 0]), v_max=np.max(conn_temptate),
cmap='Blues', ax=ax)
ax.text(0.5, 1.02, plot_panel_titles[sbi+3], horizontalalignment='center', verticalalignment='bottom', size=20, transform=ax.transAxes)
ax.text(-0.1, 1.02, subplot_labels[sbi+3], weight = 'bold', horizontalalignment='right', verticalalignment='bottom', size=22, transform=ax.transAxes)
fig.add_subplot(ax)
# -------------- feature groups --------------
inner = gridspec.GridSpecFromSubplotSpec(1, 1,
subplot_spec=outer[0], wspace=0, hspace=0)
ax = plt.Subplot(fig, inner[0])
cnt_fg = [np.mean([props_allPC[pc][0]*100, props_allPC[pc][1]*100]), props_allPC[pc][2]*100, props_allPC[pc][3]*100, props_allPC[pc][4]*100]
ax.barh(fg_names, cnt_fg, color='navy', edgecolor = 'white')
ax.tick_params(axis='both', labelsize=20)
ax.set_xlabel('% of all features', fontsize=20)
ax.text(0.5, 1.02, plot_panel_titles[0], horizontalalignment='center', verticalalignment='bottom', size=20, transform=ax.transAxes)
ax.text(-0.1, 1.02, subplot_labels[0], weight = 'bold', horizontalalignment='right', verticalalignment='bottom', size=22, transform=ax.transAxes)
fig.add_subplot(ax)
# --------- Nonlinear features - topomap ------------
inner = gridspec.GridSpecFromSubplotSpec(1, 1,
subplot_spec=outer[3], wspace=0, hspace=0)
ax3 = plt.Subplot(fig, inner[0])
df_mg = high_loading_features.loc[high_loading_features['feat_group']=='nonlinear', :]
df_mg['abs_pcacoef'] = abs(df_mg['feature_eigenvalues'])
df_mg_p = df_mg.loc[df_mg['feature_eigenvalues'] > 0, :]
df_mg_m = df_mg.loc[df_mg['feature_eigenvalues'] < 0, :]
#best channels
channels_mean_coeff = df_mg.groupby(['channel1']).mean('abs_pcacoef').sort_values(by = ['abs_pcacoef'], ascending=False)
top_channels_nonlinear_by_medication[medication_index].append(channels_mean_coeff.index.values[0:5])
cnt_ch = np.zeros((len(channel_names), ))
nm_chp, cnt_chp = np.unique(df_mg_p['channel1'], return_counts = True)
nm_chm, cnt_chm = np.unique(df_mg_m['channel1'], return_counts = True)
for chi in range(0, len(channel_names)):
cnt_ch[chi] = 0
if channel_names[chi] in nm_chp:
cnt_ch[chi] = cnt_ch[chi] + cnt_chp[np.where(nm_chp == channel_names[chi])[0][0]]
if channel_names[chi] in nm_chm:
cnt_ch[chi] = cnt_ch[chi] - cnt_chm[np.where(nm_chm == channel_names[chi])[0][0]]
#topomap
vmax = np.max([abs(np.nanmin(cnt_ch)), abs(np.nanmax(cnt_ch))])
im, cm2 = mne.viz.plot_topomap(cnt_ch, example_raw_eeg.info, names=channel_names,
vlim=(-vmax, vmax),
cmap = 'seismic', axes=ax3, show=False)
#colorbar
cbar_ax = ax3.inset_axes([colorbar_x_start, colorbar_y_start, colorbar_x_width, colorbar_y_height])
clb = fig.colorbar(im, cax=cbar_ax)
clb.ax.tick_params(labelsize=16)
cbar_ax.text(0, 1.02, 'N', horizontalalignment='left', verticalalignment='bottom', size=18, transform=cbar_ax.transAxes)
ax3.text(0.5, 1.02, plot_panel_titles[3], horizontalalignment='center', verticalalignment='bottom', size=20, transform=ax3.transAxes)
ax3.text(-0.02, 1.02, subplot_labels[3], weight = 'bold', horizontalalignment='right', verticalalignment='bottom', size=22, transform=ax3.transAxes)
fig.add_subplot(ax3)
#increase fontsize for channelnames
for tt in plt.findobj(fig, matplotlib.text.Text):
if tt.get_text() in channel_names:
tt.set_fontsize(16)
evar = round(pca_explained_variance[pc], 2)
fig.suptitle('PC' + str(pc+1) + ' (explained variance ' + str(evar) + '%)', fontsize=22, fontweight='bold')
plt.subplots_adjust(left=0.15, top = 0.91, hspace = 0.1)
plt.savefig(os.path.join(significant_pc_figure_directory, medication_name , 'PCA' + str(pc+1) + '.svg'), dpi=900, format="svg")
plt.savefig(os.path.join(significant_pc_figure_directory, medication_name , 'PCA' + str(pc+1) + '.png'), dpi=900)
plt.show()
plt.close()
# %% [markdown]
# ## 10. Single-feature graphs
#
# Plots to visualise differences between medication groups and comparison groups in single features.
# As the chosen viualisations often involve mean features from the several electrodes, the effects are calculated here, instead of using matlab output.
# %%
#Function that calculates the required group values from zdata
#input:
#feat_idx - numeric indexes of features in the data. can be multidimensional, but channels (or other dimension to be averaged) schold be in the last axis
#matched_idx - indexes of persons from given groups
#transform_log - boolean, whether to calculate log10 of the data values, defaul-no
#average_chennels - boolean, whether the result should be averaged over channels, default-yes
def get_meanfeaturevals_and_sem(feat_idx, zdata, metadata, matched_idx_ovr, matched_idx_dn, transform_log=False, average_channels=True):
psd_nodrug = []
psd_drug = []
psd_drug_naive = []
semI_drug = []
semI_nodrug = []
semI_drug_naive = []
if average_channels:
#mean of feature over channels
psd_meanchann = np.mean(zdata[:, feat_idx], axis=-1)
else:
psd_meanchann = zdata[:, feat_idx]
if transform_log:
#rescale the mean psd so it above 0, to be able to calculate log10 of it
psd_meanchann = psd_meanchann + abs(np.min(psd_meanchann)) + 0.001
#preloacte p matrices
dim_feat = 1
if len(np.shape(psd_meanchann)) > 1:
dim_feat = np.shape(psd_meanchann)[1]
p_ovr = np.zeros((len(matched_idx_ovr), dim_feat))
p_dn = np.zeros((len(matched_idx_ovr), dim_feat))
#Create feature matrix by adressing zdata
for i in range(len(matched_idx_ovr)):
feature_drug = psd_meanchann[matched_idx_ovr[i][:, 0] - 1] #-1 converts the matlab index to python
feature_nodrug = psd_meanchann[matched_idx_ovr[i][:, 1] - 1]
feature_drug_naive = psd_meanchann[matched_idx_dn[i][:, 1] - 1]
if transform_log==1:
feature_drug = np.log10(feature_drug)
feature_nodrug = np.log10(feature_nodrug)
feature_drug_naive = np.log10(feature_drug_naive)
#Mean over participants
psd_drug.append(np.mean(feature_drug, axis=0))
psd_nodrug.append(np.mean(feature_nodrug, axis=0))
psd_drug_naive.append(np.mean(feature_drug_naive, axis=0))
#standard error of the mean for the matrices
semI_drug.append(scipy.stats.sem(feature_drug))
semI_nodrug.append(scipy.stats.sem(feature_nodrug))
semI_drug_naive.append(scipy.stats.sem(feature_drug_naive))
# --------------- Significance for the mean of channels:
#Create dataframe for hierarchical regression model for each frequency
df_drug = pd.DataFrame(feature_drug)
df_drug["year_c"] = metadata[matched_idx_ovr[i][:, 0] - 1, 3]
df_drug["site"] = metadata[matched_idx_ovr[i][:, 0] - 1, 4]
df_drug["group"] = "drugY"
df_nodrug = pd.DataFrame(feature_nodrug)
df_nodrug["year_c"] = metadata[matched_idx_ovr[i][:, 1] - 1, 3]
df_nodrug["site"] = metadata[matched_idx_ovr[i][:, 1] - 1, 4]
df_nodrug["group"] = "drugN"
df_drugnai = pd.DataFrame(feature_drug_naive)
df_drugnai["year_c"] = metadata[matched_idx_dn[i][:, 1] - 1, 3]
df_drugnai["site"] = metadata[matched_idx_dn[i][:, 1] - 1, 4]
df_drugnai["group"] = "drug_naive"
dfi_ovr = pd.concat([df_drug, df_nodrug])
dfi_dn = pd.concat([df_drug, df_drugnai])
#Calculate singnificance for each frequency with linear mixed model:
for fqi in range(dim_feat):
#OVR
dfi_ovr_fqi = dfi_ovr.iloc[:, -3:]
dfi_ovr_fqi["psd"] = dfi_ovr.iloc[:, fqi]
md = smf.mixedlm("psd ~ group + year_c", dfi_ovr_fqi, groups=dfi_ovr_fqi["site"])
mdf = md.fit()
p_ovr[i, fqi] = mdf.pvalues['group[T.drugY]']
#Drug-naive
dfi_dn_fqi = dfi_dn.iloc[:, -3:]
dfi_dn_fqi["psd"] = dfi_dn.iloc[:, fqi]
md = smf.mixedlm("psd ~ group + year_c", dfi_dn_fqi, groups=dfi_dn_fqi["site"])
mdf = md.fit()
p_dn[i, fqi] = mdf.pvalues['group[T.drug_naive]']
#Mean over n_matching_repetitions:
M_psd_drug = np.mean(psd_drug, axis=0)
M_psd_nodrug = np.mean(psd_nodrug, axis=0)
M_psd_drug_naive = np.mean(psd_drug_naive, axis=0)
M_semI_drug = np.mean(semI_drug, axis=0)
M_semI_nodrug = np.mean(semI_nodrug, axis=0)
M_semI_drug_naive = np.mean(semI_drug_naive, axis=0)
#Stack above matrices into one for mean features and one for SEM:
M_psd_all = np.vstack([M_psd_drug, M_psd_nodrug, M_psd_drug_naive])
M_sem_all = np.vstack([M_semI_drug, M_semI_nodrug, M_semI_drug_naive])
#corrected p-vales
rej, p_ovr_corr, sth, sth = multipletests(p_ovr.flatten(), method='holm')
p_ovr_corr = p_ovr_corr.reshape(np.shape(p_ovr))
rej, p_dn_corr, sth, sth = multipletests(p_dn.flatten(), method='holm')
p_dn_corr = p_dn_corr.reshape(np.shape(p_dn))
print('Returns: \n', "feature matrix of shape:", np.shape(M_psd_all), ', groups: drug, others, drug-naive), \n',
"SEM matrix of shape: ", np.shape(M_sem_all), '\n',
"p-values for ovr and drug-naive comparison for each repetition, shape:", np.shape(p_dn))
return(M_psd_all, M_sem_all, p_ovr_corr, p_dn_corr)
# %% [markdown]
# ## PSD
# %%
freq_grid_nolinenoise = list(np.arange(4, 48.5, 0.5))
freq_grid_nolinenoise.extend(np.arange(52, 60.5, 0.5))
# %%
len(freq_grid_nolinenoise)
# %%
#Find indexes in zdata for PSD features
#widmo indexy
PSD_idx = {}
for ev in event_names:
PSD_idx[ev] = np.zeros((np.shape(freq_grid_nolinenoise)[0], np.shape(channel_names)[0]))
for ch in range(len(channel_names)):
for fq in range(len(freq_grid_nolinenoise)):
PSD_idx[ev][fq, ch] = int(feature_metadata.loc[((feature_metadata['measures']=='PSD') & (feature_metadata['events']==ev) & (feature_metadata['channel1']==channel_names[ch]) & (feature_metadata['freq']==freq_grid_nolinenoise[fq]))].index[0])
# %%
#Set which electrodes / events / medicines to plot. Here, we will use an example from Figure 4 in our manuscript.
savepath_psd = os.path.join(figure_directory, 'single_features_comparisons')
if not os.path.isdir(savepath_psd):
os.mkdir(savepath_psd)
chpsd = ['F3', 'F4', 'F7', 'F8'] #best in PCA 10 in psd
ev = 'OZ'
medname = "BDZ"
colors = ['forestgreen', 'royalblue', 'orange']
#Convert names of required channels to indices in data
chpsd_idx = [list(channel_names).index(ch) for ch in chpsd]
feat_idx = PSD_idx[ev][:, chpsd_idx].astype(int)
#find matched indices for given drug
matched_idx_ovr = matched_indices_ovr[medication_names.index(medname)]
matched_idx_dn = matched_indices_drug_naive[medication_names.index(medname)]
#Create matrices of mean log PSD values and their SEM for each group, and p-values for mean psd (mean over channels)
#Note: it can take a while as new linear mixed model is calculated for each feature.
M_psd_all, M_sem_all, p_ovr_corr, p_dn_corr = get_meanfeaturevals_and_sem(feat_idx, data, metadata, matched_idx_ovr, matched_idx_dn, transform_log=1, average_channels=1)
# %%
#PLOT
fig, (ax1, ax15, ax2) = plt.subplots(3, 1, sharex=True, height_ratios=[9,0.5,0.5], figsize=(10*0.85,6*0.85))# hspace=0.2)
for group in range(np.shape(M_psd_all)[0]):
ax1.plot(freq_grid_nolinenoise, M_psd_all[group], color=colors[group])
ax1.legend([medname, 'other', 'drug-naive'], fontsize="16")
for group in range(np.shape(M_psd_all)[0]):
ax1.fill_between(freq_grid_nolinenoise, M_psd_all[group]+2*M_sem_all[group], M_psd_all[group]-2*M_sem_all[group], alpha=.4, color=colors[group])
#remove upper and right box edges
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)
#set y limits so the graph is readable
xmax = 40
plt.xlim((4, xmax))
ax1.set_ylim(min(M_psd_all[:, 2*xmax-2]), max(ax1.get_ylim()))
ax1.set_ylabel('Z-scored log(power)', fontsize=20) #unit of power: [V$^2$]
ax1.tick_params(axis='both', labelsize=18)
tit = "logPSD " + str(chpsd) + " during " + ev
ax15.fill_between(x=freq_grid_nolinenoise, y1=0, y2=1, where = np.sum(p_dn_corr<alpha, axis=0) >= min_significant_matches, color=colors[2], alpha=0.5)
ax15.text(-0.02, 0.5, "p<0.05 vs drug-naive", horizontalalignment='right', verticalalignment='center', size=14, transform=ax15.transAxes)
ax2.fill_between(x=freq_grid_nolinenoise, y1=0, y2=1, where = np.sum(p_ovr_corr<alpha, axis=0) >= min_significant_matches, color=colors[1], alpha=0.5)
ax2.text(-0.02, 0.5, "p<0.05 vs others", horizontalalignment='right', verticalalignment='center', size=14, transform=ax2.transAxes)
ax2.set_xlabel("Frequency [Hz]", fontsize=20)
ax2.tick_params(axis='x', labelsize=18)
ax2.set_yticklabels([])
ax15.set_yticklabels([])
plt.tight_layout()
plt.subplots_adjust(hspace=0.02)
plt.savefig(os.path.join(savepath_psd, medname + '_' + tit + '.png'))
plt.savefig(os.path.join(savepath_psd, medname + '_' + tit + '.svg'))
plt.show()
plt.close()
# %% [markdown]
# ## Nonlinear
# %%
nonlin_idx = {}
for ev in event_names:
nonlin_idx[ev] = {}
for m in np.unique(feature_metadata.loc[feature_metadata["feat_group"]=='nonlinear', "measures"]):
aa = np.zeros(np.shape(channel_names))
for ch in range(0, len(channel_names)):
aa[ch] = feature_metadata.loc[((feature_metadata['measures']==m) & (feature_metadata['events']==ev) & (feature_metadata['channel1']==channel_names[ch]))].index.values[0]
nonlin_idx[ev][m] = aa
# %%
nonlin_idx[ev].keys()
# %%
chpsd = ['Fz', 'F3', 'Fp1', 'F4', 'Fp2']
ev = 'OZ'
medname = "BDZ"
measure= "HFD_"
#colors = ['forestgreen', 'royalblue', 'orange']
#Convert names of required channels to indices in data
chpsd_idx = [list(channel_names).index(ch) for ch in chpsd]
feat_idx = nonlin_idx[ev][measure][chpsd_idx].astype(int)
#find matched indices for given drug
matched_idx_ovr = matched_indices_ovr[medication_names.index(medname)]
matched_idx_dn = matched_indices_drug_naive[medication_names.index(medname)]
#Create matrices of mean log PSD values and their SEM for each group, and p-values for mean psd (mean over channels)
#Note: it can take a while as new linear mixed model is calculated for each feature.
M_hfd_all, M_sem_all, p_ovr_corr, p_dn_corr = get_meanfeaturevals_and_sem(feat_idx, zdata, metadata, matched_idx_ovr, matched_idx_dn, transform_log=0, average_channels=1)
# %%
#plot
medname_short = medication_short_names[medication_names.index(medname)]
fig = plt.figure(figsize=(6,6))
plt.bar([medname_short, "Others", "Drug \n naive"], M_hfd_all.flatten(), width=0.4)
plt.errorbar([0,1,2], M_hfd_all.flatten(), yerr=M_sem_all.flatten(), fmt="o")
plt.ylabel(measure.replace('_', ' ') + "\n (z-scored)", fontsize=25)
plt.xticks(fontsize=24)
plt.yticks(fontsize=24)
plt.tight_layout()
tit_png = medname + '_' + measure + '_' + str(chpsd) + ev + '.png'
tit_svg = medname + '_' + measure + '_' + str(chpsd) + ev +'.svg'
plt.savefig(os.path.join(savepath_psd, tit_svg))
plt.savefig(os.path.join(savepath_psd, tit_png))
plt.show()
plt.close()
# %%
np.unique(feature_metadata['measures'])
# %%
#Indexes for connectivity measures with 0.5 Hz resolution
freq_all = list(np.unique(feature_metadata.loc[feature_metadata['feat_group']=='coh_fq', 'freq']))
measures = np.unique(feature_metadata.loc[feature_metadata['feat_group']=='coh_fq', 'measures'])
COH_idx = {}
#channel pair indexes
channel_pairs = {}
ichp = 0
for ch1 in range(0, len(channel_names)):
for ch2 in range(ch1+1, len(channel_names)):
channel_pairs[ichp] = (channel_names[ch1], channel_names[ch2])
ichp+=1
for ev in event_names:
COH_idx[ev] = {}
for m in measures:
COH_idx[ev][m] = np.zeros((len(channel_pairs), len(freq_all)))
for fqz in freq_all:
for ichp in channel_pairs:
COH_idx[ev][m][ichp, freq_all.index(fqz)] = feature_metadata.loc[((feature_metadata['measures']==m) & (feature_metadata['events']==ev) & (feature_metadata['channel1']==channel_pairs[ichp][0]) & (feature_metadata['channel2']==channel_pairs[ichp][1]) & (feature_metadata['freq']== fqz))].index[0]
print(m, 'index shape:', np.shape(COH_idx[ev][m]))
# %%
#Indexes for binned frequency measures
freq_all_bin = list(np.unique(feature_metadata.loc[feature_metadata['feat_group']=='coh_bins', 'freq']))
measures = np.unique(feature_metadata.loc[feature_metadata['feat_group']=='coh_bins', 'measures'])
COH_bin_idx = {}
for ev in event_names:
COH_bin_idx[ev] = {}
for m in measures:
COH_bin_idx[ev][m] = np.zeros((len(channel_pairs), len(freq_all_bin)))
for fqz in freq_all_bin:
for ichp in channel_pairs:
COH_bin_idx[ev][m][ichp, freq_all_bin.index(fqz)] = feature_metadata.loc[((feature_metadata['measures']==m) & (feature_metadata['events']==ev) & (feature_metadata['channel1']==channel_pairs[ichp][0]) & (feature_metadata['channel2']==channel_pairs[ichp][1]) & (feature_metadata['freq']== fqz))].index[0]
print(m, 'index shape:', np.shape(COH_bin_idx[ev][m]))
# %%
ev = 'OZ'
medname = "BDZ"
measure= "coh"
freq=813. #alpha
#to see the example of the same plot but for 0.5-Hz resolution coherence measures, uncomment the following two lines:
#measure= "Coherence_scipy"
#freq=np.arange(8, 13.5, 0.5) #alpha range
#Convert names of required channels to indices in data
if measure in np.unique(feature_metadata.loc[feature_metadata['feat_group']=='coh_bins', 'measures']):
feat_idx = COH_bin_idx[ev][measure][:, freq_all_bin.index(freq)].astype(int)
av_chan = 0
elif measure in np.unique(feature_metadata.loc[feature_metadata['feat_group']=='coh_fq', 'measures']):
feat_idx = COH_idx[ev][measure][:, freq_all.index(freq[0]):freq_all.index(freq[-1])].astype(int)
av_chan = 1
else:
print("ERROR: Measure not in allowed connectivity measures")
#find matched indices for given drug
matched_idx_ovr = matched_indices_ovr[medication_names.index(medname)]
matched_idx_dn = matched_indices_drug_naive[medication_names.index(medname)]
#Create matrices of mean log PSD values and their SEM for each group, and p-values for mean psd (mean over channels)
#Note: it can take a while as new linear mixed model is calculated for each feature.
M_coh_all, M_sem_all, p_ovr_corr, p_dn_corr = get_meanfeaturevals_and_sem(feat_idx, zdata, metadata, matched_idx_ovr, matched_idx_dn, transform_log=0, average_channels=av_chan)
# %%
#Reshape the mean coherence matrices and p_values to n_channels x n_channels shape:
def reshape_matrices_chanelpairs(matrix_to_reshape, channel_names, channel_pairs):
channel_pairs_r = {v: k for k, v in channel_pairs.items()}
reshaped_matrix = np.zeros((np.shape(matrix_to_reshape)[0], len(channel_names), len(channel_names)))
for g in range(np.shape(matrix_to_reshape)[0]):
for ch1 in range(0, len(channel_names)):
for ch2 in range(ch1+1, len(channel_names)):
reshaped_matrix[g, ch1, ch2] = matrix_to_reshape[g, channel_pairs_r[(channel_names[ch1], channel_names[ch2])]]
reshaped_matrix[g, ch2, ch1] = matrix_to_reshape[g, channel_pairs_r[(channel_names[ch1], channel_names[ch2])]]
return reshaped_matrix
# %%
M_coh_all = reshape_matrices_chanelpairs(M_coh_all, channel_names, channel_pairs)
p_ovr_corr = reshape_matrices_chanelpairs(p_ovr_corr, channel_names, channel_pairs)
p_dn_corr = reshape_matrices_chanelpairs(p_dn_corr, channel_names, channel_pairs)
# %%
np.shape(M_coh_all), np.shape(p_ovr_corr)
# %%
#Plot three groups connectivity side by side
groups = [medname, "other", "drug-naive"]
fig, ax = plt.subplots(1,3, figsize=(18, 5))
for g in range(np.shape(M_coh_all)[0]):
plot_conn = plot_sensor_connectivity_topomap(example_raw_eeg.info, M_coh_all[g], cmap="Reds", cbar_label=None, ax=ax[g]) #, top_n=30, v_min=0.015, v_max=0.050, cmap='viridis',ax=None, save_path=None):
ax[g].text(0.5, 1.02, groups[g], horizontalalignment='center', verticalalignment='bottom', size=16, transform=ax[g].transAxes)
# %%
#Plot differences between the medicine group and comparison group - OVR
comparison2 = "other"
g2 = groups.index(comparison2)
M_coh_diff = M_coh_all[0] - M_coh_all[g2]
#remove insignificant connections
not_signif = np.sum(p_ovr_corr<alpha, axis=0) < min_significant_matches
M_coh_diff[not_signif] = 0
# %%
v_max = np.max(np.abs(M_coh_diff))
if measure in np.unique(feature_metadata.loc[feature_metadata['feat_group']=='coh_fq', 'measures']):
freq2 = str(freq[0]) + '-' + str(freq[-1]) + 'Hz'
else:
freq2 = str(freq) + 'Hz'
tit = measure + '_' + medname + '_vs_' + comparison2 + '_at_' + freq2 + '_during_' + ev
plot_sensor_connectivity_topomap(example_raw_eeg.info,
M_coh_diff,
cmap="RdBu_r", cbar_label=None,
v_min=-v_max, v_max=v_max,
top_n = min(30, np.sum(~not_signif)),
#save_path=os.path.join(savepath_psd, tit + '.png'))
save_path=None) #change to a path if you want the graph saved instead of showing

Brainwaves_Under_Medication_Visualisation_and_PCA_analysis_code_clean.ipynb at commit 99898aa, no license · at the source

Overview

Authors: Magdalena Szponar¹, Patrycja Dzianok^1,2, Bartłomiej Gmaj³, Wojciech Jernajczyk⁴, Jan Kamiński^5,1

ORCID iDs: Bartłomiej Gmaj

Institutions: Instytut Biologii Doświadczalnej im. Marcelego Nenckiego (Poland); International Institute of Molecular and Cell Biology (Poland); Medical University of Warsaw (Poland); Institute of Psychiatry and Neurology (Poland); SUNY Upstate Medical University (United States)

Journal: EBioMedicine, volume 130, article 106375

Dates: received 8 December 2025; accepted 23 June 2026; published online 9 July 2026; in print August 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1016/j.ebiom.2026.106375 · PMID 42424703 · PMCID PMC13380497 · OpenAlex W7167808942

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: EEG (modality), human (organism), clinical / translational (subfield)

Methods: Connectivity, Spectral & time-frequency, Statistics, Smoothing, state filtering, decompositions, Preprocessing, Complexity, Physiology & signal measures

Keywords: Psychotropic drugs, Psychiatry, Pharmacotherapy, Mental disorders, Electroencephalography (EEG), PharmacoEEG

MeSH: Brain*, Brain Waves*, Electroencephalography*, Mental Disorders*, Psychotropic Drugs*, Cross-Sectional Studies, Female, Humans, Male (* major topic)

Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Funding: European Funds for Smart Economy; This publication was financially supported by the Foundation for Polish Science (FNP) within the BRAINCITY IRAP project, funded by the EU under the European Funds for Smart Economy programme

Citations: cited by 1 paper (Europe PMC); 99 references in the paper

Abstract

Background: Psychotropic medications remain foundational in psychiatric care, yet the neurophysiological mechanisms through which they exert therapeutic and adverse effects are still poorly characterised, limiting the field's ability to optimise treatment selection and monitoring. Electroencephalography (EEG) offers a non-invasive, real-time window into brain function that could support more precise, mechanism-informed prescribing; however, progress has been constrained by the absence of sufficiently large and systematically analysed pharmaco-EEG datasets.

Methods: In this cross-sectional observational study, we analysed over 24,000 clinical EEG recordings (∼6000 h of data) obtained across a wide range of psychiatric diagnoses and medication regimens. We compared more than 75,000 spectral, connectivity, and nonlinear EEG features across major drug classes, including benzodiazepines, SSRIs, antipsychotics, and anticonvulsants.

Findings: Dimensionality-reduced analyses revealed robust, class-specific neurophysiological signatures that can be linked to psychotropic drugs' mechanisms of action: benzodiazepines increased beta and decreased theta–alpha power; SSRIs enhanced gamma-band coherence; and antipsychotics and anticonvulsants produced marked slow-wave amplification and reductions in signal complexity. All results are made publicly accessible through an interactive resource (BrainwavesRX), enabling clinicians and researchers to explore medication-specific EEG effects at multiple levels of granularity.

Interpretation: By establishing a population-level reference atlas of psychotropic medication effects on human neural dynamics, this study provides an important foundation for future studies leveraging EEG to predict treatment response, detect insufficient or excessive pharmacological effects, and ultimately advance the development of individualised, data-driven psychiatric care.

Funding: The publication was prepared as part of Foundation of Polish Science's Proof of Concept (FENG.02.01-IP.05-0010/24) and BRAINCITY IRAP (FENG.02.07-IP.05-0179/23) projects.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 8 matches between paragraphs and lines of code.

labianca/EEG-psychotropic-medications

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 99898aadd758379da6d9f23b98763479c731c631, 7 July 2026

Languages: MATLAB (5), Jupyter (1)

Size: 7 files, 6 scripts

Software Heritage: not archived

Found in: “Data sharing statement”

Holds: README, 1 notebook

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: EEGLAB (1 file), h5py (1 file), ICLabel (1 file), Statistics and Machine Learning Toolbox (1 file), Matplotlib (1 file), MNE-Python (1 file), NumPy (1 file), pandas (1 file), SciPy (1 file), seaborn (1 file), statsmodels (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

7 files

Brainwaves_Under_Medication_Visualisation_and_PCA_analysis_code_clean.ipynb, Jupyter, 1,338 lines, 3 matches
Brainwaves_Under_Medication_code_commented.m, MATLAB, 417 lines, 2 matches
balance_groups_meds.m, MATLAB, 346 lines
balance_groups_meds_DN.m, MATLAB, 350 lines, 1 match
bonf_holm.m, MATLAB, 37 lines
preprocessing_REST_ASR_flexible_commented.m, MATLAB, 652 lines, 2 matches
README.md, Text, 7 lines

At the source: github.com/labianca/EEG-psychotropic-medications

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
6 scripts, each with its path and the digest of its content;
8 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data sharing statement

The results from all comparisons between all drug classes are available on the interactive website, at https://brainwavesrx.nencki.edu.pl/.

The code used for data preprocessing and statistical analysis is available on GitHub at https://github.com/labianca/EEG-psychotropic-medications.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 5 authors, 6 keywords, 9 MeSH terms, 2 funders, 91 references.

Cite

This paper

Szponar, M., Dzianok, P., Gmaj, B., Jernajczyk, W., & Kamiński, J. (2026). Brainwaves under medication: revealing class-specific neural signatures of psychotropic medication from 24,000 EEGs. EBioMedicine, 130, 106375. https://doi.org/10.1016/j.ebiom.2026.106375

BibTeX

@article{szponar2026brainwaves, author = {Szponar, Magdalena and Dzianok, Patrycja and Gmaj, Bartłomiej and Jernajczyk, Wojciech and Kamiński, Jan}, title = {{Brainwaves under medication: revealing class-specific neural signatures of psychotropic medication from 24,000 EEGs}}, journal = {EBioMedicine}, year = {2026}, month = jul, volume = {130}, pages = {106375}, publisher = {Elsevier}, issn = {2352-3964}, doi = {10.1016/j.ebiom.2026.106375}, url = {https://doi.org/10.1016/j.ebiom.2026.106375}, pmid = {42424703}, pmcid = {PMC13380497} }

RIS

TY - JOUR AU - Szponar, Magdalena AU - Dzianok, Patrycja AU - Gmaj, Bartłomiej AU - Jernajczyk, Wojciech AU - Kamiński, Jan TI - Brainwaves under medication: revealing class-specific neural signatures of psychotropic medication from 24,000 EEGs T2 - EBioMedicine J2 - eBioMedicine PY - 2026 DA - 2026/07/09 VL - 130 SP - 106375 SN - 2352-3964 PB - Elsevier DO - 10.1016/j.ebiom.2026.106375 UR - https://doi.org/10.1016/j.ebiom.2026.106375 LA - en ER -

CSL-JSON

{ "id": "10.1016/j.ebiom.2026.106375", "type": "article-journal", "title": "Brainwaves under medication: revealing class-specific neural signatures of psychotropic medication from 24,000 EEGs", "container-title": "EBioMedicine", "author": [ { "family": "Szponar", "given": "Magdalena" }, { "family": "Dzianok", "given": "Patrycja" }, { "family": "Gmaj", "given": "Bartłomiej" }, { "family": "Jernajczyk", "given": "Wojciech" }, { "family": "Kamiński", "given": "Jan" } ], "container-title-short": "eBioMedicine", "volume": "130", "page": "106375", "DOI": "10.1016/j.ebiom.2026.106375", "PMID": "42424703", "PMCID": "PMC13380497", "ISSN": "2352-3964", "publisher": "Elsevier", "URL": "https://doi.org/10.1016/j.ebiom.2026.106375", "language": "en", "issued": { "date-parts": [ [ 2026, 7, 9 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.3389/fpsyt.2026.1737357 [code]: A computational pipeline for a neurotransmitter-centric analysis of the effects of psychiatric medication on EEG spectral power.
Journal: Frontiers in psychiatry
In common: EEGLAB, Statistics and Machine Learning Toolbox, EEG, 6 references
[2] doi:10.1007/s10548-026-01238-y [code]: Topographic Reorganization of EEG Complexity During Visual Mental Imagery: Insights from Lempel-Ziv Complexity in High-Density EEG.
Journal: Brain topography
In common: ICLabel, MNE-Python, statsmodels, 4 other tools, EEG, 2 references
[3] doi:10.1002/mds.70348 [code]: Electroencephalography-Based Clustering Reveals Robust Neurophysiological Subtypes in Parkinson's Disease.
Journal: Movement disorders : official journal of the Movement Disorder Society
In common: ICLabel, EEGLAB, MNE-Python, 6 other tools, EEG
[4] doi:10.7554/elife.107088 [code]: Development of auditory and spontaneous movement responses to music over the first postnatal year.
Journal: eLife
In common: ICLabel, EEGLAB, MNE-Python, 6 other tools, EEG
[5] doi:10.1162/imag.a.1229 [code]: 40 Hz audiovisual stimulation improves sustained attention and related brain oscillations.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: ICLabel, EEGLAB, Statistics and Machine Learning Toolbox, 5 other tools, EEG, 1 reference
[6] doi:10.2196/80286 [code]: At-Home Sleep Electroencephalography Assessment in Young and Older Adults Using a Novel Wireless Soft Electronics Sleep Monitoring System: Experimental Study.
Journal: JMIR formative research
In common: EEGLAB, MNE-Python, Statistics and Machine Learning Toolbox, 5 other tools, EEG, 1 reference
[7] doi:10.1162/imag.a.1245 [code]: Towards precision EEG connectomics: Evaluating the benefits of dense sampling.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: ICLabel, MNE-Python, statsmodels, 6 other tools, EEG
[8] doi:10.1162/imag.a.105 [code]: Right posterior theta reflects human parahippocampal phase resetting by salient cues during goal-directed navigation
Journal: n/a
In common: EEGLAB, MNE-Python, statsmodels, 6 other tools, EEG
[9] doi:10.1038/s41598-026-56070-y [code]: SSDLabeler: realistic semi-synthetic data generation for multi-label artifact classification in EEG.
Journal: Scientific reports
In common: ICLabel, EEGLAB, statsmodels, 5 other tools, EEG
[10] doi:10.3389/fncom.2026.1786996 [code]: Schumann-anchored golden ratio organization of human neural oscillations.
Journal: Frontiers in computational neuroscience
In common: MNE-Python, h5py, statsmodels, 5 other tools, EEG, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.