Divergent disruption of brain networks following total and chronic sleep loss: a longitudinal fMRI study.

Code ↔ Paper

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

[1] § Materials and Methods › Time series extraction ↔ Scripts/quality_check_FD_code.ipynb, lines 228–338 · score 0.68 · framewise displacement, Quality check, BOLD signal, Power, ROI, brain
[2] § Materials and Methods › Subjective data analysis ↔ Scripts/subjective_behavior_state_traits.ipynb, lines 1–142 · score 0.67 · subjective sleepiness, global graph metrics, global efficiency, graph distance, modularity, mixed
[3] § Materials and Methods › Subjective data analysis ↔ Scripts/subjective_behavior_state_traits.ipynb, lines 1–142 · score 0.63 · subjective behavioral, subjective sleepiness, objective graph, trait, baseline
[4] § Materials and Methods › Statistical analyses › Global and nodal metrics comparisons ↔ Scripts/Graphs_nodal_global_metrics_HDI_CCML.ipynb, lines 714–810 · score 0.62 · unpaired permutation, global metrics, robustness, shuffled, LMMs, FDR
[5] § Materials and Methods › Covariate-constraint manifold learning ↔ Scripts/Graphs_nodal_global_metrics_HDI_CCML.ipynb, lines 2931–3051 · score 0.62 · classical ISOMAP, global metrics, embeddings, manifold, CCML, covariates
[6] § Results › Nodal graph metrics ↔ Scripts/Graphs_nodal_global_metrics_HDI_CCML.ipynb, lines 3399–3440 · score 0.60 · III VI, Thalamus, Vermis, VII, FPN, limbic
[7] § Materials and Methods › Subjective data analysis ↔ Scripts/Graphs_nodal_global_metrics_HDI_CCML.ipynb, lines 527–586 · score 0.58 · linear mixed, global efficiency, graph distance, modularity, model, clustering
[8] § Materials and Methods › Subjective data analysis ↔ Scripts/subjective_behavior_state_traits.ipynb, lines 145–255 · score 0.58 · global graph metric, HC3, scored, OLS, trait, PSQI
[9] § Results › Nodal graph metrics ↔ Scripts/Graphs_nodal_global_metrics_HDI_CCML.ipynb, lines 3399–3440 · score 0.57 · III VI, Heschl, Precuneus, parietal, Angular, DMN
[10] § Materials and Methods › Time series extraction ↔ Scripts/quality_check_FD_code.ipynb, lines 228–338 · score 0.56 · temporal derivatives, motion parameters, outlier, signal, zero, regression

Paper

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

The paper is loaded when this pane is shown.

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

The authors' code

Jupyter notebook · 3,559 lines · 126 KB · no license · 5 matches

