OSCR

Feasibility of precision functional mapping in youth multi-echo fMRI data.

Code ↔ Paper

6 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 6 matches
  1. [1] § Methods › Individualized networks ↔ ToShareCode.zip/ToShareCode/FC_Reliability/makeFC.m, lines 165–227 · score 1.00 · Somato Cognitive Action, Cingulo opercular, Somatomotor Foot, Medial Parietal, Visual V5, Default Anterolateral
  2. [2] § Methods › Reliability and stability measures ↔ ToShareCode.zip/ToShareCode/Spatial_Topology/DiceScript.m, lines 62–145 · score 0.69 · Normalized Mutual Information, perfect overlap, Dice coefficients, NMI, Network
  3. [3] § Results › Network topology ↔ ToShareCode.zip/ToShareCode/FC_Reliability/makeFC.m, lines 165–227 · score 0.58 · visual lateral, default retrosplenial, streams, language, salience, Networks
  4. [4] § Methods › Data preprocessing and denoising ↔ ToShareCode.zip/ToShareCode/Miscellaneous/plot_sizes.m, lines 140–178 · score 0.57 · Signal complexity, Head motion, displacement, correlation, networks
  5. [5] § Methods › Reliability and stability measures ↔ ToShareCode.zip/ToShareCode/Spatial_Topology/Patch_Pipeline/diceForPatches.m, lines 515–601 · score 0.56 · Normalized Mutual Information, Dice coefficients, NMI, Network
  6. [6] § Results › Network topology ↔ ToShareCode.zip/ToShareCode/Miscellaneous/plot_sizes.m, lines 140–178 · score 0.52 · signal complexity, head motion, FDR, correlations, Networks

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 · 797 lines · 28 KB · no license · 2 matches

  1. % CREATE FC MATRICES - first three cells
  2. % silhouette calculations on FC matrices - fourth cell
  3. % FC properties - fifth cell
  4. %% (1) set up
  5. addpath(genpath('~/Documents/Randy/APM_Imaging/processingPFM/'))
  6. addpath('/Users/isaactreves/MIT Dropbox/Isaac Treves/BACKUP_LOCAL_MIT_2/MATLAB/GraphVar_2.03/GraphVar_2.03a/src/ext/BCT');
  7. addpath(genpath('~/Documents/Homologous-Functional-Regions/HFR_ai/'))
  8. ProgramPath = '~/Documents/Homologous-Functional-Regions/HFR_ai';
  9. % add freesurfer code to matlab path
  10. addpath(genpath('/Applications/MATLAB_R2024a.app/bin'))
  11. addpath(genpath('/Users/isaactreves/Documents/CBIG'))
  13. % ADJUST HERE
  14. sub='sub-abd1113';
  15. pfmpath=['/Users/isaactreves/Documents/Randy/APM_Imaging/' sub '-pfm-sparse/pfmsplitcomparison/'];
  17. columns = 21;
  18. load("~/Documents/PFM-Depression/PFM-Tutorial/Utilities/priors.mat")
  19. labels=Priors.NetworkLabels(1:21);
  20. Templates=Priors.Spatial;
  21. ArbitraryCutoff=0.25;
  22. Templates(Templates>=0.25)=1;
  23. Templates(Templates<0.25)=0;
  25. tc = ft_read_cifti_mod([pfmpath 'concatenated_regressed_smoothed2.55_32k_fsLR.dtseries.nii']);
  26. tc.brainstructure(tc.brainstructure==-1)=[];
  27. tc.data = tc.data(tc.brainstructure<3,:);
  28. %tc_corr=fisherz(corr(tc.data'));
  29. %%
  30. %networksload = ft_read_cifti_mod([pfmpath 'Bipartite_PhysicalCommunities+FinalLabeling.dlabel.nii']);
  31. % adjustment to use networks only from small amount
  32. pfmpath=['/Users/isaactreves/Documents/Randy/APM_Imaging/' sub '-pfm-sparse/pfmsplit400/'];
  33. networksload = ft_read_cifti_mod([pfmpath 'Bipartite_PhysicalCommunities+AlgorithmicLabeling.dlabel.nii']);
  34. %networksload = niftiread(['/Users/isaactreves/Documents/Randy/APM_Imaging/bids/derivatives/sub-abd1006/CIFTI/pfm_charles/Bipartite_PhysicalCommunities+AlgorithmicLabeling.dlabel.nii']);
  36. networks=struct;
  37. networks.data=networksload.data;
  38. % networks.data=networks.data +1 ;
  40. %Ic = ft_read_cifti_mod(['/Users/isaactreves/Documents/Randy/APM_Imaging/' sub{1} '/pfm/Bipartite_PhysicalCommunities+SpatialFiltering.dtseries.nii']);
  41. communities=ft_read_cifti_mod([pfmpath 'Bipartite_PhysicalCommunities+AlgorithmicLabeling_InfoMapCommunities.dlabel.nii']);
  42. if isfile([pfmpath 'Bipartite_PhysicalCommunities+FinalLabeling_NetworkLabels.xls'])
  43. [assignment]=readtable([pfmpath 'Bipartite_PhysicalCommunities+FinalLabeling_NetworkLabels.xls']);
  44. "manual found"
  45. else
  46. [assignment]=readtable([pfmpath 'Bipartite_PhysicalCommunities+AlgorithmicLabeling_NetworkLabels.xls']);
  47. end
  48. assignment=table2cell(assignment); % fine
  49. communities.brainstructure(communities.brainstructure==-1)=[];
  50. communities.data=communities.data(communities.brainstructure<3,:);
  51. %% (2)
  52. maps_data=communities.data';
  53. [labeledVector, uniqueLabels] = niftiMapsToLabeledVector(maps_data,assignment);
  55. %% (3)
  56. [refFC,ref_network]=createStructuredFCMatrix(tc,labeledVector,Priors.NetworkColors,0);
  59. %% (4) silhouette - adjust output
  60. outputpath=['/Users/isaactreves/Documents/Randy/APM_Imaging/outputs/' sub '/whole/Silhouette.csv'];
  61. FC = refFC;
  62. network_mapping= ref_network;
  63. % FC is an n x n matrix of correlations between patches
  64. % network_mapping is a structure with fields corresponding to network names
  65. % Each field contains a vector of patch indices belonging to that network
  67. % Get all network names
  68. network_names = network_mapping.keys;
  69. num_networks = length(network_names);
  71. % Initialize output array for network silhouette indices
  72. network_silhouette_indices = zeros(num_networks, 1);
  73. network_si_std = zeros(num_networks, 1);
  75. % not FIsher z transformed, so this works
  76. % Convert correlation matrix to distance matrix (1 - correlation)
  77. % Assuming correlation ranges from -1 to 1, distance will range from 0 to 2
  78. distance_matrix = 1 - FC;
  80. % For each network, calculate silhouette index
  81. for net_idx = 1:num_networks
  82. current_network = network_names{net_idx};
  83. current_patches = network_mapping(current_network).patches;
  84. startIdx = network_mapping(current_network).startIdx
  86. % Initialize array to store silhouette indices for each patch in current network
  87. patch_si = zeros(length(current_patches), 1);
  89. % Calculate silhouette index for each patch in the network
  90. for i = 1:length(current_patches)
  91. patch = current_patches(i)+startIdx-1;
  93. % Calculate mean within-network distance (a)
  94. current_patches_adjusted=current_patches+startIdx-1;
  95. within_patches = current_patches_adjusted(current_patches_adjusted ~= patch); % All patches in same network except current
  96. if isempty(within_patches)
  97. a = 0; % If only one patch in network, set a to 0
  98. else
  99. a = mean(distance_matrix(patch, within_patches));
  100. end
  102. % Initialize array to store mean distances to other networks
  103. b_values = zeros(num_networks - 1, 1);
  104. b_idx = 1;
  106. % Calculate mean distance to each other network
  107. % for other_net_idx = 1:num_networks
  108. % if other_net_idx ~= net_idx
  109. % other_network = network_names{other_net_idx};
  110. % other_patches = network_mapping(other_network).patches;
  111. % other_idx=network_mapping(other_network).startIdx;
  112. % if ~isempty(other_patches)
  113. % b_values(b_idx) = mean(distance_matrix(patch, other_patches+other_idx-1));
  114. % b_idx = b_idx + 1;
  115. % end
  116. % end
  117. % end
  119. % b is the minimum mean distance to other networks - seems
  120. % restrictive - this is like the closest other network
  121. % b = min(b_values);
  122. % Collect all patches not in the current network
  123. all_other_patches = [];
  124. for other_net_idx = 1:num_networks
  125. if other_net_idx ~= net_idx
  126. other_network = network_names{other_net_idx};
  127. other_patches_relative = network_mapping(other_network).patches;
  128. other_start_idx = network_mapping(other_network).startIdx;
  130. % Convert relative indices to absolute indices
  131. other_patches_absolute = other_patches_relative + other_start_idx - 1;
  132. all_other_patches = [all_other_patches; other_patches_absolute];
  133. end
  134. end
  136. % Calculate mean distance to all patches outside the current network
  137. if isempty(all_other_patches)
  138. b = 0 % Should not happen if there are multiple networks
  139. else
  140. b = mean(distance_matrix(patch, all_other_patches));
  141. end
  143. % Calculate silhouette index
  144. patch_si(i) = (b - a) / max(a, b);
  145. end
  147. % Store mean silhouette index for network
  148. network_silhouette_indices(net_idx) = mean(patch_si);
  149. if length(current_patches)<2
  150. network_silhouette_indices(net_idx)=0;
  151. end
  152. % Store standard deviation of silhouette indices for network (extra credit)
  153. network_si_std(net_idx) = std(patch_si);
  154. end
  156. % Create table for output
  157. result_table = table(network_names', network_silhouette_indices, network_si_std, ...
  158. 'VariableNames', {'Network', 'SilhouetteIndex', 'SilhouetteStd'});
  160. % Display results
  161. disp('Network Silhouette Indices:');
  162. disp(result_table);
  163. writetable(result_table,outputpath)
  165. %% save to scalar
  166. custom_order = {...
  167. 'Default_Parietal', ...
  168. 'Default_Anterolateral', ...
  169. 'Default_Dorsolateral', ...
  170. 'Default_Retrosplenial', ...
  171. 'Visual_Lateral', ...
  172. 'Visual_Dorsal/VentralStream', ...
  173. 'Visual_V5', ...
  174. 'Visual_V1', ...
  175. 'Frontoparietal', ...
  176. 'DorsalAttention', ...
  177. 'Premotor/DorsalAttentionII', ...
  178. 'Language', ...
  179. 'Salience', ...
  180. 'CinguloOpercular/Action-mode', ...
  181. 'MedialParietal', ...
  182. 'Somatomotor_Hand', ...
  183. 'Somatomotor_Face', ...
  184. 'Somatomotor_Foot', ...
  185. 'Auditory', ...
  186. 'SomatoCognitiveAction', ...
  187. 'Noise'};
  189. % Create table for output
  190. result_table = table(network_names', network_silhouette_indices, network_si_std, ...
  191. 'VariableNames', {'Network', 'SilhouetteIndex', 'SilhouetteStd'});
  193. % Sort the table according to the custom order
  194. [~, sort_idx] = ismember(custom_order, result_table.Network);
  195. valid_idx = sort_idx > 0; % In case some networks in custom_order are not in the data
  196. sort_idx = sort_idx(valid_idx);
  197. sorted_table = result_table(sort_idx, :);
  199. % Display results
  200. disp('Network Silhouette Indices:');
  201. disp(sorted_table);
  204. %
  205. network_silhouette_indices = sorted_table.SilhouetteIndex;
  207. networkscopy=networksload;
  209. for i =1:length(network_silhouette_indices)
  210. networkscopy.data(networkscopy.data==i)=network_silhouette_indices(i);
  211. end
  212. networkscopy.data(networkscopy.data==21)=0;
  213. networkscopy.mapname={'bananas'};
  214. ft_write_cifti_mod([pfmpath 'SI_mean.dscalar.nii'],networkscopy)
  217. % If you need the indices in the sorted order as return value
  218. network_silhouette_indices = sorted_table.SilhouetteStd;
  220. networkscopy=networksload;
  222. for i =1:length(network_silhouette_indices)
  223. networkscopy.data(networkscopy.data==i)=network_silhouette_indices(i);
  224. end
  225. networkscopy.data(networkscopy.data==21)=0;
  226. networkscopy.mapname={'bananas'};
  227. ft_write_cifti_mod([pfmpath 'SI_std.dscalar.nii'],networkscopy)
  228. %% loop through - deprecated
  229. Stopping=100:100:1600;
  230. pfm_ref=['/Users/isaactreves/Documents/Randy/APM_Imaging/' sub{1} '/CIFTIS/pfm/'];
  231. tc = ft_read_cifti_mod([pfm_ref 'concatenated_regressed_smoothed2.55_32k_fsLR.dtseries.nii']);
  232. tc.brainstructure(tc.brainstructure==-1)=[];
  233. tc.data = tc.data(tc.brainstructure<3,:);
  234. %tc_corr=fisherz(corr(tc.data'));
  235. metrics_all=[];
  236. for i = Stopping
  237. pfmpath=['/Users/isaactreves/Documents/Randy/APM_Imaging/' sub{1} '/CIFTIS/pfm' num2str(i) '/'];
  238. [labeledVector]=prepare_pfm(pfmpath);
  239. [targetFC,target_network]=createStructuredFCMatrix(tc,labeledVector,Priors.NetworkColors,1);
  240. [metrics,networksBoth]=calculateFCSimilarity(refFC,ref_network,targetFC,target_network);
  241. close all;
  242. metrics_all=[metrics_all;metrics];
  243. end
  244. %% (5) actual method - run reference FC stuff first
  245. Stopping=100:100:1600;
  246. Stop_end=1600; %for 1065
  247. % remember for split you have to use other data like split1600
  248. % for rest, you should use pfm rest ,same data
  249. pfm_target=['/Users/isaactreves/Documents/Randy/APM_Imaging/' sub '-pfm-sparse/pfmsplit1600/'];
  250. tc = ft_read_cifti_mod([pfm_target 'concatenated_regressed_smoothed2.55_32k_fsLR.dtseries.nii']);
  251. [labeledVector]=prepare_pfm(pfm_target); % I think we only need to do this once - unless you want to get fancy and use different
  252. %%
  253. % parcellations for each data chunk
  254. tc.brainstructure(tc.brainstructure==-1)=[];
  255. tc.data = tc.data(tc.brainstructure<3,:);
  256. %tc_corr=fisherz(corr(tc.data'));
  257. metrics_all=nan(length(Stopping),4);
  258. shuffles=100;
  259. delete(gcp('nocreate'))
  260. parpool(3) % reminder this is RAM intensive
  261. % not in parallel does 5-8 shuffles / min
  262. parfor i = 1:length(Stopping)
  263. i
  264. if Stopping(i)<Stop_end
  266. allshufflemetrics=[];
  267. for j = 1:shuffles
  268. max_start = size(tc.data,2) - Stopping(i) + 1;
  269. start_col = randi(max_start);
  270. tc_copy=tc;
  271. tc_copy.data=tc_copy.data(:,start_col:(start_col+Stopping(i)-1));
  272. [targetFC,target_network]=createStructuredFCMatrix(tc_copy,labeledVector,Priors.NetworkColors,1);
  273. [metrics,networksBoth]=calculateFCSimilarity(refFC,ref_network,targetFC,target_network);
  274. allshufflemetrics=[allshufflemetrics;metrics];
  275. end
  276. close all;
  277. metrics=mean(allshufflemetrics,1);
  278. metrics_all(i,:)=metrics;
  279. end
  280. end
  283. %%
  284. metricsAll_save=array2table(metrics_all,"VariableNames",{'FC_sim','PartCoef_Sim','Global Eff Diff','Modularity Diff'})
  285. writetable(metricsAll_save,['/Users/isaactreves/Documents/Randy/APM_Imaging/outputs/' sub '/metrics_FC_split_10min.csv'])
  286. %%
  287. function [labeledVector]=prepare_pfm(pfmpath)
  288. communities=ft_read_cifti_mod([pfmpath 'Bipartite_PhysicalCommunities+AlgorithmicLabeling_InfoMapCommunities.dlabel.nii']);
  289. if isfile([pfmpath 'Bipartite_PhysicalCommunities+FinalLabeling_NetworkLabels.xls'])
  290. [assignment]=readtable([pfmpath 'Bipartite_PhysicalCommunities+FinalLabeling_NetworkLabels.xls']);
  291. "manual found prepare pfm"
  292. else
  293. [assignment]=readtable([pfmpath 'Bipartite_PhysicalCommunities+AlgorithmicLabeling_NetworkLabels.xls']);
  294. end
  295. assignment=table2cell(assignment);
  296. communities.brainstructure(communities.brainstructure==-1)=[];
  297. communities.data=communities.data(communities.brainstructure<3,:);
  298. %%
  299. maps_data=communities.data';
  300. [labeledVector, uniqueLabels] = niftiMapsToLabeledVector(maps_data,assignment);
  301. end
  303. function [labeledVector, uniqueLabels] = niftiMapsToLabeledVector(data, labelCell)
  304. % Converts a set of NIFTI maps to a 1D labeled vector where vertices with
  305. % value > 0 in a map receive the corresponding label from a cell array
  306. %
  307. % Inputs:
  308. % data - NIFTI data with maps in dim 5 and vertices in dim 6
  309. % labelCell - Cell array where rows 2 to (2+numMaps-1) contain labels:
  310. % column 2 has primary labels, column 3 (if exists)
  311. % takes precedence when non-empty
  312. %
  313. % Outputs:
  314. % labeledVector - Cell array with two columns:
  315. % Column 1: Network label (e.g., 'Frontoparietal', 'Default')
  316. % Column 2: Community number for that network (e.g., 1, 2, ...)
  317. % uniqueLabels - List of unique network labels used in the resulting vector
  319. % Get dimensions
  320. numMaps = size(data, 1)
  321. numVertices = size(data, 2)
  323. % Check if labelCell has sufficient rows
  324. if size(labelCell, 1) < numMaps
  325. error('Label cell array must have at least %d rows', numMaps + 1);
  326. end
  328. % Initialize output vector (empty cells will represent unlabeled vertices)
  329. labeledVector = cell(numVertices, 2);
  330. for i = 1:numVertices
  331. labeledVector{i,1} = '';
  332. labeledVector{i,2} = 0;
  333. end
  335. % Track communities for each label
  336. labelCommunities = containers.Map('KeyType', 'char', 'ValueType', 'any');
  338. % Process each map
  339. for mapIdx = 1:numMaps
  340. % Get label for this map from cell array
  341. labelRowIdx = mapIdx; % Skip first row
  343. % Determine label - column 3 takes precedence if it exists and is non-empty
  344. if size(labelCell, 2) >= 3 && ~isempty(labelCell{labelRowIdx, 3}) && all(~isnan(labelCell{labelRowIdx, 3}))
  345. mapLabel = labelCell{labelRowIdx, 3};
  346. else
  347. mapLabel = labelCell{labelRowIdx, 2};
  348. end
  350. % Skip if label is empty
  351. if isempty(mapLabel)
  352. continue;
  353. end
  355. % Convert to string if it's a number
  356. if isnumeric(mapLabel)
  357. mapLabel = num2str(mapLabel);
  358. end
  360. % Update community counter for this label
  361. if ~isKey(labelCommunities, mapLabel)
  362. labelCommunities(mapLabel) = 1;
  363. else
  364. labelCommunities(mapLabel) = labelCommunities(mapLabel) + 1;
  365. end
  367. % Get current community number for this label
  368. communityNumber = labelCommunities(mapLabel);
  370. % Find vertices where the current map has value > 0
  371. activeVertices = data(mapIdx,:) > 0;
  373. % If data is high-dimensional, ensure we get correct vertex indices
  374. if ndims(activeVertices) > 1
  375. activeVertices = reshape(activeVertices, [], numVertices);
  376. activeVertices = any(activeVertices, 1)';
  377. end
  379. % Assign label and community to active vertices
  380. for i = find(activeVertices)'
  381. labeledVector{i,1} = mapLabel;
  382. labeledVector{i,2} = communityNumber;
  383. end
  384. end
  386. % Get list of unique network labels (excluding empty/unlabeled)
  387. networkLabels = cellfun(@(x) num2str(x), labeledVector(:,1), 'UniformOutput', false);
  388. validNetworks = ~cellfun(@isempty, networkLabels);
  390. if any(validNetworks)
  391. uniqueLabels = unique(networkLabels(validNetworks));
  392. else
  393. uniqueLabels = {};
  394. end
  396. % Remove empty string from uniqueLabels if present
  397. if ~isempty(uniqueLabels) && any(cellfun(@isempty, uniqueLabels))
  398. uniqueLabels(cellfun(@isempty, uniqueLabels)) = [];
  399. end
  400. end
  402. function numericalVector = convertLabelsToNumerical(labeledVector, networkLabels)
  403. % Converts a cell array of network names to a numerical vector using a
  404. % reference cell array of network labels
  405. %
  406. % Inputs:
  407. % labeledVector - Cell array of network names (e.g., {'Frontoparietal', 'Default', ...})
  408. % networkLabels - Cell array containing network names indexed from 1 to N
  409. %
  410. % Output:
  411. % numericalVector - Vector where each element is the numerical index
  412. % of the corresponding network in networkLabels
  414. % Initialize output vector with zeros (for unmatched labels)
  415. numericalVector = zeros(size(labeledVector));
  417. % Convert to cell array if input is not already a cell array
  418. if ~iscell(labeledVector)
  419. error('labeledVector must be a cell array of network names');
  420. end
  422. % Ensure networkLabels is a cell array
  423. if ~iscell(networkLabels)
  424. error('networkLabels must be a cell array');
  425. end
  427. % Process each label in the input vector
  428. for i = 1:length(labeledVector)
  429. % Skip empty cells
  430. if isempty(labeledVector{i})
  431. continue;
  432. end
  434. % Find the index of the current label in the networkLabels array
  435. idx = find(strcmpi(labeledVector{i}, networkLabels));
  437. % If found, assign the index to the output vector
  438. if ~isempty(idx)
  439. numericalVector(i) = idx(1); % Use first match if multiple found
  440. else
  441. warning('Network "%s" not found in networkLabels', labeledVector{i});
  442. end
  443. end
  444. end
  446. function [FC,networkMap]=createStructuredFCMatrix(data, labeledVector, networkColors,noPlot)
  447. % data: struct with field 'data' of size N vertices x T timepoints
  448. % labeledVector: cell array of size N x 2 with network and patch info
  449. % networkColors: matrix of size 21 x 3 with RGB colors (first entry to be ignored)
  451. % Extract necessary data
  452. vertexData = data.data; %
  453. networkIDs = labeledVector(:, 1);
  454. patchIDs = labeledVector(:, 2);
  456. % Get unique networks
  457. % Handle potential issues with cell arrays
  458. if iscell(networkIDs)
  459. uniqueNetworks = unique(networkIDs);
  460. else
  461. uniqueNetworks = unique(num2cell(networkIDs));
  462. end
  463. load("~/Documents/PFM-Depression/PFM-Tutorial/Utilities/priors.mat")
  464. labels=Priors.NetworkLabels(1:21);
  465. [~, indices] = ismember(uniqueNetworks, labels);
  467. % Sort foundNetworks based on these indices
  468. [~, sortOrder] = sort(indices);
  469. uniqueNetworks = uniqueNetworks(sortOrder);
  471. numNetworks = length(uniqueNetworks);
  473. % Create mapping of each vertex to its network and patch
  474. networkMap = containers.Map('KeyType', 'char', 'ValueType', 'any');
  475. patchMap = containers.Map('KeyType', 'char', 'ValueType', 'double');
  476. patchCount = 0;
  478. % First, identify all patches and count them
  479. for i = 1:numNetworks
  480. currentNetwork = uniqueNetworks{i};
  482. % Handle different data types for comparison
  483. if iscell(networkIDs)
  484. networkVertices = strcmp(networkIDs, currentNetwork);
  485. else
  486. networkVertices = (networkIDs == currentNetwork);
  487. end
  489. % Extract patches for this network safely
  490. networkPatchIDs = patchIDs(networkVertices);
  492. % Handle different data types for unique operation
  493. if iscell(networkPatchIDs)
  494. uniquePatches = unique(cell2mat(networkPatchIDs));
  495. else
  496. uniquePatches = unique(networkPatchIDs); % Convert to cell if numeric
  497. end
  499. % Store patch information for this network
  500. if ischar(currentNetwork)
  501. networkKey = currentNetwork;
  502. else
  503. networkKey = num2str(currentNetwork);
  504. end
  506. networkMap(networkKey) = struct('patches', uniquePatches, 'startIdx', patchCount+1);
  508. % Map each patch to an index for the FC matrix
  509. for j = 1:length(uniquePatches)
  510. patchCount = patchCount + 1;
  511. if iscell(uniquePatches)
  512. patchID = uniquePatches{j};
  513. else
  514. patchID = uniquePatches(j);
  515. end
  517. % Create a string key regardless of input type
  518. if ischar(patchID)
  519. patchKey = patchID;
  520. else
  521. patchKey = num2str(patchID);
  522. end
  524. patchMap([networkKey '_' patchKey]) = patchCount;
  525. end
  526. end
  527. % debugged up to here
  528. % Create FC matrix (correlation matrix)
  529. FC = zeros(patchCount, patchCount);
  531. % For each patch, compute average time series
  532. patchTimeSeries = zeros(patchCount, size(vertexData, 2));
  534. for i = 1:numNetworks
  535. currentNetwork = uniqueNetworks{i};
  537. % Convert to string key for map lookup
  538. if ischar(currentNetwork)
  539. networkKey = currentNetwork;
  540. else
  541. networkKey = num2str(currentNetwork);
  542. end
  544. patches = networkMap(networkKey).patches;
  546. for j = 1:length(patches)
  547. if iscell(patches)
  548. currentPatch = patches{j};
  549. else
  550. currentPatch = patches(j);
  551. end
  554. % Handle different data types for comparison
  555. patchIDs_forcompare=cell2mat(patchIDs);
  556. if iscell(networkIDs) && iscell(patchIDs_forcompare)
  557. % doesn't work
  558. patchVertices = strcmp(networkIDs, currentNetwork) & strcmp(patchIDs_forcompare, currentPatch);
  559. elseif ~iscell(networkIDs) && ~iscell(patchIDs_forcompare)
  560. patchVertices = (networkIDs == currentNetwork) & (patchIDs_forcompare == currentPatch);
  561. elseif iscell(networkIDs) && ~iscell(patchIDs_forcompare)
  562. patchVertices = strcmp(networkIDs, currentNetwork) & (patchIDs_forcompare == currentPatch);
  563. else
  564. patchVertices = (networkIDs == currentNetwork) & strcmp(patchIDs_forcompare, currentPatch);
  565. end
  567. % Average time series for this patch
  568. if any(patchVertices)
  569. % Convert to string for map lookup
  570. if ischar(currentPatch)
  571. patchKey = currentPatch;
  572. else
  573. patchKey = num2str(currentPatch);
  574. end
  576. patchIdx = patchMap([networkKey '_' patchKey]);
  577. patchTimeSeries(patchIdx, :) = mean(vertexData(patchVertices, :), 1);
  578. end
  579. end
  580. end
  582. % Compute correlation matrix (FC)
  583. FC = corrcoef(patchTimeSeries');
  585. if ~noPlot
  586. % Create figure
  587. figure('Position', [100, 100, 1200, 1000]);
  589. % Plot the FC matrix
  590. axes('Position', [0.1, 0.1, 0.8, 0.8]);
  591. imagesc(FC, [-1, 1]);
  592. %colormap(redbluecmap);
  593. colorbar;
  594. end
  595. % Create labels for the axes
  596. axisLabels = cell(patchCount, 1);
  597. networkBoundaries = zeros(numNetworks+1, 1);
  598. networkBoundaries(1) = 0.5;
  600. % Tick positions and labels
  601. tickPositions = [];
  602. networkCenters = [];
  604. for i = 1:numNetworks
  605. currentNetwork = uniqueNetworks{i};
  607. % Convert to string key for map lookup
  608. if ischar(currentNetwork)
  609. networkKey = currentNetwork;
  610. else
  611. networkKey = num2str(currentNetwork);
  612. end
  614. patches = networkMap(networkKey).patches;
  615. startIdx = networkMap(networkKey).startIdx;
  617. % Calculate center position for this network's label
  618. networkCenter = startIdx + (length(patches)-1)/2;
  619. networkCenters = [networkCenters; networkCenter];
  621. for j = 1:length(patches)
  622. if iscell(patches)
  623. currentPatch = patches{j};
  624. else
  625. currentPatch = patches(j);
  626. end
  628. % Convert to string for map lookup
  629. if ischar(currentPatch)
  630. patchKey = currentPatch;
  631. else
  632. patchKey = num2str(currentPatch);
  633. end
  635. patchIdx = patchMap([networkKey '_' patchKey]);
  637. % Format label based on type
  638. if ischar(currentPatch)
  639. axisLabels{patchIdx} = currentPatch;
  640. else
  641. axisLabels{patchIdx} = num2str(currentPatch);
  642. end
  644. tickPositions = [tickPositions; patchIdx];
  645. end
  647. networkBoundaries(i+1) = startIdx + length(patches) - 0.5;
  648. end
  649. if ~noPlot
  650. % Set axis ticks and labels
  651. set(gca, 'XTick', tickPositions, 'XTickLabel', axisLabels, 'XTickLabelRotation', 90,'FontSize',8);
  652. set(gca, 'YTick', tickPositions, 'YTickLabel', axisLabels);
  654. % Draw network boundaries
  655. hold on;
  656. for i = 1:length(networkBoundaries)
  657. % Horizontal lines
  658. if i > 1
  659. line([0.5, patchCount+0.5], [networkBoundaries(i), networkBoundaries(i)], 'Color', 'k', 'LineWidth', 2);
  660. end
  661. % Vertical lines
  662. if i > 1
  663. line([networkBoundaries(i), networkBoundaries(i)], [0.5, patchCount+0.5], 'Color', 'k', 'LineWidth', 2);
  664. end
  665. end
  667. % Add network colors on the left and top using the provided RGB values
  668. axisWidth = 0.05;
  670. % Colors on the left
  671. axes('Position', [0.05, 0.1, axisWidth, 0.8]);
  672. colorMap = zeros(patchCount, 3);
  673. end
  674. % problem is here
  675. for i = 1:numNetworks
  676. currentNetwork = uniqueNetworks{i};
  678. % Convert to string key for map lookup
  679. if ischar(currentNetwork)
  680. networkKey = currentNetwork;
  681. else
  682. networkKey = num2str(currentNetwork);
  683. end
  685. startIdx = networkMap(networkKey).startIdx;
  686. patches = networkMap(networkKey).patches;
  687. numPatches = length(patches);
  689. % Get the color for this network
  690. colorIdx=find(strcmp(labels,networkKey));
  691. %colorIdx = min(i, size(networkColors, 1));
  692. networkColor = networkColors(colorIdx, :);
  694. % Fill the color block for this network
  695. colorMap(startIdx:(startIdx+numPatches-1), :) = repmat(networkColor, numPatches, 1);
  696. end
  697. if ~noPlot
  698. imagesc((1:patchCount)', 'CDataMapping', 'direct');
  699. colormap(gca, colorMap);
  700. axis off;
  702. % Colors on the top
  703. axes('Position', [0.1, 0.9, 0.76, axisWidth]);
  704. imagesc(1:patchCount, 'CDataMapping', 'direct');
  705. colormap(gca, colorMap);
  706. axis off;
  708. % Return to main axes and label
  709. axes('Position', [0.1, 0.1, 0.8, 0.8], 'Visible', 'off');
  710. title('Functional Connectivity Matrix by Network and Patch', 'FontSize', 16);
  712. % Save as PDF
  713. print('-dpdf', 'FunctionalConnectivity.pdf','-bestfit');
  714. end
  715. end
  717. function [metrics,networksInBoth] = calculateFCSimilarity(refFC, refNetworkMap, testFC, testNetworkMap)
  718. % Calculate similarity between functional connectivity matrices
  719. % across different data lengths, comparing only networks that exist in both
  720. %
  721. % Inputs:
  722. % refFC - Reference FC matrix (NxN, where N is the total number of patches)
  723. % refNetworkMap - Structure with fields for each network containing:
  724. % - patches: indices of patches belonging to the network
  725. % - startIdx: starting index for the network
  726. % testFC - Test FC matrix to compare against reference
  727. % testNetworkMap - Network map for the test data
  728. %
  729. % Outputs:
  730. % similarity - Correlation similarity between average FC values (upper triangle only)
  731. % networksInBoth - Cell array of networks found in both datasets
  733. % Get the network names from both maps
  734. refNetworks = refNetworkMap.keys;
  735. testNetworks = testNetworkMap.keys;
  737. % Find networks that exist in both datasets
  738. networksInBoth = intersect(refNetworks, testNetworks);
  740. % Create average FC matrices for networks in both datasets
  741. numNetworks = length(networksInBoth);
  742. refNetFC_avg = zeros(numNetworks, numNetworks);
  743. testNetFC_avg = zeros(numNetworks, numNetworks);
  744. %% quick do modularity and efficiency
  745. geff1=efficiency_wei(refFC);
  746. geff2=efficiency_wei(testFC);
  747. geff_dif=abs(geff2-geff1);
  748. %
  749. [~,Q1]=modularity_louvain_und_sign(refFC);
  750. [~,Q2]=modularity_louvain_und_sign(testFC);
  751. modularity_dif=abs(Q1-Q2);
  752. % For each pair of networks, calculate the average FC
  753. for i = 1:numNetworks
  754. netI = networksInBoth{i};
  756. % Get true patch indices for network i in reference data
  757. refPatchesI = refNetworkMap(netI).patches + (refNetworkMap(netI).startIdx - 1);
  759. % Get true patch indices for network i in test data
  760. testPatchesI = testNetworkMap(netI).patches + (testNetworkMap(netI).startIdx - 1);
  762. for j = 1:numNetworks
  763. netJ = networksInBoth{j};
  765. % Get true patch indices for network j in reference data
  766. refPatchesJ = refNetworkMap(netJ).patches + (refNetworkMap(netJ).startIdx - 1);
  768. % Get true patch indices for network j in test data
  769. testPatchesJ = testNetworkMap(netJ).patches + (testNetworkMap(netJ).startIdx - 1);
  771. % Calculate average FC between networks i and j
  772. refNetFC_avg(i,j) = mean(mean(refFC(refPatchesI, refPatchesJ)));
  773. testNetFC_avg(i,j) = mean(mean(testFC(testPatchesI, testPatchesJ)));
  774. end
  775. end
  777. % figure;
  778. % imagesc(refNetFC_avg)
  779. % figure;
  780. % imagesc(testNetFC_avg)
  782. % Extract only the upper triangular part (including diagonal)
  783. % to avoid counting symmetric pairs twice
  784. upperTri_mask = triu(true(numNetworks));
  786. % Get values from upper triangular part
  787. refValues = refNetFC_avg(upperTri_mask);
  788. testValues = testNetFC_avg(upperTri_mask);
  790. % Calculate correlation between these values
  791. similarity = corr(refValues, testValues);
  792. participation_coef1=participation_coef_sign(refNetFC_avg,1:length(networksInBoth));
  793. participation_coef2=participation_coef_sign(testNetFC_avg,1:length(networksInBoth));
  794. participation_sim=corrcoef(participation_coef1,participation_coef2,'rows','complete');
  795. participation_sim=participation_sim(1,2);
  796. metrics=[similarity,participation_sim,modularity_dif,geff_dif];
  797. end

makeFC.m, no license · at the source

Overview

Authors: Isaac N Treves1,2, David Pagliaccio1,2, Gaurav H Patel1,3, Reem Tamimi4, Jihoon A Kimerty1,5, Randy P Auerbach1,2, Hilary A Marusak4,6,7
ORCID iDs: Isaac N Treves
  1. Department of Psychiatry, Columbia University Irving Medical Center, New York, NY 10032, USA
  2. Division of Child and Adolescent Psychiatry, New York State Psychiatric Institute, New York, NY 10032, USA
  3. Division of Experimental Therapeutics, New York State Psychiatric Institute, New York, NY 10032, USA
  4. Department of Psychiatry and Behavioral Neurosciences, Wayne State University School of Medicine, 3901 Chrysler Service Dr., Detroit, MI 48201, USA
  5. Department of Psychiatry, Weill Cornell Medicine, New York, NY 10065, USA
  6. Merrill Palmer Skillman Institute for Child and Family Development, Wayne State University, 71 East Ferry Street, Detroit, MI 48202, USA
  7. Department of Pharmacology, Wayne State University School of Medicine, 540 E Canfield St, Detroit, MI 48201, USA
Institutions: Columbia University Irving Medical Center (United States); New York State Psychiatric Institute; Columbia University (United States); Wayne State University (United States); Weill Cornell Medicine (United States)
Journal: Developmental cognitive neuroscience, volume 81, article 101789
Dates: received 26 October 2025; accepted 23 July 2026; published online 24 July 2026; in print July 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1016/j.dcn.2026.101789 · PMID 42508279 · PMCID PMC13446092 · OpenAlex W7170783909
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: fMRI (modality), human (organism), developmental (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Preprocessing, Connectivity, fMRI & imaging
Keywords: Precision functional mapping, PFM, Adolescents, FMRI, Topology, Individualized network, Feasibility
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: NIMH NIH HHS (K01 MH119241, R01 MH132830); NICHD NIH HHS (R21 HD105882)
Citations: not cited yet (Europe PMC); 51 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.

Repository

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

OSF kgvsb

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 6 files
Software Heritage: not checked
Found in: “Data availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: FieldTrip (11 files), Parallel Computing Toolbox (7 files), FreeSurfer (5 files), Statistics and Machine Learning Toolbox (5 files), Connectome Workbench (5 files), Brain Connectivity Toolbox (3 files), NumPy (2 files), pandas (2 files), fdr_bh (Benjamini-Hochberg FDR) (1 file), Image Processing Toolbox (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
29 files
At the source: osf.io/kgvsb/overview?view_only=73be63e38b18407da85a6bb85bd38602

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;
  • 29 scripts, each with its path and the digest of its content;
  • 6 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.

Code and data availability statement

The paper has a code and data 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.1016/j.dcn.2026.101789.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 7 authors, 7 keywords, 2 funders, 49 references.

Cite

This paper

Treves, I. N., Pagliaccio, D., Patel, G. H., Tamimi, R., Kimerty, J. A., Auerbach, R. P., & Marusak, H. A. (2026). Feasibility of precision functional mapping in youth multi-echo fMRI data. Developmental cognitive neuroscience, 81, 101789. https://doi.org/10.1016/j.dcn.2026.101789

BibTeX

@article{treves2026feasibility,
author = {Treves, Isaac N and Pagliaccio, David and Patel, Gaurav H and Tamimi, Reem and Kimerty, Jihoon A and Auerbach, Randy P and Marusak, Hilary A},
title = {{Feasibility of precision functional mapping in youth multi-echo fMRI data}},
journal = {Developmental cognitive neuroscience},
year = {2026},
month = jul,
volume = {81},
pages = {101789},
publisher = {Elsevier},
issn = {1878-9293},
doi = {10.1016/j.dcn.2026.101789},
url = {https://doi.org/10.1016/j.dcn.2026.101789},
pmid = {42508279},
pmcid = {PMC13446092}
}

RIS

TY - JOUR
AU - Treves, Isaac N
AU - Pagliaccio, David
AU - Patel, Gaurav H
AU - Tamimi, Reem
AU - Kimerty, Jihoon A
AU - Auerbach, Randy P
AU - Marusak, Hilary A
TI - Feasibility of precision functional mapping in youth multi-echo fMRI data
T2 - Developmental cognitive neuroscience
J2 - Dev Cogn Neurosci
PY - 2026
DA - 2026/07/24
VL - 81
SP - 101789
SN - 1878-9293
PB - Elsevier
DO - 10.1016/j.dcn.2026.101789
UR - https://doi.org/10.1016/j.dcn.2026.101789
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.dcn.2026.101789",
"type": "article-journal",
"title": "Feasibility of precision functional mapping in youth multi-echo fMRI data",
"container-title": "Developmental cognitive neuroscience",
"author": [
{
"family": "Treves",
"given": "Isaac N"
},
{
"family": "Pagliaccio",
"given": "David"
},
{
"family": "Patel",
"given": "Gaurav H"
},
{
"family": "Tamimi",
"given": "Reem"
},
{
"family": "Kimerty",
"given": "Jihoon A"
},
{
"family": "Auerbach",
"given": "Randy P"
},
{
"family": "Marusak",
"given": "Hilary A"
}
],
"container-title-short": "Dev Cogn Neurosci",
"volume": "81",
"page": "101789",
"DOI": "10.1016/j.dcn.2026.101789",
"PMID": "42508279",
"PMCID": "PMC13446092",
"ISSN": "1878-9293",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.dcn.2026.101789",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
24
]
]
}
}

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.1162/imag.a.1285 [code]
Individualized mapping of functional brain networks in older adulthood.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Connectome Workbench, FieldTrip, Parallel Computing Toolbox, 1 other tool, fMRI, 13 references
[2] doi:10.1016/j.neuron.2026.04.011 [code]
Precision fMRI reveals densely interdigitated network patches with conserved motifs in the lateral prefrontal cortex.
Journal: Neuron
In common: Connectome Workbench, FieldTrip, Parallel Computing Toolbox, 4 other tools, fMRI, 10 references
[3] doi:10.1093/psyrad/kkag013 [code]
Convergent and divergent spatial topographies of individualized brain functional networks and their developmental origins.
Journal: Psychoradiology
In common: Connectome Workbench, FieldTrip, Parallel Computing Toolbox, 2 other tools, developmental, fMRI, 11 references
[4] doi:10.1162/imag.a.1262 [code]
Frame-wise multi-echo distortion correction for superior functional MRI.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Connectome Workbench, FieldTrip, Parallel Computing Toolbox, 5 other tools, fMRI, 7 references
[5] doi:10.1162/imag.a.1222 [code]
Network-based near-scalp personalized brain stimulation targets.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Connectome Workbench, FieldTrip, FreeSurfer, 3 other tools, 9 references
[6] doi:10.21203/rs.3.rs-9326213/v1 [code]
Multi-task fMRI outperforms resting-state fMRI for revealing task-invariant organization of the human brain
Journal: Research Square (preprint)
In common: Connectome Workbench, FreeSurfer, Image Processing Toolbox, 2 other tools, fMRI, 9 references
[7] doi:10.64898/2026.03.09.710558 [code]
Multi-task fMRI outperforms resting-state fMRI for revealing task-invariant organization of the human brain
Journal: bioRxiv (preprint)
In common: Connectome Workbench, FreeSurfer, Image Processing Toolbox, 2 other tools, fMRI, 9 references
[8] doi:10.1038/s41467-026-73668-y [code]
Convergent and divergent brain-cognition development in early adolescence.
Journal: Nature communications
In common: fdr_bh (Benjamini-Hochberg FDR), Connectome Workbench, FieldTrip, 5 other tools, fMRI, 4 references
[9] doi:10.1002/advs.202523009 [code]
Personalized Network-Guided Neuromodulation Enhances Human Working Memory.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: Connectome Workbench, FieldTrip, FreeSurfer, 4 other tools, 6 references
[10] doi:10.1038/s41467-026-74565-0 [code]
The functional neurobiology of dispositions towards negative emotions.
Journal: Nature communications
In common: Brain Connectivity Toolbox, fdr_bh (Benjamini-Hochberg FDR), Connectome Workbench, 7 other tools, 1 reference

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.