OSCR

Spatiotemporal asymmetries on brain energy landscape uncover system entrapment related to depression severity.

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. [1] § Results › Asymmetric state switching relates to anhedonia and rumination ↔ tutorial.ipynb, lines 491–514 · score 0.76 · medication status, RRS Brooding, RRS Depression, RRS Reflection, MASQ, QIDS
  2. [2] § Results › Asymmetric state switching relates to anhedonia and rumination ↔ tutorial.ipynb, lines 491–514 · score 0.71 · MASQ GD scores, MASQ AD scores, MASQ AA, QIDS, depression, MDD
  3. [3] § Results › Canonical RSNs described by spontaneous coactivation patterns ↔ tutorial.ipynb, lines 377–404 · score 0.67 · medoid silhouette coefficients, correlation distance, cluster variance, iterations
  4. [4] § Methods › Data acquisition ↔ meica.libs/nibabel/parrec.py, lines 29–67 · score 0.67 · phase encoding, AP, gradient, EPI, repetition, resolution
  5. [5] § Methods › Clustering of fMRI volumes ↔ tutorial.ipynb, lines 377–404 · score 0.66 · medoid silhouette coefficients, pairwise correlation, iteration, variance, clusters
  6. [6] § Methods › Network control theory and dynamics on networks ↔ tutorial.ipynb, lines 241–288 · score 0.60 · state trajectories, control signals, structural connectivity, energies, matrix, transition
  7. [7] § Results › Structural connectivity modulates empirical state transitions ↔ tutorial.ipynb, lines 241–288 · score 0.59 · state trajectories, control signal, Control energy, structural connectomes, transitions
  8. [8] § Methods › Time series extraction and connectome construction ↔ tutorial.ipynb, lines 601–635 · score 0.52 · nearest neighbor interpolation, parcellation, affine, mapped, space, brain
  9. [9] § Methods › Statistical inference ↔ tutorial.ipynb, lines 1099–1160 · score 0.52 · clinical score, transition probability, sex, medication, age
  10. [10] § Results › Brain state dynamics are associated with depression severity ↔ tutorial.ipynb, lines 853–892 · score 0.51 · MASQ AD, clinical scores, medication, age, HC, QIDS

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 authors' code