# %%
import networkx as nx
import matplotlib.pyplot as plt
import string
import numpy as np
from networkx.algorithms import community
import random
import copy
import scipy as sp
from scipy.stats import ttest_ind
from scipy.stats import mannwhitneyu
import statsmodels.formula.api as smf
from scipy.stats import ttest_rel
from scipy.stats import ttest_1samp
import statsmodels.api as sm
from statsmodels.formula.api import mixedlm
from statsmodels.stats.anova import AnovaRM
from statsmodels.stats.multitest import multipletests
from itertools import combinations
import statsmodels.stats.multitest as smm
from statsmodels.stats.multitest import fdrcorrection
import scipy.optimize as opt
import sklearn as sk
from sklearn import manifold, datasets
from sklearn.utils import resample
from sklearn.linear_model import LinearRegression
from scipy.interpolate import griddata
from scipy.spatial.distance import cdist
from matplotlib import ticker
import pandas as pd
import seaborn as sns
from matplotlib.lines import Line2D
import matplotlib.cm as cm
from netplotbrain import plot
from multipy.fwer import bonferroni
from multipy.fdr import lsu
import os
import glob
import re
import warnings
# %% [markdown]
# ## In this code, we will use the adjacency matrix and calculate nodal graph metrics (e.g., closeness, degree, betweenness).
# %%
#load coordinates once - the coordinates are always the same, because they are obtained directly from the atlas
coord_path = "/path/to/coordinates/from/AAL/atlas.txt"
coord = pd.read_csv(coord_path, sep="\s+", header=None).to_numpy()
pos = {k: (coord[k, 1], coord[k, 2]) for k in range(89)}
#loop through all sessions and subjects
base_path = "/path/to/mai/folder/with/data"
sessions = ["ses-1", "ses-2", "ses-3"]
for ses in sessions:
ses_path = os.path.join(base_path, ses)
subjects = glob.glob(os.path.join(ses_path, "sub-*"))
for sub_path in subjects:
sub = os.path.basename(sub_path) # e.g., "sub-01"
graph_path = os.path.join(sub_path, "Graphs/wAALours")
#only match files ending with "_400.txt" — this ensures we use adjacency matrices containing 400 edges.
#you can modify the number (e.g., "_200.txt", "_600.txt") to explore matrices with a different number of edges.
adj_files = glob.glob(os.path.join(graph_path, "*_400.txt"))
for adj_file in adj_files:
try:
#load adjacency matrix
df = pd.read_csv(adj_file, sep="\s+", header=None)
X = df.to_numpy()
G = nx.from_numpy_array(X)
#plot
plt.figure(figsize=(8, 8))
nx.draw(G, pos, node_size=50, with_labels=False)
#output directory
metrics_dir = os.path.join(graph_path, "graph_metrics")
os.makedirs(metrics_dir, exist_ok=True)
#file names follow the pattern "ses-1_sub-01_Adj_mat_..._400.txt".
#including "ses" (session) and "sub" (subject) helps verify that each file corresponds to the correct participant and session.
file_id = f"{ses}_{sub}_{os.path.splitext(os.path.basename(adj_file))[0]}"
#paths
image_path = os.path.join(metrics_dir, f"{file_id}.png")
metrics_path = os.path.join(metrics_dir, f"{file_id}_metrics.csv")
plt.savefig(image_path)
plt.close()
#compute metrics using NetworkX
closeness = nx.closeness_centrality(G)
betweenness = nx.betweenness_centrality(G)
clustering = nx.clustering(G)
degree = nx.degree_centrality(G)
#save to csv
metrics_df = pd.DataFrame({
"Node": list(closeness.keys()),
"Closeness": list(closeness.values()),
"Betweenness": list(betweenness.values()),
"Clustering": list(clustering.values()),
"Degree": list(degree.values())
})
metrics_df.to_csv(metrics_path, index=False)
print(f"Processed: {file_id}")
except Exception as e:
print(f"Error processing {adj_file}: {e}")
# %% [markdown]
# ## In this code, we will use the adjacency matrix and calculate global graph metrics (e.g., global efficiency).
# %% [markdown]
# GLOBAL EFFICIENCY
# %%
#sessions
session_labels = {
"ses-1": "baseline",
"ses-2": "acute",
"ses-3": "chronic"
}
base_path = "/path/to/mai/folder/with/data"
#result list
efficiency_results = []
#loop through each session and subject
for ses_id in session_labels.keys():
ses_path = os.path.join(base_path, ses_id)
subject_paths = glob.glob(os.path.join(ses_path, "sub-*"))
for sub_path in subject_paths:
sub_id = os.path.basename(sub_path)
adj_path = os.path.join(sub_path, "Graphs/wAALours", "*cost_400.txt")
adj_files = glob.glob(adj_path)
for adj_file in adj_files:
try:
df = pd.read_csv(adj_file, sep="\s+", header=None)
X = df.to_numpy()
G = nx.from_numpy_array(X)
efficiency = nx.global_efficiency(G)
efficiency_results.append({
"subject_id": sub_id,
"session_id": ses_id,
"global_efficiency": efficiency
})
print(f"{sub_id} | {ses_id} | {efficiency:.4f}")
except Exception as e:
print(f"Error processing {adj_file}: {e}")
#convert to df
df_out = pd.DataFrame(efficiency_results)
#pivot the dataframe so that each session becomes a column, subject_id forms the rows,
#session_id serves as the column names, and global_efficiency values fill the table.
glob_eff = df_out.pivot(index="subject_id", columns="session_id", values="global_efficiency").reset_index()
glob_eff.columns.name = None # remove column index name, e.g. "session_id""
#save to CSV
output_file = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/global_efficiency_all_subjects.csv"
os.makedirs(os.path.dirname(output_file), exist_ok=True)
glob_eff.to_csv(output_file, index=False)
print(f"Final results saved to: {output_file}")
# %% [markdown]
# AVERAGE CLUSTERING
# %%
#sessions
session_labels = {
"ses-1": "baseline",
"ses-2": "acute",
"ses-3": "chronic"
}
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
clustering_results = []
#loop through each session and subject
for ses_id in session_labels.keys():
ses_path = os.path.join(base_path, ses_id)
subject_paths = glob.glob(os.path.join(ses_path, "sub-*"))
for sub_path in subject_paths:
sub_id = os.path.basename(sub_path)
adj_path = os.path.join(sub_path, "Graphs/wAALours", "*cost_400.txt")
adj_files = glob.glob(adj_path)
for adj_file in adj_files:
try:
df = pd.read_csv(adj_file, sep="\s+", header=None)
X = df.to_numpy()
G = nx.from_numpy_array(X)
avg_clustering = nx.average_clustering(G)
clustering_results.append({
"subject_id": sub_id,
"session_id": ses_id,
"average_clustering": avg_clustering
})
print(f"{sub_id} | {ses_id}: {avg_clustering:.4f}")
except Exception as e:
print(f"Error processing {adj_file}: {e}")
#convert to df
df_out = pd.DataFrame(clustering_results)
#pivot the dataframe so that each session becomes a column, subject_id forms the rows,
#session_id serves as the column names, and global_efficiency values fill the table.
clustering_df = df_out.pivot(index="subject_id", columns="session_id", values="average_clustering").reset_index()
clustering_df.columns.name = None # remove column index name
#save to CSV
output_file = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/average_clustering_all_subjects.csv"
os.makedirs(os.path.dirname(output_file), exist_ok=True)
clustering_df.to_csv(output_file, index=False)
print(f"Final results saved to: {output_file}")
# %% [markdown]
# PATH LENGTH
# %%
#sessions
session_labels = {
"ses-1": "baseline",
"ses-2": "acute",
"ses-3": "chronic"
}
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
path_length_results = []
#loop through each session and subject
for ses_id in session_labels.keys():
ses_path = os.path.join(base_path, ses_id)
subject_paths = glob.glob(os.path.join(ses_path, "sub-*"))
for sub_path in subject_paths:
sub_id = os.path.basename(sub_path)
adj_path = os.path.join(sub_path, "Graphs/wAALours", "*cost_400.txt")
adj_files = glob.glob(adj_path)
for adj_file in adj_files:
try:
df = pd.read_csv(adj_file, sep="\s+", header=None)
X = df.to_numpy()
G = nx.from_numpy_array(X)
if nx.is_connected(G):
avg_path_length = nx.average_shortest_path_length(G)
else:
largest_cc = max(nx.connected_components(G), key=len)
G_largest_cc = G.subgraph(largest_cc)
avg_path_length = nx.average_shortest_path_length(G_largest_cc)
path_length_results.append({
"subject_id": sub_id,
"session_id": ses_id,
"average_path_length": avg_path_length
})
print(f"{sub_id} | {ses_id}: {avg_path_length:.4f}")
except Exception as e:
print(f"Error processing {adj_file}: {e}")
#convert to df
df_out = pd.DataFrame(path_length_results)
#pivot the dataframe so that each session becomes a column, subject_id forms the rows,
#session_id serves as the column names, and global_efficiency values fill the table.
path_length_df = df_out.pivot(index="subject_id", columns="session_id", values="average_path_length").reset_index()
path_length_df.columns.name = None # remove column index name
#save to CSV
output_file = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/average_path_length_all_subjects.csv"
os.makedirs(os.path.dirname(output_file), exist_ok=True)
path_length_df.to_csv(output_file, index=False)
print(f"Final results saved to: {output_file}")
# %% [markdown]
# MODULARITY
# %%
#sessions
session_labels = {
"ses-1": "baseline",
"ses-2": "acute",
"ses-3": "chronic"
}
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
modularity_results = []
#loop through each session and subject
for ses_id in session_labels.keys():
ses_path = os.path.join(base_path, ses_id)
subject_paths = glob.glob(os.path.join(ses_path, "sub-*"))
for sub_path in subject_paths:
sub_id = os.path.basename(sub_path)
adj_path = os.path.join(sub_path, "Graphs/wAALours", "*cost_400.txt")
adj_files = glob.glob(adj_path)
for adj_file in adj_files:
try:
df = pd.read_csv(adj_file, sep="\s+", header=None)
X = df.to_numpy()
G = nx.from_numpy_array(X)
communities = community.greedy_modularity_communities(G)
modularity = community.modularity(G, communities)
modularity_results.append({
"subject_id": sub_id,
"session_id": ses_id,
"modularity": modularity
})
print(f"{sub_id} | {ses_id}: {modularity:.4f}")
except Exception as e:
print(f"Error processing {adj_file}: {e}")
#convert to df
df_out = pd.DataFrame(modularity_results)
#pivot the dataframe so that each session becomes a column, subject_id forms the rows,
#session_id serves as the column names, and global_efficiency values fill the table.
modularity_df = df_out.pivot(index="subject_id", columns="session_id", values="modularity").reset_index()
modularity_df.columns.name = None # remove column index name
#save to CSV
output_file = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/modularity_all_subjects.csv"
os.makedirs(os.path.dirname(output_file), exist_ok=True)
modularity_df.to_csv(output_file, index=False)
print(f"Final results saved to: {output_file}")
# %% [markdown]
# MODULARITY, COMMUNITY DISTANCES
# %%
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
coord_path = os.path.join(base_path, "ses-1", "coord_AALours.txt")
coord = pd.read_csv(coord_path, sep="\s+", header=None).to_numpy()
session_labels = {
"ses-1": "baseline",
"ses-2": "acute",
"ses-3": "chronic"
}
community_summary = []
# --- DISTANCE FUNCTIONS ---
def compute_spatial_distances(communities, coord):
comm_centroids = []
for comm in communities:
comm_coords = coord[list(comm)]
centroid = comm_coords.mean(axis=0)
comm_centroids.append(centroid)
return cdist(comm_centroids, comm_centroids, metric='euclidean')
def compute_graph_distances(G, communities):
n = len(communities)
dist_matrix = np.zeros((n, n))
for i in range(n):
for j in range(i + 1, n):
paths = []
for node_i in communities[i]:
for node_j in communities[j]:
try:
length = nx.shortest_path_length(G, source=node_i, target=node_j)
paths.append(length)
except nx.NetworkXNoPath:
continue
avg_dist = np.mean(paths) if paths else np.nan
dist_matrix[i, j] = dist_matrix[j, i] = avg_dist
return dist_matrix
for ses_id in session_labels.keys():
ses_path = os.path.join(base_path, ses_id)
subject_paths = glob.glob(os.path.join(ses_path, "sub-*"))
for sub_path in subject_paths:
sub_id = os.path.basename(sub_path)
adj_path_pattern = os.path.join(sub_path, "Graphs/wAALours", "*cost_400.txt")
adj_files = glob.glob(adj_path_pattern)
for adj_file in adj_files:
try:
df = pd.read_csv(adj_file, sep="\s+", header=None)
adj_matrix = df.to_numpy()
G = nx.from_numpy_array(adj_matrix)
communities = community.greedy_modularity_communities(G)
modularity_value = community.modularity(G, communities)
num_communities = len(communities)
node_community_map = {node: idx for idx, comm in enumerate(communities) for node in comm}
cmap = cm.get_cmap('tab20', num_communities)
node_colors = [cmap(node_community_map[n]) for n in G.nodes()]
#views
views = {
'sagittal': {k: (coord[k, 0], coord[k, 2]) for k in range(len(G.nodes))},
'axial': {k: (coord[k, 0], coord[k, 1]) for k in range(len(G.nodes))},
'coronal': {k: (coord[k, 1], coord[k, 2]) for k in range(len(G.nodes))}
}
#plot
fig, axs = plt.subplots(1, 3, figsize=(18, 6))
for i, (view_name, pos) in enumerate(views.items()):
nx.draw(G, pos, node_color=node_colors, with_labels=False, node_size=100, ax=axs[i])
axs[i].set_title(f'{view_name.capitalize()} view')
adj_basename = os.path.basename(adj_file).replace(".txt", "")
fig.suptitle(f"{sub_id} | {ses_id} | Community Structure\n{adj_basename}", fontsize=16)
plt.tight_layout()
output_dir = os.path.join(sub_path, "Graphs/wAALours", "graph_metrics")
os.makedirs(output_dir, exist_ok=True)
plot_file = os.path.join(output_dir, f"{adj_basename}_communities.png")
plt.savefig(plot_file)
plt.close()
#distances
spatial_distances = compute_spatial_distances(communities, coord)
graph_distances = compute_graph_distances(G, communities)
#save matrices
spatial_file = os.path.join(output_dir, f"{sub_id}_{ses_id}_communities_spatial_distances.csv")
graph_file = os.path.join(output_dir, f"{sub_id}_{ses_id}_communities_graph_distances.csv")
pd.DataFrame(spatial_distances).to_csv(spatial_file, index=False)
pd.DataFrame(graph_distances).to_csv(graph_file, index=False)
#add to global summary
avg_spatial_dist = np.nanmean(spatial_distances[np.triu_indices(num_communities, k=1)])
avg_graph_dist = np.nanmean(graph_distances[np.triu_indices(num_communities, k=1)])
community_summary.append({
"subject_id": sub_id,
"session_id": ses_id,
"modularity": modularity_value,
"num_communities": num_communities,
"avg_spatial_distance": avg_spatial_dist,
"avg_graph_distance": avg_graph_dist
})
print(f"{sub_id} {ses_id}: saved plot and metrics to {output_dir}")
except Exception as e:
print(f"Error processing {adj_file}: {e}")
output_base_dir = os.path.join(base_path, "graph_metrics")
os.makedirs(output_base_dir, exist_ok=True)
#convert summary to df
summary_df = pd.DataFrame(community_summary)
#create and save pivoted dfs
modularity_df = summary_df.pivot(index="subject_id", columns="session_id", values="modularity").reset_index()
num_communities_df = summary_df.pivot(index="subject_id", columns="session_id", values="num_communities").reset_index()
avg_spatial_distance_df = summary_df.pivot(index="subject_id", columns="session_id", values="avg_spatial_distance").reset_index()
avg_graph_distance_df = summary_df.pivot(index="subject_id", columns="session_id", values="avg_graph_distance").reset_index()
for df in [modularity_df, num_communities_df, avg_spatial_distance_df, avg_graph_distance_df]:
df.columns.name = None
#save to CSV
modularity_df.to_csv(os.path.join(output_base_dir, "modularity_all_subjects.csv"), index=False)
num_communities_df.to_csv(os.path.join(output_base_dir, "num_communities_all_subjects.csv"), index=False)
avg_spatial_distance_df.to_csv(os.path.join(output_base_dir, "avg_spatial_distance_all_subjects.csv"), index=False)
avg_graph_distance_df.to_csv(os.path.join(output_base_dir, "avg_graph_distance_all_subjects.csv"), index=False)
print(f"All community metrics saved to: {output_base_dir}")
# %% [markdown]
# ## Statistical comparisons of global metrics
# %%
glob_eff = pd.read_csv("/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/global_efficiency_all_subjects.csv")
clustering_df = pd.read_csv("/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/average_clustering_all_subjects.csv")
path_length_df = pd.read_csv("/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/average_path_length_all_subjects.csv")
modularity_df = pd.read_csv("/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/modularity_all_subjects.csv")
num_communities_df = pd.read_csv("/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/num_communities_all_subjects.csv")
avg_spatial_distance_df = pd.read_csv("/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/avg_spatial_distance_all_subjects.csv")
avg_graph_distance_df = pd.read_csv("/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/avg_graph_distance_all_subjects.csv")
# %%
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics"
metrics = {
"Global efficiency": glob_eff,
"Average clustering": clustering_df,
"Average path length": path_length_df,
"Modularity": modularity_df,
"Average graph distance": avg_graph_distance_df #between communities
}
lmm_results = []
#linear mixed model
def run_lmm_with_ref(df_long, metric_name, ref_session):
#session have to be "categories", otherwise they would be treated as continous variables
df_long["session"] = pd.Categorical(df_long["session"], categories=["ses-1", "ses-2", "ses-3"])
#by default LMM compares everything to ses-1 (alphabetic order), but we want to have the ses-2 vs ses-3 comparison as well, so here we will
# build a new category order with the reference session first. Reference session is defined lower in the loop
df_long["session"] = df_long["session"].cat.reorder_categories([ref_session,
*(s for s in ["ses-1", "ses-2", "ses-3"] if s != ref_session)])
try:
model = smf.mixedlm("value ~ session", data=df_long, groups=df_long["subject_id"])
result = model.fit()
conf_int = result.conf_int()
print(f"\n{metric_name} (reference: {ref_session}):\n", result.summary())
for param in result.fe_params.index:
lmm_results.append({
"Metric": metric_name,
"Reference Session": ref_session,
"Effect": param,
"Estimate": result.fe_params[param],
"CI Lower Bound": conf_int.loc[param][0],
"CI Upper Bound": conf_int.loc[param][1],
"p-value": result.pvalues[param]
})
except Exception as e:
print(f"Error fitting model for {metric_name} with ref {ref_session}: {e}")
for metric_name, df_metric in metrics.items():
df_long = df_metric.melt(id_vars="subject_id",
var_name="session",
value_name="value")
#run lmm with ses-1 as reference (default)
run_lmm_with_ref(df_long.copy(), metric_name, ref_session="ses-1")
#run lmm with ses-2 as reference to get ses-3 vs ses-2 - this "ref_session" is used in our function "run_lmm_with_ref"
run_lmm_with_ref(df_long.copy(), metric_name, ref_session="ses-2")
#save
lmm_df = pd.DataFrame(lmm_results)
output_file = os.path.join(base_path, "lmm_results_all_metrics_with_ref.csv")
lmm_df.to_csv(output_file, index=False)
print(f"\nLinear mixed model results saved to: {output_file}")
# %%
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.stats.multitest import multipletests
# LMM results with CIs
file_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/lmm_results_all_metrics_with_ref.csv"
lmm_df = pd.read_csv(file_path)
selected_metrics = [
"Global efficiency",
"Average clustering",
"Average path length",
"Modularity",
"Average graph distance"
]
custom_titles = [
"TSD vs RW", # ses-2 vs ses-1 (control)
"CSR vs RW", # ses-3 vs ses-1 (control)
"CSR vs TSD" # ses-3 vs ses-2 (control)
]
#remove intercepts and filter metrics
lmm_df = lmm_df[~lmm_df["Effect"].str.contains("Intercept")]
lmm_df = lmm_df[lmm_df["Metric"].isin(selected_metrics)].copy()
#extract comparison info
lmm_df["Compared Session"] = lmm_df["Effect"].str.extract(r'session\[T\.(.*?)\]')
lmm_df["Comparison"] = lmm_df["Compared Session"] + " vs " + lmm_df["Reference Session"]
#add new columns for FDR correction
lmm_df["p-value (FDR)"] = None
lmm_df["Significant (FDR)"] = None
#FDR correction
for (comp_ses, ref_ses) in [
("ses-2", "ses-1"),
("ses-3", "ses-1"),
("ses-3", "ses-2")
]:
mask = (lmm_df["Compared Session"] == comp_ses) & (lmm_df["Reference Session"] == ref_ses)
pvals = lmm_df.loc[mask, "p-value"].values
if len(pvals) > 0:
rejected, pvals_corrected, _, _ = multipletests(pvals, method='fdr_bh', alpha=0.05)
#corrected p-values and rejection flags - decisions
lmm_df.loc[mask, "p-value (FDR)"] = pvals_corrected
lmm_df.loc[mask, "Significant (FDR)"] = rejected
#convert to correct types - numeric for p value and boolean for decision
lmm_df["p-value (FDR)"] = pd.to_numeric(lmm_df["p-value (FDR)"])
lmm_df["Significant (FDR)"] = lmm_df["Significant (FDR)"].astype(bool)
#significance markers and colors
def get_sig_marker_fdr(p):
if p <= 0.001: return '***'
elif p <= 0.01: return '**'
elif p <= 0.05: return '*'
elif p <= 0.1: return '#'
else: return ''
def get_color(p):
if p < 0.05: return 'darkturquoise'
elif p < 0.1: return 'black'
else: return 'gray'
lmm_df["Significance (FDR)"] = lmm_df["p-value (FDR)"].apply(get_sig_marker_fdr)
lmm_df["Color"] = lmm_df["p-value (FDR)"].apply(get_color)
#plot
comparison_order = [
("ses-2", "ses-1"),
("ses-3", "ses-1"),
("ses-3", "ses-2")
]
fig, axes = plt.subplots(1, 3, figsize=(18, 6), sharey=True)
for i, (comp_ses, ref_ses) in enumerate(comparison_order):
df_comp = lmm_df[(lmm_df["Compared Session"] == comp_ses) & (lmm_df["Reference Session"] == ref_ses)].copy()
metric_order = selected_metrics[::-1]
df_comp["Metric"] = pd.Categorical(df_comp["Metric"], categories=metric_order, ordered=True)
df_comp = df_comp.sort_values("Metric")
ax = axes[i]
ax.axvline(x=0, color='gray', linestyle='--')
ax.set_xlim(-0.6, 0.6)
for _, row in df_comp.iterrows():
y = row["Metric"]
x = row["Estimate"]
x_low = row["CI Lower Bound"]
x_high = row["CI Upper Bound"]
color = row["Color"]
sig = row["Significance (FDR)"]
label = f'{x:.4f}{sig}'
ax.plot([x_low, x_high], [y, y], color=color, linewidth=4)
ax.plot(x, y, 'o', color=color, markersize=10)
ax.text(min(x_high + 0.03, 0.48), y, label, va='center', fontsize=13)
ax.set_title(custom_titles[i], fontsize=18)
ax.set_xlabel("Estimate", fontsize=16)
if i == 0:
ax.set_ylabel("Graph metric", fontsize=16)
else:
ax.set_ylabel("")
ax.tick_params(axis='y', labelsize=14)
# plt.suptitle("LMM coefficients and CI (FDR corrected)", fontsize=15)
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
#save to CSV
output_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/lmm_results_with_FDR.csv"
lmm_df.to_csv(output_path, index=False)
# %% [markdown]
# ## Unpaired permutations (broken pairs)
# %%
############### permutations for checking the robustness of LMM for global metrics - broken pairs ##############
#metrics
metrics = {
"Global efficiency": glob_eff,
"Average clustering": clustering_df,
"Average path length": path_length_df,
"Modularity": modularity_df,
"Average graph distance": avg_graph_distance_df
}
n_iterations = 10000
n_subjects = 28
session_names = ['ses-1', 'ses-2', 'ses-3']
session_pairs = [
('ses-1', 'ses-2'),
('ses-1', 'ses-3'),
('ses-2', 'ses-3')
]
comparison_labels = [
"Acute vs Control",
"Chronic vs Control",
"Chronic vs Acute"
]
# Loop over each metric
for metric_name, df_metric in metrics.items():
print(f"\n=== {metric_name} ===")
# stack all session data together in one big list
list_all = np.hstack([df_metric['ses-1'], df_metric['ses-2'], df_metric['ses-3']])
stats = np.zeros((n_iterations, 3)) #place for storing t-stats for 3 comparisons
for i in range(n_iterations):
#now we want to randomly shuffle the values, with replacement, using resample
new_sample = resample(list_all)
#here we define the "new" groups - numbers are defined automatically based on the number of participants in the original groups
new_ses_1 = new_sample[0:n_subjects]
new_ses_2 = new_sample[n_subjects:2*n_subjects]
new_ses_3 = new_sample[2*n_subjects:3*n_subjects]
#pairwise t-tests
stat12, _ = ttest_rel(new_ses_1, new_ses_2)
stat13, _ = ttest_rel(new_ses_1, new_ses_3)
stat23, _ = ttest_rel(new_ses_2, new_ses_3)
stats[i, 0] = stat12
stats[i, 1] = stat13
stats[i, 2] = stat23
#real t-tests (non-permuted, without resampling)
stat12_real, p12_real = ttest_rel(df_metric['ses-1'], df_metric['ses-2'])
stat13_real, p13_real = ttest_rel(df_metric['ses-1'], df_metric['ses-3'])
stat23_real, p23_real = ttest_rel(df_metric['ses-2'], df_metric['ses-3'])
real_stats = [stat12_real, stat13_real, stat23_real]
real_pvals = [p12_real, p13_real, p23_real]
#FDR correction
rejected, real_pvals_fdr, _, _ = multipletests(real_pvals, alpha=0.05, method='fdr_bh')
#plot
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
for j, (ax, (ses_a, ses_b)) in enumerate(zip(axes, session_pairs)):
ax.hist(stats[:, j], bins=30, alpha=0.7)
ax.axvline(real_stats[j], color='red', linestyle='dashed', linewidth=2, label=f'Real t-stat\nt = {real_stats[j]:.2f}\np (FDR)= {real_pvals_fdr[j]:.4f}')
ax.set_title(f'{metric_name}\n{comparison_labels[j]}\n(Unpaired permutation)', fontsize=16)
ax.set_xlabel('t-statistic', fontsize=14)
ax.set_ylabel('Frequency', fontsize=14)
# ax.text(
# 0.85, 0.85,
# f"t = {real_stats[j]:.2f}\n"
# f"p = {real_pvals[j]:.4f}\n"
# f"FDR p = {real_pvals_fdr[j]:.4f}",
# ha='right', va='top',
# transform=ax.transAxes,
# fontsize=10,
# bbox=dict(boxstyle="round,pad=0.3", edgecolor="black", facecolor="white")
# )
ax.legend(fontsize=12)
plt.tight_layout()
plt.show()
#print real and FDR corrected values
print('Real t-statistics, raw p-values, and FDR-corrected p-values:')
for label, t_val, p_val, p_fdr, sig in zip(comparison_labels, real_stats, real_pvals, real_pvals_fdr, rejected):
significance = "SIGNIFICANT" if sig else "NON-SIGNIFICANT"
print(f"{label}: t = {t_val:.3f}, p = {p_val:.4f}, FDR p = {p_fdr:.4f} [{significance}]")
# %% [markdown]
# ## paired permutations (pairs kept, sleep deprivation labels randomly changed)
# %%
############### permutations for checking the robustness of LMM for global metrics - paired permutations ##############
metrics = {
"Global efficiency": glob_eff,
"Average clustering": clustering_df,
"Average path length": path_length_df,
"Modularity": modularity_df,
"Average graph distance": avg_graph_distance_df
}
n_iterations = 10000
session_pairs = [
('ses-1', 'ses-2'),
('ses-1', 'ses-3'),
('ses-2', 'ses-3')
]
comparison_labels = [
"Acute vs Control",
"Chronic vs Control",
"Chronic vs Acute"
]
for metric_name, df_metric in metrics.items():
print(f"\n=== {metric_name} ===")
n_subjects = len(df_metric)
#extract values for each session
ses1 = df_metric['ses-1'].values
ses2 = df_metric['ses-2'].values
ses3 = df_metric['ses-3'].values
#store session values in a dictionary
session_data = {
'ses-1': ses1,
'ses-2': ses2,
'ses-3': ses3
}
#prepare array for storing permutation t-stats and real stats
perm_stats = np.zeros((n_iterations, 3))
real_stats = []
real_pvals = []
#loop over comparisons, pairwise
for idx, (ses_a, ses_b) in enumerate(session_pairs):
x = session_data[ses_a]
y = session_data[ses_b]
#real t-statistic from "session_pairs"
t_real, p_real = ttest_rel(x, y)
real_stats.append(t_real)
real_pvals.append(p_real)
#permutation testing: subject-wise label flipping - random fliping of ses_a and ses_b, but within pair, 10000 times we will calculate
#the differences (ttest_dep) based on mixed groups (randomly flipped labels)
for i in range(n_iterations):
flip = np.random.choice([True, False], size=n_subjects) #this will generate list of "decisions" about flipping, lenght of this list
#will be equal to the number of participants (in my case 28), so each participant will be or will not be swapped
x_perm = np.where(flip, x, y)
y_perm = np.where(flip, y, x)
#tstat after permutations
t_perm, _ = ttest_rel(x_perm, y_perm)
perm_stats[i, idx] = t_perm
# FDR Correction
rejected, real_pvals_fdr, _, _ = multipletests(real_pvals, alpha=0.05, method='fdr_bh')
#plot
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
for j, ax in enumerate(axes):
ax.hist(perm_stats[:, j], bins=30, alpha=0.7)
ax.axvline(real_stats[j], color='red', linestyle='dashed', linewidth=2, label=f'Real t-stat\nt = {real_stats[j]:.2f}\np (FDR)= {real_pvals_fdr[j]:.4f}')
ax.set_title(f'{metric_name}\n{comparison_labels[j]}\n(Paired permutation)', fontsize=16)
ax.set_xlabel('t-statistic', fontsize=14)
ax.set_ylabel('Frequency', fontsize=14)
# ax.text(
# 0.85, 0.85,
# f"t = {real_stats[j]:.2f}\n"
# f"p = {real_pvals[j]:.4f}\n"
# f"FDR p = {real_pvals_fdr[j]:.4f}",
# ha='right', va='top',
# transform=ax.transAxes,
# fontsize=10,
# bbox=dict(boxstyle="round,pad=0.3", edgecolor="black", facecolor="white")
# )
ax.legend(fontsize=12)
plt.tight_layout()
plt.show()
#print real and FDR corrected values
print('Real t-statistics, raw p-values, and FDR-corrected p-values:')
for label, t_val, p_val, p_fdr, sig in zip(comparison_labels, real_stats, real_pvals, real_pvals_fdr, rejected):
significance = "SIGNIFICANT" if sig else "NON-SIGNIFICANT"
print(f"{label}: t = {t_val:.3f}, p = {p_val:.4f}, FDR p = {p_fdr:.4f} [{significance}]")
# %% [markdown]
# # nodal metrics
# %%
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
session_labels = {
"ses-1": "baseline",
"ses-2": "acute",
"ses-3": "chronic"
}
metrics_list = ["Closeness", "Betweenness", "Clustering", "Degree"]
n_regions = 89
output_dir = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/nodal_metrics_LMM"
os.makedirs(output_dir, exist_ok=True)
warnings.filterwarnings("ignore")
lmm_results = []
def run_lmm_with_ref(df_long, metric_name, ref_session, region_index):
df_long["session"] = pd.Categorical(df_long["session"], categories=["ses-1", "ses-2", "ses-3"])
df_long["session"] = df_long["session"].cat.reorder_categories([ref_session,
*(s for s in ["ses-1", "ses-2", "ses-3"] if s != ref_session)])
try:
model = mixedlm("value ~ session", data=df_long, groups=df_long["subject_id"])
result = model.fit()
for param in result.fe_params.index:
if param == "Intercept":
continue
pval = result.pvalues[param]
estimate = result.fe_params[param]
lmm_results.append({
"Region": region_index,
"Metric": metric_name,
"Reference Session": ref_session,
"Effect": param,
"Estimate": estimate,
"p-value": pval
})
except Exception as e:
print(f"Error fitting model for {metric_name} with ref {ref_session}, region {region_index}: {e}")
#loop over metrics
for metric in metrics_list:
print(f"\nProcessing metric: {metric}")
pvals_per_region = np.full((n_regions, 3), np.nan)
estimates_per_region = np.full((n_regions, 3), np.nan)
for region in range(n_regions):
data_long = []
for ses_id in session_labels.keys():
ses_path = os.path.join(base_path, ses_id)
subject_paths = glob.glob(os.path.join(ses_path, "sub-*"))
for sub_path in subject_paths:
sub_id = os.path.basename(sub_path)
graph_path = os.path.join(sub_path, "Graphs/wAALours", "graph_metrics")
metric_files = glob.glob(os.path.join(graph_path, "*_metrics.csv"))
for metrics_file in metric_files:
try:
df_metrics = pd.read_csv(metrics_file)
value = df_metrics.loc[df_metrics['Node'] == region, metric].values[0]
data_long.append({
'subject_id': sub_id,
'session': ses_id,
'value': value
})
except Exception as e:
print(f"Error processing {metrics_file}: {e}")
continue
df_long = pd.DataFrame(data_long)
if df_long['session'].nunique() < 2 or df_long['subject_id'].nunique() < 2:
continue
#reference sessions 1 and 2!!!
run_lmm_with_ref(df_long.copy(), metric, ref_session="ses-1", region_index=region)
run_lmm_with_ref(df_long.copy(), metric, ref_session="ses-2", region_index=region)
df_lmm = pd.DataFrame(lmm_results)
for i in range(n_regions):
# p-values and estimates
for j, (ref, contrast, col) in enumerate([
("ses-1", "session[T.ses-2]", 0),
("ses-1", "session[T.ses-3]", 1),
("ses-2", "session[T.ses-3]", 2)
]):
match = df_lmm[(df_lmm["Metric"] == metric) & (df_lmm["Region"] == i) &
(df_lmm["Reference Session"] == ref) &
(df_lmm["Effect"] == contrast)]
if not match.empty:
pvals_per_region[i, col] = match["p-value"].values[0]
estimates_per_region[i, col] = match["Estimate"].values[0]
#save both results
df_pvals = pd.DataFrame(pvals_per_region, columns=["ses1_vs_ses2", "ses1_vs_ses3", "ses2_vs_ses3"])
df_ests = pd.DataFrame(estimates_per_region, columns=["ses1_vs_ses2", "ses1_vs_ses3", "ses2_vs_ses3"])
df_pvals.to_csv(os.path.join(output_dir, f"{metric.lower()}_pvals_per_region_LMM.csv"), index=False)
df_ests.to_csv(os.path.join(output_dir, f"{metric.lower()}_estimates_per_region_LMM.csv"), index=False)
print(f"Saved: {metric.lower()}_pvals_per_region_LMM.csv")
print(f"Saved: {metric.lower()}_estimates_per_region_LMM.csv")
# %%
#folder where the per-region p-value CSVs are saved
output_dir = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics"
#list of metrics
metrics_list = ["closeness", "betweenness", "clustering", "degree"]
#load all p-values into a dictionary
pvals_all = {}
for metric in metrics_list:
file_path = os.path.join(output_dir, f"{metric}_pvals_per_region_LMM.csv")
if os.path.exists(file_path):
df = pd.read_csv(file_path)
pvals_all[metric.capitalize()] = df.values
else:
print(f"File not found: {file_path}")
# %%
#plot the histograms
session_pairs = ['Ses-1 (control) vs Ses-2 (acute)', 'Ses-1 (control) vs Ses-3 (chronic)', 'Ses-2 (acute) vs Ses-3 (chronic)']
#loop over metrics, pvals_all is the dictionary from the previous code
for metric_name, pvals in pvals_all.items():
print(f"\nPlotting for {metric_name}")
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for i, ax in enumerate(axes):
ax.hist(pvals[:, i], bins=20)
ax.set_title(f'{metric_name} - {session_pairs[i]}')
ax.set_xlabel('p-value')
ax.set_ylabel('Number of regions')
plt.tight_layout()
plt.show()
# %% [markdown]
# ## Multiple comparisons correction (nodal level)
# %%
correction_results = {}
session_comparisons = ['ses1_vs_ses2', 'ses1_vs_ses3', 'ses2_vs_ses3']
for metric_name, pvals in pvals_all.items():
print(f"\n=== Multiple comparisons correction for {metric_name} ===")
correction_results[metric_name] = {}
for i_session, session_comparison in enumerate(session_comparisons):
print(f"\n--- {session_comparison } ---")
pvals_session = pvals[:, i_session] # p-values for this sessions comparison
#bonferroni correction
rej_bonf = bonferroni(pvals_session, alpha=0.05)
print('Bonferroni reject decisions:')
print(rej_bonf)
# FDR correction (Benjamini-Hochberg LSU)
rej_fdr = lsu(pvals_session, q=0.1)
print('FDR (LSU) reject decisions:')
print(rej_fdr)
#save into dictionary
correction_results[metric_name][session_comparison] = {
"bonferroni": rej_bonf,
"fdr": rej_fdr
}
# %%
######## plot significant ROIs of the graph on the brain #####
os.makedirs("/Users/patrycjascislewska/Analizy_neuro/Graphs/significant_nodes", exist_ok=True)
#ROI names
roi_table_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/AAL_ours_89_regions_list.csv"
roi_df = pd.read_csv(roi_table_path, sep=";")
roi_df.columns = roi_df.columns.str.strip()
roi_lookup = dict(zip(roi_df["Node_number"], roi_df["Region"]))
#coordinates and adjacency matrix
coord_path = "/Users/patrycjascislewska/Analizy_neuro/AAL3/AALours_coords_MNI.txt"
coord = pd.read_csv(coord_path, sep="\s+", header=None).to_numpy()
adj_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/ses-1/sub-14/Graphs/wAALours/Adj_mat_n.levels_3_n.regions_89proc_length384.cost_400.txt"
adj_matrix = pd.read_csv(adj_path, sep="\s+", header=None).to_numpy()
G = nx.from_numpy_array(adj_matrix)
#metrics and session comparisons
metrics_list = ["Closeness", "Betweenness", "Clustering", "Degree"]
session_comparisons = ['ses1_vs_ses2', 'ses1_vs_ses3', 'ses2_vs_ses3']
#comparison labels
comparison_labels_map = {
'ses1_vs_ses2': 'Acute vs Control',
'ses1_vs_ses3': 'Chronic vs Control',
'ses2_vs_ses3': 'Chronic vs Acute'
}
#file mapping for LMM coefficient paths
coef_file_map = {
"Degree": "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/nodal_metrics_LMM/degree_estimates_per_region_LMM.csv",
"Clustering": "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/nodal_metrics_LMM/clustering_estimates_per_region_LMM.csv",
"Closeness": "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/nodal_metrics_LMM/closeness_estimates_per_region_LMM.csv",
"Betweenness": "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/nodal_metrics_LMM/betweenness_estimates_per_region_LMM.csv"
}
#main plotting loop
for metric_name in metrics_list:
#load coefficients
coef_path = coef_file_map[metric_name]
coef_df = pd.read_csv(coef_path, sep=",")
coef_df.columns = coef_df.columns.str.strip()
for session_comparison in session_comparisons:
print(f"\n=== {metric_name} - {session_comparison} ===")
fdr_significant = correction_results[metric_name][session_comparison]["fdr"].astype(bool)
coef_session = coef_df[session_comparison].to_numpy()
#encode colors: -1 = ↓, 1 = ↑, 0 = n.s.
node_colors = np.zeros_like(coef_session)
node_colors[fdr_significant & (coef_session > 0)] = 1
node_colors[fdr_significant & (coef_session < 0)] = -1
#significant nodes with region names and coordinates
print("Significant nodes (FDR-corrected):")
significant_nodes_data = []
for idx, (is_sig, coef) in enumerate(zip(fdr_significant, coef_session)):
if is_sig:
direction = "↑" if coef > 0 else "↓"
roi_name = roi_lookup.get(idx, f"ROI_{idx}")
coord_str = f"({coord[idx,0]:.1f}, {coord[idx,1]:.1f}, {coord[idx,2]:.1f})"
print(f"Node {idx:02d} ({roi_name}): {direction} coef = {coef:.4f} | Coord: {coord_str}")
significant_nodes_data.append({
"Node": idx,
"Region": roi_name,
"X": round(coord[idx, 0], 1),
"Y": round(coord[idx, 1], 1),
"Z": round(coord[idx, 2], 1),
"Coefficient": round(coef, 4),
"Direction": "increase" if coef > 0 else "decrease"
})
#save significant nodes to CSV
if significant_nodes_data:
sig_df = pd.DataFrame(significant_nodes_data)
out_path = f"/Users/patrycjascislewska/Analizy_neuro/Graphs/significant_nodes/{metric_name}_{session_comparison}_sig_nodes.csv"
sig_df.to_csv(out_path, index=False)
print(f"\nSaved significant node data to:\n{out_path}")
#plotting setup
color_map = { -1: '#0080FF', 0: '#A0A0A0', 1: '#FF8000' }
node_color_list = [color_map[val] for val in node_colors]
edge_color = '0.8'
fig, axs = plt.subplots(1, 3, figsize=(18, 6))
# Sagittal view (X-Z)
pos_sagittal = {k: (coord[k, 1], coord[k, 2]) for k in range(89)}
nx.draw(G, pos_sagittal, node_color=node_color_list, edge_color=edge_color,
with_labels=False, node_size=100, ax=axs[0])
axs[0].set_title('Sagittal view (Y-Z)')
# Axial view (X-Y)
pos_axial = {k: (coord[k, 1], coord[k, 0]) for k in range(89)}
nx.draw(G, pos_axial, node_color=node_color_list, edge_color=edge_color,
with_labels=False, node_size=100, ax=axs[1])
axs[1].set_title('Axial view (X-Y)')
# Coronal view (Y-Z)
pos_coronal = {k: (coord[k, 0], coord[k, 2]) for k in range(89)}
nx.draw(G, pos_coronal, node_color=node_color_list, edge_color=edge_color,
with_labels=False, node_size=100, ax=axs[2])
axs[2].set_title('Coronal view (X-Z)')
# Legend
legend_elements = [
Line2D([0], [0], marker='o', color='w', label='Significant increase',
markerfacecolor='#FF8000', markersize=10),
Line2D([0], [0], marker='o', color='w', label='Significant decrease',
markerfacecolor='#0080FF', markersize=10),
Line2D([0], [0], marker='o', color='w', label='Not significant',
markerfacecolor='#A0A0A0', markersize=10)
]
axs[2].legend(handles=legend_elements, loc='upper right', fontsize=10)
fig.suptitle(
f"{metric_name} - {comparison_labels_map.get(session_comparison, session_comparison)}",
fontsize=16
)
plt.tight_layout()
plt.show()
# %% [markdown]
# ## Nodal results on the glass brain template
# %%
##### plotting significant nodes (only) on the glass brain using netplotbrain package
from netplotbrain import plot
sig_nodes_dir = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/significant_nodes"
output_dir = os.path.join(sig_nodes_dir, "brain_plots")
os.makedirs(output_dir, exist_ok=True)
# all significant node CSVs
sig_files = glob.glob(os.path.join(sig_nodes_dir, "*_sig_nodes.csv"))
# color mapping for directions
direction_color_map = {
"increase": "#FF8000", # orange
"decrease": "#0080FF" # blue
}
#loop through all "siginificant_nodes" files and and plot
for file in sig_files:
df = pd.read_csv(file)
if df.empty:
continue
#extract metric and comparison from filename
basename = os.path.basename(file)
parts = basename.replace("_sig_nodes.csv", "").split("_")
metric = parts[0]
comparison = "_".join(parts[1:])
#prepare plotting df to fit the netplotbrain requirements
df_plot = df.copy()
df_plot["x"] = df_plot["X"]
df_plot["y"] = df_plot["Y"]
df_plot["z"] = df_plot["Z"]
df_plot["label"] = df_plot["Region"]
df_plot["color"] = df_plot["Direction"].map(direction_color_map)
fig, ax = plot(
nodes=df_plot,
node_color="color",
node_size=20,
node_labels=False,
view="APLRISs", # here we can change views
template='MNI152NLin6Asym',
template_style='glass', # there is also e.g. surface
title=f"{metric} - {comparison}"
)
#save plot
output_path = os.path.join(output_dir, f"{metric}_{comparison}_brain.png")
fig.savefig(output_path, dpi=300)
plt.close(fig)
# %% [markdown]
# ## Within-Subject Hub Disruption Index
# %% [markdown]
# degree centrality
# %%
#here I use the degree from previously saved files with nodal metrics
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
ses_paths = {
"control": os.path.join(base_path, "ses-1"),
"acute": os.path.join(base_path, "ses-2"),
"chronic": os.path.join(base_path, "ses-3")
}
metrics_pattern = "Graphs/wAALours/graph_metrics/*_metrics.csv"
#subject IDs in all 3 sessions
ids_ses1 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["control"], "sub-*"))}
ids_ses2 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["acute"], "sub-*"))}
ids_ses3 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["chronic"], "sub-*"))}
common_ids = sorted(list(ids_ses1 & ids_ses2 & ids_ses3))
print(f"Subjects present in all 3 sessions: {len(common_ids)}")
#load degree vectors from metrics CSV (previously saved)
def load_degree_vector(subject_id, session):
folder = os.path.join(ses_paths[session], subject_id)
files = glob.glob(os.path.join(folder, metrics_pattern))
if not files:
raise FileNotFoundError(f"No metrics file found for {subject_id} in {session}")
df = pd.read_csv(files[0])
return df["Degree_centrality"].values
#plot multiple histograms - one per subject per comparison
fig, axes = plt.subplots(4, len(common_ids), figsize=(4 * len(common_ids), 12), sharey=True)
#kappas storage
kappa_acute_vs_control = []
kappa_chronic_vs_control = []
kappa_chronic_vs_acute = []
kappa_acute_vs_chronic = []
#main loop
for i, subject_id in enumerate(common_ids):
deg_control = load_degree_vector(subject_id, "control")
deg_acute = load_degree_vector(subject_id, "acute")
deg_chronic = load_degree_vector(subject_id, "chronic")
model = LinearRegression()
# Acute vs control
x = deg_control.reshape((-1, 1))
y = deg_acute - deg_control
model.fit(x, y)
kappa_acute_vs_control.append(model.coef_[0])
axes[0, i].scatter(deg_control, y, alpha=0.7)
axes[0, i].plot(deg_control, model.predict(x), color='red')
axes[0, i].set_title(f"{subject_id}\nκ1={model.coef_[0]:.3f}")
if i == 0:
axes[0, i].set_ylabel("Δ Degree (Acute - Control)")
# Chronic vs control
y = deg_chronic - deg_control
model.fit(x, y)
kappa_chronic_vs_control.append(model.coef_[0])
axes[1, i].scatter(deg_control, y, alpha=0.7)
axes[1, i].plot(deg_control, model.predict(x), color='red')
axes[1, i].set_title(f"κ2={model.coef_[0]:.3f}")
if i == 0:
axes[1, i].set_ylabel("Δ Degree (Chronic - Control)")
# Chronic vs Acute
x = deg_acute.reshape((-1, 1))
y = deg_chronic - deg_acute
model.fit(x, y)
kappa_chronic_vs_acute.append(model.coef_[0])
axes[2, i].scatter(deg_acute, y, alpha=0.7)
axes[2, i].plot(deg_acute, model.predict(x), color='red')
axes[2, i].set_title(f"κ3={model.coef_[0]:.3f}")
if i == 0:
axes[2, i].set_ylabel("Δ Degree (Chronic - Acute)")
# Acute vs Chronic
x = deg_chronic.reshape((-1, 1))
y = deg_acute - deg_chronic
model.fit(x, y)
kappa_acute_vs_chronic.append(model.coef_[0])
axes[3, i].scatter(deg_chronic, y, alpha=0.7)
axes[3, i].plot(deg_chronic, model.predict(x), color='red')
axes[3, i].set_title(f"κ4={model.coef_[0]:.3f}")
if i == 0:
axes[3, i].set_ylabel("Δ Degree (Acute - Chronic)")
#plot
plt.suptitle("Within subject 'HDI' scatterplots using degree", fontsize=16)
plt.tight_layout(rect=[0, 0, 1, 0.95])
plt.show()
# %%
# output_dir = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/HDI_kappa_results_within_subject"
# os.makedirs(output_dir, exist_ok=True)
# output_path = os.path.join(output_dir, "kappas_degree.csv")
# df_kappa.to_csv(output_path, index=False)
# %% [markdown]
# closenness centrality
# %%
#here I use the closeness values from previously saved files with nodal metrics
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
ses_paths = {
"control": os.path.join(base_path, "ses-1"),
"acute": os.path.join(base_path, "ses-2"),
"chronic": os.path.join(base_path, "ses-3")
}
metrics_pattern = "Graphs/wAALours/graph_metrics/*_metrics.csv"
#subject IDs in all 3 sessions
ids_ses1 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["control"], "sub-*"))}
ids_ses2 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["acute"], "sub-*"))}
ids_ses3 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["chronic"], "sub-*"))}
common_ids = sorted(list(ids_ses1 & ids_ses2 & ids_ses3))
print(f"Subjects present in all 3 sessions: {len(common_ids)}")
#load closeness vectors from metrics CSV (previously saved)
def load_Closeness_vector(subject_id, session):
folder = os.path.join(ses_paths[session], subject_id)
files = glob.glob(os.path.join(folder, metrics_pattern))
if not files:
raise FileNotFoundError(f"No metrics file found for {subject_id} in {session}")
df = pd.read_csv(files[0])
return df["Closeness"].values
#plot multiple histograms - one per subject per comparison
fig, axes = plt.subplots(4, len(common_ids), figsize=(4 * len(common_ids), 12), sharey=True)
#kappas storage
kappa_acute_vs_control = []
kappa_chronic_vs_control = []
kappa_chronic_vs_acute = []
kappa_acute_vs_chronic = []
#main loop
for i, subject_id in enumerate(common_ids):
Closeness_control = load_Closeness_vector(subject_id, "control")
Closeness_acute = load_Closeness_vector(subject_id, "acute")
Closeness_chronic = load_Closeness_vector(subject_id, "chronic")
model = LinearRegression()
# Acute vs control
x = Closeness_control.reshape((-1, 1))
y = Closeness_acute - Closeness_control
kappa_acute_vs_control.append(model.coef_[0])
axes[0, i].scatter(Closeness_control, y, alpha=0.7)
axes[0, i].plot(Closeness_control, model.predict(x), color='red')
axes[0, i].set_title(f"{subject_id}\nκ1={model.coef_[0]:.3f}")
if i == 0:
axes[0, i].set_ylabel("Δ Closeness (Acute - Control)")
# Chronic vs control
y = Closeness_chronic - Closeness_control
model.fit(x, y)
kappa_chronic_vs_control.append(model.coef_[0])
axes[1, i].scatter(Closeness_control, y, alpha=0.7)
axes[1, i].plot(Closeness_control, model.predict(x), color='red')
axes[1, i].set_title(f"κ2={model.coef_[0]:.3f}")
if i == 0:
axes[1, i].set_ylabel("Δ Closeness (Chronic - Control)")
# Chronic vs Acute
x = Closeness_acute.reshape((-1, 1))
y = Closeness_chronic - Closeness_acute
model.fit(x, y)
kappa_chronic_vs_acute.append(model.coef_[0])
axes[2, i].scatter(Closeness_acute, y, alpha=0.7)
axes[2, i].plot(Closeness_acute, model.predict(x), color='red')
axes[2, i].set_title(f"κ3={model.coef_[0]:.3f}")
if i == 0:
axes[2, i].set_ylabel("Δ Closeness (Chronic - Acute)")
# Acute vs Chronic
x = Closeness_chronic.reshape((-1, 1))
y = Closeness_acute - Closeness_chronic
model.fit(x, y)
kappa_acute_vs_chronic.append(model.coef_[0])
axes[3, i].scatter(Closeness_chronic, y, alpha=0.7)
axes[3, i].plot(Closeness_chronic, model.predict(x), color='red')
axes[3, i].set_title(f"κ4={model.coef_[0]:.3f}")
if i == 0:
axes[3, i].set_ylabel("Δ Closeness (Acute - Chronic)")
#plot
plt.suptitle("Within subject 'HDI' scatterplots using Closeness", fontsize=16)
plt.tight_layout(rect=[0, 0, 1, 0.95])
plt.show()
# %%
# output_dir = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/HDI_kappa_results_within_subject"
# os.makedirs(output_dir, exist_ok=True)
# output_path = os.path.join(output_dir, "kappas_closeness.csv")
# df_kappa.to_csv(output_path, index=False)
# %% [markdown]
# Clustering coefficient
# %%
#here I use the clustering coefficients from previously saved files with nodal metrics
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
ses_paths = {
"control": os.path.join(base_path, "ses-1"),
"acute": os.path.join(base_path, "ses-2"),
"chronic": os.path.join(base_path, "ses-3")
}
metrics_pattern = "Graphs/wAALours/graph_metrics/*_metrics.csv"
#subject IDs in all 3 sessions
ids_ses1 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["control"], "sub-*"))}
ids_ses2 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["acute"], "sub-*"))}
ids_ses3 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["chronic"], "sub-*"))}
common_ids = sorted(list(ids_ses1 & ids_ses2 & ids_ses3))
print(f"Subjects present in all 3 sessions: {len(common_ids)}")
#load clustering vectors from metrics CSV (previously saved)
def load_Clustering_vector(subject_id, session):
folder = os.path.join(ses_paths[session], subject_id)
files = glob.glob(os.path.join(folder, metrics_pattern))
if not files:
raise FileNotFoundError(f"No metrics file found for {subject_id} in {session}")
df = pd.read_csv(files[0])
return df["Clustering"].values
#plot multiple histograms - one per subject per comparison
fig, axes = plt.subplots(4, len(common_ids), figsize=(4 * len(common_ids), 12), sharey=True)
#kappas storage
kappa_acute_vs_control = []
kappa_chronic_vs_control = []
kappa_chronic_vs_acute = []
kappa_acute_vs_chronic = []
#main loop
for i, subject_id in enumerate(common_ids):
Clustering_control = load_Clustering_vector(subject_id, "control")
Clustering_acute = load_Clustering_vector(subject_id, "acute")
Clustering_chronic = load_Clustering_vector(subject_id, "chronic")
model = LinearRegression()
# Acute vs control
x = Clustering_control.reshape((-1, 1))
y = Clustering_acute - Clustering_control
model.fit(x, y)
kappa_acute_vs_control.append(model.coef_[0])
axes[0, i].scatter(Clustering_control, y, alpha=0.7)
axes[0, i].plot(Clustering_control, model.predict(x), color='red')
axes[0, i].set_title(f"{subject_id}\nκ1={model.coef_[0]:.3f}")
if i == 0:
axes[0, i].set_ylabel("Δ Clustering (Acute - Control)")
# Chronic vs control
y = Clustering_chronic - Clustering_control
model.fit(x, y)
kappa_chronic_vs_control.append(model.coef_[0])
axes[1, i].scatter(Clustering_control, y, alpha=0.7)
axes[1, i].plot(Clustering_control, model.predict(x), color='red')
axes[1, i].set_title(f"κ2={model.coef_[0]:.3f}")
if i == 0:
axes[1, i].set_ylabel("Δ Clustering (Chronic - Control)")
# Chronic vs Acute
x = Clustering_acute.reshape((-1, 1))
y = Clustering_chronic - Clustering_acute
model.fit(x, y)
kappa_chronic_vs_acute.append(model.coef_[0])
axes[2, i].scatter(Clustering_acute, y, alpha=0.7)
axes[2, i].plot(Clustering_acute, model.predict(x), color='red')
axes[2, i].set_title(f"κ3={model.coef_[0]:.3f}")
if i == 0:
axes[2, i].set_ylabel("Δ Clustering (Chronic - Acute)")
# Acute vs Chronic
x = Clustering_chronic.reshape((-1, 1))
y = Clustering_acute - Clustering_chronic
model.fit(x, y)
kappa_acute_vs_chronic.append(model.coef_[0])
axes[3, i].scatter(Clustering_chronic, y, alpha=0.7)
axes[3, i].plot(Clustering_chronic, model.predict(x), color='red')
axes[3, i].set_title(f"κ4={model.coef_[0]:.3f}")
if i == 0:
axes[3, i].set_ylabel("Δ Clustering (Acute - Chronic)")
#plot
plt.suptitle("Within subject 'HDI' scatterplots using Clustering", fontsize=16)
plt.tight_layout(rect=[0, 0, 1, 0.95])
plt.show()
# %%
# output_dir = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/HDI_kappa_results_within_subject"
# os.makedirs(output_dir, exist_ok=True)
# output_path = os.path.join(output_dir, "kappas_clustering.csv")
# df_kappa.to_csv(output_path, index=False)
# %% [markdown]
# ## Validation of within-subject Hub Disruption Index (permutations)
# %%
## kappa within-subject - kappa group ###
#real degree vectors from our csv files
def load_degree_vector(subject_id, session):
folder = os.path.join(ses_paths[session], subject_id)
files = glob.glob(os.path.join(folder, metrics_pattern))
if not files:
raise FileNotFoundError(f"No metrics file found for {subject_id} in {session}")
df = pd.read_csv(files[0])
return df["Degree_centrality"].values
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
ses_paths = {
"control": os.path.join(base_path, "ses-1"),
"acute": os.path.join(base_path, "ses-2")
}
metrics_pattern = "Graphs/wAALours/graph_metrics/*_metrics.csv"
#subject IDs present in both sessions
ids_ses1 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["control"], "sub-*"))}
ids_ses2 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["acute"], "sub-*"))}
common_ids = sorted(list(ids_ses1 & ids_ses2))
print(f"Subjects in both control and acute: {len(common_ids)}")
#load degree vectors
deg_control_all = []
deg_acute_all = []
for subject_id in common_ids:
deg_control = load_degree_vector(subject_id, "control")
deg_acute = load_degree_vector(subject_id, "acute")
deg_control_all.append(deg_control)
deg_acute_all.append(deg_acute)
deg_control_all = np.array(deg_control_all)
deg_acute_all = np.array(deg_acute_all)
#compute group-averaged control degree vector
group_control_mean = np.mean(deg_control_all, axis=0)
#compute real kappa differences (within kappa - group average kappa)
real_diffs = []
model = LinearRegression()
for i in range(len(common_ids)):
subj_control = deg_control_all[i]
subj_acute = deg_acute_all[i]
# 1. Within-subject
delta_within = subj_acute - subj_control
model.fit(subj_control.reshape(-1, 1), delta_within)
kappa_within = model.coef_[0]
# 2. Group-based
delta_groupmean = subj_acute - group_control_mean
model.fit(group_control_mean.reshape(-1, 1), delta_groupmean)
kappa_groupmean = model.coef_[0]
real_diffs.append(kappa_within - kappa_groupmean)
real_diffs = np.array(real_diffs)
real_mean_diff = np.mean(real_diffs)
print(f"Real mean (κ_within - κ_groupmean): {real_mean_diff:.4f}")
# PERMUTATIONS - Shuffle node order in chronic session
n_iterations = 10000
null_diffs = np.zeros(n_iterations)
for it in range(n_iterations):
diffs = []
for i in range(len(common_ids)):
subj_control = deg_control_all[i]
subj_acute = np.random.permutation(deg_acute_all[i]) # shuffle node order
# 1. Within-subject
delta_within = subj_acute - subj_control
model.fit(subj_control.reshape(-1, 1), delta_within)
kappa_within = model.coef_[0]
# 2. Group-based
delta_groupmean = subj_acute - group_control_mean
model.fit(group_control_mean.reshape(-1, 1), delta_groupmean)
kappa_groupmean = model.coef_[0]
diffs.append(kappa_within - kappa_groupmean)
null_diffs[it] = np.mean(diffs)
#plot
plt.figure(figsize=(10, 5))
plt.hist(null_diffs, bins=30, alpha=0.7)
plt.axvline(real_mean_diff, color='red', linestyle='dashed', linewidth=2,
label=f"Real mean diff = {real_mean_diff:.3f}")
plt.xlabel("Mean κ_within − κ_groupmean (permutation distribution)")
plt.ylabel("Frequency")
plt.title("Permutations distribution of 'kappa within - kappa group' after node shuffling \n(acute nodes shuffled, control-acute pairs kept in order)")
plt.legend()
plt.tight_layout()
plt.show()
# P-value
p_value = np.mean(null_diffs <= real_mean_diff)
print(f"P-value: {p_value:.4f}")
# %%
## the "mixed WITHIN PAIRS" aproach ##
#real degree vectors from our csv files
def load_degree_vector(subject_id, session):
folder = os.path.join(ses_paths[session], subject_id)
files = glob.glob(os.path.join(folder, metrics_pattern))
if not files:
raise FileNotFoundError(f"No metrics file found for {subject_id} in {session}")
df = pd.read_csv(files[0])
return df["Degree_centrality"].values
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
ses_paths = {
"control": os.path.join(base_path, "ses-1"),
"chronic": os.path.join(base_path, "ses-3")
}
metrics_pattern = "Graphs/wAALours/graph_metrics/*_metrics.csv"
#subject IDs present in both sessions
ids_ses1 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["control"], "sub-*"))}
ids_ses3 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["chronic"], "sub-*"))}
common_ids = sorted(list(ids_ses1 & ids_ses3))
print(f"Subjects in both control and chronic: {len(common_ids)}")
#load degree vectors
deg_control_all = []
deg_chronic_all = []
for subject_id in common_ids:
deg_control = load_degree_vector(subject_id, "control")
deg_chronic = load_degree_vector(subject_id, "chronic")
deg_control_all.append(deg_control)
deg_chronic_all.append(deg_chronic)
deg_control_all = np.array(deg_control_all)
deg_chronic_all = np.array(deg_chronic_all)
#COMPUTE REAL KAPPAS
real_kappas = []
model = LinearRegression()
for i in range(len(common_ids)):
x = deg_control_all[i].reshape(-1, 1)
y = deg_chronic_all[i] - deg_control_all[i]
model.fit(x, y)
real_kappas.append(model.coef_[0])
real_kappas = np.array(real_kappas)
real_kappa_mean = np.mean(real_kappas)
print(f"Real mean kappa (chronic vs control): {real_kappa_mean:.4f}")
n_iterations = 10000
null_kappa_means = np.zeros(n_iterations)
for it in range(n_iterations):
perm_kappas = []
for i in range(len(common_ids)):
x = deg_control_all[i].reshape(-1, 1)
delta = deg_chronic_all[i] - deg_control_all[i]
if np.random.rand() < 0.5:
delta = -delta
model.fit(x, delta)
perm_kappas.append(model.coef_[0])
null_kappa_means[it] = np.mean(perm_kappas)
#plot
plt.figure(figsize=(10, 5))
plt.hist(null_kappa_means, bins=30, alpha=0.7)
plt.axvline(real_kappa_mean, color='red', linestyle='dashed', linewidth=2,
label=f"Real κ = {real_kappa_mean:.3f}")
plt.xlabel("Mean kappa (permutation distribution)")
plt.ylabel("Frequency")
plt.title("Permutations - distribution of mean kappa (chronic vs control) \n (flip signs of Δy ONLY and keep x axis the same)\n mixed WITHIN pairs")
plt.legend()
plt.tight_layout()
plt.show()
# %%
## the "mixed WITHIN PAIRS" aproach ##
#real vectors from our csv files
def load_metric_vector(subject_id, session, metric_name):
folder = os.path.join(ses_paths[session], subject_id)
files = glob.glob(os.path.join(folder, metrics_pattern))
if not files:
raise FileNotFoundError(f"No metrics file found for {subject_id} in {session}")
df = pd.read_csv(files[0])
return df[metric_name].values
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
ses_paths = {
"control": os.path.join(base_path, "ses-1"),
"acute": os.path.join(base_path, "ses-2")
}
metrics_pattern = "Graphs/wAALours/graph_metrics/*_metrics.csv"
#subject IDs present in both sessions
ids_ses1 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["control"], "sub-*"))}
ids_ses2 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["acute"], "sub-*"))}
common_ids = sorted(list(ids_ses1 & ids_ses2))
print(f"Subjects in both control and acute: {len(common_ids)}")
metric_names = ["Degree_centrality", "Closeness", "Clustering"]
real_kappa_means = {}
null_kappa_distributions = {}
pvals_uncorrected = {}
cohen_ds = {}
# loop over each metric
for metric_name in metric_names:
#load metric vectors
control_all = []
acute_all = []
for subject_id in common_ids:
control = load_metric_vector(subject_id, "control", metric_name)
acute = load_metric_vector(subject_id, "acute", metric_name)
control_all.append(control)
acute_all.append(acute)
control_all = np.array(control_all)
acute_all = np.array(acute_all)
#COMPUTE REAL KAPPAS
real_kappas = []
model = LinearRegression()
for i in range(len(common_ids)):
x = control_all[i].reshape(-1, 1)
y = acute_all[i] - control_all[i]
model.fit(x, y)
real_kappas.append(model.coef_[0])
real_kappas = np.array(real_kappas)
real_kappa_mean = np.mean(real_kappas)
real_kappa_means[metric_name] = real_kappa_mean
print(f"{metric_name} - Real mean kappa (acute vs control): {real_kappa_mean:.4f}")
n_iterations = 10000
null_kappa_means = np.zeros(n_iterations)
for it in range(n_iterations):
perm_kappas = []
for i in range(len(common_ids)):
x = control_all[i].reshape(-1, 1)
delta = acute_all[i] - control_all[i]
if np.random.rand() < 0.5:
delta = -delta
model.fit(x, delta)
perm_kappas.append(model.coef_[0])
null_kappa_means[it] = np.mean(perm_kappas)
null_kappa_distributions[metric_name] = null_kappa_means
#p-value (two-tailed), this counts how many null values are as extreme or more extreme than the absolute value of my observed effect.
p_value_two_tailed = (np.sum(np.abs(null_kappa_means) >= np.abs(real_kappa_mean)) + 1) / (n_iterations + 1)
pvals_uncorrected[metric_name] = p_value_two_tailed
#Cohen's d
mean_null = np.mean(null_kappa_means)
std_null = np.std(null_kappa_means, ddof=1)
cohen_d = (real_kappa_mean - mean_null) / std_null
cohen_ds[metric_name] = cohen_d
#FDR CORRECTION
pvals_list = [pvals_uncorrected[m] for m in metric_names]
rejects, pvals_corrected = fdrcorrection(pvals_list, alpha=0.05)
print("\n==== FDR-corrected results (acute vs control, α = 0.05) ====")
for i, metric_name in enumerate(metric_names):
print(f"{metric_name}:\n"
f" Real κ = {real_kappa_means[metric_name]:.4f}\n"
f" Uncorrected p = {pvals_uncorrected[metric_name]:.4f}\n"
f" FDR-corrected p = {pvals_corrected[i]:.4f}\n"
f" Significant after FDR: {rejects[i]}\n"
f" Cohen's d = {cohen_ds[metric_name]:.2f}\n")
#plot each metric
for i, metric_name in enumerate(metric_names):
plt.figure(figsize=(10, 5))
plt.hist(null_kappa_distributions[metric_name], bins=30, alpha=0.7, edgecolor="black")
plt.axvline(real_kappa_means[metric_name], color='red', linestyle='dashed', linewidth=4)
plt.text(real_kappa_means[metric_name], plt.ylim()[1]*0.7,
f"Real mean κ = {real_kappa_means[metric_name]:.3f}\n"
f"Cohen's d = {cohen_ds[metric_name]:.2f}",
color='black', fontsize=10, ha='left', va='top')
plt.xlabel("κ", fontsize=20)
plt.ylabel("Frequency", fontsize=20)
plt.title(f"Permutation test: {metric_name.replace('_', ' ')} (acute vs control)")
plt.tick_params(axis='x', labelsize = 20)
plt.tick_params(axis='y', labelsize = 20)
# plt.legend()
plt.tight_layout()
plt.show()
# %%
## the "mixed WITHIN PAIRS" aproach ##
#real vectors from our csv files
def load_metric_vector(subject_id, session, metric_name):
folder = os.path.join(ses_paths[session], subject_id)
files = glob.glob(os.path.join(folder, metrics_pattern))
if not files:
raise FileNotFoundError(f"No metrics file found for {subject_id} in {session}")
df = pd.read_csv(files[0])
return df[metric_name].values
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
ses_paths = {
"control": os.path.join(base_path, "ses-1"),
"chronic": os.path.join(base_path, "ses-3")
}
metrics_pattern = "Graphs/wAALours/graph_metrics/*_metrics.csv"
#subject IDs present in both sessions
ids_ses1 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["control"], "sub-*"))}
ids_ses3 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["chronic"], "sub-*"))}
common_ids = sorted(list(ids_ses1 & ids_ses3))
print(f"Subjects in both control and chronic: {len(common_ids)}")
metric_names = ["Degree_centrality", "Closeness", "Clustering"]
real_kappa_means = {}
null_kappa_distributions = {}
pvals_uncorrected = {}
cohen_ds = {}
# loop over each metric
for metric_name in metric_names:
#load metric vectors
control_all = []
chronic_all = []
for subject_id in common_ids:
control = load_metric_vector(subject_id, "control", metric_name)
chronic = load_metric_vector(subject_id, "chronic", metric_name)
control_all.append(control)
chronic_all.append(chronic)
control_all = np.array(control_all)
chronic_all = np.array(chronic_all)
#COMPUTE REAL KAPPAS
real_kappas = []
model = LinearRegression()
for i in range(len(common_ids)):
x = control_all[i].reshape(-1, 1)
y = chronic_all[i] - control_all[i]
model.fit(x, y)
real_kappas.append(model.coef_[0])
real_kappas = np.array(real_kappas)
real_kappa_mean = np.mean(real_kappas)
real_kappa_means[metric_name] = real_kappa_mean
print(f"{metric_name} - Real mean kappa (chronic vs control): {real_kappa_mean:.4f}")
#PERMUTATIONS - COMPUTE "NEW" KAPPAS
n_iterations = 10000
null_kappa_means = np.zeros(n_iterations)
for it in range(n_iterations):
perm_kappas = []
for i in range(len(common_ids)):
x = control_all[i].reshape(-1, 1)
delta = chronic_all[i] - control_all[i]
if np.random.rand() < 0.5:
delta = -delta
model.fit(x, delta)
perm_kappas.append(model.coef_[0])
null_kappa_means[it] = np.mean(perm_kappas)
null_kappa_distributions[metric_name] = null_kappa_means
p_value_two_tailed = (np.sum(np.abs(null_kappa_means) >= np.abs(real_kappa_mean)) + 1) / (n_iterations + 1)
pvals_uncorrected[metric_name] = p_value_two_tailed
#Cohen's d
mean_null = np.mean(null_kappa_means)
std_null = np.std(null_kappa_means, ddof=1)
cohen_d = (real_kappa_mean - mean_null) / std_null
cohen_ds[metric_name] = cohen_d
#FDR CORRECTION
pvals_list = [pvals_uncorrected[m] for m in metric_names]
rejects, pvals_corrected = fdrcorrection(pvals_list, alpha=0.05)
print("\n==== FDR-corrected results (chronic vs control, α = 0.05) ====")
for i, metric_name in enumerate(metric_names):
print(f"{metric_name}:\n"
f" Real κ = {real_kappa_means[metric_name]:.4f}\n"
f" Uncorrected p = {pvals_uncorrected[metric_name]:.4f}\n"
f" FDR-corrected p = {pvals_corrected[i]:.4f}\n"
f" Significant after FDR: {rejects[i]}\n"
f" Cohen's d = {cohen_ds[metric_name]:.2f}\n")
#plot each metric
for i, metric_name in enumerate(metric_names):
plt.figure(figsize=(10, 5))
plt.hist(null_kappa_distributions[metric_name], bins=30, alpha=0.7, edgecolor="black")
plt.axvline(real_kappa_means[metric_name], color='red', linestyle='dashed', linewidth=4)
plt.text(real_kappa_means[metric_name], plt.ylim()[1]*0.7,
f"Real mean κ = {real_kappa_means[metric_name]:.3f}\n"
f"Cohen's d = {cohen_ds[metric_name]:.2f}",
color='black', fontsize=10, ha='left', va='top')
plt.xlabel("κ", fontsize=20)
plt.ylabel("Frequency", fontsize=20)
plt.title(f"Permutation test: {metric_name.replace('_', ' ')} (chronic vs control)")
plt.tick_params(axis='x', labelsize = 20)
plt.tick_params(axis='y', labelsize = 20)
# plt.legend()
plt.tight_layout()
plt.show()
# %% [markdown]
# ## permutations with histogram of real values
# %% [markdown]
# Degree centrality
# %%
###the "mixed WITHIN PAIRS" aproach with the histogram of real values ##
#real degree vectors from our csv files
def load_degree_vector(subject_id, session):
folder = os.path.join(ses_paths[session], subject_id)
files = glob.glob(os.path.join(folder, metrics_pattern))
if not files:
raise FileNotFoundError(f"No metrics file found for {subject_id} in {session}")
df = pd.read_csv(files[0])
return df["Degree_centrality"].values
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
ses_paths = {
"control": os.path.join(base_path, "ses-1"),
"chronic": os.path.join(base_path, "ses-3")
}
metrics_pattern = "Graphs/wAALours/graph_metrics/*_metrics.csv"
#subject IDs present in both sessions
ids_ses1 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["control"], "sub-*"))}
ids_ses3 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["chronic"], "sub-*"))}
common_ids = sorted(list(ids_ses1 & ids_ses3))
print(f"Subjects in both control and chronic: {len(common_ids)}")
#load degree vectors
deg_control_all = []
deg_chronic_all = []
for subject_id in common_ids:
deg_control = load_degree_vector(subject_id, "control")
deg_chronic = load_degree_vector(subject_id, "chronic")
deg_control_all.append(deg_control)
deg_chronic_all.append(deg_chronic)
deg_control_all = np.array(deg_control_all)
deg_chronic_all = np.array(deg_chronic_all)
#COMPUTE REAL KAPPAS
real_kappas = []
model = LinearRegression()
for i in range(len(common_ids)):
x = deg_control_all[i].reshape(-1, 1)
y = deg_chronic_all[i] - deg_control_all[i]
model.fit(x, y)
real_kappas.append(model.coef_[0])
real_kappas = np.array(real_kappas)
real_kappa_mean = np.mean(real_kappas)
print(f"Real mean kappa (chronic vs control): {real_kappa_mean:.4f}")
#PERMUTATIONS - COMPUTE "NEW" KAPPAS
n_iterations = 10000
null_kappa_means = np.zeros(n_iterations)
for it in range(n_iterations):
perm_kappas = []
for i in range(len(common_ids)):
x = deg_control_all[i].reshape(-1, 1)
delta = deg_chronic_all[i] - deg_control_all[i]
if np.random.rand() < 0.5:
delta = -delta
model.fit(x, delta)
perm_kappas.append(model.coef_[0])
null_kappa_means[it] = np.mean(perm_kappas)
#load real participant kappa values
real_kappa_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/HDI_kappa_results_within_subject/kappas_degree.csv"
kappa_df = pd.read_csv(real_kappa_path, sep=",")
real_participant_kappas = kappa_df["κ_Chronic_vs_Control"].dropna().values # drop NAs just in case
# null_mean = np.mean(null_kappa_means)
# extreme_count = np.sum(np.abs(null_kappa_means - null_mean) >= np.abs(real_kappa_mean - null_mean))
# p_value_two_sided = extreme_count / n_iterations
# print(f"Two-sided P-value: {p_value_two_sided:.4f}")
#plot
fig, ax1 = plt.subplots(figsize=(10, 5))
#histogram of permuted kappa means (null distribution)
color1 = "darkorange"
counts1, bins1, patches1 = ax1.hist(real_participant_kappas, bins=12, alpha=0.8, color=color1, edgecolor="black")
ax1.set_ylabel("Frequency (real participants)", color=color1, fontsize=20)
ax1.tick_params(axis='y', labelcolor=color1, labelsize=20)
ax1.set_ylim(0, 8)
ax1.tick_params(axis='x', labelsize=20)
#secondary axis for real participants' kappa distribution
ax2 = ax1.twinx()
color2 = "grey"
counts2, bins2, patches2 = ax2.hist(null_kappa_means, bins=20, alpha=0.6, color=color2, edgecolor="black")
ax2.axvline(real_kappa_mean, color='red', linestyle='dashed', linewidth=2, label=f"Real κ mean = {real_kappa_mean:.3f}")
ax2.set_xlabel("κ value")
ax2.set_ylabel("Frequency (permutations)", color=color2, fontsize=20)
ax2.tick_params(axis='y', labelcolor=color2, labelsize=20)
#add legends for both histograms
lines1, labels1 = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(lines1 + lines2, labels1 + labels2, loc="upper left", fontsize=14)
# # Add p-value text to plot
# ax1.text(0.95, 0.95, f'p = {p_value_two_sided:.4f}', transform=ax1.transAxes,
# fontsize=12, verticalalignment='top', horizontalalignment='right',
# bbox=dict(facecolor='white', alpha=0.8, edgecolor='black'))
plt.title("κ values distibution using degree - chronic vs control\nPermutation null ditribution vs Real participant values")
plt.tight_layout()
plt.show()
# %% [markdown]
# Closeness centrality
# %%
## the "mixed WITHIN PAIRS" aproach with the histogram of real values ##
#real closeness vectors from our csv files
def load_closeness_vector(subject_id, session):
folder = os.path.join(ses_paths[session], subject_id)
files = glob.glob(os.path.join(folder, metrics_pattern))
if not files:
raise FileNotFoundError(f"No metrics file found for {subject_id} in {session}")
df = pd.read_csv(files[0])
return df["Closeness"].values
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
ses_paths = {
"control": os.path.join(base_path, "ses-1"),
"chronic": os.path.join(base_path, "ses-3")
}
metrics_pattern = "Graphs/wAALours/graph_metrics/*_metrics.csv"
#subject IDs present in both sessions
ids_ses1 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["control"], "sub-*"))}
ids_ses3 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["chronic"], "sub-*"))}
common_ids = sorted(list(ids_ses1 & ids_ses3))
print(f"Subjects in both control and chronic: {len(common_ids)}")
#load degree vectors
closeness_control_all = []
closeness_chronic_all = []
for subject_id in common_ids:
closeness_control = load_closeness_vector(subject_id, "control")
closeness_chronic = load_closeness_vector(subject_id, "chronic")
closeness_control_all.append(closeness_control)
closeness_chronic_all.append(closeness_chronic)
closeness_control_all = np.array(closeness_control_all)
closeness_chronic_all = np.array(closeness_chronic_all)
#COMPUTE REAL KAPPAS
real_kappas = []
model = LinearRegression()
for i in range(len(common_ids)):
x = closeness_control_all[i].reshape(-1, 1)
y = closeness_chronic_all[i] - closeness_control_all[i]
model.fit(x, y)
real_kappas.append(model.coef_[0])
real_kappas = np.array(real_kappas)
real_kappa_mean = np.mean(real_kappas)
print(f"Real mean kappa (chronic vs control): {real_kappa_mean:.4f}")
#PERMUTATIONS - COMPUTE "NEW" KAPPAS
n_iterations = 10000
null_kappa_means = np.zeros(n_iterations)
for it in range(n_iterations):
perm_kappas = []
for i in range(len(common_ids)):
x = closeness_control_all[i].reshape(-1, 1)
delta = closeness_chronic_all[i] - closeness_control_all[i]
if np.random.rand() < 0.5:
delta = -delta
model.fit(x, delta)
perm_kappas.append(model.coef_[0])
null_kappa_means[it] = np.mean(perm_kappas)
# Load real participant kappa values
real_kappa_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/HDI_kappa_results_within_subject/kappas_closeness.csv"
kappa_df = pd.read_csv(real_kappa_path, sep=",") # Assuming it's tab-separated
real_participant_kappas = kappa_df["Chronic_vs_Control"].dropna().values # drop NAs just in case
# null_mean = np.mean(null_kappa_means)
# extreme_count = np.sum(np.abs(null_kappa_means - null_mean) >= np.abs(real_kappa_mean - null_mean))
# p_value_two_sided = extreme_count / n_iterations
# print(f"Two-sided P-value: {p_value_two_sided:.4f}")
# Create the plot
fig, ax1 = plt.subplots(figsize=(10, 5))
# Histogram of permuted kappa means (null distribution)
color1 = "darkorange"
counts1, bins1, patches1 = ax1.hist(real_participant_kappas, bins=12, alpha=0.8, color=color1, edgecolor="black")
ax1.set_ylabel("Frequency (real participants)", color=color1, fontsize=20)
ax1.tick_params(axis='y', labelcolor=color1, labelsize=20)
ax1.set_ylim(0, 8)
ax1.tick_params(axis='x', labelsize=20)
# Secondary axis for real participants' kappa distribution
ax2 = ax1.twinx()
color2 = "grey"
counts2, bins2, patches2 = ax2.hist(null_kappa_means, bins=20, alpha=0.6, color=color2, edgecolor="black")
ax2.axvline(real_kappa_mean, color='red', linestyle='dashed', linewidth=2, label=f"Real κ mean = {real_kappa_mean:.3f}")
ax2.set_xlabel("κ value")
ax2.set_ylabel("Frequency (permutations)", color=color2, fontsize=20)
ax2.tick_params(axis='y', labelcolor=color2, labelsize=20)
# Add legends for both histograms
lines1, labels1 = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(lines1 + lines2, labels1 + labels2, loc="upper left", fontsize=14)
# # Add p-value text to plot
# ax1.text(0.95, 0.95, f'p = {p_value_two_sided:.4f}', transform=ax1.transAxes,
# fontsize=12, verticalalignment='top', horizontalalignment='right',
# bbox=dict(facecolor='white', alpha=0.8, edgecolor='black'))
plt.title("κ values distibution using closeness - chronic vs control\nPermutation null ditribution vs Real participant values")
plt.tight_layout()
plt.show()
# %% [markdown]
# Clustering coefficient
# %%
###the "mixed WITHIN PAIRS" aproach with the histogram of real values ##
#real clustering vectors from our csv files
def load_clustering_vector(subject_id, session):
folder = os.path.join(ses_paths[session], subject_id)
files = glob.glob(os.path.join(folder, metrics_pattern))
if not files:
raise FileNotFoundError(f"No metrics file found for {subject_id} in {session}")
df = pd.read_csv(files[0])
return df["Clustering"].values
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
ses_paths = {
"control": os.path.join(base_path, "ses-1"),
"acute": os.path.join(base_path, "ses-2")
}
metrics_pattern = "Graphs/wAALours/graph_metrics/*_metrics.csv"
#subject IDs present in both sessions
ids_ses1 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["control"], "sub-*"))}
ids_ses2 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["acute"], "sub-*"))}
common_ids = sorted(list(ids_ses1 & ids_ses2))
print(f"Subjects in both control and acute: {len(common_ids)}")
#load degree vectors
clustering_control_all = []
clustering_acute_all = []
for subject_id in common_ids:
clustering_control = load_clustering_vector(subject_id, "control")
clustering_acute = load_clustering_vector(subject_id, "acute")
clustering_control_all.append(clustering_control)
clustering_acute_all.append(clustering_acute)
clustering_control_all = np.array(clustering_control_all)
clustering_acute_all = np.array(clustering_acute_all)
#COMPUTE REAL KAPPAS
real_kappas = []
model = LinearRegression()
for i in range(len(common_ids)):
x = clustering_control_all[i].reshape(-1, 1)
y = clustering_acute_all[i] - clustering_control_all[i]
model.fit(x, y)
real_kappas.append(model.coef_[0])
real_kappas = np.array(real_kappas)
real_kappa_mean = np.mean(real_kappas)
print(f"Real mean kappa (acute vs control): {real_kappa_mean:.4f}")
#PERMUTATIONS - COMPUTE "NEW" KAPPAS
n_iterations = 10000
null_kappa_means = np.zeros(n_iterations)
for it in range(n_iterations):
perm_kappas = []
for i in range(len(common_ids)):
x = clustering_control_all[i].reshape(-1, 1)
delta = clustering_acute_all[i] - clustering_control_all[i]
if np.random.rand() < 0.5:
delta = -delta
model.fit(x, delta)
perm_kappas.append(model.coef_[0])
null_kappa_means[it] = np.mean(perm_kappas)
# Load real participant kappa values
real_kappa_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/HDI_kappa_results_within_subject/kappas_clustering.csv"
kappa_df = pd.read_csv(real_kappa_path, sep=",")
real_participant_kappas = kappa_df["Acute_vs_Control"].dropna().values
# null_mean = np.mean(null_kappa_means)
# extreme_count = np.sum(np.abs(null_kappa_means - null_mean) >= np.abs(real_kappa_mean - null_mean))
# p_value_two_sided = extreme_count / n_iterations
# print(f"Two-sided P-value: {p_value_two_sided:.4f}")
# Create the plot
fig, ax1 = plt.subplots(figsize=(10, 5))
# Histogram of permuted kappa means (null distribution)
color1 = "darkorange"
counts1, bins1, patches1 = ax1.hist(real_participant_kappas, bins=12, alpha=0.8, color=color1, edgecolor="black")
ax1.set_ylabel("Frequency (real participants)", color=color1, fontsize=20)
ax1.tick_params(axis='y', labelcolor=color1, labelsize=20)
ax1.set_ylim(0, 8)
ax1.tick_params(axis='x', labelsize=20)
# Secondary axis for real participants' kappa distribution
ax2 = ax1.twinx()
color2 = "grey"
counts2, bins2, patches2 = ax2.hist(null_kappa_means, bins=20, alpha=0.6, color=color2, edgecolor="black")
ax2.axvline(real_kappa_mean, color='red', linestyle='dashed', linewidth=2, label=f"Real κ mean = {real_kappa_mean:.3f}")
ax2.set_xlabel("κ value")
ax2.set_ylabel("Frequency (permutations)", color=color2, fontsize=20)
ax2.tick_params(axis='y', labelcolor=color2, labelsize=20)
# Add legends for both histograms
lines1, labels1 = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(lines1 + lines2, labels1 + labels2, loc="upper left", fontsize=14)
# # Add p-value text to plot
# ax1.text(0.95, 0.95, f'p = {p_value_two_sided:.4f}', transform=ax1.transAxes,
# fontsize=12, verticalalignment='top', horizontalalignment='right',
# bbox=dict(facecolor='white', alpha=0.8, edgecolor='black'))
plt.title("κ values distibution using clustering - acute vs control\nPermutation null ditribution vs Real participant values")
plt.tight_layout()
plt.show()
# %%
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
sessions = {
"control": "ses-1",
"acute": "ses-2",
"chronic": "ses-3"
}
metrics_list = ["Degree_centrality", "Closeness", "Clustering"]
session_pairs = [
("acute", "control"),
("chronic", "control"),
("chronic", "acute")
]
# Custom labels for each comparison
comparison_labels = {
("acute", "control"): "TSD vs RW",
("chronic", "control"): "CSR vs RW",
("chronic", "acute"): "CSR vs TSD"
}
n_iterations = 10000
def load_metric_vector(subject_id, session, metric):
folder = os.path.join(base_path, sessions[session], subject_id)
pattern = os.path.join(folder, "Graphs/wAALours/graph_metrics/*_metrics.csv")
files = glob.glob(pattern)
if not files:
raise FileNotFoundError(f"No metrics file found for {subject_id} in {session}")
df = pd.read_csv(files[0])
return df[metric].values
# Initialize plot with metrics as rows, comparisons as columns
fig, axes = plt.subplots(len(metrics_list), len(session_pairs), figsize=(18, 12))
fig.subplots_adjust(hspace=0.4, wspace=0.3)
# Loop over all combinations
for col_idx, (ses1, ses_ref) in enumerate(session_pairs):
ids_ses1 = {os.path.basename(p) for p in glob.glob(os.path.join(base_path, sessions[ses1], "sub-*"))}
ids_ses_ref = {os.path.basename(p) for p in glob.glob(os.path.join(base_path, sessions[ses_ref], "sub-*"))}
common_ids = sorted(list(ids_ses1 & ids_ses_ref))
print(f"{ses1} vs {ses_ref}: {len(common_ids)} subjects")
for row_idx, metric in enumerate(metrics_list):
vec_ref_all = []
vec_test_all = []
for subject_id in common_ids:
vec_ref = load_metric_vector(subject_id, ses_ref, metric)
vec_test = load_metric_vector(subject_id, ses1, metric)
vec_ref_all.append(vec_ref)
vec_test_all.append(vec_test)
vec_ref_all = np.array(vec_ref_all)
vec_test_all = np.array(vec_test_all)
group_ref_mean = np.mean(vec_ref_all, axis=0)
# Real kappa differences
model = LinearRegression()
real_diffs = []
for i in range(len(common_ids)):
delta_within = vec_test_all[i] - vec_ref_all[i]
model.fit(vec_ref_all[i].reshape(-1, 1), delta_within)
kappa_within = model.coef_[0]
delta_group = vec_test_all[i] - group_ref_mean
model.fit(group_ref_mean.reshape(-1, 1), delta_group)
kappa_group = model.coef_[0]
real_diffs.append(kappa_within - kappa_group)
real_diffs = np.array(real_diffs)
real_mean_diff = np.mean(real_diffs)
# Permutations
null_diffs = np.zeros(n_iterations)
for it in range(n_iterations):
diffs = []
for i in range(len(common_ids)):
perm_test = np.random.permutation(vec_test_all[i])
delta_within = perm_test - vec_ref_all[i]
model.fit(vec_ref_all[i].reshape(-1, 1), delta_within)
kappa_within = model.coef_[0]
delta_group = perm_test - group_ref_mean
model.fit(group_ref_mean.reshape(-1, 1), delta_group)
kappa_group = model.coef_[0]
diffs.append(kappa_within - kappa_group)
null_diffs[it] = np.mean(diffs)
# P-value
p_value = np.mean(null_diffs <= real_mean_diff)
# Plot
ax = axes[row_idx, col_idx]
ax.hist(null_diffs, bins=30, alpha=0.7, color='grey')
ax.axvline(real_mean_diff, color='red', linestyle='dashed', linewidth=2,
label=f"Real = {real_mean_diff:.3f}\nP = {p_value:.4f}")
ax.set_title(f"{comparison_labels[(ses1, ses_ref)]}", fontsize=10)
if col_idx == 0:
ax.set_ylabel(f"{metric}\nFrequency")
else:
ax.set_ylabel("")
ax.set_xlabel("κ_within − κ_group")
ax.legend()
plt.suptitle("Permutation Test: κ_within − κ_groupmean", fontsize=16, y=1.02)
plt.tight_layout()
plt.show()
# %%
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
sessions = {
"control": "ses-1",
"acute": "ses-2",
"chronic": "ses-3"
}
metrics_list = ["Degree_centrality", "Closeness", "Clustering"]
session_pairs = [
("acute", "control"),
("chronic", "control"),
("chronic", "acute")
]
# Custom labels
comparison_labels = {
("acute", "control"): "TSD vs RW",
("chronic", "control"): "CSR vs RW",
("chronic", "acute"): "CSR vs TSD"
}
n_iterations = 10000
# Load metric vector
def load_metric_vector(subject_id, session, metric):
folder = os.path.join(base_path, sessions[session], subject_id)
pattern = os.path.join(folder, "Graphs/wAALours/graph_metrics/*_metrics.csv")
files = glob.glob(pattern)
if not files:
raise FileNotFoundError(f"No metrics file found for {subject_id} in {session}")
df = pd.read_csv(files[0])
return df[metric].values
# Create subplot: metrics = rows, comparisons = columns
fig, axes = plt.subplots(len(metrics_list), len(session_pairs), figsize=(18, 12))
fig.subplots_adjust(hspace=0.4, wspace=0.3)
# Loop over each comparison and metric
for col_idx, (ses1, ses_ref) in enumerate(session_pairs):
ids_ses1 = {os.path.basename(p) for p in glob.glob(os.path.join(base_path, sessions[ses1], "sub-*"))}
ids_ses_ref = {os.path.basename(p) for p in glob.glob(os.path.join(base_path, sessions[ses_ref], "sub-*"))}
common_ids = sorted(list(ids_ses1 & ids_ses_ref))
print(f"{ses1} vs {ses_ref}: {len(common_ids)} subjects")
for row_idx, metric in enumerate(metrics_list):
vec_ref_all = []
vec_test_all = []
for subject_id in common_ids:
vec_ref = load_metric_vector(subject_id, ses_ref, metric)
vec_test = load_metric_vector(subject_id, ses1, metric)
vec_ref_all.append(vec_ref)
vec_test_all.append(vec_test)
vec_ref_all = np.array(vec_ref_all)
vec_test_all = np.array(vec_test_all)
# Compute real kappas
model = LinearRegression()
real_kappas = []
for i in range(len(common_ids)):
x = vec_ref_all[i].reshape(-1, 1)
y = vec_test_all[i] - vec_ref_all[i]
model.fit(x, y)
real_kappas.append(model.coef_[0])
real_kappas = np.array(real_kappas)
real_kappa_mean = np.mean(real_kappas)
# Permutation: mix subjects
null_kappa_means = np.zeros(n_iterations)
for it in range(n_iterations):
perm_kappas = []
perm_indices = np.random.permutation(len(common_ids))
for i in range(len(common_ids)):
x = vec_ref_all[i].reshape(-1, 1)
y = vec_test_all[perm_indices[i]] - vec_ref_all[i]
model.fit(x, y)
perm_kappas.append(model.coef_[0])
null_kappa_means[it] = np.mean(perm_kappas)
# P-value
# p_value = np.mean(null_kappa_means <= real_kappa_mean)
# Two-tailed p-value
null_mean = np.mean(null_kappa_means)
p_value = np.mean(np.abs(null_kappa_means - null_mean) >= np.abs(real_kappa_mean - null_mean))
# Plotting
ax = axes[row_idx, col_idx]
ax.hist(null_kappa_means, bins=30, alpha=0.7, color='grey')
ax.axvline(real_kappa_mean, color='red', linestyle='dashed', linewidth=2,
label=f"Real κ = {real_kappa_mean:.3f}\nP = {p_value:.4f}")
ax.set_title(f"{comparison_labels[(ses1, ses_ref)]}", fontsize=10)
if col_idx == 0:
ax.set_ylabel(f"{metric}\nFrequency")
else:
ax.set_ylabel("")
ax.set_xlabel("Mean κ (shuffled)")
ax.legend()
plt.suptitle("Permutation Test: Mean κ (unpaired)", fontsize=16, y=1.02)
plt.tight_layout()
plt.show()
# %%
## the "mix the order of patients" aproach ##
#real degree vectors from our csv files
def load_degree_vector(subject_id, session):
folder = os.path.join(ses_paths[session], subject_id)
files = glob.glob(os.path.join(folder, metrics_pattern))
if not files:
raise FileNotFoundError(f"No metrics file found for {subject_id} in {session}")
df = pd.read_csv(files[0])
return df["Degree_centrality"].values
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
ses_paths = {
"control": os.path.join(base_path, "ses-1"),
"chronic": os.path.join(base_path, "ses-3")
}
metrics_pattern = "Graphs/wAALours/graph_metrics/*_metrics.csv"
#subject IDs present in both sessions
ids_ses1 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["control"], "sub-*"))}
ids_ses3 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["chronic"], "sub-*"))}
common_ids = sorted(list(ids_ses1 & ids_ses3))
print(f"Subjects in both control and chronic: {len(common_ids)}")
#load degree vectors
deg_control_all = []
deg_chronic_all = []
for subject_id in common_ids:
deg_control = load_degree_vector(subject_id, "control")
deg_chronic = load_degree_vector(subject_id, "chronic")
deg_control_all.append(deg_control)
deg_chronic_all.append(deg_chronic)
deg_control_all = np.array(deg_control_all)
deg_chronic_all = np.array(deg_chronic_all)
#COMPUTE REAL KAPPAS
real_kappas = []
model = LinearRegression()
for i in range(len(common_ids)):
x = deg_control_all[i].reshape(-1, 1)
y = deg_chronic_all[i] - deg_control_all[i]
model.fit(x, y)
real_kappas.append(model.coef_[0])
real_kappas = np.array(real_kappas)
real_kappa_mean = np.mean(real_kappas)
print(f"Real mean kappa (chronic vs control): {real_kappa_mean:.4f}")
#PERMUTATIONS - RANDOMLY BREAK SUBJECT PAIRINGS BETWEEN CHRONIC AND CONTROL
n_iterations = 10000
null_kappa_means = np.zeros(n_iterations)
for it in range(n_iterations):
perm_kappas = []
# random permutation of chronic subjects (break the subject pairing)
permuted_chronic_indices = np.random.permutation(len(common_ids))
for i in range(len(common_ids)):
x = deg_control_all[i].reshape(-1, 1)
y = deg_chronic_all[permuted_chronic_indices[i]] - deg_control_all[i]
model.fit(x, y)
perm_kappas.append(model.coef_[0])
null_kappa_means[it] = np.mean(perm_kappas)
#plot
plt.figure(figsize=(10, 5))
plt.hist(null_kappa_means, bins=30, alpha=0.7)
plt.axvline(real_kappa_mean, color='red', linestyle='dashed', linewidth=2,
label=f"Real mean κ = {real_kappa_mean:.3f}")
plt.xlabel("Mean kappa (permutation distribution)")
plt.ylabel("Frequency")
plt.title("Permutations - distribution of mean kappa (chronic vs control) using degree\n (controls in correct order, chronics random)")
plt.legend()
plt.tight_layout()
plt.show()
# %%
## the "mixed nodes" aproach ##
#real degree vectors from our csv files
def load_degree_vector(subject_id, session):
folder = os.path.join(ses_paths[session], subject_id)
files = glob.glob(os.path.join(folder, metrics_pattern))
if not files:
raise FileNotFoundError(f"No metrics file found for {subject_id} in {session}")
df = pd.read_csv(files[0])
return df["Degree_centrality"].values
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
ses_paths = {
"control": os.path.join(base_path, "ses-1"),
"chronic": os.path.join(base_path, "ses-3")
}
metrics_pattern = "Graphs/wAALours/graph_metrics/*_metrics.csv"
# subject IDs present in both sessions
ids_ses1 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["control"], "sub-*"))}
ids_ses3 = {os.path.basename(p) for p in glob.glob(os.path.join(ses_paths["chronic"], "sub-*"))}
common_ids = sorted(list(ids_ses1 & ids_ses3))
print(f"Subjects in both control and chronic: {len(common_ids)}")
# load degree vectors
deg_control_all = []
deg_chronic_all = []
for subject_id in common_ids:
deg_control = load_degree_vector(subject_id, "control")
deg_chronic = load_degree_vector(subject_id, "chronic")
deg_control_all.append(deg_control)
deg_chronic_all.append(deg_chronic)
deg_control_all = np.array(deg_control_all)
deg_chronic_all = np.array(deg_chronic_all)
# COMPUTE REAL KAPPAS
real_kappas = []
model = LinearRegression()
for i in range(len(common_ids)):
x = deg_control_all[i].reshape(-1, 1)
y = deg_chronic_all[i] - deg_control_all[i]
model.fit(x, y)
real_kappas.append(model.coef_[0])
real_kappas = np.array(real_kappas)
real_kappa_mean = np.mean(real_kappas)
print(f"Real mean kappa (chronic vs control): {real_kappa_mean:.4f}")
# PERMUTATIONS - shuffle node order in ONE of the vectors (e.g., chronic)
n_iterations = 10000
null_kappa_means = np.zeros(n_iterations)
for it in range(n_iterations):
perm_kappas = []
for i in range(len(common_ids)):
x = deg_control_all[i] # do not shuffle
y = deg_chronic_all[i].copy()
#shuffle node order in chronic session
y_shuffled = np.random.permutation(y)
delta = y_shuffled - x #mismatch in node correspondence, break the pairs
model.fit(x.reshape(-1, 1), delta)
perm_kappas.append(model.coef_[0])
null_kappa_means[it] = np.mean(perm_kappas)
# plot
plt.figure(figsize=(10, 5))
plt.hist(null_kappa_means, bins=30, alpha=0.7)
plt.axvline(real_kappa_mean, color='red', linestyle='dashed', linewidth=2,
label=f"Real κ = {real_kappa_mean:.3f}")
plt.xlabel("Mean kappa (permutation distribution)")
plt.ylabel("Frequency")
plt.title("Permutations - distribution of mean kappa (chronic vs control) using degree\n(node-wise shuffling WITHIN pairs)")
plt.legend()
plt.tight_layout()
plt.show()
# %% [markdown]
# ## Covariate constraint manifold learning
# %% [markdown]
# CCML degree centrality
# %%
metric = "Degree_centrality"
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
metric_filename_pattern = "*_metrics.csv"
#kappas file
kappa_all = pd.read_csv(os.path.join(base_path, "HDI_kappa_results_within_subject", "kappas_degree.csv"))
kappa_all.set_index("Subject_ID", inplace=True)
#subject directories
control_dir = os.path.join(base_path, "ses-2")
patient_dir = os.path.join(base_path, "ses-3")
control_subjects = sorted(glob.glob(os.path.join(control_dir, "sub-*")))
patient_subjects = sorted(glob.glob(os.path.join(patient_dir, "sub-*")))
#extract the metric for each subject
def extract_metric(subject_paths):
data = []
ids = []
for subj_path in subject_paths:
pattern = os.path.join(subj_path, "Graphs/wAALours/graph_metrics", metric_filename_pattern)
files = glob.glob(pattern)
if not files:
continue
df = pd.read_csv(files[0])
if metric not in df.columns:
continue
data.append(df[metric].values)
ids.append(os.path.basename(subj_path))
return np.array(data), ids
#extract data
control_data, control_ids = extract_metric(control_subjects)
patient_data, patient_ids = extract_metric(patient_subjects)
#just in case - clean subject IDs to match kappa file format
control_ids_cleaned = [cid.replace("sub-", "") for cid in control_ids]
patient_ids_cleaned = [pid.replace("sub-", "") for pid in patient_ids]
#extract matching kappa values
kappa_control_values = kappa_all.loc[["sub-" + cid for cid in control_ids_cleaned], "κ_Acute_vs_Control"].values
kappa_patient_values = kappa_all.loc[["sub-" + pid for pid in patient_ids_cleaned], "κ_Chronic_vs_Control"].values
#construct covariate vector
cov = np.concatenate([kappa_control_values, kappa_patient_values])
#group labels: 0 = control, 1 = patient
labels = [0] * len(control_ids) + [1] * len(patient_ids)
#subject names
subject_names = [f"T{i}" for i in range(len(control_ids))] + [f"C{i}" for i in range(len(patient_ids))]
#combine metric data
X = np.vstack([control_data, patient_data])
#save outputs for CCML (exactly the same format as in Sophie's tutorial code)
output_dir = os.path.join(base_path, "CCML_inputs_within_subjects")
os.makedirs(output_dir, exist_ok=True)
df_X = pd.DataFrame(X, index=subject_names)
df_X.to_csv(os.path.join(output_dir, f"{metric}_chronic_and_acute_vs_control.csv"))
df_cov = pd.DataFrame(cov, index=subject_names, columns=["HDI_kappa"])
df_cov.to_csv(os.path.join(output_dir, f"HDI_kappa_{metric}_chronic_and_acute_vs_control.csv"))
df_labels = pd.DataFrame(labels, index=subject_names, columns=["Group"])
df_labels.to_csv(os.path.join(output_dir, f"group_labels_{metric}_chronic_and_acute_vs_control.csv"))
print("Files for CCML saved in:", output_dir)
# %% [markdown]
# %%
def f_coord1(x,X2,cov2,alpha2,dist2):
X2 = np.vstack([X2,x])
cov2 = cov2.reshape(X2.shape)
X_tmp2 = np.hstack((alpha2*cov2,X2))
D2 = sk.metrics.pairwise_distances(X_tmp2)
return np.linalg.norm(dist2-D2)
def f_glob(X2,cov2,alpha2,dist2):
cov2 = cov2.reshape(X2.shape)
X_tmp2 = np.vstack((alpha2*cov2,X2)).T
D2 = sk.metrics.pairwise_distances(X_tmp2)
return np.linalg.norm(dist2-D2)
def f_alpha(alpha,X2,cov2,dist2):
cov2 = cov2.reshape(X2.shape)
X_tmp = np.vstack((alpha*cov2,X2))
D2 = sk.metrics.pairwise_distances(X_tmp.T)
return np.linalg.norm(dist2-D2)
dfX = pd.read_csv("/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/CCML_inputs_within_subjects/Degree_centrality_chronic_and_acute_vs_control.csv", index_col=0)
X = dfX.to_numpy()
dfcov = pd.read_csv("/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/CCML_inputs_within_subjects/HDI_kappa_Degree_centrality_chronic_and_acute_vs_control.csv", index_col=0)
cov = dfcov.to_numpy()
labels_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/CCML_inputs_within_subjects/group_labels_Degree_centrality_chronic_and_acute_vs_control.csv"
labels_df = pd.read_csv(labels_path, index_col=0)
Indiv = labels_df.index.tolist()
l = labels_df["Group"].values
#compute Euclidean distance matrix
D = sk.metrics.pairwise.pairwise_distances(X)
#initialisation aleatoire
#list_extract = sk.utils.shuffle(np.arange(X.shape[0]))
#initialize with the most remote couple of points
ll = []
ll.append(np.where(D==np.max(D))[0][0])
ll.append(np.where(D==np.max(D))[1][0])
# Add sequentially the furthest point of the previous set
while(len(ll)!=D.shape[0]):
D_tmp = D[ll]
D_tmp[:,ll] = 0
ll.append(np.where(D_tmp == np.max(D_tmp))[1][0])
list_extract = np.array(ll)
# Initialize alpha parameter with the classic Isomap
iso = sk.manifold.Isomap(n_components=2,n_neighbors=4)
iso.fit(X)
tmp = iso.embedding_[:,0]
dist = iso.dist_matrix_
delta_1 = np.max(tmp) -np.min(tmp)
delta_2 = np.max(cov) -np.min(cov)
alpha = delta_1/delta_2
#Initialize the first point at 0
X_tmp = np.zeros([1])
for i in range(1,X.shape[0]):
print("one coord",i)
cov_tmp = cov[list_extract[0:i+1]]
dist_tmp = dist[list_extract[0:i+1]][:,list_extract[0:i+1]]
opt_test = 10000
#Multi-start opitmisation
for j in range(-5,5):
xtmptmp,B,tmp,tmp1,tmp3 = opt.fmin(f_coord1,j,(X_tmp,cov_tmp,alpha,dist_tmp),xtol=0.000000001, ftol=0.000000001,maxiter=1000000,maxfun=1000000,disp = 0 , full_output = 1)
if opt_test>B:
xtmp = xtmptmp
opt_test = B
print(j)
X_tmp = np.vstack([X_tmp, xtmp])
# We can add an additional global optimisation before the estimation of the alpha coefficient
# X_tmp2,B,tmp,tmp1,tmp3 = opt.fmin(f_glob,X_tmp,(cov_tmp,alpha,dist_tmp), xtol=0.000000001, ftol=0.000000001,maxiter=1000000,maxfun=1000000,disp = 1 , full_output = 1)
# print "all_coord"
# X_tmp = X_tmp2.reshape(-1,1)
alpha,B,tmp,tmp1,tmp3 = opt.fmin(f_alpha,alpha,(X_tmp,cov_tmp,dist_tmp), xtol=0.000000001, ftol=0.00000001,disp = 1 , full_output = 1)
X_tmp2,B,tmp,tmp1,tmp3 = opt.fmin(f_glob,X_tmp,(cov_tmp,alpha,dist_tmp), xtol=0.000000000001, ftol=0.00000000001,maxiter=1000000,maxfun=1000000,disp = 1 , full_output = 1)
print("all_coord")
X_tmp = X_tmp2.reshape(-1,1)
# Merge the covariate modulated by alpha and the first component
cov2 = cov.reshape(X_tmp.shape)
X_iso = np.hstack((alpha*cov2[list_extract],X_tmp))
# Re-order the rows of the matrix to correspond to the initial order of the subject
I = np.argsort(list_extract)
X_iso_f = X_iso[I]
import sys
sys.path.append('/Users/patrycjascislewska/Analizy_neuro/Graphs/tutorials/TP2_Graphs_ILCB_summer_School')
import PlotFunction as pf
###### Visualization
Xt = cov
scaling = 5
grid_x , grid_y = np.mgrid[-1:1:1000j,-1:1:1000j]*scaling
grid_lin = griddata(X_iso_f,Xt,(grid_x,grid_y),method='linear')
grid_lin = grid_lin.reshape(1000,1000)
pf.scatter_2D(X_iso_f, l,Indiv)
#plt.imshow(grid_lin.T, extent=(-1,1,-1,1), origin='lower')
plt.imshow(grid_lin.T,extent=(-scaling,scaling,-scaling,scaling), origin='lower')
plt.colorbar().ax.tick_params(labelsize=20)
plt.title('Degree: TSD and CSR compared to RW', fontsize=16)
plt.tick_params(labelsize=20)
# %%
# calculation of difference between two clusters in CCML - for acute and chronic
group_acute = labels_df["Group"] == 0
group_chronic = labels_df["Group"] == 1
#extract group coordinates
X_acute = X_iso_f[group_acute.values]
X_chronic = X_iso_f[group_chronic.values]
#calculate observed centroid distance - distance between two clusters
d_obs = np.linalg.norm(X_acute.mean(axis=0) - X_chronic.mean(axis=0))
#permutations
n_perm = 10000
d_perm = []
group_labels = labels_df["Group"].values
for _ in range(n_perm):
perm_labels = np.random.permutation(group_labels) #randomly mix the labels of acutes and chronics
perm_acute = X_iso_f[perm_labels == 0]
perm_chronic = X_iso_f[perm_labels == 1]
d = np.linalg.norm(perm_acute.mean(axis=0) - perm_chronic.mean(axis=0))
d_perm.append(d)
#p-value
p_val = np.mean(np.array(d_perm) >= d_obs) #greater than or equal to the observed distance
print(f"Observed centroid distance = {d_obs:.4f}")
print(f"P-value (permutation test) = {p_val:.4f}")
#histogram of permuted distances
plt.figure(figsize=(8, 5))
plt.hist(d_perm, bins=50, edgecolor='black', alpha=0.7, color='grey')
plt.axvline(d_obs, color='red', linestyle='--', linewidth=2, label=f'Observed distance = {d_obs:.3f}\np-value = {p_val:.4f}')
plt.title("Permutation test - centroid distance", fontsize=20)
plt.xlabel("Centroid distance", fontsize=18)
plt.ylabel("Frequency", fontsize=18)
plt.legend(fontsize=14)
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()
# %% [markdown]
# CCML closeness centrality
# %%
metric = "Closeness"
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
metric_filename_pattern = "*_metrics.csv"
#kappas file
kappa_all = pd.read_csv(os.path.join(base_path, "HDI_kappa_results_within_subject", "kappas_closeness.csv"))
kappa_all.set_index("Subject_ID", inplace=True)
#subject directories
control_dir = os.path.join(base_path, "ses-2")
patient_dir = os.path.join(base_path, "ses-3")
control_subjects = sorted(glob.glob(os.path.join(control_dir, "sub-*")))
patient_subjects = sorted(glob.glob(os.path.join(patient_dir, "sub-*")))
#extract the metric for each subject
def extract_metric(subject_paths):
data = []
ids = []
for subj_path in subject_paths:
pattern = os.path.join(subj_path, "Graphs/wAALours/graph_metrics", metric_filename_pattern)
files = glob.glob(pattern)
if not files:
continue
df = pd.read_csv(files[0])
if metric not in df.columns:
continue
data.append(df[metric].values)
ids.append(os.path.basename(subj_path))
return np.array(data), ids
#extract data
control_data, control_ids = extract_metric(control_subjects)
patient_data, patient_ids = extract_metric(patient_subjects)
#just in case - clean subject IDs to match kappa file format
control_ids_cleaned = [cid.replace("sub-", "") for cid in control_ids]
patient_ids_cleaned = [pid.replace("sub-", "") for pid in patient_ids]
#extract matching kappa values
kappa_control_values = kappa_all.loc[["sub-" + cid for cid in control_ids_cleaned], "Acute_vs_Control"].values
kappa_patient_values = kappa_all.loc[["sub-" + pid for pid in patient_ids_cleaned], "Chronic_vs_Control"].values
#construct covariate vector
cov = np.concatenate([kappa_control_values, kappa_patient_values])
#group labels: 0 = control, 1 = patient
labels = [0] * len(control_ids) + [1] * len(patient_ids)
#subject names
subject_names = [f"T{i}" for i in range(len(control_ids))] + [f"C{i}" for i in range(len(patient_ids))]
#combine metric data
X = np.vstack([control_data, patient_data])
#save outputs for CCML (exactly the same format as in Sophie's tutorial code)
output_dir = os.path.join(base_path, "CCML_inputs_within_subjects")
os.makedirs(output_dir, exist_ok=True)
df_X = pd.DataFrame(X, index=subject_names)
df_X.to_csv(os.path.join(output_dir, f"{metric}_chronic_and_acute_vs_control.csv"))
df_cov = pd.DataFrame(cov, index=subject_names, columns=["HDI_kappa"])
df_cov.to_csv(os.path.join(output_dir, f"HDI_kappa_{metric}_chronic_and_acute_vs_control.csv"))
df_labels = pd.DataFrame(labels, index=subject_names, columns=["Group"])
df_labels.to_csv(os.path.join(output_dir, f"group_labels_{metric}_chronic_and_acute_vs_control.csv"))
print("Files for CCML saved in:", output_dir)
# %%
def f_coord1(x,X2,cov2,alpha2,dist2):
X2 = np.vstack([X2,x])
cov2 = cov2.reshape(X2.shape)
X_tmp2 = np.hstack((alpha2*cov2,X2))
D2 = sk.metrics.pairwise_distances(X_tmp2)
return np.linalg.norm(dist2-D2)
def f_glob(X2,cov2,alpha2,dist2):
cov2 = cov2.reshape(X2.shape)
X_tmp2 = np.vstack((alpha2*cov2,X2)).T
D2 = sk.metrics.pairwise_distances(X_tmp2)
return np.linalg.norm(dist2-D2)
def f_alpha(alpha,X2,cov2,dist2):
cov2 = cov2.reshape(X2.shape)
X_tmp = np.vstack((alpha*cov2,X2))
D2 = sk.metrics.pairwise_distances(X_tmp.T)
return np.linalg.norm(dist2-D2)
dfX = pd.read_csv("/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/CCML_inputs_within_subjects/Closeness_chronic_and_acute_vs_control.csv", index_col=0)
X = dfX.to_numpy()
dfcov = pd.read_csv("/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/CCML_inputs_within_subjects/HDI_kappa_Closeness_chronic_and_acute_vs_control.csv", index_col=0)
cov = dfcov.to_numpy()
labels_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/CCML_inputs_within_subjects/group_labels_Closeness_chronic_and_acute_vs_control.csv"
labels_df = pd.read_csv(labels_path, index_col=0)
Indiv = labels_df.index.tolist()
l = labels_df["Group"].values
#compute Euclidean distance matrix
D = sk.metrics.pairwise.pairwise_distances(X)
#initialisation aleatoire
#list_extract = sk.utils.shuffle(np.arange(X.shape[0]))
#initialize with the most remote couple of points
ll = []
ll.append(np.where(D==np.max(D))[0][0])
ll.append(np.where(D==np.max(D))[1][0])
# Add sequentially the furthest point of the previous set
while(len(ll)!=D.shape[0]):
D_tmp = D[ll]
D_tmp[:,ll] = 0
ll.append(np.where(D_tmp == np.max(D_tmp))[1][0])
list_extract = np.array(ll)
# Initialize alpha parameter with the classic Isomap
iso = sk.manifold.Isomap(n_components=2,n_neighbors=4)
iso.fit(X)
tmp = iso.embedding_[:,0]
dist = iso.dist_matrix_
delta_1 = np.max(tmp) -np.min(tmp)
delta_2 = np.max(cov) -np.min(cov)
alpha = delta_1/delta_2
#Initialize the first point at 0
X_tmp = np.zeros([1])
for i in range(1,X.shape[0]):
print("one coord",i)
cov_tmp = cov[list_extract[0:i+1]]
dist_tmp = dist[list_extract[0:i+1]][:,list_extract[0:i+1]]
opt_test = 10000
#Multi-start opitmisation
for j in range(-5,5):
xtmptmp,B,tmp,tmp1,tmp3 = opt.fmin(f_coord1,j,(X_tmp,cov_tmp,alpha,dist_tmp),xtol=0.000000001, ftol=0.000000001,maxiter=1000000,maxfun=1000000,disp = 0 , full_output = 1)
if opt_test>B:
xtmp = xtmptmp
opt_test = B
print(j)
X_tmp = np.vstack([X_tmp, xtmp])
# We ca add an additional global optimisation before the estimation of the alpha coefficient
# X_tmp2,B,tmp,tmp1,tmp3 = opt.fmin(f_glob,X_tmp,(cov_tmp,alpha,dist_tmp), xtol=0.000000001, ftol=0.000000001,maxiter=1000000,maxfun=1000000,disp = 1 , full_output = 1)
# print "all_coord"
# X_tmp = X_tmp2.reshape(-1,1)
alpha,B,tmp,tmp1,tmp3 = opt.fmin(f_alpha,alpha,(X_tmp,cov_tmp,dist_tmp), xtol=0.000000001, ftol=0.00000001,disp = 1 , full_output = 1)
X_tmp2,B,tmp,tmp1,tmp3 = opt.fmin(f_glob,X_tmp,(cov_tmp,alpha,dist_tmp), xtol=0.000000000001, ftol=0.00000000001,maxiter=1000000,maxfun=1000000,disp = 1 , full_output = 1)
print("all_coord")
X_tmp = X_tmp2.reshape(-1,1)
# Merge the covariate modulated by alpha and the first component
cov2 = cov.reshape(X_tmp.shape)
X_iso = np.hstack((alpha*cov2[list_extract],X_tmp))
# Re-order the rows of the matrix to correspond to the initial order of the subject
I = np.argsort(list_extract)
X_iso_f = X_iso[I]
import sys
sys.path.append('/Users/patrycjascislewska/Analizy_neuro/Graphs/tutorials/TP2_Graphs_ILCB_summer_School')
import PlotFunction as pf
###### Visualization
Xt = cov
scaling = 5
grid_x , grid_y = np.mgrid[-1:1:1000j,-1:1:1000j]*scaling
grid_lin = griddata(X_iso_f,Xt,(grid_x,grid_y),method='linear')
grid_lin = grid_lin.reshape(1000,1000)
pf.scatter_2D(X_iso_f, l,Indiv)
#plt.imshow(grid_lin.T, extent=(-1,1,-1,1), origin='lower')
plt.imshow(grid_lin.T,extent=(-scaling,scaling,-scaling,scaling), origin='lower')
plt.colorbar().ax.tick_params(labelsize=20)
plt.title('Closeness: CSR and TSD compared to RW')
plt.tick_params(labelsize=20)
# %% [markdown]
# CCML clustering coefficient
# %%
metric = "Clustering"
base_path = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
metric_filename_pattern = "*_metrics.csv"
#kappas file
kappa_all = pd.read_csv(os.path.join(base_path, "HDI_kappa_results_within_subject", "kappas_clustering.csv"))
kappa_all.set_index("Subject_ID", inplace=True)
#subject directories
control_dir = os.path.join(base_path, "ses-2")
patient_dir = os.path.join(base_path, "ses-3")
control_subjects = sorted(glob.glob(os.path.join(control_dir, "sub-*")))
patient_subjects = sorted(glob.glob(os.path.join(patient_dir, "sub-*")))
#extract the metric for each subject
def extract_metric(subject_paths):
data = []
ids = []
for subj_path in subject_paths:
pattern = os.path.join(subj_path, "Graphs/wAALours/graph_metrics", metric_filename_pattern)
files = glob.glob(pattern)
if not files:
continue
df = pd.read_csv(files[0])
if metric not in df.columns:
continue
data.append(df[metric].values)
ids.append(os.path.basename(subj_path))
return np.array(data), ids
#extract data
control_data, control_ids = extract_metric(control_subjects)
patient_data, patient_ids = extract_metric(patient_subjects)
#just in case - clean subject IDs to match kappa file format
control_ids_cleaned = [cid.replace("sub-", "") for cid in control_ids]
patient_ids_cleaned = [pid.replace("sub-", "") for pid in patient_ids]
#extract matching kappa values
kappa_control_values = kappa_all.loc[["sub-" + cid for cid in control_ids_cleaned], "Acute_vs_Control"].values
kappa_patient_values = kappa_all.loc[["sub-" + pid for pid in patient_ids_cleaned], "Chronic_vs_Control"].values
#construct covariate vector
cov = np.concatenate([kappa_control_values, kappa_patient_values])
#group labels: 0 = control, 1 = patient
labels = [0] * len(control_ids) + [1] * len(patient_ids)
#subject names
subject_names = [f"A{i}" for i in range(len(control_ids))] + [f"Ch{i}" for i in range(len(patient_ids))]
#combine metric data
X = np.vstack([control_data, patient_data])
#save outputs for CCML (exactly the same format as in Sophie's tutorial code)
output_dir = os.path.join(base_path, "CCML_inputs_within_subjects")
os.makedirs(output_dir, exist_ok=True)
df_X = pd.DataFrame(X, index=subject_names)
df_X.to_csv(os.path.join(output_dir, f"{metric}_chronic_and_acute_vs_control.csv"))
df_cov = pd.DataFrame(cov, index=subject_names, columns=["HDI_kappa"])
df_cov.to_csv(os.path.join(output_dir, f"HDI_kappa_{metric}_chronic_and_acute_vs_control.csv"))
df_labels = pd.DataFrame(labels, index=subject_names, columns=["Group"])
df_labels.to_csv(os.path.join(output_dir, f"group_labels_{metric}_chronic_and_acute_vs_control.csv"))
print("Files for CCML saved in:", output_dir)
# %% [markdown]
# ## Chord plots
# %%
#### individual subjects, 200 edges (from adjacency matrix), size of node corresponds to the number of connecitons
#### color coding = networks
import holoviews as hv
from holoviews import opts
hv.extension('bokeh')
base_dir = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA"
roi_info_csv = "/Users/patrycjascislewska/Analizy_neuro/Graphs/AAL_ours_89_regions_list.csv"
out_dir = "/Users/patrycjascislewska/Analizy_neuro/Graphs/DATA/graph_metrics/chord_diagrams"
sessions = {
"control": "ses-1_AICHA_atlas",
"acute": "ses-2_AICHA_atlas",
"chronic": "ses-3_AICHA_atlas"
}
##ROI TO NETWORK MAPPING
roi_networks = {
'PreGy_L': 'SMN', 'PreGy_R': 'SMN', 'SMA_L': 'SMN', 'SMA_R': 'SMN',
'RolandOperc_L': 'SMN', 'RolandOperc_R': 'SMN', 'Postcentral_L': 'SMN', 'Postcentral_R': 'SMN',
'ParacentralLob_L': 'SMN', 'ParacentralLob_R': 'SMN',
'FrontSup_L': 'FPN', 'FrontSup_R': 'FPN', 'FrontMid_L': 'FPN', 'FrontMid_R': 'FPN',
'FrontInfOperc_L': 'FPN', 'FrontInfOperc_R': 'FPN', 'FrontInfTri_L': 'FPN', 'FrontInfTri_R': 'FPN',
'FrontInfOrb_L': 'FPN', 'FrontInfOrb_R': 'FPN', 'FrontSupMed_L': 'FPN', 'FrontSupMed_R': 'FPN',
'CingAnt_L': 'DMN', 'CingAnt_R': 'DMN', 'CingMid_L': 'DMN', 'CingMid_R': 'DMN',
'CingPost_L': 'DMN', 'CingPost_R': 'DMN', 'Precuneus_L': 'DMN', 'Precuneus_R': 'DMN',
'Angular_L': 'DMN', 'Angular_R': 'DMN', 'FrontMedOrb_L': 'DMN', 'FrontMedOrb_R': 'DMN',
'TempMid_L': 'DMN', 'TempMid_R': 'DMN',
'Insula_L': 'Salience', 'Insula_R': 'Salience', 'Olfactory_L': 'Salience', 'Olfactory_R': 'Salience',
'Calcarine_L': 'VN', 'Calcarine_R': 'VN', 'Cuneus_L': 'VN', 'Cuneus_R': 'VN',
'Lingual_L': 'VN', 'Lingual_R': 'VN', 'Occipital_L': 'VN', 'Occipial_R': 'VN',
'Fusiform_L': 'VN', 'Fusiform_R': 'VN',
'Amygdala_L': 'Limbic', 'Amygdala_R': 'Limbic', 'Hippocampus_L': 'Limbic', 'Hippocampus_R': 'Limbic',
'ParaHippoc_L': 'Limbic', 'ParaHippoc_R': 'Limbic', 'TempPole_L': 'Limbic', 'TempPole_R': 'Limbic',
'Heschl_L': 'Auditory', 'Heschl_R': 'Auditory', 'TempSup_L': 'Auditory', 'TempSup_R': 'Auditory',
'TempInf_L': 'Language', 'TempInf_R': 'Language',
'Caudate_L': 'Subcortical', 'Caudate_R': 'Subcortical', 'Putamen_L': 'Subcortical',
'Putamen_R': 'Subcortical', 'Pallidum_L': 'Subcortical', 'Pallidum_R': 'Subcortical',
'Thalamus_L': 'Subcortical', 'Thalamus_R': 'Subcortical',
'ParietalSup_L': 'DAN', 'ParietalSup_R': 'DAN', 'ParietalInf_L': 'DAN', 'ParietalInf_R': 'DAN',
'SupraMarginal_L': 'DAN', 'SupraMarginal_R': 'DAN',
'Cereb_I_II_L': 'Cerebellum', 'Cereb_I_II_R': 'Cerebellum',
'Cereb_III_VI_L': 'Cerebellum', 'Cereb_III_VI_R': 'Cerebellum',
'Cereb_VII_X_L': 'Cerebellum', 'Cereb_VII_X_R': 'Cerebellum',
'Vermis': 'Cerebellum',
'FrontSupOrb_L': 'Other', 'FrontSupOrb_R': 'Other',
'FrontMidOrb_L': 'Other', 'FrontMidOrb_R': 'Other'
}
#ROI names and sort by atlas Node_number
roi_df = pd.read_csv(roi_info_csv, sep=";")
roi_df = roi_df.sort_values("Node_number")
original_region_names = roi_df["Region"].tolist()
nodes_df = pd.DataFrame({
'name': original_region_names,
'network': [roi_networks.get(roi, 'Other') for roi in original_region_names]
})
nodes_df = nodes_df.sort_values(by=['network', 'name']).reset_index(drop=True)
region_names = nodes_df['name'].tolist()
n_rois = len(region_names)
def find_all_participants():
participants = set()
for ses_folder in sessions.values():
ses_path = os.path.join(base_dir, ses_folder)
for root, dirs, files in os.walk(ses_path):
for file in files:
if file.startswith("Adj_mat") and file.endswith("200.txt"):
parts = root.split(os.sep)
for p in parts:
if p.startswith("sub-"):
participants.add(p)
return sorted(participants)
#I used adj matrix with 200 edges from R script from Veronica, because chart with 400 edges was to dense
def get_matrix(sub_id, session_folder):
path = os.path.join(base_dir, session_folder, sub_id, "Graphs", "wAALours")
try:
files = [f for f in os.listdir(path) if f.startswith("Adj_mat") and f.endswith("200.txt")]
if not files:
return None
matrix = np.loadtxt(os.path.join(path, files[0]))
if matrix.shape != (n_rois, n_rois):
return None
return matrix.astype(int)
except:
return None
def create_chord(matrix, title):
reordered_indices = [original_region_names.index(name) for name in region_names]
matrix = matrix[np.ix_(reordered_indices, reordered_indices)]
edges = []
connected_nodes = set()
for i in range(matrix.shape[0]):
for j in range(i + 1, matrix.shape[1]):
if matrix[i, j] == 1:
source = region_names[i]
target = region_names[j]
network_src = roi_networks.get(source, 'Other')
network_tgt = roi_networks.get(target, 'Other')
edges.append((source, target, 1, network_src))
connected_nodes.add(source)
connected_nodes.add(target)
for name in region_names:
if name not in connected_nodes:
edges.append((name, name, 0.01, roi_networks.get(name, 'Other')))
edges_df = pd.DataFrame(edges, columns=["source", "target", "value", "network_src"])
connection_counts = edges_df[['source', 'target']].stack().value_counts()
nodes_with_size = nodes_df.copy()
nodes_with_size['connections'] = nodes_with_size['name'].map(connection_counts).fillna(0)
nodes_with_size['size'] = nodes_with_size['connections'] * 4
chord = hv.Chord((edges_df, hv.Dataset(nodes_with_size, kdims='name')))
return chord.opts(
opts.Chord(
labels='name',
node_color='network',
edge_color='network_src',
cmap='Set1',
edge_cmap='Set1',
edge_alpha=0.8,
node_size='size',
width=650,
height=650,
title=title,
colorbar=False,
tools=['hover']
)
)
#MAIN
participants = find_all_participants()
print(f"Found {len(participants)} participants")
os.makedirs(out_dir, exist_ok=True)
for sub_id in participants:
chords = []
for condition, session_folder in sessions.items():
mat = get_matrix(sub_id, session_folder)
if mat is not None:
title = f"{condition.capitalize()} – {sub_id}"
chords.append(create_chord(mat, title))
else:
print(f"Missing or invalid matrix for {sub_id} in {session_folder}")
if chords:
layout = hv.Layout(chords).cols(3)
filename = os.path.join(out_dir, f"{sub_id}_chord_diagram_all_sessions_200_sized_colored.html")
hv.save(layout, filename, backend='bokeh')
print(f"Saved {filename}")
else:
print(f"No valid matrices for {sub_id}, skipping.")
print("All participant-level multi-session chord diagrams saved.")

