Unveiling the Distinctive Brain Functional Dynamics Between Parkinson's Disease and Progressive Supranuclear Palsy.

Code ↔ Paper

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

[1] § Methods › Characterization of Brain States ↔ hmmleida/clustering/_clustering.py, lines 985–1055 · score 0.79 · Davies Bouldin, state sequence, model fit, distortion, concatenated, decoding
[2] § Methods › Discriminative Power of Brain State Features ↔ hmmleida/plotting/_plotting.py, lines 707–789 · score 0.71 · static functional connectivity, Pearson correlation, static FC, rows, edge, matrices
[3] § Methods › Characterization of Brain States ↔ hmmleida/signal_tools/_signal_tools.py, lines 66–133 · score 0.68 · unit variance, BOLD signals, fMRI, zero, sensitivity
[4] § Methods › Discriminative Power of Brain State Features ↔ hmmleida/_data_load.py, lines 396–503 · score 0.62 · static functional connectivity, static FC, diagonal, Pearson, Class, matrices

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,055 lines · 35 KB · no license · 1 match

import pandas as pd
import numpy as np
import os
import matplotlib.pyplot as plt
import seaborn as sns
from hmmlearn import hmm
from scipy.spatial.distance import cdist
from sklearn.metrics.pairwise import manhattan_distances
from sklearn.metrics import (
silhouette_score,
davies_bouldin_score,
adjusted_mutual_info_score,
adjusted_rand_score
)
from sklearn.model_selection import (
train_test_split,
KFold,
StratifiedKFold
)
from sklearn.decomposition import PCA
from scipy.spatial import Voronoi,voronoi_plot_2d
import pickle
from ..data_utils.validation import _check_isint
class KMeansLeida():
"""
K-Means algorithm with cosine or euclidean
distance-based optimization.
Params:
-------
k : int.
The number of clusters to form as well as
the number of centroids to generate.
metric : str.
Whether to use 'cosine' or 'euclidean' distance
to optimize the centroids.
n_init : int.
Number of time the k-means algorithm will be
run with different centroid seeds. The final
results will be the best output of n_init
consecutive runs in terms of inertia.
n_iter : int.
Maximum number of iterations of the k-means
algorithm for a single run.
Attributes:
-----------
cluster_centers_ : ndarray of shape (n_centroids,n_ROIs).
Coordinates of cluster centers.
labels_ : ndarray of shape (n_samples,).
Labels of each point.
inertia_ : float.
Sum of squared distances of samples to their closest
cluster center (aka distortion).
"""
def __init__(self,k=2, metric='cosine',n_init=10,n_iter=1_000):
#validation of input data
if not isinstance(metric,str):
raise TypeError("'metric' must be a string!")
if metric not in ['euclidean','cosine']:
raise ValueError("'metric' must be either 'cosine' or 'euclidean'")
_check_isint({
'k':k,
'n_init':n_init,
'n_iter':n_iter
})
if k<2:
raise ValueError("'k' must be > 1")
self._k_ = k
self._metric_ = metric
self._n_iter_ = n_iter
self._n_init_ = n_init
def fit(self,y,random_state=None):
"""
Compute k-means clustering.
Params:
--------
y : ndarray of shape (n_samples,n_features).
Training instances to cluster.
random_state : int or None.
Determines random number generation
for centroid initialization. Use an
int to make the randomness deterministic.
Returns:
--------
self : object.
Fitted estimator.
"""
for init_idx in range(self._n_init_):
#print(f'\tCurrent initialization: {init_idx+1}')
if random_state is not None:
np.random.seed(random_state)
idx = np.random.choice(len(y), self._k_, replace=False)
#Step 1. Randomly assigning first centroids positions
centroids = y[idx, :]
#Step 2. Finding the distance between centroids and all the data labels
distances = cdist(y, centroids ,self._metric_)
#Step 3. Predict each observation label based on the minimum Distance
labels = np.array([np.argmin(i) for i in distances]) #Step 3
#Step 4.Repeating the above steps for a defined number of iterations
labels_all = []
for iter in range(self._n_iter_):
centroids = []
for centroid_idx in range(self._k_):
#Updating centroids by computing the mean of the observations in each cluster
temp_cent = y[labels==centroid_idx].mean(axis=0)
centroids.append(temp_cent)
centroids = np.vstack(centroids) #Updated centroids
distances = cdist(y, centroids ,self._metric_)
labels = np.array([np.argmin(i) for i in distances])
if iter!=0:
if np.array_equal(labels,labels_all[-1]):
#print(f'\tConverged at iteration n° {iter+1}')
break
labels_all.append(labels)
#compute distortion on final centroids of current centroids initialization
distortion = self._distortion(y,centroids)
#update centroids,distortion and labels keeping always the best updated
if init_idx == 0:
best_distortion = distortion
optimum_centroids = centroids
optimum_labels = labels.copy()
if distortion < best_distortion:
best_distortion = distortion.copy()
optimum_centroids = centroids.copy()
optimum_labels = labels.copy()
self.cluster_centers_ = optimum_centroids
self.labels_ = optimum_labels
self.inertia_ = best_distortion
self._remap_labels()
self._is_fitted = True
def _check_is_fitted(self):
"""
Check if the k-means model had been already fitted.
"""
if not hasattr(self,"_is_fitted"):
raise Exception("You have to fit the model first by using the 'fit' method.")
def transform(self,y,closest=False):
"""
Computes distances between each observation
or data point and each centroid. Differs from
sklearn transform method in that here we can select
if retrieve all the distances or only the distances
to closest centroid. In the new space, each dimension
is the distance to the cluster centers.
Params:
--------
y : ndarray of shape (n_samples,n_features).
Data to transform.
Returns:
--------
distances : ndarray of shape (n_samples,n_centroids).
y transformed in the new space.
"""
self._check_is_fitted()
distances = cdist(y, self.cluster_centers_ ,self._metric_)
if not closest:
return distances
else:
return distances.min(1)
def predict(self,y):
"""
Assign cluster label to each data point in 'y'.
Predict the closest cluster each sample in y belongs to.
Params:
-------
y : ndarray of shape (n_samples,n_features).
Data to predict.
Returns:
--------
labels : ndarray of shape (n_samples,).
Index of the cluster each sample belongs to.
"""
self._check_is_fitted()
labels = np.array([np.argmin(i) for i in self.transform(y)])
return labels
def _distortion(self,y,centroids):
"""
Sum of squared distances of samples to their
closest cluster center.
Params:
------
y : ndarray of shape (n_samples,n_features).
Samples clustered.
centroids : ndarray of shape (n_centroids,n_features).
Computed centroids/prototypes.
Returns:
--------
sum_sqr_dist : float.
Computed distortion.
"""
#distances = self.transform(y,closest=True)
distances = cdist(y, centroids ,self._metric_).min(1)
sqr_dist = distances**2
sum_sqr_dist = np.sum(sqr_dist)
return sum_sqr_dist
def _remap_labels(self):
""""
After fitting the k-means model and making the
labelling of each sample, this function computes
the frequency of each cluster label across all samples,
and relabels the original labels so that centroid 1 is
always the cluster with higher number of observations.
"""
label,freq = np.unique(self.labels_,return_counts=True)
sorted_idxs = np.argsort(freq)[::-1]
dct = {k:v for k,v in zip(sorted_idxs,label)}
new_labels = np.array([dct[n] for n in self.labels_])
self.labels_ = new_labels
self.cluster_centers_ = self.cluster_centers_[sorted_idxs,:]
def identify_states(eigens_dataset,K_min=2,K_max=20,n_init=15,random_state=None,plot=True,save_results=False,path=None):
"""
Perform k-means clustering for each value of 'k' to identify the
different (discrete) number of clusters (i.e.,the phase-locking
states). At eack 'k', a k-means model is fitted and a cluster is
assign to each eigenvector. These 'predictions' are located in a
new column of the provided 'eigens_dataset'.
In addition, each 'k' is evaluated by means of the Dunn score or
Dunn index, distortion (aka sum squared errors), silhouete score,
and Davies Bouldin score.
Params:
-------
eigens : ndarray with shape (n_eigevectors, n_ROIs).
Contains the extracted eigenvectors from the dFC matrices.
K_max, K_min : int.
Max and min number of clusters to fit.
n_init : int.
Number of times the clustering algorithm
will be run with different centroid seeds.
random_state : None or int.
Determines random number generation for
centroid initialization.
Use an int to make the randomness deterministic.
plot : bool.
Whether to create plots showing the number
of k in the X axis, and the Dunn mean score,
Silhoutte score, distortion, and Davis Bouldin
score for each k in the Y axis.
save_results : bool.
Whether to save results in '{path}/clustering'.
path : str.
Where to create the 'clustering' folder if
'save_results' was set to True.
Returns:
--------
eigens_predictions : pd.dataframe.
Contains the metadata (subjects_ids and conditions),
and a column for each K partition is added containing
the predicted label for each observation (row).
clustering_performance : pd.dataframe.
Contains the Dunn score, Silhouette score and distorion
values for each K explored.
models : dict.
Contain the fitted k-means models.
"""
if save_results:
try:
results_path = f'{path}/clustering'
if not os.path.exists(results_path):
os.makedirs(results_path)
print(f"-Creating folder to save results: './{results_path}'")
except:
print("Warning: the folder to save the results could't be created.")
#Execute the k-means models fitting process
eigens_predictions = eigens_dataset[['subject_id','condition']] #copy provided data (to avoid overwriting)
X = np.array(eigens_dataset.iloc[:,2:],dtype=np.float32) #keep array with the eigenvectors (remove 'subject_ids' and 'condition' columns)
N_samples = X.shape[0]
if N_samples>30_000:
random_samples_idx = np.random.choice(np.arange(N_samples),size=30_000,replace=False)
clustering_performance = []
models = {} #dict to save each fitted model
for k in range(K_min,K_max+1):
print(f'k = {k}')
kmeans = KMeansLeida(k=k,n_init=n_init,n_iter=1_000)
kmeans.fit(X,random_state=random_state)
models[f'k_{k}'] = kmeans
clustering_performance.append({
'k':k,
'Dunn_score':dunn_fast(X,kmeans.labels_) if N_samples<30_000
else dunn_fast(X[random_samples_idx,:],kmeans.labels_[random_samples_idx]),
'distortion':kmeans.inertia_,
'silhouette_score':silhouette_score(X=X, labels=kmeans.labels_, metric='cosine') if N_samples<30_000
else silhouette_score(X=X[random_samples_idx,:], labels=kmeans.labels_[random_samples_idx], metric='cosine'),
'Davies-Bouldin_score':davies_bouldin_score(X, kmeans.labels_)
})
#predict each volume label
eigens_predictions.loc[:,f'k_{k}'] = kmeans.labels_ #adding column with clusters assignment
#save model to disk
if save_results:
models_folder = f'{results_path}/models'
if not os.path.exists(models_folder):
os.mkdir(models_folder)
model_filename = f'model_k_{k}.pkl'
pickle.dump(kmeans, open(f'{models_folder}/{model_filename}', 'wb'))
clustering_performance = pd.DataFrame(clustering_performance)
if plot:
plt.ioff()
plot_clustering_scores(clustering_performance)
if save_results:
try:
plt.savefig(f'{results_path}/clustering_performance.png',dpi=300)
except:
print('The figures were not saved on local folder.')
plt.close()
#save results
if save_results:
try:
clustering_performance.to_csv(f'{results_path}/clustering_performance.csv',sep='\t',index=False)
eigens_predictions.to_csv(f'{results_path}/predictions.csv',sep='\t',index=False)
except:
print('The .csv files were not saved on disk. '
'You can still save them by yourself using the returned variables.')
print('\n*The clustering process has succesfully ended.')
return eigens_predictions,clustering_performance,models
def centroid2matrix(centroids,plot=True,cmap='jet',darkstyle=False):
"""
Take the controids resulting from the k-means
clustering (i.e., the phase-locking states) and
reconstruct the connectivity pattern in matrix
format.
Note: see Cabral et al. 2017 [p.4 and p.6].
Params:
-------
centroids : ndarray with shape (n_centroids,n_rois).
Clusters centroids of a specific K partition.
plot : bool.
Whether to create plot showing
the matrices.
cmap : str. Default = 'jet'.
Select the colormap to use.
darkstyle : bool.
Whether to use a black background.
Returns:
--------
states : ndarray with shape (n_rois,n_rois,n_centroids).
PL states in matrix format.
"""
if not isinstance(centroids,np.ndarray) or centroids.ndim!=2:
raise TypeError("'centroids' must be a 2D array (n_centroids,n_rois)!")
N_states, N_regions = centroids.shape
states = np.empty((N_regions,N_regions,N_states))
for state_idx in range(N_states):
#scale centroid by its maximum value and transpose the matrix
centroid = centroids[state_idx,:]/np.max(np.abs(centroids[state_idx,:]))
states[:,:,state_idx] = np.outer(centroid,centroid.T)
if plot:
plt.ion()
with plt.style.context('dark_background' if darkstyle else 'default'):
_, axs = plt.subplots(
nrows=1,
ncols=N_states,
figsize=(N_states*2, 2 if N_states<=10 else 1),
edgecolor='black',
subplot_kw=dict(box_aspect=1)
)
axs = axs.ravel()
ticks_size = 6 if N_states<=10 else 5 if 10<N_states<=15 else 4
for state_idx in range(N_states):
sns.heatmap(
states[:,:,state_idx],
ax=axs[state_idx],
vmax=1,
vmin=-1,
center=0,
cmap=cmap,
square=True,
cbar=False
)
axs[state_idx].set_title(f'PL state {state_idx+1}',fontsize=7 if N_states<15 else 5)
#axs[state_idx].set_xlabel(f'ROIs',fontsize=6.5)
#axs[state_idx].set_ylabel(f'ROIs',fontsize=6.5)
axs[state_idx].set_xticks(
np.arange(20,N_regions,20),
np.arange(20,N_regions,20).tolist(),
rotation=0
)
axs[state_idx].set_yticks(
np.arange(20,N_regions,20),
np.arange(20,N_regions,20).tolist(),
)
axs[state_idx].tick_params(
axis='both',
which='both',
bottom=False,
left=False,
top=False,
#labelsize=5 if N_states<15 else 3.5,
labelsize=ticks_size,
pad=0
#labelbottom=False,
#labelleft=False
)
plt.tight_layout(pad=1,w_pad=1)
plt.ioff()
plt.show()
return states
def patterns_stability(X,y=None,n_clusters=None,folds=5,metric='ari',plot=True,darkstyle=False):
"""
Run a stratified KFold cross-validation to explore
the stability of assigned clusters across folds.
The provided data in 'X' and is splitted in a train
set and a test set. Then, a cross-validation is
performed by splitting the train set into n number
of 'folds'. At each iteration, a given fold is used to
fit a KMeans model, which is used to assign a cluster
label to each observation of the test set. Finally,
the similarity between clusters assignments of the test
set across folds is evaluated by using either the ARI
(adjusted Rand index) score ('ari') or the AMI (adjusted
mutual information) score ('ami').
Params:
-------
X : ndarray with shape (n_volumes,n_rois).
Contains the eigenvectors.
y : ndarray with shape (n_volumes,). Optional.
Clusters assignment of each observation in 'X'.
If 'None', you must provide 'n_clusters'.
n_clusters : int. Optional.
Select the numbers of clusters to fit the
K-Means models. If 'y' is not None, then
the number of clusters is inferred from
the provided array.
folds : int.
Select the number of folds for the
cross-validation.
metric : str.
Metric to compute the similarity between
clusters assignments across folds.
plot : bool.
Whether to create a heatmap showing the
scores obtained in the cross-validation.
darkstyle : bool.
Whether to use a dark background for plotting.
Returns:
--------
scores : ndarray with shape (n_folds,n_folds).
Array containing the computed values of the
selected metric for each pair of folds.
"""
if not isinstance(metric,str) or metric not in ['ari','ami']:
raise Exception("'metric' must be either 'ari' or 'ami'")
if not isinstance(X,np.ndarray):
raise TypeError("'X' must be a 2D array!")
if y is None and n_clusters is None:
raise Exception("You must specify either predicted labels in 'y' "
"or a number of clusters ('n_clusters').")
elif y is not None and n_clusters is not None:
print("Warning: when 'y' is provided, the 'n_clusters' "
"is inferred from 'y'.")
if y is not None:
if not isinstance(y,np.ndarray):
raise TypeError("'y' must be an array!")
n_clusters = np.unique(y).size
X_tr,X_ts,y_tr,_ = train_test_split(X,y,stratify=y,train_size=.5)
else:
X_tr,X_ts = train_test_split(X,train_size=.5)
#compute centroids on different folds
kfold = StratifiedKFold(n_splits=folds) if y is not None else KFold(n_splits=folds)
predictions_ = np.empty((X_ts.shape[0],folds))
if y is not None:
for i,(tr_idx,_) in enumerate(kfold.split(X_tr,y_tr)):
print(f'Current fold: {i+1}')
km = KMeansLeida(k=n_clusters)
km.fit(X_tr[tr_idx])
predictions_[:,i] = km.predict(X_ts)
else:
for i,(tr_idx,_) in enumerate(kfold.split(X_tr)):
print(f'Current fold: {i+1}')
km = KMeansLeida(k=n_clusters)
km.fit(X_tr[tr_idx])
predictions_[:,i] = km.predict(X_ts)
#compute scores
scores = np.empty((folds,folds))
for fold_idx in range(folds):
for fold_idx2 in range(folds):
if metric=='ari':
scores[fold_idx,fold_idx2] = adjusted_rand_score(
predictions_[:,fold_idx],
predictions_[:,fold_idx2]
)
else:
scores[fold_idx,fold_idx2] = adjusted_mutual_info_score(
predictions_[:,fold_idx],
predictions_[:,fold_idx2]
)
font_sizes = {'folds':np.arange(2,21,1),'fsize':np.linspace(6,15,19)[::-1]}
if plot:
plt.ion()
with plt.style.context('dark_background' if darkstyle else 'default'):
matrix_ = np.triu(scores)
plt.figure()
sns.heatmap(
scores[1:,:-1],
vmin=0,
center=.5,
vmax=1,
cmap='viridis',
square=True,
linecolor='black' if darkstyle else 'white',
linewidths=.4,
annot=True,
annot_kws={'size':font_sizes['fsize'][folds-2]},
mask=matrix_[1:,:-1],
yticklabels=[f'{i+2}' for i in range(folds-1)],
xticklabels=[f'{i+1}' for i in range(folds-1)],
fmt='.2f',
cbar_kws={'label': 'ARI score' if metric=='ari' else 'AMI score','shrink': 0.5}
)
plt.xlabel('Fold',fontsize=16)
plt.ylabel('Fold',fontsize=16)
plt.yticks(rotation=0)
plt.tight_layout()
return scores
# Plotting functions
def plot_clustering_scores(data):
"""
Create a 2x2 panel with lineplots showing
the clustering evaluation metrics for each
explored k value (Dunn score, distortion,
silhouette score, and Davis-Bouldin score).
Params:
-------
data : pandas.dataframe.
Contain k in the 1st col, with the
rest of columns containing the values
of the clustering evaluation metrics.
"""
_,axs = plt.subplots(nrows=2,ncols=2,figsize=(11,6))
axs = np.ravel(axs)
for fig_idx,metric in enumerate(data.columns[1:].values):
sns.lineplot(x='k',y=metric,data=data,ax=axs[fig_idx],linewidth=3)
axs[fig_idx].set_xticks([i+2 for i in range(np.max(data.k)-1)])
axs[fig_idx].set_xlabel('Number of clusters',fontsize=15)
axs[fig_idx].set_ylabel(metric.replace('_',' '),fontsize=15)
axs[fig_idx].grid(False)
plt.tight_layout()
def barplot_eig(eig,features_list):
"""
Create barplot showing the values of the
eigenvector for each brain region. This
eigenvector could be either an eigenvector
of a particular time point, or a cluster
centroid.
Params:
-------
eig : ndarray with shape (N_ROIs,).
Contains the eigenvector to plot.
features_list : list.
Contains the names of the ROIs.
Must be in the same order as in 'eig'.
"""
# generating list of colors.
# Positive values get red, negative values get blue
cols = ['mediumblue' if i<0 else 'firebrick' for i in eig]
plt.ion()
plt.figure(figsize=(4,15))
sns.barplot(y=features_list,x=eig,orient='h',palette=cols)
plt.axvline(0,color='black')
if np.max(np.abs(eig))>0.15:
plt.xlim(-1.05,1.05)
else:
plt.xlim(-0.11,0.11)
plt.xticks(fontweight='regular')
plt.yticks(fontweight='regular',fontsize=8)
plt.ylabel('Brain areas',rotation='vertical',fontsize=24,labelpad=40)
plt.tight_layout()
def barplot_states(centroids,rois_labels,save_path=None):
"""
Create subplots with barplots showing
the values of each cluster centroid.
Params:
------
centroids : ndarray with shape (n_centroids,n_rois).
Contains the centroids.
rois_labels : list.
Contains the names of each ROI.
"""
if not isinstance(centroids,np.ndarray):
raise TypeError("'centroids' must be a numpy 2D array!")
if centroids.shape[1]!=len(rois_labels):
raise Exception("The number of brain regions in 'centroids' and 'rois_labels' must be the same!")
N_centroids = centroids.shape[0]
plt.ion()
_,axs = plt.subplots(
ncols=N_centroids,
nrows=1,
sharey=True,
figsize=(N_centroids if N_centroids>4 else 4,8)
)
axs = np.ravel(axs)
for fig_idx in range(N_centroids):
# generating list of colors.
# Positive values get red, negative values get blue
cols = ['mediumblue' if i<0 else 'firebrick' for i in centroids[fig_idx,:]]
sns.barplot(
y=rois_labels,
x=centroids[fig_idx,:],
orient='h',
palette=cols,
ax=axs[fig_idx]
)
axs[fig_idx].axvline(0,color='black')
if np.max(np.abs(centroids[fig_idx,:]))>0.15:
axs[fig_idx].set_xlim(-1.05,1.05)
else:
axs[fig_idx].set_xlim(-0.11,0.11)
#axs[fig_idx].set_xlim(-0.11,0.11)
axs[fig_idx].tick_params(axis='y',labelsize=5)
axs[fig_idx].tick_params(axis='x',labelsize=6 if N_centroids>5 else 3)
axs[fig_idx].set_title(
f'PL pattern {fig_idx+1}',
fontweight='regular',
fontsize=4,
pad=10
)
axs[0].set_ylabel('Brain areas',rotation='vertical',fontsize=24,labelpad=25)
plt.tight_layout()
# plt.ioff()
plt.show()
if save_path is not None:
plt.savefig(save_path)
print(f"Figure saved to: {save_path}")
plt.close()
def plot_clusters3D(eigens,labels,clusters_colors=None,grid=True,alpha=.7,dot_size=3,edgecolor=None):
"""
Visualize the identified clusters (brain states) in a
3D scatter plot, which constitutes a low-dimensional
representation of the 'state space'.
Method : take the eigenvectors and extract the first
three principal components to reduce the dimensionality
of the data to a 3D space. Each dot in the plot thus
represents a single eigenvector, and is coloured according
to the cluster it belongs to.
Params:
-------
eigens : ndarray with shape (n_samples, n_rois).
Contains the eigenvectors of each
subject for each time point.
labels : ndarray with shape (n_samples,).
Contains the predicted cluster assignment
for each eigenvector in 'eigens'.
clusters_colors : list [optional].
Provide a list with the desired color
of each cluster. If not provided, then
a predefined set of colors will be used.
grid : bool.
Whether to show grid or not.
alpha : float.
Set transparency of dots.
dot_size : float or int.
Set dot size.
edge_color : None | str.
Specify an edge color to use on dots.
Returns:
--------
fig : generated plot.
"""
if clusters_colors is not None:
if np.unique(labels).size!=len(clusters_colors):
raise ValueError('The number of clusters and the number of colours provided must match')
else:
clusters_colors = [
'royalblue',
'grey',
'tomato',
'orange',
'cyan',
'violet',
'yellow',
'purple',
'firebrick',
'teal',
'orchid',
'red',
'green',
'steelblue',
'indigo',
'gold',
'sienna',
'coral',
'olive',
'salmon'
]
pca = PCA(n_components=3)
pcs = pca.fit_transform(eigens)
x_pcs = pd.concat((pd.DataFrame(pcs),pd.Series(labels)),axis=1)
x_pcs.columns = np.array(['PC_1','PC_2','PC_3','Cluster'])
#plotting
fig = plt.figure(figsize=(8,8))
ax = plt.axes(projection ="3d")
for i,color in zip(np.unique(labels),clusters_colors):
ax.scatter3D(
x_pcs[x_pcs.Cluster==i]['PC_3'],
x_pcs[x_pcs.Cluster==i]['PC_2'],
x_pcs[x_pcs.Cluster==i]['PC_1'],
c=color,
edgecolors=edgecolor,
s=dot_size,
alpha=alpha
)
ax.w_xaxis.pane.fill = False
ax.w_yaxis.pane.fill = False
ax.w_zaxis.pane.fill = False
ax.set_xlabel('3rd PC', fontsize=15, fontweight='regular',labelpad=20)
ax.set_ylabel('2nd PC', fontsize=15, fontweight='regular',labelpad=20)
ax.set_zlabel('1st PC', fontsize=15, fontweight='regular',labelpad=20)
if not grid:
ax.grid(False)
plt.tight_layout()
plt.show()
return fig
def plot_voronoi(centroids):
"""
Plot the clusters centroids in a 2D Voronoi
cells space. Performs a PCA to reduce the
dimensionality of the original centroid space
to a 2D space.
Params:
--------
centroids : ndarray with shape (n_centroids, n_regions).
Centroids of a particular K partition.
"""
pcs = PCA(n_components=2).fit_transform(centroids)
vor = Voronoi(pcs)
plt.ion()
voronoi_plot_2d(vor,point_size=10,show_vertices=False)
plt.title(f'K = {centroids.shape[0]}')
plt.tight_layout()
def plot_clusters_boundaries(y,n_clusters=2,alpha=0.05):
"""
Plot cluster centroids decision boundaries in
2D space after applying PCA.
Params:
-------
y : ndarray with shape (n_samples, n_rois).
Contains the eigenvectors.
n_clusters : int.
Specify the k number of clusters.
alpha : float.
Specify the transparency of the
dots (observations).
"""
ỹ = PCA(n_components=2).fit_transform(y) #embedding
kmeans = KMeansLeida(k=n_clusters, n_init=4)
kmeans.fit(ỹ)
# Step size of the mesh. Decrease to increase the quality of the VQ.
h = 0.0007 # point in the mesh [x_min, x_max]x[y_min, y_max].
# Plot the decision boundary. For that, we will assign a color to each
x_min, x_max = ỹ[:, 1].min(), ỹ[:, 1].max()
y_min, y_max = ỹ[:, 0].min(), ỹ[:, 0].max()
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
# Obtain labels for each point in mesh. Use last trained model.
Z = kmeans.predict(np.c_[xx.ravel(), yy.ravel()])
# Put the result into a color plot
Z = Z.reshape(xx.shape)
plt.ion()
plt.figure(figsize=(8,8))
#plt.clf()
plt.imshow(
Z,
interpolation="nearest",
extent=(xx.min(), xx.max(), yy.min(), yy.max()),
cmap=plt.cm.Paired,
aspect="auto",
origin="lower",
alpha=.9
)
# Plot dots (eigenvectors/samples)
plt.plot(
ỹ[:, 1],
ỹ[:, 0],
'o',
markersize=4,
markerfacecolor='black',
markeredgecolor='grey',
markeredgewidth=.5,
alpha=alpha
)
# Plot the centroids as a white X
centroids = kmeans.cluster_centers_
plt.scatter(
centroids[:, 1],
centroids[:, 0],
marker="X",
s=400,
linewidths=3,
color="black",
zorder=10,
edgecolors='w'
)
plt.title(
f"K-means clustering (K={n_clusters}) on the eigenvectors (PCA-reduced data)"
)
plt.xlim(x_min, x_max)
plt.xlabel("Centroids are marked with white cross.\nObservations are marked with black dots")
plt.ylim(y_min, y_max)
plt.xticks(())
plt.yticks(())
plt.tight_layout()
# Compute Dunn Score functions
def dunn_fast(points, labels):
"""
Compute the Dunn index.
Params:
----------
points : ndarray with shape (N_samples,N_features).
Observations/samples.
labels : ndarray with shape (N_samples).
Labels of each observation in 'points'.
"""
distances = manhattan_distances(points)
ks = np.sort(np.unique(labels))
deltas = np.ones([len(ks), len(ks)])*1_000_000
big_deltas = np.zeros([len(ks), 1])
l_range = list(range(0, len(ks)))
for k in l_range:
for l in (l_range[0:k]+l_range[k+1:]):
deltas[k, l] = _delta_fast((labels == ks[k]), (labels == ks[l]), distances)
big_deltas[k] = _big_delta_fast((labels == ks[k]), distances)
di = np.min(deltas)/np.max(big_deltas)
return di
def _delta_fast(ck, cl, distances):
values = distances[np.where(ck)][:, np.where(cl)]
values = values[np.nonzero(values)]
return np.min(values)
def _big_delta_fast(ci, distances):
values = distances[np.where(ci)][:, np.where(ci)]
#values = values[np.nonzero(values)]
return np.max(values)
def hmm_analysis(concatenated_dataset,K_min=2,K_max=20,n_init=20,random_state=None,plot=True,save_results=False,path=None):
if save_results:
try:
results_path = f'{path}/clustering'
if not os.path.exists(results_path):
os.makedirs(results_path)
print(f"-Creating folder to save results: './{results_path}'")
except:
print("Warning: the folder to save the results could't be created.")
# Execute the k-means models fitting process
predictions = concatenated_dataset[['subject_id', 'condition']] # copy provided data (to avoid overwriting)
X = np.array(concatenated_dataset.iloc[:, 2:],
dtype=np.float32) # keep array with the eigenvectors (remove 'subject_ids' and 'condition' columns)
clustering_performance = []
models = {} # dict to save each fitted model
for k in range(K_min, K_max + 1):
print(f'k = {k}')
model = hmm.GaussianHMM(n_components=k, covariance_type='diag', n_iter=1000, tol=1e-6, random_state=seed)
model.fit(X)
models[f'k_{k}'] = model
_, state_sequences = model.decode(X)
means = model.means_
distances = cdist(X, means, metric='euclidean')
nearest_center_idx = np.argmin(distances, axis=1)
distance_sum = np.sum(distances[np.arange(len(X)), nearest_center_idx] ** 2)
clustering_performance.append({
'k': k,
'Dunn_score': dunn_fast(X, state_sequences),
'distortion': distance_sum,
'silhouette_score': silhouette_score(X=X, labels=state_sequences, metric='cosine'),
'Davies-Bouldin_score': davies_bouldin_score(X, state_sequences)
})
# predict each volume label
predictions.loc[:, f'k_{k}'] = state_sequences # adding column with clusters assignment
# save model to disk
if save_results:
models_folder = f'{results_path}/models'
if not os.path.exists(models_folder):
os.mkdir(models_folder)
model_filename = f'model_k_{k}.pkl'
pickle.dump(model, open(f'{models_folder}/{model_filename}', 'wb'))
clustering_performance = pd.DataFrame(clustering_performance)
if plot:
plt.ioff()
plot_clustering_scores(clustering_performance)
if save_results:
try:
plt.savefig(f'{results_path}/clustering_performance.png', dpi=300)
except:
print('The figures were not saved on local folder.')
plt.close()
# save results
if save_results:
try:
clustering_performance.to_csv(f'{results_path}/clustering_performance.csv', sep='\t', index=False)
predictions.to_csv(f'{results_path}/predictions.csv', sep='\t', index=False)
except:
print('The .csv files were not saved on disk. '
'You can still save them by yourself using the returned variables.')
print('\n*The clustering process has succesfully ended.')
return predictions, clustering_performance, models

