Dynamic neural states underpin motor symptom severity in Parkinson's disease: a longitudinal analysis of chronic cortico-subthalamic nucleus recordings.

Code ↔ Paper

9 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 9 matches

[1] § Methods › Modelling symptom severity using static spectral features ↔ eBioMed_postAnalysis_bradykinesia.py, lines 525–542 · score 0.77 · 20–35 Hz, 40–70 Hz, 70–100 Hz, delta alpha, low gamma, low beta
[2] § Methods › Modelling symptom severity using HMM derived spatio-spectral features ↔ eBioMed_postAnalysis_bradykinesia.py, lines 525–542 · score 0.77 · 20–35 Hz, 40–70 Hz, 70–100 Hz, delta alpha, low gamma, low beta
[3] § Methods › Modelling symptom severity using HMM derived temporal features ↔ eBioMed_postAnalysis_HMM_dyskinesia.py, lines 266–379 · score 0.66 · avoid multicollinearity, independent variable, temporal properties, fit, dyskinesia, GLM
[4] § Methods › Modelling symptom severity using HMM derived temporal features ↔ eBioMed_postAnalysis_HMM_tremor.py, lines 206–332 · score 0.66 · avoid multicollinearity, independent variable, temporal properties, fit, tremor, GLM
[5] § Methods › Hidden Markov Model fit ↔ SpectralUpdateHMM.py, lines 208–248 · score 0.57 · OSL dynamics toolbox, Python, dimensional, matrix, window, channel
[6] § Methods › Static spectral analysis ↔ eBioMed_DualEstimation_HMM.py, lines 604–658 · score 0.56 · 2–100 Hz, power spectral densities, window, PSDs, coherence, HMM
[7] § Results › Spectral fingerprints of Cortico-Subthalamic states ↔ eBioMed_postAnalysis_spectral_and_temporal_prop_plots_quantiles.py, lines 14–64 · score 0.56 · Hidden Markov Model, power spectral density, bradykinesia score, HMM state, wirelessly, Parkinson
[8] § Results › Spectral fingerprints of Cortico-Subthalamic states ↔ eBioMed_postAnalysis_spectral_and_temporal_prop_plots_quantiles.py, lines 14–64 · score 0.55 · hidden Markov model, power spectral density, Spectral properties, bradykinesia scores, standard error, shifts
[9] § Methods › Modelling symptom severity using static spectral features ↔ eBioMed_postAnalysis_HMM_dyskinesia.py, lines 439–451 · score 0.51 · delta alpha, static spectral, low gamma, low beta, high beta, high gamma

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

Python · 1,827 lines · 71 KB · AGPL-3.0 · 2 matches