Graphs_nodal_global_metrics_HDI_CCML.ipynb at commit af790e1, no license · at the source

Overview

Authors: Patrycja Scislewska¹, Arturo Cabrera Vazquez^2,3,4, Iwona Szatkowska⁵, Halszka Kontrymowicz-Ogińska⁶, Sophie Achard², Aleksandra Domagalik⁷

ORCID iDs: Patrycja Scislewska, Arturo Cabrera Vazquez, Iwona Szatkowska, Sophie Achard, Aleksandra Domagalik

Institutions: University of Warsaw (Poland); Institut polytechnique de Grenoble (France); Centre National de la Recherche Scientifique (France); Inserm (France); Grenoble Institute of Neurosciences (France); Toulouse NeuroImaging Center (France); Université Grenoble Alpes (France); Instytut Biologii Doświadczalnej im. Marcelego Nenckiego (Poland); Polish Academy of Sciences (Poland); Jagiellonian University (Poland)

Journal: Sleep, volume 49, issue 5, article zsag030

Dates: received 20 October 2025; accepted 29 January 2026; published online 3 February 2026; in print May 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1093/sleep/zsag030 · PMID 41631633 · PMCID PMC13163182 · OpenAlex W7127334101

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: fMRI (modality), human (organism), systems (subfield)

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