_clustering.py at commit 09ba8f9, no license · at the source

Overview

Authors: Chenfei Ye¹, Jiahui Shi², Jiliang Fang³, Tao Wu^4,5, Piu Chan^6,7,8,9, Jinghong Ma¹⁰, Ting Ma¹

ORCID iDs: Chenfei Ye, Piu Chan, Ting Ma

School of Biomedical Engineering Harbin Institute of Technology (Shenzhen) Shenzhen China
Department of Electronic and Information Engineering Harbin Institute of Technology (Shenzhen) Shenzhen China
Department of Radiology Guang'anmen Hospital, China Academy of Chinese Medical Sciences Beijing China
Department of Neurology, China National Clinical Research Center for Neurological Diseases Beijing Tiantan Hospital, Capital Medical University Beijing China
Center for Movement Disorders, Department of Neurology Beijing Tiantan Hospital, Capital Medical University Beijing China
Department of Neurology and Neurobiology Xuanwu Hospital Capital Medical University Beijing China
Key Laboratory on Neurodegenerative Disorders of Ministry of Education Xuanwu Hospital Capital Medical University Beijing China
Key Laboratory on Parkinson's Disease of Beijing Xuanwu Hospital Capital Medical University Beijing China
Parkinson's Disease Center of Beijing Institute of Brain Disorders, Collaborative Innovation Center for Brain Disorders Xuanwu Hospital Capital Medical University Beijing China
Department of Neurology Xuanwu Hospital of Capital Medical University Beijing China