Jupyter notebook · 1,306 lines · 60 KB · GPL-3.0 · 9 matches

  1. # %%
  2. import numpy as np
  3. import pandas as pd
  5. import statsmodels.api as sm
  6. from statsmodels.stats.multitest import fdrcorrection
  7. from scipy.spatial.distance import squareform, pdist, cdist, cosine
  8. from scipy.stats import zscore, ttest_ind, ttest_rel
  9. import scipy.integrate as spint
  11. if not hasattr(spint, 'simps'):
  12. spint.simps = spint.simpson
  14. from sklearn_extra.cluster import KMedoids
  15. from itertools import groupby
  17. import nibabel as nib
  18. from nilearn.image import resample_img, resample_to_img
  19. from nilearn.plotting import plot_stat_map
  21. from nctpy.energies import get_control_inputs, integrate_u
  22. from nctpy.utils import normalize_state, matrix_normalization
  24. import matplotlib.pyplot as plt
  25. from mpl_toolkits.axes_grid1 import make_axes_locatable
  26. import matplotlib.cm as cm
  27. from matplotlib.patches import RegularPolygon
  28. from matplotlib.transforms import Affine2D
  29. from matplotlib.path import Path
  30. from matplotlib.projections import register_projection
  31. from matplotlib.projections.polar import PolarAxes
  32. from matplotlib.spines import Spine
  34. import seaborn as sns
  36. # %%
  37. def r2(x,y):
  38. return(cdist(np.array([x]), np.array([y]),'correlation')[0][0])
  40. def medoid_silhouette(diss, ind):
  41. medoid_silhouette_values = np.zeros(diss.shape[0])
  42. for ind1 in range(diss.shape[0]):
  43. silh = sorted([diss[ind1, ind2] for ind2 in ind])
  44. num = silh[0]#closest cluster center
  45. if silh[1] == 0:
  46. medoid_silhouette_values[ind1] = 1
  47. else:
  48. medoid_silhouette_values[ind1] = 1-num/silh[1]
  49. return(medoid_silhouette_values)
  51. def silhouette_plot(values, clusters, K, it, ax):
  52. y_lower = 10
  53. for k in range(K):
  54. # Aggregate the silhouette scores for samples belonging to
  55. # cluster i, and sort them
  56. ith_cluster_silhouette_values = values[K,it][clusters[K,it] == k]
  57. ith_cluster_silhouette_values.sort()
  59. size_cluster_i = ith_cluster_silhouette_values.shape[0]
  60. y_upper = y_lower + size_cluster_i
  62. color = cm.nipy_spectral(float(k) / K)
  63. ax.fill_betweenx(np.arange(y_lower, y_upper),0,ith_cluster_silhouette_values,facecolor=color,edgecolor=color,alpha=0.7)
  65. # Label the silhouette plots with their cluster numbers at the middle
  66. ax.text(-0.05, y_lower + 0.5 * size_cluster_i, str(k))
  68. # Compute the new y_lower for next plot
  69. y_lower = y_upper + 10 # 10 for the 0 samples
  70. # The vertical line for average silhouette score of all the values
  71. ax.set_title('Silhouette Plot', fontsize = 16)
  72. ax.set_xlabel('Silhouette Coefficient Values', fontsize = 14)
  73. ax.set_ylabel('Cluster label', fontsize = 14)
  74. ax.axvline(x=values[K,it].mean(), color="black", linestyle="--", label = 'Mean')
  75. ax.set_yticks([]) # Clear the yaxis labels / ticks
  76. ax.set_xticks([-0.1, 0, 0.2, 0.4, 0.6, 0.8, 1])
  77. ax.legend()
  79. def fractional_occs(labels, K):#calculate fractional occupancy for a given subject's state sequence and number of clusters
  80. uniq, cts = np.unique(labels, return_counts=True)
  81. frq = np.zeros(K)
  82. for _, k in enumerate(uniq):
  83. frq[k] = cts[_]/len(labels)#percentage
  84. return(frq)
  86. def dwell_times(labels, K, t_r, thresh = 1):
  87. #calculate dwell times for a given subject's state sequence and number of clusters. thresh variable adjusts how many repetitions counts as dwelling i.e., if thresh ==1, at least 2 consecutive repetitions counts as dwelling
  88. count_dups = [(_,sum(1 for _ in group)) for _, group in groupby(labels)]
  89. dwell_t = [[] for _ in range(K)]
  90. for st, ct in count_dups:
  91. if ct > thresh:
  92. dwell_t[st].append(ct)
  93. dwelltime = []
  94. for dwell in dwell_t:
  95. if dwell:
  96. dwelltime.append(np.mean(dwell) * t_r)
  97. else:
  98. dwelltime.append(0 * t_r)
  99. return(dwelltime)
  102. def get_yeo_networks(states, yeo8, resampled_aparcaseg, size, K):
  103. all_cos = []
  104. for k in range(K):
  105. yeo_array = yeo8.get_fdata()[:,:,:,0]
  106. aparc_array = resampled_aparcaseg.get_fdata()
  108. cos_sim = []
  109. for y in range(1,len(np.unique(yeo_array))):
  110. yeo_networkx = np.where(yeo_array == y)
  111. aparc_network = aparc_array[yeo_networkx]
  113. uniq,cts = np.unique(aparc_network[np.nonzero(aparc_network)], return_counts=True)
  114. roi_counts = {int(roi): count for roi, count in zip(uniq, cts)}
  115. yeo_vector = np.zeros(size)
  116. for i,roi in enumerate(np.unique(aparc_array)[1:]): #ignore 0
  117. if int(roi) in roi_counts:
  118. total_roi = np.sum(aparc_array == int(roi))
  119. yeo_vector[i] = roi_counts[int(roi)]/total_roi if total_roi > 0 else 0
  120. cos_sim.append(float(1-cosine(yeo_vector,states[k])))
  121. all_cos.append(cos_sim)
  122. return(all_cos)
  125. def radar_factory(num_vars, frame='circle'):
  126. """
  127. Create a radar chart with `num_vars` Axes.
  129. This function creates a RadarAxes projection and registers it.
  131. Parameters
  132. ----------
  133. num_vars : int
  134. Number of variables for radar chart.
  135. frame : {'circle', 'polygon'}
  136. Shape of frame surrounding Axes.
  138. """
  139. # calculate evenly-spaced axis angles
  140. theta = np.linspace(0, 2*np.pi, num_vars, endpoint=False)
  142. class RadarTransform(PolarAxes.PolarTransform):
  144. def transform_path_non_affine(self, path):
  145. # Paths with non-unit interpolation steps correspond to gridlines,
  146. # in which case we force interpolation (to defeat PolarTransform's
  147. # autoconversion to circular arcs).
  148. if path._interpolation_steps > 1:
  149. path = path.interpolated(num_vars)
  150. return Path(self.transform(path.vertices), path.codes)
  152. class RadarAxes(PolarAxes):
  154. name = 'radar'
  155. PolarTransform = RadarTransform
  157. def __init__(self, *args, **kwargs):
  158. super().__init__(*args, **kwargs)
  159. # rotate plot such that the first axis is at the top
  160. self.set_theta_zero_location('N')
  162. def fill(self, *args, closed=True, **kwargs):
  163. """Override fill so that line is closed by default"""
  164. return super().fill(closed=closed, *args, **kwargs)
  166. def plot(self, *args, **kwargs):
  167. """Override plot so that line is closed by default"""
  168. lines = super().plot(*args, **kwargs)
  169. for line in lines:
  170. self._close_line(line)
  172. def _close_line(self, line):
  173. x, y = line.get_data()
  174. # FIXME: markers at x[0], y[0] get doubled-up
  175. if x[0] != x[-1]:
  176. x = np.append(x, x[0])
  177. y = np.append(y, y[0])
  178. line.set_data(x, y)
  180. def set_varlabels(self, labels, fs):
  181. self.set_thetagrids(np.degrees(theta), labels, fontsize = fs)
  183. def _gen_axes_patch(self):
  184. # The Axes patch must be centered at (0.5, 0.5) and of radius 0.5
  185. # in axes coordinates.
  186. if frame == 'circle':
  187. return Circle((0.5, 0.5), 0.5)
  188. elif frame == 'polygon':
  189. return RegularPolygon((0.5, 0.5), num_vars, radius=0.5, edgecolor="k")
  190. else:
  191. raise ValueError("Unknown value for 'frame': %s" % frame)
  193. def _gen_axes_spines(self):
  194. if frame == 'circle':
  195. return super()._gen_axes_spines()
  196. elif frame == 'polygon':
  197. # spine_type must be 'left'/'right'/'top'/'bottom'/'circle'.
  198. spine = Spine(axes=self,
  199. spine_type='circle',
  200. path=Path.unit_regular_polygon(num_vars))
  201. # unit_regular_polygon gives a polygon of radius 1 centered at
  202. # (0, 0) but we want a polygon of radius 0.5 centered at (0.5,
  203. # 0.5) in axes coordinates.
  206. spine.set_transform(Affine2D().scale(0.5).translate(.5, .5)+ self.transAxes)
  207. return {'polar': spine}
  208. else:
  209. raise ValueError("Unknown value for 'frame': %s" % frame)
  211. register_projection(RadarAxes)
  212. return theta
  214. def organize_stats(tempdf):
  215. # Remove commas from string values in the DataFrame (e.g., "1,234" → "1234")
  216. # This helps ensure numeric columns are clean for later processing.
  217. df = tempdf.map(lambda x: x.replace(',', '') if isinstance(x, str) else x)
  219. # Split the first column's values into multiple columns based on one or more spaces
  220. # (regex '\s+' matches one or more whitespace characters)
  221. df_split = df[df.keys()[0]].str.split(r'\s+', expand=True)
  223. # Assign descriptive names to the new columns
  224. # The last column is unnamed because it might contain an extra blank split
  225. df_split.columns = ["ROI", "Name", "Type", "Volume-mm3", ""]
  227. # Drop the unnecessary empty column that resulted from splitting
  228. df_cleaned = df_split.drop(columns=[''])
  230. # Return the cleaned DataFrame
  231. return df_cleaned
  234. # Function to generate a single random adjacency matrix
  235. def random_adjacency(n):
  236. mat = np.random.rand(n, n) # random values [0,1)
  237. mat = (mat + mat.T)/2 # make symmetric
  238. np.fill_diagonal(mat, 0) # set diagonal to 0
  239. return mat
  241. def NCT_multi(healthyID, mddID, DTI_matrix, states, K, size, time_h = 1, rho = 1, system = 'continuous'):
  243. s_trajectories = {idx:np.zeros((K,K,1001,size)) for idx in healthyID+mddID}
  244. c_signals= {idx:np.zeros((K,K,1001,size)) for idx in healthyID+mddID}
  245. n_energies = {idx:np.zeros((K,K,size)) for idx in healthyID+mddID}
  246. t_energies = {idx:np.zeros((K,K)) for idx in healthyID+mddID}
  248. s_trajectories_pers = {idx:np.zeros((K,K,1001,size)) for idx in healthyID+mddID}
  249. c_signals_pers = {idx:np.zeros((K,K,1001,size)) for idx in healthyID+mddID}
  250. n_energies_pers = {idx:np.zeros((K,K,size)) for idx in healthyID+mddID}
  251. t_energies_pers = {idx:np.zeros((K,K)) for idx in healthyID+mddID}
  253. control_s, trajectory_c = np.eye(size), np.eye(size)
  254. # normalize structural connectivity
  255. for idx in healthyID+mddID:
  256. norm_adjs_struct = matrix_normalization(A=DTI_matrix[idx], system=system, c=1)
  257. for k1 in range(K):
  258. for k2 in range(K):
  259. if k1 != k2:
  260. #get the state trajectory, x(t), and the control signals, u(t)
  261. s_trajectories[idx][k2,k1], c_signals[idx][k2,k1], numerical_errors = get_control_inputs(A_norm = norm_adjs_struct, T = time_h, B = control_s, x0 = states[k1], xf = states[k2], system = system, rho = rho, S = trajectory_c)
  262. # print errors
  263. thr = 1e-8
  264. if (numerical_errors[0] >= thr) or (numerical_errors[1] >= thr):
  265. # the first numerical error corresponds to the inversion error # the second numerical error corresponds to the reconstruction error
  266. print("Subject: %s, Transition: %dto%d"%(idx,k1,k2),
  267. "Inversion error = {:.2E} (<{:.2E}={:})".format(numerical_errors[0], thr, numerical_errors[0] < thr),
  268. "Reconstruction error = {:.2E} (<{:.2E}={:})".format(numerical_errors[1], thr, numerical_errors[1] < thr))
  270. # integrate control signals to get control energy
  271. n_energies[idx][k2,k1] = integrate_u(c_signals[idx][k2,k1])
  272. t_energies[idx][k2,k1] = np.sum(n_energies[idx][k2,k1]) ##from state k1 to k2
  273. else:#get persistence energies in a seperate array
  274. #get the state trajectory, x(t), and the control signals, u(t)
  275. s_trajectories_pers[idx][k2,k1], c_signals_pers[idx][k2,k1], numerical_errors = get_control_inputs(A_norm = norm_adjs_struct, T = time_h, B = control_s, x0 = states[k1], xf = states[k2], system = system, rho = rho, S = trajectory_c)
  276. # print errors
  277. thr = 1e-8
  278. if (numerical_errors[0] >= thr) or (numerical_errors[1] >= thr):
  279. # the first numerical error corresponds to the inversion error # the second numerical error corresponds to the reconstruction error
  280. print("Subject: %s, Transition: %dto%d"%(idx,k1,k2),
  281. "Inversion error = {:.2E} (<{:.2E}={:})".format(numerical_errors[0], thr, numerical_errors[0] < thr),
  282. "Reconstruction error = {:.2E} (<{:.2E}={:})".format(numerical_errors[1], thr, numerical_errors[1] < thr))
  284. # integrate control signals to get control energy
  285. n_energies_pers[idx][k2,k1] = integrate_u(c_signals_pers[idx][k2,k1])
  286. t_energies_pers[idx][k2,k1] = np.sum(n_energies_pers[idx][k2,k1]) ##from state k1 to k2
  288. return(t_energies, n_energies, s_trajectories, c_signals, t_energies_pers, n_energies_pers, s_trajectories_pers, c_signals_pers)
  290. # %%
  291. working_path = '/Path/to/Brain_states/'
  293. # Adjust your Parameters
  294. size = 85 # number of regions
  295. duration = 200 # time points per subject
  296. n_subjects = 20 # total subjects
  297. n_healthy = n_subjects // 2
  298. n_mdd = n_subjects // 2
  300. # Simulate data as dicts: each subject -> size x duration matrix
  301. np.random.seed(42) # for reproducibility
  302. healthyFMRI = {str(500+i): np.random.randn(size, duration) for i in range(n_healthy)}
  303. mddFMRI = {str(600+i): np.random.randn(size, duration) for i in range(n_mdd)}
  305. # %%
  306. # Get IDs
  307. healthyID = list(healthyFMRI.keys())
  308. mddID = list(mddFMRI.keys())
  310. # Preallocate arrays
  311. X_hcz1 = np.zeros((size, sum([healthyFMRI[idx].shape[1] for idx in healthyID])))
  312. X_mddz1 = np.zeros((size, sum([mddFMRI[idx].shape[1] for idx in mddID])))
  314. # Fill in healthy
  315. dur = 0
  316. for idx in healthyID:
  317. mat = zscore(healthyFMRI[idx], axis=1, nan_policy='raise')
  318. X_hcz1[:, dur:dur + mat.shape[1]] = mat
  319. dur += mat.shape[1]
  321. # Fill in mdd
  322. dur = 0
  323. for idx in mddID:
  324. mat = zscore(mddFMRI[idx], axis=1, nan_policy='raise')
  325. X_mddz1[:, dur:dur + mat.shape[1]] = mat
  326. dur += mat.shape[1]
  328. # Concatenate
  329. X_prez1 = np.concatenate((X_hcz1, X_mddz1), axis=1)
  330. X_prez1[np.isnan(X_prez1)] = 0
  332. # %%
  333. #Display point cloud
  334. fig, ax = plt.subplots(1,1,figsize = (21,4))
  335. im = ax.imshow(X_prez1, aspect = 'auto', interpolation = 'none', cmap = 'coolwarm')
  336. divider = make_axes_locatable(ax)
  337. cax = divider.append_axes('right', size='5%', pad=0.05)
  338. fig.colorbar(im, cax=cax)
  339. plt.tight_layout()
  341. # %%
  342. # Define range of cluster numbers to test (K = 2 to 10)
  343. Ks = np.arange(2, 10)
  345. # Number of iterations (repeated clustering runs for stability analysis)
  346. niter = 10
  348. clusters, indices, inertias = {}, {}, {}
  350. # Compute the pairwise dissimilarity matrix between all time points
  351. # - X_prez1.T: shape (n_timepoints, n_regions)
  352. # - 'correlation' metric → distance = 1 - correlation coefficient
  353. # - squareform: convert condensed distance vector into a square matrix
  354. diss = squareform(pdist(X_prez1.T, metric='correlation'))
  356. # Run clustering multiple times for each K
  357. for it in range(niter):
  358. for K in Ks:
  359. # Initialize and fit K-Medoids clustering with:
  360. # - precomputed distance matrix
  361. # - 'k-medoids++' initialization
  362. # - 'alternate' approximation method which is fast or 'pam' method which is more robust but much slower
  363. km = KMedoids(init="k-medoids++",n_clusters=K,metric='precomputed',method='alternate').fit(diss)
  365. # Store the results for this iteration and K
  366. clusters[K, it] = km.labels_
  367. indices[K, it] = km.medoid_indices_
  368. inertias[K, it] = km.inertia_
  370. # Dictionary to store the actual medoid time series
  371. medoids = {}
  372. for it in range(niter):
  373. for _, K in enumerate(Ks):
  374. # Extract medoid time series from X_prez1.T using saved medoid indices
  375. medoids[K, it] = X_prez1.T[indices[K, it]]
  377. # %%
  378. bt_vr, wt_vr = {}, {} # Initialize dictionaries to store between- and within-cluster variances
  380. for it in range(niter): # Loop over iterations
  381. between_var, within_var = [], [] # Lists to store variances for each K
  382. for K in Ks: # Loop over different numbers of clusters
  383. btw_var, wtn_var = [], [] # Lists to store variance per cluster
  384. for k1 in range(K): # Loop over each cluster
  385. # Calculate BETWEEN-cluster variance for cluster k1
  386. var = [r2(medoids[K, it][k2], medoids[K, it][k1]) for k2 in range(K)] # Pairwise r2 with all other medoids
  387. btw_var.append(np.array(var).mean()) # Average r2 across other clusters
  389. # Calculate WITHIN-cluster variance for cluster k1
  390. var = [r2(X_prez1.T[tr_label], medoids[K, it][k1]) for tr_label in np.where(clusters[K, it] == k1)[0]] # r2 with members
  391. wtn_var.append(np.array(var).mean()) # Average r2 within the cluster
  393. between_var.append(np.array(btw_var).mean()) # Average BETWEEN variance across clusters
  394. within_var.append(np.array(wtn_var).mean()) # Average WITHIN variance across clusters
  396. bt_vr[it] = between_var # Store results for this iteration
  397. wt_vr[it] = within_var
  399. # Calculate medoid silhouette coefficients
  400. medoids_silhouette1_values = {} # Initialize dictionary
  401. diss = squareform(pdist(X_prez1.T, metric='correlation')) # Compute pairwise correlation distance matrix
  402. for it in range(niter): # Loop over iterations
  403. for K in Ks: # Loop over cluster sizes
  404. medoids_silhouette1_values[K, it] = medoid_silhouette(diss, indices[K, it]) # Compute silhouette for medoids
  406. # %%
  407. explained_var, var_gain = {}, {}
  409. fig,ax = plt.subplots(1,6, figsize = (36,4))
  410. for it in range(niter):
  411. ax[0].plot(Ks, wt_vr[it], marker = '*')
  412. ax[1].plot(Ks, bt_vr[it], marker = '*')
  414. explained_var[it] = bt_vr[it]/(np.array(bt_vr[it])+np.array(wt_vr[it]))
  415. ax[2].plot(Ks, explained_var[it], marker = '>')
  417. var_gain[it] = np.diff(bt_vr[it]/(np.array(bt_vr[it])+np.array(wt_vr[it])))
  418. ax[3].plot(Ks[1:], var_gain[it], marker = 'x')
  420. ax[4].plot(Ks,[medoids_silhouette1_values[K,it].mean() for K in Ks])
  422. ax[0].plot(Ks, np.array(list(wt_vr.values())).mean(axis = 0), c = 'black', lw = 2, ls = '--', label = 'Mean')
  423. ax[1].plot(Ks, np.array(list(bt_vr.values())).mean(axis = 0), c = 'black', lw = 2, ls = '--', label = 'Mean')
  424. ax[2].plot(Ks, np.array(list(explained_var.values())).mean(axis = 0), c = 'black', lw = 2, ls = '--', label = 'Mean')
  425. ax[3].plot(Ks[1:],np.array(list(var_gain.values())).mean(axis = 0), c = 'black', lw = 2, ls = '--', label = 'Mean')
  426. ax[4].plot(Ks,np.array([[medoids_silhouette1_values[K,it].mean() for K in Ks] for it in range(niter)]).mean(axis = 0), c = 'black', lw = 2, ls = '--', label = 'Mean')
  428. K=4
  429. it = np.argmax([inertias[K,it] for it in range(niter)])
  430. silhouette_plot(medoids_silhouette1_values, clusters, K, it, ax[5])
  432. for i in range(5):
  433. ax[i].set_xticks(Ks)
  434. ax[i].legend()
  435. ax[i].set_xlabel('# of Clusters (K)', fontsize = 12)
  437. ax[0].set_title('Within Cluster Variance', fontsize = 16)
  438. ax[1].set_title('Between Cluster Variance', fontsize = 16)
  439. ax[2].set_title('Explained Variance', fontsize = 16)
  440. ax[3].set_title('Variance Gain', fontsize = 16) # variance gain from increasing from k-1 to k
  441. ax[4].set_title('Mean Silhouette Scores', fontsize = 16)
  443. ax[2].set_ylabel(r'$R^{2}$', fontsize = 12)
  444. ax[3].set_ylabel(r'$R^{2}$ Gain', fontsize = 12)
  445. ax[4].set_ylabel('Silhouette Coefficient', fontsize = 14)
  447. plt.tight_layout()
  449. # %%
  450. spatial_corr_matrices = {}
  451. fig,ax = plt.subplots(1,1,figsize = (12,4))
  452. for _K, K in enumerate(Ks):
  453. it = np.argmax([inertias[K,it] for it in range(niter)])
  454. spatial_corr_matrices[K] = 1-squareform(pdist(medoids[K,it],metric='correlation'))
  455. upper = np.triu(spatial_corr_matrices[K],1)
  456. upper[np.nonzero(upper)]
  457. ax.scatter([K]*len(upper[np.nonzero(upper)]),upper[np.nonzero(upper)], s = 25, marker = 'x')
  458. ax.set_xlabel('# of Clusters (K)', fontsize = 14)
  459. ax.set_ylabel('Correlation', fontsize = 14)
  460. ax.set_title('Spatial correlation between each pair of states', fontsize = 16)
  461. ax.errorbar(Ks, [(np.triu(spatial_corr_matrices[K],1)[np.nonzero(np.triu(spatial_corr_matrices[K],1))]).mean() for K in Ks], [(np.triu(spatial_corr_matrices[K],1)[np.nonzero(np.triu(spatial_corr_matrices[K],1))]).std() for K in Ks], c = 'black', lw = 1.5,ls = 'dashed', label = 'Mean')
  463. ax.set_xticks(Ks)
  464. ax.legend()
  466. # %%
  467. state_fraq_occ, dwellTIME = {}, {} # Initialize dictionaries to store fractional occupancy and dwell times
  469. tr = 2 # Repetition time (TR) in seconds
  471. for K in Ks: # Loop over number of clusters
  472. it = np.argmax([inertias[K, it] for it in range(niter)]) # Select iteration with maximum inertia
  473. dur = 0 # Initialize cumulative duration index
  474. for i, idx in enumerate(healthyID + mddID): # Loop over all subjects
  475. if idx in healthyID:
  476. size, duration = healthyFMRI[idx].shape # Get subject data shape
  477. else:
  478. size, duration = mddFMRI[idx].shape
  480. subj_labels = clusters[K, it][dur:dur + duration] # Extract cluster labels for this subject
  481. dur = dur + duration # Update cumulative duration
  483. # Calculate fractional occupancy and dwell times for each state
  484. frq = fractional_occs(subj_labels, K) # Fraction of time spent in each state
  485. dwell_time = dwell_times(subj_labels, K, tr) # Average dwell time in seconds
  487. state_fraq_occ[idx, K] = frq * 100 # Store fractional occupancy as percentage
  488. dwellTIME[idx, K] = dwell_time # Store dwell time in seconds
  491. # %%
  492. # Import clinical scores CSV (these ara randomly generated)
  493. clinical = pd.read_csv(working_path + 'Clinical_Scores.csv')
  495. # List of clinical measures of interest
  496. clinical_measures = ['QIDS', 'MASQ_AA', 'MASQ_AD', 'MASQ_GD', 'RRS_DR', 'RRS_BR', 'RRS_RF']
  498. # Create dictionaries mapping subject ID → demographic/medication info
  499. SEX = {idx: 0 if clinical[clinical.study_ID == int(idx)].SUB_gender.iloc[0] == "F" else 1 for idx in healthyID + mddID} # 0=F, 1=M
  501. AGE = {idx: clinical[clinical.study_ID == int(idx)]['SUB_Age '].iloc[0] for idx in healthyID + mddID} # Age in years
  503. MEDS = {idx: 0 if clinical[clinical.study_ID == int(idx)].Medication.iloc[0] == 'no' else 1 for idx in healthyID + mddID} # Medication status (0=no, 1=yes)
  505. # Create dictionaries mapping subject ID → clinical scores
  506. QIDS = {idx: clinical[clinical.study_ID == int(idx)].QIDS_score.iloc[0] for idx in healthyID + mddID}
  508. MASQ_aa = {idx: clinical[clinical.study_ID == int(idx)].MASQ_aa_score.iloc[0] for idx in healthyID + mddID}
  509. MASQ_ad = {idx: clinical[clinical.study_ID == int(idx)].MASQ_ad_score.iloc[0] for idx in healthyID + mddID}
  510. MASQ_gd = {idx: clinical[clinical.study_ID == int(idx)].MASQ_gd_score.iloc[0] for idx in healthyID + mddID}
  512. RRS_dr = {idx: clinical[clinical.study_ID == int(idx)].RRS_Depression_Related.iloc[0] for idx in healthyID + mddID}
  513. RRS_b = {idx: clinical[clinical.study_ID == int(idx)].RRS_Brooding.iloc[0] for idx in healthyID + mddID}
  514. RRS_r = {idx: clinical[clinical.study_ID == int(idx)].RRS_Reflection.iloc[0] for idx in healthyID + mddID}
  516. # %%
  517. # Initialize lists to collect subject-level, state-level, and clinical data
  518. totst_, KS_, ids_, type_, sex_, age_, meds_ = [], [], [], [], [], [], []
  519. dwellt_, fraqoc_ = [], []
  520. qids_, masqaa_, masqad_, masqgd_, rrsdr_, rrsb_, rrsr_ = [], [], [], [], [], [], []
  522. # Loop over all subjects (healthy + MDD)
  523. for i, idx in enumerate(healthyID + mddID):
  525. # Assign group label
  526. typ = 'HC' if idx in healthyID else 'MDD'
  528. # Loop over all tested K values and states within each K
  529. for K in Ks:
  530. for k in range(K):
  531. # Basic subject and state info
  532. totst_.append(K)
  533. KS_.append(k)
  534. ids_.append(idx)
  535. type_.append(typ)
  536. sex_.append(SEX[idx])
  537. age_.append(AGE[idx])
  538. meds_.append(MEDS[idx])
  540. # State-level metrics
  541. dwellt_.append(dwellTIME[idx, K][k])
  542. fraqoc_.append(state_fraq_occ[idx, K][k])
  544. # Clinical scores
  545. qids_.append(QIDS[idx])
  546. masqaa_.append(MASQ_aa[idx])
  547. masqad_.append(MASQ_ad[idx])
  548. masqgd_.append(MASQ_gd[idx])
  549. rrsdr_.append(RRS_dr[idx])
  550. rrsb_.append(RRS_b[idx])
  551. rrsr_.append(RRS_r[idx])
  553. # Create DataFrame with all collected information
  554. df = pd.DataFrame({'SubjectID': ids_, 'Morbidity': type_, 'Sex': sex_, 'Age': age_, 'Medication': meds_,
  555. 'totalstates': totst_, 'StateID': KS_,'DWELLTIME': dwellt_, 'FRAQOCC': fraqoc_,
  556. 'QIDS': qids_, 'MASQ_AA': masqaa_, 'MASQ_AD': masqad_, 'MASQ_GD': masqgd_,'RRS_DR': rrsdr_, 'RRS_BR': rrsb_, 'RRS_RF': rrsr_})
  558. # %%
  559. K = 4
  560. it = np.argmax([inertias[K,it] for it in range(niter)])
  561. my_pal = {"HC": "darkgrey", "MDD": "deepskyblue"}
  563. fig, ax = plt.subplots(1,1, figsize = (10,4))
  564. sns.violinplot(data = df[df.totalstates == K], x = 'StateID', y = 'DWELLTIME', hue = 'Morbidity', ax = ax, palette=my_pal)
  565. ax.set_title('Dwell Time', fontsize = 16)
  566. ax.set_ylabel('Seconds')
  567. dwell_p = ['t:%.3f, p:%.3f'%(ttest_ind(df[(df.Morbidity == 'HC')&(df.StateID == k)&(df.totalstates == K)].DWELLTIME.to_numpy(),df[(df.Morbidity == 'MDD')&(df.StateID == k)&(df.totalstates == K)].DWELLTIME.to_numpy())) for k in range(K)]
  568. print('K = %d'%K)
  569. print('Dwell TIME:', dwell_p, f'\n'
  570. '-----------------------------------')
  572. # %%
  573. K = 4
  574. it = np.argmax([inertias[K,it] for it in range(niter)])
  576. my_pal = {"HC": "darkgrey", "MDD": "deepskyblue"}
  578. fig, ax = plt.subplots(1,1, figsize = (10,4))
  579. sns.violinplot(data = df[df.totalstates == K], x = 'StateID', y = 'FRAQOCC', hue = 'Morbidity', ax = ax, palette=my_pal)
  580. ax.set_title('Fractional Occupancy', fontsize = 16)
  581. ax.set_ylabel('Percentage')
  583. frac_p = ['t:%.3f, p:%.3f'%(ttest_ind(df[(df.Morbidity == 'HC')&(df.StateID == k)&(df.totalstates == K)].FRAQOCC.to_numpy(),df[(df.Morbidity == 'MDD')&(df.StateID == k)&(df.totalstates == K)].FRAQOCC.to_numpy(), equal_var=False)) for k in range(K)]
  584. print('K = %d'%K)
  585. print('Frac Occupancy:', frac_p, f'\n'
  586. '-----------------------------------')
  588. # %%
  589. #ADJUST THESE PATHS, see the github repo for details
  590. brain = nib.load(working_path + 'BIDS_dir/derivatives/cmp-v3.1.0/sub-001/ses-1/anat/sub-001_ses-1_desc-cmp_T1w.nii.gz') #T1w anatomical image
  591. aparcaseg = nib.load(working_path + 'BIDS_dir/derivatives/cmp-v3.1.0/sub-001/ses-1/anat/sub-001_ses-1_atlas-L2018_res-scale1_dseg.nii.gz') # Whole-brain parcellation
  592. yeo7 = nib.load(working_path + 'BIDS_dir/code/Yeo_JNeurophysiol11_MNI152/Yeo2011_7Networks_MNI152_FreeSurferConformed1mm.nii.gz')
  593. yeo7rois = pd.read_csv(working_path + 'BIDS_dir/code/Yeo_JNeurophysiol11_MNI152/7NetworksOrderedNames.csv')
  595. df_stats = pd.read_csv(working_path+ 'BIDS_dir/derivatives/cmp-v3.1.0/sub-001/ses-1/anat/sub-001_ses-1_atlas-L2018_res-scale1_stats.tsv', delimiter= '\t')
  596. df_info = organize_stats(df_stats)
  598. # %%
  599. df_info
  601. # %%
  602. K = 4
  603. it = np.argmax([inertias[K,it] for it in range(niter)])
  605. fig,ax = plt.subplots(1,3,figsize = (16,4))
  606. # Resample aparcaseg (Lausanne2018 parcellation) to match the voxel size and shape of the 'brain' image. Use nearest-neighbor interpolation since these are label maps, not continuous data
  607. resampled_aparcaseg = resample_img(aparcaseg,target_affine=brain.affine,target_shape=brain.shape,interpolation='nearest',force_resample=True,copy_header=True)
  608. # Plot the resampled aparcaseg map on top of the brain anatomy
  609. plot_stat_map(resampled_aparcaseg, threshold=0.01, bg_img=brain, axes=ax[0])
  610. ax[0].set_title('Native brain parcellation (Lausanne2018, N = 85)', fontsize=10)
  612. # Resample the Yeo7 parcellation to match the 'brain' space
  613. resampled_yeo7 = resample_img(yeo7,target_affine=brain.affine,target_shape=brain.shape,interpolation='nearest',force_resample=True,copy_header=True)
  614. # Plot the resampled Yeo7 map with a categorical colormap
  615. plot_stat_map(resampled_yeo7, cmap='Paired', axes=ax[1], bg_img=brain, transparency=0.4)
  616. ax[1].set_title('Native brain parcellation (Yeo7, N = 7)', fontsize=10)
  618. # List of label IDs corresponding to subcortical regions in aparcaseg
  619. subcortical_inds = df_info[(df_info.Type == 'subcortical')].ROI.to_numpy()
  620. # Find voxel indices for each subcortical region in the resampled aparcaseg
  621. subcortical_args = [np.where(resampled_aparcaseg.get_fdata() == int(si)) for si in subcortical_inds]
  623. # Copy the Yeo7 data array so we can modify it
  624. yeo_array = resampled_yeo7.get_fdata().copy()
  625. # Replace all subcortical voxels with a new label value (8), creating Yeo8
  626. for i, si in enumerate(subcortical_inds):
  627. yeo_array[subcortical_args[i]] = 8 # Assign index 8 to all subcortical voxels
  628. # Convert modified array into a NIfTI image and resample to the brain space
  629. yeo8 = resample_to_img(nib.Nifti1Image(yeo_array, affine=resampled_yeo7.affine),brain,force_resample=True,copy_header=True)
  630. # Add the new "SUB" (subcortical) label to the Yeo ROIs dataframe
  631. yeo7rois.loc[len(yeo7rois)] = [8, 'SUB']
  633. # Plot the new Yeo8 parcellation (Yeo7 + subcortical regions)
  634. plot_stat_map(yeo8, cmap='Paired', axes=ax[2], bg_img=brain, transparency=0.4)
  635. ax[2].set_title('Native brain parcellation (Yeo7+Subcortical, N = 8)', fontsize=10)
  637. # %%
  638. K = 4
  639. it = np.argmax([inertias[K,it] for it in range(niter)])
  641. fig,ax = plt.subplots(1,K,figsize = (20,4))
  642. # Loop through each state (cluster) from 0 to K-1
  643. for k in range(K):
  644. # Get the aparcaseg (parcellation) data array
  645. aparcaseg_data = resampled_aparcaseg.get_fdata()
  646. # Initialize a mask for the current state with all zeros
  647. state_mask = np.zeros(aparcaseg_data.shape)
  648. # Loop through each ROI in aparcaseg (skip the first unique value, usually background = 0)
  649. for i, roi in enumerate(np.unique(aparcaseg_data)[1:]):
  650. # Assign the medoid value for this ROI in the current state k
  651. # medoids[K, it][k][i] contains the feature/activation value for ROI i in state k
  652. state_mask[np.where(aparcaseg_data == roi)] = medoids[K, it][k][i]
  653. # Convert the state mask array into a NIfTI image, preserving spatial info
  654. img = nib.Nifti1Image(state_mask, affine=resampled_aparcaseg.affine)
  655. # Plot the state map on top of the brain anatomy
  656. plot_stat_map(img, threshold=0.01, bg_img=brain, axes=ax[k])
  657. ax[k].set_title('State %s' % k, fontsize=12)
  659. # %%
  660. K = 4
  661. it = np.argmax([inertias[K, it] for it in range(niter)])
  663. # Positive amplitude states: keep only positive z-scores, set all others to 0
  664. medoid1 = medoids[K, it].copy()
  665. medoid1[medoid1 <= 0] = 0
  666. pos_amplitude = medoid1.copy()
  668. # Negative amplitude states: keep only negative z-scores, set all others to 0. Then take the absolute value so magnitudes are positive for plotting
  669. medoid2 = medoids[K, it].copy()
  670. medoid2[medoid2 > 0] = 0
  671. neg_amplitude = abs(medoid2.copy())
  673. # Dictionary to store results for each state
  674. all_data8 = {}
  676. # Map positive components to Yeo8 networks
  677. yeo8_pos = get_yeo_networks(pos_amplitude, yeo8, resampled_aparcaseg, size, K)
  679. # Map negative components to Yeo8 networks
  680. yeo8_neg = get_yeo_networks(neg_amplitude, yeo8, resampled_aparcaseg, size, K)
  682. for k in range(K):
  683. all_data8[k] = [yeo8_pos[k], yeo8_neg[k]]
  685. # %%
  686. theta = radar_factory(8, frame = 'polygon')
  688. fig, ax = plt.subplots(figsize=(60, 10), nrows=1, ncols=K,subplot_kw=dict(projection='radar'))
  689. fig.subplots_adjust(wspace=0.25, hspace=0.20, top=0.85, bottom=0.05)
  691. colors = ['b', 'r']
  692. # Plot the four cases from the example data on separate Axes
  693. for k in range(K):
  694. ax[k].set_rgrids([0.2, 0.4, 0.6, 0.8,1], fontsize = 30)
  695. ax[k].set_ylim(0, 1)
  696. for d, color in zip(all_data8[k], colors):
  697. ax[k].plot(theta, d, color='black')
  698. ax[k].fill(theta, d, facecolor=color, alpha=0.25, label='_nolegend_')
  699. ax[k].set_varlabels(yeo7rois[' Network Name'].to_list(), 30)
  700. ax[k].set_title('State %d'%k, weight='bold', fontsize = 30, position=(0.5, 1.1), horizontalalignment='center', verticalalignment='center')
  702. # add legend relative to top-left plot
  703. labels = ('High amplitude', 'Low amplitude')
  704. legend = ax[0].legend(labels, loc=(0.9, .95),labelspacing=0.6, fontsize=40)
  706. # %%
  707. K = 4
  708. it = np.argmax([inertias[K, it] for it in range(niter)])
  710. # Dictionaries for subject-level transition & persistence probabilities
  711. transition_probs, persistence_probs = {}, {}
  712. transition_hc, transition_mdd = np.zeros((len(healthyID), K, K)), np.zeros((len(mddID), K, K))
  713. persistence_hc, persistence_mdd = np.zeros((len(healthyID), K)), np.zeros((len(mddID), K))
  715. # Index to track the time offset when slicing from the concatenated cluster sequence
  716. dur = 0
  717. for i, idx in enumerate(healthyID + mddID):
  718. # Get the fMRI data dimensions for the current subject
  719. if idx in healthyID:
  720. size, duration = healthyFMRI[idx].shape
  721. else:
  722. size, duration = mddFMRI[idx].shape
  723. trnstion_prb = np.zeros((K, K)) # Transition probability matrix
  724. prsstnce_prb = np.zeros((K)) # Persistence counts per state
  726. # Extract the sequence of visited states for this subject. `clusters[K, it]` contains the state assignment for each timepoint in all subjects. We slice from dur to dur+duration to get only this subject's sequence
  727. # new_seq: state sequence without consecutive duplicates (unique transitions)
  728. # new_seq_dup: (state, length_of_consecutive_run) for each run
  729. new_seq = [_ for _, group in groupby(clusters[K, it][dur:dur + duration])]
  730. new_seq_dup = [(_, sum(1 for _ in group)) for _, group in groupby(clusters[K, it][dur:dur + duration])]
  732. # Count how many times each state appears in the deduplicated sequence
  733. uniq, cts = np.unique(new_seq, return_counts=True)
  734. # Update time offset so the next subject's sequence starts correctly
  735. dur += duration
  736. for j, (currst, ct) in enumerate(new_seq_dup):
  737. try:
  738. # Next state after the current run
  739. nextst = new_seq[j + 1]
  741. # Increase transition probability from currst → nextst.
  742. trnstion_prb[nextst, currst] += (cts / cts.sum())[currst]
  744. # Increase persistence probability for currst. (ct - 1) = number of consecutive timepoints staying in the same state
  745. prsstnce_prb[currst] += (ct - 1)
  746. except:
  747. # Last run in the sequence — no next state to transition to
  748. pass
  750. # Normalize transition and persistence probabilities
  751. transition_probs[idx] = trnstion_prb / trnstion_prb.sum()
  752. persistence_probs[idx] = prsstnce_prb / prsstnce_prb.sum()
  754. if idx in healthyID:
  755. transition_hc[i] = trnstion_prb / trnstion_prb.sum()
  756. persistence_hc[i] = prsstnce_prb / prsstnce_prb.sum()
  757. else:
  758. transition_mdd[int(i - len(healthyID))] = trnstion_prb / trnstion_prb.sum()
  759. persistence_mdd[int(i - len(healthyID))] = prsstnce_prb / prsstnce_prb.sum()
  761. # %%
  762. fig,ax = plt.subplots(1,2,figsize = (12,5))
  763. im = ax[0].imshow(transition_hc.mean(axis = 0)-transition_mdd.mean(axis = 0), aspect = 'auto', interpolation = 'none', cmap = 'coolwarm', vmin = -0.02, vmax = 0.02)
  764. divider = make_axes_locatable(ax[0])
  765. cax = divider.append_axes('right', size='5%', pad=0.05)
  766. fig.colorbar(im, cax=cax)
  767. ax[0].set_title('HC-MDD (Transition Prob.)', fontsize = 18)
  768. ax[0].set_xticks([k for k in range(K)])
  769. ax[0].set_yticks([k for k in range(K)])
  770. ax[0].set_xlabel('Current State', fontsize = 15)
  771. ax[0].set_ylabel('Next State', fontsize = 15)
  773. im = ax[1].imshow(np.diag(persistence_hc.mean(axis = 0)-persistence_mdd.mean(axis = 0)), aspect = 'auto', interpolation = 'none', cmap = 'coolwarm', vmin = -0.03, vmax = 0.03)
  774. divider = make_axes_locatable(ax[1])
  775. cax = divider.append_axes('right', size='5%', pad=0.05)
  776. fig.colorbar(im, cax=cax)
  777. ax[1].set_title('HC-MDD (Persistence Prob.)', fontsize = 18)
  778. ax[1].set_xticks([k for k in range(K)])
  779. ax[1].set_yticks([k for k in range(K)])
  780. ax[1].set_xlabel('Current State', fontsize = 15)
  781. ax[1].set_ylabel('Next State', fontsize = 15)
  783. plt.tight_layout()
  785. # %%
  786. fig,ax = plt.subplots(1,2,figsize = (12,2))
  787. im = ax[0].imshow((transition_hc.mean(axis = 0).sum(axis=0)-transition_mdd.mean(axis = 0).sum(axis=0)).reshape(1,4)/3, aspect = 'auto', interpolation = 'none', cmap = 'coolwarm')
  788. divider = make_axes_locatable(ax[0])
  789. cax = divider.append_axes('right', size='5%', pad=0.05)
  790. fig.colorbar(im, cax=cax)
  791. ax[0].set_title('HC-MDD (Exit Prob.)', fontsize = 18)
  792. ax[0].set_xticks([k for k in range(K)])
  794. im = ax[1].imshow((transition_hc.mean(axis = 0).sum(axis=1)-transition_mdd.mean(axis = 0).sum(axis=1)).reshape(1,4)/3, aspect = 'auto', interpolation = 'none', cmap = 'coolwarm')
  795. divider = make_axes_locatable(ax[1])
  796. cax = divider.append_axes('right', size='5%', pad=0.05)
  797. fig.colorbar(im, cax=cax)
  798. ax[1].set_title('HC-MDD (Enter Prob.)', fontsize = 18)
  799. ax[1].set_xticks([k for k in range(K)])
  801. plt.tight_layout()
  803. # %%
  804. # Initialize lists to collect data for each state-to-state transition and clinical score
  805. totst_, KS_, ids_, type_, sex_, age_, meds_ = [], [], [], [], [], [], []
  806. st_prob_, st_trans_ = [], []
  807. clnc_mes_, clnc_score_ = [], []
  809. # Loop over all subjects (healthy + MDD)
  810. for i, idx in enumerate(healthyID + mddID):
  811. # Assign group label and select corresponding transition matrix
  812. if idx in healthyID:
  813. typ = 'HC'
  814. trnst = transition_hc[i]
  815. else:
  816. typ = 'MDD'
  817. trnst = transition_mdd[int(i - len(healthyID))]
  819. # Loop over all clinical scores and their names
  820. for score, mes in zip([QIDS, MASQ_aa, MASQ_ad, MASQ_gd, RRS_dr, RRS_b, RRS_r],['QIDS', 'MASQ_AA', 'MASQ_AD', 'MASQ_GD', 'RRS_DR', 'RRS_BR', 'RRS_RF']):
  821. # Loop over all state-to-state transitions (excluding self-transitions)
  822. for x in range(trnst.shape[0]):
  823. for y in range(trnst.shape[0]):
  824. if x != y:
  825. # Subject and demographic info
  826. ids_.append(idx)
  827. type_.append(typ)
  828. sex_.append(SEX[idx])
  829. age_.append(AGE[idx])
  830. meds_.append(MEDS[idx])
  832. # Transition info
  833. st_trans_.append('%dto%d' % (x, y))
  834. st_prob_.append(trnst[y, x])
  836. # Clinical measure info
  837. clnc_mes_.append(mes)
  838. clnc_score_.append(score[idx])
  840. # Create DataFrame with one row per subject × transition × clinical score
  841. dfff = pd.DataFrame({'SubjectID': ids_, 'Morbidity': type_, 'Sex': sex_, 'Age': age_, 'Medication': meds_,
  842. 'STATETRANS': st_trans_, 'STATEPROB': st_prob_,'CLINICALMES': clnc_mes_, 'CLINICALSCORE': clnc_score_})
  844. # %%
  845. dfff
  847. # %%
  848. my_pal = {"HC": "darkgrey", "MDD": "deepskyblue"}
  849. st_trans = ['0to3', '3to0', '1to2', '2to1', '2to0', '3to1']
  850. fig,ax = plt.subplots(1,1,figsize = (16,5))
  851. sns.boxenplot(data = dfff[(dfff.STATETRANS.isin(st_trans))&(dfff.CLINICALMES == 'QIDS')], x = 'STATETRANS', y ='STATEPROB', hue = 'Morbidity' ,ax = ax, palette=my_pal)
  853. # %%
  854. # Initialize lists to store subject, state, and clinical data
  855. totst_, KS_, ids_, type_, sex_, age_, meds_ = [], [], [], [], [], [], []
  856. st_prob_, st_trans_ = [], []
  857. clnc_mes_, clnc_score_ = [], []
  859. # Loop over all subjects (healthy + MDD)
  860. for i, idx in enumerate(healthyID + mddID):
  861. # Assign group label and select corresponding transition matrix
  862. if idx in healthyID:
  863. typ = 'HC'
  864. trnst = transition_hc[i]
  865. else:
  866. typ = 'MDD'
  867. trnst = transition_mdd[int(i - len(healthyID))]
  868. # Loop over all clinical scores and their names
  869. for score, mes in zip([QIDS, MASQ_aa, MASQ_ad, MASQ_gd, RRS_dr, RRS_b, RRS_r],['QIDS','MASQ_AA','MASQ_AD','MASQ_GD', 'RRS_DR', 'RRS_BR', 'RRS_RF']):
  870. # Loop over each state
  871. for x in range(trnst.shape[0]):
  872. # Compute both EXIT (outgoing) and ENTER (incoming) probabilities
  873. for tty, transtype in enumerate(['EXIT','ENTER']):
  874. # Subject & demographic info
  875. ids_.append(idx)
  876. type_.append(typ)
  877. sex_.append(SEX[idx])
  878. age_.append(AGE[idx])
  879. meds_.append(MEDS[idx])
  881. # Transition type and probability
  882. st_trans_.append('%s %d' % (transtype, x))
  883. st_prob_.append(trnst.sum(axis=tty)[x] / 3) # normalized by number of states?
  885. # Clinical measure info
  886. clnc_mes_.append(mes)
  887. clnc_score_.append(score[idx])
  889. # Create a DataFrame with all subjects × states × transition types × clinical measures
  890. dfff_total = pd.DataFrame({'SubjectID': ids_, 'Morbidity': type_, 'Sex': sex_, 'Age': age_, 'Medication': meds_,
  891. 'STATETRANS': st_trans_, 'STATEPROB': st_prob_,
  892. 'CLINICALMES': clnc_mes_, 'CLINICALSCORE': clnc_score_})
  894. # %%
  895. fig,ax = plt.subplots(1,1,figsize = (24,5))
  896. sns.boxenplot(data = dfff_total, x = 'STATETRANS', y ='STATEPROB', hue = 'Morbidity' ,ax = ax, palette=my_pal)
  898. # %%
  899. fig,ax = plt.subplots(1,1,figsize = (10,4))
  900. sns.violinplot(data = dfff_total[(dfff_total.STATETRANS.isin(['EXIT 2', 'ENTER 2']))&(dfff_total.CLINICALMES=='QIDS')], x = 'STATETRANS', y ='STATEPROB', hue = 'Morbidity' ,ax = ax, palette=my_pal, density_norm='width')
  902. # %%
  905. # %%
  906. K = 4 # Number of states/clusters
  907. it = np.argmax([inertias[K, it] for it in range(niter)]) # Select iteration with maximal inertia
  908. k = 2 # Focus on state k=2
  909. my_pal = {"HC": "darkgrey", "MDD": "deepskyblue"} # Colors for plotting
  910. pvalshc = [] # Store p-values for HC
  911. pvalsmdd = [] # Store p-values for MDD
  912. fig, ax = plt.subplots(1, 4, figsize=(20, 4)) # Create figure with 1 row x 4 columns
  914. for i, mes in enumerate(['QIDS', 'MASQ_AD', 'MASQ_AA', 'MASQ_GD']): # Loop over clinical measures
  916. data = df[(df.totalstates == K) & (df.StateID == k) & (df.Morbidity == 'HC')] # Filter HC data for state k
  917. sns.regplot(data=data, x='DWELLTIME', y=mes, ax=ax[i], label='HC', color="darkgrey") # Plot regression for HC
  919. mask = ~np.isnan(data['DWELLTIME'].to_numpy()) & ~np.isnan(data[mes].to_numpy()) # Remove NaNs in case not every subject has their clinical assesment completed
  920. X = data[[mes, 'Sex', 'Age', 'Medication']][mask] # Design matrix: clinical + covariates
  921. X = sm.add_constant(X) # Add intercept
  922. y = data['DWELLTIME'][mask] # Dependent variable
  924. model = sm.OLS(y, X).fit() # Fit OLS regression
  925. print(model.summary()) # Print model summary
  926. pvalshc.append(model.pvalues.iloc[1]) # Save p-value for clinical measure coefficient
  928. data = df[(df.totalstates == K) & (df.StateID == k) & (df.Morbidity == 'MDD')] # Filter MDD data for state k
  929. sns.regplot(data=data, x='DWELLTIME', y=mes, ax=ax[i], label='MDD', color="deepskyblue") # Plot regression for MDD
  931. mask = ~np.isnan(data['DWELLTIME'].to_numpy()) & ~np.isnan(data[mes].to_numpy()) # Remove NaNs in case not every subject has their clinical assesment completed
  932. X = data[[mes, 'Sex', 'Age', 'Medication']][mask] # Design matrix
  933. X = sm.add_constant(X) # Add intercept
  934. y = data['DWELLTIME'][mask] # Dependent variable
  936. model = sm.OLS(y, X).fit() # Fit OLS regression
  937. print(model.summary()) # Print model summary
  938. pvalsmdd.append(model.pvalues.iloc[1]) # Save p-value for clinical measure coefficient
  940. ax[i].set_title(mes, fontsize=20) # Add subplot title with clinical measure
  942. plt.tight_layout() # Adjust spacing between subplots
  944. # %%
  945. fdrcorrection(pvalsmdd)
  947. # %%
  948. pvalshc = [] # Store p-values for HC
  949. pvalsmdd = [] # Store p-values for MDD
  951. st_trans = '3to0' # State transition of interest
  953. K = 4 # Number of states/clusters
  954. it = np.argmax([inertias[K, it] for it in range(niter)]) # Select iteration with maximal inertia
  955. fig, ax = plt.subplots(1, 7, figsize=(40, 4)) # Create figure with 1 row x 7 columns
  957. for i, mes in enumerate(clinical_measures): # Loop over all clinical measures
  958. # --- Healthy Controls ---
  959. data = dfff[(dfff['Morbidity'] == "HC") & (dfff['CLINICALMES'] == mes) & (dfff['STATETRANS'] == st_trans)] # Filter HC data
  960. sns.regplot(data=data, x='STATEPROB', y='CLINICALSCORE', ax=ax[i], label='HC', color="darkgrey") # Plot regression
  962. mask = ~np.isnan(data['STATEPROB'].to_numpy()) & ~np.isnan(data['CLINICALSCORE'].to_numpy()) # Remove NaNs
  963. X = data[['CLINICALSCORE', 'Sex', 'Age', 'Medication']][mask] # Design matrix with covariates
  964. X = sm.add_constant(X) # Add intercept
  965. y = data['STATEPROB'][mask] # Dependent variable
  967. model = sm.OLS(y, X).fit() # Fit OLS regression
  968. pvalshc.append(model.pvalues.iloc[1]) # Save p-value of clinical measure coefficient
  970. # --- MDD Subjects ---
  971. data = dfff[(dfff['Morbidity'] == "MDD") & (dfff['CLINICALMES'] == mes) & (dfff['STATETRANS'] == st_trans)] # Filter MDD data
  972. sns.regplot(data=data, x='STATEPROB', y='CLINICALSCORE', ax=ax[i], label='MDD', color="deepskyblue") # Plot regression
  974. mask = ~np.isnan(data['STATEPROB'].to_numpy()) & ~np.isnan(data['CLINICALSCORE'].to_numpy()) # Remove NaNs
  975. X = data[['CLINICALSCORE', 'Sex', 'Age', 'Medication']][mask] # Design matrix
  976. X = sm.add_constant(X) # Add intercept
  977. y = data['STATEPROB'][mask] # Dependent variable
  979. model = sm.OLS(y, X).fit() # Fit OLS regression
  980. pvalsmdd.append(model.pvalues.iloc[1]) # Save p-value of clinical measure coefficient
  982. ax[i].set_title(mes, fontsize=20) # Set subplot title
  983. ax[i].set_xlabel('Transition Prob. from %s' % st_trans, fontsize=15) # Label x-axis
  985. print(model.summary()) # Print model summary
  987. # %%
  988. fdrcorrection(pvalshc),pvalshc
  990. # %% [markdown]
  991. # ## Network Control Theory
  993. # %%
  994. np.random.seed(42) # reproducibility
  996. # Generate adjacency matrices for healthy subjects
  997. healthyDTI = {f"{500+i}": random_adjacency(size) for i in range(n_healthy)}
  999. # Generate adjacency matrices for MDD subjects
  1000. mddDTI = {f"{600+i}": random_adjacency(size) for i in range(n_mdd)}
  1002. allDTI = all_adj = {**healthyDTI, **mddDTI}
  1004. # %%
  1005. K = 4 # Number of states/clusters
  1006. it = np.argmax([inertias[K, it] for it in range(niter)]) # Select iteration with maximal inertia
  1008. norm_states = {} # Dictionary to store normalized states
  1009. for k in range(K):
  1010. norm_states[k] = normalize_state(medoids[K, it][k]) # Normalize the k-th state (z-score or other scaling)
  1012. # Compute network control theory metrics for all subjects
  1013. # Returns total/control energies, node-level energies, state trajectories, and control signals for both groups
  1014. totalenergies, node_energies, state_trajectory, control_signals, totalenergies_pers, node_energies_pers, state_trajectory_pers, control_signals_pers = NCT_multi(healthyID, mddID, allDTI, norm_states, K, size) # Run NCT analysis
  1016. # %%
  1017. ids_, type_, sex_, age_, meds_ = [], [], [], [], [] # Initialize subject and demographic lists
  1018. roi_, st_trans_, ne_ = [], [], [] # Initialize ROI, state transition, and node energy lists
  1019. clnc_mes_, clnc_score_ = [], [] # Initialize clinical measure lists
  1021. for idx in healthyID + mddID: # Loop over all subjects
  1022. if idx in healthyID: group = 'HC' # Assign group label
  1023. elif idx in mddID: group = 'MDD'
  1025. # Loop over all clinical measures and their names
  1026. for score, mes in zip([QIDS, MASQ_aa, MASQ_ad, MASQ_gd, RRS_dr, RRS_b, RRS_r],['QIDS', 'MASQ_AA', 'MASQ_AD', 'MASQ_GD', 'RRS_DR', 'RRS_BR', 'RRS_RF']):
  1028. # Loop over all state-to-state transitions
  1029. for k1 in range(K):
  1030. for k2 in range(K):
  1031. if k1 != k2: # Between-state transitions
  1032. for n in range(size): # Loop over all ROIs
  1033. ids_.append(idx) # Subject ID
  1034. sex_.append(SEX[idx]) # Sex
  1035. age_.append(AGE[idx]) # Age
  1036. meds_.append(MEDS[idx]) # Medication
  1038. type_.append(group) # Group
  1039. st_trans_.append('%dto%d' % (k1, k2)) # State transition label
  1040. roi_.append(n) # ROI index
  1041. ne_.append(node_energies[idx][k2, k1, n]) # Node-level energy
  1043. clnc_mes_.append(mes) # Clinical measure name
  1044. clnc_score_.append(score[idx]) # Clinical score
  1045. else: # Persistence (within-state) transitions
  1046. for n in range(size):
  1047. ids_.append(idx)
  1048. sex_.append(SEX[idx])
  1049. age_.append(AGE[idx])
  1050. meds_.append(MEDS[idx])
  1052. type_.append(group)
  1053. st_trans_.append('%dto%d' % (k1, k2))
  1054. roi_.append(n)
  1055. ne_.append(node_energies_pers[idx][k2, k1, n]) # Node energy for persistence
  1057. clnc_mes_.append(mes)
  1058. clnc_score_.append(score[idx])
  1060. # EXIT and ENTER total energies for each state
  1061. for t, total in enumerate(['EXIT', 'ENTER']):
  1062. for n in range(size):
  1063. ids_.append(idx)
  1064. sex_.append(SEX[idx])
  1065. age_.append(AGE[idx])
  1066. meds_.append(MEDS[idx])
  1068. type_.append(group)
  1069. st_trans_.append('%s%d' % (total, k1)) # Label EXIT or ENTER
  1070. roi_.append(n)
  1071. ne_.append(node_energies[idx].sum(axis=t)[k1, n] / 3) # Average node energy across states?
  1073. clnc_mes_.append(mes)
  1074. clnc_score_.append(score[idx])
  1076. # Create DataFrame with all node-level energies and associated info
  1077. dfff_nct = pd.DataFrame({'SubjectID': ids_, 'Morbidity': type_, 'Sex': sex_, 'Age': age_, 'Medication': meds_,
  1078. 'STATETRANS': st_trans_, 'ROI': roi_, 'NODEENERGIES': ne_,
  1079. 'CLINICALMES': clnc_mes_, 'CLINICALSCORE': clnc_score_})
  1081. # %%
  1082. dfff_nct
  1084. # %% [markdown]
  1085. # ## Compare control energies between groups for given ROIS
  1087. # %%
  1088. ROIs = [np.random.randint(40) for _ in range(2)]
  1090. fig, ax = plt.subplots(2,8, figsize=(24,6))
  1091. for s,st_trans in enumerate(['3to0', '2to1', 'EXIT0', 'ENTER0', 'EXIT2', 'ENTER2', 'EXIT3', 'ENTER3']):
  1092. for n,roi in enumerate(ROIs):
  1093. sns.violinplot(data = dfff_nct[(dfff_nct.ROI == roi)&(dfff_nct.CLINICALMES=='QIDS')&(dfff_nct.STATETRANS == st_trans)], y = 'NODEENERGIES', ax = ax[n][s], hue = 'Morbidity', palette=my_pal, split = True)
  1094. ax[0][s].set_title(f"Transition %s "%st_trans, fontsize = 16)
  1095. ax[n][s].tick_params(axis='both', which='major', labelsize=10)
  1096. ax[n][0].set_ylabel(f'Node Energies\n\n %s'%df_info[df_info.ROI=='%s'%roi].Name.iloc[0], fontsize = 12)
  1097. plt.tight_layout()
  1099. # %%
  1100. ids_, type_, sex_, age_, meds_ = [], [], [], [], [] # Initialize subject and demographic lists
  1101. st_trans_, te_, tp_ = [], [], [] # Initialize state transition, total energy, and transition probability lists
  1102. clnc_mes_, clnc_score_ = [], [] # Initialize clinical measure lists
  1104. for idx in healthyID + mddID: # Loop over all subjects
  1105. if idx in healthyID: group = 'HC' # Assign group label
  1106. elif idx in mddID: group = 'MDD'
  1108. # Loop over selected clinical measures
  1109. for score, mes in zip([QIDS, MASQ_aa, MASQ_ad, MASQ_gd], ['QIDS', 'MASQ_AA', 'MASQ_AD', 'MASQ_GD']):
  1111. # Loop over all state-to-state transitions
  1112. for k1 in range(K):
  1113. for k2 in range(K):
  1114. if k1 != k2: # Between-state transitions
  1115. ids_.append(idx) # Subject ID
  1116. sex_.append(SEX[idx]) # Sex
  1117. age_.append(AGE[idx]) # Age
  1118. meds_.append(MEDS[idx]) # Medication
  1120. type_.append(group) # Group label
  1121. st_trans_.append('%dto%d' % (k1, k2)) # State transition label
  1122. te_.append(totalenergies[idx][k2, k1]) # Total energy for transition
  1123. tp_.append(transition_probs[idx][k2, k1]) # Transition probability
  1125. clnc_mes_.append(mes) # Clinical measure name
  1126. clnc_score_.append(score[idx]) # Clinical score
  1127. else: # Persistence (within-state)
  1128. ids_.append(idx)
  1129. sex_.append(SEX[idx])
  1130. age_.append(AGE[idx])
  1131. meds_.append(MEDS[idx])
  1133. type_.append(group)
  1134. st_trans_.append('%dto%d' % (k1, k2))
  1135. te_.append(totalenergies_pers[idx][k2, k1]) # Energy for persistence
  1136. tp_.append(persistence_probs[idx][k1]) # Persistence probability
  1138. clnc_mes_.append(mes)
  1139. clnc_score_.append(score[idx])
  1141. # EXIT and ENTER total energies/probabilities for each state
  1142. for k1 in range(K):
  1143. for t, total in enumerate(['EXIT', 'ENTER']):
  1144. ids_.append(idx)
  1145. sex_.append(SEX[idx])
  1146. age_.append(AGE[idx])
  1147. meds_.append(MEDS[idx])
  1149. type_.append(group)
  1150. st_trans_.append('%s%d' % (total, k1)) # Label EXIT or ENTER
  1151. te_.append(totalenergies[idx].sum(axis=t)[k1] / 3) # Average total energy
  1152. tp_.append(transition_probs[idx].sum(axis=t)[k1] / 3) # Average probability
  1154. clnc_mes_.append(mes)
  1155. clnc_score_.append(score[idx])
  1157. # Create DataFrame with energies, probabilities, and clinical measures
  1158. dfff_energies = pd.DataFrame({'SubjectID': ids_, 'Morbidity': type_, 'Sex': sex_, 'Age': age_, 'Medication': meds_,
  1159. 'STATETRANS': st_trans_, 'STATEENERGY': te_, 'STATEPROB': tp_,
  1160. 'CLINICALMES': clnc_mes_, 'CLINICALSCORE': clnc_score_})
  1162. # %%
  1163. fig,ax = plt.subplots(1,2,figsize = (12,5))
  1165. arrhc = np.array([totalenergies[idx] for idx in healthyID])
  1166. arrmdd = np.array([totalenergies[idx] for idx in mddID])
  1168. im = ax[0].imshow(arrhc.mean(axis = 0) - arrmdd.mean(axis = 0), aspect = 'auto', interpolation = 'none', cmap = 'coolwarm', vmin = -5, vmax = 5)
  1169. divider = make_axes_locatable(ax[0])
  1170. cax = divider.append_axes('right', size='5%', pad=0.05)
  1171. fig.colorbar(im, cax=cax)
  1172. ax[0].set_title('HC-MDD (Transition Energy)', fontsize = 18)
  1173. ax[0].set_xticks([k for k in range(K)])
  1174. ax[0].set_yticks([k for k in range(K)])
  1175. ax[0].set_xlabel('Current State', fontsize = 15)
  1176. ax[0].set_ylabel('Next State', fontsize = 15)
  1178. arrhc_pers = np.array([totalenergies_pers[idx] for idx in healthyID])
  1179. arrmdd_pers = np.array([totalenergies_pers[idx] for idx in mddID])
  1181. im = ax[1].imshow(arrhc_pers.mean(axis = 0) - arrmdd_pers.mean(axis = 0), aspect = 'auto', interpolation = 'none', cmap = 'coolwarm', vmin = -5, vmax = 5)
  1182. divider = make_axes_locatable(ax[1])
  1183. cax = divider.append_axes('right', size='5%', pad=0.05)
  1184. fig.colorbar(im, cax=cax)
  1185. ax[1].set_title('HC-MDD (Persistence Energy)', fontsize = 18)
  1186. ax[1].set_xticks([k for k in range(K)])
  1187. ax[1].set_yticks([k for k in range(K)])
  1188. ax[1].set_xlabel('Current State', fontsize = 15)
  1189. ax[1].set_ylabel('Next State', fontsize = 15)
  1191. plt.tight_layout()
  1193. # %%
  1194. fig,ax = plt.subplots(1,2,figsize = (12,2))
  1195. arrhc = np.array([totalenergies[idx] for idx in healthyID])
  1196. arrmdd = np.array([totalenergies[idx] for idx in mddID])
  1198. im = ax[0].imshow((arrhc.mean(axis = 0).sum(axis=0)/3 - arrmdd.mean(axis = 0).sum(axis=0)/3).reshape(1,4), aspect = 'auto', interpolation = 'none', cmap = 'coolwarm', vmin = -5, vmax = 5)
  1199. divider = make_axes_locatable(ax[0])
  1200. cax = divider.append_axes('right', size='5%', pad=0.05)
  1201. fig.colorbar(im, cax=cax)
  1202. ax[0].set_title('HC-MDD (Exit Energy)', fontsize = 18)
  1203. ax[0].set_xticks([k for k in range(K)])
  1205. im = ax[1].imshow((arrhc.mean(axis = 0).sum(axis=1)/3 - arrmdd.mean(axis = 0).sum(axis=1)/3).reshape(1,4), aspect = 'auto', interpolation = 'none', cmap = 'coolwarm', vmin = -5, vmax = 5)
  1206. divider = make_axes_locatable(ax[1])
  1207. cax = divider.append_axes('right', size='5%', pad=0.05)
  1208. fig.colorbar(im, cax=cax)
  1209. ax[1].set_title('HC-MDD (Enter Energy)', fontsize = 18)
  1210. ax[1].set_xticks([k for k in range(K)])
  1212. plt.tight_layout()
  1214. # %%
  1215. st_trans = ['0to1','0to2', '0to3', '1to2', '2to1', '2to0', '0to2', '1to3', '3to1', '1to0', '0to1', '3to0','3to2', '2to3']
  1216. my_pal = {"HC": "darkgrey", "MDD": "deepskyblue"}
  1217. fig,ax = plt.subplots(1,1,figsize= (24,5))
  1218. sns.boxenplot(data = dfff_energies[(dfff_energies.STATETRANS.isin(st_trans))&(dfff_energies.CLINICALMES =='QIDS')], x = 'STATETRANS', y ='STATEENERGY', ax = ax, hue = 'Morbidity', palette=my_pal)
  1219. ax.tick_params(axis='both', which='major', labelsize=16)
  1221. # %%
  1222. my_pal = {"HC": "darkgrey", "MDD": "deepskyblue"}
  1224. st_trans1 = ['EXIT0', 'EXIT1', 'EXIT2', 'EXIT3','ENTER0','ENTER1', 'ENTER2','ENTER3']
  1225. fig,ax = plt.subplots(1,1,figsize= (24,5))
  1226. sns.boxplot(data = dfff_energies[(dfff_energies.STATETRANS.isin(st_trans1))&(dfff_energies.CLINICALMES =='QIDS')], x = 'STATETRANS', y ='STATEENERGY', ax = ax, palette=my_pal, hue = 'Morbidity')
  1228. # %%
  1229. ids_, type_ = [], [] # Initialize lists for subject IDs and group labels
  1230. st_trans_, te_, states_ = [], [], [] # Initialize lists for transition type, total energy, and state ID
  1232. for idx in healthyID + mddID: # Loop over all subjects
  1233. if idx in healthyID: group = 'HC' # Assign group label
  1234. elif idx in mddID: group = 'MDD'
  1235. else: print('WTF') # Safety check
  1237. for k1 in range(K): # Loop over each state
  1238. for t, total in enumerate(['EXIT', 'ENTER']): # Loop over EXIT and ENTER transitions
  1239. ids_.append(idx) # Append subject ID
  1240. type_.append(group) # Append group label
  1241. states_.append(k1) # Append state ID
  1242. st_trans_.append('%s' % (total)) # Append transition type
  1243. te_.append(totalenergies[idx].sum(axis=t)[k1] / 3) # Compute and append average total energy for this transition
  1245. # Create DataFrame with EXIT/ENTER total energies
  1246. dfff_energies_exitenter = pd.DataFrame({'SubjectID': ids_, 'Morbidity': type_, 'STATEID': states_,'STATETRANS': st_trans_, 'STATEENERGY': te_})
  1248. my_pal23 = {"EXIT": "darksalmon", "ENTER": "lightseagreen"} # Colors for violin plot
  1249. sns.violinplot(data=dfff_energies_exitenter, x='STATEID', y='STATEENERGY',palette=my_pal23, hue='STATETRANS', cut=True, split=True, density_norm='area') # Plot violin plot of EXIT vs ENTER energies by state
  1251. st_trans = ['EXIT%d', 'ENTER%d'] # Template strings for per-state energy comparisons
  1253. for k in range(4): # Loop over each state
  1254. print(k) # Print state index
  1255. # Compute mean difference between EXIT and ENTER energies for QIDS
  1256. print(dfff_energies[(dfff_energies.STATETRANS == 'EXIT%d' % k) & (dfff_energies.CLINICALMES == 'QIDS')].STATEENERGY.to_numpy().mean() - dfff_energies[(dfff_energies.STATETRANS == 'ENTER%d' % k) & (dfff_energies.CLINICALMES == 'QIDS')].STATEENERGY.to_numpy().mean())
  1257. # Paired t-test between EXIT and ENTER energies for this state
  1258. print(ttest_rel(dfff_energies[(dfff_energies.STATETRANS == 'EXIT%d' % k) & (dfff_energies.CLINICALMES == 'QIDS')].STATEENERGY.to_numpy(),dfff_energies[(dfff_energies.STATETRANS == 'ENTER%d' % k) & (dfff_energies.CLINICALMES == 'QIDS')].STATEENERGY.to_numpy()))
  1260. # %%
  1261. my_pal = {"HC": "darkgrey", "MDD": "deepskyblue"} # Color palette for groups
  1262. fig, ax = plt.subplots(1, 4, figsize=(22, 5)) # Create figure with 1 row x 4 columns
  1264. for st, st_trans in enumerate(['3to0', '2to1', 'EXIT2', "EXIT3"]): # Loop over selected state transitions
  1265. ax[st].set_title('%s' % st_trans, fontsize=16) # Initial subplot title
  1267. for group in ['HC', 'MDD']: # Loop over groups
  1268. # Filter data for current transition, clinical measure, and group
  1269. data = dfff_energies[(dfff_energies.STATETRANS == st_trans) & (dfff_energies.CLINICALMES == 'QIDS') & (dfff_energies.Morbidity == group)]
  1271. # Plot regression: Transition probability vs transition energy
  1272. sns.regplot(data=data, x='STATEENERGY', y='STATEPROB', ax=ax[st], label=group, color=my_pal[group])
  1274. # Prepare design matrix for linear regression
  1275. X = data[['STATEENERGY']] # Independent variable
  1276. X = sm.add_constant(X) # Add intercept
  1277. y = data['STATEPROB'] # Dependent variable
  1279. linear_model = sm.OLS(y, X).fit() # Fit OLS regression
  1280. print(linear_model.summary()) # Print regression summary
  1282. # Update subplot title with more descriptive text
  1283. ax[st].set_title('Trans. Prob. vs Trans. Energy for %s' % (st_trans), fontsize=20)
  1285. plt.tight_layout() # Adjust spacing between subplots
  1287. # %%
  1290. # %%
  1293. # %%
  1296. # %%
  1299. # %%
  1302. # %%
  1305. # %%