Keywords: fMRI, sleep deprivation, graphs, functional connectivity, hub disruption index, machine learning, sleepiness, chronotype

MeSH: Brain*, Magnetic Resonance Imaging*, Nerve Net*, Sleep Deprivation*, Adult, Brain Mapping, Female, Humans, Longitudinal Studies, Male, Wakefulness, Young Adult (* major topic)

Topic: Sleep and Work-Related Fatigue (Experimental and Cognitive Psychology, Psychology), according to OpenAlex

Funding: Excellence Initiative-Research University (BOB-IDUB-622-412/2025, 2020-2026); Inria and including CNRS; Agence Nationale de la Recherche under the France 2030 program (ANR-23-IACL-0006); French government grant; Ministry of Science and Higher Education; Polish National Science Centre (2018/29/B/HS6/01934)

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

Abstract

Study Objectives: Sleep loss significantly disrupts cognitive and emotional functioning, yet the neural consequences of different types of sleep deprivation remain unclear.

Methods: In a within-subject resting-state functional magnetic resonance imaging study, we examined how acute total sleep deprivation (TSD) and chronic sleep restriction (CSR) alter intrinsic functional brain organization in 28 healthy adults scanned under three conditions: rested wakefulness (RW), after one night of TSD, and after five nights of CSR.