Institutions: Harbin Institute of Technology (China); China Academy of Chinese Medical Sciences (China); Guang’anmen Hospital (China); Capital Medical University (China); Beijing Tian Tan Hospital (China); Xuan Wu Hospital of the Capital Medical University (China); Beijing Institute for Brain Disorders (China)

Journal: CNS neuroscience & therapeutics, volume 32, issue 9, article e71147

Dates: received 13 January 2026; accepted 31 August 2026; published online 4 September 2026; in print September 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1002/cns.71147 · PMID 42698292 · PMCID PMC13545542 · OpenAlex W7208805461

Open access: gold, a free copy (OpenAlex)

Status: code verified

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

Methods: Connectivity, Statistics, Machine learning, fMRI & imaging

Keywords: hidden Markov models, network dynamics, non‐motor impairment, Parkinson's disease, progressive supranuclear palsy

MeSH: Brain*, Parkinson Disease*, Supranuclear Palsy, Progressive*, Aged, Brain Mapping, Female, Hidden Markov Models, Humans, Magnetic Resonance Imaging, Male, Middle Aged (* major topic)

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

Funding: National Key Research and Development Program of China (2025YFF0517803); National Natural Science Foundation of China (62276081); Shenzhen Science and Technology Innovation Program (GXWD20231129121139001, JCYJ20240813110522029, JCYJ20250604145427037); Basic and Applied Basic Research Foundation of Guangdong Province (2023A1515010792, 2023B1515120065); Guangdong S&T Programme (2025B0101130004)

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