tutorial.ipynb at commit 8f93894, under GPL-3.0 · at the source

Overview

Authors: B Ülgen Kilic1,2, Jenna Jubeir1,3, Priti Balchandani3,4,5,6, James W Murrough1,3,7, Laurel S Morris1,8,9, Yael Jacob1,3,4,5
ORCID iDs: B Ülgen Kilic, Jenna Jubeir, James W Murrough, Laurel S Morris, Yael Jacob
  1. Department of Psychiatry, The Dennis S. Charney, MD, Depression and Anxiety Discovery Center, Icahn School of Medicine at Mount Sinai, New York, NY USA
  2. Department of Psychiatry, Center for Computational Psychiatry, Icahn School of Medicine at Mount Sinai, New York, NY USA
  3. Nash Family Department of Neuroscience & Friedman Brain Institute, Icahn School of Medicine at Mount Sinai, New York, NY USA
  4. Department of Radiology, BioMedical Engineering and Imaging Institute, Icahn School of Medicine at Mount Sinai, New York, NY USA
  5. Center for Engineering and Precision Medicine, Icahn School of Medicine at Mount Sinai & Rensselaer Polytechnic Institute, New York, NY USA
  6. Department of Diagnostic, Molecular and Interventional Radiology, Icahn School of Medicine at Mount Sinai, New York, NY USA
  7. VISN 2 Mental Illness Research, Education and Clinical Center (MIRECC), James J. Peters VA Medical Center, Bronx, NY USA
  8. Nuffield Department of Clinical Neurosciences, University of Oxford, Oxford, UK
  9. Department of Experimental Psychology, University of Oxford, Oxford, UK