To quantify network-level disruption, we applied graph-theoretical analyses, including a novel within-subject adaptation of the Hub Disruption Index and Covariate-Constrained Manifold Learning (CCML), an unsupervised embedding technique sensitive to subject-level covariates. Moreover, we assessed subjective sleep quality, sleepiness, and circadian traits.

Results: Both TSD and CSR were associated with a consistent reorganization of graph topology relative to RW. Furthermore, direct comparisons revealed that TSD and CSR affect different brain hubs. Regional changes in degree, closeness, and clustering coefficients were most prominent in subsystems of the default mode network, frontoparietal network, and cerebellum. These differences were also captured in CCML embeddings, supporting the hypothesis that acute and chronic sleep deprivation exert divergent effects on brain connectivity. Findings were robust across graph thresholds, brain atlases, and nodal metrics. Moreover, these results were further supported by subjective measures—sleepiness was associated with reduced network integration in RW, and circadian phenotype emerged as a key determinant of individual sensitivity to sleep loss.

Conclusions: Our results show that TSD and CSR induce divergent alterations in brain functional organization, offering new insights into their neural impact.

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

PatrycjaScislewska/sleep_deprivation_graphs

License: none: the authors keep all their rights

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: af790e173734b687d184ee821b8c2b63bf8947ab, 20 November 2025