Abstract

Background: Dynamic analysis of resting‐state fMRI (rs‐fMRI) offers a novel approach to differentiate Parkinson's disease (PD) from progressive supranuclear palsy (PSP) by capturing temporal features of brain network activity, which may shed light on the mechanisms underlying non‐motor symptoms.

Objectives: To characterize differences in brain dynamics between PD and PSP using dynamic brain metrics, evaluate their exploratory discriminative performance, and investigate associations with non‐motor symptoms.

Methods: Sixty‐nine healthy controls, 82 PD patients, and 29 PSP patients underwent standardized clinical assessment and rs‐fMRI. Hidden Markov models extracted temporal features including fraction occurrence (FO), dwell time, and transition probability. These metrics were used for group comparisons, correlation analyses, and machine learning classification.

Results: PD and PSP showed opposite trends in fraction occurrence of unimodal network‐dominant states. Compared to PD and controls, PSP exhibited prolonged duration in the dorsal attention and limbic network (DAN&LIM) state and increased occurrence in the dorsal attention and frontoparietal control network (DAN&FPCN) state. Reduced unimodal state occupancy correlated with cognitive decline in both groups, while increased DAN&LIM and unimodal persistence linked to worse mood and sleep disturbances in PSP. Machine learning with these metrics achieved moderate accuracy in differentiating PD from PSP.