# %% IMPORTS
import os
import pickle
import pandas as pd
import numpy as np
import statsmodels.api as sm
from datetime import datetime
from copy import deepcopy
import matplotlib.pyplot as plt
import seaborn as sns
from utils import list_drives_with_names, format_func_single
from SpectralUpdateHMM import SpectralUpdateHMM
from statsmodels.stats.multitest import multipletests
# %%
"""
Script documentation
Bradykinesia related analysis for the HMM model with dual estimated time windows
Author: Abhinav Sharma ([email hidden])
"""
# %% # ^ SET THE PATHS
# Autodetect the drive letter for the external hard drive with a specific name
use_external_drive = False
drive_name = 'salTnsr_X10'
drive_letter = None
drives = list_drives_with_names()
for drive in drives:
if drive[1] == drive_name:
drive_letter = drive[0]
drive_letter = drive_letter[0]
break
# %% # ^ CREATE A DATAFRAME
# *This is going to cycle through all participants and get the data
# *Walk through the folder, the folder can contain subfolders
# *We will continue to walk the folder until we find the pickle files
subjectID = ''
# side = 'left'
ipsi_contra = 'contra'
all_dicts = []
# Four state model
folder_path = r'C:\Oxford\data\neural\uscf_results\HMM_masterModelMonday_June_10_2024_08_34_15\dual_estimation'
static_spectral_folder = r'C:\Oxford\data\neural\uscf_results\HMM_masterModelMonday_June_10_2024_08_34_15\static_spectra'
ground_truth_states = 4
model_name = "HMM_masterModelMonday_June_10_2024_08_34_15"
if use_external_drive:
folder_path = f'{drive_letter}:' + folder_path[2:]
for root, dirs, files in os.walk(folder_path):
for file in files:
if file.endswith('.pkl') and ipsi_contra in file and (subjectID if 'RCS' in subjectID else '') in file:
file_path = os.path.join(root, file)
# Search for the corresponding static spectral file
static_spectral_file = file_path.replace('dual_estimation', 'static_spectra')
static_spectral_file = static_spectral_file.replace('dualProps', 'staticSpectra')
with open(static_spectral_file, 'rb') as f:
static_spectral_dict = pickle.load(f)
with open(file_path, 'rb') as f:
# try:
dicts_from_file = pickle.load(f)
# We will only select specific properties to test from
# the dictionaries
# Get only the file name from the path
file_name = os.path.basename(file_path)
print('Loaded file:', file_name)
# Get the static spectra file name
static_spectral_file_name = os.path.basename(static_spectral_file)
print('Loaded static spectral file:', static_spectral_file_name)
# Print some blank lines for separation
print('\n\n')
# Check if the number at the end of the spectral file name
# matches the number at the end of the dual estimation file name
if file_name.split('_')[-1].split('.')[0] != static_spectral_file_name.split('_')[-1].split('.')[0]:
print('File name mismatch:', file_name, static_spectral_file_name)
raise ValueError('File name mismatch')
for d, spect in zip(dicts_from_file, static_spectral_dict):
new_dict = {}
new_dict['med_status'] = d['medication_status']
# Get the static spectral properties
new_dict['static_coh'] = spect['coh']
new_dict['f'] = spect['f']
new_dict['static_psd'] = spect['psd']
# Get the spectral properties for th HMM
new_dict['HMM_coh'] = d['coh']
new_dict['HMM_psd'] = d['psd']
new_dict['HMM_f'] = d['f']
new_dict['dyskinesia_score'] = d['dyskinesia_y']
new_dict['tremor_score'] = d['Tremor_Score_y']
new_dict['bradykinesia_score'] = d['bradykinesia_y']
new_dict['bradykinesia_transformed'] = d['bradykinesia_transformed_y']
# normed scores
new_dict['norm_dyskinesia_score'] = d['norm_dyskinesia']
new_dict['norm_bradykinesia_score'] = d['norm_bradykinesia']
# status
new_dict['tremor_status'] = d['Tremor_y']
new_dict['behavior_status_dysk'] = d['behavior_status_dysk']
# temporal properties
new_dict['mean_lifetimes'] = d['mean_lifetimes']
new_dict['mean_intervals'] = d['mean_intervals']
new_dict['FO'] = d['fractional_occ']
new_dict['switching_rate'] = d['switching_rates']
# Ten minute medication status based on majority voting
new_dict['ten_minute_med_status'] = d['ten_minute_med_status']
# Sleep status
new_dict['sleep_status'] = d['sleep_status']
new_dict['ten_minute_sleep_status'] = d['ten_minute_sleep_status']
# Transition matrix
# new_dict['transition_matrix'] = d['transition_matrix']
d_one, _ = format_func_single(d['transition_matrix'], modelStates=ground_truth_states)
new_dict['transition_matrix'] = d_one[0]
# Check the shape of the transition matrix
# It should be a square matrix
# With the number of rows and columns equal to the number of ground truth states
if new_dict['transition_matrix'].shape[0] != ground_truth_states:
print('Transition matrix shape mismatch:', new_dict['transition_matrix'].shape)
raise ValueError('Transition matrix shape mismatch')
# subject ID
if subjectID:
new_dict['subjectID'] = subjectID
else:
new_dict['subjectID'] = file_name.split('_')[0]
new_dict['file_number'] = file_name.split('_')[-1].split('.')[0]
all_dicts.append(new_dict)
# except:
# print('Error loading file:', file_path)
# continue
# all_dicts.extend(dicts_from_file)
final_df = pd.DataFrame(all_dicts)
if final_df.empty:
print('No data found, maybe check your folder path')
# %% # ^ REARRANGE THE DATAFRAME
# & CREATE A NEW DATAFRAME
# Let us rearrange the dataframe above to expand the temporal properties into separate columns for each state
# This applies to the temporal properties only
# We will create a new dataframe with the temporal properties expanded
# We will have a column for each state
# We will keep the medication status and the scores as they are
clincial_df = []
for i_loc,row in final_df.iterrows():
row_dict = {}
row_dict['med_status'] = row['med_status']
row_dict['dyskinesia_score'] = row['dyskinesia_score']
row_dict['tremor_score'] = row['tremor_score']
row_dict['bradykinesia_score'] = row['bradykinesia_score']
row_dict['bradykinesia_transformed'] = row['bradykinesia_transformed']
row_dict['norm_dyskinesia_score'] = row['norm_dyskinesia_score']
row_dict['norm_bradykinesia_score'] = row['norm_bradykinesia_score']
row_dict['tremor_status'] = row['tremor_status']
row_dict['behavior_status_dysk'] = row['behavior_status_dysk']
# Static spectral properties
row_dict['static_coh'] = row['static_coh']
row_dict['static_psd'] = row['static_psd']
row_dict['f'] = row['f']
# HMM spectral properties
row_dict['HMM_coh'] = row['HMM_coh']
row_dict['HMM_psd'] = row['HMM_psd']
row_dict['HMM_f'] = row['HMM_f']
# Ten minute medication status
row_dict['ten_minute_med_status'] = row['ten_minute_med_status']
# Sleep status
row_dict['sleep_status'] = row['sleep_status']
row_dict['ten_minute_sleep_status'] = row['ten_minute_sleep_status']
# Subject ID
row_dict['subjectID'] = row['subjectID']
# File number
row_dict['file_number'] = row['file_number']
# Avoid starting from state number 0 as names
for i in range(ground_truth_states):
row_dict[f'mean_lifetimes_{i+1}'] = row['mean_lifetimes'][i]
row_dict[f'mean_intervals_{i+1}'] = row['mean_intervals'][i]
row_dict[f'FO_{i+1}'] = row['FO'][i]
row_dict[f'switching_rate_{i+1}'] = row['switching_rate'][i]
# Refashion the transition matrix
for i in range(ground_truth_states):
for j in range(ground_truth_states):
row_dict[f'transition_prob_{i+1}_to_{j+1}'] = row['transition_matrix'][i][j]
# Append the row to the list
clincial_df.append(row_dict)
# Create the final dataframe
final_clinical_df = pd.DataFrame(clincial_df)
# Delete redundant variables to free up memory
del clincial_df, all_dicts, final_df
# %% # ^ STATIC SPECTRA PLOTS
# ^ ---------------------------------STATIC SPECTRAL ANALYSIS--------------------------------
# ^ -----------------------------------------------------------------------------------------
# ^ -----------------------------------------------------------------------------------------
# ^ -----------------------------------------------------------------------------------------
# ^ -----------------------------------------------------------------------------------------
filtered_df = final_clinical_df[final_clinical_df['med_status'].isin(['med_on', 'med_off', 'med_sleep'])]
# filtered_df = final_clinical_df[final_clinical_df['med_status'].isin(['med_sleep'])]
# Re-arrange the static spectral properties to plot them
psd = filtered_df['static_psd'].tolist()
psd = np.array(psd)
med_status = filtered_df['med_status'].tolist()
med_status = np.array(med_status)
f = filtered_df['f'].tolist()
f = np.array(f)[0, :]
# Find indices corresponding to following frequencies: 2.0 and 45.0
lower_index = np.argmin(np.abs(f - 2.0))
upper_index = np.argmin(np.abs(f - 45.0))
# Slice psd according to the indices
psd = psd[:, :, lower_index:upper_index]
# Averaging the STNs together and the motor cortices together
psd1 = psd[:, 0:1, :].mean(axis=1)
psd2 = psd[:, 2:3, :].mean(axis=1)
# Concatenate the averaged PSDs
# insert a new dimension to the psd1 and psd2 arrays
psd1 = np.expand_dims(psd1, axis=1)
psd2 = np.expand_dims(psd2, axis=1)
# Concatenate the arrays
psd = np.concatenate((psd1, psd2), axis=1)
# Group PSDs based on medication status
med_on_psd = psd[med_status == 'med_on']
med_off_psd = psd[med_status == 'med_off']
# Calculate the mean and standard error across the window dimension for each group
med_on_psd_mean = med_on_psd.mean(axis=0)
med_off_psd_mean = med_off_psd.mean(axis=0)
med_on_psd_se = med_on_psd.std(axis=0) / np.sqrt(med_on_psd.shape[0])
med_off_psd_se = med_off_psd.std(axis=0) / np.sqrt(med_off_psd.shape[0])
# Sleep
med_sleep_psd = psd[med_status == 'med_sleep']
med_sleep_psd_mean = med_sleep_psd.mean(axis=0)
med_sleep_psd_se = med_sleep_psd.std(axis=0) / np.sqrt(med_sleep_psd.shape[0])
# Set font to Arial
plt.rcParams['font.family'] = 'Arial'
# Define colors and channel names
# colors = {'med_on': '0c5ac', 'med_off': 'cf560a'}
# channel_names = ['STN1', 'STN2', 'Motor Cortex 1', 'Motor Cortex 2']
channel_names = ['STN', 'Cortex']
# Create subplots
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 4))
axes = axes.flatten()
for channel in range(psd.shape[1]):
ax = axes[channel]
# Plot med_on PSD
ax.plot(f[lower_index:upper_index], med_on_psd_mean[channel, :], label='low bradykinesia', color='#cf560a', linewidth=2)
ax.fill_between(f[lower_index:upper_index],
med_on_psd_mean[channel, :] - med_on_psd_se[channel, :],
med_on_psd_mean[channel, :] + med_on_psd_se[channel, :],
color='#cf560a', alpha=0.3)
# Plot med_off PSD
ax.plot(f[lower_index:upper_index], med_off_psd_mean[channel, :], label='high bradykinesia', color='#0a69cf', linewidth=2, linestyle='--')
ax.fill_between(f[lower_index:upper_index],
med_off_psd_mean[channel, :] - med_off_psd_se[channel, :],
med_off_psd_mean[channel, :] + med_off_psd_se[channel, :],
color='#0a69cf', alpha=0.3)
# Plot med_sleep PSD
ax.plot(f[lower_index:upper_index], med_sleep_psd_mean[channel, :], label='sleep', color='#25be72', linewidth=2)
ax.fill_between(f[lower_index:upper_index],
med_sleep_psd_mean[channel, :] - med_sleep_psd_se[channel, :],
med_sleep_psd_mean[channel, :] + med_sleep_psd_se[channel, :],
color='#25be72', alpha=0.3)
# Set labels and title
ax.set_title(channel_names[channel], fontsize=16)
ax.set_xlabel('Frequency (Hz)', fontsize=16)
ax.set_ylabel('Power Spectral Density (PSD)', fontsize=16)
ax.set_yscale('log')
ax.legend(fontsize=10)
ax.grid(True, linestyle='--', alpha=0.6)
# Set xtick font and fontsize
for tick in ax.xaxis.get_major_ticks():
tick.label1.set_fontsize(16)
# Set ytick font and fontsize
for tick in ax.yaxis.get_major_ticks():
tick.label1.set_fontsize(16)
# Adjust layout
plt.tight_layout()
plt.show()
fig.savefig(r'C:\Oxford\software\Movement_disorders_revision\static_spectral_properties_avrg_with_sleep.png', dpi=600)
fig.savefig(r'C:\Oxford\software\Movement_disorders_revision\static_spectral_properties_avrg_with_sleep.svg', dpi=600)
# % Percentile based psds
filtered_df = final_clinical_df[(final_clinical_df['bradykinesia_transformed'] <= 80)]
# Define the percentiles
percentiles = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100]
f = filtered_df['f'].tolist()
f = np.array(f)[0, :]
# Find indices corresponding to following frequencies: 2.0 and 45.0
lower_index = np.argmin(np.abs(f - 2.0))
upper_index = np.argmin(np.abs(f - 45.0))
fig, ax = plt.subplots(figsize=(6, 6))
for i, percentile in enumerate(percentiles):
bradykinesia_percentile = np.percentile(filtered_df['bradykinesia_transformed'], percentile)
# Filter the dataframe based on the percentile
percentile_df = filtered_df[(filtered_df['bradykinesia_transformed'] <= bradykinesia_percentile)]
# Static spectra
psd = percentile_df['static_psd'].tolist()
psd = np.array(psd)
# Slice psd according to the indices
psd = psd[:, :, lower_index:upper_index]
# Averaging the STNs together and the motor cortices together
psd_stn = psd[:, 0:1, :].mean(axis=1)
psd_stn_mean = psd_stn.mean(axis=0)
# psd_ctx = psd[:, 2:3, :].mean(axis=1)
# psd_ctx_mean = psd_ctx.mean(axis=0)
# Plot the STN PSD
sns.lineplot(x=f[lower_index:upper_index], y=psd_stn_mean,
label = f'{percentile}th percentile',
ax=ax, linewidth=1)
# Set labels and title
ax.set_title('STN LOG', fontsize=16)
ax.set_xlabel('Frequency (Hz)', fontsize=16)
ax.set_ylabel('Power Spectral Density (PSD)', fontsize=16)
ax.set_yscale('log')
# ax.set_yscale('linear')
ax.grid(True, linestyle='--', alpha=0.6)
# Set xtick font and fontsize
for tick in ax.xaxis.get_major_ticks():
tick.label1.set_fontsize(16)
# Set ytick font and fontsize
for tick in ax.yaxis.get_major_ticks():
tick.label1.set_fontsize(16)
# Add percentile legend
percentile_labels = [f'{percentile}th percentile' for percentile in percentiles]
# Add the legend
# ax.legend(percentile_labels, fontsize=10)
# Display the legend outside the plot
# plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', fontsize=10)
# Adjust layout
plt.tight_layout()
plt.show()
# Save a very high resolution image
fig.savefig(r'C:\Oxford\software\Movement_disorders_revision\STN_LOG_static_spectral_properties_avrg_percentile.png', dpi=600)
# & Plot the static spectral properties
filtered_df = final_clinical_df[final_clinical_df['med_status'].isin(['med_on', 'med_off'])]
# We will have to re-arrange the static spectral properties to plot them
static_df = []
tsdf = []
window_number = 0
psd = filtered_df['static_psd'].tolist()
med_status = filtered_df['med_status'].tolist()
f = filtered_df['f'].tolist()
for i, p in enumerate(psd):
for j, p1 in enumerate(p):
for k, f1 in enumerate(f[j]):
temp_static_df = {}
temp_static_df['f'] = f1
temp_static_df['psd'] = p1
temp_static_df['med_status'] = med_status[i]
temp_static_df['window_number'] = i * len(p) + j + 1
temp_static_df['channel_number'] = k
tsdf.append(temp_static_df)
static_df = pd.DataFrame(tsdf)
# % Plotting the static spectral properties
fig, ax = plt.subplots(figsize=(50, 8))
sns.lineplot(data=static_df, x='f', y='psd', hue='med_status', ax=ax, errorbar=('ci', 95))
# All the x ticks should be shown
# There are 196 x ticks with unique values of f in the static_df being the labels
# We will show all the x ticks
ax.set_xticks(static_df['f'].unique())
ax.set_xticklabels(static_df['f'].unique(), rotation=90)
plt.show()
# We will create four subplots for the static spectral properties for each channel
fig, axs = plt.subplots(4, 1, figsize=(60, 40))
for i, ax in enumerate(axs.flatten()):
sns.lineplot(data=static_df[static_df['channel_number'] == i], x='f', y='psd', hue='med_status', ax=ax, errorbar=('ci', 95))
ax.set_title(f'Channel {i}')
ax.set_xticks(static_df['f'].unique())
ax.set_xticklabels(static_df['f'].unique(), rotation=90)
# ^ ---------------------------------END OF STATIC SPECTRAL ANALYSIS-------------------------
# ^ -----------------------------------------------------------------------------------------
# ^ -----------------------------------------------------------------------------------------
# ^ -----------------------------------------------------------------------------------------
# ^ -----------------------------------------------------------------------------------------
# %% # ^ LOAD THE PREVIOUSLY EXISTING DATA TO UPDATE THE SPECTRAL PROPERTIES
# & LOAD THE PREVIOUSLY SAVED DATA
# Load the final_clinical_df
model_name = 'HMM_masterModelMonday_June_10_2024_08_34_15'
final_clinical_df = pd.read_pickle(f'{model_name}_clinical_scores.pkl')
# %% # ^ FILTER THE DATA TO SAVE MEMORY
# & Filter data based on the medication status
# # Filter the final_clinical_df
# filtered_df = final_clinical_df[final_clinical_df['med_status'].isin(['med_on', 'med_off'])]
# del final_clinical_df
# & Filter data based on the sleep status (keep both awake and sleep data)
filtered_df = final_clinical_df[final_clinical_df['sleep_status'].isin(['awake', 'asleep'])]
del final_clinical_df
# %% # ^ CHANNEL GROUPING
""" Channels and spectral grouping for averaging and creating spectral band entries in the dataframe
"""
# These are bipolar channels
node_names = ['LFP1', 'LFP2', 'ECoG1', 'ECoG2']
node_identifiers = [0, 1, 2, 3]
# lfp1-ecog1, lfp1-ecog2, lfp2-ecog1, lfp2-ecog2 are all the same
# The resulting tuples for the above combinations are: (0, 2), (0, 3), (1, 2), (1, 3)
# These will all be averaged to get the coherence between the two channels: lfp-ecog
lfp_ecog = [(0, 2), (0, 3), (1, 2), (1, 3)]
# The two other unique combinations are: lfp1-lfp2, ecog1-ecog2
# The resulting tuples for the above combinations are: (0, 1), (2, 3)
lfp1_lfp2 = [(0, 1)]
ecog1_ecog2 = [(2, 3)]
# %% # ^ SPECTRAL BAND DEFINITIONS
delta_alpha_group = (2, 9)
low_beta_group = (10, 20) # 10 20
high_beta_group = (20, 35) # 20 35
low_gamma_group = (40, 70) # this is fine
high_gamma_group = (70, 100) # this is fine
groups = {'delta_alpha': delta_alpha_group,
'low_beta': low_beta_group,
'high_beta': high_beta_group,
'low_gamma': low_gamma_group,
'high_gamma': high_gamma_group}
coarse_beta_group = (12, 40)
coarse_gamma_group = (70, 100)
coarse_groups = {'coarse_beta': coarse_beta_group,
'coarse_gamma': coarse_gamma_group}
# %% # ^ UPDATE THE SPECTRAL PROPERTIES
# % Update the HMM coherene and reorganize the data
# Path to the directory containing the NNMF projectors
path_NNMF = r'C:\Oxford\software\wireless\results\HMM_masterModelMonday_June_10_2024_08_34_15\NNMF_factors'
# Instantiate the class
spectral_update = SpectralUpdateHMM(filtered_df,
node_names,
node_identifiers,
lfp_ecog,
lfp1_lfp2,
ecog1_ecog2,
groups,
path_NNMF)
# Check if the dataframe is correct
spectral_update.check_df()
spectral_update.update_coherence()
spectral_update.update_coherence_unprojected()
spectral_update.update_psd()
spectral_update.update_psd_unprojected()
del filtered_df
filtered_df = deepcopy(spectral_update.df)
# %%Now for every group find the upper and lower indices corresponding to the frequency vector
# and the group definitions above
xx = filtered_df['f'].iloc[0]
groups_indices = {}
for group_name, group in groups.items():
lower_index = np.argmin(np.abs(xx - group[0]))
upper_index = np.argmin(np.abs(xx - group[1]))
groups_indices[group_name] = (lower_index, upper_index)
coarse_groups_indices = {}
for group_name, group in coarse_groups.items():
lower_index = np.argmin(np.abs(xx - group[0]))
upper_index = np.argmin(np.abs(xx - group[1]))
coarse_groups_indices[group_name] = (lower_index, upper_index)
# %% # ^ STATIC SPECTRAL PROPERTIES UPDATE
"""
& PSD caclulations
^ The following groups are created
* 1) Refined spectral bands with all four channels stored within the same group
* 2) Coarse spectral bands with all four channels stored within the same group
"""
# Calculating average spectral band psd values for each group
for group_name, group in groups_indices.items():
filtered_df[f'{group_name}_psd'] = filtered_df['static_psd'].apply(lambda x: x[:,group[0]:group[1]].mean(axis=1))
for group_name, group in coarse_groups_indices.items():
filtered_df[f'{group_name}_psd'] = filtered_df['static_psd'].apply(lambda x: x[:,group[0]:group[1]].mean(axis=1))
"""
& PSD caclulations
^ The following groups are created
* 1) Refined spectral bands with four channels completely separated
* 2) Coarse spectral bands with four channels completely separated
* 3) Refined spectral bands with lfp and ecog channels separated (lfp1 and lfp2 are averaged, ecog1 and ecog2 are averaged)
* 4) Coarse spectral bands with lfp and ecog channels separated (lfp1 and lfp2 are averaged, ecog1 and ecog2 are averaged)
"""
for nm in node_names:
for group_name, group in groups_indices.items():
filtered_df[f'{nm}_{group_name}_psd'] = filtered_df['static_psd'].apply(lambda x: x[node_identifiers[node_names.index(nm)], group[0]:group[1]].mean())
for nm in node_names:
for group_name, group in coarse_groups_indices.items():
filtered_df[f'{nm}_{group_name}_psd'] = filtered_df['static_psd'].apply(lambda x: x[node_identifiers[node_names.index(nm)], group[0]:group[1]].mean())
# Let us also calculate the average psd for lfps and ecogs separately
for group_name, group in groups_indices.items():
filtered_df[f'lfp_{group_name}_psd'] = filtered_df['static_psd'].apply(lambda x: np.mean([x[0, group[0]:group[1]], x[1, group[0]:group[1]]]))
for group_name, group in coarse_groups_indices.items():
filtered_df[f'lfp_{group_name}_psd'] = filtered_df['static_psd'].apply(lambda x: np.mean([x[0, group[0]:group[1]], x[1, group[0]:group[1]]]))
for group_name, group in groups_indices.items():
filtered_df[f'ecog_{group_name}_psd'] = filtered_df['static_psd'].apply(lambda x: np.mean([x[2, group[0]:group[1]], x[3, group[0]:group[1]]]))
for group_name, group in coarse_groups_indices.items():
filtered_df[f'ecog_{group_name}_psd'] = filtered_df['static_psd'].apply(lambda x: np.mean([x[2, group[0]:group[1]], x[3, group[0]:group[1]]]))
# %% # ^ STATIC COHERENCE CALCULATIONS UPDATE
"""
& Coherence calculations
* Coherence is a matrix of shape (4, 4, n) where n is the number of frequency points
* We will calculate the average coherence for each entry in the matrix and
* average across based on the groups defined above for the frequency bands
* In the end for a spectral band we will have a 4x4 matrix with the average coherence values
"""
for group_name, group in groups_indices.items():
filtered_df[f'{group_name}_coh'] = filtered_df['static_coh'].apply(lambda x: x[:, :, group[0]:group[1]].mean(axis=2))
for group_name, group in coarse_groups_indices.items():
filtered_df[f'{group_name}_coh'] = filtered_df['static_coh'].apply(lambda x: x[:, :, group[0]:group[1]].mean(axis=2))
# & Coherence calculations
# * We will also calculate the spectral bands for lfp-lfp and ecog-ecog pairs
# ^ fine spectral bands
# For lfp-lfp pair
for group_name, group in groups_indices.items():
for nm in lfp1_lfp2:
filtered_df[f'lfp_lfp_{group_name}_coh'] = filtered_df['static_coh'].apply(lambda x: x[nm[0], nm[1], group[0]:group[1]].mean(axis=0))
# For the eco-ecog pair
for group_name, group in groups_indices.items():
for nm in ecog1_ecog2:
filtered_df[f'ecog_ecog_{group_name}_coh'] = filtered_df['static_coh'].apply(lambda x: x[nm[0], nm[1], group[0]:group[1]].mean(axis=0))
# ^ coarse spectral bands
# For lfp-lfp pair
for group_name, group in coarse_groups_indices.items():
for nm in lfp1_lfp2:
filtered_df[f'lfp_lfp_{group_name}_coh'] = filtered_df['static_coh'].apply(lambda x: x[nm[0], nm[1], group[0]:group[1]].mean(axis=0))
# For the eco-ecog pair
for group_name, group in coarse_groups_indices.items():
for nm in ecog1_ecog2:
filtered_df[f'ecog_ecog_{group_name}_coh'] = filtered_df['static_coh'].apply(lambda x: x[nm[0], nm[1], group[0]:group[1]].mean(axis=0))
# & Coherence calculations
# * Okay so now we will calculate spectral bands for lfp-ecog pairs and average across the pairs
# ^ fine spectral bands
for group_name, group in groups_indices.items():
# Add an empty column to the dataframe
filtered_df[f'lfp_ecog_{group_name}_coh'] = np.nan
for iloc, row in filtered_df.iterrows():
# Collect the data across all lfp-ecog pairs for the given group
ls = []
for nm in lfp_ecog:
val = row['static_coh'][nm[0], nm[1], group[0]:group[1]].mean(axis=0)
ls.append(val)
ls = np.array(ls)
assert ls.shape[0] == len(lfp_ecog)
filtered_df.at[iloc, f'lfp_ecog_{group_name}_coh'] = ls.mean()
# ^ coarse spectral bands
for group_name, group in coarse_groups_indices.items():
# Add an empty column to the dataframe
filtered_df[f'lfp_ecog_{group_name}_coh'] = np.nan
for iloc, row in filtered_df.iterrows():
# Collect the data across all lfp-ecog pairs for the given group
ls = []
for nm in lfp_ecog:
val = row['static_coh'][nm[0], nm[1], group[0]:group[1]].mean(axis=0)
ls.append(val)
ls = np.array(ls)
assert ls.shape[0] == len(lfp_ecog)
filtered_df.at[iloc, f'lfp_ecog_{group_name}_coh'] = ls.mean()
# %% # ^ SAVE THE UPDATED DATAFRAME
# Save the final med on med off dataframe to a pickle file
current_date = datetime.now().strftime('%A_%B_%d_%Y_%H_%M_%S')
# filtered_df.to_pickle(f'{model_name}_med_on_med_off_{current_date}_threshed.pkl')
filtered_df.to_pickle(f'{model_name}_awake_and_asleep_{current_date}_threshed.pkl')
# Save the final med on med off dataframe to a pickle file
# %% # ^ LOAD THE DATAFRAME
# This frame contains the temporal properties, clinical scores, and the spectral properties
# The spectral properties include both the static spectral properties and the HMM spectral properties
ground_truth_states = 4
model_name = "HMM_masterModelMonday_June_10_2024_08_34_15"
# Add date and time to file creation
# ? THESE ARE IN RECYCLE BIN
# load_date = 'Tuesday_December_10_2024_13_00_57'
# load_date = 'Monday_February_03_2025_11_10_48'
# load_date = 'Wednesday_February_12_2025_12_13_44'
# ? THIS IS THE LATEST AND CORRECT
load_date = 'Thursday_February_13_2025_22_16_39'
filtered_df = pd.read_pickle(f'{model_name}_med_on_med_off_{load_date}_threshed.pkl')
# %% # ^ Create a new column in the filtered_df where the bradykinesia_transformed is shifted as follows:
# ! CRITICAL:
# ^ This next bradykinesia score calculation and assignment has to be done before any other filtering
# ^ is applied. This is absolutely essential for correct temporal alignment.
# Get unique subject IDs
subject_ids = filtered_df['subjectID'].unique()
# The current row will have the next row's bradykinesia_transformed value and the last row will have no value
# This has to be performed within each subject and not across subjects
filtered_df = filtered_df.assign(next_bradykinesia_transformed=np.nan)
# Rearrange the dataframe columns so that bradykinesia_transformed and next_bradykinesia_transformed are next to each other
# All other columns should remain where they are
# Get the column names
cols = filtered_df.columns.tolist()
# Rearrange the columns
cols = cols[:cols.index('bradykinesia_transformed') + 1] + cols[-1:] + cols[cols.index('bradykinesia_transformed') + 1:-1]
# Now rearrange the columns
filtered_df = filtered_df[cols]
for subject in subject_ids:
# For this subject get the file numbers
file_numbers = filtered_df[filtered_df['subjectID'] == subject]['file_number'].unique()
# This ensures that the next bradykinesia scores are found within the same file number which is essentially the same session
for file_number in file_numbers:
# Get the indices of the subject and file number
indices = filtered_df[(filtered_df['subjectID'] == subject) & (filtered_df['file_number'] == file_number)].index
# Shift the values
filtered_df.loc[indices, 'next_bradykinesia_transformed'] = filtered_df.loc[indices, 'bradykinesia_transformed'].shift(-1)
# Last row should have no value assert this
assert np.isnan(filtered_df.loc[indices[-1], 'next_bradykinesia_transformed'])
# Drop the rows where the next_bradykinesia_transformed is nan
filtered_df = filtered_df.dropna(subset=['next_bradykinesia_transformed'])
# %% # ^ FILTER THE DATAFRAME
filtered_df = filtered_df[filtered_df['bradykinesia_transformed'] < 80]
# ! do not use this filter because med_on classification is not reliable
# filtered_df = filtered_df[filtered_df['med_status'].isin(['med_off', 'med_on'])]
# ! THIS IS THE MOST APPROPRIATE FILTER BECAUSE MED ON CLASSIFICATION IS NOT RELIABLE
filtered_df = filtered_df[filtered_df['med_status'].isin(['med_off'])]
# dyskinesia filter: behavior_status_dysk is either 'no-result' or 'non-dyskinetic'
filtered_df = filtered_df[filtered_df['behavior_status_dysk'].isin(['no-result', 'non-dyskinetic'])]
# Get unique subject IDs
subject_ids = filtered_df['subjectID'].unique()
# %%
# & STATISTICAL ANALYSIS
# ~ --------------------------------------------------------------------------------
# ^ ---------------------------------ANALYSIS------------------------------------- ^
# ~ --------------------------------------------------------------------------------
# %%
def benjamini_hochberg(p_values, alpha):
"""
Apply the Benjamini-Hochberg procedure to a list of p-values.
Parameters:
p_values (list or np.array): List or array of p-values to be corrected.
alpha (float): Desired false discovery rate (e.g., 0.05).
Returns:
np.array: Adjusted p-values.
"""
# Convert p_values to a numpy array for convenience
p_values = np.array(p_values)
m = len(p_values)
# Sort p-values and get the sorted indices
sorted_indices = np.argsort(p_values)
sorted_p_values = p_values[sorted_indices]
# Compute critical values
critical_values = (np.arange(1, m + 1) / m) * alpha
# Find the largest k where p-value <= critical value
below_critical = sorted_p_values <= critical_values
if np.any(below_critical):
k = np.max(np.where(below_critical)[0]) + 1
else:
k = 0
# Compute adjusted p-values
adjusted_p_values = np.minimum(1, (p_values * m) / (np.arange(1, m + 1)))
# Sort adjusted p-values according to their original order
final_p_values = np.ones_like(p_values)
if k > 0:
final_p_values[sorted_indices[:k]] = adjusted_p_values[sorted_indices[:k]]
return final_p_values
# %% # ^ GLM analysis for behavioral scores as dependent variables
def glm_analysis_med_status_full(df, clinical_score, temporal_property,
states,include_temporal_properties=True,
link_function='log',
glm_family = sm.families.Gamma(),
include_spectral_properties=False,
spectral_properties=None,
include_subjectID=False,
within_subject_centering=False,
subjectID=None,
scale_dependent_variable=False,
write_to_file=False, file_name=None):
results = {}
# Filter data for a specific subject if subject is not none
if subjectID and not include_subjectID:
df = df[df['subjectID'] == subjectID]
elif subjectID and include_subjectID:
raise ValueError('subjectID should not be included as a covariate when a single subject is being analyzed')
if include_temporal_properties is False and include_spectral_properties is False:
raise ValueError('At least one of include_temporal_properties or include_spectral_properties must be True')
if within_subject_centering and include_subjectID:
raise ValueError('Within subject centering is not possible with subjectID included as a covariate\n'
'Use only one of within_subject_centering or include_subjectID')
if include_temporal_properties:
# Directly create the X dataframe with the temporal property for all the states
X = df[[f'{temporal_property}_{s}' for s in states]]
# Center the temporal property within each subject
if within_subject_centering:
X = X - df.groupby('subjectID')[X.columns].transform('mean')
# Scale the temporal property between 0 and 1 for all states individually
for s in states:
X[f'{temporal_property}_{s}'] = (X[f'{temporal_property}_{s}'] - X[f'{temporal_property}_{s}'].min()) / (X[f'{temporal_property}_{s}'].max() - X[f'{temporal_property}_{s}'].min())
# Add static spectral properties as independent variables if needed
if include_spectral_properties:
if include_temporal_properties:
X = pd.concat([X, df[spectral_properties]], axis=1)
# For each static spectral property, scale between 0 and 1
for prop in spectral_properties:
if within_subject_centering:
X[prop] = X[prop] - df.groupby('subjectID')[prop].transform('mean')
# Global scaling for stable model fitting
X[prop] = (X[prop] - X[prop].min()) / (X[prop].max() - X[prop].min())
else:
X = df[spectral_properties]
for prop in spectral_properties:
if within_subject_centering:
X[prop] = X[prop] - df.groupby('subjectID')[prop].transform('mean')
# Global scaling for stable model fitting
X[prop] = (X[prop] - X[prop].min()) / (X[prop].max() - X[prop].min())
if include_subjectID:
# ! Always drop the first column to avoid multicollinearity !!!
X = pd.concat([X, pd.get_dummies(df['subjectID'], drop_first=True, dtype = 'int8')], axis=1)
# ! Always drop the first column to avoid multicollinearity !!!
# ! BUT DO NOT DROP ANY ROWS FROM X BASED ON NAN VALUES !!
# ! THERE SHOULD BE NO NAN VALUES IN X !!
y = df.loc[X.index, clinical_score].dropna()
# Scale the clinical score between 0 and 1
if scale_dependent_variable:
y = (y - y.min()) / (y.max() - y.min())
# Ensure X and y are aligned
X = X.loc[y.index]
# Add a constant term for the intercept
X = sm.add_constant(X)
# Ensure all columns in X are numeric
X = X.apply(pd.to_numeric, errors='coerce')
# Ensure y is numeric
y = pd.to_numeric(y, errors='coerce')
# Define the Tweedie family and GLM
if link_function == 'log':
link = sm.families.links.Log()
elif link_function == 'identity':
link = sm.families.links.Identity()
glm_model = sm.GLM(y, X, family=glm_family, link=link)
glm_results = glm_model.fit(cov_type='HC3', use_t=True)
# Predict the clinical score using the model
y_pred = glm_results.predict(X)
# Add subject ID to the predictions
y_pred = pd.concat([y_pred, df.loc[X.index, 'subjectID']], axis=1)
# Add the actual clinical score to the predictions
y_pred = pd.concat([y_pred, y], axis=1)
# Name the columns
y_pred.columns = ['Predicted', 'subjectID', 'Actual']
return results, glm_results, glm_model, y_pred
# return glm_results
# %% # ^ FEATURES TO INCLUDE IN THE ANALYSIS
# property_to_analyze = 'switching_rate'
# property_to_analyze = 'FO'
property_to_analyze = 'mean_lifetimes'
# property_to_analyze = 'mean_intervals'
# %%
gaussian_family = sm.families.Gaussian(link=sm.families.links.Identity())
# %% # ^ Independent variables 2
# ^ Only unprojects properties: stn and ecog psds and lfp-ecog coherence
STN_psd_props = []
STN_psd_props = ['state_1_delta_alpha_lfp_psd_unprojected', 'state_2_delta_alpha_lfp_psd_unprojected',
'state_3_delta_alpha_lfp_psd_unprojected', 'state_4_delta_alpha_lfp_psd_unprojected',
'state_1_low_beta_lfp_psd_unprojected', 'state_2_low_beta_lfp_psd_unprojected',
'state_3_low_beta_lfp_psd_unprojected', 'state_4_low_beta_lfp_psd_unprojected',
'state_1_high_beta_lfp_psd_unprojected', 'state_2_high_beta_lfp_psd_unprojected',
'state_3_high_beta_lfp_psd_unprojected', 'state_4_high_beta_lfp_psd_unprojected',
'state_1_low_gamma_lfp_psd_unprojected', 'state_2_low_gamma_lfp_psd_unprojected',
'state_3_low_gamma_lfp_psd_unprojected', 'state_4_low_gamma_lfp_psd_unprojected',
'state_1_high_gamma_lfp_psd_unprojected', 'state_2_high_gamma_lfp_psd_unprojected',
'state_3_high_gamma_lfp_psd_unprojected', 'state_4_high_gamma_lfp_psd_unprojected',]
ECOG_psd_props = []
ECOG_psd_props = ['state_1_delta_alpha_ecog_psd_unprojected', 'state_2_delta_alpha_ecog_psd_unprojected',
'state_3_delta_alpha_ecog_psd_unprojected', 'state_4_delta_alpha_ecog_psd_unprojected',
'state_1_low_beta_ecog_psd_unprojected', 'state_2_low_beta_ecog_psd_unprojected',
'state_3_low_beta_ecog_psd_unprojected', 'state_4_low_beta_ecog_psd_unprojected',
'state_1_high_beta_ecog_psd_unprojected', 'state_2_high_beta_ecog_psd_unprojected',
'state_3_high_beta_ecog_psd_unprojected', 'state_4_high_beta_ecog_psd_unprojected',
'state_1_low_gamma_ecog_psd_unprojected', 'state_2_low_gamma_ecog_psd_unprojected',
'state_3_low_gamma_ecog_psd_unprojected', 'state_4_low_gamma_ecog_psd_unprojected',
'state_1_high_gamma_ecog_psd_unprojected', 'state_2_high_gamma_ecog_psd_unprojected',
'state_3_high_gamma_ecog_psd_unprojected', 'state_4_high_gamma_ecog_psd_unprojected',]
LFP_ECOG_coherence_props = []
LFP_ECOG_coherence_props = ['state_1_delta_alpha_lfp_ecog_coherence_unprojected', 'state_2_delta_alpha_lfp_ecog_coherence_unprojected',
'state_3_delta_alpha_lfp_ecog_coherence_unprojected', 'state_4_delta_alpha_lfp_ecog_coherence_unprojected',
'state_1_low_beta_lfp_ecog_coherence_unprojected', 'state_2_low_beta_lfp_ecog_coherence_unprojected',
'state_3_low_beta_lfp_ecog_coherence_unprojected', 'state_4_low_beta_lfp_ecog_coherence_unprojected',
'state_1_high_beta_lfp_ecog_coherence_unprojected', 'state_2_high_beta_lfp_ecog_coherence_unprojected',
'state_3_high_beta_lfp_ecog_coherence_unprojected', 'state_4_high_beta_lfp_ecog_coherence_unprojected',
'state_1_low_gamma_lfp_ecog_coherence_unprojected', 'state_2_low_gamma_lfp_ecog_coherence_unprojected',
'state_3_low_gamma_lfp_ecog_coherence_unprojected', 'state_4_low_gamma_lfp_ecog_coherence_unprojected',
'state_1_high_gamma_lfp_ecog_coherence_unprojected', 'state_2_high_gamma_lfp_ecog_coherence_unprojected',
'state_3_high_gamma_lfp_ecog_coherence_unprojected', 'state_4_high_gamma_lfp_ecog_coherence_unprojected']
# Combine all the unprojected properties in a single list
final_spectral_props = (STN_psd_props + ECOG_psd_props + LFP_ECOG_coherence_props)
# %% # ^ Independent variables 3: Only static spectral props > local psds and coherence
static_spectral_properties = ['lfp_low_beta_psd', 'ecog_low_beta_psd',
'lfp_high_beta_psd', 'ecog_high_beta_psd',
'lfp_low_gamma_psd', 'ecog_low_gamma_psd',
'lfp_high_gamma_psd', 'ecog_high_gamma_psd',
'lfp_delta_alpha_psd', 'ecog_delta_alpha_psd']
static_lfp_ecog_coherence = ['lfp_ecog_delta_alpha_coh',
'lfp_ecog_low_beta_coh',
'lfp_ecog_high_beta_coh',
'lfp_ecog_low_gamma_coh',
'lfp_ecog_high_gamma_coh',]
final_spectral_props = (static_spectral_properties + static_lfp_ecog_coherence)
# %% # ^ Check if the final_spectral_props are in the dataframe
# Check if final_spectral_props are in the dataframe
for prop in final_spectral_props:
if prop not in filtered_df.columns:
print(f'{prop} not in the dataframe')
# %%
# & GLM Analysis
single_subjectID = False
if single_subjectID:
single_subjectID_name = 'RCS05L'
include_subjectID = False
else:
single_subjectID_name = None
include_subjectID = True
# ^ For FO, we have to use one state at a time due to the compositional nature of FO
# ^ Refer to compositional analysis for more details (hint: FO sums to 1)
include_temporal_properties = False
include_spectral_properties = True
if property_to_analyze == 'FO' and include_temporal_properties:
states = [1] # Change this to the state you want to analyze, repeat for each state
else:
states = [1, 2, 3, 4]
_, results, model, y_pred = glm_analysis_med_status_full(filtered_df,
'next_bradykinesia_transformed',
temporal_property=property_to_analyze,
states = states,
include_temporal_properties=include_temporal_properties,
glm_family=gaussian_family,
link_function='identity',
include_spectral_properties = include_spectral_properties,
include_subjectID=include_subjectID,
within_subject_centering=False,
subjectID=single_subjectID_name,
spectral_properties=final_spectral_props,
scale_dependent_variable=True,
write_to_file=False)
print(results.summary())
# Get the degrees of freedom for residuals , the model and the number of observations from the results
results.df_resid
results.df_model
results.nobs
# %% # ^ TEMPORAL PROPERTY SIMPLE OLS MODELS
dependent_variable_name = 'next_bradykinesia_transformed'
# independent_variable = ['FO_2', 'FO_4']
# independent_variable = [ 'mean_lifetimes_2', 'mean_lifetimes_4']
# independent_variable = ['mean_lifetimes_1', 'mean_lifetimes_2', 'mean_lifetimes_3' ,'mean_lifetimes_4']
# independent_variable = ['switching_rate_1', 'switching_rate_2', 'switching_rate_3', 'switching_rate_4']
# independent_variable = ['FO_1', 'FO_2', 'FO_3', 'FO_4']
independent_variable = ['mean_intervals_1', 'mean_intervals_2', 'mean_intervals_3', 'mean_intervals_4']
single_subjectID_name = 'RCS08R'
single_subjectID = True
# Simple model
subject_filtered_df = filtered_df[filtered_df['subjectID'] == single_subjectID_name]
independent_variable = subject_filtered_df[independent_variable]
# Add a constant term for the intercept
independent_variable = sm.add_constant(independent_variable)
dependent_variable = subject_filtered_df[dependent_variable_name]
# Scale the dependent variable between 0 and 1 # This is optional t-values are not affected by scaling
max_dep = np.max(dependent_variable)
min_dep = np.min(dependent_variable)
dependent_variable = (dependent_variable - min_dep)/(max_dep-min_dep)
# Fit a simple linear regression model
model = sm.OLS(dependent_variable, independent_variable)
results = model.fit(cov_type='HC3', use_t=True)
# display the results
print(results.summary())
# %% # ^ SPECTRAL PROPERTY
# Calculate the R2 value for the model based on the y_pred
# calculate the R2 value
R2 = 1 - np.sum((y_pred['Actual'] - y_pred['Predicted'])**2) / np.sum((y_pred['Actual'] - np.mean(y_pred['Actual']))**2)
print('The R2 value is: ' + str(R2))
# number of observations
n = len(y_pred)
# Get the number of independent variables
p = len(model.exog_names) - 1
# calculate the adjusted R2 value
adj_R2 = 1 - (1 - R2) * (n - 1) / (n - p - 1)
print('The adjusted R2 value is: ' + str(adj_R2))
# %% # ^ SPECTRAL PROPERTY Let us plot the predicted vs actual clinical score using the y pred dataframe
# Get the unique subject IDs
subject_ids = y_pred['subjectID'].unique()
# Cycle through the subject IDs and plot the predicted vs actual clinical score along with the best fit line
# We will use a unique color for each subject ID
# Assign a unique color to each subject ID
colors = sns.color_palette('husl', n_colors=len(subject_ids))
for i, subject_id in enumerate(subject_ids):
# Get the data for the current subject ID
curr_data = y_pred[y_pred['subjectID'] == subject_id]
# Create a new figure for each subject ID
plt.figure(figsize=(10, 6))
# Fit a best fit line to the plot
sns.regplot(data=curr_data, x='Actual',
y='Predicted', lowess=True, color=colors[i],
scatter_kws={'color': colors[i], 'alpha': 0.1, 's': 50},
line_kws={'color': 'black', 'linewidth': 2})
# Set the title for the plot
plt.title(f'Subject ID: {subject_id}')
# Show the plot
# plt.show()
# %% # ^ SPECTRAL PROPERTY DATAFRAMES
# & Extract the results
all_spectral_bands = ['high_gamma', 'low_gamma', 'high_beta', 'low_beta', 'delta_alpha']
all_channel_types_psd = ['lfp_psd', 'ecog_psd', 'lfp_ecog']
coef_df = pd.DataFrame({
"Coefficient": results.params,
"Standard Error": results.bse,
"t-value": results.tvalues,
"p-value": results.pvalues,
"Confidence Interval (0.025)": results.conf_int()[0],
"Confidence Interval (0.975)": results.conf_int()[1]
})
significant_results = deepcopy(coef_df[coef_df['p-value'] < 1.00])
significant_results = significant_results.assign(State=None, Spectral_Band=None, Channel_Type=None)
for index, row in significant_results.iterrows():
# Split the predictor name
predictor_name = index
# Split the predictor name
split_name = predictor_name.split('_')
# Check if the predictor name contains the state
if 'state' in split_name[0]:
state = split_name[0] + '_' + split_name[1]
# Check if the predictor name contains the spectral band
spectral_band = None
channel_type = None
if 'coh' in split_name[-1]:
spectral_band = split_name[2] + '_' + split_name[3]
channel_type = split_name[4] + '_' + split_name[5]
if 'coh' in split_name and 'unprojected' in split_name:
spectral_band = split_name[2] + '_' + split_name[3]
channel_type = split_name[4] + '_' + split_name[5] # dont add unprojected
if 'coherence' in split_name and 'unprojected' in split_name:
spectral_band = split_name[2] + '_' + split_name[3]
channel_type = split_name[4] + '_' + split_name[5]
if 'psd' in split_name[-1]:
spectral_band = split_name[2] + '_' + split_name[3]
channel_type = split_name[4] + '_' + split_name[5]
if 'psd' in split_name and 'unprojected' in split_name:
spectral_band = split_name[2] + '_' + split_name[3]
channel_type = split_name[4] + '_' + split_name[5] # dont add unprojected
else:
state = 'None'
spectral_band = 'None'
channel_type = 'None'
if 'lfp' in split_name and 'ecog' in split_name:
# If the predictor name contains lfp and ecog, we will set the channel type to lfp_ecog
channel_type = 'lfp_ecog'
# Remove lfp, ecog and coh/coherence from the split name
split_name = [name for name in split_name if name not in ['lfp', 'ecog', 'coh', 'coherence', 'psd', 'unprojected']]
# Now the join the remaining parts of the split name to get the spectral band
if len(split_name) > 1:
spectral_band = '_'.join(split_name[0:])
elif 'lfp' in split_name and 'psd' in split_name:
channel_type = 'lfp_psd'
# Remove lfp and psd from the split name
split_name = [name for name in split_name if name not in ['lfp', 'psd', 'unprojected']]
# Now the join the remaining parts of the split name to get the spectral band
if len(split_name) > 1:
spectral_band = '_'.join(split_name[0:])
elif 'ecog' in split_name and 'psd' in split_name:
channel_type = 'ecog_psd'
# Remove ecog and psd from the split name
split_name = [name for name in split_name if name not in ['ecog', 'psd', 'unprojected']]
# Now the join the remaining parts of the split name to get the spectral band
if len(split_name) > 1:
spectral_band = '_'.join(split_name[0:])
# Update the dataframe
significant_results.at[index, 'State'] = state
significant_results.at[index, 'Spectral Band'] = spectral_band
significant_results.at[index, 'Channel Type'] = channel_type
# Drop the following columns from the significant results
# significant_results.drop(columns=['Standard Error',
# 'Confidence Interval (0.025)',
# 'Confidence Interval (0.975)'],
# inplace=True)
# # Reorder the columns
significant_results = significant_results[[ 'State', 'Spectral Band', 'Channel Type',
'Coefficient','t-value', 'p-value', 'Standard Error',
'Confidence Interval (0.025)',
'Confidence Interval (0.975)']]
significant_results = significant_results[
(significant_results['Spectral Band'] != 'None') &
(significant_results['Channel Type'] != 'None')]
if single_subjectID:
# Add a column for the subject ID in the significant results dataframe for all rows
significant_results['Subject ID'] = single_subjectID_name
# ^ ULTRA CORRECTION: Correcting across states, spectral bands and channel types
# # Filter the results based on the corrected p-values
# significant_results_new = significant_results[significant_results['p-value_corr'] < 0.001]
if not single_subjectID:
# Correcting every p-value found
# Correct the p-values using fdr
p_values = significant_results['p-value']
p_values_corrected = multipletests(p_values, method='fdr_bh', alpha = 0.05)[1]
# Update the p-values in the dataframe
significant_results['p-value_corr'] = p_values_corrected
# Add a new column called: corrected_significance: True if the corrected p-value is less than 0.01
significant_results['corrected_significance'] = significant_results['p-value_corr'] < 0.01
# Add a new row at the end of the dataframe
significant_results.loc[len(significant_results)] = {
'State': 'DF_residuals',
'Spectral Band': results.df_resid,
'Channel Type': 'DF_model',
'Coefficient': results.df_model,
't-value': 'Num obs',
'p-value': results.nobs,
'Standard Error': 0,
'Confidence Interval (0.025)': 0,
'Confidence Interval (0.975)': 0
}
# Filter the results based on the corrected p-values
significant_results_ULTRA_CORRECTED = significant_results[significant_results['corrected_significance']]
# Add a new row at the end of the dataframe
significant_results_ULTRA_CORRECTED.loc[len(significant_results_ULTRA_CORRECTED)] = {
'State': 'DF_residuals',
'Spectral Band': results.df_resid,
'Channel Type': 'DF_model',
'Coefficient': results.df_model,
't-value': 'Num obs',
'p-value': results.nobs,
'Standard Error': 0,
'Confidence Interval (0.025)': 0,
'Confidence Interval (0.975)': 0
}
# %% # ^ SPECTRAL PROPERTIES ULTRA CORRECTION PLOTS
posthoc_ULTRA_glm_figure_dir = r'C:\Oxford\software\wireless\results\panel_figures\posthocs-ULTRA_NEWNEW'
spectral_band_colors = {
'delta_alpha': '#00A36C', # green
'low_beta': '#FFB336', # ligher yellow
'high_beta': '#E18500', # darker yellow
'low_gamma': '#15A8FF', # light blue
'high_gamma': '#004FB4' # dark blue
}
spectral_band_tuple_dict = {
'delta_alpha': (1.8, 1.9),
'low_beta': (1.4, 1.5),
'high_beta': (1.0, 1.1),
'low_gamma': (0.6, 0.7),
'high_gamma': (0.2, 0.3)
}
all_spectral_bands = ['high_gamma', 'low_gamma', 'high_beta', 'low_beta', 'delta_alpha']
all_spectral_band_names = ['high gamma', 'low gamma', 'high beta', 'low beta', 'delta/alpha']
bar_height = 0.2
bar_offset = 0.2 # Offset to ensure consistent spacing between positive and negative bars
left = 0
bar_alpha = 0.9
global_min = -15.0
global_max = 15.0
# Add stars based on siginificance level
# if p_values_corrected < 0.001:
# stars = '***'
# if p_values_corrected < 0.01:
# stars = '**'
# if p_values_corrected < 0.05:
# stars = '*'
for st in significant_results['State'].unique():
for ch in significant_results['Channel Type'].unique():
st_ch_data = significant_results[(significant_results['State'] == st) &
(significant_results['Channel Type'] == ch)]
fig, ax = plt.subplots(figsize=(6, 4))
ax.axvline(x=0, color='black', linewidth=1)
ax.set_xlim(global_min, global_max)
for i, row in st_ch_data.iterrows():
row_tval = row['t-value']
corrected_significance = row['corrected_significance']
if corrected_significance:
bar_color = spectral_band_colors[row['Spectral Band']]
else:
bar_color = 'darkgrey'
y_tuple = spectral_band_tuple_dict[row['Spectral Band']]
if row_tval > 0:
ax.barh(
y=y_tuple[0],
width=row_tval,
xerr=row['Standard Error'],
color=bar_color,
edgecolor='black',
left=left,
height=bar_height,
alpha=bar_alpha
)
ax.barh(
y=y_tuple[1],
width=-0.1,
xerr=None,
color='white',
edgecolor=bar_color,
left=left,
height=bar_height,
alpha=bar_alpha
)
if row_tval < 0:
ax.barh(
y=y_tuple[0],
width=0.1,
xerr=None,
color='white',
edgecolor=bar_color,
left=left,
height=bar_height,
alpha=bar_alpha
)
ax.barh(
y=y_tuple[1],
width=row_tval,
xerr=row['Standard Error'],
color=bar_color,
edgecolor='black',
left=left,
height=bar_height,
alpha=bar_alpha
)
# Add stars for significance
if row['p-value_corr'] < 0.001:
stars = '***'
elif row['p-value_corr'] < 0.01:
stars = '**'
elif row['p-value_corr'] < 0.05:
stars = '*'
else:
stars = ''
if stars:
ax.text(
x=row_tval + (0.5 if row_tval > 0 else -0.5), # Adjust position based on bar direction
y=y_tuple[0] + bar_height / 2,
s=stars,
fontsize=30,
ha='center',
va='center',
color='black'
)
# Set y tick labels as the spectral bands
ax.set_yticks([0.25, 0.65, 1.05, 1.45, 1.85])
ax.set_yticklabels(all_spectral_band_names, fontsize=16)
ax.set_xticks([-10, 0, 10])
ax.set_xlabel('t-value', fontsize=16, fontname='Arial')
ax.set_ylabel('Spectral Band', fontsize=16, fontname='Arial')
ax.tick_params(axis='x', labelsize=18)
ax.tick_params(axis='y', labelsize=18)
# Increase the size of the actual ticks and make them inward
ax.tick_params(axis='both', which='major', length=15, direction='in')
# Switch off the frame
ax.grid(True, axis='x', linestyle='-', alpha=1.0)
ax.grid(True, axis='y', linestyle='--', alpha=0.7)
# Set the title
ax.set_title(f'State {st} {ch}', fontsize=18, fontname='Arial', pad=20)
plt.show()
# Save the figure
fig.savefig(os.path.join(posthoc_ULTRA_glm_figure_dir, f'STARRED_ULTRA_{st}_{ch}.png'), dpi=600, bbox_inches='tight')
# %% # ^ TEMPORAL PROPERTY DATAFRAMES
coef_df = pd.DataFrame({
"Coefficient": results.params,
"Standard Error": results.bse,
"t-value": results.tvalues,
"p-value": results.pvalues,
"Confidence Interval (0.025)": results.conf_int()[0],
"Confidence Interval (0.975)": results.conf_int()[1]
})
significant_results = deepcopy(coef_df[coef_df['p-value'] < 1.00])
significant_results = significant_results.assign(State=None, Temporal_Property=None)
for index, row in significant_results.iterrows():
# Split the predictor name
predictor_name = index
# Split the predictor name
split_name = predictor_name.split('_')
if len(split_name) == 2:
significant_results.at[index, 'State'] = split_name[1]
significant_results.at[index, 'Temporal_Property'] = split_name[0]
elif len(split_name) == 3:
significant_results.at[index, 'State'] = split_name[2]
significant_results.at[index, 'Temporal_Property'] = split_name[0] + '_' + split_name[1]
else:
significant_results.at[index, 'State'] = 'None'
significant_results.at[index, 'Temporal_Property'] = 'None'
significant_results = significant_results[[ 'State', 'Temporal_Property',
'Coefficient','t-value', 'p-value', 'Standard Error',
'Confidence Interval (0.025)',
'Confidence Interval (0.975)']]
significant_results = significant_results[(significant_results['State'] != 'None')]
if single_subjectID:
# Add a column for the subject ID in the significant results dataframe for all rows
significant_results['Subject ID'] = single_subjectID_name
# ^ ULTRA CORRECTION: For a specific temporal property, correcting across all states
# ^ This is not applicable to Fractional Occupancy (FO): For FO correct manually
if not single_subjectID and len(states) > 1:
# Correcting every p-value found
# Correct the p-values using fdr
p_values = significant_results['p-value']
p_values_corrected = multipletests(p_values, method='fdr_tsbh', alpha = 0.00001)[1]
# Update the p-values in the dataframe
significant_results['p-value_corr'] = p_values_corrected
# Add a new column called: corrected_significance: True if the corrected p-value is less than 0.01
significant_results['corrected_significance'] = significant_results['p-value_corr'] < 0.01
if property_to_analyze == 'FO':
# Prompt the user to enter 4 p-values for the states
p_values = []
for state in range(1, 5):
while True:
try:
p_value = float(input(f"Enter the p-value for State {state}: "))
if 0 <= p_value <= 1:
p_values.append(p_value)
break
else:
print("Please enter a valid p-value between 0 and 1.")
except ValueError:
print("Invalid input. Please enter a numeric value.")
# Correct the p-values using the Benjamini-Hochberg procedure
alpha = 0.05
corrected_p_values = multipletests(p_values, method='fdr_bh', alpha=alpha)[1]
# Display the corrected p-values
for state, corrected_p in enumerate(corrected_p_values, start=1):
print(f"Corrected p-value for State {state}: {corrected_p}")
# %% # ^ ULTRA CORRECTION TEMPORAL PROPERTIES PLOTS
posthoc_ULTRA_glm_figure_dir = r'C:\Oxford\software\wireless\results\panel_figures\posthocs_temporal-ULTRA'
temporal_property = 'FO'
model_name = 'HMM_masterModelMonday_June_10_2024_08_34_15'
# Load the data
temporal_property_data_dir = r'C:\Oxford\software\wireless\results\HMM_masterModelMonday_June_10_2024_08_34_15\GLMs\Temporal_properties_May_2025'
file_name = f'None_{model_name}_ALLSTATES_MODEL_TEMPORAL_PROP_{temporal_property}_2025-05-14.csv'
# Load the data
temporal_property_data = pd.read_csv(os.path.join(temporal_property_data_dir, file_name))
# Get the unique states
unique_states = temporal_property_data['State'].unique()
# State colors = orange, blue, green, lilac
state_colors = {
1: '#FFB336', # orange
2: '#004FB4', # blue
3: '#00A36C', # green
4: '#E18500' # lilac
}
state_tuple_dict = {
1: (1.8, 1.9),
2: (1.4, 1.5),
3: (1.0, 1.1),
4: (0.6, 0.7)
}
bar_height = 0.2
bar_offset = 0.2 # Offset to ensure consistent spacing between positive and negative bars
left = 0
bar_alpha = 0.9
global_min = -6.0
global_max = 6.0
all_state_names = ['4', '3', '2', '1']
fig, ax = plt.subplots(figsize=(6, 4))
for st in temporal_property_data['State'].unique():
st_ch_data = temporal_property_data[(temporal_property_data['State'] == st)]
ax.axvline(x=0, color='black', linewidth=1)
ax.set_xlim(global_min, global_max)
for i, row in st_ch_data.iterrows():
row_tval = row['t-value']
corrected_significance = row['corrected_significance']
if corrected_significance:
bar_color = state_colors[row['State']]
else:
bar_color = 'darkgrey'
y_tuple = state_tuple_dict[row['State']]
if row_tval > 0:
ax.barh(
y=y_tuple[0],
width=row_tval,
xerr=row['Standard Error'],
color=bar_color,
edgecolor='black',
left=left,
height=bar_height,
alpha=bar_alpha
)
ax.barh(
y=y_tuple[1],
width=-0.1,
xerr=None,
color='white',
edgecolor=bar_color,
left=left,
height=bar_height,
alpha=bar_alpha
)
if row_tval < 0:
ax.barh(
y=y_tuple[0],
width=0.1,
xerr=None,
color='white',
edgecolor=bar_color,
left=left,
height=bar_height,
alpha=bar_alpha
)
ax.barh(
y=y_tuple[1],
width=row_tval,
xerr=row['Standard Error'],
color=bar_color,
edgecolor='black',
left=left,
height=bar_height,
alpha=bar_alpha
)
# Add stars for significance
# if row['p-value_corr'] < 0.001:
# stars = '***'
# elif row['p-value_corr'] < 0.01:
# stars = '**'
# elif row['p-value_corr'] < 0.05:
# stars = '*'
# else:
# stars = ''
# if stars:
# ax.text(
# x=row_tval + (0.5 if row_tval > 0 else -0.5) if abs(row_tval) > 1 else row_tval + (0.2 if row_tval > 0 else -0.2), # Adjust position based on bar length
# y=y_tuple[0] + bar_height / 2,
# s=stars,
# fontsize=30,
# ha='center',
# va='center',
# color='black'
# )
# Set y tick labels as the spectral bands
ax.set_yticks([0.65, 1.05, 1.45, 1.85])
ax.set_yticklabels(all_state_names, fontsize=16)
ax.set_xticks([-6, 0, 6])
ax.set_xlabel('t-value', fontsize=16, fontname='Arial')
ax.set_ylabel('States', fontsize=20, fontname='Arial')
ax.tick_params(axis='x', labelsize=18)
ax.tick_params(axis='y', labelsize=18)
# Increase the size of the actual ticks and make them inward
ax.tick_params(axis='both', which='major', length=15, direction='in')
# Switch off the frame
ax.grid(True, axis='x', linestyle='-', alpha=1.0)
ax.grid(True, axis='y', linestyle='--', alpha=0.7)
# Set the title
# Change temporal property case to first letter capital
temporal_property = temporal_property.capitalize()
ax.set_title(f'{temporal_property}', fontsize=18, fontname='Arial', pad=20)
plt.show()
# Save the figure
posthoc_ULTRA_temporals = r'C:\Oxford\software\wireless\results\panel_figures\posthocs_temporal-ULTRA'
fig.savefig(os.path.join(posthoc_ULTRA_temporals, f'ULTRA_{temporal_property}.png'), dpi=600, bbox_inches='tight')
# %% # ^ FOR BOTH SPECTRAL AND TEMPORAL PROPERTIES
# Save the results to a csv file drop the index
# significant_results.to_csv(f'significant_glm_results_bradykinesia_HMM_coh_spectra_without_med_control_within_subjects.csv', index=False)
# Add date of creation to the file name
# Get the current date
current_date = datetime.now().strftime('%Y-%m-%d')
# Get current time in hours, minutes and seconds am or pm
# current_time = datetime.now().strftime('%I-%M-%p')
# Add the date and time to the file name
# ^ PLEASE CHECK THE FILE NAME AND DIRECTORY
# ! DO NOT OVERWRITE THE FILE NAME
# GLM_result_directory = r'C:\Oxford\software\wireless\results\HMM_masterModelMonday_June_10_2024_08_34_15\GLMs\Temporal_properties_May_2025'
# Ask the user to enter the directory
GLM_result_directory = input('Enter the directory to save the results: ')
# Check if the directory exists, if not create it
if not os.path.exists(GLM_result_directory):
os.makedirs(GLM_result_directory)
if not(property_to_analyze == 'FO') and include_temporal_properties:
significant_results.to_csv(os.path.join(GLM_result_directory,
f'{single_subjectID_name}_{model_name}_ALLSTATES_MODEL_TEMPORAL_PROP_INTERVALS_{current_date}.csv'),
index=False)
if property_to_analyze == 'FO' and include_temporal_properties:
significant_results.to_csv(os.path.join(GLM_result_directory,
f'{single_subjectID_name}_{model_name}_STATE_{states[0]}_MODEL_TEMPORAL_PROP_FO_{current_date}.csv'),
index=False)
if not include_temporal_properties:
significant_results.to_csv(os.path.join(GLM_result_directory,
f'{single_subjectID_name}_{model_name}_ALLSTATES_MODEL_SPECTRAL_PROP_{current_date}.csv'),
index=False)
# If the variable significant_results_ULTRA_CORRECTED exists, save it to a csv file
if 'significant_results_ULTRA_CORRECTED' in locals():
significant_results_ULTRA_CORRECTED.to_csv(os.path.join(GLM_result_directory,
f'{single_subjectID_name}_{model_name}_ALLSTATES_MODEL_SPECTRAL_PROP_ULTRA_CORRECTED_{current_date}.csv'),
index=False)
# %% # ^ Collect model information and diagnostics for GLM
#
model_info = {
"Dependent Variable": results.model.endog_names,
"Model Type": type(model).__name__,
"Family": results.family.__class__.__name__,
"Link Function": results.family.link.__class__.__name__,
"Observations": results.nobs,
"Degrees of Freedom (Residual)": results.df_resid,
"Degrees of Freedom (Model)": results.df_model,
"Pearson Chi2": results.pearson_chi2,
"Scale": results.scale,
"Log-Likelihood": results.llf,
"Deviance": results.deviance,
"Null Deviance": results.null_deviance,
"AIC": results.aic,
"BIC": results.bic,
"Pseudo R-squared (CS)": results.pseudo_rsquared() # Pseudo R-squared (only applicable for certain GLM types)
}
# Convert to DataFrame
model_info_df = pd.DataFrame(list(model_info.items()), columns=["Metric", "Value"])
# %% # ^ Collect model information and diagnostics for GLM
# Perform the Wald test for the model
wald_test = results.wald_test_terms()
# Extract the datafraem
wald_test_df = wald_test.summary_frame()
# Only extract the dataframe for results where p-value is less than 0.01
significant_wald_test = wald_test_df[wald_test_df['P>F'] < 0.01]
#%% # ^ Let us segegrate the significant wald test based on states
# Will use this to only evaluate results if spectral properties are being studied
significant_spectral_wald_state = {}
for state in [1, 2, 3, 4]:
state_wald_test = significant_wald_test[significant_wald_test.index.str.contains(f'{state}')]
curr_state = {}
for row in state_wald_test.iterrows():
# Get the row index
index = row[0]
# Get the coefficient value
curr_state[index] = significant_results.loc[index].loc['Coefficient']
significant_spectral_wald_state[f'State {state}'] = curr_state
# %%