Languages: Jupyter (3)

Size: 183 files, 3 scripts

Software Heritage: not archived

Found in: “Data availability”

Holds: README, 3 notebooks

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

Tools: NumPy (3 files), statsmodels (3 files), Matplotlib (2 files), pandas (2 files), SciPy (2 files), NetworkX (1 file), scikit-learn (1 file), seaborn (1 file)

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

28 September 2026: the link answers

5 files

Scripts/Graphs_nodal_global_metrics_HDI_CCML.ipynb, Jupyter, 3,559 lines, 5 matches
Scripts/quality_check_FD_code.ipynb, Jupyter, 772 lines, 2 matches
Scripts/subjective_behavior_state_traits.ipynb, Jupyter, 377 lines, 3 matches
README.md, Text, 36 lines
readme.txt, Text, 1 line

At the source: github.com/PatrycjaScislewska/sleep_deprivation_graphs/

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

Time series, atlases and code used in this study are available in the following GitHub repository: https://github.com/PatrycjaScislewska/sleep_deprivation_graphs/

All preprocessing steps can be fully reproduced using the following GitHub repository: https://github.com/veronicamunoz/rs_graph_processing

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, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 6 authors, 8 keywords, 12 MeSH terms, 6 funders, 80 references.

Cite

This paper

Scislewska, P., Cabrera Vazquez, A., Szatkowska, I., Kontrymowicz-Ogińska, H., Achard, S., & Domagalik, A. (2026). Divergent disruption of brain networks following total and chronic sleep loss: a longitudinal fMRI study. Sleep, 49(5), zsag030. https://doi.org/10.1093/sleep/zsag030