Conclusions: Temporal features of brain network dynamics may provide candidate imaging markers for distinguishing PD and PSP while offering mechanistic insights into non‐motor symptomatology. However, their diagnostic applicability requires validation in independent external multicenter cohorts.

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

Repositories

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

netneurolab/neuromaps

License: CC-BY-NC-SA-4.0

State: the link answers, verified on 26 September 2026

Evidence: files inventoried

Commit: ffcc2e0f657943ce00a1b6a968396f32250e495c, 4 June 2026

Languages: Python (45), Shell (3)

Size: 90 files, 48 scripts

Software Heritage: archived

Found in: “Data Availability Statement”

Holds: README, license file, environment (Dockerfile, pyproject.toml, requirements.txt, setup.py, docs/requirements.txt), tests, continuous integration, documentation

Not found: CITATION.cff

Tools: neuromaps (29 files), NumPy (19 files), NiBabel (8 files), SciPy (7 files), Nilearn (5 files), scikit-learn (3 files), Matplotlib (2 files), BrainSMASH (1 file), BrainSpace (1 file), pandas (1 file)

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

26 September 2026: the link answers

50 files

docs/conf.py, Python, 126 lines
examples/plot_fetch_datasets.py, Python, 108 lines
examples/plot_spatial_nulls.py, Python, 97 lines
neuromaps/__init__.py, Python, 2 lines
neuromaps/_version.py, Python, 683 lines
neuromaps/caret.py, Python, 164 lines
neuromaps/civet.py, Python, 120 lines
neuromaps/datasets/__init__.py, Python, 17 lines
neuromaps/datasets/annotations.py, Python, 416 lines
neuromaps/datasets/atlases.py, Python, 364 lines
neuromaps/datasets/tests/__init__.py, Python, 1 line
neuromaps/datasets/tests/test_annotations.py, Python, 43 lines
neuromaps/datasets/tests/test_atlases.py, Python, 76 lines
neuromaps/datasets/tests/test_utils.py, Python, 78 lines
neuromaps/datasets/utils.py, Python, 602 lines
neuromaps/images.py, Python, 624 lines
neuromaps/nulls/__init__.py, Python, 11 lines
neuromaps/nulls/burt.py, Python, 190 lines
neuromaps/nulls/nulls.py, Python, 739 lines
neuromaps/nulls/spins.py, Python, 687 lines
neuromaps/nulls/tests/__init__.py, Python, 1 line
neuromaps/nulls/tests/test_burt.py, Python, 34 lines
neuromaps/nulls/tests/test_nulls.py, Python, 64 lines
neuromaps/nulls/tests/test_spins.py, Python, 57 lines
neuromaps/parcellate.py, Python, 230 lines
neuromaps/plotting.py, Python, 137 lines
neuromaps/points.py, Python, 420 lines
neuromaps/resampling.py, Python, 317 lines
neuromaps/stats.py, Python, 275 lines
neuromaps/tests/__init__.py, Python, 1 line
neuromaps/tests/conftest.py, Python, 20 lines
neuromaps/tests/test_caret.py, Python, 34 lines
neuromaps/tests/test_civet.py, Python, 22 lines
neuromaps/tests/test_images.py, Python, 161 lines
neuromaps/tests/test_parcellate.py, Python, 22 lines
neuromaps/tests/test_plotting.py, Python, 20 lines
neuromaps/tests/test_points.py, Python, 59 lines
neuromaps/tests/test_resampling.py, Python, 63 lines
neuromaps/tests/test_stats.py, Python, 51 lines
neuromaps/tests/test_transforms.py, Python, 112 lines
neuromaps/tests/test_utils.py, Python, 26 lines
neuromaps/transforms.py, Python, 621 lines
neuromaps/utils.py, Python, 152 lines
setup.py, Python, 7 lines
tools/install_dependencies.sh, Shell, 33 lines
tools/install_package.sh, Shell, 21 lines
tools/run_checks.sh, Shell, 22 lines
versioneer.py, Python, 2,277 lines
LICENSE, License, 437 lines
README.rst, Text, 90 lines

