OSCR

Prenatal and early postnatal periods differentially shape the maturation of human cortical microstructure and myelin.

Code ↔ Paper

4 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 4 matches · 1 of them tie a paragraph to a whole file without given lines: a weak match, whose lines are not tinted
  1. [1] § Methods › MRI acquisition and preprocessing ↔ code/infant_micro_figures.m, lines 93 to 142 · score 0.67 · parcellated cortical, von Economo, cortical thickness, hemisphere, validation, Schaefer · Report a mismatch
  2. [2] § Methods › Microstructure profiles and their central moments ↔ code/infant_micro_figures.m, lines 93 to 142 · score 0.57 · left hemisphere, microstructure profiles, von Economo, scans, vertex, week · Report a mismatch
  3. [3] § Methods › Participants and demographics ↔ code/infant_micro_figures.m, lines 144 to 212 · score 0.53 · postmenstrual age, gestational age, singleton, scans, birth, infants · Report a mismatch
  4. [4] § Methods › Statistical analysis ↔ code/moment_age_model.m, the whole file · a weak match · score 0.51 · response variable, linear model, validating, predictors, metrics, thickness · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

MATLAB · 1,969 lines · 104 KB · no license · 3 matches

  1. for init_project = 1
  3. GH = '/Users/ttsigaras/Documents/GitHub'; % GitHub directory
  4. homeDir = [GH '/InfantMicro'];
  5. % directory of each participant's native microstructure profile csv
  6. MPs_dir = [homeDir '/data/native_MPs'];
  7. % directory of vertex-wise Von Economo parcel indices
  8. labels_dir_economo = [homeDir '/data/native_economo_ind'];
  9. % directory of vertex-wise Schaefer-200 parcel indices
  10. labels_dir_schaefer200 = [homeDir '/data/native_schaefer200_ind'];
  11. % directory of parcel-wise (Von Economo) cortical thickness
  12. thick_dir = [homeDir '/data/native_thick_economo'];
  14. wk40_tpl_dir = [homeDir '/tpl-week-40']; % week-40 template directory
  16. % add paths for GitHub directories (gifti toolbox and scripts)
  17. addpath(genpath([homeDir '/code']))
  18. addpath(genpath([GH '/gifti']))
  20. % load demographic data
  21. demographics = readtable([homeDir '/data/nnsi01.txt']); % load demographics file
  23. depths = 1:12; % number of depths per profile
  24. % defining valid parcels for each parcellation. Von Economo's
  25. % parcellation includes Corpus Callosum parcels (#2 and #47), which are
  26. % not assigned to any vertex, therefore although the official parcellation
  27. % includes 90 parcels, there are only 88 unique parcel indices when
  28. % extracting the unique indices from the label files. The below defined
  29. % parcels indices are actually the positions of the parcels within the
  30. % total number of parcels and not the actual indices
  31. total_num_parcels_economo = 88;
  32. total_num_parcels_schaefer200 = 201;
  33. excluded_indices_economo = [1, 15, 16, 17, 20, 21, 22, 23, 24, 25, 26, 45, ...
  34. 59, 60, 61, 64, 65, 66, 67, 68, 69, 70]; % non-cortical and limbic parcels
  35. valid_parcels_economo = setdiff(1:total_num_parcels_economo, ...
  36. excluded_indices_economo); % cortical and non-limbic parcel
  37. valid_parcels_schaefer200 = 2:201;
  39. % load week-40 pial surface
  40. wk40_lh_pial = gifti([wk40_tpl_dir '/week-40_hemi-left_space-dhcpSym_dens-32k_pial.surf.gii']);
  41. wk40_rh_pial = gifti([wk40_tpl_dir '/week-40_hemi-right_space-dhcpSym_dens-32k_pial.surf.gii']);
  42. wk40.tri = [wk40_lh_pial.faces; wk40_rh_pial.faces + size(wk40_lh_pial.vertices, 1)];
  43. wk40.coord = [wk40_lh_pial.vertices; wk40_rh_pial.vertices]';
  44. % load Von Economo's week-40 tpl vertex-wise parcel indices
  45. wk40_eco_lh_ind = gifti([wk40_tpl_dir '/week-40_hemi-left_desc-economo.label.gii']);
  46. wk40_eco_rh_ind = gifti([wk40_tpl_dir '/week-40_hemi-right_desc-economo.label.gii']);
  47. wk40_economo_ind = [wk40_eco_lh_ind.cdata + 1; wk40_eco_rh_ind.cdata + 46];
  48. % load Schaefer-200's week-40 tpl vertex-wise parcel indices
  49. wk40_sch200_lh_ind = gifti([wk40_tpl_dir '/week-40_hemi-left_desc-Schaefer2018_7Networks_200.32k.label.gii']);
  50. wk40_sch200_rh_ind = gifti([wk40_tpl_dir '/week-40_hemi-right_desc-Schaefer2018_7Networks_200.32k.label.gii']);
  51. wk40_schaefer200_ind = [wk40_sch200_lh_ind.cdata; wk40_sch200_rh_ind.cdata];
  53. % load colour maps
  54. load([homeDir '/resources/colourmaps/colorbrewer.mat'])
  55. scicol = [homeDir '/resources/colourmaps/ScientificColourMaps7'];
  56. scicol_names = ["devon", "lajolla", "lapaz", "roma", "vik"];
  57. for ii = 1:length(scicol_names)
  58. load(strcat(scicol, '/', scicol_names(ii), '.mat'));
  59. end
  61. % load eigenmodes 2, 3 and 4
  62. eigenmodes = readmatrix([homeDir '/resources/week-40_pial_hemi-left_emode_2-4.txt']);
  63. % load previously used permutation indices
  64. load([homeDir '/resources/spin_test_ind_wk40.mat'], 'perm_id')
  66. end
  68. for preparing_data = 1
  70. variables = {'src_subject_id', 'scan_validation', 'nscan_ga_at_birth_weeks', ...
  71. 'nscan_ga_at_scan_weeks', 'sex'};
  72. filt_demographics = table();
  74. for sub = 1:height(demographics)
  75. sub_id = demographics.src_subject_id{sub}; % subject ID
  76. ses_id = num2str(demographics.scan_validation(sub)); % session ID
  78. MP_file_dir = fullfile(MPs_dir, [sub_id '_ses-' ses_id '_desc-MP.csv']);
  79. left_label_file_dir_economo = fullfile(labels_dir_economo, [sub_id ...
  80. '_ses-' ses_id '_hemi-left_desc-economo.label.gii']);
  81. right_label_file_dir_economo = fullfile(labels_dir_economo, [sub_id ...
  82. '_ses-' ses_id '_hemi-right_desc-economo.label.gii']);
  84. % if all 3 files exist, add subject to the filtered table
  85. if isfile(MP_file_dir) && isfile(left_label_file_dir_economo) && ...
  86. isfile(right_label_file_dir_economo)
  87. filt_demographics = [filt_demographics; demographics(sub, variables)];
  88. else
  89. fprintf('Missing files for sub-%s_ses-%s.\n', sub_id, ses_id);
  90. end
  91. end
  93. % remove subsequent scans of longitudinal participants
  94. filt_demographics.scan_validation = double(filt_demographics.scan_validation);
  95. filt_demographics = sortrows(filt_demographics, ...
  96. {'src_subject_id', 'scan_validation'});
  97. [~, first_idx] = unique(filt_demographics.src_subject_id, 'stable');
  98. filt_demographics = filt_demographics(first_idx, :);
  99. filt_demographics.postnatal_age = filt_demographics.nscan_ga_at_scan_weeks ...
  100. - filt_demographics.nscan_ga_at_birth_weeks; % calculate postnatal age at scan
  101. % excluding participants with PNA > 7
  102. filt_demographics(filt_demographics.postnatal_age > 7, :) = [];
  104. MP_parc_economo = nan(length(depths), total_num_parcels_economo, height(filt_demographics));
  105. MP_parc_schaefer200 = nan(length(depths), total_num_parcels_schaefer200, height(filt_demographics));
  106. thick_parc_economo = nan(total_num_parcels_economo, height(filt_demographics));
  107. for sub = 1:height(filt_demographics)
  108. sub_id = filt_demographics.src_subject_id{sub}; % subject ID
  109. ses_id = num2str(filt_demographics.scan_validation(sub)); % session ID
  111. MP_file_dir = fullfile(MPs_dir, [sub_id '_ses-' ses_id '_desc-MP.csv']);
  112. thick_file_dir = fullfile(thick_dir, [sub_id '_ses-' ses_id '_desc-thickness_parc-economo.csv']);
  113. left_label_file_dir_economo = fullfile(labels_dir_economo, ...
  114. [sub_id '_ses-' ses_id '_hemi-left_desc-economo.label.gii']);
  115. right_label_file_dir_economo = fullfile(labels_dir_economo, ...
  116. [sub_id '_ses-' ses_id '_hemi-right_desc-economo.label.gii']);
  117. left_label_file_dir_schaefer200 = fullfile(labels_dir_schaefer200, ...
  118. [sub_id '_ses-' ses_id '_hemi-left_desc-Schaefer2018_7Networks_200.32k.label.gii']);
  119. right_label_file_dir_schaefer200 = fullfile(labels_dir_schaefer200, ...
  120. [sub_id '_ses-' ses_id '_hemi-right_desc-Schaefer2018_7Networks_200.32k.label.gii']);
  122. MP = readmatrix(MP_file_dir); % read microstructure profiles
  123. sub_thick_parc = readmatrix(thick_file_dir); % read parcellated cortical thickness
  124. left_labels_economo = gifti(left_label_file_dir_economo);
  125. right_labels_economo = gifti(right_label_file_dir_economo);
  126. left_labels_schaefer200 = gifti(left_label_file_dir_schaefer200);
  127. right_labels_schaefer200 = gifti(right_label_file_dir_schaefer200);
  128. left_ind_economo = left_labels_economo.cdata;
  129. % add 45 to the right hemisphere parcel indices to differentiate them from the left hemisphere ones
  130. right_ind_economo = right_labels_economo.cdata + 45;
  131. left_ind_schaefer200 = left_labels_schaefer200.cdata;
  132. right_ind_schaefer200 = right_labels_schaefer200.cdata;
  133. economo_ind = cat(1, left_ind_economo, right_ind_economo); % vertex-wise Von Economo parcel indices
  134. schaefer200_ind = cat(1, left_ind_schaefer200, right_ind_schaefer200); % vertex-wise Schaefer-200 parcel indices
  136. num_of_parc_economo = unique(economo_ind);
  137. num_of_parc_schaefer200 = unique(schaefer200_ind);
  138. if length(num_of_parc_economo) == 88
  139. if length(num_of_parc_schaefer200) == 201
  140. sub_MP_parc_economo = zeros(size(MP, 1), length(num_of_parc_economo));
  141. sub_MP_parc_schaefer200 = zeros(size(MP, 1), length(num_of_parc_schaefer200));
  142. MP(isinf(MP)) = NaN;
  144. % calculate the average intensity per parcel at each depth
  145. for depth = 1:size(MP, 1)
  146. sub_MP_parc_economo(depth,:) = grpstats(MP(depth,:) , economo_ind);
  147. sub_MP_parc_schaefer200(depth,:) = grpstats(MP(depth,:) , schaefer200_ind);
  148. end
  150. MP_parc_economo(:,:,sub) = sub_MP_parc_economo;
  151. MP_parc_schaefer200(:,:,sub) = sub_MP_parc_schaefer200;
  152. thick_parc_economo(:,sub) = sub_thick_parc;
  153. else
  154. fprintf('Subject %s session %s has wrong number of Schaefer parcels.\n', sub_id, ses_id);
  155. end
  156. else
  157. fprintf('Subject %s session %s has wrong number of Von Economo parcels.\n', sub_id, ses_id);
  158. end
  159. end
  161. moments_parc_economo = NaN(5, total_num_parcels_economo, size(MP_parc_economo, 3));
  162. moments_parc_schaefer200 = NaN(5, total_num_parcels_schaefer200, size(MP_parc_schaefer200, 3));
  163. for sub = 1:height(filt_demographics)
  164. MP_economo = MP_parc_economo(:, :, sub);
  165. MP_schaefer200 = MP_parc_schaefer200(:, :, sub);
  166. MP_economo = round(MP_economo*100); % required for the calculate_moments() function to work
  167. MP_schaefer200= round(MP_schaefer200*100); % required for the calculate_moments() function to work
  168. % calculate the moments of this profile
  169. moments_economo = calculate_moments(MP_economo);
  170. moments_schaefer200 = calculate_moments(MP_schaefer200);
  171. moments_parc_economo(:, :, sub) = moments_economo;
  172. moments_parc_schaefer200(:, :, sub) = moments_schaefer200;
  173. end
  174. moments_parc_economo(1,:,:) = moments_parc_economo(1,:,:)/100;
  175. moments_parc_schaefer200(1,:,:) = moments_parc_schaefer200(1,:,:)/100;
  176. % exclude mean intensity, skewness and kurtosis from the moments
  177. moments_parc_economo([1 4 5], :, :) = [];
  178. moments_parc_schaefer200([1 4 5], :, :) = [];
  180. % exclude twins. Twins/Triples have the string "AN", "BN", "CN"
  181. % following the number in their subject ID, whereas singleton infants
  182. % have the string "XX"
  183. ids = string(filt_demographics.src_subject_id);
  184. expression = "CC(?<group>\d{5})(?<twin_id>[A-Z]{2})";
  185. tokens = regexp(ids, expression, 'names');
  186. tokens = vertcat(tokens{:});
  187. twin_id = string({tokens.twin_id}');
  188. is_twin = twin_id ~= "XX";
  189. filt_demographics_no_twins = filt_demographics(~is_twin,:);
  190. moments_parc_no_twins = moments_parc_economo(:,:,~is_twin);
  191. MP_parc_no_twins = MP_parc_economo(:,:,~is_twin);
  193. end
  195. for creating_gen_table = 1
  197. GA = filt_demographics.nscan_ga_at_birth_weeks; % gestational age [weeks]
  198. PNA = filt_demographics.postnatal_age; % postnatal age at scan [weeks]
  199. PMA = filt_demographics.nscan_ga_at_scan_weeks; % postmenstrual age at scan [weeks]
  200. sex = filt_demographics.sex; % sex
  201. sub_id = filt_demographics.src_subject_id; % subject ID
  202. ses_id = filt_demographics.scan_validation; % session ID
  204. GA_no_twins = filt_demographics_no_twins.nscan_ga_at_birth_weeks; % gestational age [weeks]
  205. PNA_no_twins = filt_demographics_no_twins.postnatal_age; % postnatal age at scan [weeks]
  206. PMA_no_twins = filt_demographics_no_twins.nscan_ga_at_scan_weeks; % postmenstrual age at scan [weeks]
  207. sex_no_twins = filt_demographics_no_twins.sex; % sex
  208. sub_id_no_twins = filt_demographics_no_twins.src_subject_id; % subject ID
  209. ses_id_no_twins = filt_demographics_no_twins.scan_validation; % session ID
  211. num_subjects = size(filt_demographics, 1);
  212. num_subjects_no_twins = size(filt_demographics_no_twins, 1);
  214. % repeating entries, in order to match the number of parcels per subject
  215. sub_id_rep_economo = repelem(sub_id, total_num_parcels_economo, 1);
  216. ses_id_rep_economo = repelem(ses_id, total_num_parcels_economo, 1);
  217. GA_rep_economo = repelem(GA, total_num_parcels_economo, 1);
  218. PNA_rep_economo = repelem(PNA, total_num_parcels_economo, 1);
  219. PMA_rep_economo = repelem(PMA, total_num_parcels_economo, 1);
  220. sex_rep_economo = repelem(sex, total_num_parcels_economo, 1);
  221. sub_id_rep_no_twins = repelem(sub_id_no_twins, total_num_parcels_economo, 1);
  222. ses_id_rep_no_twins = repelem(ses_id_no_twins, total_num_parcels_economo, 1);
  223. GA_rep_no_twins = repelem(GA_no_twins, total_num_parcels_economo, 1);
  224. PNA_rep_no_twins = repelem(PNA_no_twins, total_num_parcels_economo, 1);
  225. PMA_rep_no_twins = repelem(PMA_no_twins, total_num_parcels_economo, 1);
  226. sex_rep_no_twins = repelem(sex_no_twins, total_num_parcels_economo, 1);
  227. sub_id_rep_schaefer200 = repelem(sub_id, total_num_parcels_schaefer200, 1);
  228. ses_id_rep_schaefer200 = repelem(ses_id, total_num_parcels_schaefer200, 1);
  229. GA_rep_schaefer200 = repelem(GA, total_num_parcels_schaefer200, 1);
  230. PNA_rep_schaefer200 = repelem(PNA, total_num_parcels_schaefer200, 1);
  231. PMA_rep_schaefer200 = repelem(PMA, total_num_parcels_schaefer200, 1);
  232. sex_rep_schaefer200 = repelem(sex, total_num_parcels_schaefer200, 1);
  234. % multiplying the parcel indices with the number of subjects
  235. parcel_id_rep_economo = repmat(unique(wk40_economo_ind), num_subjects, 1);
  236. parcel_id_rep_no_twins = repmat(unique(wk40_economo_ind), num_subjects_no_twins, 1);
  237. parcel_id_rep_schaefer200 = repmat(unique(wk40_schaefer200_ind), num_subjects, 1);
  239. % reshaping the moments and cortical thickness to match the table's structure
  240. moments_rep_economo = reshape(moments_parc_economo, [2, total_num_parcels_economo * num_subjects])';
  241. thickness_rep = reshape(thick_parc_economo, [total_num_parcels_economo * num_subjects,1]);
  242. moments_rep_no_twins = reshape(moments_parc_no_twins, [2, total_num_parcels_economo * num_subjects_no_twins])';
  243. moments_rep_schaefer200 = reshape(moments_parc_schaefer200, [2, total_num_parcels_schaefer200 * num_subjects])';
  245. % creating the tables
  246. gen_table_economo = table(sub_id_rep_economo, ses_id_rep_economo, GA_rep_economo, ...
  247. PMA_rep_economo, PNA_rep_economo, sex_rep_economo, parcel_id_rep_economo, ...
  248. thickness_rep, moments_rep_economo(:,1), moments_rep_economo(:,2), ...
  249. 'VariableNames', {'subject_id', 'session_id', 'GA', 'PMA', 'PNA', ...
  250. 'sex', 'parcel', 'thick', 'cog', 'variance'});
  251. gen_table_no_twins = table(sub_id_rep_no_twins, ses_id_rep_no_twins, ...
  252. GA_rep_no_twins, PMA_rep_no_twins, PNA_rep_no_twins, sex_rep_no_twins, ...
  253. parcel_id_rep_no_twins, moments_rep_no_twins(:,1), moments_rep_no_twins(:,2), ...
  254. 'VariableNames', {'subject_id', 'session_id', 'GA', 'PMA', 'PNA', ...
  255. 'sex', 'parcel', 'cog', 'variance'});
  256. gen_table_schaefer200 = table(sub_id_rep_schaefer200, ses_id_rep_schaefer200, ...
  257. GA_rep_schaefer200, PMA_rep_schaefer200, PNA_rep_schaefer200, sex_rep_schaefer200, ...
  258. parcel_id_rep_schaefer200, moments_rep_schaefer200(:,1), moments_rep_schaefer200(:,2), ...
  259. 'VariableNames', {'subject_id', 'session_id', 'GA', 'PMA', 'PNA', ...
  260. 'sex', 'parcel', 'cog', 'variance'});
  262. % parcels to exclude from the Von Economo general table
  263. % (i.e. cortical wall, limbic and hippocampal lobes)
  264. exclude_parcels = [1, 16, 17, 18, 21, 22, 23, 24, 25, 26, 27, ...
  265. 46, 61, 62, 63, 66, 67, 68, 69, 70, 71, 72];
  267. % excluding subjects with a PNA lower than 7 weeks and also
  268. % parcels belonging to the cortical wall, hippocampal and limbic lobes
  269. gen_table_economo = gen_table_economo(gen_table_economo.PNA <= 7 & ...
  270. ~ismember(gen_table_economo.parcel, exclude_parcels), :);
  271. gen_table_economo.sex = categorical(gen_table_economo.sex); % making sex a categorical variable
  272. gen_table_no_twins = gen_table_no_twins(gen_table_no_twins.PNA <= 7 & ...
  273. ~ismember(gen_table_no_twins.parcel, exclude_parcels), :);
  274. gen_table_no_twins.sex = categorical(gen_table_no_twins.sex); % making sex a categorical variable
  275. gen_table_schaefer200.sex = categorical(gen_table_schaefer200.sex); % making sex a categorical variable
  277. % averaging central moments across subjects
  278. group_vars = {'subject_id', 'session_id', 'GA', 'PMA', 'PNA', 'sex'};
  279. avg_vars1 = {'thick', 'cog', 'variance'};
  280. avg_vars2 = {'cog', 'variance'};
  281. gen_table_avg_economo = groupsummary(gen_table_economo, group_vars, 'mean', avg_vars1);
  282. gen_table_avg_no_twins = groupsummary(gen_table_no_twins, group_vars, 'mean', avg_vars2);
  283. gen_table_avg_schaefer200 = groupsummary(gen_table_schaefer200, group_vars, 'mean', avg_vars2);
  284. for i = 1:numel(avg_vars1)
  285. old_name = ['mean_' avg_vars1{i}];
  286. newName = avg_vars1{i};
  287. gen_table_avg_economo.Properties.VariableNames{old_name} = newName;
  288. end
  289. for i = 1:numel(avg_vars2)
  290. old_name = ['mean_' avg_vars2{i}];
  291. newName = avg_vars2{i};
  292. gen_table_avg_no_twins.Properties.VariableNames{old_name} = newName;
  293. gen_table_avg_schaefer200.Properties.VariableNames{old_name} = newName;
  294. end
  295. gen_table_avg_economo.GroupCount = [];
  296. gen_table_avg_no_twins.GroupCount = [];
  297. gen_table_avg_schaefer200.GroupCount = [];
  299. uparc_economo = unique(gen_table_economo.parcel); % unique Von Economo parcel indices
  300. uparc_schaefer200 = unique(gen_table_schaefer200.parcel); % unique Schaefer-200 parcel indices
  301. uparc_schaefer200(1) = [];
  303. end
  305. for figure_1 = 1
  307. for panel_A = 1
  309. GA = gen_table_avg_economo.GA;
  310. PNA = gen_table_avg_economo.PNA;
  311. PMA = gen_table_avg_economo.PMA;
  312. sex = gen_table_avg_economo.sex;
  314. [~, sort_idx] = sortrows([PMA, PNA]); % sort first by PMA, then by PNA
  315. GA_sorted = GA(sort_idx);
  316. PMA_sorted = PMA(sort_idx);
  317. sex_sorted = sex(sort_idx);
  318. y_vals = 1:length(PMA_sorted);
  320. % plot each participant's age as a line spanning from their time of
  321. % birth to the day of the scan, separately-coloured for each sex
  322. figure('units','centimeters','outerposition',[0 0 56 20]);
  323. hold on;
  324. for i = 1:length(GA_sorted)
  325. if sex_sorted(i) == 'M'
  326. plot([GA_sorted(i), PMA_sorted(i)], [y_vals(i), y_vals(i)], 'Color', lapaz(70,:), 'LineWidth', 0.5);
  327. scatter(GA_sorted(i), y_vals(i), 10, lapaz(70,:), 'filled', 'MarkerEdgeColor', lapaz(70,:));
  328. scatter(PMA_sorted(i), y_vals(i), 10, lapaz(70,:), 'filled', 'MarkerEdgeColor', lapaz(70,:));
  329. elseif sex_sorted(i) == 'F'
  330. plot([GA_sorted(i), PMA_sorted(i)], [y_vals(i), y_vals(i)], 'Color', lajolla(75,:), 'LineWidth', 0.5);
  331. scatter(GA_sorted(i), y_vals(i), 10, lajolla(75,:), 'filled', 'MarkerEdgeColor', lajolla(75,:));
  332. scatter(PMA_sorted(i), y_vals(i), 10, lajolla(75,:), 'filled', 'MarkerEdgeColor', lajolla(75,:));
  333. end
  334. end
  335. set(gca, 'YTick', []);
  336. set(gca, 'FontSize', 14);
  337. xlim([24 max(PMA_sorted)]);
  338. grid off;
  339. hold off;
  341. end
  343. for panel_B = 1
  345. % plot a scatter plot of PNA against GA, coloured by sex
  346. figure('units','centimeters','outerposition',[0 0 21 20]);
  347. hold on;
  348. scatter(gen_table_avg_economo.GA(gen_table_avg_economo.sex == 'F'), ...
  349. gen_table_avg_economo.PNA(gen_table_avg_economo.sex == 'F'), ...
  350. 50, 'filled', 'MarkerFaceColor', lajolla(75,:));
  351. scatter(gen_table_avg_economo.GA(gen_table_avg_economo.sex == 'M'), ...
  352. gen_table_avg_economo.PNA(gen_table_avg_economo.sex == 'M'), ...
  353. 50, 'filled', 'MarkerFaceColor', devon(80,:));
  354. set(gca, 'FontSize', 14);
  355. ylim([0 max(gen_table_avg_economo.PNA)]);
  356. xlim([min(gen_table_avg_economo.GA) max(gen_table_avg_economo.GA)])
  357. grid off;
  358. hold off;
  360. % plot smoothed histograms of PNA and GA distributions
  361. variables = {'GA', 'PNA'};
  362. for i = 1:length(variables)
  363. var = variables{i};
  364. [f1, xi1] = ksdensity(gen_table_avg_economo.(var)(gen_table_avg_economo.sex == 'M'), 'Bandwidth', 0.4);
  365. [f2, xi2] = ksdensity(gen_table_avg_economo.(var)(gen_table_avg_economo.sex == 'F'), 'Bandwidth', 0.4);
  366. figure('units','centimeters','outerposition',[0 0 28 8]);
  367. hold on;
  368. plot(xi1, f1, 'Color', devon(80,:), 'LineWidth', 2);
  369. fill(xi1, f1, devon(80,:), 'FaceAlpha', 0.3);
  370. plot(xi2, f2, 'Color', lajolla(75,:), 'LineWidth', 2);
  371. fill(xi2, f2, lajolla(75,:), 'FaceAlpha', 0.3);
  372. set(gca, 'XTick', [], 'YTick', []);
  373. xlim([min(gen_table_avg_economo.(var)) max(gen_table_avg_economo.(var))]);
  374. grid on;
  375. hold off;
  376. end
  378. end
  380. end
  382. for figure_2 = 1
  384. for panel_A = 1
  386. % calculate and plot the average profile across all participants
  387. % and parcels
  388. depth = 1:12;
  389. avg_profile = mean(MP_parc_economo(:, :, :), [2, 3]);
  390. figure('units','centimeters','outerposition',[0 0 10 22]);
  391. plot(avg_profile, -depth, 'Color', 'k', 'LineWidth', 1.5);
  392. set(gca, 'YTick', []);
  393. set(gca, 'FontSize', 14);
  394. grid off;
  396. end
  398. for panel_B = 1
  400. cmap = flipud(roma);
  401. nsubs = size(MP_parc_economo, 3);
  402. nbins = 100;
  404. for m = 1:2
  405. profiles_reshaped = reshape(MP_parc_economo, [12, ...
  406. total_num_parcels_economo * nsubs]); % reshape profile array
  407. moments_reshaped = reshape(moments_parc_economo(m,:,:), ...
  408. [total_num_parcels_economo * nsubs, 1]); % flatten moment array
  410. if m == 1
  411. % mask profiles based on CoG value
  412. profiles_reshaped = profiles_reshaped(:, ...
  413. moments_reshaped > 5.8); % arbitrary limit for CoG
  414. % assign each profile to a bin based on their CoG value
  415. mom_bins = discretize(moments_reshaped(moments_reshaped > 5.8), ...
  416. nbins); % arbitrary limit for CoG
  417. else
  418. % mask profiles based on variance value
  419. profiles_reshaped = profiles_reshaped(:, ...
  420. moments_reshaped < 3.4); % arbitrary limit for variance
  421. % assign each profile to a bin based on their variance value
  422. mom_bins = discretize(moments_reshaped(moments_reshaped < 3.4), ...
  423. nbins); % arbitrary limit for variance
  424. end
  426. % plot the average profile within each bin
  427. figure('units','centimeters','outerposition',[0 0 10 22]);
  428. hold on;
  429. for ii = 1:nbins
  430. plot(mean(profiles_reshaped(:,mom_bins==ii),2), -depths, ...
  431. 'Color', cmap(round(ii*height(cmap)/100), :), 'LineWidth', 1)
  432. end
  433. colormap(cmap);
  434. set(gca, 'FontSize', 14);
  435. set(gca, 'YTick', []);
  436. grid off;
  437. hold off;
  439. end
  441. end
  443. for panel_C = 1
  445. moment = ["Centre of gravity", "Variance"];
  446. % calculate average moments per parcel
  447. avg_moments = mean(moments_parc_economo, 3);
  448. % replace non-cortical and limbic parcel values with -100 for
  449. % visualisation purposes
  450. avg_moments(:, ~ismember(1:88, valid_parcels_economo)) = -100;
  452. for mom = 1:length(moment)
  453. min_val = min(avg_moments(mom, avg_moments(mom, :) > -100));
  454. min_val = min_val - min_val/100;
  455. vert_avg_moments = BoSurfStatMakeParcelData(avg_moments(mom,:), wk40, wk40_economo_ind);
  456. figure; SurfStatViewData(vert_avg_moments, wk40, strcat(moment(mom)));
  457. colormap([0.7 0.7 0.7; flipud(lajolla)]);
  458. SurfStatColLim([min_val, max(avg_moments(mom,:))]);
  459. end
  461. end
  463. end
  465. for figure_3 = 1
  467. for panel_A =1
  469. moments = ["Centre of gravity", "Variance"];
  471. for mom = 1:length(moments)
  472. if mom == 1
  473. cog_sign_parcel_data_PMA = moment_age_model("cog", "PMA", ...
  474. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 1, 0);
  475. else
  476. var_sign_parcel_data_PMA = moment_age_model("variance", "PMA", ...
  477. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 1, 0);
  478. end
  479. end
  481. % calculate maximum absolute value for limits
  482. min_val = min([min(cog_sign_parcel_data_PMA(cog_sign_parcel_data_PMA > -100)); ...
  483. min(var_sign_parcel_data_PMA(var_sign_parcel_data_PMA > -100))]);
  484. max_val = max([max(cog_sign_parcel_data_PMA(cog_sign_parcel_data_PMA < 100)); ...
  485. max(var_sign_parcel_data_PMA(var_sign_parcel_data_PMA < 100))]);
  486. abs_lim = max(abs(min_val), abs(max_val));
  487. abs_lim = abs_lim + abs_lim/100;
  489. % plot the significant effects of PMA on each moment
  490. for mom = 1:length(moments)
  491. if mom == 1
  492. sign_parc_data_PMA = cog_sign_parcel_data_PMA;
  493. else
  494. sign_parc_data_PMA = var_sign_parcel_data_PMA;
  495. end
  496. sign_vert_data_PMA = BoSurfStatMakeParcelData(sign_parc_data_PMA, wk40, wk40_economo_ind);
  497. figure; SurfStatViewData(sign_vert_data_PMA, wk40, strcat(moments(mom), ' sign. PMA'));
  498. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  499. % SurfStatColLim([-abs_lim, abs_lim]);
  500. % in the manuscript, we use the same colourscale limits as in
  501. % figure 4 for the purpose of comparison. If you want to use
  502. % limits based only on this data, use the line above
  503. SurfStatColLim([-30.8513, 30.8513]); % limits of figure 4
  504. end
  506. end
  508. for panel_B_profile_plots = 1
  510. moments = ["Centre of gravity", "Variance"];
  511. PMA = gen_table_avg_economo.PMA;
  512. min_PMA = min(PMA);
  513. max_PMA = max(PMA);
  514. % define a window size of 1.5 weeks and 50% overlap
  515. window_size = 1.5;
  516. step = 0.75;
  517. ages = min_PMA:step:max_PMA;
  519. for mom = 1:length(moments)
  520. if mom == 1
  521. sign_parcel_data_PMA = moment_age_model("cog", "PMA", ...
  522. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 1, 0);
  523. else
  524. sign_parcel_data_PMA = moment_age_model("variance", "PMA", ...
  525. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 1, 0);
  526. end
  528. % isolate the parcel-indices with the maximum/minimum significant effects
  529. max_change_parc = find(sign_parcel_data_PMA == ...
  530. max(sign_parcel_data_PMA(sign_parcel_data_PMA < 100)));
  531. min_change_parc = find(sign_parcel_data_PMA == ...
  532. min(sign_parcel_data_PMA(sign_parcel_data_PMA > -100)));
  533. % isolate the parcel-indices with the second maximum/minimum significnat effects
  534. sec_max_change_parc = find(sign_parcel_data_PMA == ...
  535. max(sign_parcel_data_PMA(sign_parcel_data_PMA ...
  536. < max(sign_parcel_data_PMA(sign_parcel_data_PMA<100)))));
  537. sec_min_change_parc = find(sign_parcel_data_PMA == ...
  538. min(sign_parcel_data_PMA(sign_parcel_data_PMA ...
  539. > min(sign_parcel_data_PMA(sign_parcel_data_PMA>-100)))));
  541. % isolate the MPs of the parcels with the first/second
  542. % maximum/minimum significant effects
  543. max_change_profs = MP_parc_economo(:,[max_change_parc ...
  544. sec_max_change_parc],:);
  545. min_change_profs = MP_parc_economo(:,[min_change_parc ...
  546. sec_min_change_parc],:);
  547. profiles_to_plot = [max_change_profs min_change_profs];
  549. % plot the profiles of the parcels with the first/second
  550. % maximum/minimum significant effects on CoG
  551. for i = 1:2
  552. profiles = profiles_to_plot(:,i*2-1:i*2,:);
  553. figure('units','centimeters','outerposition',[0 0 13 17]);
  554. hold on;
  555. cmap = flipud(roma(round(linspace(1, 256, length(ages))), :));
  556. colormap(cmap);
  557. for j = 1:length(ages)-1
  558. age_min = ages(j);
  559. age_max = ages(j) + window_size;
  560. idx = PMA >= age_min & PMA < age_max;
  561. if sum(idx) == 0 % if no profiles within age-range, skip
  562. continue;
  563. end
  564. avg_profile = mean(profiles(:, :, idx), [2, 3]);
  565. plot(avg_profile, -depth, 'Color', cmap(j, :), 'LineWidth', 1.5);
  566. end
  567. colorbar;
  568. clim([min_PMA max_PMA]);
  569. set(gca, 'YTick', []);
  570. set(gca, 'FontSize', 14);
  571. if i==1
  572. if mom == 1
  573. xlim([0.6 1.5]) % arbitrary limits for max cog
  574. else
  575. xlim([0.7 1.9]) % arbitrary limits for max variance
  576. end
  577. else
  578. if mom == 1
  579. xlim([0.6 1.8]) % arbitrary limits for min cog
  580. else
  581. xlim([0.8 1.7]) % arbitrary limits for min variance
  582. end
  583. end
  584. grid off;
  585. hold off;
  586. end
  588. % store the profiles of the parcels with the maximum/minimum
  589. % effects separately for each moment to use in the scatter plots
  590. if mom == 1
  591. max_cog_change_profs = max_change_profs;
  592. min_cog_change_profs = min_change_profs;
  593. else
  594. max_var_change_profs = max_change_profs;
  595. min_var_change_profs = min_change_profs;
  596. end
  597. end
  599. end
  601. for panel_B_scatter_plots = 1
  603. for mom = 1:length(moments)
  604. for i = 1:2
  605. if mom == 1 && i == 1
  606. profiles = max_cog_change_profs;
  607. elseif mom == 1
  608. profiles = min_cog_change_profs;
  609. elseif mom == 2 && i == 1
  610. profiles = max_var_change_profs;
  611. else
  612. profiles = min_var_change_profs;
  613. end
  615. % initialise array to store correlations at each depth
  616. depth_corrs = zeros(length(depths),1);
  617. % calculate correlation of intensities with PMA at each depth
  618. for depth = 1:length(depths)
  619. depth_corrs(depth) = corr(squeeze(mean(profiles(depth,:,:),2)), gen_table_avg_economo.PMA);
  620. end
  622. % plot depth-wise correlations as a scatter plot
  623. figure('units','centimeters','outerposition',[0 0 13 17]); hold on
  624. scatter(depth_corrs, depths, 90, 'filled', 'MarkerFaceColor', vik(80,:),...
  625. 'MarkerEdgeColor', 'k', 'LineWidth', 1.8);
  626. set(gca, 'YDir','reverse', 'YTick', [], 'YLim',[0 13], ...
  627. 'XLim', [0.18 0.75], 'TickLength',[0 0], 'FontSize', 14);
  628. grid off
  629. end
  630. end
  632. end
  634. end
  636. for figure_4 = 1
  638. for controlling_effects = 1
  639. moments = ["Centre of gravity", "Variance"];
  640. % run models which contain both GA and PNA as predictors (effects
  641. % controlling for the effects of the other)
  642. for mom = 1:length(moments)
  643. if mom == 1
  644. cog_sign_parcel_data = moment_age_model("cog", ["GA", "PNA"], ...
  645. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 1, 0);
  646. else
  647. var_sign_parcel_data = moment_age_model("variance", ["GA", "PNA"], ...
  648. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 1, 0);
  649. end
  650. end
  652. % calculate maximum absolute value for limits
  653. min_val = min([min(cog_sign_parcel_data(cog_sign_parcel_data > -100)); ...
  654. min(var_sign_parcel_data(var_sign_parcel_data > -100));]);
  655. max_val = max([max(cog_sign_parcel_data(cog_sign_parcel_data < 100)); ...
  656. max(var_sign_parcel_data(var_sign_parcel_data < 100))]);
  657. abs_lim = max([abs(min_val), abs(max_val)]);
  658. abs_lim = abs_lim + abs_lim/100;
  660. % plot the significant effects of GA and PNA on each moment
  661. for mom = 1:length(moments)
  662. if mom == 1
  663. sign_parcel_data_GA = cog_sign_parcel_data(:,1);
  664. sign_parcel_data_PNA = cog_sign_parcel_data(:,2);
  665. else
  666. sign_parcel_data_GA = var_sign_parcel_data(:,1);
  667. sign_parcel_data_PNA = var_sign_parcel_data(:,2);
  668. end
  670. sign_vert_data_PNA = BoSurfStatMakeParcelData(sign_parcel_data_PNA, wk40, wk40_economo_ind);
  671. figure; SurfStatViewData(sign_vert_data_PNA, wk40, strcat(moments(mom), ' sign. PNA'));
  672. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  673. SurfStatColLim([-abs_lim, abs_lim]);
  674. sign_vert_data_GA = BoSurfStatMakeParcelData(sign_parcel_data_GA, wk40, wk40_economo_ind);
  675. figure; SurfStatViewData(sign_vert_data_GA, wk40, strcat(moments(mom), ' sign. GA'));
  676. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  677. % SurfStatColLim([-abs_lim, abs_lim]);
  678. % in the manuscript, we use the same colourscale limits as the
  679. % figures of the independent models for the purpose of comparison.
  680. % If you want to use limits based only on this data, use the line above
  681. SurfStatColLim([-30.8513, 30.8513]); % limits of independent models
  682. end
  683. end
  685. for no_control = 1
  686. moments = ["Centre of gravity", "Variance"];
  687. % run individual models for GA and PNA as predictors (raw effects)
  688. for mom = 1:length(moments)
  689. if mom == 1
  690. cog_sign_parcel_data_GA = moment_age_model("cog", "GA", ...
  691. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 1, 0);
  692. cog_sign_parcel_data_PNA = moment_age_model("cog", "PNA", ...
  693. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 1, 0);
  694. else
  695. var_sign_parcel_data_GA = moment_age_model("variance", "GA", ...
  696. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 1, 0);
  697. var_sign_parcel_data_PNA = moment_age_model("variance", "PNA", ...
  698. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 1, 0);
  699. end
  700. end
  702. % calculate maximum absolute value for limits
  703. min_val = min([min(cog_sign_parcel_data_GA(cog_sign_parcel_data_GA > -100)); ...
  704. min(cog_sign_parcel_data_PNA(cog_sign_parcel_data_PNA > -100));
  705. min(var_sign_parcel_data_GA(var_sign_parcel_data_GA > -100)); ...
  706. min(var_sign_parcel_data_PNA(var_sign_parcel_data_PNA > -100))]);
  707. max_val = max([max(cog_sign_parcel_data_GA(cog_sign_parcel_data_GA < 100)); ...
  708. max(cog_sign_parcel_data_PNA(cog_sign_parcel_data_PNA < 100));
  709. max(var_sign_parcel_data_GA(var_sign_parcel_data_GA < 100)); ...
  710. max(var_sign_parcel_data_PNA(var_sign_parcel_data_PNA < 100))]);
  711. abs_lim = max([abs(min_val), abs(max_val)]);
  712. abs_lim = abs_lim + abs_lim/100;
  714. % plot the significant effects of GA and PNA on each moment
  715. for mom = 1:length(moments)
  716. if mom == 1
  717. sign_parcel_data_GA = cog_sign_parcel_data_GA;
  718. sign_parcel_data_PNA = cog_sign_parcel_data_PNA;
  719. else
  720. sign_parcel_data_GA = var_sign_parcel_data_GA;
  721. sign_parcel_data_PNA = var_sign_parcel_data_PNA;
  722. end
  724. sign_vert_data_PNA = BoSurfStatMakeParcelData(sign_parcel_data_PNA, wk40, wk40_economo_ind);
  725. figure; SurfStatViewData(sign_vert_data_PNA, wk40, strcat(moments(mom), ' sign. PNA'));
  726. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  727. SurfStatColLim([-abs_lim, abs_lim]);
  728. sign_vert_data_GA = BoSurfStatMakeParcelData(sign_parcel_data_GA, wk40, wk40_economo_ind);
  729. figure; SurfStatViewData(sign_vert_data_GA, wk40, strcat(moments(mom), ' sign. GA'));
  730. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  731. SurfStatColLim([-abs_lim, abs_lim]);
  732. end
  733. end
  735. end
  737. for figure_5 = 1
  739. for panel_A = 1
  741. % plotting first eigenmode
  742. eigenmode1 = eigenmodes(:,1);
  743. eigenmode1(eigenmode1~=0) = 1;
  744. figure; SurfStatViewData(eigenmode1, wk40);
  745. colormap([1 1 1; 0.1 0.01 0.56])
  747. % plotting vertex-wise eigenmodes 2-4
  748. for i = 1:size(eigenmodes, 2)
  749. figure; SurfStatViewData(eigenmodes(:,i), wk40);
  750. colormap(vik);
  751. end
  753. end
  755. for panel_B = 1
  757. % Note: the brain plots under the scatter plots are the same brain
  758. % plots generated in figures 3 and 4, but without applying significance
  759. % threshold on the t-maps and without a common colourscale. The code
  760. % to generate these is not repeated in this section
  761. n_perm = 10000; % number of permutations
  762. moments = ["Centre of gravity", "Variance"];
  764. % run linear models and concatenate t-values
  765. for mom = 1:length(moments)
  766. if mom == 1
  767. cog_parcel_data_PMA = moment_age_model("cog", "PMA", ...
  768. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  769. cog_parcel_data_GA_PNA = moment_age_model("cog", ["GA", "PNA"], ...
  770. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  771. else
  772. var_parcel_data_PMA = moment_age_model("variance", "PMA", ...
  773. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  774. var_parcel_data_GA_PNA = moment_age_model("variance", ["GA", "PNA"], ...
  775. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  776. end
  777. end
  778. t_maps = {cog_parcel_data_PMA, var_parcel_data_PMA, ...
  779. cog_parcel_data_GA_PNA(:,1), var_parcel_data_GA_PNA(:,1), ...
  780. cog_parcel_data_GA_PNA(:,2), var_parcel_data_GA_PNA(:,2)};
  781. for i = 1:numel(t_maps)
  782. t_maps{i}(t_maps{i} == -100) = NaN;
  783. end
  785. % spin the maps using the precomputed permutation indices
  786. spun_map = zeros(total_num_parcels_economo, n_perm, length(t_maps));
  787. for i = 1:length(t_maps)
  788. spun_map(:,:,i) = t_maps{i}(perm_id);
  789. end
  791. r_emp = zeros(length(t_maps), size(eigenmodes, 2));
  792. r_perm = zeros(n_perm, size(eigenmodes, 2), length(t_maps));
  793. spun_p = zeros(length(t_maps), size(eigenmodes, 2));
  794. null_prctiles = zeros(2, size(eigenmodes, 2), length(t_maps));
  795. for ii = 1:length(t_maps)
  796. spun_map = zeros(size(eigenmodes, 1), n_perm);
  797. % upsample parcel-wise effect maps to vertex-wise maps
  798. vert_t_map = BoSurfStatMakeParcelData(t_maps{ii}, wk40, wk40_economo_ind);
  799. for jj = 1:size(eigenmodes, 2)
  800. % correlate original effect map with eigenmode
  801. r_emp(ii, jj) = corr(vert_t_map', eigenmodes(:, jj), 'Rows', 'complete', 'Type', 'Pearson');
  802. % spin the parcel-wise effect map, upsample to vertex-wise
  803. % maps and correlate with eigenmode
  804. for kk = 1:n_perm
  805. spun_map(:,kk) = BoSurfStatMakeParcelData(t_maps{ii}(perm_id(:,kk)), ...
  806. wk40, wk40_economo_ind);
  807. r_perm(kk, jj, ii) = corr(spun_map(:, kk), eigenmodes(:, jj), ...
  808. 'Rows', 'complete', 'Type', 'Pearson');
  809. end
  810. % calculate two-tailed permutation p-value
  811. spun_p(ii, jj) = sum(abs(r_perm(:, jj, ii)) >= abs(r_emp(ii, jj))) / n_perm;
  812. % 95% null distribution interval (2.5th and 97.5th percentiles)
  813. null_prctiles(:, jj, ii) = prctile(r_perm(:, jj, ii), [2.5, 97.5]);
  814. end
  815. end
  817. % plot correlations with each eigenmode for each age-metric separately
  818. start_x = 30;
  819. side_gap = 15;
  820. gap_within = 5;
  821. gap_between = 20;
  822. % positions to plot values
  823. positions = [start_x, start_x+gap_within, start_x+2*gap_within, ...
  824. start_x+2*gap_within+gap_between, start_x+3*gap_within+gap_between, ...
  825. start_x+4*gap_within+gap_between];
  827. groupped_by = 'moments'; % choose between 'moments' or 'eigenmodes'
  828. orientation = 'vertical'; % choose between 'horizontal' or 'vertical'
  830. for age_idx = 1:3 % {for age = PMA, GA, PNA}
  831. figure;
  832. hold on;
  834. if strcmp(groupped_by, 'moments')
  835. correlation = r_emp(age_idx*2-1:age_idx*2, :);
  836. correlation = reshape(correlation', 1, []);
  837. whiskers = null_prctiles(:, :, age_idx*2-1:age_idx*2);
  838. whiskers = reshape(whiskers, 2, []);
  839. elseif strcmp(groupped_by, 'eigenmodes')
  840. correlation = reshape(correlation, 1, []);
  841. whiskers = permute(whiskers, [1, 3, 2]);
  842. whiskers = reshape(whiskers, 2, 6);
  843. end
  845. col_scaled = round(1 + ((correlation + max(abs(correlation))) ./ ...
  846. (2 * max(abs(correlation)))) * 255);
  848. if strcmp(orientation, 'horizontal')
  849. xline(0, 'k--', 'LineWidth', 1);
  850. scatter(correlation, -positions, 120, vik(col_scaled,:), 'filled');
  851. for i = 1:6
  852. plot([whiskers(1,i), whiskers(2,i)], [-positions(i), -positions(i)], ...
  853. 'Color', 'k', 'LineWidth', 1);
  854. end
  855. xlim([-(max(abs(correlation))+0.1), max(abs(correlation))+0.1]);
  856. ylim([-(max(positions)+side_gap), -(min(positions)-side_gap)]);
  857. set(gca, 'YTick', []);
  858. elseif strcmp(orientation, 'vertical')
  859. yline(0, 'k--', 'LineWidth', 1);
  860. scatter(positions, correlation, 120, vik(col_scaled,:), 'filled');
  861. for i = 1:6
  862. plot([positions(i), positions(i)], [whiskers(1,i), whiskers(2,i)], ...
  863. 'Color', 'k', 'LineWidth', 1);
  864. end
  865. ylim([-(max(abs(correlation))+0.1), max(abs(correlation))+0.1]);
  866. xlim([min(positions)-side_gap, max(positions)+side_gap]);
  867. set(gca, 'XTick', []);
  868. end
  870. set(gca, 'FontSize', 14);
  871. grid off;
  872. hold off;
  873. end
  875. end
  877. for panel_C = 1
  879. % upsample parcel-wise effects maps to vertex-wise maps
  880. vertex_maps = cell(1, numel(t_maps));
  881. for i = 1:numel(t_maps)
  882. vertex_maps{i} = BoSurfStatMakeParcelData(t_maps{i}, ...
  883. wk40, wk40_economo_ind);
  884. end
  885. % concatenate vertex-wise effect maps and eigenmodes in a single table
  886. vertex_matrix = cat(2, vertex_maps{:});
  887. vertex_matrix = reshape(vertex_matrix, [], numel(vertex_maps));
  888. vertex_matrix = [vertex_matrix, eigenmodes(:,1), eigenmodes(:,2), eigenmodes(:,3)];
  889. tmaps_emodes_vert = array2table(vertex_matrix, 'VariableNames', ...
  890. {'cog_PMA','var_PMA', 'cog_GA','var_GA', ...
  891. 'cog_PNA','var_PNA', 'EM1','EM2','EM3'});
  893. response_vars = ["cog_PMA","var_PMA","cog_GA","var_GA","cog_PNA","var_PNA"];
  895. regression_results = zeros(7,3,length(response_vars));
  896. % run linear models with each of the effect maps as response
  897. % variables and all 7 different combinations of eigenmodes as predictors
  898. for var = 1:length(response_vars)
  899. mdl_EM1 = fitlm(tmaps_emodes_vert, strcat(response_vars(var), ' ~ EM1'));
  900. regression_results(1,1,var) = mdl_EM1.Rsquared.Adjusted;
  901. mdl_EM2 = fitlm(tmaps_emodes_vert, strcat(response_vars(var), ' ~ EM2'));
  902. regression_results(2,1,var) = mdl_EM2.Rsquared.Adjusted;
  903. mdl_EM3 = fitlm(tmaps_emodes_vert, strcat(response_vars(var), ' ~ EM3'));
  904. regression_results(3,1,var) = mdl_EM3.Rsquared.Adjusted;
  905. mdl_EM1_2 = fitlm(tmaps_emodes_vert, strcat(response_vars(var), ' ~ EM1 + EM2'));
  906. regression_results(4,1,var) = mdl_EM1_2.Rsquared.Adjusted;
  907. mdl_EM1_3 = fitlm(tmaps_emodes_vert, strcat(response_vars(var), ' ~ EM1 + EM3'));
  908. regression_results(5,1,var) = mdl_EM1_3.Rsquared.Adjusted;
  909. mdl_EM2_3 = fitlm(tmaps_emodes_vert, strcat(response_vars(var), ' ~ EM2 + EM3'));
  910. regression_results(6,1,var) = mdl_EM2_3.Rsquared.Adjusted;
  911. mdl_EM1_2_3 = fitlm(tmaps_emodes_vert, strcat(response_vars(var), ' ~ EM1 + EM2 + EM3'));
  912. regression_results(7,1,var) = mdl_EM1_2_3.Rsquared.Adjusted;
  913. end
  915. r_sqrd_vals = squeeze(regression_results(:,1,:))';
  917. tmap_names = ["cog-PMA","var-PMA","cog-GA","var-GA","cog-PNA","var-PNA"];
  918. model_names = ["EM2","EM3","EM4","EM2+EM3","EM2+EM4","EM3+EM4","EM2+EM3+EM4"];
  920. % plot the adjusted R-square values of each model in a bar plot
  921. cols = [vik(145,:); vik(195,:)];
  922. gap_between = 1;
  923. block_width = size(r_sqrd_vals,2) + gap_between;
  924. for age = 1:3
  925. figure('units','centimeters','outerposition',[0 0 20 12]); hold on;
  926. data_to_plot = r_sqrd_vals((age*2-1):age*2,:);
  927. for mom = 1:2
  928. x_pos = (mom-1)*block_width + (1:size(r_sqrd_vals,2));
  929. bar(x_pos,data_to_plot(mom,:), 0.8, 'FaceColor', cols(mom,:), ...
  930. 'EdgeColor', 'k', 'LineWidth', 1);
  931. end
  932. set(gca, 'FontSize', 14);
  933. ylim([0 max(r_sqrd_vals(:))*1.1]);
  934. box on; grid off
  935. end
  937. end
  939. end
  941. for supp_figure_1 = 1
  943. moments = ["Centre of gravity", "Variance"];
  944. % initialise table to store the correlations between the effect maps of
  945. % the original dataset and the dataset excluding twins
  946. r_emp = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  947. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  949. for mom = 1:length(moments)
  950. if mom == 1
  951. cog_parcel_data_PMA = moment_age_model("cog", "PMA", ...
  952. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  953. cog_parcel_data_PMA_no_twins = moment_age_model("cog", "PMA", ...
  954. gen_table_no_twins, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  955. else
  956. var_parcel_data_PMA = moment_age_model("variance", "PMA", ...
  957. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  958. var_parcel_data_PMA_no_twins = moment_age_model("variance", "PMA", ...
  959. gen_table_no_twins, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  960. end
  961. end
  963. % calculate maximum absolute value for limits
  964. min_val = min([min(cog_parcel_data_PMA(cog_parcel_data_PMA > -100)); ...
  965. min(cog_parcel_data_PMA_no_twins(cog_parcel_data_PMA_no_twins > -100));
  966. min(var_parcel_data_PMA(var_parcel_data_PMA > -100)); ...
  967. min(var_parcel_data_PMA_no_twins(var_parcel_data_PMA_no_twins > -100))]);
  968. max_val = max([max(cog_parcel_data_PMA(cog_parcel_data_PMA < 100)); ...
  969. max(cog_parcel_data_PMA_no_twins(cog_parcel_data_PMA_no_twins < 100)); ...
  970. max(var_parcel_data_PMA(var_parcel_data_PMA < 100)); ...
  971. max(var_parcel_data_PMA_no_twins(var_parcel_data_PMA_no_twins < 100))]);
  972. abs_lim = max(abs(min_val), abs(max_val));
  973. abs_lim = abs_lim + abs_lim/100;
  975. % plot the significant effects of PMA on each moment
  976. for mom = 1:length(moments)
  977. if mom == 1
  978. parc_data_PMA = cog_parcel_data_PMA;
  979. parc_data_PMA_no_twins = cog_parcel_data_PMA_no_twins;
  980. else
  981. parc_data_PMA = var_parcel_data_PMA;
  982. parc_data_PMA_no_twins = var_parcel_data_PMA_no_twins;
  983. end
  985. vert_data_PMA = BoSurfStatMakeParcelData(parc_data_PMA, wk40, wk40_economo_ind);
  986. figure; SurfStatViewData(vert_data_PMA, wk40, strcat(moments(mom), ' PMA'));
  987. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  988. SurfStatColLim([-abs_lim, abs_lim]);
  989. vert_data_PMA_no_twins = BoSurfStatMakeParcelData(parc_data_PMA_no_twins, wk40, wk40_economo_ind);
  990. figure; SurfStatViewData(vert_data_PMA_no_twins, wk40, strcat(moments(mom), ' PMA - no twins'));
  991. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  992. SurfStatColLim([-abs_lim, abs_lim]);
  994. % replace -100 values with NaNs for the correlation
  995. parc_data_PMA(parc_data_PMA==-100)=NaN;
  996. parc_data_PMA_no_twins(parc_data_PMA_no_twins==-100)=NaN;
  997. spun_map = parc_data_PMA(perm_id); % spin one of the maps of effects
  998. r_emp.Correlation(mom)=corr(parc_data_PMA, parc_data_PMA_no_twins, ...
  999. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the effects
  1000. % run permutation test
  1001. r_perm = zeros(n_perm,1);
  1002. for kk=1:n_perm
  1003. r_perm(kk) = corr(spun_map(:, kk), parc_data_PMA_no_twins, ...
  1004. 'Rows', 'complete', 'Type', 'Pearson');
  1005. end
  1006. r_emp.p_value(mom) = sum(abs(r_perm(:)) >= abs(r_emp.Correlation(mom))) / n_perm; % p-value after spin test
  1007. end
  1009. end
  1011. for supp_figures_2_3 = 1
  1013. moments = ["Centre of gravity", "Variance"];
  1014. % initialise tables to store the correlations between the effect maps of
  1015. % the original dataset and the dataset excluding twins
  1016. r_emp_GA = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1017. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1018. r_emp_PNA = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1019. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1021. % run models which contain both GA and PNA as predictors (effects
  1022. % controlling for the effects of the other)
  1023. for mom = 1:length(moments)
  1024. if mom == 1
  1025. cog_parcel_data = moment_age_model("cog", ["GA", "PNA"], ...
  1026. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1027. cog_parcel_data_no_twins = moment_age_model("cog", ["GA", "PNA"], ...
  1028. gen_table_no_twins, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1029. else
  1030. var_parcel_data = moment_age_model("variance", ["GA", "PNA"], ...
  1031. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1032. var_parcel_data_no_twins = moment_age_model("variance", ["GA", "PNA"], ...
  1033. gen_table_no_twins, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1034. end
  1035. end
  1037. % calculate maximum absolute value for limits for each age-metric separately
  1038. min_val_GA = min([min(cog_parcel_data(cog_parcel_data(:,1) > -100, 1)); ...
  1039. min(cog_parcel_data_no_twins(cog_parcel_data_no_twins(:,1) > -100, 1)); ...
  1040. min(var_parcel_data(var_parcel_data(:,1) > -100, 1)); ...
  1041. min(var_parcel_data_no_twins(var_parcel_data_no_twins(:,1) > -100, 1))]);
  1042. max_val_GA = max([max(cog_parcel_data(cog_parcel_data(:,1) < 100, 1)); ...
  1043. max(cog_parcel_data_no_twins(cog_parcel_data_no_twins(:,1) < 100, 1)); ...
  1044. max(var_parcel_data(var_parcel_data(:,1) < 100, 1)); ...
  1045. max(var_parcel_data_no_twins(var_parcel_data_no_twins(:,1) < 100, 1))]);
  1046. abs_lim_GA = max([abs(min_val_GA), abs(max_val_GA)]);
  1047. abs_lim_GA = abs_lim_GA + abs_lim_GA/100;
  1049. min_val_PNA = min([min(cog_parcel_data(cog_parcel_data(:,2) > -100, 2)); ...
  1050. min(cog_parcel_data_no_twins(cog_parcel_data_no_twins(:,2) > -100, 2)); ...
  1051. min(var_parcel_data(var_parcel_data(:,2) > -100, 2)); ...
  1052. min(var_parcel_data_no_twins(var_parcel_data_no_twins(:,2) > -100, 2))]);
  1053. max_val_PNA = max([max(cog_parcel_data(cog_parcel_data(:,2) < 100, 2)); ...
  1054. max(cog_parcel_data_no_twins(cog_parcel_data_no_twins(:,2) < 100, 2)); ...
  1055. max(var_parcel_data(var_parcel_data(:,2) < 100, 2)); ...
  1056. max(var_parcel_data_no_twins(var_parcel_data_no_twins(:,2) < 100, 2))]);
  1057. abs_lim_PNA = max([abs(min_val_PNA), abs(max_val_PNA)]);
  1058. abs_lim_PNA = abs_lim_PNA + abs_lim_PNA/100;
  1060. % plot the significant effects of GA and PNA on each moment
  1061. for mom = 1:length(moments)
  1062. if mom == 1
  1063. parcel_data_GA = cog_parcel_data(:,1);
  1064. parcel_data_GA_no_twins = cog_parcel_data_no_twins(:,1);
  1065. parcel_data_PNA = cog_parcel_data(:,2);
  1066. parcel_data_PNA_no_twins = cog_parcel_data_no_twins(:,2);
  1067. else
  1068. parcel_data_GA = var_parcel_data(:,1);
  1069. parcel_data_GA_no_twins = var_parcel_data_no_twins(:,1);
  1070. parcel_data_PNA = var_parcel_data(:,2);
  1071. parcel_data_PNA_no_twins = var_parcel_data_no_twins(:,2);
  1072. end
  1074. vert_data_PNA = BoSurfStatMakeParcelData(parcel_data_PNA, wk40, wk40_economo_ind);
  1075. figure; SurfStatViewData(vert_data_PNA, wk40, strcat(moments(mom), ' PNA'));
  1076. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1077. SurfStatColLim([-abs_lim_PNA, abs_lim_PNA]);
  1078. vert_data_PNA_no_twins = BoSurfStatMakeParcelData(parcel_data_PNA_no_twins, wk40, wk40_economo_ind);
  1079. figure; SurfStatViewData(vert_data_PNA_no_twins, wk40, strcat(moments(mom), ' PNA - no twins'));
  1080. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1081. SurfStatColLim([-abs_lim_PNA, abs_lim_PNA]);
  1082. vert_data_GA = BoSurfStatMakeParcelData(parcel_data_GA, wk40, wk40_economo_ind);
  1083. figure; SurfStatViewData(vert_data_GA, wk40, strcat(moments(mom), ' GA'));
  1084. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1085. SurfStatColLim([-abs_lim_GA, abs_lim_GA]);
  1086. vert_data_GA_no_twins = BoSurfStatMakeParcelData(parcel_data_GA_no_twins, wk40, wk40_economo_ind);
  1087. figure; SurfStatViewData(vert_data_GA_no_twins, wk40, strcat(moments(mom), ' GA - no twins'));
  1088. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1089. SurfStatColLim([-abs_lim_GA, abs_lim_GA]);
  1091. % replace -100 values with NaNs for the correlation
  1092. parcel_data_GA(parcel_data_GA==-100)=NaN;
  1093. parcel_data_GA_no_twins(parcel_data_GA_no_twins==-100)=NaN;
  1094. parcel_data_PNA(parcel_data_PNA==-100)=NaN;
  1095. parcel_data_PNA_no_twins(parcel_data_PNA_no_twins==-100)=NaN;
  1096. spun_map_GA = parcel_data_GA(perm_id); % spin one of the GA maps of effects
  1097. spun_map_PNA = parcel_data_PNA(perm_id); % spin one of the PNA maps of effects
  1098. r_emp_GA.Correlation(mom)=corr(parcel_data_GA, parcel_data_GA_no_twins, ...
  1099. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the GA effects
  1100. r_emp_PNA.Correlation(mom)=corr(parcel_data_PNA, parcel_data_PNA_no_twins, ...
  1101. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the PNA effects
  1102. % run permutation test
  1103. r_perm_GA = zeros(n_perm,1);
  1104. r_perm_PNA = zeros(n_perm,1);
  1105. for kk=1:n_perm
  1106. r_perm_GA(kk) = corr(spun_map_GA(:, kk), parcel_data_GA_no_twins, ...
  1107. 'Rows', 'complete', 'Type', 'Pearson');
  1108. r_perm_PNA(kk) = corr(spun_map_PNA(:, kk), parcel_data_PNA_no_twins, ...
  1109. 'Rows', 'complete', 'Type', 'Pearson');
  1110. end
  1111. r_emp_GA.p_value(mom) = sum(abs(r_perm_GA(:)) ...
  1112. >= abs(r_emp_GA.Correlation(mom))) / n_perm; % p-value after spin test
  1113. r_emp_PNA.p_value(mom) = sum(abs(r_perm_PNA(:)) ...
  1114. >= abs(r_emp_PNA.Correlation(mom))) / n_perm; % p-value after spin test
  1115. end
  1117. end
  1119. for supp_figure_4 = 1
  1121. moments = ["Centre of gravity", "Variance"];
  1122. % initialise table to store the correlations between the effect maps of the
  1123. % original dataset and the dataset excluding extremely preterm (EPT) infants
  1124. r_emp_PMA = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1125. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1126. r_emp_GA = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1127. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1128. r_emp_PNA = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1129. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1131. for mom = 1:length(moments)
  1132. if mom == 1
  1133. cog_parcel_data_PMA = moment_age_model("cog", "PMA", ...
  1134. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1135. cog_parcel_data_PMA_no_EPT = moment_age_model("cog", "PMA", ...
  1136. gen_table_economo(gen_table_economo.GA >= 28, :), uparc_economo, ...
  1137. total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1138. cog_parcel_data_GA_PNA = moment_age_model("cog", ["GA", "PNA"], ...
  1139. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1140. cog_parcel_data_GA_PNA_no_EPT = moment_age_model("cog", ["GA", "PNA"], ...
  1141. gen_table_economo(gen_table_economo.GA >= 28, :), uparc_economo, ...
  1142. total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1143. else
  1144. var_parcel_data_PMA = moment_age_model("variance", "PMA", ...
  1145. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1146. var_parcel_data_PMA_no_EPT = moment_age_model("variance", "PMA", ...
  1147. gen_table_economo(gen_table_economo.GA >= 28, :), uparc_economo, ...
  1148. total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1149. var_parcel_data_GA_PNA = moment_age_model("variance", ["GA", "PNA"], ...
  1150. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1151. var_parcel_data_GA_PNA_no_EPT = moment_age_model("variance", ["GA", "PNA"], ...
  1152. gen_table_economo(gen_table_economo.GA >= 28, :), uparc_economo, ...
  1153. total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1154. end
  1155. end
  1157. % calculate maximum absolute value for limits
  1158. min_val = min([min(cog_parcel_data_PMA(cog_parcel_data_PMA > -100)); ...
  1159. min(cog_parcel_data_PMA_no_EPT(cog_parcel_data_PMA_no_EPT > -100)); ...
  1160. min(cog_parcel_data_GA_PNA(cog_parcel_data_GA_PNA > -100)); ...
  1161. min(cog_parcel_data_GA_PNA_no_EPT(cog_parcel_data_GA_PNA_no_EPT > -100)); ...
  1162. min(var_parcel_data_PMA(var_parcel_data_PMA > -100)); ...
  1163. min(var_parcel_data_PMA_no_EPT(var_parcel_data_PMA_no_EPT > -100)); ...
  1164. min(var_parcel_data_GA_PNA(var_parcel_data_GA_PNA > -100)); ...
  1165. min(var_parcel_data_GA_PNA_no_EPT(var_parcel_data_GA_PNA_no_EPT > -100))]);
  1166. max_val = max([max(cog_parcel_data_PMA(cog_parcel_data_PMA < 100)); ...
  1167. max(cog_parcel_data_PMA_no_EPT(cog_parcel_data_PMA_no_EPT < 100)); ...
  1168. max(cog_parcel_data_GA_PNA(cog_parcel_data_GA_PNA < 100)); ...
  1169. max(cog_parcel_data_GA_PNA_no_EPT(cog_parcel_data_GA_PNA_no_EPT < 100)); ...
  1170. max(var_parcel_data_PMA(var_parcel_data_PMA < 100)); ...
  1171. max(var_parcel_data_PMA_no_EPT(var_parcel_data_PMA_no_EPT < 100)); ...
  1172. max(var_parcel_data_GA_PNA(var_parcel_data_GA_PNA < 100)); ...
  1173. max(var_parcel_data_GA_PNA_no_EPT(var_parcel_data_GA_PNA_no_EPT < 100))]);
  1174. abs_lim = max(abs(min_val), abs(max_val));
  1175. abs_lim = abs_lim + abs_lim/100;
  1177. % plot the significant effects of PMA on each moment
  1178. for mom = 1:length(moments)
  1179. if mom == 1
  1180. parcel_data_PMA = cog_parcel_data_PMA;
  1181. parcel_data_PMA_no_EPT = cog_parcel_data_PMA_no_EPT;
  1182. parcel_data_GA = cog_parcel_data_GA_PNA(:,1);
  1183. parcel_data_GA_no_EPT = cog_parcel_data_GA_PNA_no_EPT(:,1);
  1184. parcel_data_PNA = cog_parcel_data_GA_PNA(:,2);
  1185. parcel_data_PNA_no_EPT = cog_parcel_data_GA_PNA_no_EPT(:,2);
  1186. else
  1187. parcel_data_PMA = var_parcel_data_PMA;
  1188. parcel_data_PMA_no_EPT = var_parcel_data_PMA_no_EPT;
  1189. parcel_data_GA = var_parcel_data_GA_PNA(:,1);
  1190. parcel_data_GA_no_EPT = var_parcel_data_GA_PNA_no_EPT(:,1);
  1191. parcel_data_PNA = var_parcel_data_GA_PNA(:,2);
  1192. parcel_data_PNA_no_EPT = var_parcel_data_GA_PNA_no_EPT(:,2);
  1193. end
  1195. vert_data_PMA = BoSurfStatMakeParcelData(parcel_data_PMA, wk40, wk40_economo_ind);
  1196. figure; SurfStatViewData(vert_data_PMA, wk40, strcat(moments(mom), ' PMA'));
  1197. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1198. SurfStatColLim([-abs_lim, abs_lim]);
  1199. vert_data_PMA_no_EPT = BoSurfStatMakeParcelData(parcel_data_PMA_no_EPT, wk40, wk40_economo_ind);
  1200. figure; SurfStatViewData(vert_data_PMA_no_EPT, wk40, strcat(moments(mom), ' PMA - no EPT'));
  1201. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1202. SurfStatColLim([-abs_lim, abs_lim]);
  1203. vert_data_GA = BoSurfStatMakeParcelData(parcel_data_GA, wk40, wk40_economo_ind);
  1204. figure; SurfStatViewData(vert_data_GA, wk40, strcat(moments(mom), ' GA'));
  1205. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1206. SurfStatColLim([-abs_lim, abs_lim]);
  1207. vert_data_GA_no_EPT = BoSurfStatMakeParcelData(parcel_data_GA_no_EPT, wk40, wk40_economo_ind);
  1208. figure; SurfStatViewData(vert_data_GA_no_EPT, wk40, strcat(moments(mom), ' GA - no EPT'));
  1209. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1210. SurfStatColLim([-abs_lim, abs_lim]);
  1211. vert_data_PNA = BoSurfStatMakeParcelData(parcel_data_PNA, wk40, wk40_economo_ind);
  1212. figure; SurfStatViewData(vert_data_PNA, wk40, strcat(moments(mom), ' PNA'));
  1213. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1214. SurfStatColLim([-abs_lim, abs_lim]);
  1215. vert_data_PNA_no_EPT = BoSurfStatMakeParcelData(parcel_data_PNA_no_EPT, wk40, wk40_economo_ind);
  1216. figure; SurfStatViewData(vert_data_PNA_no_EPT, wk40, strcat(moments(mom), ' PNA - no EPT'));
  1217. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1218. SurfStatColLim([-abs_lim, abs_lim]);
  1220. % replace -100 values with NaNs for the correlation
  1221. parcel_data_PMA(parcel_data_PMA==-100)=NaN;
  1222. parcel_data_PMA_no_EPT(parcel_data_PMA_no_EPT==-100)=NaN;
  1223. parcel_data_GA(parcel_data_GA==-100)=NaN;
  1224. parcel_data_GA_no_EPT(parcel_data_GA_no_EPT==-100)=NaN;
  1225. parcel_data_PNA(parcel_data_PNA==-100)=NaN;
  1226. parcel_data_PNA_no_EPT(parcel_data_PNA_no_EPT==-100)=NaN;
  1227. spun_map_PMA = parcel_data_PMA(perm_id); % spin one of the PMA maps of effects
  1228. spun_map_GA = parcel_data_GA(perm_id); % spin one of the GA maps of effects
  1229. spun_map_PNA = parcel_data_PNA(perm_id); % spin one of the PNA maps of effects
  1230. r_emp_PMA.Correlation(mom)=corr(parcel_data_PMA, parcel_data_PMA_no_EPT, ...
  1231. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the PMA effects
  1232. r_emp_GA.Correlation(mom)=corr(parcel_data_GA, parcel_data_GA_no_EPT, ...
  1233. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the GA effects
  1234. r_emp_PNA.Correlation(mom)=corr(parcel_data_PNA, parcel_data_PNA_no_EPT, ...
  1235. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the PNA effects
  1236. % run permutation test
  1237. r_perm_PMA = zeros(n_perm,1);
  1238. r_perm_GA = zeros(n_perm,1);
  1239. r_perm_PNA = zeros(n_perm,1);
  1240. for kk=1:n_perm
  1241. r_perm_PMA(kk) = corr(spun_map_PMA(:, kk), parcel_data_PMA_no_EPT, ...
  1242. 'Rows', 'complete', 'Type', 'Pearson');
  1243. r_perm_GA(kk) = corr(spun_map_GA(:, kk), parcel_data_GA_no_EPT, ...
  1244. 'Rows', 'complete', 'Type', 'Pearson');
  1245. r_perm_PNA(kk) = corr(spun_map_PNA(:, kk), parcel_data_PNA_no_EPT, ...
  1246. 'Rows', 'complete', 'Type', 'Pearson');
  1247. end
  1248. r_emp_PMA.p_value(mom) = sum(abs(r_perm(:)) ...
  1249. >= abs(r_emp_PMA.Correlation(mom))) / n_perm; % p-value after spin test
  1250. r_emp_GA.p_value(mom) = sum(abs(r_perm(:)) ...
  1251. >= abs(r_emp_GA.Correlation(mom))) / n_perm; % p-value after spin test
  1252. r_emp_PNA.p_value(mom) = sum(abs(r_perm(:)) ...
  1253. >= abs(r_emp_PNA.Correlation(mom))) / n_perm; % p-value after spin test
  1254. end
  1256. end
  1258. for supp_figure_5 = 1
  1260. moments = ["Centre of gravity", "Variance"];
  1261. % initialise table to store the correlations between the effect maps of the
  1262. % original dataset and the dataset excluding very preterm (VPT) infants
  1263. r_emp_PMA = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1264. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1265. r_emp_GA = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1266. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1267. r_emp_PNA = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1268. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1270. for mom = 1:length(moments)
  1271. if mom == 1
  1272. cog_parcel_data_PMA = moment_age_model("cog", "PMA", ...
  1273. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1274. cog_parcel_data_PMA_no_VPT = moment_age_model("cog", "PMA", ...
  1275. gen_table_economo(gen_table_economo.GA >= 32, :), uparc_economo, ...
  1276. total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1277. cog_parcel_data_GA_PNA = moment_age_model("cog", ["GA", "PNA"], ...
  1278. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1279. cog_parcel_data_GA_PNA_no_VPT = moment_age_model("cog", ["GA", "PNA"], ...
  1280. gen_table_economo(gen_table_economo.GA >= 32, :), uparc_economo, ...
  1281. total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1282. else
  1283. var_parcel_data_PMA = moment_age_model("variance", "PMA", ...
  1284. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1285. var_parcel_data_PMA_no_VPT = moment_age_model("variance", "PMA", ...
  1286. gen_table_economo(gen_table_economo.GA >= 32, :), uparc_economo, ...
  1287. total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1288. var_parcel_data_GA_PNA = moment_age_model("variance", ["GA", "PNA"], ...
  1289. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1290. var_parcel_data_GA_PNA_no_VPT = moment_age_model("variance", ["GA", "PNA"], ...
  1291. gen_table_economo(gen_table_economo.GA >= 32, :), uparc_economo, ...
  1292. total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1293. end
  1294. end
  1296. % calculate maximum absolute value for limits
  1297. min_val = min([min(cog_parcel_data_PMA(cog_parcel_data_PMA > -100)); ...
  1298. min(cog_parcel_data_PMA_no_VPT(cog_parcel_data_PMA_no_VPT > -100)); ...
  1299. min(cog_parcel_data_GA_PNA(cog_parcel_data_GA_PNA > -100)); ...
  1300. min(cog_parcel_data_GA_PNA_no_VPT(cog_parcel_data_GA_PNA_no_VPT > -100)); ...
  1301. min(var_parcel_data_PMA(var_parcel_data_PMA > -100)); ...
  1302. min(var_parcel_data_PMA_no_VPT(var_parcel_data_PMA_no_VPT > -100)); ...
  1303. min(var_parcel_data_GA_PNA(var_parcel_data_GA_PNA > -100)); ...
  1304. min(var_parcel_data_GA_PNA_no_VPT(var_parcel_data_GA_PNA_no_VPT > -100))]);
  1305. max_val = max([max(cog_parcel_data_PMA(cog_parcel_data_PMA < 100)); ...
  1306. max(cog_parcel_data_PMA_no_VPT(cog_parcel_data_PMA_no_VPT < 100)); ...
  1307. max(cog_parcel_data_GA_PNA(cog_parcel_data_GA_PNA < 100)); ...
  1308. max(cog_parcel_data_GA_PNA_no_VPT(cog_parcel_data_GA_PNA_no_VPT < 100)); ...
  1309. max(var_parcel_data_PMA(var_parcel_data_PMA < 100)); ...
  1310. max(var_parcel_data_PMA_no_VPT(var_parcel_data_PMA_no_VPT < 100)); ...
  1311. max(var_parcel_data_GA_PNA(var_parcel_data_GA_PNA < 100)); ...
  1312. max(var_parcel_data_GA_PNA_no_VPT(var_parcel_data_GA_PNA_no_VPT < 100))]);
  1313. abs_lim = max(abs(min_val), abs(max_val));
  1314. abs_lim = abs_lim + abs_lim/100;
  1316. % plot the significant effects of PMA on each moment
  1317. for mom = 1:length(moments)
  1318. if mom == 1
  1319. parcel_data_PMA = cog_parcel_data_PMA;
  1320. parcel_data_PMA_no_VPT = cog_parcel_data_PMA_no_VPT;
  1321. parcel_data_GA = cog_parcel_data_GA_PNA(:,1);
  1322. parcel_data_GA_no_VPT = cog_parcel_data_GA_PNA_no_VPT(:,1);
  1323. parcel_data_PNA = cog_parcel_data_GA_PNA(:,2);
  1324. parcel_data_PNA_no_VPT = cog_parcel_data_GA_PNA_no_VPT(:,2);
  1325. else
  1326. parcel_data_PMA = var_parcel_data_PMA;
  1327. parcel_data_PMA_no_VPT = var_parcel_data_PMA_no_VPT;
  1328. parcel_data_GA = var_parcel_data_GA_PNA(:,1);
  1329. parcel_data_GA_no_VPT = var_parcel_data_GA_PNA_no_VPT(:,1);
  1330. parcel_data_PNA = var_parcel_data_GA_PNA(:,2);
  1331. parcel_data_PNA_no_VPT = var_parcel_data_GA_PNA_no_VPT(:,2);
  1332. end
  1334. vert_data_PMA = BoSurfStatMakeParcelData(parcel_data_PMA, wk40, wk40_economo_ind);
  1335. figure; SurfStatViewData(vert_data_PMA, wk40, strcat(moments(mom), ' PMA'));
  1336. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1337. SurfStatColLim([-abs_lim, abs_lim]);
  1338. vert_data_PMA_no_VPT = BoSurfStatMakeParcelData(parcel_data_PMA_no_VPT, wk40, wk40_economo_ind);
  1339. figure; SurfStatViewData(vert_data_PMA_no_VPT, wk40, strcat(moments(mom), ' PMA - no VPT'));
  1340. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1341. SurfStatColLim([-abs_lim, abs_lim]);
  1342. vert_data_GA = BoSurfStatMakeParcelData(parcel_data_GA, wk40, wk40_economo_ind);
  1343. figure; SurfStatViewData(vert_data_GA, wk40, strcat(moments(mom), ' GA'));
  1344. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1345. SurfStatColLim([-abs_lim, abs_lim]);
  1346. vert_data_GA_no_VPT = BoSurfStatMakeParcelData(parcel_data_GA_no_VPT, wk40, wk40_economo_ind);
  1347. figure; SurfStatViewData(vert_data_GA_no_VPT, wk40, strcat(moments(mom), ' GA - no VPT'));
  1348. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1349. SurfStatColLim([-abs_lim, abs_lim]);
  1350. vert_data_PNA = BoSurfStatMakeParcelData(parcel_data_PNA, wk40, wk40_economo_ind);
  1351. figure; SurfStatViewData(vert_data_PNA, wk40, strcat(moments(mom), ' PNA'));
  1352. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1353. SurfStatColLim([-abs_lim, abs_lim]);
  1354. vert_data_PNA_no_VPT = BoSurfStatMakeParcelData(parcel_data_PNA_no_VPT, wk40, wk40_economo_ind);
  1355. figure; SurfStatViewData(vert_data_PNA_no_VPT, wk40, strcat(moments(mom), ' PNA - no VPT'));
  1356. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1357. SurfStatColLim([-abs_lim, abs_lim]);
  1359. % replace -100 values with NaNs for the correlation
  1360. parcel_data_PMA(parcel_data_PMA==-100)=NaN;
  1361. parcel_data_PMA_no_VPT(parcel_data_PMA_no_VPT==-100)=NaN;
  1362. parcel_data_GA(parcel_data_GA==-100)=NaN;
  1363. parcel_data_GA_no_VPT(parcel_data_GA_no_VPT==-100)=NaN;
  1364. parcel_data_PNA(parcel_data_PNA==-100)=NaN;
  1365. parcel_data_PNA_no_VPT(parcel_data_PNA_no_VPT==-100)=NaN;
  1366. spun_map_PMA = parcel_data_PMA(perm_id); % spin one of the PMA maps of effects
  1367. spun_map_GA = parcel_data_GA(perm_id); % spin one of the GA maps of effects
  1368. spun_map_PNA = parcel_data_PNA(perm_id); % spin one of the PNA maps of effects
  1369. r_emp_PMA.Correlation(mom)=corr(parcel_data_PMA, parcel_data_PMA_no_VPT, ...
  1370. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the PMA effects
  1371. r_emp_GA.Correlation(mom)=corr(parcel_data_GA, parcel_data_GA_no_VPT, ...
  1372. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the GA effects
  1373. r_emp_PNA.Correlation(mom)=corr(parcel_data_PNA, parcel_data_PNA_no_VPT, ...
  1374. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the PNA effects
  1375. % run permutation test
  1376. r_perm_PMA = zeros(n_perm,1);
  1377. r_perm_GA = zeros(n_perm,1);
  1378. r_perm_PNA = zeros(n_perm,1);
  1379. for kk=1:n_perm
  1380. r_perm_PMA(kk) = corr(spun_map_PMA(:, kk), parcel_data_PMA_no_VPT, ...
  1381. 'Rows', 'complete', 'Type', 'Pearson');
  1382. r_perm_GA(kk) = corr(spun_map_GA(:, kk), parcel_data_GA_no_VPT, ...
  1383. 'Rows', 'complete', 'Type', 'Pearson');
  1384. r_perm_PNA(kk) = corr(spun_map_PNA(:, kk), parcel_data_PNA_no_VPT, ...
  1385. 'Rows', 'complete', 'Type', 'Pearson');
  1386. end
  1387. r_emp_PMA.p_value(mom) = sum(abs(r_perm(:)) ...
  1388. >= abs(r_emp_PMA.Correlation(mom))) / n_perm; % p-value after spin test
  1389. r_emp_GA.p_value(mom) = sum(abs(r_perm(:)) ...
  1390. >= abs(r_emp_GA.Correlation(mom))) / n_perm; % p-value after spin test
  1391. r_emp_PNA.p_value(mom) = sum(abs(r_perm(:)) ...
  1392. >= abs(r_emp_PNA.Correlation(mom))) / n_perm; % p-value after spin test
  1393. end
  1395. end
  1397. for supp_figure_6 = 1
  1399. moments = ["Centre of gravity", "Variance"];
  1400. % initialise table to store the correlations between the effect maps of the
  1401. % original dataset and the dataset excluding preterm (PT) infants
  1402. r_emp_PMA = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1403. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1404. r_emp_GA = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1405. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1406. r_emp_PNA = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1407. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1409. for mom = 1:length(moments)
  1410. if mom == 1
  1411. cog_parcel_data_PMA = moment_age_model("cog", "PMA", ...
  1412. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1413. cog_parcel_data_PMA_no_PT = moment_age_model("cog", "PMA", ...
  1414. gen_table_economo(gen_table_economo.GA >= 37, :), uparc_economo, ...
  1415. total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1416. cog_parcel_data_GA_PNA = moment_age_model("cog", ["GA", "PNA"], ...
  1417. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1418. cog_parcel_data_GA_PNA_no_PT = moment_age_model("cog", ["GA", "PNA"], ...
  1419. gen_table_economo(gen_table_economo.GA >= 37, :), uparc_economo, ...
  1420. total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1421. else
  1422. var_parcel_data_PMA = moment_age_model("variance", "PMA", ...
  1423. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1424. var_parcel_data_PMA_no_PT = moment_age_model("variance", "PMA", ...
  1425. gen_table_economo(gen_table_economo.GA >= 37, :), uparc_economo, ...
  1426. total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1427. var_parcel_data_GA_PNA = moment_age_model("variance", ["GA", "PNA"], ...
  1428. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1429. var_parcel_data_GA_PNA_no_PT = moment_age_model("variance", ["GA", "PNA"], ...
  1430. gen_table_economo(gen_table_economo.GA >= 37, :), uparc_economo, ...
  1431. total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1432. end
  1433. end
  1435. % calculate maximum absolute value for limits
  1436. min_val = min([min(cog_parcel_data_PMA(cog_parcel_data_PMA > -100)); ...
  1437. min(cog_parcel_data_PMA_no_PT(cog_parcel_data_PMA_no_PT > -100)); ...
  1438. min(cog_parcel_data_GA_PNA(cog_parcel_data_GA_PNA > -100)); ...
  1439. min(cog_parcel_data_GA_PNA_no_PT(cog_parcel_data_GA_PNA_no_PT > -100)); ...
  1440. min(var_parcel_data_PMA(var_parcel_data_PMA > -100)); ...
  1441. min(var_parcel_data_PMA_no_PT(var_parcel_data_PMA_no_PT > -100)); ...
  1442. min(var_parcel_data_GA_PNA(var_parcel_data_GA_PNA > -100)); ...
  1443. min(var_parcel_data_GA_PNA_no_PT(var_parcel_data_GA_PNA_no_PT > -100))]);
  1444. max_val = max([max(cog_parcel_data_PMA(cog_parcel_data_PMA < 100)); ...
  1445. max(cog_parcel_data_PMA_no_PT(cog_parcel_data_PMA_no_PT < 100)); ...
  1446. max(cog_parcel_data_GA_PNA(cog_parcel_data_GA_PNA < 100)); ...
  1447. max(cog_parcel_data_GA_PNA_no_PT(cog_parcel_data_GA_PNA_no_PT < 100)); ...
  1448. max(var_parcel_data_PMA(var_parcel_data_PMA < 100)); ...
  1449. max(var_parcel_data_PMA_no_PT(var_parcel_data_PMA_no_PT < 100)); ...
  1450. max(var_parcel_data_GA_PNA(var_parcel_data_GA_PNA < 100)); ...
  1451. max(var_parcel_data_GA_PNA_no_PT(var_parcel_data_GA_PNA_no_PT < 100))]);
  1452. abs_lim = max(abs(min_val), abs(max_val));
  1453. abs_lim = abs_lim + abs_lim/100;
  1455. % plot the significant effects of PMA on each moment
  1456. for mom = 1:length(moments)
  1457. if mom == 1
  1458. parcel_data_PMA = cog_parcel_data_PMA;
  1459. parcel_data_PMA_no_PT = cog_parcel_data_PMA_no_PT;
  1460. parcel_data_GA = cog_parcel_data_GA_PNA(:,1);
  1461. parcel_data_GA_no_PT = cog_parcel_data_GA_PNA_no_PT(:,1);
  1462. parcel_data_PNA = cog_parcel_data_GA_PNA(:,2);
  1463. parcel_data_PNA_no_PT = cog_parcel_data_GA_PNA_no_PT(:,2);
  1464. else
  1465. parcel_data_PMA = var_parcel_data_PMA;
  1466. parcel_data_PMA_no_PT = var_parcel_data_PMA_no_PT;
  1467. parcel_data_GA = var_parcel_data_GA_PNA(:,1);
  1468. parcel_data_GA_no_PT = var_parcel_data_GA_PNA_no_PT(:,1);
  1469. parcel_data_PNA = var_parcel_data_GA_PNA(:,2);
  1470. parcel_data_PNA_no_PT = var_parcel_data_GA_PNA_no_PT(:,2);
  1471. end
  1473. vert_data_PMA = BoSurfStatMakeParcelData(parcel_data_PMA, wk40, wk40_economo_ind);
  1474. figure; SurfStatViewData(vert_data_PMA, wk40, strcat(moments(mom), ' PMA'));
  1475. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1476. SurfStatColLim([-abs_lim, abs_lim]);
  1477. vert_data_PMA_no_PT = BoSurfStatMakeParcelData(parcel_data_PMA_no_PT, wk40, wk40_economo_ind);
  1478. figure; SurfStatViewData(vert_data_PMA_no_PT, wk40, strcat(moments(mom), ' PMA - no PT'));
  1479. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1480. SurfStatColLim([-abs_lim, abs_lim]);
  1481. vert_data_GA = BoSurfStatMakeParcelData(parcel_data_GA, wk40, wk40_economo_ind);
  1482. figure; SurfStatViewData(vert_data_GA, wk40, strcat(moments(mom), ' GA'));
  1483. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1484. SurfStatColLim([-abs_lim, abs_lim]);
  1485. vert_data_GA_no_PT = BoSurfStatMakeParcelData(parcel_data_GA_no_PT, wk40, wk40_economo_ind);
  1486. figure; SurfStatViewData(vert_data_GA_no_PT, wk40, strcat(moments(mom), ' GA - no PT'));
  1487. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1488. SurfStatColLim([-abs_lim, abs_lim]);
  1489. vert_data_PNA = BoSurfStatMakeParcelData(parcel_data_PNA, wk40, wk40_economo_ind);
  1490. figure; SurfStatViewData(vert_data_PNA, wk40, strcat(moments(mom), ' PNA'));
  1491. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1492. SurfStatColLim([-abs_lim, abs_lim]);
  1493. vert_data_PNA_no_PT = BoSurfStatMakeParcelData(parcel_data_PNA_no_PT, wk40, wk40_economo_ind);
  1494. figure; SurfStatViewData(vert_data_PNA_no_PT, wk40, strcat(moments(mom), ' PNA - no PT'));
  1495. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1496. SurfStatColLim([-abs_lim, abs_lim]);
  1498. % replace -100 values with NaNs for the correlation
  1499. parcel_data_PMA(parcel_data_PMA==-100)=NaN;
  1500. parcel_data_PMA_no_PT(parcel_data_PMA_no_PT==-100)=NaN;
  1501. parcel_data_GA(parcel_data_GA==-100)=NaN;
  1502. parcel_data_GA_no_PT(parcel_data_GA_no_PT==-100)=NaN;
  1503. parcel_data_PNA(parcel_data_PNA==-100)=NaN;
  1504. parcel_data_PNA_no_PT(parcel_data_PNA_no_PT==-100)=NaN;
  1505. spun_map_PMA = parcel_data_PMA(perm_id); % spin one of the PMA maps of effects
  1506. spun_map_GA = parcel_data_GA(perm_id); % spin one of the GA maps of effects
  1507. spun_map_PNA = parcel_data_PNA(perm_id); % spin one of the PNA maps of effects
  1508. r_emp_PMA.Correlation(mom)=corr(parcel_data_PMA, parcel_data_PMA_no_PT, ...
  1509. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the PMA effects
  1510. r_emp_GA.Correlation(mom)=corr(parcel_data_GA, parcel_data_GA_no_PT, ...
  1511. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the GA effects
  1512. r_emp_PNA.Correlation(mom)=corr(parcel_data_PNA, parcel_data_PNA_no_PT, ...
  1513. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the PNA effects
  1514. % run permutation test
  1515. r_perm_PMA = zeros(n_perm,1);
  1516. r_perm_GA = zeros(n_perm,1);
  1517. r_perm_PNA = zeros(n_perm,1);
  1518. for kk=1:n_perm
  1519. r_perm_PMA(kk) = corr(spun_map_PMA(:, kk), parcel_data_PMA_no_PT, ...
  1520. 'Rows', 'complete', 'Type', 'Pearson');
  1521. r_perm_GA(kk) = corr(spun_map_GA(:, kk), parcel_data_GA_no_PT, ...
  1522. 'Rows', 'complete', 'Type', 'Pearson');
  1523. r_perm_PNA(kk) = corr(spun_map_PNA(:, kk), parcel_data_PNA_no_PT, ...
  1524. 'Rows', 'complete', 'Type', 'Pearson');
  1525. end
  1526. r_emp_PMA.p_value(mom) = sum(abs(r_perm(:)) ...
  1527. >= abs(r_emp_PMA.Correlation(mom))) / n_perm; % p-value after spin test
  1528. r_emp_GA.p_value(mom) = sum(abs(r_perm(:)) ...
  1529. >= abs(r_emp_GA.Correlation(mom))) / n_perm; % p-value after spin test
  1530. r_emp_PNA.p_value(mom) = sum(abs(r_perm(:)) ...
  1531. >= abs(r_emp_PNA.Correlation(mom))) / n_perm; % p-value after spin test
  1532. end
  1534. end
  1536. for supp_figure_7 = 1
  1538. moments = ["Centre of gravity", "Variance"];
  1539. % initialise table to store the correlations between the effect maps on
  1540. % Von Economo's parcels with the effect maps on Schaefer-200 parcels
  1541. r_emp = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1542. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1544. for mom = 1:length(moments)
  1545. if mom == 1
  1546. cog_parcel_data_PMA = moment_age_model("cog", "PMA", ...
  1547. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1548. cog_parcel_data_PMA_schaefer200 = moment_age_model("cog", "PMA", ...
  1549. gen_table_schaefer200, uparc_schaefer200, ...
  1550. total_num_parcels_schaefer200, valid_parcels_schaefer200, 0, 0);
  1551. else
  1552. var_parcel_data_PMA = moment_age_model("variance", "PMA", ...
  1553. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1554. var_parcel_data_PMA_schaefer200 = moment_age_model("variance", "PMA", ...
  1555. gen_table_schaefer200, uparc_schaefer200, ...
  1556. total_num_parcels_schaefer200, valid_parcels_schaefer200, 0, 0);
  1557. end
  1558. end
  1560. % calculate maximum absolute value for limits
  1561. min_val = min([min(cog_parcel_data_PMA(cog_parcel_data_PMA > -100)); ...
  1562. min(cog_parcel_data_PMA_schaefer200(cog_parcel_data_PMA_schaefer200 > -100));
  1563. min(var_parcel_data_PMA(var_parcel_data_PMA > -100)); ...
  1564. min(var_parcel_data_PMA_schaefer200(var_parcel_data_PMA_schaefer200 > -100))]);
  1565. max_val = max([max(cog_parcel_data_PMA(cog_parcel_data_PMA < 100)); ...
  1566. max(cog_parcel_data_PMA_schaefer200(cog_parcel_data_PMA_schaefer200 < 100)); ...
  1567. max(var_parcel_data_PMA(var_parcel_data_PMA < 100)); ...
  1568. max(var_parcel_data_PMA_schaefer200(var_parcel_data_PMA_schaefer200 < 100))]);
  1569. abs_lim = max(abs(min_val), abs(max_val));
  1570. abs_lim = abs_lim + abs_lim/100;
  1572. % plot the significant effects of PMA on each moment
  1573. for mom = 1:length(moments)
  1574. if mom == 1
  1575. parc_data_PMA = cog_parcel_data_PMA;
  1576. parc_data_PMA_schaefer200 = cog_parcel_data_PMA_schaefer200;
  1577. else
  1578. parc_data_PMA = var_parcel_data_PMA;
  1579. parc_data_PMA_schaefer200 = var_parcel_data_PMA_schaefer200;
  1580. end
  1582. vert_data_PMA = BoSurfStatMakeParcelData(parc_data_PMA, wk40, wk40_economo_ind);
  1583. figure; SurfStatViewData(vert_data_PMA, wk40, strcat(moments(mom), ' PMA'));
  1584. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1585. SurfStatColLim([-abs_lim, abs_lim]);
  1586. vert_data_PMA_schaefer200 = BoSurfStatMakeParcelData(parc_data_PMA_schaefer200, wk40, wk40_schaefer200_ind);
  1587. vert_data_PMA_schaefer200(vert_data_PMA == -100) = -100;
  1588. figure; SurfStatViewData(vert_data_PMA_schaefer200, wk40, strcat(moments(mom), ' PMA - Schaefer-200'));
  1589. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1590. SurfStatColLim([-abs_lim, abs_lim]);
  1592. % replace -100 values with NaNs for the correlation
  1593. parc_data_PMA(parc_data_PMA==-100)=NaN;
  1594. vert_data_PMA(vert_data_PMA==-100)=NaN;
  1595. vert_data_PMA_schaefer200(vert_data_PMA_schaefer200==-100)=NaN;
  1596. spun_map = parc_data_PMA(perm_id); % spin one of the maps of effects
  1597. r_emp.Correlation(mom)=corr(vert_data_PMA', vert_data_PMA_schaefer200', ...
  1598. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the effects
  1599. % run permutation test
  1600. r_perm = zeros(n_perm,1);
  1601. for kk=1:n_perm
  1602. vert_spun_map = BoSurfStatMakeParcelData(spun_map(:, kk), wk40, wk40_economo_ind);
  1603. r_perm(kk) = corr(vert_spun_map', vert_data_PMA_schaefer200', ...
  1604. 'Rows', 'complete', 'Type', 'Pearson');
  1605. end
  1606. r_emp.p_value(mom) = sum(abs(r_perm(:)) >= abs(r_emp.Correlation(mom))) / n_perm; % p-value after spin test
  1607. end
  1609. end
  1611. for supp_figure_8_9 = 1
  1613. moments = ["Centre of gravity", "Variance"];
  1614. % initialise tables to store the correlations between the effect maps on
  1615. % Von Economo's parcels with the effect maps on Schaefer-200 parcels
  1616. r_emp_GA = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1617. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1618. r_emp_PNA = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1619. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1621. for mom = 1:length(moments)
  1622. if mom == 1
  1623. cog_parcel_data = moment_age_model("cog", ["GA", "PNA"], ...
  1624. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1625. cog_parcel_data_schaefer200 = moment_age_model("cog", ["GA", "PNA"], ...
  1626. gen_table_schaefer200, uparc_schaefer200, ...
  1627. total_num_parcels_schaefer200, valid_parcels_schaefer200, 0, 0);
  1628. else
  1629. var_parcel_data = moment_age_model("variance", ["GA", "PNA"], ...
  1630. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1631. var_parcel_data_schaefer200 = moment_age_model("variance", ["GA", "PNA"], ...
  1632. gen_table_schaefer200, uparc_schaefer200, ...
  1633. total_num_parcels_schaefer200, valid_parcels_schaefer200, 0, 0);
  1634. end
  1635. end
  1637. % calculate maximum absolute value for limits
  1638. min_val_GA = min([min(cog_parcel_data(cog_parcel_data(:,1) > -100, 1)); ...
  1639. min(cog_parcel_data_schaefer200(cog_parcel_data_schaefer200(:,1) > -100, 1)); ...
  1640. min(var_parcel_data(var_parcel_data(:,1) > -100, 1)); ...
  1641. min(var_parcel_data_schaefer200(var_parcel_data_schaefer200(:,1) > -100, 1))]);
  1642. max_val_GA = max([max(cog_parcel_data(cog_parcel_data(:,1) < 100, 1)); ...
  1643. max(cog_parcel_data_schaefer200(cog_parcel_data_schaefer200(:,1) < 100, 1)); ...
  1644. max(var_parcel_data(var_parcel_data(:,1) < 100, 1)); ...
  1645. max(var_parcel_data_schaefer200(var_parcel_data_schaefer200(:,1) < 100, 1))]);
  1646. abs_lim_GA = max([abs(min_val_GA), abs(max_val_GA)]);
  1647. abs_lim_GA = abs_lim_GA + abs_lim_GA/100;
  1649. min_val_PNA = min([min(cog_parcel_data(cog_parcel_data(:,2) > -100, 2)); ...
  1650. min(cog_parcel_data_schaefer200(cog_parcel_data_schaefer200(:,2) > -100, 2)); ...
  1651. min(var_parcel_data(var_parcel_data(:,2) > -100, 2)); ...
  1652. min(var_parcel_data_schaefer200(var_parcel_data_schaefer200(:,2) > -100, 2))]);
  1653. max_val_PNA = max([max(cog_parcel_data(cog_parcel_data(:,2) < 100, 2)); ...
  1654. max(cog_parcel_data_schaefer200(cog_parcel_data_schaefer200(:,2) < 100, 2)); ...
  1655. max(var_parcel_data(var_parcel_data(:,2) < 100, 2)); ...
  1656. max(var_parcel_data_schaefer200(var_parcel_data_schaefer200(:,2) < 100, 2))]);
  1657. abs_lim_PNA = max([abs(min_val_PNA), abs(max_val_PNA)]);
  1658. abs_lim_PNA = abs_lim_PNA + abs_lim_PNA/100;
  1660. % plot the significant effects of PMA on each moment
  1661. for mom = 1:length(moments)
  1662. if mom == 1
  1663. parcel_data_GA = cog_parcel_data(:,1);
  1664. parcel_data_GA_schaefer200 = cog_parcel_data_schaefer200(:,1);
  1665. parcel_data_PNA = cog_parcel_data(:,2);
  1666. parcel_data_PNA_schaefer200 = cog_parcel_data_schaefer200(:,2);
  1667. else
  1668. parcel_data_GA = var_parcel_data(:,1);
  1669. parcel_data_GA_schaefer200 = var_parcel_data_schaefer200(:,1);
  1670. parcel_data_PNA = var_parcel_data(:,2);
  1671. parcel_data_PNA_schaefer200 = var_parcel_data_schaefer200(:,2);
  1672. end
  1674. vert_data_GA = BoSurfStatMakeParcelData(parcel_data_GA, wk40, wk40_economo_ind);
  1675. figure; SurfStatViewData(vert_data_GA, wk40, strcat(moments(mom), ' GA'));
  1676. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1677. SurfStatColLim([-abs_lim_GA, abs_lim_GA]);
  1678. vert_data_GA_schaefer200 = BoSurfStatMakeParcelData(parcel_data_GA_schaefer200, wk40, wk40_schaefer200_ind);
  1679. vert_data_GA_schaefer200(vert_data_GA == -100) = -100;
  1680. figure; SurfStatViewData(vert_data_GA_schaefer200, wk40, strcat(moments(mom), ' GA - Schaefer-200'));
  1681. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1682. SurfStatColLim([-abs_lim_GA, abs_lim_GA]);
  1683. vert_data_PNA = BoSurfStatMakeParcelData(parcel_data_PNA, wk40, wk40_economo_ind);
  1684. figure; SurfStatViewData(vert_data_PNA, wk40, strcat(moments(mom), ' PNA'));
  1685. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1686. SurfStatColLim([-abs_lim_PNA, abs_lim_PNA]);
  1687. vert_data_PNA_schaefer200 = BoSurfStatMakeParcelData(parcel_data_PNA_schaefer200, wk40, wk40_schaefer200_ind);
  1688. vert_data_PNA_schaefer200(vert_data_PNA == -100) = -100;
  1689. figure; SurfStatViewData(vert_data_PNA_schaefer200, wk40, strcat(moments(mom), ' PNA - Schaefer-200'));
  1690. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1691. SurfStatColLim([-abs_lim_PNA, abs_lim_PNA]);
  1693. % replace -100 values with NaNs for the correlation
  1694. parcel_data_GA(parcel_data_GA==-100)=NaN;
  1695. vert_data_GA(vert_data_GA==-100)=NaN;
  1696. vert_data_GA_schaefer200(vert_data_GA_schaefer200==-100)=NaN;
  1697. parcel_data_PNA(parcel_data_PNA==-100)=NaN;
  1698. vert_data_PNA(vert_data_PNA==-100)=NaN;
  1699. vert_data_PNA_schaefer200(vert_data_PNA_schaefer200==-100)=NaN;
  1700. spun_map_GA = parcel_data_GA(perm_id); % spin one of the GA maps of effects
  1701. spun_map_PNA = parcel_data_PNA(perm_id); % spin one of the PNA maps of effects
  1702. r_emp_GA.Correlation(mom)=corr(vert_data_GA', vert_data_GA_schaefer200', ...
  1703. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the effects
  1704. r_emp_PNA.Correlation(mom)=corr(vert_data_PNA', vert_data_PNA_schaefer200', ...
  1705. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the effects
  1706. % run permutation test
  1707. r_perm_GA = zeros(n_perm,1);
  1708. r_perm_PNA = zeros(n_perm,1);
  1709. for kk=1:n_perm
  1710. vert_spun_map_GA = BoSurfStatMakeParcelData(spun_map_GA(:, kk), wk40, wk40_economo_ind);
  1711. vert_spun_map_PNA = BoSurfStatMakeParcelData(spun_map_PNA(:, kk), wk40, wk40_economo_ind);
  1712. r_perm_GA(kk) = corr(vert_spun_map_GA', vert_data_GA_schaefer200', ...
  1713. 'Rows', 'complete', 'Type', 'Pearson');
  1714. r_perm_PNA(kk) = corr(vert_spun_map_PNA', vert_data_PNA_schaefer200', ...
  1715. 'Rows', 'complete', 'Type', 'Pearson');
  1716. end
  1717. r_emp_GA.p_value(mom) = sum(abs(r_perm(:)) ...
  1718. >= abs(r_emp_GA.Correlation(mom))) / n_perm; % p-value after spin test
  1719. r_emp_PNA.p_value(mom) = sum(abs(r_perm(:)) ...
  1720. >= abs(r_emp_PNA.Correlation(mom))) / n_perm; % p-value after spin test
  1721. end
  1723. end
  1725. for supp_figure_10 = 1
  1727. min_val = 0.7;
  1728. max_val = 2.1;
  1729. for depth = 1:height(MP_parc_economo)
  1730. % average intensities per parcel for each depth
  1731. avg_intensities = mean(MP_parc_economo(depth, :, :), 3);
  1732. % mask non-cortical and limbic parcels
  1733. avg_intensities(~ismember(1:88, valid_parcels_economo)) = -100;
  1734. vert_avg_intensities = BoSurfStatMakeParcelData(avg_intensities, wk40, wk40_economo_ind);
  1735. figure; SurfStatViewData(vert_avg_intensities, wk40, strcat('Depth: ', string(depth)));
  1736. colormap([0.7 0.7 0.7; flipud(lajolla)]);
  1737. SurfStatColLim([min_val, max_val]);
  1738. end
  1740. end
  1742. for supp_figure_11 = 1
  1744. moments = ["Centre of gravity", "Variance"];
  1745. % run models for GA and PNA as predictors separately
  1746. for mom = 1:length(moments)
  1747. if mom == 1
  1748. cog_parcel_data_GA = moment_age_model("cog", "GA", ...
  1749. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1750. cog_parcel_data_PNA = moment_age_model("cog", "PNA", ...
  1751. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1752. else
  1753. var_parcel_data_GA = moment_age_model("variance", "GA", ...
  1754. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1755. var_parcel_data_PNA = moment_age_model("variance", "PNA", ...
  1756. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1757. end
  1758. end
  1760. % calculate maximum absolute value for limits for each age-metric separately
  1761. min_val = min([min(cog_parcel_data_GA(cog_parcel_data_GA > -100, 1)); ...
  1762. min(cog_parcel_data_PNA(cog_parcel_data_PNA > -100, 1)); ...
  1763. min(var_parcel_data_GA(var_parcel_data_GA > -100, 1)); ...
  1764. min(var_parcel_data_PNA(var_parcel_data_PNA > -100, 1))]);
  1765. max_val = max([max(cog_parcel_data_GA(cog_parcel_data_GA(:,1) < 100, 1)); ...
  1766. max(cog_parcel_data_PNA(cog_parcel_data_PNA(:,1) < 100, 1)); ...
  1767. max(var_parcel_data_GA(var_parcel_data_GA(:,1) < 100, 1)); ...
  1768. max(var_parcel_data_PNA(var_parcel_data_PNA(:,1) < 100, 1))]);
  1769. abs_lim = max([abs(min_val), abs(max_val)]);
  1770. abs_lim = abs_lim + abs_lim/100;
  1772. % plot the significant effects of GA and PNA on each moment
  1773. for mom = 1:length(moments)
  1774. if mom == 1
  1775. parcel_data_GA = cog_parcel_data_GA;
  1776. parcel_data_PNA = cog_parcel_data_PNA;
  1777. else
  1778. parcel_data_GA = var_parcel_data_GA;
  1779. parcel_data_PNA = var_parcel_data_PNA;
  1780. end
  1782. vert_data_GA = BoSurfStatMakeParcelData(parcel_data_GA, wk40, wk40_economo_ind);
  1783. figure; SurfStatViewData(vert_data_GA, wk40, strcat(moments(mom), ' GA'));
  1784. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1785. SurfStatColLim([-abs_lim, abs_lim]);
  1786. vert_data_PNA = BoSurfStatMakeParcelData(parcel_data_PNA, wk40, wk40_economo_ind);
  1787. figure; SurfStatViewData(vert_data_PNA, wk40, strcat(moments(mom), ' PNA'));
  1788. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1789. SurfStatColLim([-abs_lim, abs_lim]);
  1790. end
  1792. end
  1794. for supp_figure_12 = 1
  1796. moments = ["Centre of gravity", "Variance"];
  1797. % initialise table to store the correlations between the effect maps
  1798. % with and without correcting for cortical thickness
  1799. r_emp = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1800. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1802. for mom = 1:length(moments)
  1803. if mom == 1
  1804. cog_parcel_data = moment_age_model("cog", "PMA", ...
  1805. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1806. cog_parcel_data_thick_corr = moment_age_model("cog", "PMA", ...
  1807. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 1);
  1808. else
  1809. var_parcel_data = moment_age_model("variance", "PMA", ...
  1810. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1811. var_parcel_data_thick_corr = moment_age_model("variance", "PMA", ...
  1812. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 1);
  1813. end
  1814. end
  1816. % calculate maximum absolute value for limits
  1817. min_val = min([min(cog_parcel_data(cog_parcel_data > -100)); ...
  1818. min(cog_parcel_data_thick_corr(cog_parcel_data_thick_corr > -100));
  1819. min(var_parcel_data(var_parcel_data > -100)); ...
  1820. min(var_parcel_data_thick_corr(var_parcel_data_thick_corr > -100))]);
  1821. max_val = max([max(cog_parcel_data(cog_parcel_data < 100)); ...
  1822. max(cog_parcel_data_thick_corr(cog_parcel_data_thick_corr < 100)); ...
  1823. max(var_parcel_data(var_parcel_data < 100)); ...
  1824. max(var_parcel_data_thick_corr(var_parcel_data_thick_corr < 100))]);
  1825. abs_lim = max(abs(min_val), abs(max_val));
  1826. abs_lim = abs_lim + abs_lim/100;
  1828. % plot the significant effects of PMA on each moment
  1829. for mom = 1:length(moments)
  1830. if mom == 1
  1831. parcel_data = cog_parcel_data;
  1832. parcel_data_thick_corr = cog_parcel_data_thick_corr;
  1833. else
  1834. parcel_data = var_parcel_data;
  1835. parcel_data_thick_corr = var_parcel_data_thick_corr;
  1836. end
  1838. vert_data = BoSurfStatMakeParcelData(parcel_data, wk40, wk40_economo_ind);
  1839. figure; SurfStatViewData(vert_data, wk40, strcat(moments(mom), ' PMA'));
  1840. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1841. SurfStatColLim([-abs_lim, abs_lim]);
  1842. vert_data_thick_corr = BoSurfStatMakeParcelData(parcel_data_thick_corr, wk40, wk40_economo_ind);
  1843. figure; SurfStatViewData(vert_data_thick_corr, wk40, strcat(moments(mom), ' PMA - thickness correction'));
  1844. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1845. SurfStatColLim([-abs_lim, abs_lim]);
  1847. % replace -100 values with NaNs for the correlation
  1848. parcel_data(parcel_data==-100)=NaN;
  1849. parcel_data_thick_corr(parcel_data_thick_corr==-100)=NaN;
  1850. spun_map = parcel_data(perm_id); % spin one of the maps of effects
  1851. r_emp.Correlation(mom)=corr(parcel_data, parcel_data_thick_corr, ...
  1852. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the effects
  1853. % run permutation test
  1854. r_perm = zeros(n_perm,1);
  1855. for kk=1:n_perm
  1856. r_perm(kk) = corr(spun_map(:, kk), parcel_data_thick_corr, ...
  1857. 'Rows', 'complete', 'Type', 'Pearson');
  1858. end
  1859. r_emp.p_value(mom) = sum(abs(r_perm(:)) >= abs(r_emp.Correlation(mom))) / n_perm; % p-value after spin test
  1860. end
  1862. end
  1864. for supp_figure_13_14 = 1
  1866. moments = ["Centre of gravity", "Variance"];
  1867. % initialise table to store the correlations between the effect maps
  1868. % with and without correcting for cortical thickness
  1869. r_emp_GA = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1870. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1871. r_emp_PNA = table(zeros(numel(moments),1), zeros(numel(moments),1), ...
  1872. 'RowNames', moments, 'VariableNames', ["Correlation", "p_value"]);
  1873. % run models which contain both GA and PNA as predictors (effects
  1874. % controlling for the effects of the other)
  1875. for mom = 1:length(moments)
  1876. if mom == 1
  1877. cog_parcel_data = moment_age_model("cog", ["GA", "PNA"], ...
  1878. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1879. cog_parcel_data_thick_corr = moment_age_model("cog", ["GA", "PNA"], ...
  1880. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 1);
  1881. else
  1882. var_parcel_data = moment_age_model("variance", ["GA", "PNA"], ...
  1883. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 0);
  1884. var_parcel_data_thick_corr = moment_age_model("variance", ["GA", "PNA"], ...
  1885. gen_table_economo, uparc_economo, total_num_parcels_economo, valid_parcels_economo, 0, 1);
  1886. end
  1887. end
  1889. % calculate maximum absolute value for limits for each age-metric separately
  1890. min_val_GA = min([min(cog_parcel_data(cog_parcel_data(:,1) > -100, 1)); ...
  1891. min(cog_parcel_data_thick_corr(cog_parcel_data_thick_corr(:,1) > -100, 1)); ...
  1892. min(var_parcel_data(var_parcel_data(:,1) > -100, 1)); ...
  1893. min(var_parcel_data_thick_corr(var_parcel_data_thick_corr(:,1) > -100, 1))]);
  1894. max_val_GA = max([max(cog_parcel_data(cog_parcel_data(:,1) < 100, 1)); ...
  1895. max(cog_parcel_data_thick_corr(cog_parcel_data_thick_corr(:,1) < 100, 1)); ...
  1896. max(var_parcel_data(var_parcel_data(:,1) < 100, 1)); ...
  1897. max(var_parcel_data_thick_corr(var_parcel_data_thick_corr(:,1) < 100, 1))]);
  1898. abs_lim_GA = max([abs(min_val_GA), abs(max_val_GA)]);
  1899. abs_lim_GA = abs_lim_GA + abs_lim_GA/100;
  1901. min_val_PNA = min([min(cog_parcel_data(cog_parcel_data(:,2) > -100, 2)); ...
  1902. min(cog_parcel_data_thick_corr(cog_parcel_data_thick_corr(:,2) > -100, 2)); ...
  1903. min(var_parcel_data(var_parcel_data(:,2) > -100, 2)); ...
  1904. min(var_parcel_data_thick_corr(var_parcel_data_thick_corr(:,2) > -100, 2))]);
  1905. max_val_PNA = max([max(cog_parcel_data(cog_parcel_data(:,2) < 100, 2)); ...
  1906. max(cog_parcel_data_thick_corr(cog_parcel_data_thick_corr(:,2) < 100, 2)); ...
  1907. max(var_parcel_data(var_parcel_data(:,2) < 100, 2)); ...
  1908. max(var_parcel_data_thick_corr(var_parcel_data_thick_corr(:,2) < 100, 2))]);
  1909. abs_lim_PNA = max([abs(min_val_PNA), abs(max_val_PNA)]);
  1910. abs_lim_PNA = abs_lim_PNA + abs_lim_PNA/100;
  1912. % plot the significant effects of GA and PNA on each moment
  1913. for mom = 1:length(moments)
  1914. if mom == 1
  1915. parcel_data_GA = cog_parcel_data(:,1);
  1916. parcel_data_GA_thick_corr = cog_parcel_data_thick_corr(:,1);
  1917. parcel_data_PNA = cog_parcel_data(:,2);
  1918. parcel_data_PNA_thick_corr = cog_parcel_data_thick_corr(:,2);
  1919. else
  1920. parcel_data_GA = var_parcel_data(:,1);
  1921. parcel_data_GA_thick_corr = var_parcel_data_thick_corr(:,1);
  1922. parcel_data_PNA = var_parcel_data(:,2);
  1923. parcel_data_PNA_thick_corr = var_parcel_data_thick_corr(:,2);
  1924. end
  1926. vert_data_PNA = BoSurfStatMakeParcelData(parcel_data_PNA, wk40, wk40_economo_ind);
  1927. figure; SurfStatViewData(vert_data_PNA, wk40, strcat(moments(mom), ' PNA'));
  1928. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1929. SurfStatColLim([-abs_lim_PNA, abs_lim_PNA]);
  1930. vert_data_PNA_thick_corr = BoSurfStatMakeParcelData(parcel_data_PNA_thick_corr, wk40, wk40_economo_ind);
  1931. figure; SurfStatViewData(vert_data_PNA_thick_corr, wk40, strcat(moments(mom), ' PNA - thickness correction'));
  1932. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1933. SurfStatColLim([-abs_lim_PNA, abs_lim_PNA]);
  1934. vert_data_GA = BoSurfStatMakeParcelData(parcel_data_GA, wk40, wk40_economo_ind);
  1935. figure; SurfStatViewData(vert_data_GA, wk40, strcat(moments(mom), ' GA'));
  1936. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1937. SurfStatColLim([-abs_lim_GA, abs_lim_GA]);
  1938. vert_data_GA_thick_corr = BoSurfStatMakeParcelData(parcel_data_GA_thick_corr, wk40, wk40_economo_ind);
  1939. figure; SurfStatViewData(vert_data_GA_thick_corr, wk40, strcat(moments(mom), ' GA - thickness correction'));
  1940. colormap([0.5 0.5 0.5; vik; 0.5 0.5 0.5]);
  1941. SurfStatColLim([-abs_lim_GA, abs_lim_GA]);
  1943. % replace -100 values with NaNs for the correlation
  1944. parcel_data_GA(parcel_data_GA==-100)=NaN;
  1945. parcel_data_GA_thick_corr(parcel_data_GA_thick_corr==-100)=NaN;
  1946. parcel_data_PNA(parcel_data_PNA==-100)=NaN;
  1947. parcel_data_PNA_thick_corr(parcel_data_PNA_thick_corr==-100)=NaN;
  1948. spun_map_GA = parcel_data_GA(perm_id); % spin one of the GA maps of effects
  1949. spun_map_PNA = parcel_data_PNA(perm_id); % spin one of the PNA maps of effects
  1950. r_emp_GA.Correlation(mom)=corr(parcel_data_GA, parcel_data_GA_thick_corr, ...
  1951. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the GA effects
  1952. r_emp_PNA.Correlation(mom)=corr(parcel_data_PNA, parcel_data_PNA_thick_corr, ...
  1953. 'Rows', 'complete', 'Type', 'Pearson'); % correlate the PNA effects
  1954. % run permutation test
  1955. r_perm_GA = zeros(n_perm,1);
  1956. r_perm_PNA = zeros(n_perm,1);
  1957. for kk=1:n_perm
  1958. r_perm_GA(kk) = corr(spun_map_GA(:, kk), parcel_data_GA_thick_corr, ...
  1959. 'Rows', 'complete', 'Type', 'Pearson');
  1960. r_perm_PNA(kk) = corr(spun_map_PNA(:, kk), parcel_data_PNA_thick_corr, ...
  1961. 'Rows', 'complete', 'Type', 'Pearson');
  1962. end
  1963. r_emp_GA.p_value(mom) = sum(abs(r_perm_GA(:)) ...
  1964. >= abs(r_emp_GA.Correlation(mom))) / n_perm; % p-value after spin test
  1965. r_emp_PNA.p_value(mom) = sum(abs(r_perm_PNA(:)) ...
  1966. >= abs(r_emp_PNA.Correlation(mom))) / n_perm; % p-value after spin test
  1967. end
  1969. end

infant_micro_figures.m at commit 5c4f80c, no license · at the source

Overview

Authors: Thanos Tsigaras1,2, Juergen Dukart1,2, Stuart Oldham3,4, Simon B. Eickhoff1,2, Casey Paquola1
ORCID iDs: Casey Paquola
  1. Institute for Neuroscience and Medicine (INM-7), Forschungszentrum Jülich, Jülich, Germany
  2. Institute for Systems Neuroscience, Heinrich Heine University Düsseldorf, Düsseldorf, Germany
  3. Developmental Imaging, Murdoch Children’s Research Institute, Parkville, Australia
  4. Turner Institute for Brain and Mental Health, Monash University, Melbourne, Australia
Institutions: Forschungszentrum Jülich (Germany); Heinrich Heine University Düsseldorf (Germany); Murdoch Children's Research Institute (Australia); Monash University (Australia)
Journal: PLoS biology, volume 24, issue 3, article e3003722
Dates: received 12 August 2025; accepted 9 March 2026; published online 26 March 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pbio.3003722 · PMID 41886498 · PMCID PMC13046243 · OpenAlex W7140585155
Open access: gold, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: structural MRI / diffusion (modality), human (organism), developmental (subfield)
Methods: Connectivity, Statistics, Preprocessing, Machine learning, fMRI & imaging
MeSH: Cerebral Cortex*, Myelin Sheath*, Female, Gestational Age, Humans, Infant, Newborn, Magnetic Resonance Imaging, Male, Neurodevelopment, Pregnancy (* major topic)
Journal subjects: Physical Sciences, Materials Science, Materials Physics, Microstructure, Physics, Medicine and Health Sciences, Diagnostic Medicine, Diagnostic Radiology, Magnetic Resonance Imaging, Research and Analysis Methods, Imaging Techniques, Radiology and Imaging, Mathematical and Statistical Techniques, Statistical Methods, Regression Analysis, Linear Regression Analysis, Mathematics, Statistics, Discrete Mathematics, Combinatorics, Permutation, Biology and Life Sciences, Neuroscience, Cognitive Science, Cognitive Psychology, Perception, Sensory Perception, Psychology, Social Sciences, Developmental Biology, Twins, Women's Health, Maternal Health, Pregnancy, Obstetrics and Gynecology, People and Places, Population Groupings, Age Groups, Children, Infants, Families
Topic: Advanced Neuroimaging Techniques and Applications (Radiology, Nuclear Medicine and Imaging, Medicine), according to OpenAlex
Funding: Deutsche Forschungsgemeinschaft (524408221); German Scholars Organisation (KT40)
Citations: cited by 1 paper (Europe PMC); 70 references in the paper

Abstract

During late gestation and early postnatal development a combination of intrinsic and extrinsic factors drive the maturation of the human cortex. This process is regionally heterogeneous, with cortical areas developing at different paces and trajectories. Leveraging submillimetre T1-weighted/T2w-weighted (T1w/T2w) magnetic resonance imaging (MRI) from pre- and full-term neonates (n = 599, 26–44 weeks), we sampled intracortical profiles across the cortex and characterized the profiles’ shapes according to their central moments. We found that gestational age at birth dominated the effects on early cortical development, with significant, global increases in intracortical homogeneity and a bimodal change in the balance of myelin-sensitive signal between superficial and deeper cortical layers. On the other hand, weeks since birth (i.e., postnatal age) exhibited different effects on myelin, with increasing intracortical heterogeneity and intracortical balance only shifting towards deeper layers in posterior temporal, occipital, medial parietal areas and some prefrontal areas. These effects align with low spatial-frequency geometric eigenmodes of the human cortex, specifically the anterior–posterior and superior-inferior axes. Our findings demonstrate that separating prenatal from postnatal influences, and analyzing intracortical profiles rather than macroscale features, provides finer-grained insights into how human cortical myelin changes during perinatal development and lays the groundwork for investigating the biological underpinnings that govern normative cortical maturation.

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

Repositories

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

ThanosTsig/InfantMicro

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 5c4f80c56e321abb1081aa388d2cbac5ecd67fb4, 11 March 2026
Languages: MATLAB (7)
Size: 45 files, 7 scripts
Software Heritage: not archived
Found in: “Data Availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Statistics and Machine Learning Toolbox (2 files), fdr_bh (Benjamini-Hochberg FDR) (1 file), GIfTI library for MATLAB (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
8 files
At the source: github.com/ThanosTsig/InfantMicro

Zenodo 18910237

License: CC-BY-4.0
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the text, “Alignment of developmental effects on myelin to ”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Statistics and Machine Learning Toolbox (2 files), fdr_bh (Benjamini-Hochberg FDR) (1 file), GIfTI library for MATLAB (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers (HTTP 200)
  • 30 September 2026: the link answers (HTTP 200)
8 files
At the source: doi.org/10.5281/zenodo.18910237

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

Tracing map

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

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 14 scripts, each with its path and the digest of its content;
  • 4 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data Availability

Data were obtained from the developing Human Connectome Project (https://www.developingconnectome.org/), which is openly available to download for research purposes via the NIMH data archive (https://nda.nih.gov/edit_collection.html?id=3955). Code and additional data for this study are openly available at a Github repository (https://github.com/ThanosTsig/InfantMicro) which is archived in Zenodo (https://doi.org/10.5281/zenodo.18910237).

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

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 10 MeSH terms, 2 funders, 69 references.

Cite

This paper

Tsigaras, T., Dukart, J., Oldham, S., Eickhoff, S. B., & Paquola, C. (2026). Prenatal and early postnatal periods differentially shape the maturation of human cortical microstructure and myelin. PLoS biology, 24(3), e3003722. https://doi.org/10.1371/journal.pbio.3003722

BibTeX

@article{tsigaras2026prenatal,
author = {Tsigaras, Thanos and Dukart, Juergen and Oldham, Stuart and Eickhoff, Simon B. and Paquola, Casey},
title = {{Prenatal and early postnatal periods differentially shape the maturation of human cortical microstructure and myelin}},
journal = {PLoS biology},
year = {2026},
month = mar,
volume = {24},
number = {3},
pages = {e3003722},
publisher = {PLOS},
issn = {1544-9173},
doi = {10.1371/journal.pbio.3003722},
url = {https://doi.org/10.1371/journal.pbio.3003722},
pmid = {41886498},
pmcid = {PMC13046243}
}

RIS

TY - JOUR
AU - Tsigaras, Thanos
AU - Dukart, Juergen
AU - Oldham, Stuart
AU - Eickhoff, Simon B.
AU - Paquola, Casey
TI - Prenatal and early postnatal periods differentially shape the maturation of human cortical microstructure and myelin
T2 - PLoS biology
J2 - PLoS Biol
PY - 2026
DA - 2026/03/26
VL - 24
IS - 3
SP - e3003722
SN - 1544-9173
PB - PLOS
DO - 10.1371/journal.pbio.3003722
UR - https://doi.org/10.1371/journal.pbio.3003722
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pbio.3003722",
"type": "article-journal",
"title": "Prenatal and early postnatal periods differentially shape the maturation of human cortical microstructure and myelin",
"container-title": "PLoS biology",
"author": [
{
"family": "Tsigaras",
"given": "Thanos"
},
{
"family": "Dukart",
"given": "Juergen"
},
{
"family": "Oldham",
"given": "Stuart"
},
{
"family": "Eickhoff",
"given": "Simon B."
},
{
"family": "Paquola",
"given": "Casey"
}
],
"container-title-short": "PLoS Biol",
"volume": "24",
"issue": "3",
"page": "e3003722",
"DOI": "10.1371/journal.pbio.3003722",
"PMID": "41886498",
"PMCID": "PMC13046243",
"ISSN": "1544-9173",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pbio.3003722",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
26
]
]
}
}

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.1002/hbm.26660 [code]
Feature similarity gradients detect alterations in the neonatal cortex associated with preterm birth.
Journal: Human brain mapping
In common: biomedia.github.io/dhcp-release-notes, developmental, structural MRI / diffusion, 12 references
[2] doi:10.1162/imag.a.1212 [code]
Intracortical microstructure profiling: A cross-modal method for indexing cortical lamination.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: structural MRI / diffusion, 9 references, author Casey Paquola
[3] doi:10.1371/journal.pbio.3003378 [code]
Humans have a longer period of cortical maturation across depth and hierarchy than macaques.
Journal: PLoS biology
In common: developmental, structural MRI / diffusion, 12 references
[4] doi:10.1038/s41593-024-01868-0 [code]
The architecture of the human default mode network explored through cytoarchitecture, wiring and signal flow.
Journal: Nature neuroscience
In common: GIfTI library for MATLAB, Statistics and Machine Learning Toolbox, structural MRI / diffusion, 8 references, author Casey Paquola
[5] doi:10.1038/s41467-024-51459-7 [code]
Relating sex-bias in human cortical and hippocampal microstructure to sex hormones.
Journal: Nature communications
In common: fdr_bh (Benjamini-Hochberg FDR), GIfTI library for MATLAB, Statistics and Machine Learning Toolbox, 10 references
[6] doi:10.7554/elife.78733 [code]
Development of visual cortex in human neonates is selectively modified by postnatal experience.
Journal: eLife
In common: GIfTI library for MATLAB, Statistics and Machine Learning Toolbox, structural MRI / diffusion, 10 references
[7] doi:10.1162/imag_a_00418 [code]
Perinatal development of structural thalamocortical connectivity.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: GIfTI library for MATLAB, Statistics and Machine Learning Toolbox, structural MRI / diffusion, 9 references
[8] doi:10.1371/journal.pbio.3003927 [code]
Brain structural and functional connectivity converge prenatally but diverge after birth.
Journal: PLoS biology
In common: GIfTI library for MATLAB, Statistics and Machine Learning Toolbox, biomedia.github.io/dhcp-release-notes, developmental, structural MRI / diffusion, 5 references
[9] doi:10.1038/s41467-026-70660-4 [code]
Dual-axis myelination covariance drives the functional connectivity emergence during infancy.
Journal: Nature communications
In common: Statistics and Machine Learning Toolbox, structural MRI / diffusion, 10 references
[10] doi:10.1016/j.nicl.2022.103286 [code]
Aberrant claustrum structure in preterm-born neonates: an MRI study.
Journal: NeuroImage. Clinical
In common: developmental, structural MRI / diffusion, 10 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1371/journal.pbio.3003722.svg)](https://openscicode.org/badge/doi/10.1371/journal.pbio.3003722)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1371/journal.pbio.3003722.svg
   :target: https://openscicode.org/badge/doi/10.1371/journal.pbio.3003722
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1371/journal.pbio.3003722"><img src="https://openscicode.org/badge/doi/10.1371/journal.pbio.3003722.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge 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: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.