Institutions: Icahn School of Medicine at Mount Sinai (United States); Rensselaer Polytechnic Institute (United States); James J. Peters VA Medical Center (United States); Mental Illness Research, Education and Clinical Centers (United States); University of Oxford (United Kingdom)
Journal: Nature communications, volume 17, issue 1, article 5662
Dates: received 6 August 2025; accepted 3 April 2026; published online 23 April 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-71961-4 · PMID 42026043 · PMCID PMC13315732 · OpenAlex W7155399019
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism), depression (population), computational (subfield)
Methods: Connectivity, Statistics, Machine learning, Smoothing, state filtering, decompositions, fMRI & imaging
Keywords: Dynamical systems, Depression, Control theory, Network models
MeSH: Brain*, Major Depressive Disorder*, Brain Mapping, Diffusion Magnetic Resonance Imaging, Female, Humans, Magnetic Resonance Imaging, Male, Severity of Illness Index (* major topic)
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: NIMH NIH HHS (K01 MH131855, R01 MH109544, K01 MH120433); NIH HHS (S10 OD026880, S10 OD030463); U.S. Department of Health &amp; Human Services | National Institutes of Health (R01 MH109544); U.S. Department of Health &amp; Human Services | NIH | National Institute of Mental Health (K01MH12043); NCATS NIH HHS (UL1 TR004419); U.S. Department of Health & Human Services | NIH | National Institute of Mental Health (NIMH) (K01MH12043)
Citations: cited by 1 paper (Europe PMC); 88 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

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