At the source: github.com/netneurolab/neuromaps

Shijh123456/PD_PSP

License: none: the authors keep all their rights

State: the link answers, verified on 26 September 2026

Evidence: files inventoried

Commit: 09ba8f96ac5c2945f7dcf34e4eb8c03183d5d9b3, 30 July 2026

Languages: Python (19)

Size: 55 files, 19 scripts

Software Heritage: not archived

Found in: “Data Availability Statement”

Holds: README

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

Tools: NumPy (12 files), Matplotlib (9 files), pandas (9 files), seaborn (9 files), SciPy (4 files), Nilearn (3 files), imageio (2 files), scikit-learn (2 files), NiBabel (1 file)

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

26 September 2026: the link answers

20 files

hmmleida/__init__.py, Python, 37 lines
hmmleida/_data_load.py, Python, 1,593 lines, 1 match
hmmleida/_hmm_leida.py, Python, 236 lines
hmmleida/clustering/__init__.py, Python, 34 lines
hmmleida/clustering/_clustering.py, Python, 1,055 lines, 1 match
hmmleida/clustering/rsnets_overlap.py, Python, 235 lines
hmmleida/data_utils/__init__.py, Python, 31 lines
hmmleida/data_utils/_data_utils.py, Python, 324 lines
hmmleida/data_utils/validation.py, Python, 69 lines
hmmleida/dynamics_metrics/__init__.py, Python, 26 lines
hmmleida/dynamics_metrics/_dynamics_metrics.py, Python, 686 lines
hmmleida/plotting/__init__.py, Python, 36 lines
hmmleida/plotting/_plotting.py, Python, 1,458 lines, 1 match
hmmleida/plotting/dev_surfer.py, Python, 221 lines
hmmleida/signal_tools/__init__.py, Python, 18 lines
hmmleida/signal_tools/_signal_tools.py, Python, 325 lines, 1 match
hmmleida/stats/__init__.py, Python, 21 lines
hmmleida/stats/_stats.py, Python, 467 lines
run.py, Python, 7 lines
README.md, Text, 9 lines