eBioMed_postAnalysis_bradykinesia.py at commit b4b7158, under AGPL-3.0 · at the source

Overview

Authors: Abhinav Sharma^1,2, Tao Liu^1,2, Bahman Abdi-Sargezeh^1,2, Amelia Hahn³, Maria Shcherbakova³, Wolf-Julian Neumann⁴, Simon Little³, Philip Starr³, Ashwini Oswal^1,2

ORCID iDs: Abhinav Sharma, Tao Liu, Maria Shcherbakova, Wolf-Julian Neumann, Ashwini Oswal

Institutions: University of Oxford (United Kingdom); MRC Brain Network Dynamics Unit (United Kingdom); University of California, San Francisco (United States); Stiftung Charité; Charité - Universitätsmedizin Berlin (Germany)

Journal: EBioMedicine, volume 128, article 106293

Dates: received 31 July 2025; accepted 28 April 2026; published online 12 May 2026; in print June 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1016/j.ebiom.2026.106293 · PMID 42119302 · PMCID PMC13191631 · OpenAlex W7161000153

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: human (organism), Parkinson's (population)

Methods: Spectral & time-frequency, Statistics, Connectivity, Machine learning, fMRI & imaging

Keywords: Neural states, Bradykinesia, Parkinson's disease

MeSH: Motor Cortex*, Parkinson Disease*, Subthalamic Nucleus*, Adult, Female, Humans, Hypokinesia, Longitudinal Studies, Male, Middle Aged, Severity of Illness Index, Tremor (* major topic)