prantikk/me-ica

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 8cc47cfed0203b3d6d187935ad3c2823b3e36a88, 20 January 2018
Languages: Python (240)
Size: 266 files, 240 scripts
Software Heritage: archived
Found in: the text, “MRI data preprocessing”
Holds: README, tests
Not found: license file, CITATION.cff, environment file, continuous integration, documentation
Tools: NumPy (77 files), SciPy (10 files), NiBabel (4 files), AFNI (3 files), pydicom (2 files), scikit-learn (2 files), Pillow (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
241 files
At the source: bitbucket.org/prantikk/me-ica

ulgenklc/Brain_states

License: GPL-3.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 8f93894ad16451c3f0668a8d7b1d6bd38b38c529, 16 May 2026
Languages: Jupyter (1)
Size: 63 files, 1 script
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (environment.yml), 1 notebook
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: Matplotlib (1 file), NiBabel (1 file), Nilearn (1 file), NumPy (1 file), pandas (1 file), SciPy (1 file), seaborn (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
3 files
At the source: github.com/ulgenklc/Brain_states

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41467-026-71961-4.

Tracing map

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

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 241 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 statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • it says that the data are available on request

Read it in the paper: doi.org/10.1038/s41467-026-71961-4.

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

Recorded: type, language, journal, volume, issue, pages, dates, 6 authors, 4 keywords, 9 MeSH terms, 6 funders, 83 references.

Cite

This paper

Kilic, B. Ü., Jubeir, J., Balchandani, P., Murrough, J. W., Morris, L. S., & Jacob, Y. (2026). Spatiotemporal asymmetries on brain energy landscape uncover system entrapment related to depression severity. Nature communications, 17(1), 5662. https://doi.org/10.1038/s41467-026-71961-4

BibTeX

@article{kilic2026spatiotemporal,
author = {Kilic, B Ülgen and Jubeir, Jenna and Balchandani, Priti and Murrough, James W and Morris, Laurel S and Jacob, Yael},
title = {{Spatiotemporal asymmetries on brain energy landscape uncover system entrapment related to depression severity}},
journal = {Nature communications},
year = {2026},
month = apr,
volume = {17},
number = {1},
pages = {5662},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-71961-4},
url = {https://doi.org/10.1038/s41467-026-71961-4},
pmid = {42026043},
pmcid = {PMC13315732}
}

RIS

TY - JOUR
AU - Kilic, B Ülgen
AU - Jubeir, Jenna
AU - Balchandani, Priti
AU - Murrough, James W
AU - Morris, Laurel S
AU - Jacob, Yael
TI - Spatiotemporal asymmetries on brain energy landscape uncover system entrapment related to depression severity
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/04/23
VL - 17
IS - 1
SP - 5662
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-71961-4
UR - https://doi.org/10.1038/s41467-026-71961-4
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-71961-4",
"type": "article-journal",
"title": "Spatiotemporal asymmetries on brain energy landscape uncover system entrapment related to depression severity",
"container-title": "Nature communications",
"author": [
{
"family": "Kilic",
"given": "B Ülgen"
},
{
"family": "Jubeir",
"given": "Jenna"
},
{
"family": "Balchandani",
"given": "Priti"
},
{
"family": "Murrough",
"given": "James W"
},
{
"family": "Morris",
"given": "Laurel S"
},
{
"family": "Jacob",
"given": "Yael"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "5662",
"DOI": "10.1038/s41467-026-71961-4",
"PMID": "42026043",
"PMCID": "PMC13315732",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-71961-4",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
23
]
]
}
}

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/s41398-026-04025-2 [code]
Brain energetic landscapes shape state dysregulation in major depressive disorder: a morphological network controllability perspective.
Journal: Translational psychiatry
In common: AFNI, Nilearn, NiBabel, 6 other tools, depression, 13 references
[2] doi:10.1162/imag.a.1282 [code]
Metabolic syndrome severity and the energetic cost of brain network transitions: A normative modeling study of accelerated brain aging.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Nilearn, NiBabel, statsmodels, 6 other tools, computational, structural MRI / diffusion, 9 references
[3] doi:10.1002/hbm.70600 [code]
A Data-Driven Closed-Loop Control Approach to Drive Neural State Transitions for Mechanistic Insight.
Journal: Human brain mapping
In common: Nilearn, NiBabel, statsmodels, 5 other tools, computational, depression, 6 references
[4] doi:10.1002/hbm.70485 [code]
Exploring the Role of the Rich Club in Network Control of Neurocognitive States.
Journal: Human brain mapping
In common: Nilearn, NiBabel, scikit-learn, 3 other tools, 8 references
[5] doi:10.1038/s41467-026-76011-7 [code]
Human cortex organizes dynamic co-fluctuations along the sensorimotor-association axis.
Journal: Nature communications
In common: AFNI, Pillow, NiBabel, 3 other tools, 8 references
[6] doi:10.1038/s41467-026-72931-6 [code]
Three parsimonious spatiotemporal patterns in cerebellum reveal individual traits in function and behavior.
Journal: Nature communications
In common: Nilearn, NiBabel, statsmodels, 6 other tools, 4 references
[7] doi:10.1038/s41467-026-75585-6 [code]
Brain network dynamics reflect psychiatric illness status and transdiagnostic symptom profiles across health and disease.
Journal: Nature communications
In common: seaborn, scikit-learn, pandas, 3 other tools, depression, 6 references
[8] doi:10.1038/s41467-026-71151-2 [code]
Common and distinct neural correlates of social interaction processing and theory of mind in narratives.
Journal: Nature communications
In common: AFNI, Nilearn, Pillow, 8 other tools, 2 references
[9] doi:10.1038/s41467-026-75745-8 [code]
A language network in the individualized functional connectomes of 1199 human brains doing arbitrary tasks.
Journal: Nature communications
In common: Nilearn, Pillow, NiBabel, 6 other tools, 4 references
[10] doi:10.1038/s41531-026-01354-3 [code]
Neuromodulation-induced normalization of cortical metastable dynamics signatures in Parkinson's disease.
Journal: NPJ Parkinson's disease
In common: Nilearn, NiBabel, statsmodels, 6 other tools, 4 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.