At the source: github.com/Shijh123456/PD_PSP

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:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
67 scripts, each with its path and the digest of its content;
4 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 Availability Statement

The data that support the findings of this study are available on request from the corresponding author. The data are not publicly available due to privacy or ethical restrictions. The codes for dynamic brain state analyses were conducted using the hmmlearn package (https://github.com/hmmlearn/hmmlearn) and modified from pyleida (https://psychomark.github.io/leida‐python/ (https://psychomark.github.io/leida-python/)). Codes for machine learning used the scikit‐learn package (https://github.com/scikit‐learn/scikit‐learn (https://github.com/scikit-learn/scikit-learn)). Sensitivity analyses used neuromaps (https://github.com/netneurolab/neuromaps). The complete codes are publicly available at https://github.com/Shijh123456/PD_PSP.

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, issue, pages, dates, 7 authors, 5 keywords, 11 MeSH terms, 5 funders, 38 references.

Cite

This paper

Ye, C., Shi, J., Fang, J., Wu, T., Chan, P., Ma, J., & Ma, T. (2026). Unveiling the Distinctive Brain Functional Dynamics Between Parkinson's Disease and Progressive Supranuclear Palsy. CNS neuroscience & therapeutics, 32(9), e71147. https://doi.org/10.1002/cns.71147

BibTeX

@article{ye2026unveiling, author = {Ye, Chenfei and Shi, Jiahui and Fang, Jiliang and Wu, Tao and Chan, Piu and Ma, Jinghong and Ma, Ting}, title = {{Unveiling the Distinctive Brain Functional Dynamics Between Parkinson's Disease and Progressive Supranuclear Palsy}}, journal = {CNS neuroscience \& therapeutics}, year = {2026}, month = sep, volume = {32}, number = {9}, pages = {e71147}, publisher = {Wiley}, issn = {1755-5930}, doi = {10.1002/cns.71147}, url = {https://doi.org/10.1002/cns.71147}, pmid = {42698292}, pmcid = {PMC13545542} }

RIS

TY - JOUR AU - Ye, Chenfei AU - Shi, Jiahui AU - Fang, Jiliang AU - Wu, Tao AU - Chan, Piu AU - Ma, Jinghong AU - Ma, Ting TI - Unveiling the Distinctive Brain Functional Dynamics Between Parkinson's Disease and Progressive Supranuclear Palsy T2 - CNS neuroscience & therapeutics J2 - CNS Neurosci Ther PY - 2026 DA - 2026/09/01 VL - 32 IS - 9 SP - e71147 SN - 1755-5930 PB - Wiley DO - 10.1002/cns.71147 UR - https://doi.org/10.1002/cns.71147 LA - en ER -

CSL-JSON

{ "id": "10.1002/cns.71147", "type": "article-journal", "title": "Unveiling the Distinctive Brain Functional Dynamics Between Parkinson's Disease and Progressive Supranuclear Palsy", "container-title": "CNS neuroscience & therapeutics", "author": [ { "family": "Ye", "given": "Chenfei" }, { "family": "Shi", "given": "Jiahui" }, { "family": "Fang", "given": "Jiliang" }, { "family": "Wu", "given": "Tao" }, { "family": "Chan", "given": "Piu" }, { "family": "Ma", "given": "Jinghong" }, { "family": "Ma", "given": "Ting" } ], "container-title-short": "CNS Neurosci Ther", "volume": "32", "issue": "9", "page": "e71147", "DOI": "10.1002/cns.71147", "PMID": "42698292", "PMCID": "PMC13545542", "ISSN": "1755-5930", "publisher": "Wiley", "URL": "https://doi.org/10.1002/cns.71147", "language": "en", "issued": { "date-parts": [ [ 2026, 9, 1 ] ] } }

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.