OSCR

Anterior insular co-activation patterns associated with stress markers in chronic primary pain.

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] § Materials and methods › Resting-state functional dynamics › Co-activation patterns ↔ 02_CAPs/01_Analysis/Script_twopop_PCA_HC_ref_SH.m, lines 824–916 · score 0.87 · consensus clustering, activation patterns, reference population, events, PAC, retained frames
  2. [2] § Materials and methods › Demographic and clinical characteristics ↔ 01_Clinical_Measures/Demographics_Table.R, lines 87–127 · score 0.70 · hormonal contraceptives, psychotropic medications, corticosteroids, analgesics, sex, demographic
  3. [3] § Results › Identified aIC-based CAPs in HC ↔ 02_CAPs/01_Analysis/Script_twopop_PCA_HC_ref_SH.m, lines 824–916 · score 0.69 · MNI space, co activation, reference population, retained frames, fraction, maps
  4. [4] § Materials and methods › Statistical analysis ↔ 02_CAPs/02_Statistics/02_CAPs_Stat.Rmd, lines 315–334 · score 0.68 · Post hoc, linear mixed, interaction term, ID, FDR, BDI
  5. [5] § Materials and methods › Statistical analysis ↔ 01_Clinical_Measures/Demographics_Table.R, lines 129–167 · score 0.64 · Wilcoxon rank sum, Shapiro Wilk, Welch, variances, demographic, clinical
  6. [6] § Materials and methods › Demographic and clinical characteristics › Pain provocation test ↔ 01_Clinical_Measures/Demographics_Table.R, lines 211–249 · score 0.60 · peg algometry, pain intensity, ear, Algopeg
  7. [7] § Materials and methods › Statistical analysis ↔ 01_Clinical_Measures/Demographics_Table.R, lines 211–249 · score 0.60 · BPI severity, Peg algometry, Alpha, pain, CPP
  8. [8] § Materials and methods › Demographic and clinical characteristics › Salivary cortisol and alpha-amylase ↔ 01_Clinical_Measures/Cortisol_AUC_Analysis.R, lines 69–156 · score 0.58 · post awakening, cortisol, concentrations, curve, AUCI, 45 min
  9. [9] § Materials and methods › MRI data acquisition and preprocessing ↔ 02_CAPs/01_Analysis/Script_twopop_PCA_HC_ref_SH.m, lines 68–134 · score 0.56 · motion parameters, SPM, FD, voxel, brain
  10. [10] § Materials and methods › Resting-state functional dynamics › Co-activation patterns ↔ 02_CAPs/01_Analysis/Script_twopop_PCA_HC_ref_SH.m, lines 136–178 · score 0.51 · consensus clustering, selected frames, algorithm, threshold, activations, HC

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