Topic: Neurological disorders and treatments (Neurology, Medicine), according to OpenAlex

Funding: Rosetrees Trust; China Scholarship Council; Jon Moulton Charity Trust; Medical Research Council

Citations: not cited yet (Europe PMC); 64 references in the paper

Abstract

Background: Motor symptoms in Parkinson's disease (PD) may arise due to transient, network-wide neural dynamics that extend beyond beta-band oscillatory activity within the motor cortical-subthalamic nucleus (STN) circuit.

Methods: We applied a four-state Hidden Markov Model (HMM) to identify states of local and interregional oscillatory synchrony from chronic motor cortical and STN recordings (1046 h from 10 hemispheres) in five patients with PD (mean age 49 years), with concurrent measurements of bradykinesia, dyskinesia and tremor quantified using wearable sensors.

Findings: Neural states exhibited distinct spectral and temporal features relating to symptom severity. Two states exhibited spectral signatures—particularly STN low and high gamma, STN delta/alpha, cortical beta, and cortico-STN beta coherence—that predicted worsening bradykinesia. STN beta oscillations were not consistent predictors of bradykinesia (p = 0.52), but did predict worsening tremor (p < 0.01) and also improvements in dyskinesia severity (p < 0.001), in a state specific manner. These states also displayed compensatory features associated with bradykinesia amelioration, including cortical delta/alpha activity, cortical high gamma, and cortico-STN high gamma coherence. Additionally, we identified a state, marked by STN beta without cortico-STN beta coherence, whose increased lifetimes and occurrence improved motor function (p < 0.001).