BibTeX

@article{scislewska2026divergent, author = {Scislewska, Patrycja and Cabrera Vazquez, Arturo and Szatkowska, Iwona and Kontrymowicz-Ogińska, Halszka and Achard, Sophie and Domagalik, Aleksandra}, title = {{Divergent disruption of brain networks following total and chronic sleep loss: a longitudinal fMRI study}}, journal = {Sleep}, year = {2026}, month = may, volume = {49}, number = {5}, pages = {zsag030}, publisher = {Oxford University Press}, issn = {0161-8105}, doi = {10.1093/sleep/zsag030}, url = {https://doi.org/10.1093/sleep/zsag030}, pmid = {41631633}, pmcid = {PMC13163182} }

RIS

TY - JOUR AU - Scislewska, Patrycja AU - Cabrera Vazquez, Arturo AU - Szatkowska, Iwona AU - Kontrymowicz-Ogińska, Halszka AU - Achard, Sophie AU - Domagalik, Aleksandra TI - Divergent disruption of brain networks following total and chronic sleep loss: a longitudinal fMRI study T2 - Sleep J2 - Sleep PY - 2026 DA - 2026/05/01 VL - 49 IS - 5 SP - zsag030 SN - 0161-8105 PB - Oxford University Press DO - 10.1093/sleep/zsag030 UR - https://doi.org/10.1093/sleep/zsag030 LA - en ER -

CSL-JSON

{ "id": "10.1093/sleep/zsag030", "type": "article-journal", "title": "Divergent disruption of brain networks following total and chronic sleep loss: a longitudinal fMRI study", "container-title": "Sleep", "author": [ { "family": "Scislewska", "given": "Patrycja" }, { "family": "Cabrera Vazquez", "given": "Arturo" }, { "family": "Szatkowska", "given": "Iwona" }, { "family": "Kontrymowicz-Ogińska", "given": "Halszka" }, { "family": "Achard", "given": "Sophie" }, { "family": "Domagalik", "given": "Aleksandra" } ], "container-title-short": "Sleep", "volume": "49", "issue": "5", "page": "zsag030", "DOI": "10.1093/sleep/zsag030", "PMID": "41631633", "PMCID": "PMC13163182", "ISSN": "0161-8105", "publisher": "Oxford University Press", "URL": "https://doi.org/10.1093/sleep/zsag030", "language": "en", "issued": { "date-parts": [ [ 2026, 5, 1 ] ] } }

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

Similar papers

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

[1] doi:10.1038/s41467-026-75661-x [code]: A neural signature of sleep deprivation in the human brain.
Journal: Nature communications
In common: statsmodels, seaborn, scikit-learn, 4 other tools, 5 references
[2] doi:10.1162/imag.a.1278 [code]: Young and old adult brains experience opposite effects of acute sleep restriction on the functional connectivity network.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: NetworkX, statsmodels, seaborn, 4 other tools, fMRI, 4 references
[3] doi:10.1002/hbm.70522 [code]: Investigating Emotional Reactivity in Experienced Users of Psychedelics: A Cross-Sectional fMRI Study.
Journal: Human brain mapping
In common: fMRI, 4 references, author Aleksandra Domagalik
[4] doi:10.1073/pnas.2531706123 [code]: Metabolism-weighted brain connectome reveals synaptic integration and vulnerability to neurodegeneration.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: NetworkX, statsmodels, seaborn, 5 other tools, 3 references
[5] doi:10.1016/j.isci.2026.116173 [code]: Hierarchical sparse spatiotemporal graph neural network for brain graph classification.
Journal: iScience
In common: NetworkX, seaborn, scikit-learn, 4 other tools, fMRI, 3 references
[6] doi:10.1016/j.patter.2026.101560 [code]: Automating region selection with genetic algorithms for energy landscape analyses of brain dynamics.
Journal: Patterns (New York, N.Y.)
In common: NetworkX, statsmodels, seaborn, 4 other tools, fMRI, 3 references
[7] doi:10.3389/fnins.2026.1873417 [code]: Functional brain network alterations in a rat model of dental malocclusion.
Journal: Frontiers in neuroscience
In common: seaborn, scikit-learn, pandas, 3 other tools, fMRI, systems, 4 references
[8] doi:10.1038/s41467-026-75959-w [code]: Charting higher-order models of brain function beyond pairwise interactions.
Journal: Nature communications
In common: NetworkX, statsmodels, seaborn, 5 other tools, 2 references
[9] doi:10.1016/j.nicl.2026.104012 [code]: Structural-functional multilayer brain network properties and outcome of combined repetitive transcranial magnetic stimulation and psychotherapy for obsessive-compulsive disorder.
Journal: NeuroImage. Clinical
In common: NetworkX, statsmodels, seaborn, 5 other tools, fMRI, 2 references
[10] doi:10.1002/alz.71365 [code]: Benchmarking speech biomarkers of Alzheimer's against cognitive and neural measures.
Journal: Alzheimer's & dementia : the journal of the Alzheimer's Association
In common: statsmodels, seaborn, scikit-learn, 4 other tools, fMRI, 3 references

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.