MATLAB · 922 lines · 33 KB · no license · 4 matches

  1. %% Script to run the CAPs analyses without the GUI
  2. % In this script, we assume HC as a reference population
  3. % PCA included
  5. % Adapt:
  6. % - Seed
  7. % - The parameters in section 2 can be adaped.
  8. % - K_range in section 4, as well as in the section 7
  9. % "Parameters.KMeansClustering.MaxClusterNumber = Kmax"
  11. % Funtion generated by Samantha Weber:
  12. % - createMask_CAP_SW
  13. % - loadData_SW
  14. % - loadseed_SW
  15. % - MakeViolin_SW
  16. % Remaining Functions are from the TbCAPs Toolbox (https://github.com/MIPLabCH/TbCAPs)
  17. % Adaptation and additions for the CAPs in CPP Project by Salome Häuselmann
  18. % ([email hidden]) 2024/2025
  20. clear; clc;close;
  21. addpath(genpath(pwd));
  23. %% 1. Loading the data files
  24. RootPath = '';
  25. SavePath = '';
  26. spmPath = '';
  27. mkdir(SavePath);
  29. %% 1.1 Define your variables
  30. prefix = 's5w*';
  31. selectMask = 'Population'; % make GM Population Mask based on data
  32. seedName = {'Anterior_Insula_BNA_mask.nii'};
  33. MasksPath = '/';
  34. fancyName = '';
  35. Group = {'HC', 'CPP'};
  36. n_datasets = 2;
  38. % load HC folder;
  39. filelist=dir(fullfile(RootPath,Group{1}, 'P*'));
  40. N_HC{1}=size(filelist,1);
  41. myHC=cell(1,N_HC{1});
  42. for f=1:N_HC{1}
  43. myHC{f}=filelist(f).name;
  44. end
  45. clear filelist f
  47. % load CPP folder;
  48. filelist=dir(fullfile(RootPath,Group{2}, 'P*'));
  49. N_CPP{1}=size(filelist,1);
  50. myCPP=cell(1,N_CPP{1});
  51. for f=1:N_CPP{1}
  52. myCPP{f}=filelist(f).name;
  53. end
  54. clear filelist f
  56. mySubj = [myHC, myCPP];
  57. N_Subj{1} = N_HC{1} + N_CPP{1};
  59. % Define the folders HC
  60. for i = 1:N_HC{1}
  61. functdir1{i} = fullfile(RootPath, Group{1}, myHC{i});
  62. end
  64. % Define the folders CPP
  65. for i = 1 : N_CPP{1}
  66. functdir2{i} = fullfile(RootPath, Group{2}, myCPP{i});
  67. end
  68. %% 1. Loading the data files
  69. mask ={};
  70. brain_info ={};
  71. n_dataset = 0;
  73. % Data: cell array, each cell of size n_TP x n_masked_voxels
  74. % Mask: n_voxels x 1 logical vector
  75. % Header: the header (obtained by spm_vol) of one NIFTI file with proper
  76. % data dimension and .mat information
  77. [mask, brain_info] = createMask_CAP_SW(functdir1, prefix, selectMask);
  78. Tstart = clock;
  79. % Data HC: cell array, each cell of size n_TP x n_masked_voxels
  80. if exist(fullfile(SavePath,fancyName,['HCloaded_' fancyName '.mat']))
  81. disp(['Data HC has been loaded already !'])
  82. load(fullfile(SavePath,fancyName,['HCloaded_' fancyName '.mat']));
  83. else
  84. TC ={};
  85. FD ={};
  86. [TC,brain,FD] = loadData_SW(functdir1, prefix, mask, brain_info); % here the motion parameters rp_*-txt file are loaded in
  87. FD1 = FD{:,:};
  88. TC1 = TC{:,:};
  89. % Save intermediate steps
  91. if ~exist (fullfile(SavePath,fancyName))
  92. mkdir(fullfile(SavePath,fancyName));
  93. end
  94. save(fullfile(SavePath,fancyName,['HCloaded_' fancyName '.mat']),'TC1','brain','FD1','-v7.3')
  95. end
  97. % Data CPP: cell array, each cell of size n_TP x n_masked_voxels
  98. if exist(fullfile(SavePath,fancyName,['CPPloaded_' fancyName '.mat']))
  99. disp(['Data CPP has been loaded already !'])
  100. load(fullfile(SavePath,fancyName,['CPPloaded_' fancyName '.mat']));
  101. else
  102. TC ={};
  103. FD={};
  104. [TC,brain,FD] = loadData_SW(functdir2, prefix, mask, brain_info); % here the motion parameters rp_*-txt file are loaded in
  105. FD2 = FD{:,:};
  106. TC2 = TC{:,:};
  107. % Save intermediate steps
  109. if ~exist (fullfile(SavePath,fancyName))
  110. mkdir(fullfile(SavePath,fancyName));
  111. end
  112. save(fullfile(SavePath,fancyName,['CPPloaded_' fancyName '.mat']),'TC2','brain','FD2','-v7.3')
  113. end
  115. % Seed: a n_masked_voxels x n_seed logical vector with seed information
  117. if exist(fullfile(SavePath,fancyName,['Seed_' fancyName '.mat']))
  118. disp(' ');
  119. disp('-----------------------------------------------------');
  120. disp(['Seed has been loaded already !'])
  121. load(fullfile(SavePath,fancyName,['Seed_' fancyName '.mat']));
  122. else
  123. disp(' ');
  124. disp('-----------------------------------------------------');
  125. disp('Load seed...');
  126. [Seed] = loadseed_SW(functdir1, prefix, mask, brain_info, seedName,MasksPath);
  127. save(fullfile(SavePath,fancyName,['Seed_' fancyName '.mat']),'Seed');
  128. disp('Seed successfully loaded!');
  129. end
  131. % Computes seed maps for each subject and for the population, using the
  132. % data from the chosen reference population
  133. TC = [TC1, TC2];
  134. [~,AvgSeedMap] = CAP_Compute_SeedMap(TC1,Seed,1);
  136. %% 2. Specifying the main parameters
  138. % Threshold above which to select frames
  139. T = 0.84;
  141. % Selection mode ('Threshold' or 'Percentage')
  142. SelMode = 'Threshold';
  144. % Threshold of FD above which to scrub out the frame and also the t-1 and
  145. % t+1 frames (if you want another scrubbing setting, directly edit the
  146. % code)
  147. Tmot = 0.5;
  149. % Type of used seed information: select between 'Average','Union' or
  150. % 'Intersection'
  151. SeedType = 'Average';
  153. % Contains the information, for each seed (each row), about whether to
  154. % retain activation (1 0) or deactivation (0 1) time points
  155. switch SeedType
  156. case 'Union'
  157. Activation = [1,0];
  158. for s = 1:size(seedName,2)
  159. SignMatrix(s,:) = Activation;
  160. end
  161. case 'Average'
  162. SignMatrix = [1 0];
  163. end
  165. % Percentage of frames to use in each fold of consensus clustering
  166. Pcc = 80;
  168. % Number of folds we run consensus clustering for
  169. N = 200;%50; %120;
  171. % Percentage of positive-valued voxels to retain for clustering
  172. Pp = 100;
  174. % Percentage of negative-valued voxels to retain for clustering
  175. Pn = 100;
  177. % Number of repetitions of the K-means clustering algorithm
  178. n_rep = 50;%50;
  182. %% 3. Selecting the frames to analyse
  184. % Xon will contain the retained frames, and Indices will tag the time
  185. % points associated to these frames, for each subject (it contains a
  186. % subfield for retained frames and a subfield for scrubbed frames)
  187. [Xon1,p1,Indices1,idx_sep_seeds1,Xonp_scrub1] = CAP_find_activity(TC1,Seed,T,FD1,Tmot,SelMode,SeedType,SignMatrix);
  188. [Xon2,p2,Indices2,idx_sep_seeds2,Xonp_scrub2] = CAP_find_activity(TC2,Seed,T,FD2,Tmot,SelMode,SeedType,SignMatrix);
  190. % Percentage of retained frames across subjects
  191. RetainedPercentage{1} = p1(3,:);
  192. RetainedPercentage{2} = p2(3,:);
  194. % Indices of the frames that have been retained (used later for metrics
  195. % computations)
  196. FrameIndices{1} = Indices1;
  197. FrameIndices{2} = Indices2;
  199. tmp_toplot1 = ConcatMat(RetainedPercentage(1),1,1,N_HC,'FD');
  200. tmp_toplot2 = ConcatMat(RetainedPercentage(2),1,1,N_CPP,'FD');
  202. tmp_toplot = zeros(2,N_CPP{1});
  203. tmp_toplot(1,1:N_HC{1}) = tmp_toplot1;
  204. tmp_toplot(2,1:N_CPP{1}) = tmp_toplot2;
  205. tmp_toplot(tmp_toplot == 0) = NaN;
  207. %% Quality Checks 1 - Visualization
  209. % 1. Percentage of Retained Frames
  211. % Perform two-sample t-test between the two groups
  212. %[~, p_value, ci, stats] = ttest2(tmp_toplot1, tmp_toplot2, 'Vartype', 'unequal');
  214. % Perform unpaired non-parametric t-test (Mann-Withney U)(since in the
  215. % graph the do not look nomrmally distributed)
  216. [p_value, h, stats] = ranksum(tmp_toplot1, tmp_toplot2);
  218. % Displays the violin plot of subject scrubbing percentage for the
  219. % reference population
  220. TPViolin = figure;
  221. axes1 = axes('Parent',TPViolin);
  223. % Colors used in plotting of all populations
  224. PopColor{1} = [255,255,180; 219,224,252; 188,252,188; 230,230,230]/255;
  225. PopColor{2} = [130,48,48; 51,75,163; 59,113,86; 0, 0, 0]/255;
  227. %Create plot
  228. [~,~,TPViolin] = MakeViolin_SW(tmp_toplot,axes1,{'HC' 'CPP'},'Frames ret. [%]',PopColor,1,2);
  229. set(TPViolin,'Visible','on');
  231. % Add t-test results as text on the plot
  232. xPosition = mean(xlim); % Centered on x-axis
  233. yPosition = max(ylim) - 0.05 * range(ylim); % Slightly below the top of the y-axis
  234. text(xPosition, yPosition, sprintf('W = %.2f, p = %.4f', stats.ranksum, p_value), ...
  235. 'FontSize', 12, 'FontWeight', 'bold', 'HorizontalAlignment', 'center', 'Color', 'k');
  237. %Save plot as .jpg and .fig
  238. saveas(gcf,fullfile(SavePath,fancyName,'RetainedFrames.jpg'));
  239. saveas(gcf,fullfile(SavePath,fancyName,'RetainedFrames.fig'));
  240. disp('-----------------------------------------------------');
  241. disp(' ');
  242. disp(['Figure Frames retained has been saved in ' SavePath]);
  243. close;
  245. % 2. Distribution of retained frames
  246. % Test of data distribution homogeneity (inter-subject variability)
  247. % Test normal distribution around mean % of retained frames value using Kolmogorov–Smirnov test.
  249. % Define the mean and standard deviation for both datasets
  250. mu_HC = mean(tmp_toplot1);
  251. sigma_HC = std(tmp_toplot1);
  252. mu_CPP = mean(tmp_toplot2);
  253. sigma_CPP = std(tmp_toplot2);
  255. % Standardize the data (convert to Z-scores)
  256. z_HC = (tmp_toplot1 - mu_HC) / sigma_HC;
  257. z_CPP = (tmp_toplot2 - mu_CPP) / sigma_CPP;
  259. % Perform Kolmogorov-Smirnov Test against the standard normal distribution
  260. [h_HC, p_HC] = kstest(z_HC);
  261. [h_CPP, p_CPP] = kstest(z_CPP);
  263. % Coefficient of Variation (CV)
  264. cv_HC = sigma_HC / mu_HC;
  265. cv_CPP = sigma_CPP / mu_CPP;
  267. % Define the normal distribution curves for plotting
  268. x_HC = linspace(mu_HC - 4*sigma_HC, mu_HC + 4*sigma_HC, 1000);
  269. x_CPP = linspace(mu_CPP - 4*sigma_CPP, mu_CPP + 4*sigma_CPP, 1000);
  270. pdf_HC = normpdf(x_HC, mu_HC, sigma_HC);
  271. pdf_CPP = normpdf(x_CPP, mu_CPP, sigma_CPP);
  273. % Plot histogram and fitted normal distribution curves
  274. figure;
  275. hold on;
  276. histogram(tmp_toplot1, 20, 'Normalization', 'pdf', 'FaceAlpha', 0.5, 'DisplayName', 'HC (Histogram)');
  277. histogram(tmp_toplot2, 20, 'Normalization', 'pdf', 'FaceAlpha', 0.5, 'DisplayName', 'CPP (Histogram)');
  278. plot(x_HC, pdf_HC, 'LineWidth', 2, 'DisplayName', sprintf('Normal Fit HC (KS p=%.3f)', p_HC));
  279. plot(x_CPP, pdf_CPP, 'LineWidth', 2, 'DisplayName', sprintf('Normal Fit CPP (KS p=%.3f)', p_CPP));
  280. xlabel('Retained Frames Value');
  281. ylabel('Density');
  282. title('Data Distribution and Normal Fit (HC vs. CPP)');
  283. legend;
  284. hold off;
  286. % Save plot as .jpg and .fig with a new name
  287. saveas(gcf, fullfile(SavePath, fancyName, 'HC_CPP_Distribution.jpg'));
  288. saveas(gcf, fullfile(SavePath, fancyName, 'HC_CPP_Distribution.fig'));
  290. % Display message confirming the save
  291. disp('-----------------------------------------------------');
  292. disp(' ');
  293. disp(['Figure "HC_CPP_Distribution" has been saved in ' SavePath]);
  295. % Close the figure
  296. close;
  299. %% 4. Consensus clustering (if wished to determine the optimum K)
  301. % The input to PCA should have a dimensionality n_dimensions x
  302. % n_datapoints; since we consider a population of subjects in a cell array
  303. % (each cell with size n_dims x n_samples_persubj), we want to concatenate
  304. % these cells into one giant data matrix, which we feed to the PCA
  305. % function. I will denote the number of dimensions (or voxels) by V, and
  306. % the number of total time points across subjects by T.
  308. % Change dimension of Xon1 to feed into PCA
  309. Xon1_pca = [];
  311. for i = 1:N_HC{1,1}
  313. Xon1_pca = [Xon1_pca, Xon1{1,i}];
  315. end
  317. [U, W,Eigenvals,mu] = ComputePCA(Xon1_pca);
  318. %save(fullfile(SavePath,fancyName,'PCA_Outputs.mat'), 'U', 'W', 'Eigenvals');
  320. %[U,W] = ComputePCA(cell2mat(Xon1_pca));
  322. % After the call, U will have size V x T, and W will have size T x T. Be
  323. % careful that in our function, row i of W contain the weights, for all
  324. % T frames, associated to the principal direction i (contained as the i-th
  325. % column in U).
  327. % Also, notice that we want to compute only one PCA on the
  328. % population-wise data, not one per subject! Otherwise, we would have a
  329. % different dimensionality reduction for each subject, which would make our
  330. % life complicated...
  332. % You want to input "{W}" (rows = dimensions, columns = data points) to
  333. % CAP_ConsensusClustering. The brackets are because the function wants a
  334. % cell input...
  336. % After the consensus clustering step, you will want to feed W to your
  337. % k-means clustering as well instead of "cell2mat(Xon1)". Your output CAP
  338. % matrix will have size K x T, with K the number of CAPs. In order to go
  339. % back to the original space, you will apply:
  340. % CAP_original = (U*CAP')'
  341. % This will give you CAPs with a dimensionality K x V, i.e., what you would
  342. % have obtained using CAP analysis without PCA
  345. % This specifies the range of values over which to perform consensus
  346. % clustering: if you want to run parallel consensus clustering processes,
  347. % you should feed in different ranges to each call of the function
  348. K_range = 2:6;
  349. if exist(fullfile(SavePath,fancyName,['ConsensusClustering_Range' num2str(K_range(1)) '_to_' num2str(K_range(end)) '_' fancyName '.mat']))
  350. disp('-----------------------------------------------------');
  351. disp('Consensus Clustering was performed already !')
  352. disp('Load data...');
  353. load(fullfile(SavePath,fancyName,['ConsensusClustering_Range' num2str(K_range(1)) '_to_' num2str(K_range(end)) '_' fancyName '.mat']))
  354. disp(' ');
  355. disp('Data loaded successfully!');
  356. else
  357. % Have each of these run in a separate process on the server =)
  358. disp('Running Consensus Clustering...')
  359. % Umcomment if no PCA is performed
  360. %[Consensus] = CAP_ConsensusClustering(Xon1,K_range,'items',Pcc/100,N,'correlation');
  362. %Consensus clustering using PCA reduced dataset
  363. [Consensus] = CAP_ConsensusClustering({W},K_range,'items',Pcc/100,N,'correlation');
  365. %Plot Consensus Matrix
  366. ConsensusMatrixPlot(Consensus,SavePath,fancyName, K_range)
  368. % Calculates the quality metrics
  369. [CDF,PAC] = ComputeClusteringQuality(Consensus,K_range);
  371. % % Calculates the quality metrics
  372. % [~,Qual] = ComputeClusteringQuality(Consensus,[]);
  373. save(fullfile(SavePath,fancyName,['ConsensusClustering_Range' num2str(K_range(1)) '_to_' num2str(K_range(end)) '_' fancyName '.mat']))
  374. disp(' ');
  375. disp('Consensus Cluster successfully performed!');
  376. end
  378. Ttotal=etime(clock, Tstart);
  379. disp(['** Consensus clustering completed. Total time: ' num2str(Ttotal/60,'%3.1f') ' min.']);
  381. % Qual should be inspected to determine the best cluster number(s)
  382. CCPlot = figure;
  383. ax1 = axes(CCPlot);
  385. %set(CCPlot,'Visible','on');
  386. tmp_plot = bar(2:K_range(end),1-PAC);
  387. xlabel(get(tmp_plot(1),'Parent'),'Cluster number K');
  388. ylabel(get(tmp_plot(1),'Parent'),'Stability');
  389. xlim(get(tmp_plot(1),'Parent'),[2-0.6,K_range(end)+0.6]);
  390. ylim(get(tmp_plot(1),'Parent'),[0,1]);
  391. set(get(tmp_plot(1),'Parent'),'Box','off');
  392. custom_cm = cbrewer('seq','Reds',25);
  393. colormap(CCPlot,custom_cm(6:25,:));
  395. saveas(gcf,fullfile(SavePath,fancyName,'Stability.jpg'));
  396. saveas(gcf,fullfile(SavePath,fancyName,'Stability.fig'));
  397. disp(' ');
  398. disp('-----------------------------------------------------');
  399. disp(['Figure stability measures has been saved in ' SavePath]);
  400. disp(' ');
  402. disp('Please check stability measures');
  404. % You should fill this with the actual value
  405. %K_opt = 3;
  406. prompt = 'What is the optimal cluster size? ';
  407. K_opt = input(prompt);
  408. close;
  410. %% 5. Clustering into CAPs
  412. %[CAP,Disp,Std_Clusters,idx1,CorrDist,sfrac] = Run_Clustering(cell2mat(Xon1),...
  413. % K_opt,mask{1},brain_info{1},Pp,Pn,n_rep,idx_sep_seeds1,SeedType);
  414. % CAP_raw = CAP;
  415. [CAP,Disp,Std_Clusters,idx1,CorrDist,sfrac] = Run_Clustering_PCA(W,...
  416. K_opt,mask{1},brain_info{1},Pp,Pn,n_rep,idx_sep_seeds1,SeedType);
  418. %Plot
  419. % Computation of the similarity
  420. SimMat = corr(CAP',CAP');
  421. SimMat(isnan(SimMat))=0;
  423. % Similarity Plot
  424. imagesc(SimMat);
  426. tmp_cb2 = cbrewer('div','RdBu',1000,'spline'); % had to ad 'spline', otherwise error message (SH)
  427. tmp_cb2(tmp_cb2 < 0) = 0;
  429. colormap(flipud(tmp_cb2));
  431. % Arbitrary setting of probability scale
  432. caxis([-1 1]);
  433. axis('square','on');
  434. axis('on');
  436. % Add colorbar to the figure
  437. colorbar; % This adds the color scale
  439. saveas(gcf,fullfile(SavePath,fancyName,'Similarity.jpg'));
  440. saveas(gcf,fullfile(SavePath,fancyName,'Similarity.fig'));
  441. disp(' ');
  442. disp('-----------------------------------------------------');
  443. disp(['Figure Similarity matrix has been saved in ' SavePath]);
  444. disp(' ');
  446. close;
  449. %% 5.1. PCA Reconstruction
  451. CAP_original = (U*CAP') + mu;
  452. CAP = CAP_original';
  454. %% 6. Assignment of the frames from population 2
  456. % Parameter that governs the stringency of assignment: if Ap = 5%, we
  457. % assign a frame to a CAP if spatial correlation exceeds the 5th percentile
  458. % of the distribution of spatial correlations between the CAP, and its
  459. % constituting frames
  460. Ap = 5;
  462. idx2 = CAP_AssignFrames(CAP,cell2mat(Xon2),CorrDist,Ap)';
  466. %% 7. Computing metrics
  468. % The TR of your data in seconds
  469. TR = 1.3;
  471. [ExpressionMap1,Counts1,Entries1,Avg_Duration1,Duration1,TransitionProbabilities1,...
  472. From_Baseline1,To_Baseline1,Baseline_resilience1,Resilience1,Betweenness1,...
  473. InDegree1,OutDegree1,SubjectEntries1] = Compute_Metrics_simpler(idx1,...
  474. Indices1.kept.active,Indices1.scrubbedandactive,K_opt,TR);
  476. [ExpressionMap2,Counts2,Entries2,Avg_Duration2,Duration2,TransitionProbabilities2,...
  477. From_Baseline2,To_Baseline2,Baseline_resilience2,Resilience2,Betweenness2,...
  478. InDegree2,OutDegree2,SubjectEntries2] = Compute_Metrics_simpler(idx2,...
  479. Indices2.kept.active,Indices2.scrubbedandactive,K_opt,TR);
  481. %% Transition Probability (Visualization and Statistics)
  483. % 1. Plot metrics Transition matrix for all HC
  485. tmp_toplot = squeeze(mean(TransitionProbabilities1,3));
  486. tmp_toplot = tmp_toplot(3:end-1,3:end-1);
  488. % Make graph visible and plotting
  489. TMGraph=imagesc(tmp_toplot);
  490. tmp_cb = cbrewer('seq','Greys',1000);
  491. colormap(flipud(tmp_cb));
  493. clear tmp_toplot
  495. % Arbitrary setting of probability scale from 0 to 0.03
  496. caxis([0 0.03]);
  497. axis('square','on');
  498. axis('on');
  500. % Add colorbar to the figure
  501. colorbar; % This adds the color scale
  503. saveas(gcf,fullfile(SavePath,fancyName,'TransitionProbabilityMatrix_HC.jpg'));
  504. saveas(gcf,fullfile(SavePath,fancyName,'TransitionProbabilityMatrix_HC.fig'));
  505. close;
  507. % 2. Plot metrics Transition matrix for all CPP
  508. tmp_toplot = squeeze(mean(TransitionProbabilities2,3));
  509. tmp_toplot = tmp_toplot(3:end-1,3:end-1);
  511. % Make graph visible and plotting
  512. TMGraph=imagesc(tmp_toplot);
  513. tmp_cb = cbrewer('seq','Greys',1000);
  514. colormap(flipud(tmp_cb));
  516. clear tmp_toplot
  518. % Arbitrary setting of probability scale from 0 to 0.03
  519. caxis([0 0.03]);
  520. axis('square','on');
  521. axis('on');
  523. % Add colorbar to the figure
  524. colorbar; % This adds the color scale
  526. saveas(gcf,fullfile(SavePath,fancyName,'TransitionProbabilityMatrix_CPP.jpg'));
  527. saveas(gcf,fullfile(SavePath,fancyName,'TransitionProbabilityMatrix_CPP.fig'));
  528. close;
  530. disp(' ');
  531. disp('-----------------------------------------------------');
  532. disp(['Figure Transition Probability Matrix has been saved in ' SavePath]);
  533. disp(' ');
  535. %% Visualization of CAP and Frame Dynamics
  537. % 1. Dynamic state plotting
  539. % Makes the graph visible
  540. % Concatenates information from the different datasets
  541. tmp_toplot = [];
  542. ExpressionMap{1}=ExpressionMap1;
  543. ExpressionMap{2}=ExpressionMap2;
  544. for i = 1:n_datasets
  545. tmp_toplot = [tmp_toplot; ExpressionMap{i}; 0*ones(5,size(TC{1},1))];
  546. end
  547. tmp_toplot = tmp_toplot(1:end-5,:); % combines HC and CPPs and but between the dataframes 5 rows with 0, that they can be visually seperated.
  549. % Define custom colormap: [black for scrubbed (-1), white for non-retained (0), then CAP colors]
  550. custom_cm = cbrewer('qual','Set1',K_opt+1);
  551. custom_cm = [0.05,0.05,0.05;1,1,1;custom_cm];
  553. %Plot it
  554. imagesc(tmp_toplot);
  555. colormap((custom_cm));
  556. xlabel('Time [s]');
  557. ylabel('Subjects [-]');
  558. caxis([-1,K_opt+1]); %Adjust color axis to include scrubbed, non-retained, CAPs, and unassigned
  560. clear tmp_toplot
  562. % Generate line objects for the legend and assign colors
  563. L = gobjects(1, K_opt + 3); % Preallocate for legend lines (CAPs + Scrubbed + Non-retained + Unassigned)
  565. % Loop to create legend lines for CAPs
  566. for cap_idx = 1:K_opt
  567. % Create line for CAPs
  568. L(cap_idx) = line(nan(4,1), nan(4,1), 'LineWidth', 2, 'Color', custom_cm(cap_idx + 2, :)); % +2 skips black and white
  569. end
  571. % Add lines for 'Scrubbed' (-1), 'Non-retained' (0), and 'Unassigned'
  572. L(K_opt + 1) = line(nan(4,1), nan(4,1), 'LineWidth', 2, 'Color', custom_cm(1, :)); % Black for scrubbed (-1)
  573. L(K_opt + 2) = line(nan(4,1), nan(4,1), 'LineWidth', 2, 'Color', custom_cm(2, :)); % White for non-retained (0)
  574. L(K_opt + 3) = line(nan(4,1), nan(4,1), 'LineWidth', 2, 'Color', custom_cm(end, :)); % Last color for unassigned
  576. % Create the legend dynamically based on K_opt
  577. legend_entries = cell(1, K_opt + 3); % Preallocate cell array for legend entries (CAPs + 3 categories)
  578. for cap_idx = 1:K_opt
  579. legend_entries{cap_idx} = sprintf('CAP%d', cap_idx); % Generate CAP labels
  580. end
  581. legend_entries{K_opt + 1} = 'Scrubbed'; % Black for scrubbed (-1)
  582. legend_entries{K_opt + 2} = 'Non-retained'; % White for non-retained (0)
  583. legend_entries{K_opt + 3} = 'Unassigned'; % Color from colormap for unassigned (K_opt + 1)
  585. % Create the legend
  586. lgd = legend(L, legend_entries{:}, 'Location', 'best');
  587. lgd.FontSize = 10;
  588. lgd.Position(1) = 0.65;
  589. lgd.Position(2) = 0.5;
  591. % Ensure the legend is shown
  592. legend show;
  594. saveas(gcf,fullfile(SavePath,fancyName,'DynamicStates.jpg'));
  595. saveas(gcf,fullfile(SavePath,fancyName,'DynamicStates.fig'));
  596. close;
  598. disp(' ');
  599. disp('-----------------------------------------------------');
  600. disp(['Figure Dynamic States has been saved in ' SavePath]);
  601. disp(' ');
  603. % 2. Averaged by CAP HC
  605. % Define custom colormap using cbrewer (for example, Set1 colormap)
  606. custom_cm = cbrewer('qual', 'Set1', K_opt + 1); % Add 1 for an additional color if needed
  608. % Find unique CAPs in the matrix
  609. CAPs = unique(ExpressionMap1);
  611. % Preallocate a matrix to store the normalized sums for each column
  612. normalizedSums = zeros(length(CAPs), size(ExpressionMap1, 2)); % CAPs x Columns
  614. % Loop through each column and calculate the normalized counts
  615. for col = 1:size(ExpressionMap1, 2)
  616. for i = 1:length(CAPs)
  617. CAP_value = CAPs(i);
  618. count = sum(ExpressionMap1(:, col) == CAP_value); % Count occurrences of CAP in this column
  619. normalizedSums(i, col) = count / size(ExpressionMap1, 1); % Normalize by dividing by the number of rows
  620. end
  621. end
  623. % Exclude CAP0 and CAP-1 (if they exist)
  624. excludedCAPs = [-1, 0];
  625. validCAPs = ~ismember(CAPs, excludedCAPs); % Find CAPs that are not -1 or 0
  627. % Filter out rows and CAPs
  628. filteredNormalizedSums = normalizedSums(validCAPs, :);
  629. filteredCAPs = CAPs(validCAPs); % Corresponding CAP values
  631. % Visualization: Line plot (Only lines, no markers)
  632. figure;
  633. hold on;
  635. % Use custom colormap for the lines
  636. for i = 1:length(filteredCAPs)
  637. plot(1:size(ExpressionMap1, 2), filteredNormalizedSums(i, :), '-', 'LineWidth', 1.5, ...
  638. 'Color', custom_cm(i + 1, :), 'DisplayName', sprintf('CAP%d', filteredCAPs(i))); % +1 to skip first color if needed
  639. end
  640. hold off;
  642. % Customize the plot
  643. xlabel('Frames');
  644. ylabel('Average CAP expression');
  645. title('HC');
  646. legend('Location', 'best');
  647. grid on;
  648. set(gca, 'FontSize', 12);
  650. saveas(gcf,fullfile(SavePath,fancyName,'AverageCAPExpressionHC.jpg'));
  651. saveas(gcf,fullfile(SavePath,fancyName,'AverageCAPExpressionHC.fig'));
  652. close;
  654. disp(' ');
  655. disp('-----------------------------------------------------');
  656. disp(['Figure Average CAP eypression HC has been saved in ' SavePath]);
  657. disp(' ');
  660. % 3. Averaged by CAP CPP
  662. custom_cm2 = cbrewer('qual', 'Set1', K_opt + 1); % Add 1 for an additional color if needed
  663. % Find unique CAPs in the matrix
  664. CAPs2 = unique(ExpressionMap2);
  666. % Preallocate a matrix to store the normalized sums for each column
  667. normalizedSums2 = zeros(length(CAPs2), size(ExpressionMap2, 2)); % CAPs x Columns
  669. % Loop through each column and calculate the normalized counts
  670. for col = 1:size(ExpressionMap2, 2)
  671. for i = 1:length(CAPs2)
  672. CAP_value2 = CAPs2(i);
  673. count2 = sum(ExpressionMap2(:, col) == CAP_value2); % Count occurrences of CAP in this column
  674. normalizedSums2(i, col) = count2 / size(ExpressionMap2, 1); % Normalize by dividing by the number of rows
  675. end
  676. end
  678. % Exclude CAP0 and CAP-1 (if they exist)
  679. excludedCAPs2 = [-1, 0, K_opt+1];
  680. validCAPs2 = ~ismember(CAPs2, excludedCAPs2); % Find CAPs that are not -1 or 0
  682. % Filter out rows and CAPs
  683. filteredNormalizedSums2 = normalizedSums2(validCAPs2, :);
  684. filteredCAPs2 = CAPs2(validCAPs2); % Corresponding CAP values
  686. % Visualization: Line plot (Only lines, no markers)
  687. figure;
  688. hold on;
  690. % Use custom colormap for the lines
  691. for i = 1:length(filteredCAPs2)
  692. plot(1:size(ExpressionMap2, 2), filteredNormalizedSums2(i, :), '-', 'LineWidth', 1.5, ...
  693. 'Color', custom_cm2(i + 1, :), 'DisplayName', sprintf('CAP%d', filteredCAPs2(i))); % +1 to skip first color if needed
  694. end
  695. hold off;
  697. % Customize the plot
  698. xlabel('Frames');
  699. ylabel('Average CAP expression');
  700. title('CPP');
  701. legend('Location', 'best');
  702. grid on;
  703. set(gca, 'FontSize', 12);
  705. saveas(gcf,fullfile(SavePath,fancyName,'AverageCAPExpressionCPP.jpg'));
  706. saveas(gcf,fullfile(SavePath,fancyName,'AverageCAPExpressionCPP.fig'));
  707. close;
  709. disp(' ');
  710. disp('-----------------------------------------------------');
  711. disp(['Figure Average CAP eypression CPP has been saved in ' SavePath]);
  712. disp(' ');
  714. %% Quality Check 2 - Visualization
  716. % Number of subjects in each group
  717. nHC = size(ExpressionMap{1}, 1); % Number of HC subjects
  718. nCPP = size(ExpressionMap{2}, 1); % Number of CPP subjects
  720. % Initialize vectors for subject-wise ratios
  721. ratio_HC = zeros(nHC, 1);
  722. ratio_CPP = zeros(nCPP, 1);
  725. % 1. Ration unnassigned frames to assigned frames (CAPs) in Population
  726. % Compute ratios for CPP subjects
  727. for i = 1:nCPP
  728. assigned_frames = sum(ismember(ExpressionMap{2}(i, :), 1:K_opt)); % Assigned frames
  729. unassigned_frames = sum(ExpressionMap{2}(i, :) == (K_opt+1)); % Unassigned frames
  730. ratio_CPP(i) = (unassigned_frames / assigned_frames)
  731. end
  733. % Perform a one-sample t-test to check if the mean differs from zero
  734. [h, p_value1, ci, stats] = ttest(ratio_CPP);
  736. % Create figure
  737. figure;
  738. hold on;
  740. % Set X-axis position to center
  741. boxplot(ratio_CPP, 'Positions', 1.5, 'Labels', {'CPP'}, 'Whisker', 1.5, 'Colors', 'r', 'Symbol', 'ro', 'Widths', 0.6);
  742. set(findobj(gca,'Type','Line','Tag','Median'), 'Color', 'k', 'LineWidth', 2); % Highlight median
  744. % Scatter individual points (improved jittering for better visibility)
  745. x_CPP = 1.5 + 0.08 * randn(nCPP, 1); % Adjust X position for centering
  746. scatter(x_CPP, ratio_CPP, 90, 'r', 'filled', 'MarkerFaceAlpha', 0.6, 'MarkerEdgeColor', 'k'); % Improved dot styling
  748. % Formatting improvements
  749. ylabel('Ratio of Unassigned to Assigned Frames', 'FontSize', 14, 'FontWeight', 'bold');
  750. title('Subject-Level Ratio of Unassigned to Assigned Frames (CPP)', 'FontSize', 16, 'FontWeight', 'bold');
  751. ylim([0 1]); % Set Y-axis limit for better readability
  753. % Center the X-axis
  754. xlim([1 2]); % Set range to center the boxplot
  756. % Annotate p-value on the plot
  757. text(1.5, 9, sprintf('p = %.4f', p_value1), 'HorizontalAlignment', 'center', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'blue');
  759. % Remove grid for cleaner look
  760. grid off;
  762. % Save figure
  763. saveas(gcf, fullfile(SavePath, fancyName, 'Unassigned_to_assigned_frames_CPP.jpg'));
  764. saveas(gcf, fullfile(SavePath, fancyName, 'Unassigned_to_assigned_frames_CPP.fig'));
  765. close;
  767. % 2. Ration scrubbed (but active) frames to assigned frames (CAPs)
  769. % Compute ratios for HC subjects
  770. for i = 1:nHC
  771. assigned_frames = sum(ismember(ExpressionMap{1}(i, :), 1:K_opt)); % Assigned frames
  772. scrubbed_frames = sum(ExpressionMap{1}(i, :) == -1); % Scrubbed frames
  773. ratio_HC(i) = (scrubbed_frames / assigned_frames) * 100; % Convert to percentage
  774. end
  776. % Compute ratios for CPP subjects
  777. for i = 1:nCPP
  778. assigned_frames = sum(ismember(ExpressionMap{2}(i, :), 1:K_opt)); % Assigned frames
  779. scrubbed_frames = sum(ExpressionMap{2}(i, :) == -1); % Scrubbed frames
  780. ratio_CPP(i) = (scrubbed_frames / assigned_frames) * 100; % Convert to percentage
  781. end
  783. % Combine data and create group labels
  784. all_ratios = [ratio_HC; ratio_CPP]; % Merge HC and CPP ratios
  785. group_labels = [ones(nHC, 1); 2 * ones(nCPP, 1)]; % 1 for HC, 2 for CPP
  787. % Non-parametric test (Mann-Whitney U test)
  788. [p_value2, h, stats] = ranksum(ratio_HC, ratio_CPP);
  790. % Create figure
  791. figure;
  792. hold on;
  794. % Corrected boxplot without manual 'positions' parameter
  795. boxplot(all_ratios, group_labels, 'Labels', {'HC', 'CPP'}, ...
  796. 'Whisker', 1.5, 'Colors', 'rb', 'Symbol', 'ro', 'Widths', 0.6);
  797. set(findobj(gca,'Type','Line','Tag','Median'), 'Color', 'k', 'LineWidth', 2); % Highlight median
  799. % Scatter individual points (jittered for better visibility)
  800. x_HC = 1 + 0.08 * randn(nHC, 1); % Jitter HC points
  801. x_CPP = 2 + 0.08 * randn(nCPP, 1); % Jitter CPP points
  802. scatter(x_HC, ratio_HC, 90, 'b', 'filled', 'MarkerFaceAlpha', 0.6, 'MarkerEdgeColor', 'k'); % HC points
  803. scatter(x_CPP, ratio_CPP, 90, 'r', 'filled', 'MarkerFaceAlpha', 0.6, 'MarkerEdgeColor', 'k'); % CPP points
  805. % Formatting improvements
  806. ylabel('Ratio of Scrubbed to Assigned Frames (%)', 'FontSize', 14, 'FontWeight', 'bold');
  807. title('Scrubbed (but active) to Assigned Frames (HC vs. CPP)', 'FontSize', 16, 'FontWeight', 'bold');
  808. ylim([0 100]); % Set Y-axis limit for better readability
  809. yticks(0:10:100); % Set ticks at intervals of 10%
  811. % Annotate p-value on the plot
  812. text(1.5, 90, sprintf('p = %.4f', p_value2), 'HorizontalAlignment', 'center', ...
  813. 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'blue');
  815. % Remove grid for cleaner look
  816. grid off;
  818. % Save figure
  819. saveas(gcf, fullfile(SavePath, fancyName, 'Scrubbed_to_assigned_frames_HC_CPP.jpg'));
  820. saveas(gcf, fullfile(SavePath, fancyName, 'Scrubbed_to_assigned_frames_HC_CPP.fig'));
  821. close;
  824. %% 7. Save Metrics
  825. % General information on the project
  826. OverallInfo.ProjectTitle = fancyName;
  827. OverallInfo.NumberTimePoints = size(TC{1},1);%handles.SubjSize.TP;
  828. OverallInfo.NumberVoxels = size(TC{1},2); %handles.SubjSize.VOX;
  829. OverallInfo.NumberSubjects = N_Subj; %handles.n_subjects;
  830. OverallInfo.TR = TR; %handles.TR;
  832. Parameters.Inputs.DataHeader = brain_info;%handles.brain_info;
  833. Parameters.Inputs.Mask = mask;
  834. Parameters.Inputs.Seeds = Seed;
  835. Parameters.SpatioTemporalSelection.IsSeedFree = 0; %handles.is_seed_free;
  836. Parameters.SpatioTemporalSelection.NumberSeeds = size(seedName,2); %handles.n_seed;
  837. Parameters.SpatioTemporalSelection.TypeEventRetainedPerSeed = SignMatrix;
  838. Parameters.SpatioTemporalSelection.SeedType = SeedType;
  839. Parameters.SpatioTemporalSelection.MotionThreshold = Tmot;
  840. Parameters.SpatioTemporalSelection.SelectionMode = SelMode;
  841. Parameters.SpatioTemporalSelection.FrameSelectionParameter = T;
  842. Parameters.KMeansClustering.IsConsensusRun = 1;% handles.is_consensus_clustering;
  843. Parameters.KMeansClustering.MaxClusterNumber = 6; %handles.Kmax; change this to the respective Kmax.
  844. Parameters.KMeansClustering.PercentageDataPerFold = Pcc; %handles.PCC;
  845. Parameters.KMeansClustering.NumberRepetitions = n_rep;
  846. Parameters.KMeansClustering.NumberClusters = K_opt; %handles.K;
  847. Parameters.KMeansClustering.PercentagePositiveValuedVoxelsClustered = Pp;
  848. Parameters.KMeansClustering.PercentageNegativeValuedVoxelsClustered = Pn;
  850. Outputs.SpatioTemporalSelection.RetainedFramesPerSeed{1} = idx_sep_seeds1;
  851. Outputs.SpatioTemporalSelection.RetainedFramesPerSeed{2} = idx_sep_seeds2;
  852. Outputs.SpatioTemporalSelection.Indices{1}=Indices1;
  853. Outputs.SpatioTemporalSelection.Indices{2}=Indices2; %hier war vorher eine 1
  854. Outputs.SpatioTemporalSelection.PercentageRetainedFrames = RetainedPercentage;
  855. Outputs.SpatioTemporalSelection.AverageCorrelationMap = AvgSeedMap;
  856. Outputs.KMeansClustering.ConsensusQuality = PAC; %handles.ConsensusQuality; CHECK IF THAT'S TRUE
  857. Outputs.KMeansClustering.CoActivationPatternsDispersion = Disp;
  858. Outputs.KMeansClustering.CoActivationPatterns = CAP;
  859. Outputs.KMeansClustering.CoActivationPatternsZScored = CAP_Zscore(CAP);
  860. Outputs.KMeansClustering.CoActivationPatternsSTD = Std_Clusters; %handles.STDCAP;
  861. Outputs.KMeansClustering.AssignmentsToCAPs{1} = idx1;
  862. Outputs.KMeansClustering.AssignmentsToCAPs{2} = idx2;
  863. Outputs.Metrics.CAPExpressionIndices{1} = ExpressionMap1;
  864. Outputs.Metrics.CAPExpressionIndices{2} = ExpressionMap2;
  865. Outputs.Metrics.Occurrences{1} = Counts1;
  866. Outputs.Metrics.Occurrences{2} = Counts2;
  867. Outputs.Metrics.NumberEntries{1} = Entries1;
  868. Outputs.Metrics.NumberEntries{2} = Entries2;
  869. Outputs.Metrics.AverageExpressionDuration{1} = Avg_Duration1;
  870. Outputs.Metrics.AverageExpressionDuration{2} = Avg_Duration2;
  871. Outputs.Metrics.AllExpressionDurations{1} = Duration1;
  872. Outputs.Metrics.AllExpressionDurations{2} = Duration2;
  873. Outputs.Metrics.TransitionProbabilities{1} = TransitionProbabilities1;
  874. Outputs.Metrics.TransitionProbabilities{2} = TransitionProbabilities2;
  875. Outputs.Metrics.FractionCAPFramesPerSeedCombination = sfrac;
  876. Outputs.Metrics.CAPEntriesFromBaseline{1} = From_Baseline1;
  877. Outputs.Metrics.CAPEntriesFromBaseline{2} = From_Baseline2;
  878. Outputs.Metrics.CAPExitsToBaseline{1} = To_Baseline1;
  879. Outputs.Metrics.CAPExitsToBaseline{2} = To_Baseline2;
  880. Outputs.Metrics.CAPResilience{1} = Resilience1;
  881. Outputs.Metrics.CAPResilience{2} = Resilience2;
  882. Outputs.Metrics.BaselineResilience{1} = Baseline_resilience1;
  883. Outputs.Metrics.BaselineResilience{2} = Baseline_resilience2;
  884. Outputs.Metrics.BetweennessCentrality{1} = Betweenness1;
  885. Outputs.Metrics.BetweennessCentrality{2} = Betweenness2;
  886. Outputs.Metrics.CAPInDegree{1} = InDegree1;
  887. Outputs.Metrics.CAPInDegree{2} = InDegree2;
  888. Outputs.Metrics.CAPOutDegree{1} = OutDegree1;
  889. Outputs.Metrics.CAPOutDegree{2} = OutDegree2;
  890. Outputs.Metrics.SubjectCounts{1} = SubjectEntries1;
  891. Outputs.Metrics.SubjectCounts{2} = SubjectEntries2;
  893. HeavyOutputs.SpatioTemporalSelection.ClusteredFrames{1} = Xon1; % Are this 3 lines necessary to be saved?
  894. HeavyOutputs.SpatioTemporalSelection.ClusteredFrames{2} = Xon2;%Xonp?
  895. HeavyOutputs.KMeansClustering.Consensus = Consensus;
  897. % Saves NIFTI files storing the CAPs in MNI space
  898. ReferencePopulation = 1;
  899. CAPToNIFTI(CAP,...
  900. mask{ReferencePopulation},brain_info{ReferencePopulation},...
  901. fullfile(SavePath,fancyName),['CAP_NIFTI_',fancyName]);
  903. CAPToNIFTI(CAP_Zscore(CAP),...
  904. mask{ReferencePopulation},brain_info{ReferencePopulation},...
  905. fullfile(SavePath,fancyName),['CAP_NIFTI_ZScored_',fancyName]);
  907. % Saves the different variables from the program
  908. save(fullfile(SavePath,fancyName),'OverallInfo','Parameters','Outputs','HeavyOutputs','brain','-v7.3');
  909. %save(fullfile(SavePath,fancyName),'OverallInfo','Parameters','Outputs','brain','-v7.3');
  912. disp('-----------------------------------------------------');
  913. disp('CAPs computed and saved successfully!');
  914. disp(' ');
  915. disp('Results are saved in: ');
  916. disp([SavePath]);

Script_twopop_PCA_HC_ref_SH.m at commit 1240c94, no license · at the source

Overview

Authors: Salome Häuselmann1,2,3, Anna Wyss1,3, Samantha Weber1,4,5, Nicolas Gninenko1,6, Cristina Concetti6, Eliane Müller1,2,6, Rupert Bruckmaier7, Josef Gross7, Nina Bischoff1, Chantal Berna8, Martin grosse Holtforth1,9, Selma Aybek6
ORCID iDs: Salome Häuselmann, Samantha Weber, Nicolas Gninenko, Rupert Bruckmaier, Josef Gross, Selma Aybek
  1. Psychosomatic Medicine, Department of Neurology, Inselspital, Bern University Hospital, University of Bern, 3010 Bern, Switzerland
  2. Graduate School of Cellular and Biomedical Sciences (GCB), University of Bern, 3012 Bern, Switzerland
  3. Translational Imaging Center (TIC), Swiss Institute for Translational and Entrepreneurial Medicine, 3010 Bern, Switzerland
  4. Department of Adult Psychiatry and Psychotherapy, University Hospital of Psychiatry Zurich, University of Zurich, 8032 Zurich, Switzerland
  5. Faculty of Medicine, University of Zurich, 8032 Zurich, Switzerland
  6. Department of Neurology, Faculty of Science and Medicine, University of Fribourg, 1700 Fribourg, Switzerland
  7. Veterinary Physiology, Vetsuisse Faculty, University of Bern, 3012 Bern, Switzerland
  8. Center for Integrative and Complementary Medicine, Department of Anesthesiology, Lausanne University Hospital, 1011 Lausanne, Switzerland
  9. Institute of Psychology, University of Bern, 3012 Bern, Switzerland
Institutions: University of Bern (Switzerland); University Hospital of Bern (Switzerland); University of Zurich (Switzerland); University Hospital Zurich (Switzerland); Psychiatrische Universitätsklinik Zürich (Switzerland); University of Fribourg (Switzerland); Centre Hospitalier Universitaire Vaudois (Switzerland); University of Lausanne (Switzerland)
Journal: Brain communications, volume 8, issue 2, article fcag121
Dates: received 25 October 2025; accepted 2 April 2026; published online 3 April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1093/braincomms/fcag121 · PMID 41978789 · PMCID PMC13070617 · OpenAlex W7149391556
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: fMRI (modality), human (organism), pain (population), clinical / translational (subfield)
Methods: Connectivity, Spectral & time-frequency, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, fMRI & imaging
Keywords: chronic primary pain, functional MRI, co-activation patterns, brain dynamics, stress response
Topic: Pain Mechanisms and Treatments (Physiology, Medicine), according to OpenAlex
Funding: Swiss National Science Foundation (176985, PP00P3_176985)
Citations: cited by 1 paper (Europe PMC); 74 references in the paper

Abstract

Chronic primary pain occurs without an identifiable causal disease and is marked by persistent pain, emotional distress and functional disability. The anterior insular cortex, involved in salience processing and integration of sensory, emotional and cognitive aspects of pain, has been implicated in neural processes linking pain and stress responses. This study investigates whether specific brain state dynamics, using the anterior insula as a seed region, are associated with chronic primary pain and examines their associations with pain- and stress-related measures. Resting-state functional MRI, stress biomarkers (cortisol and alpha-amylase), a pain sensitivity test, as well as subjective measures of stress and pain were collected from patients with chronic primary pain (N = 30) and healthy controls (N = 30). Co-activation pattern analysis was used to identify brain states with the anterior insula as the seed region and to assess group differences in the temporal characteristics of these brain states. Partial least squares analysis was applied to investigate multivariate associations between specific temporal brain state characteristics and pain- and stress-related measures. Three anterior insula–seeded co-activation patterns (brain states) were identified in healthy controls. In the first co-activation pattern, the anterior insula co-activated with the default mode network; in the second, with the salience-somatomotor network; and in the third, with the visual network. Chronic primary pain patients and healthy controls differed significantly in temporal brain state characteristics, namely in the relative number of entries into co-activation pattern one (pFDR = 0.002) and two (pFDR = 0.022), and the relative occurrence of co-activation patterns one (pFDR = 0.022) and two (pFDR = 0.022). Furthermore, in chronic primary pain patients, perceived stress scores and cortisol were significantly associated with these specific temporal brain state characteristics (P = 0.002), whereas no associations were found with pain-related measures. Together, these findings suggest that in chronic primary pain, reduced coupling of the anterior insula with the default mode network and increased coupling of the anterior insula with salience-related networks are associated with psychophysiological stress markers. These brain state dynamics may potentially represent a neural correlate of altered stress processing in chronic primary pain.

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.

FND-ResearchGroup/CAP_in_CPP

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 1240c94a8fdfea5a597d5aab584649140b820c0b, 30 March 2026
Languages: R (3), MATLAB (2), Jupyter (1)
Size: 10 files, 6 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, 2 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: car (3 files), cowplot (3 files), ggplot2 (3 files), ggpubr (3 files), lme4 (3 files), lmerTest (3 files), multcomp (3 files), reshape2 (3 files), rstatix (3 files), tidyverse (3 files), easystats (1 file), Statistics and Machine Learning Toolbox (1 file), NiBabel (1 file), NumPy (1 file), pandas (1 file), Plotly (1 file), psych (1 file), SciPy (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
7 files
At the source: github.com/FND-ResearchGroup/CAP_in_CPP

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

Tracing map

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

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 6 scripts, each with its path and the digest of its content;
  • 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

The data are not publicly available but can be shared upon request. The TbCAPs toolbox is publicly available at https://github.com/MIPLabCH/TbCAPs, and the PLS Toolbox can be accessed at https://github.com/FND-ResearchGroup/myPLS_SL. Additional codes used for data analysis are available at https://github.com/FND-ResearchGroup/CAP_in_CPP.

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, 12 authors, 5 keywords, 1 funder, 69 references.

Cite

This paper

Häuselmann, S., Wyss, A., Weber, S., Gninenko, N., Concetti, C., Müller, E., Bruckmaier, R., Gross, J., Bischoff, N., Berna, C., grosse Holtforth, M., & Aybek, S. (2026). Anterior insular co-activation patterns associated with stress markers in chronic primary pain. Brain communications, 8(2), fcag121. https://doi.org/10.1093/braincomms/fcag121

BibTeX

@article{hauselmann2026anterior,
author = {Häuselmann, Salome and Wyss, Anna and Weber, Samantha and Gninenko, Nicolas and Concetti, Cristina and Müller, Eliane and Bruckmaier, Rupert and Gross, Josef and Bischoff, Nina and Berna, Chantal and grosse Holtforth, Martin and Aybek, Selma},
title = {{Anterior insular co-activation patterns associated with stress markers in chronic primary pain}},
journal = {Brain communications},
year = {2026},
month = apr,
volume = {8},
number = {2},
pages = {fcag121},
publisher = {Oxford University Press},
issn = {2632-1297},
doi = {10.1093/braincomms/fcag121},
url = {https://doi.org/10.1093/braincomms/fcag121},
pmid = {41978789},
pmcid = {PMC13070617}
}

RIS

TY - JOUR
AU - Häuselmann, Salome
AU - Wyss, Anna
AU - Weber, Samantha
AU - Gninenko, Nicolas
AU - Concetti, Cristina
AU - Müller, Eliane
AU - Bruckmaier, Rupert
AU - Gross, Josef
AU - Bischoff, Nina
AU - Berna, Chantal
AU - grosse Holtforth, Martin
AU - Aybek, Selma
TI - Anterior insular co-activation patterns associated with stress markers in chronic primary pain
T2 - Brain communications
J2 - Brain Commun
PY - 2026
DA - 2026/04/03
VL - 8
IS - 2
SP - fcag121
SN - 2632-1297
PB - Oxford University Press
DO - 10.1093/braincomms/fcag121
UR - https://doi.org/10.1093/braincomms/fcag121
LA - en
ER -

CSL-JSON

{
"id": "10.1093/braincomms/fcag121",
"type": "article-journal",
"title": "Anterior insular co-activation patterns associated with stress markers in chronic primary pain",
"container-title": "Brain communications",
"author": [
{
"family": "Häuselmann",
"given": "Salome"
},
{
"family": "Wyss",
"given": "Anna"
},
{
"family": "Weber",
"given": "Samantha"
},
{
"family": "Gninenko",
"given": "Nicolas"
},
{
"family": "Concetti",
"given": "Cristina"
},
{
"family": "Müller",
"given": "Eliane"
},
{
"family": "Bruckmaier",
"given": "Rupert"
},
{
"family": "Gross",
"given": "Josef"
},
{
"family": "Bischoff",
"given": "Nina"
},
{
"family": "Berna",
"given": "Chantal"
},
{
"family": "grosse Holtforth",
"given": "Martin"
},
{
"family": "Aybek",
"given": "Selma"
}
],
"container-title-short": "Brain Commun",
"volume": "8",
"issue": "2",
"page": "fcag121",
"DOI": "10.1093/braincomms/fcag121",
"PMID": "41978789",
"PMCID": "PMC13070617",
"ISSN": "2632-1297",
"publisher": "Oxford University Press",
"URL": "https://doi.org/10.1093/braincomms/fcag121",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
3
]
]
}
}

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.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: psych, rstatix, car, 11 other tools, fMRI, 1 reference
[2] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: multcomp, rstatix, car, 10 other tools
[3] doi:10.1162/imag.a.1245 [code]
Towards precision EEG connectomics: Evaluating the benefits of dense sampling.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: psych, rstatix, lmerTest, 11 other tools
[4] doi:10.1093/braincomms/fcag146 [code]
Convergent structural brain alterations in chronic pain: a multi-metric individual participant data meta-analysis.
Journal: Brain communications
In common: multcomp, psych, car, 5 other tools, pain, 2 references
[5] doi:10.1523/eneuro.0076-26.2026 [code]
Exogenously Driven Neural Reactivation of Spatially Matching Visual Working-Memory Contents.
Journal: eNeuro
In common: multcomp, rstatix, car, 9 other tools
[6] doi:10.1093/cercor/bhag077 [code]
The longitudinal development of intrinsic timescales in infancy and their relation to alpha brain rhythm.
Journal: Cerebral cortex (New York, N.Y. : 1991)
In common: multcomp, psych, easystats, 9 other tools
[7] doi:10.1016/j.neuroimage.2026.122115 [code]
Midfrontal theta power relates to response speeding following frustrative nonreward.
Journal: NeuroImage
In common: psych, rstatix, easystats, 8 other tools
[8] doi:10.1093/brain/awag080 [code]
Early glymphatic failure in AppNL-F knock-in mice is linked to parenchymal border macrophages loss.
Journal: Brain : a journal of neurology
In common: multcomp, psych, rstatix, 7 other tools
[9] doi:10.1162/imag.a.1321 [code]
Phase similarity between similar objects indicates representational merging across retrieval training but not sleep.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: rstatix, easystats, car, 9 other tools
[10] doi:10.1038/s41380-026-03694-1 [code]
Targeting cortico-striatal-amygdalar networks via theta-band frontoparietal synchronization in opioid use disorder: a randomized tACS-fMRI Trial.
Journal: Molecular psychiatry
In common: multcomp, lmerTest, lme4, 7 other tools, pain, fMRI, clinical / translational

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.