Interpretation: Our findings highlight the multidimensional nature of motor impairments in PD and suggest that adaptive interventions targeting state features—rather than single frequency bands—may offer new opportunities for personalised neuromodulation.

Funding: AO: MRC Clinician Scientist Fellowship (MR/W024810/1), Rosetrees Trust/Race Against Dementia Team award, Oxford Hospitals Charity, and the Jon Moulton Charity Trust. TL: China Scholarship Council.

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 9 matches between paragraphs and lines of code.

saltwater-tensor/PD-neural-dynamics-motor-states

License: AGPL-3.0

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: b4b7158751cd6c6f5dc44bda2d8353066f90589f, 21 April 2026

Languages: Python (13)

Size: 14 files, 13 scripts

Software Heritage: not archived

Found in: “Data sharing statement”

Holds: license file

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

Tools: pandas (11 files), NumPy (10 files), Matplotlib (7 files), seaborn (6 files), statsmodels (4 files), SciPy (2 files), TensorFlow (2 files), scikit-posthocs (1 file)

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

28 September 2026: the link answers

14 files

BehaviorAnalysis.py, Python, 683 lines
BehaviorClassification.py, Python, 48 lines
CoherenceAnalyzer.py, Python, 382 lines
SpectralUpdateHMM.py, Python, 736 lines, 1 match
eBioMed_DualEstimation_HMM.py, Python, 870 lines, 1 match
eBioMed_Train_HMM.py, Python, 426 lines
eBioMed_postAnalysis_HMM_dyskinesia.py, Python, 1,401 lines, 2 matches
eBioMed_postAnalysis_HMM_tremor.py, Python, 1,361 lines, 1 match
eBioMed_postAnalysis_bradykinesia.py, Python, 1,827 lines, 2 matches
eBioMed_postAnalysis_dualEstim_HMM_sleep_plots.py, Python, 679 lines
eBioMed_postAnalysis_spectral_and_temporal_prop_plots_quantiles.py, Python, 1,434 lines, 2 matches
eBioMed_postAnalysis_subgrouping_data.py, Python, 163 lines
utils.py, Python, 287 lines
LICENSE, License, 661 lines

At the source: github.com/saltwater-tensor/PD-neural-dynamics-motor-states

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;
13 scripts, each with its path and the digest of its content;
9 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

Due to the sensitive nature of the data, patient-level datasets cannot be made publicly available. Fully anonymised data generated in this study will be made available to qualified researchers upon reasonable request, subject to a data sharing agreement, via the MRC CoRE in Restorative Neural Dynamics data platform (https://data.mrc.ox.ac.uk/).

Code for implementation of this manuscript is available at: https://github.com/saltwater-tensor/PD-neural-dynamics-motor-states.

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 2, 28 September 2026

Funding: added Rosetrees Trust; China Scholarship Council; Jon Moulton Charity Trust; Medical Research Council

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 9 authors, 3 keywords, 12 MeSH terms, 64 references.

Cite

This paper

Sharma, A., Liu, T., Abdi-Sargezeh, B., Hahn, A., Shcherbakova, M., Neumann, W.-J., Little, S., Starr, P., & Oswal, A. (2026). Dynamic neural states underpin motor symptom severity in Parkinson's disease: a longitudinal analysis of chronic cortico-subthalamic nucleus recordings. EBioMedicine, 128, 106293. https://doi.org/10.1016/j.ebiom.2026.106293

BibTeX

@article{sharma2026dynamic, author = {Sharma, Abhinav and Liu, Tao and Abdi-Sargezeh, Bahman and Hahn, Amelia and Shcherbakova, Maria and Neumann, Wolf-Julian and Little, Simon and Starr, Philip and Oswal, Ashwini}, title = {{Dynamic neural states underpin motor symptom severity in Parkinson's disease: a longitudinal analysis of chronic cortico-subthalamic nucleus recordings}}, journal = {EBioMedicine}, year = {2026}, month = may, volume = {128}, pages = {106293}, publisher = {Elsevier}, issn = {2352-3964}, doi = {10.1016/j.ebiom.2026.106293}, url = {https://doi.org/10.1016/j.ebiom.2026.106293}, pmid = {42119302}, pmcid = {PMC13191631} }

RIS

TY - JOUR AU - Sharma, Abhinav AU - Liu, Tao AU - Abdi-Sargezeh, Bahman AU - Hahn, Amelia AU - Shcherbakova, Maria AU - Neumann, Wolf-Julian AU - Little, Simon AU - Starr, Philip AU - Oswal, Ashwini TI - Dynamic neural states underpin motor symptom severity in Parkinson's disease: a longitudinal analysis of chronic cortico-subthalamic nucleus recordings T2 - EBioMedicine J2 - EBioMedicine PY - 2026 DA - 2026/05/12 VL - 128 SP - 106293 SN - 2352-3964 PB - Elsevier DO - 10.1016/j.ebiom.2026.106293 UR - https://doi.org/10.1016/j.ebiom.2026.106293 LA - en ER -

CSL-JSON

{ "id": "10.1016/j.ebiom.2026.106293", "type": "article-journal", "title": "Dynamic neural states underpin motor symptom severity in Parkinson's disease: a longitudinal analysis of chronic cortico-subthalamic nucleus recordings", "container-title": "EBioMedicine", "author": [ { "family": "Sharma", "given": "Abhinav" }, { "family": "Liu", "given": "Tao" }, { "family": "Abdi-Sargezeh", "given": "Bahman" }, { "family": "Hahn", "given": "Amelia" }, { "family": "Shcherbakova", "given": "Maria" }, { "family": "Neumann", "given": "Wolf-Julian" }, { "family": "Little", "given": "Simon" }, { "family": "Starr", "given": "Philip" }, { "family": "Oswal", "given": "Ashwini" } ], "container-title-short": "EBioMedicine", "volume": "128", "page": "106293", "DOI": "10.1016/j.ebiom.2026.106293", "PMID": "42119302", "PMCID": "PMC13191631", "ISSN": "2352-3964", "publisher": "Elsevier", "URL": "https://doi.org/10.1016/j.ebiom.2026.106293", "language": "en", "issued": { "date-parts": [ [ 2026, 5, 12 ] ] } }

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

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.